Estimates of genetic parameters for test day milk yields of a Holstein Friesian herd in Turkey with random regression models

Genetic parameters for test day milk yields of Holstein Friesian cows have been estimated using a random regression model (RRM). Data consisted of 1487 monthly test day milk yield records of cows calving between 1987 and 1993 in Sarmısaklı Farm, Turkey. Data were restricted to have at least 150d and at maximum 308d length of first lactations. Additive genetic and permanent environmental (co)variances were modeled with the same order polynomial regressions. Residual (measurement) error variance was assumed to be constant throughout lactation. The quadratic (k=3) order orthogonal polynomial regression was found to be sufficient. Heritability estimates for test day milk yields were high at the middle of the lactation and ranged from 0.07 to 0.32. Genetic correlations of milk yields between consecutive test days were high, but decreased as the interval between tests days increased. Genetic correlations ranged from 0.51 to 0.99. Residual error variance was estimated 13.77 kg.


Introduction
Test day records are expressions of a trait that changes over time (SWALVE, 1995a;Van der WERF et al., 1998).These records are used to predict total 305-d yields which are required to evaluate the additive genetic merit of sires and cows in traditional evaluation (ALI and SCHAEFFER, 1987).For the genetic evaluation of dairy cows using individual test day yields rather than total lactation production has a number of advantages.Test day models (TDM) allow: 1-Direct correction for environmental effects on the test day, 2-Better accounting for variation in number of tests recorded per animal, 3-Accounting for variation in the shape of the lactation curve (Van der WERF et al., 1998;SWALVE and GUO, 1999).A common approach to investigate genetic associations between test day yields is to consider every yield at each time period as a separate trait and then to estimate the genetic correlations between these traits.This approach has some disadvantages when large numbers of test day yields are considered.Biological interpretation of a large number of correlations is furthermore often difficult (VEERKAMP and THOMPSON, 1999).At the same time, heritability estimates for test day yields are usually less than for 305 day milk yields (SWALVE, 1995a;BAFFOUR-AWUAH et al., 1996;STRABEL and MISZTAL, 1999).ALI and SCHAEFFER (1987), obtained heritability estimates of first lactation 305 day milk yields of 0.28, 0.31 and 0.30 with three different TDM.PTAK and SCHAEFFER (1993) used repeatability models for first lactation test day milk yields.They assumed constant genetic variance and genetic correlations among test day records.SWALVE (1995b) estimated the heritability of 0.39 for 305 day milk yield while they were changed 0.18 to 0.36 from several test day milk records.These heritability estimates were higher than the estimates reported by ALI and SCHAEFFER (1987) with repeatability models.Alternatively, some researchers (ALI and SCHAEFFER, 1987;KIRKPATRICK et al., 1994;Van der WERF et al., 1998;GUO and SWALVE, 1997;JAMROZIK, 1997;MEYER and HILL, 1997;HORSTICK and DISTL, 2002;AMIN, 2003) have utilized covariances among all test day yields to improve the accuracy of predictions.The covariance structure is described by a covariance function, estimated by fitting a set of orthogonal polynomials or other defined covariables as random regressions on time of repeated records (OLORI et al., 1999b).A random regression model (RRM) allows different shapes of lactation curves for each cow by the inclusion of random regression coefficients for each animal (SCHAEFFER and DEKKERS, 1994).Using RRM, lactation curve for an individual cow is described by two sets of regressions on days in milk (DIM).Fixed regressions for all cows describe the general shape of lactation for cows belonging to the same subclass, for example regions, age at calving and season of calving, and the random parts of regressions for each cow describe the genetic deviation of individual regression from the fixed regressions, which allows each cow to have a genetically different shape of lactation curve (JAMROZIK et al., 1997;SWALVE and GUO, 1999).Random regression coefficients have been suggested by Henderson (1984), but SCHAEFFER andDEKKERS (1994) firstly improved this model to a RRM.In the literature many studies have been recently used RRM to estimate genetic parameters for production traits (Van der WERF et al., 1998;STRABEL and MISZTAL, 1999;VEERKAMP and THOMPSON, 1999;LIU et al., 2000;RÖHE et al., 2000;HORSTICK and DISTL, 2002).In these studies, third order RRM was used.Van der WERF et al. (1998) found high estimates of heritabilities at the periphery of the trajectory opposite to STRABEL and MISZTAL (1999); VEERKAMP andTHOMPSON (1999) andLIU et al. (2000).On the other hand, in most studies, different orders RRM were used.High estimates of heritability (0.59-0.40) for daily milk yield have been reported by JAMROZIK and SCHAEFFER (1997) and KETTUNEN et al. (1997, 1998) when using function of ALI and SCHAEFFER (1987).Conversely, the researches (REKAYA et al., 1999;OLORI et al., 1999b;LIU et al., 2000;POOL et al., 2000;ROMERO and CARABANO, 2003) have found less extreme heritability estimates, varied from 0.20 to 0.46, at the beginning and end of lactation from different RRM.Considering with previous studies, there was no any study on estimation of genetic parameters for test day milk yields of Turkish Holstein Friesian by using a RRM.In this study, the first goal was to decide the order of Legendre polynomial which gives the best fit of RRM with different orders of fit.The second aim was to estimate the additive genetic and permanent environmental (co)variances and heritability values for first lactation test day milk yields of Holstein Friesian cows using a RRM.

