https://doi.org/10.22319/rmcp.v15i3.6528
Article
Modeling lactation curves for milk production, fat and protein, and evaluation of factors that affect them in Holstein cattle in Mexico
Luis Enrique Trejo-Díaz a
Felipe De Jesús Ruiz-López b
Hugo Oswaldo Toledo-Alvarado a
Marina Durán-Aguilar c
Adriana García-Ruiz b
a Universidad Nacional Autónoma de México. Posgrado en Ciencias de la Producción y de la Salud Animal. Circuito de Posgrados, Edificio B, 1er Piso, Ciudad Universitaria, Alcaldía Coyoacán. 04510, Ciudad de México, México.
b Instituto Nacional de Investigaciones Forestales, Agrícolas y Pecuarias. Centro Nacional de Investigación Disciplinaria en Fisiología y Mejoramiento Animal. Querétaro, México.
c Universidad Autónoma de Querétaro. Facultad de Ciencias Naturales, Querétaro, México.
*Corresponding author: garcia.adriana@inifap.gob.mx
Abstract:
The knowledge and modeling of lactation curves make it possible to identify factors that help explain environmental and genetic variations that allow the implementation of a selection program. This work aimed to evaluate different models for milk production, fat, and protein curves in Holstein cattle in Mexico and some factors that affect them. The information used was from 125,982 lactations belonging to 68,804 animals born from 2000 to 2020. The effect of calving number, season of the year, and herd was evaluated. R’s Lactcurves package was employed to fit the 38 models included in the package, of which the best four (Wood, Wilmink, Ali & Schaeffer, and modified Pollot) were chosen and then used to model the individual curves through a nonlinear regression model. The parameters calculated for each model were statistically different among the number of lactations (P<0.05), as well as the number of calving, calving season, and herd (P<0.01). The modeled curves have similar shapes to those reported in other studies, except those obtained for protein in the third and fourth or more calvings with the modified Pollot model. The equation proposed by Wilmink was the one that presented the best fit for the study population according to the different evaluation criteria. Using the model that best suits the data will give a closer predictions to reality, and it can be applied to different areas, such as genetic improvement.
Keywords: Lactation curves, Milk production, Fat, Protein, Holstein.
Received: 11/07/2023
Accepted: 13/05/2024
Introduction
The lactation curve, defined as the graphical representation of milk production during the production cycle, can be described through mathematical functions explaining a biological production process subject to genetic and environmental influences(1,2). Proper modeling of lactation curves allows for a good forecast of total production from partial samples, herd planning based on reliable production prediction, and animal selection through knowledge of the different parts of the curve. Therefore, it is essential to find the mathematical function that best describes the lactation curve of animals in each system of production(2,3).
The lactation curve is usually analyzed through four consecutive sections: a) Initial production, estimated by the average production during days 4 to 6 after the colostrum period, b) Ascending or increasing production phase, which is the rate of ascent, until reaching the maximum level of production, c) Maximum point or peak of production, determined by the highest level of production that the cow reaches within the first 90 days of lactation, and d) Decline or reduction in production, also called persistence, which refers to the decrease in milk secretion from peak production(4).
The use of mathematical models has made it possible to know the lactation curves in different dairy production systems. However, not all populations and production systems adjust to a typical lactation curve, with its different parameters and phases, such as start, ascending phase, peak, and decrease. The parameters of a model that fit the lactation curve must reflect various factors, such as genetic, physiological, productive, and environmental factors, and the interactions between them(2,3). Therefore, it is possible to generate as many curves as there are lactations and sources of variation. Hence, it is essential to know the standard levels of milk production by groups of animals with similar characteristics, such as the same lactation stage, calving season, production level, or lactation number(5).
Nonlinear models to represent lactation curves were initially proposed by Wood and have been used in cattle, sheep, goats, buffaloes, and South American camelids(2). The different mathematical models proposed have presented advantages in the specific modeling of sections of the lactation curve, or they fit correctly to various production systems. For example, Wood’s model fits milk production data well, better predicts actual data during early and late lactation, and less accurately predicts data during middle lactation than other nonlinear models(2,6). Wilmink’s model is also widely used to describe lactation curves in dairy cattle, mainly used to detect environmental effects; however, it has been reported that in some populations, this model tends to underestimate the middle part of the curve and overestimate the final part. The Ali-Schaeffer model fits well for lactations that start with low production and peak earlier than usual(7).
One of the main problems with empirical models is that it has been difficult to give physiological meaning to the parameters derived from them. Several modifications have been made to some models in order to have an interpretation closer to the physiological aspects of the lactation curve, such as those proposed by Pollot(8), where the resulting parameters have a biological interpretation, based on changes in the number of cells in the mammary gland during gestation, lactation, and involution, and their effects on milk production(9).
This work aimed to evaluate different mathematical models and some factors that affect the lactation curves of milk production and its components (fat and protein) in a population of Holstein cattle in Mexico.
Material and methods
Data editing and description
The study included information on milk production in kilograms and fat and protein percentages from 68,804 Holstein cows born from 2000 to 2020, belonging to 198 herds of the intensive production system. The data comes from 17 states of the country: Aguascalientes, Baja California, Coahuila, Chihuahua, Durango, Guanajuato, Hidalgo, Jalisco, State of Mexico, Michoacán, Nayarit, Puebla, Querétaro, San Luis Potosí, Tlaxcala, Veracruz, and Zacatecas, where temperate (central zone) and semi-desert climates (northern zone) usually predominate. Querétaro, Guanajuato, Chihuahua, and the State of Mexico concentrate most information. The Holstein Association of Mexico provided the data. The database excluded lactations that did not have production weighing in the first 30 d, those greater than 500 d, and those that had double or triple peak production since it does not correspond to a standard production curve. Each lactation had information from 4 to 12 weightings, and lactations that had fewer than four useful weightings were eliminated.
The milk days of each weighing were adjusted to minimum and maximum values from 5 to 305 d. When the record was outside this range, it was not included in the analysis. Milk production in kilograms, and fat and protein in percentage were adjusted to the mean ±3 standard deviations. When no fat or protein information was available, information on both components was removed. To define the calving season variable, the animals were grouped into three categories according to the month in which they calved, which correspond to cold, hot, and rainy seasons, respectively. The first group covers from November to February, the second from March to June, and the third from July to October.
After the edition, the study included information from 68,804 Holstein cows, with information from 125,982 lactations (72,979 belonging to the first lactation, 31,371 to the second, 11,922 to the third, and 9,710 to 4 or more lactations), and there were 1,319,810 weightings in total.
RStudio(10) was used to evaluate different mathematical models to describe the representation of lactation curves. A total of 38 models included in the R(11) Lactcurves package were fitted, and the four best models were chosen according to the following selection criteria (Table 1): residual standard error (RSE), coefficient of determination (R2), adjusted coefficient of determination (R2 adjus), log-likelihood (LogL), Akaike information criterion (AIC), corrected Akaike information criterion (CAIC), Bayesian information criterion (BIC), and Durbin-Watson coefficient (DW).
The best models were adjusted to lactations per animal by means of a nonlinear (NLIN) regression model using the Statistic Analysis System(12) program. The parameters that describe the curve, persistence, days to the peak, and peak yield were obtained from each curve.
In addition, through the process of generalized linear models (PROC GLM) in SAS(12), it was evaluated whether, in each model, the number of calving, the herd, and the calving season were statistically important in milk production, with the intention of evaluating parameters that could be incorporated into the prediction model. The ggplot package of R was used to plot the curves by lactation.
The Wood(2) model used was as follows:
Where: yt= milk yield at t days in kg, a= initial yield, b= phase of increase in the curve, c= phase of decline in the curve, and t= days.
From the calculated parameters, it is possible to estimate the days to the peak (), maximum yield at the peak (), and persistence ().
Wilmink’s (13) model is described as:
Where: yt= milk yield at t days in kg, a= initial yield, b= phase of increase in the curve, c= phase of decline in the curve, k=parameter associated with the days to the peak, and t= days in production.
The calculated parameters are used to estimate the persistence (), days to the peak (), and peak yield ()(14).
The Ali-Schaeffer model(15) is:
Where: t= days in milk, a= related to peak production, b and c= related to decreased production, d and f= related to increased production.
The modified Pollot model(8) is described as:
Where: yt= milk production at day t, t= days in milk, a= maximum lactation secretion potential, b= related to milk production potential, c= relative proliferation rate of secretory cell number during early lactation, and d= relative decrease in cell number as lactation progresses.
Results
According to the selection criteria, the best-evaluated models were Wood, Wilmink, Ali-Schaeffer, and modified Pollot. Table 1 shows the results of the four models and the values of the selection criteria for estimating the milk production curves for Holstein cattle in Mexico in the intensive production system. In most of the criteria, Wilmink’s model is the one with the best results.
Table 1: Selection parameters of Wood, Wilmink, Ali-Schaeffer, and modified Pollot models in Holstein cattle in Mexico
R2= coefficient of determination, R2adj= adjusted coefficient of determination, RSE= residual standard error, LogL= log-likelihood, AIC= Akaike information criterion, CAIC = corrected Akaike information criterion, BIC= Bayesian information criterion, DW= Durbin-Watson coefficient.
Table 2 shows the results of the ANOVA and Tukey tests for the parameters of the four selected models, differentiated by the number of lactations, and the mean of all the animals. In Wood’s model, it is observed that estimators a, b, and c are statistically different between the different lactation numbers, except for the estimator c for lactations 3 and 4 or more. The values of persistence, peak production, and days to the peak for each lactation are also presented.
Regarding Wilmink’s model, parameter a of lactation 1 differed from those of lactations 2 and 3, which in turn differed from that obtained for 4 or more lactations. Regarding parameters b and k, there were no significant differences between the groups; in contrast, in parameter c, lactations 2 and 4 are the same but differ from the rest. Ali & Schaeffer’s model shows that parameters a, b, c, d, and f in lactations 1 and 4 are statistically different from those in lactations 2 and 3.
For the modified Pollot model, it is observed that parameter a is different between the first, fourth, and second-third lactations; as for parameter c, that of the first lactation differs from that of all the others; as for parameters b and d, they are different between lactations.
Table 2 also shows the estimated values of persistence, peak production in kilograms, and days to the peak by lactation for Wood, Wilmink, Ali-Schaeffer, and modified Pollot by calving number, as well as the mean for all animals. Figure 1 shows the lactation curves for each of the models.
In the evaluation of factors that affect the lactation curve presented in Table 2, it was found that the number of calving, calving season, and herd are significant (P<0.05) in the models used, except for the herd in the Ali-Schaeffer model. Table 3 shows the parameters for fat and protein with the different models used, where it can be seen that all the parameters are different between the number of calvings (P<0.05). Table 4 shows the estimates for the components of the curve with the different models, while Figures 2 and 3 show the curves calculated for fat and protein, respectively. In the modified Pollot model for protein, the parameters did not model a curve in lactations 3 and 4 or more, so it was not possible to obtain the days to the peak, peak production, and persistence.
Figure 1: Lactation curves for milk production by lactation number with Wood, Wilmink, Ali & Schaeffer and modified Pollot models
Figure 2: Lactation curves for milk fat percentage by lactation number using Wood, Wilmink, Ali & Schaeffer and modified Pollot models
Figure 3: Lactation curves for protein percentage by lactation number with Wood, Wilmink, Ali & Schaeffer and modified Pollot models
Discussion
The mean square error was similar among the different models, being slightly lower in the first lactation. The same is true for the other model selection criteria, with Wilmink’s model being slightly better.
The parameters obtained by the model proposed by Wood are different in lactation numbers, results that are far from those found by Duque et al(16) with Wood’s model in the Colombian tropics with grazing Holstein cattle. Duque et al (2018) estimated a mean of parameter b (0.12) and peak production (26.5 kg) lower than what was found in the present study; the same happens with the days to reach maximum production among the different numbers of lactations (between 28 to 32 d). It is known that milk production in the tropics is usually lower due to various factors that limit production, such as temperature, where Holstein cows do not adapt adequately to hot climates; in addition to this, being in an extensive system, grazing feed tends to vary at different times of the year. In terms of persistence, they had higher values (66 to 82 %) than reported in this study, which ranges from 15 to 19 %. This may be because they are subjected to less production stress in addition to the variation in diet depending on the time of year and heat stress.
In a study conducted by Vázquez et al(17), where they evaluated cows mostly of the Holstein breed under an intensive system in Lima, Peru, the values of parameters a and c of Wood’s model among the different lactations (16.41 to 18.11, and 0.0023 to 0.004, respectively) and the peak production (31.13 to 43.91 kg) are similar to those found in the Holstein population of Mexico. In both studies, the animals were subjected to intensive production systems, and the climatic conditions were similar, corresponding to a subtropical desert climate. For parameter b and days to peak production, Vázquez et al(17) show lower values (0.1880 to 0.3043 and 66 to 82 d, respectively). This difference could be attributed to the fact that the cows in the Peruvian study were milked 3 times a day, so the amount of milk produced in the first stage of the curve increases compared to cows that are milked 2 times a day, as is the case of the majority of the Mexican population.
In a study by Boujenane & Btissam(18) in semi-intensive production herds in Morocco with Holstein animals, the results show some differences compared to this study. The values by lactation reported in Morocco for parameter a in the first three calvings (15.9, 16.9, and 17.2, respectively) present a higher value, especially for first-calving animals. The mean of parameter b (0.1039) is the one with the biggest difference, which is reflected in the results of each lactation, where they are also higher (0.073, 0.091, and 0.096 for the first, second, and third lactations). In the same study, the parameter c is slightly higher in each lactation.
In terms of lactation components, there are differences in both studies. The three components shown by Boujenane & Btissam(18) (41.4 for days to the peak, 23.6 for peak production, and 6.56 for persistence) are lower than those found in the present study (Table 2), especially in days to the peak and peak yield. This is possibly caused by the semi-intensive production system and the high temperatures of the African country. In general, production in the animals in the Moroccan study is lower.
Regarding fat percentage, Gołębiewski et al(19) conducted a study on Holstein cattle in Poland, reporting values of 3.05, -0.07, and 0.04 for parameters a, b, and c, respectively, with Wood’s model; and for protein percentage, they report values of 4.59, -0.19, and 0.04, for parameters a, b, and c, values similar to those found in the Holstein population of Mexico (Table 3). Both the Polish and Mexican Holstein animals in both studies were under an intensive system, so the environmental conditions are similar.
Regarding Wilmink’s model, the results presented by Bouallega et al(20) in Holstein cows in Tunisia, the value of the parameters differs from that obtained in the Mexican Holstein population. The authors show values close to 28 and -7 for parameters a and b. The parameter c was similar to what was obtained (-0.3), while the value of k was set to 0.05. The values calculated in the present study for peak production and days to the peak (Table 4) were higher than those presented by the authors (26 kg and 48 d). In terms of persistence, they report values around 94 %. The number of animals used in the study was small (5,649), where the authors mention that a larger amount of data is suggested; the main difference between the cows in the Tunisian research and the Mexican population was temperature. The former were subjected to heat stress due to the climate in Tunisia, which can reduce production in Holstein cattle since animals of this breed do not usually adapt well to this type of climatic conditions.
Regarding parturition, the results found in this study are similar to those reported by Bouallega et al(20), which reiterates that the number of lactations is a significant source of variation, showing differences in animals with 1, 2, and 3 or more births, because first-calving animals have not completed the mammary gland maturation process; therefore, their production is usually lower than in subsequent lactations. In addition, Bouallega et al(20) recommend using the age at calving as a source of variation. Regarding herd as a factor affecting lactation curves, the aforementioned authors found that it is significant, attributing 30 % of the variation in milk production to it. This highlights the importance of environmental conditions in modeling the production curve, which are different between herds.
The protein percentage parameters show little similarity compared to what was reported in the present study (Table 4). An example is the peak and days-to-peak results shown for Tunisian animals (2.84 % and 53.4 d, respectively), which are lower than those of the Mexican population. In terms of fat percentage, the biggest differences in comparison to this study are observed in parameter b (1.19), where the authors present slightly higher values, while the days to the peak (50.63) they show are much lower than those found in this study (Table 4). The differences in some results for protein and fat may be due to the high Mediterranean temperatures to which the animals were exposed; however, the percentage of fat tends to vary less due to environmental conditions and during lactation than the percentage of protein.
Torshizi et al(7) found that, in first-calving Holstein cows under intensive production systems in Iran, herd and calving season are sources of significant variation using Wilmink’s model to model lactation curves, similar to what was found in the Mexican population. In addition, they used 4 fixed values for parameter k in their analyses (0.05, 0.065, 0.61, 0.10), the first being the one that yielded the highest correlation between the observed and predicted production values.
Regarding the other parameters, they are also very different, being most evident in parameters b and c (-20.227 and -0.036). Peak day and peak production (66 d and 32 kg, respectively) are lower in Iranian cows. These differences may be due to the length of lactations since Torshizi et al(7) study included animals with production cycles adjusted to 200 d, contrary to the usual in Mexico, which is to adjust to 305 d, and minimum productions of 3 kg of milk. This is due to the variation in production due to climatic issues. The authors mention that the best model for first-calving animals is Wood’s.
Results presented by Gök et al(21), with first-calving Holstein cows in the Turkish province of Konya, where they used the Ali-Schaeffer model, show values similar to those found in the present study for parameter b (138). Regarding parameters a and f, these authors show higher values (-51.92 and -3.62); in contrast, the estimates for parameters c and d (-648.66 and 32.68) are lower compared to those found in the Mexican Holstein population (Table 4). In the same way, days to the peak and peak production are lower in animals from Konya. The main difference with this study was the production system. Turkish animals were under a grazing system, and where the climate is usually extreme in the different seasons of the year.
In their study of Holstein cattle from Turkey, where they classified the animals by lactation number (from 1 to 3), Koçak & Ekiz(15), by using the Ali-Schaeffer model, reported values for parameter b (165.3, 259.3, and 280.9) and c (-101.3, -121.1, and -127.0, respectively) similar to those found in this study (Table 2). In contrast, the rest of the parameters were higher in the research carried out in Turkey (-49.0, -55.7, and -50.7 for a; 103.36, 126.97, and 41.58 for c; -0.10, -6.91, and -14.71 for f). These differences in parameters are reflected in the days to the peak, where Turkish cows have their maximum production (74.94, 47.62, and 39.62 for the first, second, and third lactations, respectively) earlier than Mexican cows, despite the fact that productions at this stage are similar. The animals in Koçak & Ekiz’s(15) study belonged to semi-intensive production systems and were milked three times a day; in said study, despite the high environmental temperatures, the houses had temperature regulation systems, unlike what happened in Mexico, where the environmental conditions of the animals are not usually controlled.
On the other hand, Nanda et al(22) carried out a study on a housed herd in Indonesia, where they modeled the curves using Ali & Schaeffer’s model by calving numbers (from 1 to 4). The parameters a, c, and f in each of the lactations of the cows in Java were higher than the parameters calculated in the present study (-40.79, -16.19, -20.86, and -26.89 for a; -16.50, -7.06, -14.74, and -25.00 for c; -6.83, -4.59, -4.52, and -4.59 for f); with respect to parameters b and d, they report lower values (68.32, 32.87, 44.25, and 59.15 for b; 38.85, 24.43, 25.52, and 26.83 for d). These values may indicate that the curves of the animals in the Indonesian study show peak production in a shorter time and yield lower than in the Mexican population. The hot and humid climate of the island of Java may be an important factor explaining the differences in the curves, as these animals were not housed in places where temperature was controlled.
In Holstein cows in Australia under a grazing system in a warm climate, Adediran et al(23) used the modified Pollot model and found that parameter a (13.36) was lower than that obtained in the present work, while parameter b was higher (1.23). Parameters c and d (2.80 and 0.0012) were similar in both studies. The main difference with the Mexican population is the type of production system. In the study in Australia, as the animals are grazing, there is less control of the environmental conditions, coupled with the intense heat reported by the authors, which may limit milk production, unlike the population in this study, which was in a housed system.
Information on the parameters of the fat and protein production curves with the different models is scarce, especially for Ali-Schaeffer and Pollot modified, so the results were compared with other models; however, it should be noted that the results presented in this study will serve as a reference for these characteristics in intensive production systems.
Conclusions and implications
Of the models evaluated in this study, the one proposed by Wilmink was the one that best fit the data of the Holstein population of Mexico. The importance of choosing a model that best suits the information lies in obtaining more accurate predictions, which translates into values that are closer to reality. In addition, the study evaluated environmental factors such as calving number, calving season, and herd; they were significant for the modeling of lactation curves, so it is essential to consider them as a source of variation in the predictions made with the different models. Future research could investigate other environmental factors that may affect the curves. The practical application of lactation curve modeling is extensive, including genetic improvement, so having identified environmental sources of variation and choosing the most appropriate model will allow the selection of animals with the highest genetic value.
Table 2: Lactation curves parameters, days to the peak (Dpeak), peak production (peak), and persistence of the curve (persistence) for Wood, Wilmink, Ali-Schaeffer, and modified Pollot by lactation number.
Dpeak= days to the peak, Peak= peak production, MSE= mean square error.
abcd Significant differences at P<0.05.
Table 3: Parameters of lactation models for protein and fat percentages with Wood, Wilmink, Ali-Schaeffer, and modified Pollot models by lactation for the models with better goodness of fit
MSE= mean square error. abcd Significant differences at P<0.05.
Table 4: Days to the peak, peak yield, and persistence for milk protein and fat percentages with Wood, Wilmink, Ali-Schaeffer, and modified Pollot models.
Literature cited: