An accurate characterization of the lactation curve in dairy cattle has become essential for herd management and animal breeding (Gipson and Grossman, 1990). Moreover, the identification of biological and environmental factors influencing its trend and slope must be of special interest for dairy farmers (Masselin et al., 1987), as these sources of variation remarkably modulate the economic profitability inherent to each cow. The statistical characterization of the lactation curves allows us to predict total milk yield from one or a few test-day records (Goodall and Sprevak, 1984), and even milk composition (Stanton et al., 1992).
The first endeavors to mathematically characterize lactation curves date from the beginning of the 20th century (Grossman and Koops, 1988), although a plethora of parametric and non-parametric designs were developed from then on (Delange et al., 1953; Pérochon et al., 1996; Grossman and Koops, 2003). Nevertheless, substantial limitations were reported in tropical dairy cattle (Val-Arreola et al., 2004) where some parameterizations such as Wood's curve (Wood, 1967) poorly fit monthly data (Papajcsik and Bodero, 1988; Val-Arreola et al., 2004). These departures from non-tropical and intensive dairy production could be linked to limited energy balance and poor forage quality.
In Colombia, where milk is produced in tropical conditions using different dairy cattle breeds and their crosses, the first genetic evaluations on a larger scale were made from 1997 with the participation of the Colombian Association of Zebu Cattle Producers (Asocebú), Colombian Agricultural Research Corporation (Corpoica) and the Universidad Nacional de Colombia in Bogota. However, there are still many limiting factors in this country to be overcome for the development and establishment of breeding programs (Quijano, 2002) and frequently the selection is not made based on the genetic evaluations of the herds, but on their phenotype. Several lactation curve works have been done, but only some of them applying adjustment models as for example those of Cañas et al. (2011, 2012) and Vargas et al. (2016).
Within this context, this research evaluated six different lactation curve models in a Colombian Holstein dairy herd data, not previously evaluated, in order to assess their performance and goodness of fit in its tropical management system and recommend the most appropriate to be used for animal breeding purposes.
MATERIALS AND METHODS
Field data source
This research focused on the experimental dairy herd of the Centro Paysandú (Santa Elena, Medellín, Colombia, 6°12'37"N and 75°30'11"W) of the Universidad Nacional de Colombia, tropical region located 2500 m of altitude, with average temperature, rainfall and humidity of 17°C, 1861 mm and 79%, respectively. This Holstein cows herd was founded in 1956 donated by the W. K. Kellogg foundation (http://www.wkkf.org), maintaining a herd size between 80 and 90. Cow feeding is based on rotational grazing of Pennisetum clandestinum meadows supplemented with 3 to 10 kg of standard concentrate (16.8% protein; 33.6% NDF) depending on production level. Colombian (40%) and international (60%; i.e. USA, Canada and Europe) Holsteins bulls were used for AI, with manual (up to 1988) or machine milking (after 1988) done twice a day at 4 am and 2 pm.
Lactation data
The data set included 425 full lactations and 244,876 test day records taken during fourteen years (1986-1999) from the Colombian Holstein herd. Data were available for the first 9 lactations although the sixth and following lactations were grouped due to their reduced incidence (first lactation: 91; second: 80; third: 76; fourth: 66; fifth: 51; sixth and following: 61). Lactation records were restricted to 305 days and records of average daily milk production of each week (AWP) were calculated by accumulating appropriate test day records; Missing records (13%) for shorter-than-305 days lactations were augmented by applying lactation number-specific persistence coefficients estimated as the average decrease between two successive production weeks for each group of lactation number. 305 days adjusted milk yield (MY305d) was calculated for each lactation by adding AWP1x7 + AWP2x7 + … + AWP42x7 + AWP43x10, where AWPk indicated the AWP record for the kth week of lactation. Lactation length (LL) was defined as the number of days elapsed between calving date and drying date. Production peak (PP) was determined as the AWP with maximum milk production, whereas the time to the PP (TPP) was the number of weeks between calving date and PP date. Lactation persistence (LP) was evaluated for a 6 month period, calculated as the percentage of milk yield still producing 27 weeks after PP, when comparing with milk yield two weeks after PP.
Statistical analysis of dairy traits
Records of actual dairy traits (i.e., MY305d, LL, PP, TPP and LP) were analyzed under the following univariate mixed linear model with ANOVA properties:
where Yijkl was the phenotypic record, eijkl was the residual, CNi was the effect of the calving number with 6 levels (1st, 2nd, 3rd, 4th, 5th, and 6th and subsequent calving), CSj was the effect of the calving season with two levels (Dry: December to March, July and August; Rainy: April to June, September to November), YCk was the effect of the year of calving with 14 levels (from 1986 to 1999), and vl was the random contribution of the lth cow. All analyses were performed with the Mixed Model procedure of SAS.
Lactation curve models
Six non-linear lactation curve models being commonly implemented for dairy cattle lactation curves worldwide were preliminary selected and evaluated on the basis of their goodness of fit on test day records. These parametric curves were the exponential model developed by Brody et al. (1924), Nelder's (1966) inverse polynomial model, Wood's (1967) incomplete gamma curve, Cobby and Le Du's (1978) model, and both single-phase, and two-phase curves, previously reported by Grossman and Koops (1988). Their Mathematical characterization is shown in Table 1, where yt was a test day record at time t (t = 0 being the calving date), and a, b, c, d, g and h where unknown parameters of the lactation curve models.
Comparing lactation curves and systematic sources of variation
Each lactation curve model was fitted on a lactation basis as suggested by Tekerli et al. (2000), Macciotta et al. (2005) and Dematawewa et al. (2007). Analyses were performed on AWP records and with the non-linear mixed models (NLMIXED) procedure of the SAS v.9.2 software (SAS Inst. Inc., Cary, NC). In order to compare the goodness-of-fit for each parameterization of the lactation curve, the means square error (MSE) was computed for each model. Additionally, model fit was compared by the Mixed Models procedure of SAS. The statistical model included the fixed effects of calving number (CNi) with 6 levels (1st, 2nd, 3rd, 4th, 5th, and 6th and subsequent calvings), calving season (CSj) with two levels (dry: December to March, July and August; rainy: April to June, September to November), their first-degree interaction, the lactation curve parameterization (LCk), and the random effect of each cow (vl) as follows:
Lactation curve parameters were also analyzed by the Generalized Linear Models procedure of SAS (SAS Institute, Inc.) and under the following model,
where θijk was the lactation-specific estimate for a given lactation parameter.
RESULTS AND DISCUSSION
Systematic sources of variation on dairy traits
Data from Colombian Holsteins involved in this study averaged 5830 ± 59 kg, 330 ± 3 d, 27.7 ± 0.3 kg, 4.7 ± 0.1 wk and 63.1 ± 0.6% for MY305d, LL, PP, TPP, and LP, respectively, being the age at first calving 30.0 ± 3 mo, and the average daily milk yield 19.3 ± 0.1 kg.
Results from the analyses of systematic sources of variation on MY305d, LL, PP, TPP and LP are shown in Table 2. Number of calving, calving season and year of calving significantly influenced all dairy traits (P<0.05) excepting number of calving on TPP (P=0.44) and calving season on LL (P=0.79) and LP (P=0.35). More specifically, MY305d, LL and PP increased (P<0.05) with the number of calving (from 1 to 5), although suggesting a slight decrease for older cows with more than 5 lactations (Table 3); In contrast, LP was higher (P<0.05) for first lactation cows than for the remaining categories.
The original data set of Centro Paysandú farm was obtained between 1986-1999, but if a more recent work was done, probably the results would be similar to those presented here, since in the last 30 years there have not been any notable changes in this farm from the environmental and/or genetic point of view. Peak of production observed during the first 2 lactations was higher (+ 21% and + 9%, respectively) than indicated by Cañas et al. (2009) in Colombian Holstein cows. Both MY305d and PP, increased during the first five lactations, showing a stationary pattern afterwards, agreeing with results previously reported by Cerón-Muñoz et al. (2003) and Restrepo et al. (2008) in tropical dairy production systems. The MY305d was slightly lower (5830 vs 6212 kg) than in the case of Cañas et al. (2012) in Colombian herds, while the average PP (27.7 ± 0.3 kg) was similar to the value indicated by the same authors. On the other hand, our persistency (63.1%) can be considered normal or even higher than usual (i.e. 90% monthly) because was calculated for a period of six months.
In accordance with Macciotta et al. (2011), first calving cows showed a lower peak of production and a higher persistency, compared to later parities (Figure 1). MY305d in the fifth lactation was 27% higher than in the first one, matching the 25%-increase suggested by Hurley (2003) and attributed to the increase in body weight (~20%) and development of the mammary gland (~80%). The trends observed in first calving cows, with lower total and peak milk production, an increase in peak time and greater persistence, could be explained, according to Dijkstra et al. (1997), for the lower rate of proliferation of mammary secretory cells in primiparous, which can affect the speed to reach the peak, and, according to Stanton et al. (1992), counteract the normal decline of the lactation curve, generating greater persistence. In addition, the greatest productive increase between the first and the second may be due to a very important increase in the cells of the mammary parenchyma between both lactations (Pollot, 2000); this shape could be compatible with the phenotypic patterns shown in Colombian Holstein. On the other hand, TPP was not influenced by the lactation number, agreeing with Cañas et al. (2009) and contradicting previous results from Macciotta et al. (2004, 2011).
Differences between calving seasons evidenced a small but significant increase (P<0.05) in MY305d and PP traits for cows calving during the rainy season, and this peak of production was obtained earlier during the lactation period (Table 3). The departure between dry and rainy seasons was small for both MY305d and PP trait, increasing less than 3% during the rainy season. On the contrary, TPP was reduced by 13.6% during the rainy season. These results are comparable to the ones provided by Cerón-Muñoz et al. (2003) in Colombian herds, where cows that calved during rainy months reached higher production level.
Adjustment of lactation curve models
Average lactation curves were characterized by the estimates provided in Table 4, and corresponding average time to the production peak (wk), milk yield in he production peak (kg d-1) and persistence (%) were summarized in Table 5. Comparing with field data estimates in Table 3, time to the production peak (with values varying from 4.5 to 5.1, depending of number of calving) was clearly underestimated by all lactation curve models, with the only exception of the two-phase lactation curve (Grossman and Koops, 1988) that averaged 6.0 ± 0.2 wk (Table 5).
On the contrary, departures for milk yield in the production peak were small, these estimates ranging between 25.2 kg d-1 (single phase model; Grossman and Koops (1988)) and 27.4 - 27.6 kg d-1 (inverse polynomial model, Nelder (1966), and exponential model, Brody et al. (1924), respectively). Persistence showed a moderate dispersion pattern with minimum and maximum estimates provided by the two-phase model (Grossman and Koops (1988) and Wood's (1967), averaging 60.0 ± 0.5% and 74.0 ± 1.0%, respectively.
Although all models evidenced different patterns (Figure 2), some of them, particularly Single-phase (Grossman and Koops, 1988; Wood,1967), provided a bad prediction of the PP, a phenomenon previously reported by Grossman and Koops (1988) and Olori et al. (1999). These authors suggested that three-parameter exponential models suffered from important limitations during the first 90 days of lactation. On the contrary, Scherchand et al. (1995) and Vargas et al. (2000) provided evidences about a reasonable fit for the two-phase lactation curve between days 60 and 90 of lactation, whereas this model overestimated TPP and underestimated the PP in our Colombian Holstein data set.
Goodness of fit for each lactation curve model was evaluated by the MSE (Table 6) of estimated average daily milk yield. This statistic characterized the magnitude of the unaccounted variability for the average daily milk production of each week or AWP. On the basis of its MSE, goodness of fit for lactation curve models was maximum (P<0.05) with the two-phase model (Grossman and Koops, 1988), and progressively decreased with the incomplete gamma (Wood, 1967; Cobby and Le Du's, 1978) and the remaining models. Lactation curve parameters also departed across systematic effects such as calving number as shown for the best fitting models in Table 7.
Our results agreed with Vargas et al. (2000) who favored the two-phase curve in Costa Rican Holstein dairy herds, and discarded others. In our case, the worst goodness of fit was achieved by the exponential model of Brody et al. (1924) as reported by Quinn et al. (2005), whereas Vargas et al. (2000) reported worse fits with Wood's curve. Limitations of the Wood model exist, but this model possesses the ability to fit four different shapes including the atypical curves (Macciotta et al., 2011). In addition, Dematawewa et al. (2007) and Macciotta et al. (2011) indicated the good fitting performances of the multiphasic approach but also the computational problems related to the estimation of a high number of parameters. These authors concluded in favor of simpler models, as the Wood function, particularly for routine use. Similarly, several Colombian authors Cañas et al. (2011, 2012); Vargas et al. (2016), under different circumstances, obtained the best results using the Wood model. In our case the biphasic function gave the best results, followed by the Wood model. Therefore, for research and technical purposes the two-phase model of Grossman and Koops (1988) would be recommended, but for common livestock use the simpler model of Wood would probably be the best choice.
CONCLUSIONS
Under the circumstances of this tropical Colombian Herd, and according to the MSE criteria, these results provide substantial evidence favoring the two-phase lactation curve of Grossman and Koops followed by the Wood model. The biphasic model would be recommended for research and technical purposes, while the Wood model would be the best for livestock use.