IntroductionTop

Buffaloes (Bubalus bubalis) are used in many countries for mozzarella cheese production because of the higher percentage of milk protein and fat. Milk yield can be defined as a longitudinal trait whose modeling in lactation curves is of great importance for the development of breeding programs for dairy species (Akers, 2000). Several models have been proposed for the longitudinal analysis of milk production, such as factor analysis (Aspilcueta-Borquis et al., 2012; Pereira et al., 2013), random regression (Pereira et al., 2010; Aspilcueta-Borquis et al., 2013) and nonlinear model (Hernández et al., 2014; Bangar & Verma, 2017) in dairy species.

The Wood’s incomplete gamma function is probably the most popular empirical model to generate lactation curves. It generates the standard shape of the lactation curve as the product of a constant, a power function and an exponential decay function (Macciota et al., 2011). The Wood’s function provides a straight forward way of calculating parameters with a direct biological interpretation and economical interest (persistency, peak yield, time to peak yield). Studies on dairy cattle have indicated the possibility of selecting animals to change the lactation curve (Rekaya et al., 2000; Chang et al., 2001). Therefore, in view of the importance of selecting for lactation curves in buffaloes, this study was carried out aimed to estimate parameters of the Wood’s curves for milk yield in buffaloes. This information was used to estimate genetic parameters (heritability and genetic correlations) for peak yield (PY), peak time (PT) and lactation persistency (LP).

Material and methodsTop

We used the information of milking control of primiparous buffalo cows, raised on properties located in the states of São Paulo, Ceará and Rio Grande do Norte, Brazil. The research was conducted in three steps: a) estimation of the individual parameters of the Wood’s curve; b) estimation of parameters defining the shape of the lactation curve; and c) estimation of the genetic p arameters.

For the first stage, information from 1470 animals with 9 monthly controls was used. Wood’s model was applied to fit the lactation curves:

where y_t represents the test-day milk yield in the month t; θ₁ has a scale meaning and tends to increase with parity, i.e. according to the production level of animals; θ₂ controls the rate of increase to the lactation peak; finally, θ₃ expresses the rate of decline after the peak and its absolute value tends to increase with parity (Macciota et al., 2011). The parameters of the nonlinear mixed models were estimated using the Saemix package (Comets et al., 2017) in which the Expectation Maximization (EM) algorithm with a stochastic approximation is implemented. A heterogeneous residual structure was used and the average curve was the only fixed effect considered. The nonlinear mixed model used in this case can be described as:

where y_ij, x_ij and ϵ_ij represent, in this order, milk yield, lactation month and error related to the i_th animal and j_th test-day; θ_i is the vector of individual parameters of Wood's function; f(.) and g (.) are functions that define the structure of the individual effects and the residuals, respectively (Comets et al., 2017). For g (.) was considered a function composed of a constant and a parameter dependent on θ_i and x_ij. In this step, discrepant records were eliminated by residual analysis. The final database structure is shown in Fig. 1.

Using the individual parameters estimated in the first step (θ₁, θ₂, and θ₃), peak yield (PY), peak time (PT) and lactation persistency (LP) were calculated for each animal:

For estimation of the genetic parameters, total milk production up to 270 days (MY) was included in the analyses. The MY was obtained using interpolation between controls (method used in the genetic evaluation of these herds). Contemporary groups were formed using the variables farm, year and season of calving. Contem-porary groups with at least three animals were kept and outlier records were eliminated (considering 3 standard deviations of the mean of the contemporary groups). The variance components were estimated by a Bayesian approach with Gibbs sampling under a multi-trait animal model, using the GIBBS2F90 software (Misztal et al., 2019). The model can be represented as follows:

where y, β, a and ε are the vectors of observations, systematic effects (contemporary groups and age at calving), additive genetic random effects and residual random effects, respectively; X and Z are incidences matrices that associate vectors β and a with vector of observations y.

As prior information for Bayesian inference, uniform, normal multivariate and inverted Wishart distributions were considered for the systematic effects (β), genetic effects (a) and (co)variance matrices, respectively. A total of 650,000 samples were generated. The first 150,000 iterations were discarded as burn-in and one sample was collected every 50 iterations (thinning interval). The POSTGIBBS1F90 software (Misztal et al., 2019) was used for post-Gibbs analysis.

Results and discussionTop

