Advertisement

Genotype by environment interactions in fertility traits in New Zealand dairy cows

Open AccessPublished:September 19, 2018DOI:https://doi.org/10.3168/jds.2017-14195

      ABSTRACT

      New Zealand's seasonal dairy farming system entails a condensed calving pattern with cows required to conceive within approximately 12 wk of the planned start of calving. This has resulted in strong selection for fertility through culling of nonpregnant cows and relatively strong emphasis on fertility in Breeding Worth, the national breeding objective that drives sire selection. Despite this, average herd-level fertility is highly variable across New Zealand dairy farms. We studied genotype by environment interaction in fertility-related traits, with the goal of improving selection decisions in different fertility environments. We used data from the New Zealand national dairy database, which contains records on 3,743,862 animals. Herds were classified into high-, mid-, or low-fertility categories or environments based on herd average fertility performance, and data were analyzed in 2 different ways. First, we estimated genetic parameters when the fertility trait was defined specifically for each fertility environment to determine the extent to which genetic correlations between high- and low-fertility environments differed from 1 and the extent of changes in genetic variance across environments. Second, we used simple regression to evaluate the impact of ancestral genetic merit for fertility on cow fertility phenotypes to compare the effect of changes in genetic merit on phenotypic performance between fertility environments. The genetic standard deviations of fertility-related traits were 1.5 to 3.6 times higher in low-fertility herds than in high-fertility herds, and the genetic correlations between the same fertility-related traits between the high- and low-fertility environments were moderate to high, albeit with high standard errors. The high standard errors of the correlations reflected the low heritabilities of the traits and potential problems of culling bias, particularly for traits expressed in later parities. Regression analysis revealed that the bottom 30% of herds (in terms of fertility) could achieve more than twice the benefit from selection for fertility than the top 30% of herds. Although our analyses do not support separate genetic evaluations of fertility in the different environments, they indicate that low-fertility herds could benefit more from targeted selection of sires with higher fertility estimated breeding values than from selection based solely on the multitrait national index. Conversely, high-fertility herds could focus their sire selection on traits other than fertility, provided they avoid very low fertility sires.

      Key words

      INTRODUCTION

      New Zealand's seasonal dairy system requires a condensed calving pattern, with a calving interval of 365 d. Consequently, a cow must be mated and conceive within approximately 12 wk of the commencement of the seasonal mating period or be culled for infertility. This system aims to maximize pasture utilization by aligning high feed demand during lactation with high seasonal pasture growth. High variability in gestation and calving patterns negatively affects the profitability and efficiency of the production system (
      • MacMillan K.L.
      Calving patterns and herd production in seasonal dairy herds.
      ;
      • MacMillan K.L.
      • Taufa V.K.
      • Phillips P.
      Recent trends in conception rates and return patterns in AB herds and their effects on calving patterns.
      ).
      Unfavorable genetic correlations between production and fertility have been reported in many dairy cattle populations (
      • Berry D.P.
      • Buckley F.
      • Dillon P.
      • Evans R.D.
      • Rath M.
      • Veerkamp R.F.
      Genetic relationships among body condition score, body weight, milk yield, and fertility in dairy cows.
      ;
      • Windig J.J.
      • Calus M.P.L.
      • Veerkamp R.F.
      Influence of herd environment on health and fertility and their relationship with milk production.
      ;
      • Holtsmark M.
      • Heringstad B.
      • Madsen P.
      • Ødegård J.
      Genetic relationship between culling, milk production, fertility, and health traits in Norwegian red cows.
      ). As a consequence of this antagonistic relationship, intensive selection for production traits has decreased cow fertility (
      • Grosshans T.
      • Xu Z.Z.
      • Burton L.J.
      • Johnson D.L.
      • Macmillan K.L.
      Performance and genetic parameters for fertility of seasonal dairy cows in New Zealand.
      ;

      Harris, B. L., A. M. Winkelman, and L. J. Burton. 2000. Comparisons of fertility measures in strains of Holstein-Friesian cows and their crosses. Massey Dairy Farming Annu. Proc. Massey Dairy Farmers Conference.

      ;
      • Harris B.L.
      • Pryce J.E.
      Genetic and phenotypic relationships between milk protien percentage, reproductive performance and body condition score in New Zealand dairy cattle.
      ). In New Zealand,
      • Macdonald K.A.
      • Thorrold B.S.
      • Glassey C.B.
      • Holmes C.W.
      • Pryce J.E.
      Impact of farm management decision rules on the production and profit of different strains of Holstein-Friesian dairy cows.
      reported that the 1990s-type Holstein-Friesian with predominantly New Zealand ancestry had reduced fertility compared with the 1970s type. This aligns with research from the United States, where conception rates at first insemination have decreased by 0.45% per year (
      • Beam S.W.
      • Butler W.R.
      Effects of energy balance on follicular development and first ovulation in postpartum dairy cows.
      ), and the average number of inseminations required for conception has increased from 1.75 to more than 3 over a 20-yr period (
      • Lucy M.C.
      Reproductive loss in high-producing dairy cattle: Where will it end?.
      ).
      Costs related to reduced fertility extend beyond the increased expenses for mating, mating management, and lost milk production due to later calving dates. Fertility also affects cow survival, or longevity, due to fertility-based culling (
      • Philipsson J.
      • Lindhé B.
      Experiences of including reproduction and health traits in Scandinavian dairy cattle breeding programmes.
      ). Lengthening the productive lives of dairy cattle has become a priority globally because selection for longer life and improved fertility reduces wastage and replacement rates in dairy herds, with corresponding economic impacts (
      • van Arendonk J.A.M.
      Use of profit equations to determine relative economic value of dairy cattle herd life and production from field data.
      ).
      • Pritchard T.
      • Coffey M.
      • Mrode R.
      • Wall E.
      Understanding the genetics of survival in dairy cows.
      argued that cows with longer lives are expected to be more profitable because fewer replacement heifers are required to achieve the same herd output, more superior offspring are available as replacements, and a greater proportion of potential lactation milk yield is captured.
      In New Zealand, fertility has been included in the Breeding Worth, the national breeding objective, since 2001 (
      • Harris B.L.
      • Montgomerie W.A.
      Fertility breeding values for seasonal dairying.
      ), and selection over the past decade has reversed an unfavorable genetic trend in fertility for the Holstein Friesian breed associated with the integration of overseas Holstein genetics into the New Zealand population. Industry data indicates the phenotypic performance of 3-wk submission rate has improved from 78.5% in the 2008–2009 season to 81.1% in the 2014–2015 dairy season (
      • DairyNZ-Livestock Improvement Corporation
      New Zealand dairy statistics 2014–15.
      ). Similarly, the mean 6-wk in-calf rate improved from 63.4 to 66.8% over the same period. In contrast, the average number of inseminations per cow is still trending upward, with an increase from 1.25 inseminations per cow in the 1997–1998 season to 1.36 in the 2014–2015 season. Although 6-wk in-calf rate has improved, it is a long way from the industry target of 78%, the average of the top quartile of farms presented by
      • Xu Z.
      • Burton L.J.
      Reproductive Performance of Dairy Cows in New Zealand: Monitoring Fertility.
      .
      Despite these recent improvements, herd-level average fertility phenotypes are highly variable in New Zealand. Management is important but genotype by fertility environment interaction (G×E) may also contribute to this variation. Genotype by fertility environment interaction interactions can take 2 forms (
      • Falconer D.S.
      • Mackay T.F.C.
      Introduction to Quantitative Genetics.
      ): a scaling effect, in which environmental factors effect some genotypes more than others but there is no reranking, or a crossing-type interaction with reranking of genotypes in different environments.
      This study had 2 objectives: (1) complete parameter estimation to quantify the effect of G×E in fertility traits, as measured by a genetic correlation <1 between environments, and (2) to quantify variation in expression of genetic merit for fertility across a range of environments by regression of phenotypic fertility performance on a pedigree index of fertility genetic merit. We will discuss how the results could be deployed within the national genetic evaluation system, and how they could be used for selective breeding decisions by commercial farmers.

      MATERIALS AND METHODS

      We obtained fertility and milk production data from the New Zealand national dairy database and added herd-level environmental descriptors based on fertility and milk production level derived from reproductive and milk test-day records. We conducted 2 analyses; first, conventional genetic parameter estimation with trait definition split according to herd average fertility performance; and second, regression analysis to understand how fertility genetic merit translates differently into cow fertility phenotypic performance in different fertility environments.

      Data

      We identified and retained animals born between 1997 and 2009. This interval ensured that all included animals would have expressed fertility traits in up to 4 parities, and provided a minimum of 10 years of records representing modern genotypes.
      We only included fertility phenotypes for cows with a known date of birth (i.e., not recorded as an estimate), born in New Zealand, no missing parent records, and a valid breed code of Jersey (JE), Holstein Friesian (HF), or JE × HF crosses.
      We generated contemporary group identification for each animal by concatenating their unique farm ID (map reference in the database), farmer reference (varies if the cow is managed by different farmers on the same property), and year of calving. We excluded cows from the data set that were either transferred between farms or sold to new farmers but remained on the same property. This ensured that all animals in the data set were on the same farm and managed by the same farmer for all parities, and therefore were under similar management for all phenotypic recording. Herds with fewer than 40 animal records per year were removed from the data set to ensure that outlier herds with small numbers of animals did not skew the distribution of fertility phenotypes.
      We used mating and calving records to derive the binary fertility traits used in the New Zealand genetic evaluation system: calving rate at 42 d (CR42), which was defined as success (1) or failure (0) to calve in the first 42 d after the planned start of calving (PSC) for cows calving in the same contemporary group (herd). Mating dates and planned start of mating dates were used to derive the binary mating trait, percentage mated at 21 d (PM21). Cows received a score of 1 if they presented for mating within the first 21 d of mating for their contemporary group, and 0 if they were mated later or not at all. Although both CR42 and PM21 were scored phenotypically as binomial traits taking values of 0 or 1, we have reported all of their results as percentages by multiplying means and coefficients by 100 before their inclusion in tables of results. A new continuous trait called calving season day (CSD) was also included. This trait treats calving records as a continuous variable (
      • Bowley F.E.
      • Green R.E.
      • Amer P.R.
      • Meier S.
      Novel approaches to genetic analysis of fertility traits in New Zealand dairy cattle.
      ;
      • Amer P.R.
      • Stachowicz K.
      • Jenkins G.M.
      • Meier S.
      Short communication: Estimates of genetic parameters for dairy fertility in New Zealand.
      ) and is the deviation, in days, from the PSC for a contemporary group (herd) and the cow's calving date. For all traits, we added a number (0–3) representing the parity in which it was measured, where 0 is heifer calving and 3 is the animal's fourth calving (see Table 1 for trait abbreviations and definitions).
      Table 1Trait definitions and corresponding trait abbreviations
      TraitTrait descriptionAbbreviation
      Calving rate at 42 d (CR42)CR42 is a binary trait, representing success (1) or failure (0) for calving in the first 42 d of a herd's calving period. This trait is recorded across 4 consecutive calving events, commencing as a first-calving heifer and ending at the fourth calving.CR42_0 (first/heifer-calving trait)
      CR42_1 (second-calving trait)
      CR42_2 (third-calving trait)
      CR42_3 (fourth-calving trait)
      Calving season day (CSD)CSD is a continuous calving trait, which is a deviation, in days, that the animal calved from the planned start of calving.CSD_0 (first/heifer-calving trait)
      CSD_1 (second-calving trait)
      CSD_2 (third-calving trait)
      CSD_3 (fourth-calving trait)
      Percentage mated in the first 21 d (PM21)PM21 is a binary trait, representing success (1) or failure (0) of presenting for mating in the first 21 d of a herd's mating period. This trait is recorded across 3 consecutive mating events, commencing at second mating and ending at the fourth mating event.PM21_1 (second-mating trait)
      PM21_2 (third-mating trait)
      PM21_3 (fourth-mating trait)

      Herd Fertility Descriptor

      To enable G×E analysis, we partitioned herds into herd fertility environments, by calculating a herd fertility environment descriptor, being the herd-level mean of the cows' second calving event (CSD_1). This constitutes the largest cohort of animals of a single age group and those most vulnerable to poor management due to the combined stresses of simultaneously recovering from their first calving, growing, and lactating. The mean CSD_1 performance of the herd was then assigned to all animals in the herd.

      High- and Low-Fertility Herd Status for Parameter Estimation

      We classified low- and high-fertility herds as having a herd average CSD_1 <24 d, and >37 d respectively, which corresponded approximately to the top and bottom 10% for fertility performance; we assigned all other herds to the mid-fertility group. After all data edits, 3,743,862 animals remained in the data set (326,679 in high-, 3,017,535 in mid-, and 399,647 in low-fertility environments). To facilitate timely computation of variance component estimates, we reduced the data further by randomly sampling herds from the 3 fertility environments. This resulted in 73,321, 667,491, and 67,117 animals in high-, mid-, and low-fertility environments, respectively.

      Decile Partitioning for Regression Analysis

      We applied a second partitioning approach using the full data set of approximately 3,750,000 animals, which assigned herds to decile groups 1 to 10 (where 1 represented favorable fertility performance and 10 was unfavorable), based on their mean herd CSD_1. To aid timely completion of regression analysis, we reduced these data by randomly sampling herds in each decile group. The sampling criteria balanced each decile to contain records of approximately 80,000 animals from 500 herds, resulting in a total of 810,000 animals included in the final data set.

      Milk Production Partitioning

      We also partitioned the data set based on milk production by grouping herds as either “high production” or “other,” based on their herd-level milk solids production. We defined milk solids as the combined fat and protein percentage multiplied by milk volume using individual herd test records over each cow's lactation. Herd test records from seasons 2000 to 2014 were available, but we only used records relating to animals in the fertility data set. To reduce computation time and memory requirements, the data set for production-level grouping was filtered to include only second-lactation records for cows that had calved from July to October and daily milk solids yields ranging from 0.5 to 3.0 kg. We also discarded herd tests recorded as being under once-a-day milking regimens. Furthermore, we retained only records between 50 and 200 DIM, when the slope of the lactation curve was linear. Finally, herd-seasons with fewer than 30 unique records were discarded, and a simple random sample (n = 100 animals, no replacement) was taken for herd-years with more than 100 records. For the remaining records, the relationships between daily milk solids (ms) and DIM (dim), percentage Friesian (fr), and herd identifier (herd) were modeled as a linear regression equation separately for each year using the statistical model
      ms=β0+β1dim+β2fr+β3herd+ε,


      where β0, …, β3 are the regression coefficients for the predictor variables, and ε is the error term. All effects were highly significant (P < 0.001) and the coefficient of determination of the 14 models, one for each year, ranged from 0.53 to 0.58 with a mean of 0.56.
      The means of herd coefficients (β3) were computed for sliding 5-yr periods (i.e., 1998–2002, 1999–2003, 2000–2004, and so on). Herds that did not have at least 2 coefficients in a period were excluded. Using 5-yr herd means, as opposed to individual herd-year values, minimized the effect of cows and herds switching between high- and low-production groups over time. Further, the 5-yr periods were chosen so that fertility traits of an individual expressed over up to 5 yr of their lifetime would fall in a unique grouping period, starting with the first and ending with the fifth year of the individual's milk production. Then, average herd coefficients of the 2008–2012 period were arbitrarily set as a baseline. High production systems were defined as those herd-periods in which production was equal to, or above, the top 20% of the baseline. Remaining herd-periods were grouped into the “other” production category.
      Animals in the fertility data set were matched with corresponding herd test records. This resulted in approximately 520,000 animals remaining in the data set for analysis, of which 144,000 animals were identified in high-production herds (from 636 herds) and the remaining 376,000 (from 2,657 herds) belonging to the “other” production group.

      Parameter Estimation of Fertility-Related Traits in High- and Low-Fertility Environments

      For these analyses, we excluded the mid-fertility group from the data and analyzed the fertility traits CSD, CR42, and PM21 in pair-wise bivariate animal models for high- and low-fertility environments using ASREML (
      • Gilmour A.R.
      • Cullis B.R.
      • Welham S.J.
      • Thompson R.
      ASREML Reference Manual. NSW Agric. Biometric Bull.
      ). This approach enables the estimation of across-environment genetic correlations in addition to genetic variances, heritabilities, breed effects, and heterosis effects that are specific to each environment (
      • Falconer D.S.
      • Mackay T.F.C.
      Introduction to Quantitative Genetics.
      ). The following statistical model was used:
      y = CG + age × breed + HFFR + HFNZ + HETJE×NZ + HETJE×FR + HETJE×NZ + RECJE×NZ + RECJE×FR + RECNZ×FR + a + e,


      where y is the phenotypic record, CG is the fixed contemporary group effect corresponding to herd-year combinations, age × breed is the fixed linear regression of age at calving in days nested within breed, HFNZ and HFFR are fixed linear regressions of New Zealand and foreign Holstein-Friesian breed proportions, respectively, HET and REC are fixed linear regressions of breed-specific (heterosis and recombination effects, respectively, a is a random animal effect, and e is the random error term (
      • Harris B.L.
      • Pryce J.E.
      • Xu Z.Z.
      • Montgomerie W.A.
      Fertility breeding values in New Zealand, the next generation.
      ). A pedigree file containing 3 generations of sire and dam records was used to generate the inverse of the numerator relationship matrix for the animal effect required to solve the mixed model equations.
      Breed proportions for each animal were calculated and used to account for heterosis (het) and recombination (rec), using the methodology presented by (
      • Dickerson G.E.
      Inbreeding and heterosis in animals.
      ) as follows:
      het = Ps(A)Pd(B) + Ps(B)Pd(A) and


      rec = Ps(A)Ps(B) + Pd(B)Pd(A),


      where s is sire, d is dam, A is breed 1, B is breed 2, and P is the proportion of breed A or B in sire or dam breed composition.
      Variance components were estimated by analyzing traits in bivariate animal models, where each trait in the low-fertility environment is modeled with the corresponding trait in the high-fertility environment (e.g., CSD1_low with CSD1_high). This approach was then expanded to produce bivariate results among all pairwise trait combinations between and within environments. To evaluate potential biases in parameter estimates for later fertility traits due to selection on fertility in early parities, we completed trivariate analyses, in which traits within environment were fitted in separate models (e.g., CSD_0_low to CSD_2_low in a single model). The fourth-calving traits (i.e., CSD_3 and CR42_3), in both low- and high-fertility herds, were not included to make full 4-trait analyses because of excessive computational time that would have been required.

      Heritability for Binomial Traits on the Underlying Continuous Scale

      Heritability of the binomial fertility traits was estimated on the observable binomial scale. Heritabilities on the observable scales were transformed to the underlying continuous scale (U) using standard formulas as described by
      • Dempster E.R.
      • Lerner I.M.
      Heritability of threshold characters.
      as follows:
      hU2=h21pi2p


      where h2 is heritability on the binomial scale, p is the proportion of animals that did not calve or present for mating within the allotted time (i.e., scoring 0 for CR42 or PM21), and i is the mean for the proportion of the animals to the right of the truncation point that results in p proportion of animals to the right of the truncation point on a standardized normal distribution.

      Pedigree Index Fertility Breeding Value Analysis

      Sire fertility EBV were obtained from New Zealand Animal Evaluation Ltd. (Hamilton) following the October 2014 national genetic evaluation run. The breeding values were generated using the statistical model outlined by
      • Harris B.L.
      • Pryce J.E.
      • Xu Z.Z.
      • Montgomerie W.A.
      Fertility breeding values in New Zealand, the next generation.
      . The multiple-trait animal model contains 8 traits: 270-d milk yield, BCS, PM21_1, CR42_1, PM21_2, CR42_2, PM21_3, and CR42_3. It should be noted that no heifer fertility traits are currently included in the national genetic evaluation, and no use is made of the first calving date record of cows.
      We calculated a pedigree index fertility breeding value (fertility-PI) for each cow by compiling 3 generations of sire EBV, based on each animal's pedigree as follows:
      Fertility-PI = SEBV × 0.5 + SDEBV × 0.52 + SDDEBV × 0.53,


      where SEBV is the fertility EBV of the sire, SDEBV is the fertility EBV of the sire of the dam (i.e., maternal grandsire), and SDDEBV is the fertility EBV of the sire of the dam of the dam. Using earlier generations to calculate an animal's pedigree index breeding value has minimal effect on its accuracy (
      • Wiggans G.R.
      • Powell R.L.
      Increasing pedigree contribution to dairy sire evaluation.
      ). We used a pedigree index instead of a cow's own EBV to avoid the confounding effect of the cow's own phenotype for fertility traits and the fertility EBV against which her records were compared. Although the cow's own individual fertility records contributed to her sire's EBV for fertility, this effect was substantially diluted by the fertility records of her siblings because most sires had very large numbers of daughters. Pedigree index averages across different herd phenotypic classifications for both fertility and milk production were tested using a simple single-factor ANOVA to determine if there were any consistent patterns, particularly to determine whether there was a tendency for cows in low-fertility herds to have genetically lower fertility genetic backgrounds.

      Linear Regression Analysis

      We used linear regression to estimate the influence of fertility-PI on fertility phenotypes within fertility decile and milk production groups. These partitioning variables will be referred to as “environmental variables” hereafter. We separately fit each fertility phenotype—CR42_0 to CR42_3, CSD_0 to CSD_3, and PM21_1 to PM21_3—as the response variable with the herd-level environmental variable of interest and contemporary group (herd-year) as fixed effects. These models also included a continuous effect of fertility-PI fitted and its interaction with the categorical environmental variable. We tested for significance effects of the environmental variables and interaction between the environmental variables and the fertility-PI using the “lm” function in R, followed by an ANOVA procedure and Tukey multiple comparisons test to determine significance between group means (
      • R Core Team
      R: A language and environment for statistical computing.
      ).
      The expectation of the regression slope was a value of 1 for CR42_1, because the phenotype in this case had an identical definition to that used in the national fertility evaluation. We also added a quadratic term to the regression models to test for a nonlinear relationship between fertility-PI and the cow's fertility phenotype. This analysis was applied to the data set in which herds were partitioned into high-, mid-, and low-fertility herds, the same data set compiled for parameter estimation.

      RESULTS

      Overall Data

      The phenotypic means for all fertility-related traits with herds partitioned into high-, mid-, and low-fertility groups are presented in Table 2. The value of CSD_0 was 3 d greater and CR42_0 was l4.8 percentage points lower in low-fertility herds than in high-fertility herds. In later parities, the divergence was much greater. The difference in CSD was 16, 14, and 10 d greater, from PSC, in low-fertility herds for CSD_1, CSD_2, and CSD_3, respectively. Similarly, for CR42, the difference in percentage calved at 42 d was 24, 20, and 16 d lower in low-fertility herds compared with the respective traits in high-fertility herds. The difference in fertility-PI was only 0.7% between high- and low-fertility herd categories.
      Table 2Phenotypic means (phenotypic SD in parentheses) for the fertility traits CSD (calving season day), CR42 (calved within 42 d from planned start of calving), and PM21 (mated within 21 d of the planned start of mating)
      Trait
      CSD is presented in days and CR42 and PM21 as percentages. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      Environment
      High fertilityMid fertilityLow fertility
      CSD_020.03 (16.38)21.34 (17.78)23.17 (20.02)
      CSD_122.61 (18.51)29.20 (24.26)38.84 (31.00)
      CSD_223.81 (18.78)29.27 (23.55)36.98 (28.68)
      CSD_323.79 (18.67)27.80 (22.63)34.07 (27.37)
      CR42_090.69 (29.06)88.87 (31.44)86.12 (34.58)
      CR42_185.56 (35.14)75.24 (43.16)61.93 (48.55)
      CR42_284.12 (36.55)75.64 (42.92)64.10 (47.97)
      CR42_384.45 (36.24)78.04 (41.40)68.46 (46.47)
      PM21_188.26 (32.20)81.77 (38.61)68.91 (46.28)
      PM21_287.59 (32.97)80.71 (39.46)66.52 (47.19)
      PM21_388.57 (31.82)83.46 (37.15)71.41 (45.18)
      1 CSD is presented in days and CR42 and PM21 as percentages. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      Figure 1 illustrates the spread of CSD_1 for cows in each fertility environment. The distributions of CSD_1 in the mid- and low-fertility environments have an extended right tail, and the low-fertility group also has 2 shoulders, which reflects estrous cycles. The shoulders appear at approximately 21-d intervals and, after a 21-d mating cycle, a significant proportion of exposed cows will be already pregnant, so fewer unmated cows are available in subsequent cycles. The transition points are not discrete because of variation in gestation length and calving dates, and because some cows only commence their first postpartum estrous cycle following the first 21 d of mating, but in decreasing numbers over time. Similarly, the density curve for the mid-fertility group has an extended tail with a slight shoulder starting 17 to 20 d from the peak in the density (data not shown).
      Figure thumbnail gr1
      Figure 1Density plot of second-parity calving trait calving season day (CSD1) for cows categorized as being in high (red), mid (green), and low (blue) fertility environments.
      The distribution of the cow's breeding values for high-, mid-, and low-fertility environments is presented in Figure 2. This illustrates that the fertility EBV distribution, across environments, is similar, although cows from low-fertility herds have a slightly extended tail toward unfavorable fertility EBV.
      Figure thumbnail gr2
      Figure 2Distribution of cow fertility EBV for cows in high (red), mid (green), and low (blue) fertility environments.

      Parameter Estimates for High- and Low-Fertility Herd Categories

      Heritabilities for and genetic correlations among the fertility-related traits with data partitioned into low- and high-fertility environments are presented in Table 3. These genetic parameter estimates were generated from the pairwise bivariate genetic analyses with the same trait measured in divergent environments treated as separate traits (i.e., CSD_1_high with CSD_1_low, and so on). These correlations were high, and, in most cases, declined from parity 1 through parity 3 with the exception of CSD, where the correlation increased from 0.68 for first calving (CSD_0) to 0.88 in the final (fourth) calving (CSD_3). Moreover, CR42_3 is the only trait for which we found a negative correlation between low- and high-fertility environments, and this estimate had a very high standard error. The standard errors for genetic correlation estimates were all high, and tended to increase for fertility traits expressed at later ages, when the number of cows with records had declined.
      Table 3Heritabilities (SE in parentheses) in high- and low-fertility environments (with genetic correlations between environments estimated from pairwise bivariate analyses) for the fertility traits PM21 (mated within 21 d of the planned start of mating), CR42 (calved within 42 d from planned start of calving), and CSD (calving season day)
      Trait
      CSD is presented in days and CR42 and PM21 as percentages. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      HeritabilityGenetic correlation
      Low fertility herdsHigh fertility herds
      PM21_10.038 (0.010)0.023 (0.007)0.704 (0.181)
      PM21_20.034 (0.010)0.020 (0.007)0.803 (0.196)
      PM21_30.050 (0.015)0.007 (0.006)0.260 (0.423)
      CR42_00.013 (0.006)0.008 (0.004)0.568 (0.341)
      CR42_10.017 (0.007)0.004 (0.004)0.853 (0.405)
      CR42_20.028 (0.010)0.008 (0.005)0.165 (0.347)
      CR42_30.009 (0.007)0.005 (0.005)−0.136 (0.649)
      CSD_00.012 (0.006)0.012 (0.005)0.679 (0.293)
      CSD_10.025 (0.008)0.016 (0.006)0.630 (0.227)
      CSD_20.021 (0.008)0.018 (0.007)0.572 (0.243)
      CSD_30.021 (0.009)0.014 (0.007)0.888 (0.276)
      1 CSD is presented in days and CR42 and PM21 as percentages. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      Genetic standard deviations between high- and low-fertility herds (results not shown), for all traits, were much greater in low-fertility herds, with the ratio between them ranging between 1.2 and 3.6. Similarly, we obtained higher genetic standard deviation estimates in low-fertility herds using trivariate models, which should better account for preselection issues. Table 4 presents genetic correlations derived from the trivariate analysis. This illustrates the variation in correlations observed between the trivariate and bivariate results, especially in the calving-based fertility traits. Correlations between high- and low-fertility herds for PM21 were more consistent between bivariate and trivariate analyses than those observed for calving date–based traits; however, the heritability estimates for PM21 in the trivariate analysis traits in low-fertility environments were unrealistically high.
      Table 4Heritabilities (diagonals) and genetic correlations (off diagonals) for the fertility traits PM21 (mated within 21 d of the planned start of mating), CR42 (calved within 42 d from planned start of calving), and CSD (calving season day)
      Estimated from trivariate analysis where data for low versus high herd fertility environments were analyzed separately. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      TraitPM21_1_lowPM21_2_lowPM21_3_lowPM21_1_highPM21_2_highPM21_3_high
      PM21_1_low0.089 (0.018)0.441 (0.120)0.271 (0.150)
      PM21_2_low0.135 (0.022)0.768 (0.080)
      PM21_3_low0.135 (0.026)
      PM21_1_high0.142 (0.272)0.696 (0.100)0.567 (0.15)0
      PM21_2_high0.060 (0.013)0.830 (0.100)
      PM21_3_high0.004 (0.014)
      CR42_0_lowCR42_1_lowCR42_2_lowCR42_0_highCR42_1_highCR42_2_high
      CR42_0_low0.012 (0.006)0.108 (0.310)0.069 (0.290)
      CR42_1_low0.017 (0.007)0.780 (0.180)
      CR42_2_low0.029 (0.008)
      CR42_0_high0.006 (0.004)0.448 (0.690)0.359 (0.460)
      CR42_1_high0.003 (0.003)0.939 (0.680)
      CR42_2_high0.008 (0.005)
      CSD_0_lowCSD_1_lowCSD_2_lowCSD_0_highCSD_1_highCSD_2_high
      CSD_0_low0.012 (0.006)0.043 (0.310)0.044 (0.310)
      CSD_1_low0.023 (0.008)0.634 (0.200)
      CSD_2_low0.020 (0.008)
      CSD_0_high0.010 (0.005)0.213 (0.300)0.369 (0.290)
      CSD_1_high0.016 (0.006)0.927 (0.160)
      CSD_2_high0.020 (0.007)
      1 Estimated from trivariate analysis where data for low versus high herd fertility environments were analyzed separately. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      Heritabilities were also provided for binomial fertility traits after transformation to the underlying scale (Table 5). On this scale, there was no evidence of higher heritability in low-fertility herds. This implies that the higher genetic standard deviations for low-fertility herds are due to the greater phenotypic variation for fertility-related traits, which have lower means (Table 2), rather than a proportionally stronger genetic signal.
      Table 5Heritabilities for binomial fertility traits PM21 (percentage mated within 21 d of the planned start of mating) and CR42 (percentage calved within 42 d from planned start of calving) when converted to the underlying scale
      Trait
      The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      Low fertility herdsHigh fertility herds
      PM21_10.0600.066
      PM21_20.0500.058
      PM21_30.0180.091
      CR42_00.0320.032
      CR42_10.0100.028
      CR42_20.0180.047
      CR42_30.0110.016
      1 The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.

      Effect of Sire Fertility Genetic Merit on Cow Phenotype

      We partitioned the data into deciles to investigate changes in regression coefficients between fertility phenotypes and fertility-PI across the range of fertility environments. The decile means for herd average CR42_1 were 85.2, 81.3, 79.1, 77.1, 76.3, 75.1, 73.0, 71.6, 68.6, and 63.0% for deciles 1 to 10, where decile 1 is high herd fertility and decile 10 is low herd fertility.
      We present the regression coefficients for the binary fertility traits CR42 and PM21 on fertility-PI for each herd fertility decile group in Table 6. The main effect of decile group and the interaction between fertility-PI and decile group were statistically significant (P < 0.001) for all traits except heifer calving season day (CSD_0). The regression coefficients indicated that the expression of sire genetic merit for fertility increased in higher decile groups, which have lower herd average fertility. The expression of sire genetic merit for fertility in deciles 9 and 10 was 2.3 and 2.2 times greater than that of decile 1 for CR42. Coefficients also tended to decrease in successive parities, probably because cows with low genetic merit for fertility were culled from the herd in early parities. This culling bias was lower in low-fertility environments, as evidenced by the modest decrease in the coefficient between second and fourth parity for decile groups 1 to 8 (average 0.16), and the more substantial decreases in deciles 9 and 10 (0.45 and 0.24). Calving season day followed a similar trend, when considering that a negative sign for CSD is beneficial (a decreasing CSD means the cow is calving earlier). The average coefficient of deciles 8 to 10 (i.e., lowest performing 30% of herds based on herd average fertility) for CSD was −1.04, compared with an average of −0.61 across deciles 1 to 3 (i.e., top 30%).
      Table 6Regression coefficients (P-values in parentheses) for phenotypic values of binary fertility traits
      CR42 = percentage calved within 42 d from planned start of calving; PM21 = percentage mated within 21 d of the planned start of mating. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      CR42 and PM21, when regressed against a pedigree index of cow breeding values for fertility when herds are partitioned into decile groupings according to their mean phenotypic level of fertility
      Decile group
      Decile 1 contains the highest fertility herds.
      CR42_0CR42_1CR42_2CR42_3PM21_1PM21_2PM21_3
      10.180 (0.048)0.719 (0.059)0.767 (0.068)0.581 (0.076)0.856 (0.053)0.952 (0.059)0.645 (0.064)
      20.127 (0.050)0.966 (0.064)1.048 (0.072)0.783 (0.078)0.978 (0.055)1.104 (0.060)0.744 (0.063)
      30.223 (0.050)1.108 (0.066)1.085 (0.074)0.998 (0.079)1.063 (0.057)1.194 (0.063)1.005 (0.065)
      40.148 (0.051)1.210 (0.069)1.116 (0.076)1.030 (0.081)1.100 (0.059)1.240 (0.066)0.952 (0.069)
      50.219 (0.050)1.222 (0.067)1.189 (0.074)1.011 (0.079)1.193 (0.058)1.162 (0.063)1.116 (0.067)
      60.165 (0.050)1.249 (0.069)1.193 (0.076)1.092 (0.081)1.092 (0.059)1.377 (0.066)1.095 (0.069)
      70.118 (0.048)1.357 (0.068)1.332 (0.074)1.222 (0.079)1.231 (0.059)1.414 (0.066)1.221 (0.069)
      80.221 (0.049)1.326 (0.069)1.314 (0.077)1.183 (0.082)1.363 (0.060)1.318 (0.067)1.151 (0.071)
      90.222 (0.050)1.667 (0.071)1.372 (0.078)1.220 (0.084)1.472 (0.063)1.603 (0.071)1.490 (0.074)
      100.313 (0.049)1.593 (0.069)1.725 (0.078)1.350 (0.084)1.398 (0.062)1.611 (0.072)1.451 (0.077)
      1 CR42 = percentage calved within 42 d from planned start of calving; PM21 = percentage mated within 21 d of the planned start of mating. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      2 Decile 1 contains the highest fertility herds.
      The regressions of fertility-PI on heifer calving traits CSD_0 and CR42_0 had a very shallow slope, indicating that the fertility-PI is only a weak predictor of heifer fertility. Neither first calving dates nor any other indicators of heifer fertility are currently included in the genetic evaluation of fertility. Consequently, any predictive ability of a sire's EBV for heifer reproductive performance relies entirely on the genetic correlations between mating and calving traits at later parities and the fertility of heifers at first mating.

      Nonlinear Model

      A quadratic term for fertility-PI was found to be statistically significant (P < 0.05) when added to the regression model for the second-parity traits CSD_1 and CR42_1 and highly significant (P < 0.001) for CSD and CR42 for all later parities and the 3 PM21 traits. The quadratic terms for heifer calving traits were not significant. Figure 3 illustrates that the regression of fertility-PI on CR42_1 had a diminishing slope in high- and mid-fertility herds. In contrast, low-fertility herds showed an increasing slope. Similar patterns were observed for CSD and PM21 (results not shown), with the pattern becoming less consistent for all traits as the cows moved through successive parities.
      Figure thumbnail gr3
      Figure 3Nonlinear regression of sire EBV on fertility phenotype for percentage calved in the first 42 d from the planned start of calving for second parity (CR42_1).

      Production Partitioning

      The means of fertility traits and regression coefficients for fertility traits on fertility-PI for each herd production level are presented in Table 7. Unsurprisingly, daily milk solids yield in peak season for second-parity cows was 15.9% higher in high-production herds than in “other” production herds. The differences in mean fertility phenotypes indicated that animals in high-production herds had slightly higher fertility than those in less productive herds. The interaction between fertility-PI and production level was significant (P < 0.05) only for CSD_2, CR42_2, PM21_2, and PM21_3. The greatest level of significance (P < 0.001) was observed for PM21_2, with a difference between slopes of 0.27.
      Table 7Regression coefficients (Reg. Coeff.) for fertility traits CSD, CR42, and PM21
      The cow's pedigree breeding value is regressed against the cow's fertility phenotype, when animals were partitioned into high production and other production level groups based on second-lactation daily milk solids production after accounting for breed composition of the herd.
      Trait
      CSD = calving season day (in days); CR42 = percentage calved within 42 d from planned start of calving; PM21 = percentage mated within 21 d of the planned start of mating; Fertility-PI = pedigree index fertility breeding value. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      High production herdOther production herds
      MeanReg. Coeff.MeanReg. Coeff.
      CSD_021.43 (0.05)−0.180 (0.021)21.49 (0.05)−0.168 (0.013)
      CSD_128.37 (0.07)−0.783 (0.028)30.12 (0.04)−0.861 (0.018)
      CSD_2
      P < 0.01,
      28.11 (0.04)−0.745 (0.029)30.11 (0.04)−0.877 (0.019)
      CSD_327.44 (0.05)−0.659 (0.032)28.44 (0.05)−0.732 (0.020)
      CR42_088.61 (0.08)0.232 (0.037)88.74 (0.05)0.189 (0.023)
      CR42_177.31 (0.11)1.200 (0.048)73.77 (0.07)1.296 (0.031)
      CR42_2
      P < 0.05,
      77.73 (0.12)1.163 (0.053)74.32 (0.08)1.310 (0.035)
      CR42_378.93 (0.13)1.008 (0.058)76.97 (0.08)1.110 (0.037)
      PM21_183.83 (0.10)1.071 (0.041)80.27 (0.07)1.196 (0.027)
      PM21_2
      P < 0.001: significance of interaction between fertility-PI and production level.
      82.69 (0.11)1.100 (0.046)79.13 (0.07)1.379 (0.030)
      PM21_3
      P < 0.05,
      84.43 (0.11)0.914 (0.048)82.14 (0.07)1.184 (0.032)
      Fertility-PI0.265 (0.0067) NA
      Not applicable.
      0.743 (0.004)NA
      Herd production2.129 (0.0002)NA1.790 (0.0003)NA
      Total animals375,488NA143,959NA
      Total number of herds2,657NA636NA
      1 The cow's pedigree breeding value is regressed against the cow's fertility phenotype, when animals were partitioned into high production and other production level groups based on second-lactation daily milk solids production after accounting for breed composition of the herd.
      2 CSD = calving season day (in days); CR42 = percentage calved within 42 d from planned start of calving; PM21 = percentage mated within 21 d of the planned start of mating; Fertility-PI = pedigree index fertility breeding value. The parity that each trait relates to is denoted by 0 to 3, where 0 is the first calving and 3 is the fourth calving.
      3 Not applicable.
      * P < 0.05,
      ** P < 0.01,
      *** P < 0.001: significance of interaction between fertility-PI and production level.

      DISCUSSION

      Identification and a deeper understanding of the impact of a G×E interaction can enable improved decision making to optimize sire selection, refine breeding objectives, and design improved genetic evaluation systems. In most cases, we found genetic correlations >0.61, suggesting that a separate trait definition or breeding program for high- and low-fertility environments is not warranted (
      • Mulder H.A.
      • Veerkamp R.F.
      • Ducro B.J.
      • van Arendonk J.A.M.
      • Bijma P.
      Optimization of dairy cattle breeding programs for different environments with genotype by environment interaction.
      ). Several studies have estimated genetic correlations that indicate differences, in terms of sire reranking, between environment (
      • Windig J.J.
      • Calus M.P.L.
      • Beerda B.
      • Veerkamp R.F.
      Genetic correlations between milk production and health and fertility depending on herd environment.
      ;
      • Haile-Mariam M.
      • Carrick M.J.
      • Goddard M.E.
      Genotype by environment interaction for fertility, survival, and milk production traits in Australian dairy cattle.
      ), and others have found no effect of G×E (
      • Kearney J.F.
      • Schutz M.M.
      • Boettcher P.J.
      Genotype x environment interaction for grazing vs. confinement. II. Health and reproduction traits.
      ;
      • van der Laak M.
      • van Pelt M.
      • de Jong G.
      • Mulder H.
      Genotype by environment interaction for production, somatic cell score, workability, and conformation traits in Dutch Holstein-Friesian cows between farms with or without grazing.
      ).
      In our study, we found no conclusive evidence to indicate G×E for fertility traits in terms of sire reranking for individual fertility traits. It is possible that the inconsistent results we found from bivariate and trivariate analyses, and across similar trait definitions, stem from the scale of the data used, because fertility traits are of low heritability and the number of records required to generate accurate genetic parameters would have made computation difficult to complete using an animal model for traits that are not normally distributed, as applied in this study. This highlights the difficulty of trying to prove the existence of a sire reranking G×E for low heritability traits.
      We found that low-fertility herds had a 1.6 to 3.2 times greater genetic standard deviation across all fertility traits compared with high-fertility herds. This suggests that response to selection for fertility should be greater in low- than in high-fertility herds, all else being equal. Given unreliable estimates from other aspects of this analysis, these results should be viewed with caution. However, findings are similar to recent research investigating G×E for livability (survival) of dairy calves from primiparous cows (
      • Ouweltjes W.
      • Windig J.J.
      • van Pelt M.L.
      • Calus M.P.L.
      Genotype by environment interaction for livability of dairy calves from first parity cows.
      ). Those researchers found across-environment genetic correlations close to 1, indicating no G×E. However, they did find a consistent trend of increasing heritability with decreasing herd levels of livability, similar to our results, and concluded that sire selection for improved livability would be twice as effective in low-livability herds as in high-livability herds. Although a high genetic correlation indicates an absence of G×E at the single-trait level, scaling effects for individual traits can still lead to a crossing-type G×E when aggregating merit across multiple economically important traits.
      Our results indicate that the expected response to selection on genetic merit for fertility was higher in low-fertility herds, with the response for the lowest 30% of herds being more than twice that of the same relationship in the top 30% of herds. This pattern was strongest when we included a quadratic term in analyses of the high, mid, and low partitioned data, as illustrated by CR42_1 in Figure 3, which shows a clear increasing response to fertility-PI in low-fertility herds. In mid-fertility herds, the response was more linear, and in high-fertility herds, a leveling-off in response was observed once fertility-PI moved past the average EBV for fertility. This indicates that high-fertility herds would gain less from using bulls with high fertility EBV. In contrast, low-fertility herds would experience the opposite—a substantially higher proportional response to using high-fertility bulls.
      Heifer calving traits are currently not included in generating fertility EBV, making it important to note that, in all cases, the heifer calving traits CSD_0 and CR42_0 had consistently much lower regression coefficients on fertility-PI. This implies that fertility-PI is not a good predictor of heifer calving, although the slope was typically in the favorable direction. However,
      • Bowley F.E.
      • Green R.E.
      • Amer P.R.
      • Meier S.
      Novel approaches to genetic analysis of fertility traits in New Zealand dairy cattle.
      , using similar fertility data from the New Zealand dairy industry, estimated a genetic correlation between heifer CSD and later-parity CSD of 0.7. Thus, the modest regression slopes found for heifer CSD on fertility-PI may reflect the fact that heifers are more fertile than lactating cows, and there is therefore less phenotypic variation in first-parity than in later-parity fertility traits. This reduced phenotypic variation would partly explain the lower slope observed for CSD_0 on fertility-PI.
      Production partitioning revealed that high-production herds had slightly higher fertility. These results are inconsistent with expectations from intensive feeding systems with much higher mean milk production per cow (
      • Pryce J.E.
      • Woolaston R.
      • Berry D.P.
      • Wall E.
      • Winters M.
      • Butler R.
      • Shaffer M.
      World trends in dairy cow fertility.
      ), which would be expected to have lower fertility. In New Zealand's more pastoral and seasonal dairy production systems, higher fertility in high-production herds may result from better management or increased management contact improving heat detection. We found a significant interaction between fertility-PI and production group only for 4 of the 11 fertility-related traits. Higher regression slopes in low-production herds for some fertility-related traits may simply reflect the lower mean fertility level in these herds, indicating that mean herd fertility level is probably the driver of differences in the slopes, rather than production level per se.
      Our genetic parameter estimates and regression analyses illustrate that low-fertility herds have greater variation for fertility traits. The nature of this means that greater gains can be made in low-fertility herds than in high-fertility herds. Results of parameter estimation indicate a greater genetic standard deviation in low-fertility herds due largely to greater phenotypic variation with the same heritability when heritability was transformed to the trait on the underlying scale.
      Expression of sire genetic merit for fertility was magnified in herds with low fertility performance. Therefore, herds with a fertility problem would benefit more from utilizing sires with higher fertility EBV than would high-fertility herds. Understanding the performance difference of sire genetics between high- and low-fertility herds can provide industry with clear, quantifiable messages that can be sent to farmers, supporting improved sire selection to increase the fertility performance of their herds (e.g., encourage farmers with fertility issues in their herds to favor bulls with high Breeding Worth that also have high fertility EBV to rectify long-term fertility problems using genetics). Alternatively, a fertility-specific index could be developed with more emphasis on fertility, for use by farmers operating a herd with low fertility performance. However, the results of this study do not provide sufficient evidence to justify a separate trait evaluation for fertility either in high- versus low-production herds or in high- versus low-fertility herds.

      ACKNOWLEDGMENTS

      This work was funded by dairy farmers through DairyNZ (Hamilton, New Zealand; RD1402, RD1404), by the Ministry for Primary Industries as part of the Transforming the Dairy Value Chain Primary Growth Partnership and the Ministry of Business, Innovation and Employment (DRCX1302).

      REFERENCES

        • Amer P.R.
        • Stachowicz K.
        • Jenkins G.M.
        • Meier S.
        Short communication: Estimates of genetic parameters for dairy fertility in New Zealand.
        J. Dairy Sci. 2016; 99 (27448853): 8227-8230
        • Beam S.W.
        • Butler W.R.
        Effects of energy balance on follicular development and first ovulation in postpartum dairy cows.
        J. Reprod. Fertil. Suppl. 1999; 54 (10692872): 411-424
        • Berry D.P.
        • Buckley F.
        • Dillon P.
        • Evans R.D.
        • Rath M.
        • Veerkamp R.F.
        Genetic relationships among body condition score, body weight, milk yield, and fertility in dairy cows.
        J. Dairy Sci. 2003; 86 (12836956): 2193-2204
        • Bowley F.E.
        • Green R.E.
        • Amer P.R.
        • Meier S.
        Novel approaches to genetic analysis of fertility traits in New Zealand dairy cattle.
        J. Dairy Sci. 2015; 98 (25597971): 2005-2012
        • DairyNZ-Livestock Improvement Corporation
        New Zealand dairy statistics 2014–15.
        DairyNZ and LIC, Hamilton New Zealand2015
        • Dempster E.R.
        • Lerner I.M.
        Heritability of threshold characters.
        Genetics. 1950; 35 (17247344): 212-236
        • Dickerson G.E.
        Inbreeding and heterosis in animals.
        in: Proc. Anim. Breed. Genet. Symp. Honor Dr. Jay L. Lush. American Society of Animal Science, American Dairy Science Association, and Poultry Science Association, Champaign, IL1973: 54-77
        • Falconer D.S.
        • Mackay T.F.C.
        Introduction to Quantitative Genetics.
        4th ed. Addison Wesley Longman, Harlow, Essex, UK1996
        • Gilmour A.R.
        • Cullis B.R.
        • Welham S.J.
        • Thompson R.
        ASREML Reference Manual. NSW Agric. Biometric Bull.
        Department of Primary Industries and Fisheries, Brisbane, Australia1999
        • Grosshans T.
        • Xu Z.Z.
        • Burton L.J.
        • Johnson D.L.
        • Macmillan K.L.
        Performance and genetic parameters for fertility of seasonal dairy cows in New Zealand.
        Livest. Prod. Sci. 1997; 51: 41-51
        • Haile-Mariam M.
        • Carrick M.J.
        • Goddard M.E.
        Genotype by environment interaction for fertility, survival, and milk production traits in Australian dairy cattle.
        J. Dairy Sci. 2008; 91 (19038960): 4840-4853
        • Harris B.L.
        • Montgomerie W.A.
        Fertility breeding values for seasonal dairying.
        Interbull Bull. 2001; 27: 139-142
        • Harris B.L.
        • Pryce J.E.
        Genetic and phenotypic relationships between milk protien percentage, reproductive performance and body condition score in New Zealand dairy cattle.
        Proc. N.Z. Soc. Anim. Prod. 2004; 64: 127-131
        • Harris B.L.
        • Pryce J.E.
        • Xu Z.Z.
        • Montgomerie W.A.
        Fertility breeding values in New Zealand, the next generation.
        Interbull Bull. 2005; 33: 47-50
      1. Harris, B. L., A. M. Winkelman, and L. J. Burton. 2000. Comparisons of fertility measures in strains of Holstein-Friesian cows and their crosses. Massey Dairy Farming Annu. Proc. Massey Dairy Farmers Conference.

        • Holtsmark M.
        • Heringstad B.
        • Madsen P.
        • Ødegård J.
        Genetic relationship between culling, milk production, fertility, and health traits in Norwegian red cows.
        J. Dairy Sci. 2008; 91 (18832226): 4006-4012
        • Kearney J.F.
        • Schutz M.M.
        • Boettcher P.J.
        Genotype x environment interaction for grazing vs. confinement. II. Health and reproduction traits.
        J. Dairy Sci. 2004; 87 (14762094): 510-516
        • Lucy M.C.
        Reproductive loss in high-producing dairy cattle: Where will it end?.
        J. Dairy Sci. 2001; 84 (11417685): 1277-1293
        • Macdonald K.A.
        • Thorrold B.S.
        • Glassey C.B.
        • Holmes C.W.
        • Pryce J.E.
        Impact of farm management decision rules on the production and profit of different strains of Holstein-Friesian dairy cows.
        Proc. N.Z. Soc. Anim. Prod. 2005; 65: 40-45
        • MacMillan K.L.
        Calving patterns and herd production in seasonal dairy herds.
        Proc. N.Z. Soc. Anim. Prod. 1979; 39: 168-174
        • MacMillan K.L.
        • Taufa V.K.
        • Phillips P.
        Recent trends in conception rates and return patterns in AB herds and their effects on calving patterns.
        Proc. N.Z. Soc. Anim. Prod. 1984; 44: 63-65
        • Mulder H.A.
        • Veerkamp R.F.
        • Ducro B.J.
        • van Arendonk J.A.M.
        • Bijma P.
        Optimization of dairy cattle breeding programs for different environments with genotype by environment interaction.
        J. Dairy Sci. 2006; 89 (16606745): 1740-1752
        • Ouweltjes W.
        • Windig J.J.
        • van Pelt M.L.
        • Calus M.P.L.
        Genotype by environment interaction for livability of dairy calves from first parity cows.
        Animal. 2015; 9 (26123138): 1617-1623
        • Philipsson J.
        • Lindhé B.
        Experiences of including reproduction and health traits in Scandinavian dairy cattle breeding programmes.
        Livest. Prod. Sci. 2003; 83: 99-112
        • Pritchard T.
        • Coffey M.
        • Mrode R.
        • Wall E.
        Understanding the genetics of survival in dairy cows.
        J. Dairy Sci. 2013; 96 (23477814): 3296-3309
        • Pryce J.E.
        • Woolaston R.
        • Berry D.P.
        • Wall E.
        • Winters M.
        • Butler R.
        • Shaffer M.
        World trends in dairy cow fertility.
        in: Proc. World Congr. Genet. Appl. Livest. Prod. 10. Vancouver, Canada. American Society of Animal Science, Champaign, IL2014: 1-10
        • R Core Team
        R: A language and environment for statistical computing.
        R Foundation for Staistical Computing, Vienna, Austria2015
        • van Arendonk J.A.M.
        Use of profit equations to determine relative economic value of dairy cattle herd life and production from field data.
        J. Dairy Sci. 1991; 74: 1101-1107
        • van der Laak M.
        • van Pelt M.
        • de Jong G.
        • Mulder H.
        Genotype by environment interaction for production, somatic cell score, workability, and conformation traits in Dutch Holstein-Friesian cows between farms with or without grazing.
        J. Dairy Sci. 2016; 99 (27040792): 4496-4503
        • Wiggans G.R.
        • Powell R.L.
        Increasing pedigree contribution to dairy sire evaluation.
        J. Dairy Sci. 1984; 67 (6725733): 893-896
        • Windig J.J.
        • Calus M.P.L.
        • Beerda B.
        • Veerkamp R.F.
        Genetic correlations between milk production and health and fertility depending on herd environment.
        J. Dairy Sci. 2006; 89 (16606748): 1765-1775
        • Windig J.J.
        • Calus M.P.L.
        • Veerkamp R.F.
        Influence of herd environment on health and fertility and their relationship with milk production.
        J. Dairy Sci. 2005; 88 (15591398): 335-347
        • Xu Z.
        • Burton L.J.
        Reproductive Performance of Dairy Cows in New Zealand: Monitoring Fertility.
        in: Proceedings of a Seminar of the Society of Dairy Cattle Veterinarians of New Zealand Veterinary Association. vol. 17. 2003: 23-41