Data
The complete data set comprised monthly 1506 test day milk yields of Holstein Friesian cows obtained from Sarmısaklı Farm, in the Northwest region of Turkey.The cows were daughters of 56 sires and 119 dams, calved from 1987 through 1993.Total of 139 animals evaluated and there were 184 test date subclasses.Test day milk yields were recorded at successive 28-d periods throughout lactation and these periods were considered as monthly intervals (TD1-TD11).Test day months are used as the time variable rather than days in milk.Lactations were used to have at least 150d and at maximum 308d long of first lactations.

Model
The following RRM was used in the analysis: where ij y is the i th test day record of the j th cow, i TD is the fixed effect of test day (month) of recording i, jm β is the m th fixed regression coefficients, jm α and jm p are the m th random regression coefficients for additive genetic and permanent environmental effects of cow j, B k , A k and p k are the order of fitted fixed, random additive and random permanent regression coefficients, ij t is the i th standardized lactation month of the j th animal, m φ is the m th polynomial evaluated for the age ij t , and ij e is the random residual effect.
The matrix notation of the model can be written as, y= Xb+Za+Wp+e Where, vector b includes fixed regression coefficients jm β and TD effects, vector a and p include random regression coefficients for additive genetic and permanent environmental effects, e is the vector of residual effects and X, Z and W are the incidence matrices.The (co)variance structure for random effects in the model was defined as: G is the genetic covariance matrix of the random regression coefficients, assumed to be the same for all cows, A is the additive genetic relationship matrix among animals, 2 p σ is the variance of the permanent environmental effects, I is the identity matrix, R is the diagonal matrix of residual variance, ⊗ is the kronecker product function (SEARLE, 1982).Variance components were estimated by derivative-free REML (DFREML) using a RRM with the DXMRR statistical package (MEYER, 1997).Because of reducing the number of parameters to be estimated and reducing the dimension of the likelihood searches, residual variance was assumed to be constant throughout lactation.Additive genetic and permanent environmental (co)variances were modeled with the same order of polynomial regressions.Legendre polynomials were used because they are orthogonal, normalized and resulted in a better convergence and more accurate results as compared to conventional polynomials (KIRKPATRICK et al., 1990).Significant differences in the fit of Legendre polynomials with order from k=2 to k=6 were tested using a chi-square test ( 2χ ) of the likelihood.

