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

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.


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 et al. (2021) 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 and Marcondes (2020) 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 Increased parity is negatively associated with survival and reproduction in different production systems in 3 herds were parity 1, indicating a substantial loss of older cows (Golder et al., 2019). Percentage of cows culled in Canada from 2015 to 2020 was 30.7% (CDIC, 2021) and the Northeastern United States in 2014 was 31.4% plus a 6.2% death rate (USDA, 2018). The cows of parity ≥ 2 in a United States-based study (Golder et al., 2019) 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 et al., 2021). Given increased use of sexed semen (Holden and Butler, 2018), 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 et al., 2006). However, De Vries and Marcondes (2020) suggest that an average productive life span-that is, from time of calving to removal (Schuster et al., 2020)-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 and Young (1998) 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 and Lean (1998) and Workie et al. (2021) 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 and Lean, 1998).
Associations between increased production and reduced pregnancy per AI or extended days open were described in the 1970s by Spalding et al. (1975) and more recently by many workers (Lucy, 2001;Berry et al., 2016). Flaws in many of the studies evaluating these associations have been detailed and include use of inappropriate measures of reproduction or production (LeBlanc, 2010). 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 pasturebased or intensively fed with a mixed or componentfed 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 478 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 Tables 1 and 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 Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION Table 1. Summary 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   Golder and Lean (unpublished) that the monthly herd tests provide sufficient precision to test our hypotheses. 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 (Tables 1 and 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 et al. (2021).

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 et al. (2020). 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 Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION 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, 2021) states that in a Cox model, the data are organized as i = 1,… n groups with j = 1,… n i observations in group i. For the jth observation in the ith group, the hazard is where h ij (t) is the hazard at time t for the jth observation in the ith group and h 0 (t) is the baseline hazard, α i is the group-level frailty or random effect, β is the coefficient for x ij , the row vector of covariates for the time interval (t 0ij ; t ij ], where t 0ij is the initial time of observation and t ij 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, 2021) states that the model represents a series of M indepen-dent clusters, and is conditional on a set of random effects u j , where j = 1, … M clusters, with cluster j consisting of i = 1, … n j observations. The responses are the binary valued y ij , with the convention of treating y ij = 1 if depvar ij (the dependent variable) ≠ 0, and treating y ij = 0 otherwise. Pr is the probability for (y ij = 1|x ij , u j ). The 1 × p row vector (x ij ) 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 z ij 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 y ij = 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 and Skrondal, 2005) 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 et al., 2022(Lean et al., , 2023 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 et al., 2014)-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).
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 et al., 2023). 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 et al. (2018) 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).
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 et al., 1989;Santos et al., 2009;LeBlanc, 2010;Rearte et al., 2018). 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 Tables 6 and 9 Table 3. Associations 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 Measure 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 Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION Figure 1. Kaplan-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. 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.

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

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

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. 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 et al., 2002). 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 Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION Figure 2. Kaplan-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.

490
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 et al., 2002).
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 (2010) 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 et al. (2018) 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.
Other studies found differing and conditional associations of milk production with reproduction. LeBlanc (2010) similarly found lower reproductive performance Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION Figure 3. Kaplan-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.
in cows with lower milk production in early lactation but did not detect a quadratic relationship of milk yield with time to pregnancy. Rearte et al. (2018) also found that higher-yielding cows had shorter time to pregnancy but only in lower-yielding herds. In higheryielding 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 et al., 2010), 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 et al. (2013) 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 reproduc-tive 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 per- Lean et al.: PARITY,SURVIVAL,AND REPRODUCTION Figure 4. The 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. centage 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.
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 et al. (1999) found first service conception was greater for parity 2 to 4 cows than those in first or parity greater than 4. Lean et al. (1989) 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 et al. (1989), Rearte et al. (2018) assessed the odds of pregnancy within 100 DIM and found that par- Lean et al.: PARITY,SURVIVAL,AND REPRODUCTION Figure 5. The 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. ity 2 cows did not differ from parity 1 and that parity 3 and ≥4 had reduced odds.
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 et al. (2018) 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 and Green (2018) 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 Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION Figure 6. The 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. 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 et al., 2014), multitrait genomic prediction for reproduction has improved these measures (Veronese et al., 2019a,b;Lima et al., 2020), 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 produc- Lean et al.: PARITY,SURVIVAL,AND REPRODUCTION Figure 7. The 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.
tion. The effects of disease on fertility have been well described (Curtis et al., 1985;McDougall, 2001;Carvalho et al., 2019). Conversely, for very high-producing, presumably healthy cows, the quadratic association may reflect the competitive demands of lactation and reproduction (Macmillan et al., 1996;Lean and Rabiee, 2006;Friggens et al., 2010) as supported by associations of output of nutrients in milk protein, diet, and estimated duodenal flux of nutrients with reproductive performance (Rodney et al., 2018).
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).
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. Lean et al.: PARITY,SURVIVAL,AND REPRODUCTION Figure 8. The 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.
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 (Figures 8 and 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 sum-marized in Table 10 are similar to Moss (2001), Madouasse et al. (2010), Morton et al. (2017), and Rodney et al. (2016). Rodney et al. (2016) 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 et al., 2016). Greater protein percentage is an indication that cows are more likely to be reproductively successful; however, the response can be quadratic (Figures 6 and 9). Further, Morton et al. (2017) 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 as- Lean et al.: PARITY, SURVIVAL, AND REPRODUCTION Figure 9. The 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. sociation with reproductive outcomes, although it may only explain a small proportion of the total variance in the odds of pregnancy per insemination (Hudson and Green, 2018). Our companion paper (Lean et al., 2023) 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.

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.