Evaluation of Non-linear Models to Predict Potential Milk Yield of Beef Cows According to Parity Order Under Grazing

This study aimed to evaluate the effect of parity order on milk yield (MY) and composition over time of grazing beef cows and to evaluate non-linear models to describe the lactation curve. Thirty-six pregnant Nellore cows (12 nulliparous, 2 years; 12 primiparous, 3 years; and 12 multiparous, 4–6 years) were included in the study. With calving day assigned as day 0, milking was performed using a milking machine to estimate MY on days 7, 14, 21, 42, 63, 91, 119, 154, and 203. Dummy variable analyses were applied to estimate its effects on MY, composition (kg and percentage), afternoon/morning, and afternoon/total proportions. Since multiparous cows had higher MY than nulliparous and primiparous cows, two different groups were used for lactation curve analysis: Mult (multiparous) and Null/Prim (nulliparous and primiparous). The MY estimated by the last edition of BR-Corte (Nutrient Requirements of Zebu and Crossbred Cattle) equation was compared with the observed values from this study. Five nonlinear models proposed by Wood (WD), Jenkins & Ferrell (JF), Wilmink (WK), Henriques (HR) and Cobby & Le Du (CL) were evaluated. Models were validated using an independent dataset of multiparous and primiparous cows. The estimates for parameters a, b, and c of the CL equation were compared between groups, and the BR-Corte equation used the model identity methodology. Nulliparous and primiparous cows displayed similar MY (P > 0.05); however, multiparous cows had an average MY that is 0.70 kg/day greater than that of nulliparous and primiparous cows (P < 0.05). Milk protein and total solids were higher for multiparous cows (P < 0.05). Effect of days in milking was found for milk fat, protein, and total solids (P < 0.05). The yield of all milk components was higher for multiparous cows than for nulliparous and primiparous cows. The afternoon/morning and afternoon/total proportions of milk production were not affected by parities and days in milking (P > 0.05), with an average of 0.76 and 0.42, respectively. The BR-Corte equation did not correctly estimate the MY (P < 0.05). The equations of WD, WK, and CL had the best estimate of MY for both Mult and Null/Prim datasets. The equations had a very similar Akaike's information criterion with correction and mean square error of prediction.


INTRODUCTION
The milking ability of beef cows is one of the main factors influencing the weaning weights of calves (1); thus, many methods have been used to attempt to estimate beef cows' milk production and its influence on calf preweaning growth (2)(3)(4). The evaluation of milk production is also necessary to estimate the nutrient requirements of cows and calves since the nutritional balance of a lactating cow is important in the calculation of energy and protein requirements (5,6). Thus, if it is overestimated or underestimated, systematic errors can impair the estimation of nutrient requirements.
Milk yield can be estimated by different methods including determining differences in calf weights before and after suckling (2,7,8) and hand milking (3,4) or machine milking procedures (9)(10)(11). Those methods are possible with small numbers of animals but are not adaptable to larger herds. Such procedures, particularly those with a controlled suckling period before calf separation as described by Boggs et al. (1), require repeated, intensive animal handling in which timing is critical. Furthermore, in general, data collection and handling procedures of grazing animals, specifically pregnant Nellore cows, are laborious since animals have a poor temperament and require extreme care. This challenge would explain the low number of studies and the low number of data points to fit an equation of milk production for grazing beef cows.
In the last BR-Corte edition [third edition of the Brazilian tables of Nutrient Requirements of Zebu and Crossbred Cattle; (12)], the estimated milk yield for Nellore cows during 7 months of lactation was based on the model proposed by Cobby and Le Du (CL) (13). The equation was generated using data from feedlot cows, but the authors employed independent data from experiments with grazing Nellore cows to validate the equation. Although the equation provided a better estimate of milk yield than did previous editions (14), none of these experiments explicitly evaluated the milk yield of different parity orders under grazing conditions, and there are still limited data in this regard.
Previous studies have shown that age and parity order can influence the metabolism and milk production of dairy (15,16) and beef cows (6,17), where primiparous cows display lower milk production and a more unbalanced nutritional status compared to multiparous cows (15,17,18). Thus, different equations are needed to predict the milk yield of Nellore cows according to parity order, since they may differ in milk potential.
Therefore, this study aimed to evaluate the effect of parity order on milk yield and composition over time of grazing beef cows and to evaluate non-linear models to describe the lactation curve.