The mean and standard deviation of the individual estimates for parameters θ₁, θ₂, and θ₃ were 8.52 (2.18), 0.59 (0.21) and 0.19 (0.06), respectively (Fig. 2). Working with information from buffaloes in Iran, Hossein-Zadeh (2015) described estimates of 5.42, 0.0968 and 0.001520 for θ₁, θ₂, and θ₃, respectively, using first lactation information. Şeahin et al. (2015) reported a lower estimate for parameter θ₁ (5.36), but similar estimates for the other parameters (0.58 and 0.23), using information of Mediterranean buffaloes.

For the population studied, the values of 9.248 kg for PY, 3.125 months for PT and 2.657 for LP were obtained (Fig. 2). In this population, PT occurred in the third month (90 days), a value higher than that described by Hossein-Zadeh (2017) at 64 days but similar to that obtained by Şeahin et al. (2015), which was close to 80 days of lactation. The LP values obtained were similar to those reported by Şeahin et al. (2015) for buffaloes (2.45) and by Bangar & Verma (2017) for Gir cattle (2.60).

Regarding the type of curves obtained, animals with atypical curves (negative estimates for the parameters θ₂ and θ₃) were not observed in the herds studied. Atypical curves are presumed to be associated with buffaloes that exhibit various data-related problems, such as lack of information, erroneous data, or management problems (Şeahin et al., 2015), indicating the quality of the database used.

The heritability estimates were of low to moderate magnitudes, with values of 0.261, 0.271, 0.155, 0.167, 0.299, 0.191 and 0.187 for MY, θ₁, θ₂, θ₃, PY, PT and LP, respectively (Fig. 3). The estimates obtained in this study indicate that the highest genetic gains can be achieved for MY, θ₁ and PY, in view of the greater influence of the additive genetic effect.

The heritability estimates for MY reported in the literature for buffaloes range from 0.14 and 0.46 (Rosati & Van Vleck, 2002; Pareek & Narang, 2014; Hossein-Zadeh, 2015). The heritability obtained in this study for parameters θ₁ and θ₃, of the Wood’s curve were similar to those described by Rekaya et al. (2000), in Holstein cattle: 0.23 and 0.17, respectively. However, the value obtained for parameter θ₂, was lower than that described by these authors: 0.36. In dairy sheep, Chang et al. (2001) reported heritability estimates higher than those found in this study: 0.35, 0.35 and 0.27 for parameters θ₁, θ₂, and θ₃, respectively. The estimates obtained for PY, PT and LP using Wood's function, the estimates were higher than those described by Boujenane & Hilal (2012) in Holstein cows, with values of 0.06, 0.10 and 0.05 for PT, PY, and LP, respectively.

In view of its economic importance, several methodologies have been applied to estimate LP. Aspilcueta-Borquis et al. (2012) described lower heritability estimates for PY (0.08) and LP (0.12), using multivariate methods (latent variables) and buffalo information. Using random regression models and different definitions of LP in buffaloes, Hossein-Zadeh et al.. (2017) estimated heritabilities ranging from 0.09 to 0.26.

The genetic correlation estimates of Wood’s parameters with MY, were of low to high magnitude: 0.880, -0.533, -0.506, 0.983, -0.146 and 0.242 for θ₁, θ₂, θ₃, PY, PT and LP, respectively (Fig. 3). Rekaya et al. (2000) described similar correlation estimates between MY, PY (0.90) and LP (0.33) in Holstein cattle. Wasike et al. (2014) reported similar correlations with PY (0.976) and a lower estimate for LP (0.053). However, discrepancies were found in the literature for the correlation between MY and PT, with estimates ranging from -0.129 (Wasike et al., 2014) to 0.59 (Rekaya et al., 2000).

Considering milk production as the main selection criterion, the estimates heritability and genetic correlation estimates obtained, we verified that selection for MY will promote significant changes in θ₁ and PY. But we note that a small change will be promoted in other traits of the shape of the curve, since estimates of genetic correlation with the parameters θ₂ and θ₃ present moderate magnitudes.

PT and LP are two aspects of the lactation curve that are important for dairy farms, and the genetic correlation estimates with MY indicate a possible independence (considering the high density intervals, Fig. 3) and lower indirect genetic gains. From a bio-economic point of view, the goal should be to select animals with lactation curves with moderate PY and higher LP and, regarding PT, the selection of animals with late peaks can be applied to reduce stress and to allow body reserves to be used more slowly (Ferris et al., 1985). Thus, an alternative would be to introduce these traits as auxiliary selection criteria combined in a selection index (Ferris et al., 1985; Grosman & Koops, 2003).

The heritability estimates obtained for the population in this study indicate the possibility of using the parameters of Wood’s function in breeding programs. However, the magnitudes and direction of the genetic correlations with MY should be considered in the elaboration of selection indices.