Results
The maximum log likelihood values and changes in the log likelihoods from the models with different orders of fit and one measurement error class were presented in Table 1.Log likelihood values and their changes were decreased with increasing order of model.The changes in the log likelihood for quadratic, cubic and quartic model have been found to be significant (P<0.05).When comparisons of the models based on significant differences in the fit of models were tested using a 2 χ test of the likelihood, the quadratic polynomial had the largest changes (5.54%) of log likelihood values.The changes in likelihood for cubic, quartic and quintic models were only 2.79%, 2.26% and 0.67%, respectively (Table 1).Because the first three eigenvalues of the additive genetic coefficient matrix were large for all orders, the first three eigenvalues were shown in Table 2.In this table, proportion of each eigenvalue in total was also given to determine their importance.Since changes in the likelihood value of quadratic polynomial regression is the largest, this model was chosen as sufficient model to fit additive genetic and permanent environmental (co)variances with the fixed regression in this study.Moreover, the first, second and third eigenvalues obtained in the quadratic order of fit for the additive genetic covariance function were 13.55; 0.46 and 0.61E-03, respectively.The first eigenvalue of the coefficient matrix of the additive genetic covariance function accounted for about 97% of total eigenvalues for quadratic model.
On the contrary to the first one, second and third eigenvalues have negligible proportions of the of total eigenvalues that is related to variation in the additive genetic variance.
Eigenfunctions of additive genetic coefficient matrix of covariance function for the quadratic model are plotted in Figure 1.First eigenfunctions slightly decreased when the test days increased.Contrary to the first one, second and third eigenfunctions show an increasing pattern throughout test days.Heritabilities for test day yields were given in Table 3.Also, tendency of heritability estimates for test day milk yields from the quadratic order of fit was plotted in Figure 2. Heritabilities for test day milk yields were ranged from 0.07 to 0.32.It is obvious that the test day milk yields at beginning and end of lactation have lower heritability than the test days in the middle part.
Additive genetic and phenotypic correlations between test day milk yields were given in Table 3.The additive genetic correlations were higher than the phenotypic correlations.While the additive genetic correlations were changed from 0.51 to 0.99, the phenotypic correlations for test day milk yields varied from 0.11 to 0.60.Discussion Use of RRM makes it possible to study changes in test day records over time and a better understanding of lactation genetics (SWALVE and GUO, 1999;JAKOBSEN et al., 2001).But, feasibility of the RRM depends on the order of fit because of the computational difficulties.Therefore, order of fit should not exceed three (POOL et al., 2000), although regression models with higher orders estimate parameters more accurately (JAMROZIK et al., 1997;Van der WERF et al., 1998;STRABEL and MISZTAL, 1999).As a matter of fact, many studies show that third order polynomial RRM was found sufficient if complete lactations were used for parameter estimation (Van der WERF et al., 1998;POOL and MEUWISSEN, 1999;OLORI et al., 1999b;VEERKAMP and THOMPSON, 1999;KETTUNEN et al., 2000;STRABEL et al., 2003).In this study, to reduce computational problems, same order of fit for fixed and random effects was used.Results indicated that the maximum log likelihood values increased when the order of fit increased.The third order polynomial gave the best fit when considering the largest changes in log likelihood with the difference between number of parameters as a degrees of freedom versus table value of 2 χ .Moreover, the first three eigenvalues of the additive genetic coefficient matrix were greater in all models.These results support that the third order of fit is sufficient for modeling of test day records random regression coefficients as mentioned in literature (OLORI et al., 1999b;POOL and MEUWISSEN, 1999).
The eigenvalue is proportional to the amount of genetic variation in the population corresponding to related eigenfunction (KIRKPATRICK et al., 1990).For the third order model, the first eigenvalue of the coefficient matrix of the additive genetic covariance function accounted for about 97% of total eigenvalues.The large first eigenvalue can be concluded that selection on first eigenfunction will cause quick changes on the lactation curve.On the other hand, second and third eigenvalues accounted for about 3% of total eigenvalues and represent an unimportant proportion of the variation in additive genetic variance.Similar results were reported by Van der WERF et al. (1998) andAKBAŞ et al. (2004).Eigenfunction values from the first eigenfunctions were almost constant at early part of lactation, after that started to decrease at the second part of lactation.This means that factors are equally affect on most of the genetic variation of milk yield at early part of lactation while it was not the case for test days in the second part of lactation.
Since the first eigenfunction is related to largest eigenvalue, selection based on the factors will slightly decrease the milk yields for late test days.The second and third eigenfunctions show increasing pattern with increasing test days.This change reveal that a factor or factors have effects on milk yield different level in early and late stages of lactation (OLORI et al., 1999b).However, the second and third eigenvalues are small; selection based on the factors defined by third eigenfunctions not alters the milk yields.
Heritability estimates for test day milk yields were similar with the results reported by LIU et al. (2000); KETTUNEN et al. (2000) and DRUET et al. (2003).But, they were significantly higher as compared to results from STRABEL and MISZTAL (1999), LIDAUER and MANTYSAARI (1999), and lower as compared to results from BAFFOUR-AWUAH et al. (1996);JAMROZIK and SCHAEFFER (1997); KETTUNEN et al. (1998);OLORI et al. (1999b);POOL et al. (2000); ROMERO and CARABANO (2003).Shape of the heritability curve from the third order polynomial model resulted in decreasing patterns when the test days were increased.The heritability values at the beginning and the last part of lactation were lower than the estimates obtained for test day milk yields in the middle part as described by VEERKAMP and THOMPSON (1999);LIU et al. (2000) andPOOL et al. (2000).
The additive genetic correlations between milk yields obtained at different test days were higher than the phenotypic correlations.Genetic and phenotypic correlations between milk yields obtained at consecutive test days were positive and high but they were decreased as the interval between tests days increased.Furthermore, the phenotypic correlations of first test day milk yields with other test milk yields were relatively low.These results agreed with the results from previous works (KETTUNEN et al., 1998;KETTUNEN et al., 2000;VEERKAMP and THOMPSON, 1999).These results suggest that selection for increased milk yield at any test day, will effect positively on milk yields in late test days.
Changes in additive genetic variance over test days were similar as reported in the literature (VEERKAMP and THOMPSON, 1999;OLORI et al., 1999b;POOL et al., 2000) which was showing higher values in the middle of lactation.However OLORI et al. (1999b) found slightly higher variance estimates than the value estimated in this study for the middle part.On the other hand, opposite to additive genetic variance, phenotypic and permanent environmental variances were increased after second part of the lactation.This finding also was very similar to reported pattern (OLORI et al., 1999a, b;POOL et al., 2000).
In this study, error variance was assumed constant throughout lactation.This assumption causes residual variance in early lactation to be underestimated in RRM models but has no significant effect on the other variance components (OLORI et al., 1999a).The estimates of total variance and heritability for milk yield at early stage of lactation are affected by the constant error assumption.Also, this may partly explain the inconsistent heritability in early and late stages of lactation (OLORI et al., 1999a).Finally, to estimate the genetic structure of test day milk yields in Turkish Holstein Friesian, further investigations using larger data set seems to be required with varying orders of RRM and also accounting for heterogeneous measurement error variances in the analysis of test day milk yields.

Table 2
Eigenvalues with their relative proportions of coefficient matrix of the additive genetic covariance functions (Die ersten drei Eigenvalues der Koeffizientenmatrix und ihr unterschiedliches Verhältnis zur additiven genetischen Estimates of additive genetic, phenotypic and permanent environmental variances were ranged from 3.01 to 9.32; 4.78 to 21.84, and 26.40 to 38 61, respectively.Residual error variance was 13.77 kg 2 .Additive genetic, phenotypic, permanent environmental variances over test days were plotted in Figure3.Shapes of curves for the phenotypic and permanent environmental variances from the third order polynomial model match oscillatory patterns.Permanent environmental variances over test days show opposite changes with additive genetic and phenotypic variance.