MATERIALS AND METHODS
The Animal Care and Use Committee at the Universidade Federal de Viçosa (UFV), Brazil (protocol CEUAP-UFV 120/2018), approved all animal care and handling procedures. Animals used in this study were provided by the UFV/Beef Cattle Research and Extension Unit, Viçosa-MG, Brazil, where the study was conducted from July 2018 to May 2019.
The nomenclature for each category related to parity was set at the beginning of the experiment (late gestation period) and used throughout the manuscript. Even though after calving the parity order changed (e.g., nulliparous cows became primiparous), the nomenclature remained constant along the manuscript.
Parity classes were systematically randomized into six paddocks, assuming two cows from each parity class in the paddocks (thus characterizing sub-blocks). The animals were assigned to paddocks 15 days before the beginning of the experiment to acclimate them to the environment and the herd. The average area of the paddocks was 7 ha, covered with Urochloa decumbens grass, and cows had free access to water and feeders.
All cows were group-fed with an energy-protein supplement (1.0 kg/day) with 35% crude protein (CP) for 60 days prepartum (41.2% corn meal, 56.3% soybean meal, and 2.5% urea: ammonium sulfate). The supplement was calculated to supply approximately 40% of the cows' protein requirements, as recommended by the BR-Corte (12). We provided a linear trough space of 0.70 m per cow to ensure homogeneous supplement consumption among groups. The supplement was supplied at 12:00 pm to minimize any interference from animal grazing behavior.
After 90 days of age, calves were offered 5 g/kg BW of an energy-protein supplement formulated to contain 20% CP in a creep-feeding system until the end of the study (203 days in milk).

Milk Collection and Analyses
Calving day was assigned as day 0, and milking was performed using a milking machine to estimate milk yield on days 7, 14, 21, 42, 63, 91, 119, 154, and 203.
Milking procedures were as described by Boggs et al. (1), providing a controlled suckling period before the calf separation. To empty udders, calves were separated from their dams at 3:00 pm and then reunited at 5:45 pm and allowed to suckle until 6:00 pm, when they were once again separated. The first milking was performed at 6:00 am on the next day after an injection of 10 IU (international unit) of oxytocin (10 IU/ml; Ocitovet R , Brazil) in the cow's mammary vein, and the produced milk was weighed. The exact time when the milking of each cow ended was recorded. After morning milking, cows were kept separated from their calves until the end of the afternoon milking, which was performed at 6:00 pm. Then, the total daily production was calculated by the sum of both milking times to obtain a 24-h milk production.
From each cow, a 30-ml sample of milk was collected at morning and afternoon milking to evaluate milk composition. Samples were stored at 4 • C in a refrigerator, each receiving one bronopol tablet per sample as a preservative for further analyses. Milk samples were analyzed fresh for percentage of protein, fat, lactose, and total solids content using infrared spectroscopy (Foss MilkoScan FT120, São Paulo, Brazil). A weighted average was calculated for each component based on morning and afternoon milk yields.

Forage Sampling and Analyses
Every 30 days, grass samples were collected by hand-plucking to evaluate the forage selected by animals and collected by cutting at the ground level from five delimited areas (0.5 × 0.5 m), selected randomly in each paddock to quantify dry matter (DM) per hectare. In these circumstances, all the samples were weighed, oven-dried (55 • C), and then ground to pass through 1-and 2-mm screens in a Wiley mill (model 3, Arthur H. Thomas, Philadelphia, USA). All data from each month were combined and expressed as an average per season as follows: dry season = August (beginning of the experiment), dry-rainy transition = September to November; rainy season = December to February; and rainy-dry transition = March to June (end of the experiment).
The forage and supplement samples were analyzed following the procedures described by the Brazilian National Institute of Science and Technology in Animal Science [INCT-CA; (19)] for DM (method G-003/1), ash (method M-001/1), CP (method N-001/1), and neutral detergent fiber corrected for ash and protein (apNDF; method F-002/1). Indigestible neutral detergent fiber [iNDF; (20)] was processed at 2 mm and quantified by in situ incubation procedures with non-woven textile bags (100 g/m 2 ) for 288 h.

