Random regression models to estimate genetic parameters for test-day milk yield and composition in Iranian buffaloes

The objective of this study was to estimate genetic parameters for milk yield and milk percentages of fat and protein in Iranian buffaloes. A total of 9 278 test-day production records obtained from 1 501 first lactation buffaloes on 414 herds in Iran between 1993 and 2009 were used for the analysis. Genetic parameters for productive traits were estimated using random regression test-day models. Regression curves were modelled using legendre polynomials. Heritability estimates were low to moderate for milk production traits and ranged from 0.09 to 0.33 for milk yield, from 0.01 to 0.27 for milk protein percentage and from 0.03 to 0.24 for milk fat percentage, respectively. Genetic correlations ranged from −0.24 to 1 for milk yield between different days in milk over the lactation. Genetic correlations of milk yield at different days in milk were often higher than permanent environmental correlations. Genetic correlations for the milk protein percentage ranged from −0.89 to 1 between different days in milk. Also, genetic correlations for the milk percentage of fat ranged from −0.60 to 1 between different days in milk. The highest estimates of genetic and permanent environmental correlations for milk traits were observed in adjacent test-days. Ignoring heritability estimates for milk yield and milk protein percentage in the first and final days of lactation, these estimates were higher in the 120 days of lactation. Test-day milk yield heritability estimates were moderate in the course of the lactation, suggesting that this trait could be applied as selection criteria in Iranian milking buffaloes.


Introduction
The world buffalo milk production increased by 43 % from 1997 to 2007.This rate is superior to that reported for cow milk production (20 %) in the same period (FAO 2009).However, the buffalo milk yield is still much lower than that of cows.In 2007, the worldwide buffalo and cow milk yields were 85 and 560 million tonnes, respectively (FAO 2009).Higher average milk yields are reported in India and Italy, probably because genetic evaluations are a common practice in these countries (Moioli & Borghese 2005).
According to climate conditions, Iranian buffaloes can be classified in three main groups: i) Azari ecotype (Western and Eastern Azerbaijan); ii) North ecotype (Guilan & Mazendaran); iii) Khoozestan ecotype (Khoozestan).Iranian water buffaloes have also some similarity to Iraqi buffaloes (Tavakolian 2000).Therefore, it may be that both groups have their origin in the same ancestor.Furthermore Iranian buffaloes in the northwest of the country (West Azerbaijan) have a great similarity to Mediterranean water buffaloes.Thus, it's thought that they have descended from the same ancestor.There are about 480 000 water buffaloes in Iran.Most of these animals are kept in the south and northwest.All of the Iranian buffaloes are riverine (Naserian & Saremi 2007).
Milk, fat and protein yields are constantly monitored traits in herds integrated in milk test programs.The test-day milk yield (TDMY), defined as the total yield of a cow over a period of 24 h, replaces the milk yield in 305 days of lactation, as calculated by using formulas and extension factors (Tonhati et al. 2004).Buffalo milk is characterised by a high percentage of fat and protein, which could be used to increase milk's added value, following the current tendency of the worldwide milk markets.In the process of cheese making, fat and casein are the major milk solids incorporated into the final product; thus, these milk components are routinely used in many countries as a criterion to determine the milk price.In Western countries, most of the buffalo milk production is transformed to buffalo mozzarella cheese, an expensive fresh cheese, and some industries are paying for the milk quality (Aspilcueta-Borquis et al. 2010).
The use of the daily production with test-day models has several advantages over the traditional procedures of evaluating lactational records, such as the ability to account for environmental effects on each test-day and to model individual lactation curves (Schaeffer et al. 2000).Random regression models have been proposed as an alternative methodology for the analysis of longitudinal data or repeated measurement records.For these reasons, random regression models were recommended for the analyses of test-day models in dairy cattle (Schaeffer & Jamrozik 2008).The random regression models allow obtaining the breeding values for milk yield at any day of lactation in a continuous manner or for functions of the lactation curve, instead of finite dimensional models that only give punctual predictions of breeding values.Also, the random regression models allow estimates of covariances between coefficients of random functions or, equivalently, estimates of covariance functions.Moreover, random regression models provide estimates of breeding values with higher accuracies than the conventional finite dimensional models because all records available from lactation and short length lactation records can be used in the genetic evaluation (Jamrozik et al. 2000, Schaeffer et al. 2000).Few studies have reported about estimates of genetic parameters for milk yield in buffaloes.These estimates are generally restricted to total milk yield, with heritability estimates ranging from 0.14 to 0.40 (Rosati & Van Vleck 1998, Tonhati et al. 2000, Peeva 2002, Mourad & Khattab 2009), and TDMY, with estimates ranging from 0.01 to 0.24 (Hurtado-Lugo et al. 2006, Aspilcueta-Borquis et al. 2007).
The improvement through the selection of traits associated with milk yield and composition for milking buffaloes is dependent on the availability of reliable genetic parameter estimates for these traits.Therefore, the objective of the present study was to estimate genetic parameters for the TDMY, the milk fat percentage and the milk protein percentage in Iranian buffaloes using random regression models.

