Advertisement

Increased parity is negatively associated with survival and reproduction in different production systems

Open AccessPublished:November 29, 2022DOI:https://doi.org/10.3168/jds.2021-21672

      ABSTRACT

      We conducted a retrospective meta-analysis based on individual cow data to assess the associations of parity, level of production, and pasture-based or intensively fed systems with fertility. Our goal was to provide understandings of the role of parity in risks for removal and reproductive failure. Multilevel models were used to evaluate the fixed effects of parity, milk, milk solids, milk fat and protein percentage and yield, and production system [intensively fed (n = 28,675) or predominantly pasture fed (n = 4,108)] on reproductive outcomes. The outcomes were the hazard of not being bred (HNBRED), hazard of pregnancy (HPREG), pregnancy to first breeding (PREG1), and odds of becoming pregnant in a lactation (OPAL). The 32,783 cows were in 13 studies conducted in Australia (14.6% of cows), Canada (2.4% of cows), and the United States (83.0% of cows). There were 38.5% of cows in the sample in parity 1, 27.3% in parity 2, 16.7% in parity 3, 9.0% in parity 4, and 8.6% in parity ≥5. Compared with cows of parity 1, parity ≥5 cows had a greater HNBRED [hazard ratio (HR) = 2.45], lesser HPREG (HR = 0.73), and reduced OPAL (odds ratio = 0.36). However, the parity ≥5 cows had similar PREG1 to other parities except for parity 1. This suggests the possibility of a higher proportion of subfertile parity ≥5 cows than for other parities. Associations between parity and reproductive measures were influenced by the different milk production measures, indicating that milk yield and milk component percentages and yields modified the odds or hazards of successful reproduction. All milk production measures had quadratic associations with OPAL, indicating that either low or high production or concentration of solids within a cohort reduced OPAL. This reduced OPAL reflected a greater HNBRED for lower milk yield and milk protein and fat yielding cows. Both milk yield and milk protein percentage had quadratic associations with HPREG. When centered milk yield was categorized into quartiles, small differences in HPREG existed. A more marked association of milk protein percentage occurred with HPREG, with optimal HPREG at approximately 0.5% above group mean milk protein percentage. Milk fat percentage (HR = 0.901), fat yield (kg/d; HR = 0.78), protein yield (kg/d; HR = 0.71), and milk solids yield (kg/d; HR = 0.84) were all linearly associated with reduced HPREG. Difference in production systems did not have substantive effects on PREG1 but did for HNBRED, HPREG, and OPAL. Estimates of associations of parity with reproductive outcomes HNBRED, HPREG, and OPAL were influenced by milk and milk solids yield; older cows had markedly lower reproductive outcomes. Interestingly, for PREG1, there were few differences among parities and differences were less influenced by milk yield and constituent measures. The marked associations of parity with removal for all reasons, deaths and culling, and reductions in HNBRED, HPREG, and OPAL indicate a need to focus on the physiological changes with parity to produce better strategies to support optimal longevity of cows.

      Graphical Abstract

      Key words

      INTRODUCTION

      It is necessary to ensure that cows economically optimize their genetic potential by providing environments that permit them to be productive, healthy, and reproductively successful.
      • Dallago G.M.
      • Wade K.M.
      • Cue R.I.
      • McClure J.T.
      • Lacroix R.
      • Pellerin D.
      • Vasseur E.
      Keeping dairy cows for longer: A critical literature review on dairy cow longevity in high milk-producing countries.
      noted that the longevity of dairy cattle has decreased in most high milk-producing countries over time and that increasing longevity would have benefits including increased profitability.
      • De Vries A.
      • Marcondes M.
      Overview of factors affecting productive lifespan of dairy cows.
      suggested that the economic trade-offs between longevity and profit may be complex given that rapid genetic gain is possible. A study conducted in the United States identified that more than 40% of cows in 3 herds were parity 1, indicating a substantial loss of older cows (
      • Golder H.M.
      • Rossow H.
      • Lean I.
      Effects of in-feed enzymes on milk production and components, reproduction, and health in dairy cows.
      ). Percentage of cows culled in Canada from 2015 to 2020 was 30.7% () and the Northeastern United States in 2014 was 31.4% plus a 6.2% death rate (
      • United States Department of Agriculture (USDA)
      Health and management practices on US dairy operations 2014, Report 3.
      ). The cows of parity ≥ 2 in a United States-based study (
      • Golder H.M.
      • Rossow H.
      • Lean I.
      Effects of in-feed enzymes on milk production and components, reproduction, and health in dairy cows.
      ) and in large Australian studies had at least 2 times the risk of many diseases and had reduced odds and hazard of pregnancy (HPREG) than parity 1 cows (
      • Golder H.M.
      • McGrath J.
      • Lean I.J.
      Effect of 25-hydroxyvitamin D3 during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows.
      ). Given increased use of sexed semen (
      • Holden S.A.
      • Butler S.T.
      Review: Applications and benefits of sexed semen in dairy and beef herds.
      ), which can increase the supply of replacement animals, farmers may view older cows as a liability because of increased risk of reproductive and health disorders or may intend to accelerate genetic progress by bringing young animals into the milking herd (
      • Fetrow J.
      • Nordlund K.
      • Norman H.
      Invited review: Culling: Nomenclature, definitions, and recommendations.
      ). However,
      • De Vries A.
      • Marcondes M.
      Overview of factors affecting productive lifespan of dairy cows.
      suggest that an average productive life span—that is, from time of calving to removal (
      • Schuster J.C.
      • Barkema H.W.
      • De Vries A.
      • Kelton D.F.
      • Orsel K.
      Invited review: Academic and applied approach to evaluating longevity in dairy cows.
      )—of approximately 5 yr may increase profitability and improve social acceptability of dairy production.
      A better understanding of reasons for removal and how these vary with parity will optimize longevity. Our goal was to provide understandings of the role of parity in risks for removal and reproductive failure with the intent of examining hypotheses based on these findings in future studies. We aim to quantify associations of parity with survival and reproductive performance, statistically controlling for the level of milk outputs, in cows in intensively fed herds and in predominantly pasture-fed herds using data from studies conducted in the United States, Canada, and Australia.
      • Bascom S.S.
      • Young A.
      A summary of the reasons why farmers cull cows.
      found that reproduction was the most commonly ascribed reason for culling in a sample of herds in the New England region of the United States. Similarly,
      • Stevenson M.A.
      • Lean I.
      Descriptive epidemiological study on culling and deaths in eight dairy herds.
      and
      • Workie Z.W.
      • Gibson J.P.
      • van der Werf J.H.J.
      Analysis of culling reasons and age at culling in Australian dairy cattle.
      found that reproductive failure was the most important cause of culling in Australian herds. It is important to note that many culling decisions are based on several factors rather than a single reason for removal (
      • Stevenson M.A.
      • Lean I.
      Descriptive epidemiological study on culling and deaths in eight dairy herds.
      ).
      Associations between increased production and reduced pregnancy per AI or extended days open were described in the 1970s by
      • Spalding R.
      • Everett R.
      • Foote R.
      Fertility in New York artificially inseminated Holstein herds in dairy herd improvement.
      and more recently by many workers (
      • Lucy M.C.
      Reproductive loss in high-producing dairy cattle: Where will it end?.
      ;
      • Berry D.P.
      • Friggens N.C.
      • Lucy M.
      • Roche J.
      Milk production and fertility in cattle.
      ). Flaws in many of the studies evaluating these associations have been detailed and include use of inappropriate measures of reproduction or production (
      • LeBlanc S.
      Assessing the association of the level of milk production with reproductive performance in dairy cattle.
      ). Measures based on total lactational milk or milk solids production contain interactions between reproductive success and subsequent milk production. Briefly, lower-producing cows may not be bred or breeding will cease sooner, and nonpregnant cows may be culled before the end of standard lactations (e.g., 305 d). Consequently, milk production measures made before breeding are less influenced by reproductive success, so are more valid. Further studies are needed on reproductive measures that are less influenced by loss of cows due to selective culling and mating decisions, such as pregnancy to first mating or time to pregnancy.
      This retrospective study is a meta-analysis based on individual cow data obtained from previous studies. We examined survival of cows and reproductive performance for different parities and hypothesized that parity, production system (categorized as pasture-based or intensively fed with a mixed or component-fed ration), milk, milk fat or protein percentage, and milk solids production would influence the probability of survival, the hazard of not being bred (HNBRED), pregnancy to first breeding (PREG1), or odds of becoming pregnant in a lactation (OPAL). The focus was to explore aspects of the patho-physiological attributes of cows of different parities with the intent of obtaining a greater understanding of factors leading to loss of older cows from herds and interactions between milk, milk fat and protein, and fertility. To account for differences in yield over time and between systems, we assess milk production and milk solids for cows within a cohort group.

      MATERIALS AND METHODS

      The studies utilized to produce the database for this study were approved by the relevant Animal Care and Use Committees at the time of conduct.

      Study Selection Criteria

      The studies represent a convenience sample of data that were collected during the conduct of previous research programs. In general, dairies were purposively selected for use in the original studies on the basis that these had good record-keeping and a history of performance that suggested that they would be capable of maintaining attention to detail consistent with successful trial conduct. The study contributors had conducted large, prospective studies that allowed a rigorous evaluation of the original study hypotheses. We included only data on Holstein or Holstein-Friesian cattle as there were relatively few cows of other breeds in the databases. When an intervention was applied in a herd, treatments within the herd were considered as separate groups for statistical analytical purposes. Cows in observational studies were treated as a group for each herd studied. This provided control for any treatment effects in the studies used that influenced survival, milk, milk fat or protein percentages or yield, or reproduction.
      Nine studies examining survival conducted between 1995 and 2020 formed the first database and ended at variable times. Some had data to >300 DIM, whereas others allowed for observations to d 95 before censoring outcomes (Table 1). Some studies recorded reasons for removal. Some studies recorded reasons for removal and, specifically, culling or death, allowing an evaluation of all removals, time to death, or time to culling. There was a maximum of 30,780 cows for survival and these were in 132 groups. Table 2 provides details on the 13 prospective studies that were used to produce the database for reproduction, including references for those that are published. Not all studies were based on interventions; a brief description of the study purpose is also provided in Table 1, Table 2. Study inclusion criteria for reproduction included being an observational or randomized controlled trial that provided details on parity, milk production, and reproductive measures, including time to first insemination or pregnancy and censoring for those events. Hormonal interventions were applied to the best of our knowledge consistently within herds. Measures of milk fat and protein production were desirable but not essential. We considered that the monthly herd tests provide sufficient precision to test our hypotheses.
      Table 1Summary of the 9 prospective studies used in the survival database including year that the study commenced, country in which it was conducted, type of production system, numbers of cattle included in this database, intervention used (purpose of study), and duration of the observation period
      ReferenceStudy yearCountry
      AU = Australia; USA = United States of America; CA = Canada.
      Production system
      Other = pasture fed; TMR = intensively fed.
      n of cowsInterventionObservation period (d)
      • Duffield T.F.
      • Leslie K.
      • Sandals D.
      • Lissemore K.
      • McBride B.
      • Lumsden J.
      • Dick P.
      • Bagg R.
      Effect of a monensin-controlled release capsule on cow health and reproductive performance.
      1994CATMR996Monensin95
      • Duffield T.
      • Bagg R.
      • DesCoteaux L.
      • Bouchard E.
      • Brodeur M.
      • DuTremblay D.
      • Keefe G.
      • LeBlanc S.
      • Dick P.
      Prepartum monensin for the reduction of energy associated disease in postpartum dairy cows.
      1998CATMR152Monensin120
      • LeBlanc S.J.
      • Duffield T.F.
      • Leslie K.E.
      • Bateman K.G.
      • TenHag J.
      • Walton J.S.
      • Johnson W.H.
      The effect of prepartum injection of vitamin E on health in transition dairy cows.
      1998CATMR975Vitamin120
      • Golder H.M.
      • Rossow H.
      • Lean I.
      Effects of in-feed enzymes on milk production and components, reproduction, and health in dairy cows.
      2016USATMR6,395Enzyme210
      Lean (unpublished)2003AUOther330Micronutrient329
      • Vieira-Neto A.
      • Duarte G.A.
      • Zimpel R.
      • Thatcher W.W.
      • Santos J.E.P.
      Days in the prepartum group are associated with subsequent performance in Holstein cows.
      2012USATMR9,076Observational300
      • Pinedo P.
      • Santos J.E.P.
      • Chebel R.C.
      • Galvão K.N.
      • Schuenemann G.M.
      • Bicalho R.C.
      • Gilbert R.O.
      • Rodriguez Zas S.
      • Seabury C.M.
      • Rosa G.
      • Thatcher W.W.
      Early-lactation diseases and fertility in 2 seasons of calving across US dairy herds.
      2013USATMR11,489Observational305
      • Golder H.M.
      • McGrath J.
      • Lean I.J.
      Effect of 25-hydroxyvitamin D3 during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows.
      , experiment 1
      2017AUOther764Calcidiol305
      Golder and Lean (unpublished)2020AUTMR603Amino acid305
      Total30,780
      1 AU = Australia; USA = United States of America; CA = Canada.
      2 Other = pasture fed; TMR = intensively fed.
      Table 2Summary of the 13 studies used in the reproductive database including year that the study commenced, country in which it was conducted, type of production system, numbers of cattle included in the present analysis, whether the study included milk component data, intervention used (purpose of study), the original period in days of observation for pregnancy in each study (data were right censored at 305 DIM), and publication details
      Study yearCountryProduction system
      TMR = intensively fed.
      n of cowsMilk component data
      Y = yes; N = no.
      InterventionOriginal period of observation for pregnancy (d)Reference
      1992AustraliaPasture237NObservational368
      • Curtis M.A.
      Epidemiology of Uterine Infections in Dairy Cows: Antioxidant and Metabolic Investigations.
      • Curtis M.A.
      • Lean I.J.
      Path analysis of metabolic and antioxidant risk factors for periparturient and postparturient conditions and reproductive performance in dairy cows.
      1994AustraliaTMR81YProtein and genetic merit150
      • Westwood C.T.
      • Lean I.J.
      • Garvin J.K.
      • Wynn P.C.
      Effects of genetic merit and varying dietary protein degradability on lactating dairy cows.
      • Westwood C.T.
      • Lean I.J.
      • Garvin J.K.
      Factors influencing fertility of Holstein dairy cows: A multivariate description.
      1995CanadaTMR796NMonensin490
      • Duffield T.F.
      • Leslie K.
      • Sandals D.
      • Lissemore K.
      • McBride B.
      • Lumsden J.
      • Dick P.
      • Bagg R.
      Effect of a monensin-controlled release capsule on cow health and reproductive performance.
      • Duffield T.F.
      • Leslie K.
      • Sandals D.
      • Lissemore K.
      • McBride B.
      • Lumsden J.
      • Dick P.
      • Bagg R.
      Effect of prepartum administration of monensin in a controlled-release capsule on milk production and milk components in early lactation.
      1999AustraliaPasture765YObservational302
      • Moss N.
      The epidemiology of subfertility in Australian dairy cows.
      2002AustraliaPasture283YProtected fat131Lean (unpublished)
      2003AustraliaPasture336NMicronutrient198Lean (unpublished)
      2005AustraliaPasture610NObservational235
      • DeGaris P.J.
      • Lean I.
      • Rabiee A.
      • Heuer C.
      Effects of increasing days of exposure to prepartum transition diets on reproduction and health in dairy cows.
      2012USATMR9,075NObservational300
      • Vieira-Neto A.
      • Duarte G.A.
      • Zimpel R.
      • Thatcher W.W.
      • Santos J.E.P.
      Days in the prepartum group are associated with subsequent performance in Holstein cows.
      2013USATMR11,728NObservational305
      • Pinedo P.
      • Santos J.E.P.
      • Chebel R.C.
      • Galvão K.N.
      • Schuenemann G.M.
      • Bicalho R.C.
      • Gilbert R.O.
      • Rodriguez Zas S.
      • Seabury C.M.
      • Rosa G.
      • Thatcher W.W.
      Early-lactation diseases and fertility in 2 seasons of calving across US dairy herds.
      2016USATMR6,392YEnzyme350
      • Golder H.M.
      • Rossow H.
      • Lean I.
      Effects of in-feed enzymes on milk production and components, reproduction, and health in dairy cows.
      2017AustraliaPasture764YCalcidiol300Exp. 1,
      • Golder H.M.
      • McGrath J.
      • Lean I.J.
      Effect of 25-hydroxyvitamin D3 during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows.
      2017AustraliaPasture1,113YCalcidiol300Exp. 2,
      • Golder H.M.
      • McGrath J.
      • Lean I.J.
      Effect of 25-hydroxyvitamin D3 during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows.
      2020AustraliaTMR603YAmino acid487Golder and Lean (unpublished)
      1 TMR = intensively fed.
      2 Y = yes; N = no.
      Veterinary diagnosis was used to determine pregnancy, and Table 2 provides detail on the maximum interval to pregnancy within a study. The studies were conducted from 1992 to 2020 in Australia, United States, and Canada and included pasture- or nonpasture-fed (intensively fed) herds using seasonal, split, or all year-round breeding strategies (Table 1, Table 2) on the basis that breeding strategies may be more intensive for seasonal and split than all year-round herds and that diets would also differ between systems. Many of the intensively fed herds used fixed time insemination hormonal synchrony methods for first insemination; if so, treatments were intended to be applied equally between groups, and we have assumed this to be the case. Thirteen studies were retained for the reproduction analyses, utilizing 72 herds and containing 120 separate analysis groups. Treated and control cows were approximately equal in number for those studies, with the exception of
      • Golder H.M.
      • McGrath J.
      • Lean I.J.
      Effect of 25-hydroxyvitamin D3 during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows.
      .

      Animal Data

      The rules for censoring of cows used in the published base studies were adopted. For cattle that did not have a date for removal recorded in the data set (n = 35) or for cows that did not have pregnancy status confirmed at removal, these were treated as not pregnant and a censoring date for those cows was based on the last date for mating or pregnancy diagnosis or last recorded date for the study. For reproductive data for cows lacking a recorded date at removal or insemination date, removal date was estimated from the last recorded information, and these were right censored. We consider that the assumption of nonpregnancy at removal should not be differential and not substantially influence results. If data were lacking to allow any definition, these cows were removed. The nonbreeding cows were considered to have no intention to treat from the date of designation. Therefore, these cows would contribute to HNBRED and HPREG up to the time that they were designated nonbreeding, at which point they would be censored. They would contribute to OPAL and would only contribute to PREG1 if mated. Studies ended at variable times; some allowed for insemination to 300 or more DIM, while others recorded insemination to d 150 before censoring outcomes (Table 2). Cows were right censored at 305 d or at the termination of the observation for pregnancy for the study. The OPAL and PREG1 were dichotomized outcomes, pregnant or not.
      Milk production and milk fat and protein percentage data were based on the approximately monthly milk recording information (DHIA) for each cow that was closest to first insemination for the herds, which was 69 DIM. A herd with daily milk obtained using in-parlor meters and milk component data from YieldSense+ (LIC Automation Services) had these estimates derived as described by
      • Golder H.
      • McGrath J.
      • Lean I.
      Supplementary material: Effect of 25-hydroxyvitamin D3 [25-(OH)D3] during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows. Figshare.
      . In 1 intensively fed herd in the United States, milk yield was recorded daily for individual cows using electronic milk meters (Perfection 3000, Boumatic), and then monthly means were calculated. Not all studies recorded milk or milk fat and protein information, and data from these cows were included in the univariable evaluations only.

      Sample Size Estimations

      The unit of interest in this study is the cow. For reproduction, we targeted a difference in HPREG of 10% with a power of 0.90, and α = 0.05, consequently requiring 1,899 cows per parity. The power calculations were made using Stata Version 16.1 (StataCorp).

      Statistical Analysis

      All statistical analysis was conducted using Stata Version 17 (StataCorp). Initial data evaluation included tabulation of data by categorical outcomes, visual appraisal of milk production, milk fat and protein percentage, and yield data for normality of distribution. Milk solids (yield of fat + protein kg/d) were calculated based on milk production and milk fat and protein percentage. The milk (kg/d), fat and protein percentages, and yields (kg/d) of milk solids, fat, and protein were centered by group to provide a model that accounted for the effects of these within each group. A quadratic term was generated for each of those centered variables.
      Survival analysis models were used to evaluate time to death or being culled, time to death only, or time to culling only with observations right censored at 305 DIM. The conditions for use of Cox's proportional hazards model for time to death or being culled were not met based on the Schoenfeld residual test, so Kaplan-Meier curves and the Weibull survival model were used. For death and sales, Cox models were tested against the Weibull survival model using Akaike's information criterion and Bayesian information criterion and found to have a lesser fit and very similar estimates of effects. Consequently, a Weibull model was used for these variables. The Weibull distribution is suitable for modeling data with monotone hazard rates that either increase or decrease exponentially with time and do not assume proportionality (Stata Version 17 StataCorp).
      Survival analysis models were also used to evaluate time not bred over the first 120 DIM (HNBRED) and HPREG. For HNBRED and HPREG, the conditions for use of Cox's proportional hazards model (stcox function) were assessed based on Schoenfeld residual tests and graphs of −log[log(Survival function)] against the natural log of time to event. The visual appraisal of Schoenfeld residuals indicated a good fit. The Cox model was tested against other survival models using Akaike's information criterion (AIC) and Bayesian information criterion (BIC) and found to have a better fit.
      The Stata 17 User's Guide (
      • StataCorp LLC
      Stata user's guide release 17.
      ) states that in a Cox model, the data are organized as i = 1,… n groups with j = 1,… ni observations in group i. For the jth observation in the ith group, the hazard is
      hij(t) = h0(ti exp(xijβ),


      where hij(t) is the hazard at time t for the jth observation in the ith group and h0(t) is the baseline hazard, αi is the group-level frailty or random effect, β is the coefficient for xij, the row vector of covariates for the time interval (t0ij; tij], where t0ij is the initial time of observation and tij is the observation at time t of the jth observation in the ith group in the dataset. The frailties are unobservable positive quantities and are assumed to have mean 1 and variance θ, to be estimated from the data. The frailties equate to random effects in other statistical models. After fitting of hazards models, predicted hazards of HPREG were estimated using the predict function.
      All models for survival and reproduction evaluated the fixed effects of parity (1 to ≥5), production system (pasture or intensively fed), and centered study year, and the random effect was cow nested within group. These effects were evaluated separately in all models. The interaction between parity and system was evaluated and retained in 2 models. For reproduction, the linear and quadratic effect of the centered production variables (milk yield or fat or protein percentage, fat yield, protein yield, or milk solids yield) were also evaluated and the random or frailty effect was group. Interactions of parity with production variables were tested and included when these improved the model fit as assessed by AIC and BIC. These effects were separately evaluated in models for each reproduction variable. A backward stepwise approach was used to remove nonsignificant terms (P > 0.05) from all survival models.
      The following mixed-effects, multilevel model using the melogit function provided an evaluation of OPAL or PREG1. The models included cow within group as a random effect.
      The Stata 17 User's Guide (
      • StataCorp LLC
      Stata user's guide release 17.
      ) states that the model represents a series of M independent clusters, and is conditional on a set of random effects uj,
      Pr(yij = 1|xij, uj) = H(xijβ + zijuj),


      where j = 1, … M clusters, with cluster j consisting of i = 1, … nj observations. The responses are the binary valued yij, with the convention of treating yij = 1 if depvarij (the dependent variable) ≠ 0, and treating yij = 0 otherwise. Pr is the probability for (yij = 1|xij, uj). The 1 × p row vector (xij) are the covariates for the fixed effects, analogous to the covariates found in a standard logistic regression model, with regression coefficients (fixed effects) β. The logistic cumulative distribution function is H, and zij are the covariates corresponding to the random effects. As the logistic cumulative distribution function, H maps the linear predictor to the probability of a success; that is, when yij = 1. After fitting of models, predicted OPAL or PREG1 were estimated using the predict function. The fixed effects tested were parity (1 to ≥5), production system (pasture or intensively fed), year, and the linear and quadratic effect of the centered production variables (milk yield or fat or protein percentage, fat yield, protein yield or milk solids yield), the interaction of each parity with the linear and the quadratic centered production variables, and the random effect was group.
      The effects of group within system were examined as random effects using the gllamm function to partition the variance components of the nested model (
      • Rabe-Hesketh S.
      • Skrondal A.
      Multilevel and Longitudinal Modeling Using Stata.
      ) for the base analyses without milk components, but the effect of system did not explain significant variation (<5%) in responses above that explained by group alone for OPAL; therefore, a 2-level model was used. The effects of milk at the group level were explored as a random slope, but did not influence OPAL and was not further explored. For milk fat or protein yield, random slope models did not converge. The mean effects of group level of milk or milk constituent yields were not significant as fixed effects when included with models containing system and were not included in final modeling. Change in parity over time was tested with the same model by examining the interaction of parity and year. The potential for interaction between parity and production variables was explored for OPAL and PREG1 by modeling within parity. Interactions of parity with production variables were tested and included when these improved the model fit as assessed by AIC and BIC. The number of studies and cows in intensive feeding systems increased over time. Consequently, effects of year of study on fertility outcomes were not included in final models.
      Collinearity among variables was subsequently explored for all models developed using the collin function that provides the variance inflation (VIF) index and condition number. The mean VIF for milk variables including the quadratic terms was 1.14 for milk, 1.07 for protein percentage, and 1.09 for fat percentage, indicating little effect of collinearity.

      RESULTS AND DISCUSSION

      This study provides insights into associations between parity and survival time, and between parity, milk production, and reproductive performance using rigorously collected data from 3 countries. All the studies used for evaluating reproductive performance had this as an important outcome and almost all pregnancy determinations made were the result of veterinary diagnosis, thereby reducing the risks of exclusion bias that exists in studies that are based on herd records without a reproductive focus. The use of purposively selected herds from 1992 to 2020 may introduce bias that could reduce the potential to apply these findings to all cows and herds. However, notwithstanding this limitation, these and companion studies (
      • Lean I.J.
      • Sheedy D.B.
      • LeBlanc S.J.
      • Duffield T.
      • Santos J.E.P.
      • Golder H.M.
      Holstein dairy cows lose body condition score and gain body weight with increasing parity in both pasture-based and total mixed ration herds.
      ,
      • Lean I.J.
      • LeBlanc S.J.
      • Sheedy D.B.
      • Duffield T.
      • Santos J.E.P.
      • Golder H.M.
      Associations of parity with health disorders and blood metabolite concentrations in Holstein cows in different production systems.
      ) represent the first studies since the 1980s with a detailed focus on parity. We provide insight into associations among parity, survival, deaths and sales, production, and pregnancy, indicating quadratic effects of milk production on reproductive performance and markedly reduced OPAL and HPREG in older cows. A possible limitation in our study is that genetic effects cannot be estimated. However, the risk of genetic confounding over time for the reproductive data is low. Generation intervals have been long and heritability of reproductive measures such as pregnancy to service are very low—that is, 0 to 0.05 (
      • Berry D.P.
      • Wall E.
      • Pryce J.E.
      Genetics and genomics of reproductive performance in dairy and beef cattle.
      )—and all effects are assessed within a temporal group of cows.

      Survival

      Of the 30,780 Holstein cows from 9 studies that evaluated survival time, 1,697 or 5.5% were in Australia, 2,123 (6.9%) were in Canada, and 26,960 or 87.6% were in the United States (Table 1). Of these, 1,094 cows were predominantly pasture fed. Cows in parity 1 were 40.4% (12,436), parity 2 were 27.3% (8,403), parity 3 were 16.0% (4,935), parity 4 were 8.6% (2,650), and parity ≥5 were 2,356 cows or 7.7% of the sample population. The median day to removal was 117, for death was 41, and sale was 135. The survival analysis data are dominated by studies and numbers of cows from the United States and could not be considered fully generalizable to the Canada and Australia populations, despite the lack of significant differences by production system or country.
      The hazard of removal for all reasons was associated with parity (P < 0.001) and the interaction of parity with system (P < 0.001) but not with system as a main effect (P = 0.957). Table 3 provides detail on the outcomes, and Figure 1 shows a Kaplan-Meier plot of the loss of cows that died and were sold by parity, with system as a covariable. There was a markedly increased hazard of removal for all reasons for older cows, with parity 3, 4, and ≥5 cows differing from parity 1, parity 4 and ≥5 differing from parity 2, and parity ≥5 differing from all other parities (Table 3; Figure 1). Parity 1 cows in intensively fed systems had a lesser hazard of removal than those in pasture systems (HR = 0.50; 95% CI 0.39 – 0.64), and parity 3 cows in intensively fed systems had a greater hazard of removal than those in pasture systems (HR = 1.55; 95% CI 1.00 – 2.40).
      Table 3Associations of parity, study year, and feeding system and their interactions with time to death or being sold, time to sale only, or time to death only; hazard ratio ± SE and P-values are reported from models including parity, study year, production system, and the interaction between parity and production system
      MeasurenParity
      The referent is parity 1.
      (P-value)
      Study yearProduction system
      TMR (intensively fed) versus pasture fed; the referent is pasture fed.
      Parity × system
      234≥5Main effect
      Died or sold30,7800.86 ± 0.17
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.71 ± 0.16
      Differs from parity 1 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.24
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.83 ± 0.28
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      <0.001
      Nonsignificant and excluded from the model.
      Nonsignificant and excluded from the model.
      <0.001
      (0.387)(0.459)(0.718)(<0.001)
      Sold30,7800.66 ± 0.16
      Differs from parity 5 at P < 0.05.
      0.58 ± 0.15
      Differs from parity 5 at P < 0.05.
      0.68 ± 0.21
      Differs from parity 5 at P < 0.05.
      1.54 ± 0.26
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      <0.001
      Nonsignificant and excluded from the model.
      Nonsignificant and excluded from the model.
      <0.001
      (0.082)(0.031)(0.217)(0.012)
      Died30,7801.57 ± 0.14
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.77 ± 0.25
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      3.95 ± 0.38
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      4.27 ± 0.43
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      <0.0010.001
      Nonsignificant and excluded from the model.
      Nonsignificant and excluded from the model.
      (<0.001)(<0.001)(<0.001)(<0.001)
      a Differs from parity 1 at P < 0.05.
      b Differs from parity 2 at P < 0.05.
      c Differs from parity 3 at P < 0.05.
      d Differs from parity 4 at P < 0.05.
      e Differs from parity 5 at P < 0.05.
      1 The referent is parity 1.
      2 TMR (intensively fed) versus pasture fed; the referent is pasture fed.
      3 Nonsignificant and excluded from the model.
      Figure thumbnail gr1
      Figure 1Kaplan-Meier survival curves for days to all removals (deaths and culls; x-axis) by parity for 30,780 cows in 9 prospective studies. The y-axis is the proportion of cows remaining. Based on the model that contains parity, study year, production system, and the interaction of parity with production system, cows in parity 1 differ from all other parities except parity 2. Cows in parity ≥5 differ from all other parities. Cows in parity 4 also differ from parity 2.
      Deaths were explained by a model that accounted for parity and year of study only (Table 3), with the hazard increasing markedly with increased parity and having a small increase over time. All parity comparisons were different except for parity 4 with ≥5.
      The pattern of sales was similar to responses for the combination of sales and deaths as the hazard of removal was associated with parity (P < 0.001) and the interaction of parity with system (P < 0.001) but not system as a main effect (P = 0.609). Parity ≥5 cows had greater hazard of sale than other parities, which did not differ among themselves (Table 3). However, the interaction of system with parity indicated that parity 1 cows in intensively fed systems had a reduced hazard of removal compared with those in pasture systems (HR = 0.46; 95% CI 0.33 – 0.63); no other within parity comparisons of production systems differed. Some of the sales will have been voluntary and this needs to be considered in the context of increased loss of older cows.
      The increased deaths with parity may reflect the increased risks for clinical hypocalcemia, mastitis, and lameness in cows of increased parity (
      • Lean I.J.
      • LeBlanc S.J.
      • Sheedy D.B.
      • Duffield T.
      • Santos J.E.P.
      • Golder H.M.
      Associations of parity with health disorders and blood metabolite concentrations in Holstein cows in different production systems.
      ). The marked differences in metabolic markers with increased parity identified in the companion paper may assist in identifying pathways that contribute to increased hazard of death.

      Reproduction

      A summary of the milk measures for each production system is provided in Table 4. Table 5 provides detail on the differences in OPAL, PREG1, and milk production variables for the different parities. Differences in production systems did not, generally, have significant associations with PREG1, but HPREG estimates indicated a reduced probability of pregnancy per day, and OPAL was markedly lower for intensively fed herds when milk constituents were considered (Tables 6 to 8). The variable system acted as a proxy for herd production level, as mean group milk was 31.4 ± 8.0 kg/d in pasture-fed groups and 42.2 ± 9.7 kg/d in intensively fed groups (Table 4). Consequently, we evaluated system rather than mean group milk production as
      • Rearte R.
      • LeBlanc S.J.
      • Corva S.G.
      • de la Sota R.L.
      • Lacau-Mengido I.M.
      • Giuliodori M.J.
      Effect of milk production on reproductive performance in dairy herds.
      did. Estimates of the effect of system also may have been influenced by changes in reproductive treatments including hormonal interventions and outcomes over time as the mean year of study was 2008 (standard deviation (SD) = 8.8 yr) for pasture-fed cows and 2013 (SD = 3.8 yr) for intensively fed cows. There was no significant difference in mean parity of 2.2 (SD 1.3) over time (P = 0.13).
      Table 4Summary of univariable reproduction and milk measures by production system either predominantly pasture fed or intensively fed (TMR) with cow numbers and mean ± SD or median estimates
      MeasurePasture fedTMR
      nMean ± SDMediannMean ± SDMedian
      Parity4,0832.90 ± 1.8728,6752.20 ± 1.37
      Reproduction measure
       Proportion bred4,1030.90 ± 0.3028,6750.88 ± 0.32
       Days to first breeding4,1037172
       Odds of becoming pregnant in a lactation (OPAL)4,1080.70 ± 0.4628,6750.72 ± 0.45
       Days to pregnancy4,1089097
       Proportion pregnancy to first breeding (PREG1)2,8620.48 ± 0.5020,5290.42 ± 0.49
      Milk measure
       Production (kg/d)3,18331.4 ± 7.9824,91742.2 ± 9.71
       Protein (%)2,5613.16 ± 0.326,2042.89 ± 0.26
       Fat (%)2,5603.78 ± 0.946,2013.84 ± 0.83
       Protein yield (kg/d)2,5611.05 ± 0.236,1911.28 ± 0.30
       Fat yield (kg/d)2,5601.24 ± 0.386,1911.71 ± 0.58
       TS (kg/d)2,5602.29 ± 0.536,1912.99 ± 0.83
      Table 5Univariable summary of reproductive and milk measures by parity with cow numbers and mean ± SD or median estimates
      VariableAll paritiesParity
      1234≥5
      nMean ± SDMediannMean ± SDMediannMean ± SDMediannMean ± SDMediannMean ± SDMediannMean ± SDMedian
      Reproduction measure
       Proportion bred32,7780.89 ± 0.3212,6110.90 ± 0.318,9380.92 ± 0.275,4540.88 ± 0.322,9360.83 ± 0.372,8130.80 ± 0.40
       Days to first breeding32,7787212,611728,938725,454722,936722,81371
       Odds of becoming pregnant in a lactation (OPAL)32,7830.71 ± 0.4512,6110.76 ± 0.438,9400.74 ± 0.445,4550.69 ± 0.462,9360.62 ± 0.482,8160.56 ± 0.50
       Days to pregnancy32,7839612,611948,940965,455992,9361002,81696
       Proportion pregnancy to first breeding (PREG1)23,3910.43 ± 0.499,5850.44 ± 0.506,6490.42 ± 0.493,7470.40 ± 0.491,8310.43 ± 0.501,5640.42 ± 0.49
      Milk yield (kg/d)
       Milk28,10041.0 ± 10.110,99635.3 ± 6.897,88644.1 ± 9.394,63045.8 ± 10.82,39946.2 ± 10.62,21543.0 ± 10.4
       Protein8,7521.20 ± 0.303,0851.03 ± 0.222,3761.28 ± 0.291,6851.36 ± 0.307771.34 ± 0.298251.24 ± 0.27
       Fat8,7511.58 ± 0.573,0851.31 ± 0.392,3761.71 ± 0.591,6841.80 ± 0.637771.75 ± 0.628251.56 ± 0.51
       TS8,7512.78 ± 0.823,0852.34 ± 0.562,3762.99 ± 0.821,6843.16 ± 0.877773.09 ± 0.868252.81 ± 0.70
      Table 6Associations of parity and milk production variables at approximately 70 DIM with hazard of pregnancy (HPREG) by d 305, controlling for parity, study year, and production system
      The HPREG was assessed at 3 study population levels, depending on data availability. Hazard ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The hazard ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      MeasurenParity
      The referent is parity 1.
      (P-value)
      Production system
      Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      Effect of milk production variable
      234≥5LinearQuadratic
      Whole population
       Univariate parity32,6960.90 ± 0.01
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.86 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.79 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.73 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (<0.001)(<0.001)(<0.001)(<0.001)
       Univariate parity censored32,6961.07 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.36 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.62 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.05 ± 0.07
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (0.015)(<0.001)(<0.001)(<0.001)
      Milk yield subpopulation
       Univariate parity28,0770.90 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.86 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.79 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.72 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (<0.001)(<0.001)(<0.001)(<0.001)
       Milk yield (kg/d)28,0770.92 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.90 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.83 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      0.75 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.67 ± 0.051.00 ± 0.001.00 ± 0.00
      (<0.001)(<0.001)(<0.001)(<0.001)(<0.001)(0.003)(<0.001)
      Component subpopulation
       Univariate parity8,7471.03 ± 0.03
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.90 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.69 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (0.375)(0.004)(0.048)(<0.001)
       Milk yield
      Linear interaction between milk yield (kg/d) and parity 1; hazard ratio = 0.96 ± 0.01; P < 0.001. Linear interaction between milk yield (kg/d) and parity 2; hazard ratio = 0.97 ± 0.00; P < 0.001. Linear interaction between milk yield (kg/d) and parity 3; hazard ratio = 0.98 ± 0.01; P = 0.003. Quadratic interaction between milk yield (kg/d) and parity 1; hazard ratio = 1.00 ± 0.00; P = 0.011. Quadratic interaction between milk yield (kg/d) and parity 2; hazard ratio = 1.00 ± 0.00; P = 0.002.
      (kg/d)
      8,7471.16 ± 0.05
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.09 ± 0.05
      Differs from parity 1 at P < 0.05.
      1.00 ± 0.06
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.77 ± 0.05
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.68 ± 0.111.01 ± 0.011.00 ± 0.01
      (0.001)(0.086)(0.993)(<0.001)(0.014)(0.080)(0.036)
       Milk fat
      Interaction between milk fat percentage and parity 1; hazard ratio = 1.12 ± 0.04; P = 0.001.
      (%)
      8,7571.03 ± 0.03
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.05
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.69 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.67 ± 0.100.90 ± 0.02
      Nonsignificant and excluded from the model.
      (0.390)(0.009)(0.061)(<0.001)(0.009)(<0.001)
       Fat yield (kg/d)8,7471.14 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.04 ± 0.04
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.05 ± 0.06
      Differs from parity 5 at P < 0.05.
      0.78 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.67 ± 0.110.78 ± 0.02
      Nonsignificant and excluded from the model.
      (<0.001)(0.369)(0.396)(<0.001)(0.011)(<0.001)
       Milk protein (%)8,7611.03 ± 0.03
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.00
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.05
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.69 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.67 ± 0.111.11 ± 0.060.75 ± 0.08
      (0.361)(0.009)(0.068)(<0.001)(0.010)(0.045)(0.005)
       Protein yield (kg/d)8,7481.12 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.02 ± 0.04
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.02 ± 0.06
      Differs from parity 5 at P < 0.05.
      0.77 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.67 ± 0.110.71 ± 0.04
      Nonsignificant and excluded from the model.
      (0.001)(0.673)(0.654)(<0.001)(0.011)(<0.001)
       TS (kg/d)8,7471.16 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.07 ± 0.04
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.07 ± 0.06
      Differs from parity 5 at P < 0.05.
      0.80 ± 0.05
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.67 ± 0.110.836 ± 0.02
      Nonsignificant and excluded from the model.
      (<0.001)(0.129)(0.191)(<0.001)(0.011)(<0.001)
      a Differs from parity 1 at P < 0.05.
      b Differs from parity 2 at P < 0.05.
      c Differs from parity 3 at P < 0.05.
      d Differs from parity 4 at P < 0.05.
      e Differs from parity 5 at P < 0.05.
      1 The HPREG was assessed at 3 study population levels, depending on data availability. Hazard ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The hazard ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      2 The referent is parity 1.
      3 Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      4 Linear interaction between milk yield (kg/d) and parity 1; hazard ratio = 0.96 ± 0.01; P < 0.001. Linear interaction between milk yield (kg/d) and parity 2; hazard ratio = 0.97 ± 0.00; P < 0.001. Linear interaction between milk yield (kg/d) and parity 3; hazard ratio = 0.98 ± 0.01; P = 0.003. Quadratic interaction between milk yield (kg/d) and parity 1; hazard ratio = 1.00 ± 0.00; P = 0.011. Quadratic interaction between milk yield (kg/d) and parity 2; hazard ratio = 1.00 ± 0.00; P = 0.002.
      5 Interaction between milk fat percentage and parity 1; hazard ratio = 1.12 ± 0.04; P = 0.001.
      6 Nonsignificant and excluded from the model.
      Table 8Associations of parity and milk production variables at approximately 70 DIM with odds of becoming pregnant in a lactation (OPAL), controlling for parity, study year, and production system
      The OPAL were assessed at 3 study population levels, depending on data availability. Odds ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The odds ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      MeasurenParity
      The referent is parity 1.
      (P-value)
      Production system
      Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      Effect of milk production variable
      234≥5LinearQuadratic
      Whole population
      Univariable parity33,4120.89 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.67 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.49 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.36 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (<0.001)(<0.001)(<0.001)(<0.001)
      Milk yield subpopulation
       Univariable parity28,0960.80 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.65 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.50 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.35 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (<0.001)(<0.001)(<0.001)(<0.001)
       Milk yield
      Linear interaction between milk yield (kg/d) and parity 1; odds ratio = 0.95 ± 0.01; P < 0.001. Quadratic interaction between milk yield (kg/d) and parity 1; odds ratio = 1.00 ± 0.00; P < 0.001. Linear interaction between milk yield (kg/d) and parity 2; odds ratio = 0.97 ± 0.00; P < 0.001.
      (kg/d)
      28,0960.68 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.48 ± 0.027
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.36 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.27 ± 0.02
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.56 ± 0.121.05 ± 0.001.0 ± 0.00
      (<0.001)(<0.001)(<0.001)(<0.001)(0.007)(<0.001)(<0.001)
      Milk yield subpopulation
       Univariable parity8,7471.01 ± 0.06
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.79 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.68 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.41 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (0.896)(0.001)(<0.001)(<0.001)
       Milk yield
      Linear interaction between milk yield (kg/d) and parity 1; odds ratio = 0.94 ± 0.01; P < 0.001. Linear interaction between milk yield (kg/d) and parity 2; odds ratio = 0.96 ± 0.01; P < 0.001. Linear interaction between milk yield (kg/d) and parity 3; odds ratio = 0.98 ± 0.01; P < 0.001. Quadratic interaction between milk yield (kg/d) and parity 1; odds ratio = 1.00 ± 0.00; P < 0.001. Quadratic interaction between milk yield (kg/d) and parity 2; odds ratio = 1.00 ± 0.00; P = 0.018.
      (kg/d)
      8,7470.95 ± 0.08
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.71 ± 0.07
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.61 ± 0.07
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.37 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.41 ± 0.151.03 ± 0.011.0 ± 0.01
      (0.581)(<0.001)(<0.001)(<0.001)(0.015)(<0.001)(<0.001)
       Milk fat (%)8,7571.03 ± 0.07
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.82 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.70 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.43 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.40 ± 0.150.89 ± 0.030.95 ± 0.02
      (0.684)(0.006)(<0.001)(<0.001)(0.011)(<0.001)(0.024)
       Fat yield
      Linear interaction between fat yield (kg/d) and parity 1; odds ratio = 1.04 ± 0.15; P = 0.803. Quadratic interaction between fat yield (kg/d) and parity 1; odds ratio = 0.36 ± 0.06 P < 0.001. Linear interaction between fat yield (kg/d) and parity 2; odds ratio = 0.75 ± 0.08; P = 0.006.
      (kg/d)
      8,7471.49 ± 0.37
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.73 ± 0.17
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.63 ± 0.15
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.39 ± 0.09
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.41 ± 0.150.86 ± 0.06
      Nonsignificant and excluded from the model.
      (0.113)(0.179)(0.056)(<0.001)(0.015)(0.048)
       Milk protein (%)8,7611.01 ± 0.07
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.80 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.68 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.41 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.40 ± 0.150.95 ± 0.090.66 ± 0.12
      (0.936)(0.002)(<0.001)(<0.001)(0.012)(0.608)(0.024)
       Protein yield
      Linear interaction between protein yield (kg/d) and parity 1; odds ratio = 0.02 ± 0.02; P < 0.001. Quadratic interaction between protein yield (kg/d) and parity 1; odds ratio = 0.20 ± 0.13; P = 0.012. Linear interaction between protein yield (kg/d) and parity 2; odds ratio = 0.01 ± 0.01; P < 0.001. Quadratic interaction between protein yield (kg/d) and parity 2; odds ratio = 3.53 ± 2.15; P = 0.039. Linear interaction between protein yield (kg/d) and parity 3; odds ratio = 0.02 ± 0.02; P < 0.001. Linear interaction between protein yield (kg/d) and parity 4; odds ratio = 0.02 ± 0.03; P = < 0.001. Linear interaction between protein yield (kg/d) and parity ≥5; odds ratio = 0.03 ± 0.04; P = 0.002.
      (kg/d)
      8,7481.90 ± 0.75
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.50 ± 0.21
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.33 ± 0.16
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.13 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Nonsignificant and excluded from the model.
      70.2 ± 73.90.36 ± 0.13
      (0.105)(0.094)(0.022)(<0.001)(<0.001)(0.003)
       TS
      Linear interaction between total milk solids (kg/d) and parity 1; odds ratio = 0.53 ± 0.08; P < 0.001. Quadratic interaction between total milk solids (kg/d) and parity 1; odds ratio = 0.62 ± 0.07; P < 0.001. Linear interaction between total milk solids (kg/d) and parity 2; odds ratio = 0.57 ± 0.07; P < 0.001. Quadratic interaction between total milk solids (kg/d) and parity 1; odds ratio = 1.28 ± 0.10; P = 0.002.
      (kg/d)
      8,7470.93 ± 0.08
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.74 ± 0.07
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.64 ± 0.07
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.39 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.41 ± 0.151.22 ± 0.100.84 ± 0.04
      (0.388)(0.001)(<0.001)(<0.001)(0.015)(0.019)(0.001)
      a Differs from parity 1 at P < 0.05.
      b Differs from parity 2 at P < 0.05.
      c Differs from parity 3 at P < 0.05.
      d Differs from parity 4 at P < 0.05.
      e Differs from parity 5 at P < 0.05.
      1 The OPAL were assessed at 3 study population levels, depending on data availability. Odds ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The odds ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      2 The referent is parity 1.
      3 Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      4 Linear interaction between milk yield (kg/d) and parity 1; odds ratio = 0.95 ± 0.01; P < 0.001. Quadratic interaction between milk yield (kg/d) and parity 1; odds ratio = 1.00 ± 0.00; P < 0.001. Linear interaction between milk yield (kg/d) and parity 2; odds ratio = 0.97 ± 0.00; P < 0.001.
      5 Linear interaction between milk yield (kg/d) and parity 1; odds ratio = 0.94 ± 0.01; P < 0.001. Linear interaction between milk yield (kg/d) and parity 2; odds ratio = 0.96 ± 0.01; P < 0.001. Linear interaction between milk yield (kg/d) and parity 3; odds ratio = 0.98 ± 0.01; P < 0.001. Quadratic interaction between milk yield (kg/d) and parity 1; odds ratio = 1.00 ± 0.00; P < 0.001. Quadratic interaction between milk yield (kg/d) and parity 2; odds ratio = 1.00 ± 0.00; P = 0.018.
      6 Linear interaction between fat yield (kg/d) and parity 1; odds ratio = 1.04 ± 0.15; P = 0.803. Quadratic interaction between fat yield (kg/d) and parity 1; odds ratio = 0.36 ± 0.06 P < 0.001. Linear interaction between fat yield (kg/d) and parity 2; odds ratio = 0.75 ± 0.08; P = 0.006.
      7 Nonsignificant and excluded from the model.
      8 Linear interaction between protein yield (kg/d) and parity 1; odds ratio = 0.02 ± 0.02; P < 0.001. Quadratic interaction between protein yield (kg/d) and parity 1; odds ratio = 0.20 ± 0.13; P = 0.012. Linear interaction between protein yield (kg/d) and parity 2; odds ratio = 0.01 ± 0.01; P < 0.001. Quadratic interaction between protein yield (kg/d) and parity 2; odds ratio = 3.53 ± 2.15; P = 0.039. Linear interaction between protein yield (kg/d) and parity 3; odds ratio = 0.02 ± 0.02; P < 0.001. Linear interaction between protein yield (kg/d) and parity 4; odds ratio = 0.02 ± 0.03; P = < 0.001. Linear interaction between protein yield (kg/d) and parity ≥5; odds ratio = 0.03 ± 0.04; P = 0.002.
      9 Linear interaction between total milk solids (kg/d) and parity 1; odds ratio = 0.53 ± 0.08; P < 0.001. Quadratic interaction between total milk solids (kg/d) and parity 1; odds ratio = 0.62 ± 0.07; P < 0.001. Linear interaction between total milk solids (kg/d) and parity 2; odds ratio = 0.57 ± 0.07; P < 0.001. Quadratic interaction between total milk solids (kg/d) and parity 1; odds ratio = 1.28 ± 0.10; P = 0.002.
      There were 32,783 cows in the 13 studies conducted in Australia (14.6% of cows), Canada (2.4% of cows), and the United States (83.0% of cows). After evaluating the data, we found only 8.6% were parity ≥5, and these were combined into a single group. A plurality of cows were in parity 1 with 12,611 cows or 38.5% of the sample, 27.3% were in parity 2, 16.7% were in parity 3, and 9.0% were in parity 4, indicating considerable loss of cows with increased parity. Consequently, 30.1% of cows entering parity 1 would not enter parity 2, 39.0% of parity 2 were lost before parity 3, and 46.2% of cows entering parity 3 would not enter parity 4. There were 4,108 pasture-based cows and 28,675 intensively fed cows. The sample size requirements for the study were approached or exceeded as the smallest subgroup was parity 4 in the milk fat yield model (n = 777 cows).
      Controlling for the effects of parity reduces the potentially confounding effects of inherently greater milk production, with increasing parity influencing associations among parity, production level and composition, and reproductive outcomes. The use of milk, milk fat and protein percentages and yields, and milk solids production estimates close to first insemination (72.8 DIM on average; median 72 d) minimizes the potential for the effects of pregnancy to influence lactational persistency of production and for selective culling and breeding decisions to bias relationships (
      • Lean I.J.
      • Galland J.
      • Scott J.
      Relationships between fertility, peak milk yields and lactational persistency in dairy cows.
      ;
      • Santos J.E.P.
      • Rutigliano H.M.
      • Sá Filho M.F.
      Risk factors for resumption of postpartum estrous cycles and embryonic survival in lactating dairy cows.
      ;
      • LeBlanc S.
      Assessing the association of the level of milk production with reproductive performance in dairy cattle.
      ;
      • Rearte R.
      • LeBlanc S.J.
      • Corva S.G.
      • de la Sota R.L.
      • Lacau-Mengido I.M.
      • Giuliodori M.J.
      Effect of milk production on reproductive performance in dairy herds.
      ). However, for low-producing cows, there is still the potential for removal bias as producers may choose to remove these.

      Time to Censoring for Breeding

      The analyses presented in Table 6, Table 9 provide evaluations of not being inseminated up to d 120 (HNBRED) and HPREG, respectively, as estimated by the hazard function over time. The univariable Kaplan-Meier survival graphs for HNBRED and HPREG are presented in Figure 2, Figure 3, respectively. The median times to first insemination were 72 d for parities 1 to 4 and 71 d for parity 5. Seventy-five percent of cows in all parities were inseminated by 80 to 82 DIM.
      Table 9Associations of parity and milk production variables at approximately 70 DIM with the hazard of not being bred (HNBRED) by 120 DIM, controlling for parity, study year, and production system
      The HNBRED was assessed at 3 study population levels, depending on data availability. Hazard ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The hazard ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      MeasurenParity
      The referent is parity 1.
      (P-value)
      Production system
      Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      Effect of milk production variable
      234≥5LinearQuadratic
      Whole population
       Univariate parity32,6850.82 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.24 ± 0.06
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.84 ± 0.10
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.45 ± 0.13
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (<0.001)(<0.001)(<0.001)(<0.001)
      Milk yield subpopulation
       Univariate parity28,0710.80 ± 0.0
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.29 ± 0.13
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.39 ± 0.25
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      3.53 ± 0.36
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      (0.028)(0.013)(<0.001)(<0.001)
       Milk yield
      Linear interaction between milk yield (kg/d) and parity 1; hazard ratio = 0.97 ± 0.01; P = 0.003.
      (kg/d)
      28,0711.95 ± 0.27
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      3.42 ± 0.50
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      6.64 ± 0.98
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      8.82 ± 1.26
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      3.09 ± 1.050.95 ± 0.001.0 ± 0.00
      (<0.001)(<0.001)(<0.001)(<0.001)(0.001)(<0.001)(<0.001)
      Component subpopulation
       Univariate parity8,7470.72 ± 0.10
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.97 ± 0.14
      Differs from parity 5 at P < 0.05.
      1.35 ± 0.23
      Differs from parity 1 at P < 0.05.
      1.52 ± 0.28
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      (0.016)(0.840)(0.079)(0.020)
       Milk yield (kg/d)8,7471.33 ± 0.20
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.09 ± 0.34
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      3.03 ± 0.58
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      3.00 ± 0.58
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      6.71 ± 3.250.95 ± 0.001.0 ± 0.00
      (0.057)(<0.001)(<0.001)(<0.001)(<0.001)(<0.001)(<0.001)
       Milk fat (%)8,7570.69 ± 0.10
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.91 ± 0.13
      Differs from parity 5 at P < 0.05.
      1.28 ± 0.22
      Differs from parity 2 at P < 0.05.
      1.44 ± 0.26
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      6.71 ± 3.141.28 ± 0.08
      Nonsignificant and excluded from the model.
      (0.007)(0.503)(0.156)(0.044)(<0.001)(<0.001)
       Fat yield
      Linear interaction between milk fat yield (kg/d) and parity 1; hazard ratio = 1.13 ± 0.29; P = 0.628. Quadratic interaction between fat yield (kg/d) and parity 1; hazard ratio = 3.71 ± 0.67; P < 0.001.
      (kg/d)
      8,7471.40 ± 0.56
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.98 ± 0.81
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.77 ± 1.16
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      3.01 ± 1.24
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      6.91 ± 3.321.13 ± 0.10
      Nonsignificant and excluded from the model.
      (0.398)(0.093)(0.014)(0.008)(<0.001)(0.023)
       Milk protein (%)8,7610.71 ± 0.40
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.95 ± 0.14
      Differs from parity 5 at P < 0.05.
      1.33 ± 0.23
      Differs from parity 2 at P < 0.05.
      1.51 ± 0.27
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      6.91 ± 3.320.90 ± 0.171.96 ± 0.61
      (0.011)(0.704)(0.099)(0.022)(<0.001)(0.596)(0.032)
       Protein yield (kg/d)8,7481.27 ± 0.19
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      2.01 ± 0.33
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      2.87 ± 0.54
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      2.81 ± 0.54
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      6.94 ± 3.330.17 ± 0.034.81 ± 1.30
      (0.104)(<0.001)(<0.001)(<0.001)(<0.001)(<0.001)(<0.001)
       TS (kg/d)8,7471.14 ± 0.17
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      1.79 ± 0.29
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      2.50 ± 0.47
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      2.58 ± 0.50
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      6.81 ± 3.260.54 ± 0.041.27 ± 0.04
      (0.391)(<0.001)(<0.001)(<0.001)(<0.001)(<0.001)(<0.001)
      a Differs from parity 1 at P < 0.05.
      b Differs from parity 2 at P < 0.05.
      c Differs from parity 3 at P < 0.05.
      d Differs from parity 4 at P < 0.05.
      e Differs from parity 5 at P < 0.05.
      1 The HNBRED was assessed at 3 study population levels, depending on data availability. Hazard ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The hazard ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      2 The referent is parity 1.
      3 Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      4 Linear interaction between milk yield (kg/d) and parity 1; hazard ratio = 0.97 ± 0.01; P = 0.003.
      5 Nonsignificant and excluded from the model.
      6 Linear interaction between milk fat yield (kg/d) and parity 1; hazard ratio = 1.13 ± 0.29; P = 0.628. Quadratic interaction between fat yield (kg/d) and parity 1; hazard ratio = 3.71 ± 0.67; P < 0.001.
      Figure thumbnail gr2
      Figure 2Kaplan-Meier survival curves for the days to removal before first insemination [i.e., the hazard of not being bred (HNBRED) or censoring for the hazard of being inseminated at least once] by parity. The x-axis is DIM, and the y-axis is the proportion of cows not bred. The data are from 32,356 cows. Based on the model that contains parity, study year, and production system, cows in parity ≥5 had 2.45 times greater HNBRED than cows in parity 1. Cows in parity 2 had reduced hazard of HNBRED, and parities 3 and 4 had greater hazard of not being inseminated than parity 1 and lesser hazard than parity ≥5.
      Figure thumbnail gr3
      Figure 3Kaplan-Meier survival curves for days to pregnancy by parity. The data are from 32,969 cows in 13 studies. The x-axis is DIM, and the y-axis is the proportion of cows pregnant. Based on the model that contains parity, study year, and production system, all parity groups differed in the daily probability of pregnancy.
      Eighty-eight percent of cows had a first breeding. Of the 12% of cows that were never inseminated, 75% were removed by 57 DIM, suggesting that many of these died or were culled for disease or low production. Presumably, high-producing cows that were selected not to be rebred would be kept in the herd producing at least until after peak production. Therefore, cows removed in early lactation probably include those with a primary health disorder resulting in low production, a health risk despite acceptable production (e.g., chronically high SCC), or low production without a diagnosed health disorder The pattern of censoring in Figure 2 and Table 9 show that parity ≥5 cows had 2.5 times greater hazard of not being inseminated (removed before first insemination) than parity 1 cows. In the univariable estimates of parity associations with HNBRED in the subsample of cows with milk fat and protein data, cows with parity ≥5 had 1.5 times greater HNBRED than parity 1 and greater HNBRED than parity 2 or 3. Parity 2 cows had a lesser HNBRED than parity 1, but other comparisons among parity did not differ aside from parities 2 and 4.
      When the quadratic effect of milk yield and the effect of system were included in the model (n = 28,071), all parities had greater hazards of not being inseminated than parity 1, with the hazard being >5 times greater for parities 4 and ≥5. Figure 4 shows that only cows with milk yield approximately 15 kg/d less than group mean were less likely to be inseminated at all. The hazard of not being inseminated was 3 times greater in intensively managed cows than in the pasture-based system (Table 9). For the subsample with data on milk fat and protein percentages, a 1% point increase in milk fat percentage above group mean was associated with a 28% increase in the hazard of not being inseminated (Table 9). This finding may indicate that cows mobilizing more lipids that transferred to milk may be at more risk of not being inseminated, a finding consistent with associations between elevated blood free fatty acids and reduced HPREG (
      • Westwood C.T.
      • Lean I.J.
      • Garvin J.K.
      Factors influencing fertility of Holstein dairy cows: A multivariate description.
      ). Pairwise differences among parities are detailed in Table 9. As seen in the models without production data, production system was associated with almost 7 times greater hazard of HNBRED for intensively fed than pasture-fed cows (Table 9) in all milk constituent models. However, there was an interaction between parity 1 cows and fat yield such that parity 1 cows had more risk of being not inseminated with a low or higher fat yield compared with the group mean. This finding may reflect a more aggressive removal policy for poor producing parity 1 cows and possibly increased risk of anestrus with higher production of FCM (
      • Westwood C.T.
      • Lean I.J.
      • Garvin J.K.
      Factors influencing fertility of Holstein dairy cows: A multivariate description.
      ).
      Figure thumbnail gr4
      Figure 4The association of milk yield (kg/d) centered on group within herd at approximately 70 DIM (x-axis) with the relative predicted hazard of not being inseminated by 120 DIM on the y-axis. The data are from 28,071 cows in 7 studies. Based on the model that contains parity, study year, and production system, all parities had increased hazard of not being bred (HNBRED) to parity 1.
      There was a quadratic effect of milk protein percentage on HNBRED. Controlling for milk protein percentage and system, parity 1 cows were more likely to be inseminated than parity 4 or ≥5, less likely than parity 2, and similar to parity 3 (Table 9). Parity 4 cows had increased HNBRED to parity 2 and did not differ from other parity cows. Similarly, also accounting for system, there was a quadratic association of protein yield with HNBRED (HR linear = 0.17; P < 0.001; quadratic = 4.81; P < 0.001). Parity 1 and 2 cows did not differ in HNBRED but were more likely to be inseminated than parity 3, 4, or ≥5 (Table 9). Similar associations were observed considering milk solids yield (Table 9).
      The metabolic demands of lactation, as indicated by increased milk volume and solids yield were associated with small increases in the hazard of not being inseminated. However, lower-producing cows within a herd, cows of parity ≥5, and cows in the intensively fed herds were less likely to be bred (Table 9, Figure 4).

      Time to Pregnancy

      Univariable estimates indicated that increased parity was associated with reduced HPREG (Table 6). Associations of increased parity with reduced HPREG remained in models accounting for system and milk yield (Table 6, Figure 5).
      • LeBlanc S.
      Assessing the association of the level of milk production with reproductive performance in dairy cattle.
      reported that cows in the top quartile of milk yield at first test day had a reduced time to first service and to pregnancy compared with cows in the quartile with least milk production.
      • Rearte R.
      • LeBlanc S.J.
      • Corva S.G.
      • de la Sota R.L.
      • Lacau-Mengido I.M.
      • Giuliodori M.J.
      Effect of milk production on reproductive performance in dairy herds.
      found that cows in herds in the top quartile of milk yield that produced more than 1 SD (8 kg/d) above the herd mean daily milk at 63 DIM had 1.3% lower probability of pregnancy per day (HR = 0.99), whereas they had 14.8% greater HPREG (HR = 1.15) in herds in the lowest quartile of milk yield. In the present study, the interquartile range was from −6.0 kg/d below to 5.7 kg/d above the group mean. Figure 5 shows that HPREG was reduced both above and below the interquartile range except for parity 2 cows.
      Figure thumbnail gr5
      Figure 5The association of milk yield (kg/d) centered on group within herd at approximately 70 DIM on the x-axis with the relative predicted hazard of being pregnant on the y-axis. The data are from 20,071 cows. Based on the model that contains parity, study year, and production system, all parity groups differed in probability or hazard of pregnancy (HPREG). Cows in parity 1 had the greatest daily probability or HPREG, and parity ≥5 had the least probability of pregnancy.
      Other studies found differing and conditional associations of milk production with reproduction.
      • LeBlanc S.
      Assessing the association of the level of milk production with reproductive performance in dairy cattle.
      similarly found lower reproductive performance in cows with lower milk production in early lactation but did not detect a quadratic relationship of milk yield with time to pregnancy.
      • Rearte R.
      • LeBlanc S.J.
      • Corva S.G.
      • de la Sota R.L.
      • Lacau-Mengido I.M.
      • Giuliodori M.J.
      Effect of milk production on reproductive performance in dairy herds.
      also found that higher-yielding cows had shorter time to pregnancy but only in lower-yielding herds. In higher-yielding herds, higher-yielding cows had longer times to pregnancy, although the magnitude of the effect was small. In a large study in the United Kingdom (
      • Madouasse A.
      • Huxley J.
      • Browne W.
      • Bradley A.
      • Dryden I.
      • Green M.
      Use of individual cow milk recording data at the start of lactation to predict the calving to conception interval.
      ), the strongest predictors of pregnancy by 145 DIM were milk yield and protein percentage at second herd test 2 and lactose percentage at first herd test after calving, respectively.
      • Bello N.M.
      • Steibel J.
      • Erskine R.
      • Tempelman R.
      Cows and herds constitute distinct hierarchical levels of heterogeneity in the variability of and association between milk yield and pregnancy outcome in dairy cows.
      found no association of cow-level milk yield but an antagonistic association of herd-level milk yield with pregnancy at first AI, with heterogeneity among herds as functions of the relative level of production and the baseline probability of pregnancy at first AI. These studies and the present results underline the desirability of assessing the nature of the cow-level association of milk production with reproduction (i.e., linear vs. quadratic relationships), their herd-level associations, and conditional relationships between the 2 levels of evaluation. Specifically, greater production at the herd level need not be associated with reduced reproductive performance, but within a herd or group, there is evidence that milk and milk solids production greater or lower than the mean are associated with reduced performance.
      When milk constituent data were evaluated in the models, increased parity remained associated with reduced HPREG (Table 6). In these models, cows in intensively fed groups consistently had a lower HPREG (HR ∼0.67) than pasture-fed cows. This association of production system with HPREG represents not just the effects of feeding and management system, but also greater mean milk yield and may be confounded by effects of changes in HPREG over the time during which the studies were conducted. Increased fat percentage (HR = 0.90) was associated with a lower HPREG. Milk fat yield also was associated (HR = 0.78) with HPREG; each additional kg of milk fat was associated with a 22.6% reduction in HPREG. In the model with fat yield, the associations between parity and HPREG changed slightly as detailed in Table 6.
      Consistently among parities, milk protein percentage was associated quadratically with HPREG (HR linear = 1.11 and HR quadratic = 0.75; Figure 6). The predicted HPREG increased up to approximately 0.5 percentage units above the mean protein percentage, after which HPREG decreased. Conversely, milk protein yield had a linear negative association with HPREG (HR = 0.71; Table 6). Accounting for this, parity ≥5 cows had lesser HPREG than all other parity groups. The total milk solids model found a linear decrease in HPREG associated with increased milk solids yield (HR = 0.84); however, accounting for the influence of milk solids yield altered parity effect estimates. Parity 2 cows had greater HPREG than parity 1, 3, or ≥5, but did not differ to parity 4. Parity ≥5 had lesser HPREG than all other groups.
      Figure thumbnail gr6
      Figure 6The association of milk protein percentage centered on group within herd at approximately 70 DIM with the relative predicted hazard of pregnancy (HPREG). The x-axis is milk protein percentage centered by group, and the y-axis is the predicted hazard ratio of being pregnant. The data are from 8,761 cows in 7 studies. Based on the model that contains parity, study year, and production system, cows in parity ≥5 had lesser daily probability of pregnancy than other parities. Parity 1 did not differ from parity 2; however, cows in parity 3 and 4 had a lesser HPREG than cows in parity 1 and 2, which did not differ from each other.
      In summary, accounting for measures of milk production, associations of parity with HPREG were reasonably consistent (Table 6), with parity ≥5 cows consistently having the lowest HPREG, parity 1 and 2 being similar, and parity 3 and 4 showing some reduction in HPREG. The reduction in differences among parities 1 to 4 when yields of fat and protein were included suggest that output of fat and protein explained much of the difference in HPREG, whereas accounting for milk outputs did not obviate the reduced HPREG for parity ≥5 cows. The larger database that evaluated milk yield showed a linear decline in HPREG with age.

      Pregnancy at First Breeding

      The differences in full lactation HPREG suggested a need to evaluate PREG1; the time to first breeding (HBRED1) was similar among parities, yet there were differences in time to pregnancy, particularly for parity ≥5. The random effect of group explained only 3.7% of the variance in PREG1. In the univariable analysis, parity 1 cows had greater PREG1 than most other parity groups (Table 7). For the model including milk yield, greater milk yield was associated with lesser odds of PREG1 [odds ratio (OR) linear = 1.00; Figure 7]; few differences among parities were evident. Parity 2 cows had an interaction with milk production indicating a greater effect of milk production on odds of PREG1 (Figure 7).
      • Pryce J.E.
      • Nielsen B.L.
      • Veerkamp R.F.
      • Simm G.
      Genotype and feeding system effects and interactions for health and fertility traits in dairy cattle.
      found first service conception was greater for parity 2 to 4 cows than those in first or parity greater than 4.
      • Lean I.J.
      • Galland J.
      • Scott J.
      Relationships between fertility, peak milk yields and lactational persistency in dairy cows.
      found that cows with higher peak milk yield had a lower probability of pregnancy from the first 2 matings and that cows in parity 1 and 2 had a greater probability of pregnancy than older cows. Similar to
      • Lean I.J.
      • Galland J.
      • Scott J.
      Relationships between fertility, peak milk yields and lactational persistency in dairy cows.
      ,
      • Rearte R.
      • LeBlanc S.J.
      • Corva S.G.
      • de la Sota R.L.
      • Lacau-Mengido I.M.
      • Giuliodori M.J.
      Effect of milk production on reproductive performance in dairy herds.
      assessed the odds of pregnancy within 100 DIM and found that parity 2 cows did not differ from parity 1 and that parity 3 and ≥4 had reduced odds.
      Table 7Associations of parity and milk production variables at approximately 70 DIM with odds of pregnancy to first breeding (PREG1), controlling for parity, study year, and production system
      The PREG1 was assessed at 3 study population levels depending on data availability. Odds ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The odds ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      MeasurenParity
      The referent is parity 1.
      (P-value)
      Production system
      Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      Effect of milk production variable
      234≥5LinearQuadratic
      Whole population23,637
       Univariable parity23,6370.87 ± 0.03
      Differs from parity 1 at P < 0.05.
      0.81 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.92 ± 0.05
      Differs from parity 3 at P < 0.05.
      0.86 ± 0.05
      Differs from parity 1 at P < 0.05.
      (<0.001)(<0.001)(0.115)(0.009)
      Milk yield subpopulation
       Univariable parity21,9640.87 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 3 at P < 0.05.
      0.79 ± 0.03
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      0.94 ± 0.05
      Differs from parity 3 at P < 0.05.
      0.86 ± 0.05
      Differs from parity 1 at P < 0.05.
      (<0.001)(<0.001)(0.227)(0.014)
       Milk yield
      Linear interaction between milk yield (kg/d) and parity 2; odds ratio = 0.98 ± 0.00; P < 0.001.
      (kg/d)
      21,9640.96 ± 0.04
      Differs from parity 3 at P < 0.05.
      0.84 ± 0.04
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      1.00 ± 0.060.92 ± 0.06
      Nonsignificant and excluded from the model.
      1.00 ± 0.00
      Nonsignificant and excluded from the model.
      (0.260)(0.001)(0.988)(0.195)(0.048)
      Component subpopulation
       Univariable parity6,2771.11 ± 0.071.03 ± 0.081.19 ± 0.121.04 ± 0.11
      (0.098)(0.687)(0.085)(0.745)
       Milk yield (kg/d)6,2771.27 ± 0.09
      Differs from parity 1 at P < 0.05.
      1.17 ± 0.111.37 ± 0.16
      Differs from parity 1 at P < 0.05.
      1.18 ± 0.14
      Nonsignificant and excluded from the model.
      0.99 ± 0.000.97 ± 0.01
      (0.001)(0.083)(0.007)(0.159)(0.027)(<0.001)
       Milk fat (%)6,2841.12 ± 0.071.05 ± 0.081.20 ± 0.121.05 ± 0.11
      Nonsignificant and excluded from the model.
      0.93 ± 0.03
      Nonsignificant and excluded from the model.
      (0.073)(0.539)(0.064)(0.644)(0.022)
       Fat yield
      Linear interaction between milk fat yield (kg/d) and parity 1; odds ratio = 0.60 ± 0.07; P < 0.001. Linear interaction between milk fat yield (kg/d) and parity 2; odds ratio = 0.63 ± 0.06; P < 0.001.
      (kg/d)
      6,2771.21 ± 0.26
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.53 ± 0.09
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      0.62 ± 0.12
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      0.56 ± 0.10
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Nonsignificant and excluded from the model.
      NS
      Nonsignificant and excluded from the model.
      (0.380)(<0.001)(0.010)(0.002)
       Milk protein (%)6,2841.12 ± 0.071.04 ± 0.081.20 ± 0.121.05 ± 0.11
      Nonsignificant and excluded from the model.
      1.16 ± 0.120.56 ± 0.11
      (0.085)(0.594)(0.063)(0.673)(0.168)(0.004)
       Protein yield
      Linear interaction between protein yield (kg/d) and parity 1; odds ratio = 0.38 ± 0.11; P = 0.001. Quadratic interaction between protein yield (kg/d) and parity 1; odds ratio = 0.154 ± 0.110; P = 0.009. Linear interaction between protein yield (kg/d) and parity 2; odds ratio = 0.47 ± 0.11; P = 0.001.
      (kg/d)
      6,2770.99 ± 0.37
      Differs from parity 3 at P < 0.05.
      Differs from parity 4 at P < 0.05.
      Differs from parity 5 at P < 0.05.
      0.37 ± 0.13
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      0.44 ± 0.16
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      0.39 ± 0.14
      Differs from parity 1 at P < 0.05.
      Differs from parity 2 at P < 0.05.
      Nonsignificant and excluded from the model.
      0.76 ± 0.13
      Nonsignificant and excluded from the model.
      (0.973)(0.006)(0.021)(0.008)(0.101)
       TS
      Linear interaction between TS yield (kg/d) and parity 2; odds ratio = 0.72 ± 0.07; P < 0.001.
      (kg/d)
      6,2771.23 ± 0.09
      Differs from parity 1 at P < 0.05.
      1.16 ± 0.101.35 ± 0.15
      Differs from parity 1 at P < 0.05.
      1.17 ± 0.14
      Nonsignificant and excluded from the model.
      0.88 ± 0.05
      Nonsignificant and excluded from the model.
      (0.003)(0.089)(0.008)(0.185)(0.021)
      a Differs from parity 1 at P < 0.05.
      b Differs from parity 2 at P < 0.05.
      c Differs from parity 3 at P < 0.05.
      d Differs from parity 4 at P < 0.05.
      e Differs from parity 5 at P < 0.05.
      1 The PREG1 was assessed at 3 study population levels depending on data availability. Odds ratio ± SE and P-values are reported from models including only parity, or parity, study year, production system, and a milk production variable. Study year was not significant and was excluded from models. Interaction terms are listed in the footnotes below. The odds ratio is estimated for a 1-kg or 1-percentage point difference for each respective measure.
      2 The referent is parity 1.
      3 Production systems were pasture fed or TMR (intensively fed); the referent is pasture fed.
      4 Linear interaction between milk yield (kg/d) and parity 2; odds ratio = 0.98 ± 0.00; P < 0.001.
      5 Nonsignificant and excluded from the model.
      6 Linear interaction between milk fat yield (kg/d) and parity 1; odds ratio = 0.60 ± 0.07; P < 0.001. Linear interaction between milk fat yield (kg/d) and parity 2; odds ratio = 0.63 ± 0.06; P < 0.001.
      7 Linear interaction between protein yield (kg/d) and parity 1; odds ratio = 0.38 ± 0.11; P = 0.001. Quadratic interaction between protein yield (kg/d) and parity 1; odds ratio = 0.154 ± 0.110; P = 0.009. Linear interaction between protein yield (kg/d) and parity 2; odds ratio = 0.47 ± 0.11; P = 0.001.
      8 Linear interaction between TS yield (kg/d) and parity 2; odds ratio = 0.72 ± 0.07; P < 0.001.
      Figure thumbnail gr7
      Figure 7The association of milk yield (kg/d) centered on group at approximately 70 DIM with the predicted coefficient for odds of pregnancy at first insemination. The data are from 28,071 cows. The x-axis is milk yield centered by group, and the y-axis is the predicted coefficient of odds of pregnancy at first insemination. Based on the model that contains parity and production system, and the interaction of parity 2 with production system, odds of pregnancy at first insemination were greater in cows in parity 1 than parity 3 but did not differ from other parities. Cows in parity 4 had greater odds of pregnancy at first insemination than cows in parity 2 and 3.
      Production system did not influence the odds of PREG1 for any of the subpopulation models, which differed from HPREG and HNBRED, with which system had a marked association. The inclusion of time in outcomes HPREG and HNBRED is an important difference to PREG1. However, milk component production was associated with pregnancy at first breeding. For the subsample with data on solids percentages and yields, when the decrease in PREG1 associated with increased fat percentage was accounted for, there were no differences in pregnancy at first breeding among parities (Table 7). Fat yields had an interaction for parity 1 and parity 2 cows that influenced differences among parities with parity 1 (OR = 0.60) and 2 cows (OR = 0.63), which declined in odds of PREG1 with increased fat yield and resulted in differences in odds of PREG1 for parity 1 and 2 to the parity 3 to 5 cows. Protein percentage had a quadratic association with PREG1 and again no differences were found among parities (Table 7). In contrast to protein percentage effects, increased protein yield was associated with a reduction in the odds of PREG1 and, similarly to fat yield interactions, with protein yield for parity 1 (OR linear = 0.38 and quadratic = 0.14) and parity 2 (OR = 0.47). Accounting for that effect, parity 1 and 2 had lower odds of PREG1 (Table 7).
      • Rodney R.M.
      • Celi P.
      • Scott W.
      • Breinhild K.
      • Santos J.
      • Lean I.
      Effects of nutrition on the fertility of lactating dairy cattle.
      also found that increased milk protein yield reduced the probability of pregnancy, with a 1 SD (0.2 kg) increase in milk protein resulting in an 11% reduction in the odds of pregnancy, a result consistent with the findings in this study. The milk solids model indicated that increased solids yield was associated with lower odds of PREG1 (OR = 0.88). There was an interaction with milk solids yield for parity 2 cows, resulting in a more marked decline in PREG1 associated with increased milk solids production than for other parities. Parity differences were limited to parity 2 and 4 cows with greater odds of PREG1 for cows of parity >1 (Table 7).
      • Hudson C.D.
      • Green M.J.
      Associations between routinely collected dairy herd improvement data and insemination outcome in UK dairy herds.
      found conditional associations of milk yield and components with the pregnancy risk per AI (P/AI). Generally, yield at the first 2 tests and protein percentage at the first test were positively associated with P/AI while fat percentage at the first test and protein percentage at the second test were negatively associated. However, combined, these effects explained < 8% of the variance in P/AI; most of the variance was at the herd level with unmeasured effects. Therefore, milk yield or components are associated with measures of reproductive performance but may account for small fractions of the variance in these multifactorial outcomes. Large-scale studies with comprehensive cow- and herd-level data will be needed to accurately partition the effects of age, nutrition, level of milk yield, and health conditions.

      Odds of Pregnancy in a Lactation

      The univariable association between parity and pregnancy in a lactation was examined in 33,412 cows, which was all cows that calved. As parity increased, OPAL was reduced (Table 8). The random effect of group explained only 2.7% of the variance in OPAL. When milk yield and production system were included in the same statistical model (n = 28,096 cows), the magnitude of reduction in OPAL as parity increased was markedly reduced, although parity group remained associated with OPAL (Table 8). Intensive production system was associated with lower OPAL (OR = 0.562; Table 8). Figure 8 shows the quadratic association of milk yield with predicted OPAL and the interactions for parity 1 and 2 cows in which the quadratic effects are more marked. For parity 1 cows, the quadratic effect is more centered on the group mean milk than for older groups that have an optimal OPAL above mean milk production. The marked reduction in OPAL with increased parity, independent of production, is also demonstrated in Figure 8. While the heritability of pregnancy to first service (0.023) or days open (0.06) are low (
      • Berry D.P.
      • Wall E.
      • Pryce J.E.
      Genetics and genomics of reproductive performance in dairy and beef cattle.
      ), multitrait genomic prediction for reproduction has improved these measures (
      • Veronese A.
      • Marques O.
      • Moreira R.
      • Belli A.L.
      • Bisinotto R.S.
      • Bilby T.R.
      • Peñagaricano F.
      • Chebel R.C.
      Genomic merit for reproductive traits. I: Estrous characteristics and fertility in Holstein heifers.
      ,
      • Veronese A.
      • Marques O.
      • Peñagaricano F.
      • Bisinotto R.S.
      • Pohler K.G.
      • Bilby T.R.
      • Chebel R.C.
      Genomic merit for reproductive traits. II: Physiological responses of Holstein heifers.
      ;
      • Lima F.S.
      • Silvestre F.T.
      • Penagaricano F.
      • Thatcher W.W.
      Early genomic prediction of daughter pregnancy rate is associated with improved reproductive performance in Holstein dairy cows.
      ), suggesting the potential to improve fertility through genetic selection. The quadratic association between milk yield and OPAL may be confounded by health disorders. Presumably cows with unresolved disease problems or very low production were never inseminated. Cows with milk yield at ∼72 DIM markedly below their group mean probably reflects effects of disease in many cases, or in a smaller proportion of cases, simply the low extremes of phenotypic potential for milk production. The effects of disease on fertility have been well described (
      • Curtis C.R.
      • Erb H.N.
      • Sniffen C.J.
      • Smith R.D.
      • Kronfeld D.S.
      Path analysis of dry period nutrition, postpartum metabolic and reproductive disorders, and mastitis in Holstein cows.
      ;
      • McDougall S.
      Effects of periparturient diseases and conditions on the reproductive performance of New Zealand dairy cows.
      ;
      • Carvalho M.R.
      • Peñagaricano F.
      • Santos J.
      • DeVries T.
      • McBride B.
      • Ribeiro E.
      Long-term effects of postpartum clinical disease on milk production, reproduction, and culling of dairy cows.
      ). Conversely, for very high-producing, presumably healthy cows, the quadratic association may reflect the competitive demands of lactation and reproduction (
      • MacMillan K.L.
      • Lean I.
      • Westwood C.
      The effects of lactation on the fertility of dairy cows.
      ;

      Lean, I., and A. Rabiee. 2006. Quantitative metabolic and epidemiological approaches to the fertility of the dairy cow. Proceedings of the Dairy Cattle Reproductive Council, DCRC:115–131.

      ;
      • Friggens N.C.
      • Disenhaus C.
      • Petit H.V.
      Nutritional sub-fertility in the dairy cow: Towards improved reproductive management through a better biological understanding.
      ) as supported by associations of output of nutrients in milk protein, diet, and estimated duodenal flux of nutrients with reproductive performance (
      • Rodney R.M.
      • Celi P.
      • Scott W.
      • Breinhild K.
      • Santos J.
      • Lean I.
      Effects of nutrition on the fertility of lactating dairy cattle.
      ).
      Figure thumbnail gr8
      Figure 8The association of milk yield (kg/d) centered on group within herd at approximately 70 DIM with the predicted coefficient of odds of becoming pregnant in a lactation (OPAL). The data are from 28,071 cows. The x-axis is milk yield centered by group, and the y-axis is the predicted coefficient of odds of pregnancy at first insemination. Based on the model that contains parity and production system and the interaction of parity with production system, all parities differed, with cows in parity 1 having the greatest odds and parity ≥5 having the least odds of pregnancy in the lactation.
      The results in Table 8 support an association between the metabolic challenge from milk and milk component production and fertility. The univariable model in the sample of cows with milk composition data around peak milk yield had fewer cows (n = 8,757) than the model including milk yield alone (n = 28,096). Although parity 2 did not differ from parity 1, and parity 3 did not differ from parity 4, when milk components were included in the models, the same general trends were observed, with a reduced OPAL as parity number increased, with parity ≥5 having the lowest OPAL (Table 8).
      Milk fat percentage at peak production was also associated quadratically with OPAL in each parity (Table 8). Intensively fed systems had lesser OPAL than pasture fed (OR = 0.40). While these analyses are not directly comparable to the milk yield model with all 28,096 cows, they support the concept that a greater metabolic demand is associated with reduced OPAL. The results for fat yield were similar to milk yield, with a quadratic interaction for parity 1 and a linear interaction for parity 2 (Table 8). Parity ≥5 cows had a lower OPAL than all other parities, and parity 3 and 4 had lower OPAL than parity 2.
      Milk protein percentage was also quadratically associated with OPAL (Table 8; Figure 9). As for fat percentage, the pattern was for lower OPAL with greater parity. Cows in intensively fed systems had lower OPAL than those on pasture (Table 8).
      Figure thumbnail gr9
      Figure 9The association of milk protein percentage centered on group with herd at approximately 70 DIM with the relative predicted odds of becoming pregnant in a lactation (OPAL). The data are from 8,761 cows in 7 studies. The x-axis is milk yield centered by group, and the y-axis is the predicted coefficient of odds of pregnancy at first insemination. Based on the model that contains parity and production system and the interaction of parity with production system, all parity comparisons differed, apart from cows in parity 1 and 2, and cows in parity 3 that did not differ from parity 4.
      For milk protein yield, there were a large number of interactions. All parities had a linear interaction with protein yield, and parity 1 and 2 had quadratic interactions, as with fat and protein percentages. There were quadratic associations between milk protein yield and OPAL, suggesting that an optimum milk protein yield may exist for OPAL (Table 8). Parity differences were present for all except parity 1 with 2 and 3, and for parity 4 with 3 and 5.
      Milk solids yield also was associated quadratically with OPAL (Table 8). In the milk solids and protein yield models, system was not a significant predictor in contrast with the milk and fat yield models. The parity differences were consistent with the other models with milk constituent data, and all parity groups differed except parity 1 from 2 and parity 3 from 4.
      There is a strong consistency of findings evident for the models in Table 8 (Figure 8, Figure 9) that suggest that milk yield or milk protein yield near group means may optimize the odds of pregnancy over a lactation. Interestingly, the absolute values of these group means varied substantially among the studies included. Therefore, these results may reflect nutrition and management that are sufficient for cows close to the herd average, but possibly inadequate to support cows with the propensity to produce well above the group average.

      Overall Discussion

      The increased risk of removal with increased parity was substantial for the older cows. The quadratic associations between milk protein percentage and HPREG, PREG1, and OPAL (Tables 6 to 8; Figure 6) and summarized in Table 10 are similar to
      • Moss N.
      The epidemiology of subfertility in Australian dairy cows.
      ,
      • Madouasse A.
      • Huxley J.
      • Browne W.
      • Bradley A.
      • Dryden I.
      • Green M.
      Use of individual cow milk recording data at the start of lactation to predict the calving to conception interval.
      ,
      • Morton J.M.
      • Auldist M.J.
      • Douglas M.L.
      • Macmillan K.L.
      Milk protein concentration, estimated breeding value for fertility, and reproductive performance in lactating dairy cows.
      , and
      • Rodney R.M.
      • Hall J.K.
      • Westwood C.T.
      • Celi P.
      • Lean I.J.
      Precalving and early lactation factors that predict milk casein and fertility in the transition dairy cow.
      .
      • Rodney R.M.
      • Hall J.K.
      • Westwood C.T.
      • Celi P.
      • Lean I.J.
      Precalving and early lactation factors that predict milk casein and fertility in the transition dairy cow.
      found that HPREG in the first 150 DIM was 28% lower in cows producing milk in the lowest quartile of protein percentage compared with cows in the upper 3 quartiles. Cows in the lowest quartile of milk protein percentage produced 3.87 L/cow per day more milk on average during the first 150 DIM and had greater prepartum plasma calcium and glucose concentrations and estimated MP balance than cows producing a higher milk protein percentage (
      • Rodney R.M.
      • Hall J.K.
      • Westwood C.T.
      • Celi P.
      • Lean I.J.
      Precalving and early lactation factors that predict milk casein and fertility in the transition dairy cow.
      ). Greater protein percentage is an indication that cows are more likely to be reproductively successful; however, the response can be quadratic (Figure 6, Figure 9). Further,
      • Morton J.M.
      • Auldist M.J.
      • Douglas M.L.
      • Macmillan K.L.
      Milk protein concentration, estimated breeding value for fertility, and reproductive performance in lactating dairy cows.
      found that the relationship between milk protein percentage and the odds of pregnancy at first service were influenced by the level of milk yield. That is similar to our findings that the association of milk protein yield with OPAL, HPREG, and PREG1 were consistent with those for milk, milk fat, and milk solids yield, rather than with milk protein percentage. Therefore, yield of milk protein around peak lactation and the early breeding period looks to have a robust association with reproductive outcomes, although it may only explain a small proportion of the total variance in the odds of pregnancy per insemination (
      • Hudson C.D.
      • Green M.J.
      Associations between routinely collected dairy herd improvement data and insemination outcome in UK dairy herds.
      ). Our companion paper (
      • Lean I.J.
      • LeBlanc S.J.
      • Sheedy D.B.
      • Duffield T.
      • Santos J.E.P.
      • Golder H.M.
      Associations of parity with health disorders and blood metabolite concentrations in Holstein cows in different production systems.
      ) highlights many differences in risk of disease and differences in blood metabolites that reduce the risk of survival and may help to explain differences in reproductive success for older cows. These findings suggest a need to emphasize nutritional and other interventions that support older or higher-producing cows and further evaluate reasons for loss in lower-producing cows. Future studies should include accurate and comparable health and environmental data to refine understanding of factors that influence reproductive performance and, in particular, characterize the metabolic and other physiological changes that so dramatically increase risk of reproductive failure with parity.
      Table 10Summary of the associations of milk production, milk fat percentage and yield, milk protein percentage and yield, and yield of milk solids with reproductive outcomes
      I indicates an interaction only, −L indicates a negative linear association, and Q indicates a quadratic association.
      Milk measureHazard of not being bred (HNBRED)
      Probability or hazard of not being bred over time.
      Hazard of pregnancy (HPREG)
      Probability or hazard of pregnancy over time.
      Pregnancy to first breeding (PREG1)Odds of becoming pregnant in a lactation (OPAL)
      Production (kg/d)QQ−LQ
      Fat (%)−L−L−LQ
      Fat yield (kg/d)−L−LI−L
      Protein (%)QQQQ
      Protein yield (kg/d)Q−L−LQ
      TS yield (kg/d)Q−L−LQ
      1 I indicates an interaction only, −L indicates a negative linear association, and Q indicates a quadratic association.
      2 Probability or hazard of not being bred over time.
      3 Probability or hazard of pregnancy over time.

      CONCLUSIONS

      Cows have increased hazards of removal for all reasons and of death and culling with increased parity. This effect of parity is evident in reduced reproductive performance. Cows in intensively fed systems had lesser probability of HNBRED than those in pasture-based systems, but production system had little influence on PREG1. Relationships of all milk production variables, except fat yield, with OPAL were quadratic, and older cows had markedly lower odds of OPAL than parity 1 and 2. Milk yield was quadratically associated with all reproductive measures except PREG1, where the association was linearly negative. As yield of milk increased, OPAL improved up to approximately 10 kg/d greater than the group mean (Figure 8), after which it decreased. However, parity 1 cows had optimal OPAL at approximately group mean milk. In contrast to other milk or milk component outputs, milk protein percentage had a largely positive association with reproductive outcomes, and fat percentage had largely negative associations. The associations of protein yield that were linearly negative for HPREG and PREG1 can be contrasted with the quadratic associations of milk protein percentage with the same variables. Evaluations of the effects of parity on the odds of pregnancy are striking when considered over the whole lactation (OPAL) versus at PREG1 alone; there is much less association for parity with PREG1 than OPAL. The effect on OPAL may not just reflect reproductive failure, but other reasons for loss from the herd. Lower producing cows, for example, had increased HNBRED. Our results provide insights on the complex relationships among parity, measures of production, and measures of reproductive success.

      ACKNOWLEDGMENTS

      Thanks to all that provided their data and C. Ha for assistance in formatting the results. This project was funded by Dairy UP (Brownlow Hill, NSW, Australia), a joint project between University of Sydney (Camden, NSW, Australia), Scibus (Camden, NSW, Australia), and New South Wales Department of Primary Industries (Orange, NSW, Australia) and supported by the NSW Government, Australian Fresh Milk Holding Ltd. (Gooloogong, NSW, Australia), Bega Cheese (Bega, NSW, Australia), Dairy Australia (Southbank, Vic, Australia), Dairy Connect (Mascot, NSW, Australia), DairyNSW (Camden, NSW, Australia), Local Land Services–Hunter (Tocal, NSW, Australia), Leppington Pastoral Co. (Bringelly, NSW, Australia), Norco Dairy Co-Op (South Lismore, NSW, Australia), NSW Farmers (St Leonards, NSW, Australia), NSW Department of Primary Industries (Menangle, NSW, Australia), Scibus, and South East Local Land Services (Goulburn, NSW, Australia). The authors have not stated any conflicts of interest.

      REFERENCES

        • Bascom S.S.
        • Young A.
        A summary of the reasons why farmers cull cows.
        J. Dairy Sci. 1998; 81 (9749397): 2299-2305
        • Bello N.M.
        • Steibel J.
        • Erskine R.
        • Tempelman R.
        Cows and herds constitute distinct hierarchical levels of heterogeneity in the variability of and association between milk yield and pregnancy outcome in dairy cows.
        J. Dairy Sci. 2013; 96 (23462172): 2314-2326
        • Berry D.P.
        • Friggens N.C.
        • Lucy M.
        • Roche J.
        Milk production and fertility in cattle.
        Annu. Rev. Anim. Biosci. 2016; 4 (26526546): 269-290
        • Berry D.P.
        • Wall E.
        • Pryce J.E.
        Genetics and genomics of reproductive performance in dairy and beef cattle.
        Animal. 2014; 8 (24703258): 105-121
        • Canadian Dairy Information Centre (CDIC)
        Culling and replacement rates in dairy herds in Canada.
        • Carvalho M.R.
        • Peñagaricano F.
        • Santos J.
        • DeVries T.
        • McBride B.
        • Ribeiro E.
        Long-term effects of postpartum clinical disease on milk production, reproduction, and culling of dairy cows.
        J. Dairy Sci. 2019; 102 (31548073): 11701-11717
        • Curtis C.R.
        • Erb H.N.
        • Sniffen C.J.
        • Smith R.D.
        • Kronfeld D.S.
        Path analysis of dry period nutrition, postpartum metabolic and reproductive disorders, and mastitis in Holstein cows.
        J. Dairy Sci. 1985; 68 (4067048): 2347-2360
        • Curtis M.A.
        Epidemiology of Uterine Infections in Dairy Cows: Antioxidant and Metabolic Investigations.
        (PhD thesis) Department of Veterinary Clinical Studies, University of Sydney, 1997
        • Curtis M.A.
        • Lean I.J.
        Path analysis of metabolic and antioxidant risk factors for periparturient and postparturient conditions and reproductive performance in dairy cows.
        in: Proc. XX World Buiatrics Congress. Sydney, Australia. 1998: 809-818
        • Dallago G.M.
        • Wade K.M.
        • Cue R.I.
        • McClure J.T.
        • Lacroix R.
        • Pellerin D.
        • Vasseur E.
        Keeping dairy cows for longer: A critical literature review on dairy cow longevity in high milk-producing countries.
        Animals (Basel). 2021; 11 (33805738): 808-833
        • De Vries A.
        • Marcondes M.
        Overview of factors affecting productive lifespan of dairy cows.
        Animal. 2020; 14 (32024570): s155-s164
        • DeGaris P.J.
        • Lean I.
        • Rabiee A.
        • Heuer C.
        Effects of increasing days of exposure to prepartum transition diets on reproduction and health in dairy cows.
        Aust. Vet. J. 2010; 88 (20402690): 84-92
        • Duffield T.
        • Bagg R.
        • DesCoteaux L.
        • Bouchard E.
        • Brodeur M.
        • DuTremblay D.
        • Keefe G.
        • LeBlanc S.
        • Dick P.
        Prepartum monensin for the reduction of energy associated disease in postpartum dairy cows.
        J. Dairy Sci. 2002; 85 (11913700): 397-405
        • Duffield T.F.
        • Leslie K.
        • Sandals D.
        • Lissemore K.
        • McBride B.
        • Lumsden J.
        • Dick P.
        • Bagg R.
        Effect of a monensin-controlled release capsule on cow health and reproductive performance.
        J. Dairy Sci. 1999; 82 (10575604): 2377-2384
        • Duffield T.F.
        • Leslie K.
        • Sandals D.
        • Lissemore K.
        • McBride B.
        • Lumsden J.
        • Dick P.
        • Bagg R.
        Effect of prepartum administration of monensin in a controlled-release capsule on milk production and milk components in early lactation.
        J. Dairy Sci. 1999; 82 (10068948): 272-279
        • Fetrow J.
        • Nordlund K.
        • Norman H.
        Invited review: Culling: Nomenclature, definitions, and recommendations.
        J. Dairy Sci. 2006; 89 (16702253): 1896-1905
        • Friggens N.C.
        • Disenhaus C.
        • Petit H.V.
        Nutritional sub-fertility in the dairy cow: Towards improved reproductive management through a better biological understanding.
        Animal. 2010; 4 (22444617): 1197-1213
        • Golder H.
        • McGrath J.
        • Lean I.
        Supplementary material: Effect of 25-hydroxyvitamin D3 [25-(OH)D3] during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows. Figshare.
        • Golder H.M.
        • Rossow H.
        • Lean I.
        Effects of in-feed enzymes on milk production and components, reproduction, and health in dairy cows.
        J. Dairy Sci. 2019; 102 (31279550): 8011-8026
        • Golder H.M.
        • McGrath J.
        • Lean I.J.
        Effect of 25-hydroxyvitamin D3 during prepartum transition and lactation on production, reproduction, and health of lactating dairy cows.
        J. Dairy Sci. 2021; 104 (33663856): 5345-5374
        • Holden S.A.
        • Butler S.T.
        Review: Applications and benefits of sexed semen in dairy and beef herds.
        Animal. 2018; 12 (29631644): s97-s103
        • Hudson C.D.
        • Green M.J.
        Associations between routinely collected dairy herd improvement data and insemination outcome in UK dairy herds.
        J. Dairy Sci. 2018; 101 (30316603): 11262-11274
        • Lean I.J.
        • Galland J.
        • Scott J.
        Relationships between fertility, peak milk yields and lactational persistency in dairy cows.
        Theriogenology. 1989; 31 (16726627): 1093-1103
        • Lean I.J.
        • LeBlanc S.J.
        • Sheedy D.B.
        • Duffield T.
        • Santos J.E.P.
        • Golder H.M.
        Associations of parity with health disorders and blood metabolite concentrations in Holstein cows in different production systems.
        J. Dairy Sci. 2023; 106: 500-518
      1. Lean, I., and A. Rabiee. 2006. Quantitative metabolic and epidemiological approaches to the fertility of the dairy cow. Proceedings of the Dairy Cattle Reproductive Council, DCRC:115–131.

        • Lean I.J.
        • Sheedy D.B.
        • LeBlanc S.J.
        • Duffield T.
        • Santos J.E.P.
        • Golder H.M.
        Holstein dairy cows lose body condition score and gain body weight with increasing parity in both pasture-based and total mixed ration herds.
        JDS Commun. 2022; 3: 431-435
        • LeBlanc S.
        Assessing the association of the level of milk production with reproductive performance in dairy cattle.
        J. Reprod. Dev. 2010; 56 (20629210): S1-S7
        • LeBlanc S.J.
        • Duffield T.F.
        • Leslie K.E.
        • Bateman K.G.
        • TenHag J.
        • Walton J.S.
        • Johnson W.H.
        The effect of prepartum injection of vitamin E on health in transition dairy cows.
        J. Dairy Sci. 2002; 85 (12146472): 1416-1426
        • Lima F.S.
        • Silvestre F.T.
        • Penagaricano F.
        • Thatcher W.W.
        Early genomic prediction of daughter pregnancy rate is associated with improved reproductive performance in Holstein dairy cows.
        J. Dairy Sci. 2020; 103 (32089311): 3312-3324
        • Lucy M.C.
        Reproductive loss in high-producing dairy cattle: Where will it end?.
        J. Dairy Sci. 2001; 84 (11417685): 1277-1293
        • MacMillan K.L.
        • Lean I.
        • Westwood C.
        The effects of lactation on the fertility of dairy cows.
        Aust. Vet. J. 1996; 73 (8660229): 141-147
        • Madouasse A.
        • Huxley J.
        • Browne W.
        • Bradley A.
        • Dryden I.
        • Green M.
        Use of individual cow milk recording data at the start of lactation to predict the calving to conception interval.
        J. Dairy Sci. 2010; 93 (20855002): 4677-4690
        • McDougall S.
        Effects of periparturient diseases and conditions on the reproductive performance of New Zealand dairy cows.
        N. Z. Vet. J. 2001; 49 (16032164): 60-67
        • Morton J.M.
        • Auldist M.J.
        • Douglas M.L.
        • Macmillan K.L.
        Milk protein concentration, estimated breeding value for fertility, and reproductive performance in lactating dairy cows.
        J. Dairy Sci. 2017; 100 (28478010): 5850-5862
        • Moss N.
        The epidemiology of subfertility in Australian dairy cows.
        (PhD Thesis) Faculty of Veterinary Science, University of Sydney, 2001
        • Pinedo P.
        • Santos J.E.P.
        • Chebel R.C.
        • Galvão K.N.
        • Schuenemann G.M.
        • Bicalho R.C.
        • Gilbert R.O.
        • Rodriguez Zas S.
        • Seabury C.M.
        • Rosa G.
        • Thatcher W.W.
        Early-lactation diseases and fertility in 2 seasons of calving across US dairy herds.
        J. Dairy Sci. 2020; 103 (32896394): 10560-10576
        • Pryce J.E.
        • Nielsen B.L.
        • Veerkamp R.F.
        • Simm G.
        Genotype and feeding system effects and interactions for health and fertility traits in dairy cattle.
        Livest. Prod. Sci. 1999; 57: 193-201
        • Rabe-Hesketh S.
        • Skrondal A.
        Multilevel and Longitudinal Modeling Using Stata.
        STATA Press, 2005
        • Rearte R.
        • LeBlanc S.J.
        • Corva S.G.
        • de la Sota R.L.
        • Lacau-Mengido I.M.
        • Giuliodori M.J.
        Effect of milk production on reproductive performance in dairy herds.
        J. Dairy Sci. 2018; 101 (29803419): 7575-7584
        • Rodney R.M.
        • Celi P.
        • Scott W.
        • Breinhild K.
        • Santos J.
        • Lean I.
        Effects of nutrition on the fertility of lactating dairy cattle.
        J. Dairy Sci. 2018; 101 (29605330): 5115-5133
        • Rodney R.M.
        • Hall J.K.
        • Westwood C.T.
        • Celi P.
        • Lean I.J.
        Precalving and early lactation factors that predict milk casein and fertility in the transition dairy cow.
        J. Dairy Sci. 2016; 99 (27320662): 7554-7567
        • Santos J.E.P.
        • Rutigliano H.M.
        • Sá Filho M.F.
        Risk factors for resumption of postpartum estrous cycles and embryonic survival in lactating dairy cows.
        Anim. Reprod. Sci. 2009; 110 (18295986): 207-221
        • Schuster J.C.
        • Barkema H.W.
        • De Vries A.
        • Kelton D.F.
        • Orsel K.
        Invited review: Academic and applied approach to evaluating longevity in dairy cows.
        J. Dairy Sci. 2020; 103 (33222845): 11008-11024
        • Spalding R.
        • Everett R.
        • Foote R.
        Fertility in New York artificially inseminated Holstein herds in dairy herd improvement.
        J. Dairy Sci. 1975; 58: 718-723
        • StataCorp LLC
        Stata user's guide release 17.
        Stata Press, 2021
        https://www.stata.com/manuals/u.pdf
        Date accessed: October 3, 2022
        • Stevenson M.A.
        • Lean I.
        Descriptive epidemiological study on culling and deaths in eight dairy herds.
        Aust. Vet. J. 1998; 76 (9700403): 482-488
        • United States Department of Agriculture (USDA)
        Health and management practices on US dairy operations 2014, Report 3.
        USDA–APHIS–VS–CEAH–NAHMS, Fort Collins, CO2018: 1-202
        • Veronese A.
        • Marques O.
        • Moreira R.
        • Belli A.L.
        • Bisinotto R.S.
        • Bilby T.R.
        • Peñagaricano F.
        • Chebel R.C.
        Genomic merit for reproductive traits. I: Estrous characteristics and fertility in Holstein heifers.
        J. Dairy Sci. 2019; 102 (31030916): 6624-6638
        • Veronese A.
        • Marques O.
        • Peñagaricano F.
        • Bisinotto R.S.
        • Pohler K.G.
        • Bilby T.R.
        • Chebel R.C.
        Genomic merit for reproductive traits. II: Physiological responses of Holstein heifers.
        J. Dairy Sci. 2019; 102 (31030930): 6639-6648
        • Vieira-Neto A.
        • Duarte G.A.
        • Zimpel R.
        • Thatcher W.W.
        • Santos J.E.P.
        Days in the prepartum group are associated with subsequent performance in Holstein cows.
        J. Dairy Sci. 2021; 104 (33663839): 5964-5978
        • Vieira-Neto A.
        • Duarte G.A.
        • Zimpel R.
        • Thatcher W.W.
        • Santos J.E.P.
        Days in the prepartum group are associated with subsequent performance in Holstein cows.
        J. Dairy Sci. 2021; 104 (33663839): 5964-5978
        • Westwood C.T.
        • Lean I.J.
        • Garvin J.K.
        Factors influencing fertility of Holstein dairy cows: A multivariate description.
        J. Dairy Sci. 2002; 85 (12512596): 3225-3237
        • Westwood C.T.
        • Lean I.J.
        • Garvin J.K.
        • Wynn P.C.
        Effects of genetic merit and varying dietary protein degradability on lactating dairy cows.
        J. Dairy Sci. 2000; 83 (11132865): 2926-2940
        • Workie Z.W.
        • Gibson J.P.
        • van der Werf J.H.J.
        Analysis of culling reasons and age at culling in Australian dairy cattle.
        Anim. Prod. Sci. 2021; 61: 680-689