Statistical Analyses
Variables measured during lactation were analyzed using the following model: where Y ijkl = observation measured on animal i, belonging to parity class j, within paddock k and sub-block l; µ = overall mean; C j = parity effect j (fixed); P k(l) = paddock effect k with the nested sub-block l (random); and e ijkl = residual random effect, assumed to be independent and identically distributed (0, σ 2 e ). The covariate day in milk was added to this model when its effect was significant. Given the ordinal nature of parity class (0, 1, and 2), dummy variable analyses (21) were applied to estimate its effects on milk yield, composition, and afternoon/morning and afternoon/total proportions. For this, the mixed model presented in [1] was fitted by using PROC MIXED of SAS R software, assuming the dummy variable as covariates. The parity effect estimates were added to the overall mean to report the results on the same scale as the observed data.
Based on model [1], since multiparous cows had higher milk yield than nulliparous and primiparous cows, two different groups were used for lactation curve analysis: Mult (multiparous cows) and Null/Prim (both nulliparous and primiparous categories together).
The milk yield estimated by the BR-Corte (12) equation was compared with the observed values from this study using the Model Evaluation System (MES version 3.2.2) using the following regression model: where x = predicted values; y = observed values; and β 0 and β 1 = intercept and slope, respectively.
The regression was evaluated according to the following statistical hypothesis: Estimates were evaluated using the estimated value of the mean square error of the prediction and its components (22): where x is the predicted values; y is the observed values; MSEP is the mean squared error of prediction; SB is the squared bias; MaF is the magnitude of random fluctuation; MoF is the model random fluctuation; S x and S y are the standard deviations of the predicted and observed values, respectively; and r is the Pearson linear correlation between the predicted and observed values. For all calculations of variance and covariance, the total number of observations was used as a divisor because it was an estimate of the prediction error (22). The prediction of efficiency was determined by estimating the correlation and concordance coefficient (CCC) or reproducibility index described by Tedeschi (23).
The CCC indicates models with good accuracy and precision (when close to 1.0) or models with a problem of reproducibility (when close to 0.0). The smallest mean square error of prediction indicates the best model in the evaluation. In this study, it can indicate that the model error is associated with the squared bias (SB) or errors related to the high dispersion of data around the mean (MaF) or systematic errors concerning the direction of the curve predicted (MoF).
Five non-linear different equation forms were fitted to the datasets Mult and Null/Prim: [CL] Cobby and Le Du (13)   where Y t = milk yield (kg/day) from lactating beef cows; t = time of lactation (weeks), exp = exponential of natural logarithm; a = theoretical initial yield (production scale) (WD, CL, WK, JF, and HR); b = decrease rate of production after peak of lactation (CL) and increase rate of production up to peak of lactation (WD, HR, and WK); c = increase rate of production up to peak of lactation (CL) and decrease rate of production after peak of lactation (WD, JF, and HR); d = decrease rate of production after peak of lactation (WK); and e t = residual term, assumed as NIID (0, σ 2 e ).
WD (24) is the most widely applied equation for dairy cattle lactation curve. JF (26) proposed a similar equation to WD but without the parameter b. The HR (27) equation is based on the JF model with the addition of a parameter for adjustment for the beginning of lactation. WK (25) is a modification of the lactation curve function in CL.
The parameters of the functions were estimated through the procedure NLMIXED of SAS (version 9.3, SAS Inst. Inc., Cary, NC) and were adjusted by the Gauss-Newton method.
An independent database (Rodrigues, 2021, unpublished data) of Nellore cows from the UFV/Beef Cattle Research and Extension Unit was used for training the models, which were composed of 130 and 112 observations of multiparous (n = 16) and primiparous (n = 14) cows, respectively, with eight lactation points each (1,2,4,6,9,13,20, and 29 weeks). The milking procedures were performed using the exact same approach as the present study: milking twice (morning and afternoon) with a controlled suckling period before the calf separation. The predicted values of the alternative equations generated were compared to the external data values of milk yield using MES (version 3.2.2). Akaike's information criterion with correction (AICC), MSEP, and the distribution of the error (SB, MaF, and MoF) were used to choose the best models.
Milk production upon lactation peak and time until peak were calculated for a 30-week lactation length by the derivative of the equations equal to zero.
The estimates for the parameters a, b, and c from the CL (13) model were compared between groups and the BR-Corte equation (12) (current equation to estimate milk yield of Nellore cows) using model identity methodology based on the overlapping of asymptotic confidence intervals (21).
All statistical evaluations were performed considering 0.05 as the critical level of probability for the occurrence of type I error.