Material and methods
A total of 9 278 test-day production records obtained from 1 501 first lactation buffaloes on 414 herds in Iran between 1993 and 2009 were used for the analysis.The age of the animals at first calving varied from 24 to 60 months.For milk yield and milk fat percentage, contemporary groups were defined by the effects of herd, year, month and day of milk test and year-season of calving.For the milk protein percentage, contemporary groups were defined by herd, year, month and day of milk test.Only records of cows with at least four test days were kept in the analysis.The structure of the data set and pedigree information are shown in Table 1.Legendre polynomial functions were chosen to fit the lactation curves in the framework of a random regression test day model for estimating (co)variance components.Model specification and the choice of fixed effects to be considered in the model were based on the backward elimination method and variables which were significant at P<0.05 were included in the model.Also, various orders of polynomials were examined and the optimum set was selected based on the logarithm of the likelihood function at the point of conversion and total number of parameters to be estimated.The models of analysis included the fixed effects of contemporary groups and days in milk (DIM) as linear and quadratic covariables (only for milk yield).Additive genetic and permanent environmental effects were included in the models of analysis as random effects.The additive genetic and permanent environmental random effects were modelled by Legendre polynomials of DIM.For the milk yield, the Legendre polynomial functions of orders 3 and 5 were used to model the additive genetic and permanent environmental effects, respectively.For the milk fat percentage, the Legendre polynomial functions of orders 6 and 7 were used to model the additive genetic and permanent environmental effects, respectively.For the milk protein percentage, the Legendre polynomial functions of orders 5 and 6 were used to model the additive genetic and permanent environmental effects, respectively.The final model for the genetic analysis of the productive traits can be written as: where y is the vector of observations, b is the vector of the systematic effects and fixed regression coefficients, a is the vector of additive genetic random regression coefficients, c is the vector of permanent environmental random regression coefficients, e is the random vector of residual effects, and X, Z, and W are the incidence matrices of fixed effects and additive genetic and permanent environmental random effects, respectively.The model equation is based on the following assumptions: where k a and k c are (co)variance matrices between random regression coefficients for additive genetic and permanent environment effects, respectively, A is the relationship matrix, I is an identity matrix, ⊗ is the Kronecker product between matrices and R is a block diagonal matrix containing residual variances.Correlations between random regression coefficients for different effects were set to zero.The (co)variance components and genetic parameters for productive traits were estimated using the Average Information (AI) REML algorithm of the Wombat program (Meyer 2006).

