Short communication: Genetic analysis of lactation curves in buffaloes, using Wood’s model

Aim of study: To estimate the heritability and genetic correlations for lactation curve traits in buffaloes.Area of study: The buffalo cows were raised on properties located in the states of São Paulo, Ceará and Rio Grande do Norte, Brazil.Material and methods: The individual parameters of Wood’s model ( , , and ) were obtained using a non-linear mixed model. Peak yield (PY), peak time (PT) and lactation persistency (LP) were also calculated. These individual parameters were employed in multi-trait analysis with the milk yield (MY) using Bayesian inference.Main results: The heritability estimates were of low to moderate magnitudes, with values ranging from 0.156 ( ) to 0.299 (PY). The estimates for genetic correlation between the Wood’s parameters and MY were of low to high magnitude and ranged from -0.533 (  and MY) to 0.983 (PY and MY).Research highlights: The heritability estimates obtained indicate that the traits studied can be used in animal breeding programs.


Introduction
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 methods
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; θ 1 has a scale meaning and tends to increase with parity, i.e. according to the production level of animals; θ 2 controls the rate of increase to the lactation peak; finally, θ 3 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 = ( , ) + ( , , ) (.) are functions that define the structure of the individual effects and the residuals, respectively (Comets et al., 2017). For = ( , ) + ( , , ) (.) 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 (θ 1 , θ 2 , and θ 3 ), 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. Contemporary 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.
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 θ 2 and θ 3 ) 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 3 Short communication: Genetic analysis of lactation curve in buffaloes  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, θ 1 , θ 2 , θ 3 , PY, PT and LP, respectively (Fig. 3). The estimates obtained in this study indicate that the highest genetic gains can be achieved for MY, θ 1 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 θ 1 and θ 3 , 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 θ 2 , 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 θ 1 , θ 2 , and θ 3 , 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.
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 θ 1 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 θ 2 and θ 3 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.