RESULTS
The average DM yield was expressed per period (season) during the experiment as follows: dry season = 4.69 t/ha; dry-rainy transition = 4.33 t/ha; rainy season = 2.93 t/ha; and rainy-dry transition = 3.74 t/ha (Figure 1). Forage chemical composition by season is presented in Table 1.
Milk yield was different between parities (P < 0.05). Nulliparous and primiparous cows displayed similar milk yield (P > 0.05). However, multiparous cows had an average milk yield that is 0.70 kg/day greater than that of nulliparous and primiparous cows. The effect of days in milking was significant (P < 0.05) and was estimated as −0.008 kg/day, which is a decrease in production per day ( Table 2).
Milk fat did not differ between parities (P > 0.05); however, there was an effect of days in milking, in which milk fat increased linearly by 0.001% per day (P < 0.05; Table 2 and Figure 2A).
Milk lactose was on average 4.52% and did not differ between parities or days in milking (P > 0.05; Table 2 and Figure 2C).
Milk protein and total solids were different between parities (P < 0.05; Table 2), being higher for multiparous cows (0.11 and 0.31%, respectively). Days in milk also affected these variables, as they increased linearly by 0.001 and 0.003% per day, respectively (Figures 2B,D).
However, the yield of milk components expressed in kilogram (fat, protein, lactose, and total solids yield) was all higher for multiparous cows (P < 0.05), while nulliparous and primiparous cows displayed similar yields (P > 0.05; Table 2).
The intercepts for fat, protein, lactose, and total solids yield were 0.367, 0.242, 0.328, and 0.932 for multiparous cows and 0.317, 0.283, 0.242, and 0.866 for nulliparous and primiparous cows, respectively, decreasing linearly by 0.0003, 0.0002, 0.0004, and 0.0009 per day, respectively (P < 0.05, Figure 3). The afternoon/morning and afternoon/total proportions of milk production were not affected by parity or days of milking (P < 0.05), with an average of 0.76 and 0.42, respectively ( Table 2). Based on the milk yield differences between parities, two different groups were used for lactation curve analysis: Mult (multiparous cows) and Null/Prim (both nulliparous and primiparous categories together).
The BR-Corte equation did not estimate the milk yield of the dataset correctly because the intercept and slope did not differ from 0 and 1, respectively (P < 0.05). The BR-Corte equation demonstrated overprediction, with 53 and 67% of the MSEP related to model bias Table 3, which shows the need for the development of new equations to estimate the milk yield of grazing Nellore cows.
Therefore, five non-linear alternative equations forms were fitted to the datasets Mult and Null/Prim: CL, WD, WK, JF, and HR (Figures 4, 5).
The equation JF fitted for Mult and Null/Prim did not estimate the milk yield of the independent data correctly because the intercept and slope did not differ from 0 and 1, respectively. In addition, both equations had the highest AICC and MSEP, showing that the models did not adjust to the independent data (Tables 4, 5).
The HR equation did not estimate the milk yield of the independent data for the Null/Prim equation because the slope did not differ from 1 (P < 0.05; Table 5). Even though the HR equation fitted well to the Mult data, we did not consider it as a plausible equation form for use, as it did not adjust also to the Null/Prim dataset.
The equations of WD, WK, and CL had the best estimate of milk yield and description of the lactation curve of grazing Nellore cows for both Mult and Null/Prim datasets (Tables 4, 5). The equations had very similar AICC and MSEP, with 96% for Mult and around 92 to 96% for Null/Prim of the MSEP associated with random error (MoF). All three equations were similar in estimate mean milk yield; hence, the choice of any of those would not impact negatively the estimate of milk yield.
The estimated milk yield at the peak of production and time until peak for Mult and Null/Prim according to equations form are as follows:     (Figure 6).

DISCUSSION
Nulliparous and primiparous cows are in a different physiological state than multiparous cows since they are not physically mature at this age (28). Therefore, a difference in milk production was expected; however, the magnitude of the differences was not known. Our outcome reveals an average milk yield that was approximately 11% lower (−0.70 kg/day) for nulliparous and primiparous cows than for multiparous Nellore cows. These results are consistent with other studies evaluating Bos taurus dairy (15,16) and beef cows (6,29). However, the range of the differences in milk yield according to age or parity, as well as the lactation peak and persistence, may vary between species (such as B. taurus and Bos indicus) and breed crosses.
In the Nutrient Requirements of Beef Cattle, NASEM (2016), the equation used to predict milk yield of beef cows includes an age coefficient that represents 26 and 12% less milk yield for cows of ≤2 years and >2 but ≤3 years, respectively (6). Also, according to NASEM, the peak lactation of beef cows (based on a wide variety of B. taurus breed and breed crosses) occurs approximately at 8.5 weeks (6). In lactating dairy cows, milk production usually peaks at 4 to 8 weeks postpartum (16). In contrast to both studies, for Nellore cows, we found that the peak of lactation is substantially earlier (around 3 to 6 weeks) than that observed for B. taurus dairy and beef cows. There is no mention of differences in the milk yield and time to peak of lactation according to the parity or age of Nellore cows in the BR-Corte editions (12,14).
Greater milk production of multiparous cows may be explained by the largest balance of net energy for milk production, since this category of nutritional requirements is only for maintenance and milk production, unlike young cows that still have nutritional growth requirements. Also, as nulliparous cows still require nutrients for their continued growth, the mammary gland is not completely developed, and the capacity of milk production is reduced compared to multiparous cows. In dairy cows, although the underlying mechanisms are not well understood, several indicators suggest that the mammary gland is more metabolically active in multiparous cows (i.e., higher expression of genes related to metabolic activity) than in primiparous cows, especially at the onset and peak of lactation (30). This observation suggests, at least in part, that the lower milk production observed in primiparous cows could be related to a lower density of secretory cells. Generally, multiparous animals have a higher milk yield but lower lactation persistency than younger animals (30)(31)(32). In this experiment, besides the lower milk yield, the lactation curve of nulliparous and primiparous cows was indeed much flatter than that of multiparous cows.
The average milk composition of Nellore cows in BR-Corte (12) was 15.0% total solids, 4.58% lactose, and 5.61% fat during all lactation periods. However, milk protein increased during the lactation period from 3.57% (4th week) to 3.97% (28th week). In the present study, all milk components changed during the lactation period except for lactose. Total solids, protein, and fat increased linearly throughout lactation. However, when milk composition is expressed in kilograms per day, all components' yield decreased throughout lactation because the percentage of milk yield reduction was greater than the solids increase. The greatest fat percentage value observed in BR-Corte (12) may be related to a higher-energy diet since the cows were maintained in a feedlot system and fed corn silage, which also impacted the milk total solids values. Milk composition can vary according to breed and diet; however, in general, data from previous studies show a percentage of milk fat ranging from 4.5 to 5% for Nellore cows (10,11,33).
Reduced milk protein and total solids for nulliparous cows were expected. They may be explained by the differences in protein metabolic status since nulliparous cows have lower blood albumin and total protein than multiparous beef (17) and dairy cows (34,35). These differences are even higher when milk    components were expressed in kilograms since parities also differ in milk yield. Over the years, several models and approaches have been used to describe the lactation curve of cattle (13, 24-27, 36, 37). The equation most widely used for describing the lactation curve of dairy cows is the model proposed by WD (24); however, its use for beef cattle milk production has been limited since it demands a relatively large number of data points to fit the equation. Therefore, NASEM used a similar equation proposed by JF (26) that requires fewer data points (6).
In the study of HR (27), the best equations to estimate the milk yield of Nellore cows were based on WD (24) and a modification of the JF (26) equation. However, the method of milk collection was weigh-suckle-weigh, then multiplying the morning yield by 2 to estimate daily milk yield, which has been criticized due to data variation and wrong daily estimates. The model that fits better for describing a lactation curve for Nellore cows, according to the BR-Corte (12), is the one proposed by CL (13) in which, after peak, milk yield tends to decline linearly. Although the equation presented in BR-Corte (12) provides a better estimation than the past editions, there are still limitations on its database.
Studies regarding milk production of Nellore cows are often conducted in pens and using diets that mischaracterize the range cattle system, as well as methods of milking that can overestimate or underestimate the milk production depending on the approach. Therefore, the BR-Corte model overprediction was expected because the equation was based on a study developed in a feedlot system using corn silage as a roughage source (38). Although the experiment validated the equation with independent data of studies with grazing Nellore cows, in all of these studies (including the observed values), the cows were milked once in the morning, and then this value was doubled to estimate daily production or extrapolated the morning milk yield for a 24-h period. We believe that this overestimation is also related to the equivocated estimation of milk yield based on only one milking because the morning and afternoon yields contribute to the total daily production in different ratios (10,39). Based on this rationale, the development of new equations that correctly estimate milk production on a grazing range system was necessary. Those explanations are supported by the analysis of model identity that showed that the parameters of the CL equations created compared with the BR-Corte were indeed different, even for the multiparous equation, which used the same animal category as the BR-Corte database.
The equations of WD, WK, and CL were the best in estimating milk yield and describing the lactation curve of both Mult and Null/Prim datasets. All three equations were similar in estimates of mean milk yield and milk production upon lactation peak. The only difference between the equations is the occurrence peak of production, in which the CL equation estimates an earlier peak of production than the WD and WK equation forms. In terms of using milk yield information for energy and protein requirements, choosing any of those would not negatively impact the data estimates.
In this experiment, morning and afternoon yields indeed contributed at different ratios to the total daily production. However, the values are not consistent with those of other studies (10,39). In addition, we did not find any effect of days in milking for these ratios. Thus, it can be assumed that only one ratio for the whole lactation is needed to estimate total daily milk production, based only on the morning milking.

CONCLUSIONS
Parity influences the milk production and composition of beef cows. Milk yield can be estimated based on WD, WK, and CL equation forms, using different equations for multiparous (4-6 years) and nulliparous/primiparous (2-3 years) cows, with an estimated average peak of milk production of 7.5 and 6.3 kg, respectively, and time until peak ranging from 3 to 6 weeks.

DATA AVAILABILITY STATEMENT
The raw data supporting the conclusions of this article will be made available by the authors, without undue reservation.

ETHICS STATEMENT
The animal study was reviewed and approved by the Animal Care and Use Committee at the Universidade Federal de Viçosa, Brazil (protocol CEUAP-UFV 120/2018) approved all animal care and handling procedures.

AUTHOR CONTRIBUTIONS
MF conceived the study, carried out the experimental trial, performed the statistical analysis and chemical analyses, and wrote the manuscript. LR contributed to draft the manuscript and coordinate the research group. IR carried out the experimental trial and contributed to the draft of the manuscript. SF contributed to data interpretation and to draft the manuscript. FS, ED, and MP contributed to designing the experiment, statistical analysis, and drafting the manuscript. All authors read and approved the final manuscript.