Results and discussion
Overall phenotypic means, standard deviations and coefficients of variations for productive traits of Iranian buffaloes are presented in Table 2.The peak milk yield occurred on the 90th day of lactation.Also, the lowest milk yield was observed in the 270th day of lactation.Breda et al. (2010) reported that the peak period of milk production of Murrah buffaloes was observed around the 11th week of lactation.Hurtado-Lugo et al. (2009) observed the highest values of the TDMY towards the middle of the lactation in Colombia North buffaloes.The milk fat percentage was greater at the beginning of lactation but decreased and became lowest on the 60th day of lactation, and its peak occurred on the 270th day of lactation.In the milk protein percentage curves, the trend almost was similar to that observed for daily milk yields, with higher means in the early lactation, thereby suggesting a positive phenotypic association between milk yield and milk protein percentage.
Estimates of additive genetic, permanent environmental and phenotypic variances and heritabilities for TDMY, test-day records of fat (TDF) and protein (TDP) percentages are shown in Table 3.The additive genetic variance of the milk yield at DIM 30 was 0.41, which was decreased to the value of 0.16 in the following seven test-days and then reached the maximum estimated value of 0.72 at DIM 270.Estimates of the genetic and phenotypic variances for milk yield were within the range reported by White et al. (1999) and Olori et al. (1999) in dairy cattle.Tijani et al. (1999) used covariance functions developed by Gengler et al. (1999) to produce covariance functions for all days of lactation, and similar to the result of this study, they concluded that the genetic variance increased with DIM.A similar trend was observed for the additive genetic variances of TDP in this study.On the other hand, there was a decreasing trend for the estimates of the genetic variance for TDF.This variability in our results may be partially explained by the fact that estimates of the additive and permanent environmental variances for test-day yields in a random regression model depend on the submodel fitted (Olori et al. 1999).The other permanent environmental and phenotypic variances tend to follow oscillation at midlactation.The heritability estimates were higher for the milk yield at the beginning and at the end of lactation (Table 3).The reason for the highest heritability estimates at initial lactation might be due to the fact that the milk yield during the first test-days is critical to calf survival in terms of both volume and content and, as such, could have a large genetic component (Geetha et al. 2007).Hurtado-Lugo et al. (2006) and Tonhati et al. (2008) analysed TDMY records of milking buffaloes by finite models and the pattern of heritability estimates in their studies was similar to the current results.Similar to the results of this study, Geetha et al. (2007) and Tonhati et al. (2008) reported about the highest heritabilities of TDMY on the initial DIM.
Ignoring estimates of heritabilities for the milk yield in the first and final days of lactation, TDMY heritabilities were higher at DIM 90 and 120 due to higher genetic variance.The highest estimates of heritability for TDMY obtained in the third and fourth month of lactation by Jamrozik & Schaeffer (1997), Lidauer & Mäntysaari (1999) and Silvestre et al. (2005) in dairy cattle and it could be suggested that selection for milk yield should be made at the days of lactation 90 or 120.The selection could be anticipated by promoting a reduction in the generation interval.In the present study, the heritability for the milk yield (0.09 to 0.33) was higher than that reported by Hurtado-Lugo et al. (2006) for Murrah buffaloes in Colombia (0.01 to 0.20) and was consistent with the report of Chakraborty et al. (2010) and Breda et al. (2010) who reported that the heritability estimates for the TDMY were low to medium in Murrah buffaloes.The current heritability estimates ranged from 0.03 to 0.24 for TDF (Table 3).Similar to the results obtained for TDMY, the highest estimate of heritability for TDF observed on the 30th day of lactation due to higher genetic variation compared to other testdays.In this study, the highest estimates of heritabilities for TDP were obtained in the 120 days of lactation (0.27).In general, the heritability estimates for TDF and TDP were lower than those for the TDMY.Consistent with our results, Swalve (1995) reported that the estimates of heritabilities for test-day fat and protein traits were lower than the estimates of milk yield.
Estimates of genetic and permanent environmental correlations between different days of lactation and for each of the productive traits are separately shown in Table 4.The estimates of genetic correlations for milk yield ranged from −0.24 to 1.0 and varied from −0.85 to 0.92 for permanent environmental correlations.Genetic correlations of milk yield between different days of lactation were often higher than the permanent environmental correlations.The adjacent lactation days had generally high genetic correlations which varied from 0.73 (DIM 210-DIM 240) to 1.0 (generally at the midlactation).Negative genetic correlations for milk yield between the DIM 270 and DIM 90-180 were obtained in this study.Sesana et al. (2010) reported about unexpected negative genetic correlation estimates between TDMY records in the first weeks with records from the middle to the end of lactation.Except for DIM 30 and DIM 60, the genetic correlations were higher in the initial test-day yield and decreased when the test-day interval was increased.White et al. (1999) reported about higher estimates of genetic correlations between test-day yields throughout the lactation compared with the estimates obtained in this study.The pattern of variation of the genetic correlation in the current study has similarities with the results achieved by Brotherstone et al. (2000) and Kettunen et al. (2000) in dairy cattle.The permanent environmental correlation estimates for the milk yield were positive and high for adjacent days and permanent environmental correlations decreased from the maximum value of 0.92 on adjacent days to the value of −0.85 between days 30 and 150.In addition, the estimates of genetic correlations for milk fat and milk protein percentages throughout the lactation, ranged from −0.60 to 1.0 and −0.89 to 1.0, respectively.Also, the estimates of permanent environmental correlations for milk fat and milk protein percentages throughout the lactation varied from −0.83 to 0.79 and −0.91 to 1.0, respectively.The highest estimates of genetic correlations for TDMY, TDF and TDP were observed in adjacent test-days.The estimates of genetic correlations for productive traits in this study showed a lack of consistency between the estimates at the beginning and at the end of lactation.Also, similar to the results of this study, Brotherstone et al. (2000), Rekaya et al. (1999), Kettunen et al. (2000) and Bignardi et al. (2009) reported about negative genetic correlations for the productive traits in the course of the lactation in dairy cattle using random regression models.The post-calving cow's stress during the first days of lactation may influence the obtained because cows generally show an energy deficit

Table 2
Phenotypic mean, standard deviation and coefficient of variation for test-day milk yield, test-day milk fat percentage (Fat, %) and test-day milk protein percentage (Protein, %)