ORIGINAL

Random regression models for milk, fat and protein in Colombian Buffaloes

Modelos de regresión aleatoria para leche, grasa y proteína en búfalos de Colombia

Naudin Hurtado-Lugo,^1,2 Ph.D, Humberto Tonhati,¹ Ph.D, Raul Aspilcuelta-Borquis,¹ Ph.D, Cruz Enríquez-Valencia,¹ M.Sc, Mario Cerón-Muñoz,²* Ph.D.

¹University of São Paulo, Faculty of Agrarian and Veterinary Sciences (FCAV), State Jaboticabal (UNESP), 14884900, SP, Brazil.
²Universidad de Antioquia, Facultad de Ciencias Agrarias, Grupo de investigación GaMMA, Carrera 75 No 65-87. Medellín, Colombia.

*Correspondence: cerongamma@gmail.com

Received: June 2013; Accepted: August 2014.

ABSTRACT

Objective. Covariance functions for additive genetic and permanent environmental effects and, subsequently, genetic parameters for test-day milk (MY), fat (FY) protein (PY) yields and mozzarella cheese (MP) in buffaloes from Colombia were estimate by using Random regression models (RRM) with Legendre polynomials (LP). Materials and Methods. Test-day records of MY, FY, PY and MP from 1884 first lactations of buffalo cows from 228 sires were analyzed. The animals belonged to 14 herds in Colombia between 1995 and 2011. Ten monthly classes of days in milk were considered for test-day yields. The contemporary groups were defined as herd-year-month of milk test-day. Random additive genetic, permanent environmental and residual effects were included in the model. Fixed effects included the contemporary group, linear and quadratic effects of age at calving, and the average lactation curve of the population, which was modeled by third-order LP. Random additive genetic and permanent environmental effects were estimated by RRM using third- to- sixth-order LP. Residual variances were modeled using homogeneous and heterogeneous structures. Results. The heritabilities for MY, FY, PY and MP ranged from 0.38 to 0.05, 0.67 to 0.11, 0.50 to 0.07 and 0.50 to 0.11, respectively. Conclusions. In general, the RRM are adequate to describe the genetic variation in test-day of MY, FY, PY and MP in Colombian buffaloes.

Key words: Cattle, genetics, zootechnics (Source: EuroVoc).

RESUMEN

Objetivo. Fueron estimadas funciones de covariancia para los efectos genéticos aditivos y de ambiente permanente, y posteriormente fueron estimados los parámetros genéticos para la producción de leche en el día del control (PDC) para leche (MY), grasa (FY), proteína (PY) y producción de queso mozzarella (MP), usando modelos de regresión aleatoria (RRM) mediante Polinomios Ortogonales de Legendre (LP) en búfalos de Colombia. Materiales y Métodos. Fueron analizados 1884 registros de PDC de primeras lactancias para MY, FY, PY y MP, hijas de 228 reproductores. Los registros procedían de 14 rebaños entre 1995 y 2011. Para las PDC fueron consideradas 10 clases mensuales de días en lactancia. Los grupos contemporáneos fueron definidos como rebaño-año-mes de control. Fueron incluidos en el modelo efectos aleatorios genético aditivo, de ambiente permanente y residual. Los efectos fijos incluidos fueron grupo contemporáneo, efecto lineal y cuadrático de edad al parto, y curva promedia de la lactancia de la población, la cual fue modelada por un LP de tercera orden. Los efectos aleatorios genético aditivo y de ambiente permanente fueron estimados por RRM mediante LP de tercero a sexto orden. Las variancias residuales fueron modeladas usando estructuras homogéneas y heterogéneas. Resultados. Las heredabilidades para MY, FY, PY y MP variaron de 0.38 a 0.05, 0.67 a 0.11, 0.50 a 0.07 y 0.50 a 0.11, respectivamente. Conclusiones. En general, los RRM son adecuados para describir la variación genética para la PDC de MY, FY, PF y MP en búfalos de Colombia.

Palabras clave: bovino, genética, zootecnia (Fuente: EuroVoc).

INTRODUCTION

Buffalo milk is mainly used for the production of milk derivatives due to the composition and physicochemical quality of its components (1). The nutritional quality and economic value of mozzarella cheese does not only depend on the quantity of milk produced, but also on its relationship with fat and protein (1,2).

Random regression models (RRM) can be applied to test-day milk yield as an alternative to standard procedures used for the genetic evaluation of longitudinal traits in dairy cattle. According to El Faro et al (3), in dairy cattle of different breeds, heritabilities of 0.15 to 0.48, 0.13 to 0.28 and 0.16 to 0.28 were estimated for milk (MY), fat (FY) and protein (PY) yield, respectively, using RRM (4-7).

Studies estimating heritabilities for mozzarella cheese production (MP) with RRM are scarce. The objetive of this study was to estimate covariance functions for additive genetic and permanent environmental effects and, subsequently, genetic parameters for MY, FY, PY and MP in buffaloes from Colombia using RRM with Legendre polynomials (LP).

MATERIALS AND METHODS

Description of the data. In the study, test-day MY, FY, PY and MP from 1884 first lactations of buffalo cows aged 26 to 49 months, from 228 sires, were analyzed. The animals belonged to 14 herds that comprise the database of Faculty of Agricultural Sciences, University of Antioquia, Colombia and Colombian Buffaloes Breeders Association. Calvings occurred between 1995 and 2011. Lactations were considered when the first control was between 5 and 45 d. Test-day yields were divided into monthly classes of lactation (1 to 10 classes) and animals with at least three test-day records were included in the analysis. The contemporary groups (CG) were defined as herd-year-month of milk test, with the restriction that each group should contain at least three animals. A pedigree file containing 7.090 animals was used in all analyses.

MP was estimated using the mozzarella index proposed by Borghese and Mazzi (8), which is based on MY, %FY and %PY. This index is adopted in genetic and animal breeding programs in Italy and is calculated as follows:
MP(kg)=MY*{[(3.5*%P)+(1.23*%F)-0.88]/100}.

Statistical method. Single-trait RRM were used for the analysis of MY, FY, PY, and MP. Random additive genetic effects, permanent environmental effects of the animal, and residual effects were included in all models. The fixed effects were CG, linear and quadratic effects of the covariate age of dam at calving, and the average lactation curve of the population modeled by a third-order orthogonal polynomial. The matrix representation of the model is given as:

Where
y = vector of observations;
b = vector of fixed effects (CG and age of dam at calving and average curve of the population covariates);
a = vector of solutions for additive genetic random regression coefficients;
ap = vector of solutions for permanent environmental random regression coefficients;
e = vector of residual effects, and
X,Z,W = incidence matrices for fixed,random additive genetic and permanent environmental effects, respectively.

The dimension of vector “a” consists of k_a x N_a coefficients, where k_a corresponds to the order of the polynomial and N_a is the number of animals in the relationship matrix. The dimension of vector ap consists of k_ap x N_d coefficients, where k_ap corresponds to the order of the polynomial and N_d is the number of animals with records.

For analysis, it was assumed that the records were distributed with mean Xβ and the random effects were assumed to be normally distributed with:

Where
K_a and K_ap are covariance matrices between additive genetic and permanent environmental random regression coefficients, respectively;
A is the relationship matrix between individuals;
I_Nd is the identity matrix of dimension N_d;
is the Kronecker product between matrices,
and R is a block diagonal matrix containing residual variances. It was assumed that the residuals are independent.

The random additive genetic and permanent environmental effects were modeled using third- to sixth-order LP. A homogeneous (R) and a heterogeneous (r) structure with 10 classes (each month was considered to be a different class) were adopted for residual variances. In addition, heterogeneous variances with smaller classes were also considered, which were divided according to similarity between variances based on the 10 residual classes. Thus, test-day records for MY were divided into four classes (1^st, 2^nd to 4^th, 5^th to 8^th and 9^th to 10^th month). The residual variance structure for FY was defined by four classes (1^st, 2^nd, 3^rd to 8^th and 9^th to 10^th month). For PY, residual variance was modeled using four classes (1^st, 2^nd to 7^th, 8^th and 9^th to 10^th month). Residual variance for MP was modeled using four classes (1^st to 2^nd, 3^rd to 7^th, 8^th and 9^th to 10^th month).

The general RRM is given as LEGk_a.k_ap_R or LEGk_a.k_ap_r, where k_a and k_ap correspond to the order of the covariance function for additive genetic and permanent environmental effects, respectively, and R or r is the structure of residual variances. The covariance functions were estimated by the restricted maximum likelihood method using the WOMBAT program (9).

The different models were compared using the logarithm of the likelihood function (log L), the likelihood ratio test (LRT) at 1% probability, restricted maximum likelihood Akaike's (AIC) and Schwarz Bayesian (BIC) Information Criteria, and examination of variances and correlations estimated for the traits. The information criteria can be described as: and , where p is the number of parameters estimated, N is the number of data, r(X) is the rank of the incidence matrix of fixed effects in the model, and log L is the logarithm of the restricted maximum likelihood function.

RESULTS

The mean test-day MY data followed a typical lactation curve for buffaloes, starting with 3.69 kg and increasing until peak yield in the third month (4.14 kg), followed by a decrease until the end of lactation (2.34 kg). The mean test-day FY, PY and MP were 0.24, 0.14 and 0.75 kg, respectively, with a standard deviation of 0.11, 0.07 and 0.33 kg and coefficient of variation of 46.58, 45.89 and 44.41% (Table 1).

Milk component yields remained constant during early lactation until the third month (0.27 kg for FY and 0.17 kg for PY). Production declined after the fourth month with increasing number of days in milk until the end of lactation (0.15 kg for FY and 0.08 kg for PY). In contrast, MP increased from early lactation (0.89 kg) to the second month (0.90 kg), followed by a decline after the third month to the end of the lactation (0.46 kg). In general, a greater variation in FY, PY and MP was observed at the end of lactation.

The LEG3,3_4 model should be used for MY, PY and MP, and the LEG3,3_6 or LEG3,3_4 models should be used for FY, with the biological interpretation of this trait being easier with the LEG3,3_4 model because of the smaller number of parameters (Table 2).

RRM with LP require the definition of the most appropriate order for each random effect included in the model. For MY, FY, PY and MP (Table 3), improvement of the log L (p>0.01 by the LRT) and AIC criteria was generally observed when the order of fit was increased from three to four for additive genetic variance combined with an order of six for permanent environmental variance.

Table 4 shows the eigenvalues associated with the matrix of random regression coefficients for additive genetic and permanent environmental effects obtained for the best models for MY, FY, PY and MP according to the BIC. As can be seen, lower-order RRM permit to model most part of the additive genetic variation in MY, FY, PY and MP, with a fourth coefficient being associated with an eigenvalue of zero, in contrast to the permanent environmental effect. The first eigenvalue analyzed was responsible for at least 62% of the data variation in MRA adjusted to the traits studied.

The first three eigenvalues explained most of the additive genetic and permanent environmental variability (over 90%). On the basis of the results of the eigenvalues, the dimension of the three random effects can be reduced without the loss of information, disagreeing with the log L and AIC criteria for additive genetic and permanent environmental effects and agreeing with the BIC criterion for the additive genetic effect.

The estimates of additive genetic, permanent environmental, residual and phenotypic variances for test-day MY (LEG3,5_4), FY (LEG3,4_4), PY (LEG3,4_4), and MP (LEG3,4_4) are shown in figure 1. In general, the additive variances for the traits were higher in the first month, decreased until the eighth month, and then increased in the ninth and tenth month.

The estimates of permanent environmental variance for MY were practically constant across lactation and were higher than the additive variances. Residual variance was declined from the fifth to the tenth month (Figure 1). For PY, the additive genetic variances showed the same trend as the phenotypic variances, with the observation of higher values at the end of lactation. For FY and MP, the additive genetic variances showed the same trend as the phenotypic and residual variances and were higher at the end of lactation. These higher estimates for FY, PY and MP might be due to the number of records found or to the high fat and protein content at the end of lactation (Figure 1).

The heritabilities for MY, FY, PY and MP ranged from 0.38 to 0.05, 0.67 to 0.11, 0.50 to 0.07 and 0.50 to 0.08, respectively (Figure 2).

The estimates of genetic and phenotypic correlations for test-day MY, FY, PY and MP are shown in figure 3. The genetic correlations for these traits ranged from -0.56 to 0.96. The phenotypic correlations of the traits studied were lower as the interval between test days increased (Figure 3). In general, the phenotypic correlations were of lower magnitude than the genetic correlations.

DISCUSSION

Evaluation of the best fit for residual variance (Table 2) showed increases of log L (p<0.01), which were significant by the LRT, with increasing number of heterogeneous classes. According to the criteria used for the assessment of goodness-of-fit (AIC and BIC), the model considering homogeneity of residual variances was found to be inappropiate. This fact indicates a distinct behavior of residual variances across lactation and the need to consider heterogeneous variance structures for residuals effects. The heterogeneity of residual variances can be attributed to factors such as stage of pregnancy, body condition, duration of the lactation interval and management conditions, since these factors are not easily incorporated into models analyses because of the lack of information (6-7).

The changes in log L were of small magnitude and nonsignificant by the LRT (p>0.01), indicating that four residual classes would be sufficient to obtain a good fit of residual variance for MY, PY, FY, and MP. This would avoid the use of superparameterized models which generally present parameter estimation problems (7).

The BIC, which is more rigorous because of parameterization, indicated that LEG3,5_4, LEG3,4_4, LEG3,4_4 and LEG3,4_4 were the best fitted models for MY, FY, PY and MP, respectively. Therefore, these models would be more adequate to describe the biological variation in yield traits across lactation.

Specifically for PY and FY, the increases seen in late lactation may be related mainly to the small number of records available for this period (Table 1) or may have been due to computing artifacts rather than biology (8). In addition, fitting RRM at the beginning and end of lactation is difficult, resulting in parameter estimation problems (3-4). The results might be explained by the fact that buffalo breeding programs are only starting in Colombia and the population therefore presents small selection rates for the traits studied. In addition, the low yields in early and late lactation are related to the lack of feeding and management strategies during these periods.

The trends of additive and permanent environmental variances across lactation obtained for the traits studied (MY, PY and FY) are comparable to those reported by Jamrozik and Schaeffer (4) and Silvestre et al (5), who estimated higher additive and permanent environmental variances at the beginning and end of lactation for MY, PY and FY in dairy cows.

The heritabilities for MY ranged from 0.38 (first month) to 0.05 (seventh month). Sesana et al (7) found higher estimates, particularly at the extremes of the lactation curve, a finding that might be attributed to the fact that the authors used weekly records. The heritabilities for FY ranged from 0.67 (first month) to 0.11 (seventh month), whereas the heritabilities for PY ranged from 0.50 (first month) to 0.07 (ninth month). Similar results have been reported for dairy cattle using RRM, with high estimates for FY and PY at the beginning and end of lactation (4-6, 9-11). Aspilcueta et al (6) reported heritabilities of lower magnitude using finite dimension models for monthly records of FY and PY in dairy buffaloes. The heritabilities for MP ranged from 0.50 (first month) to 0.08 (ninth month). Rosati and Van Vleck (12) reported a heritability for MP of 0.13 for 270-day lactation in the Italian buffalo population. In general, the heritability estimates for MY, PY, FY and MP showed a similar trend across lactation, with the observation of a marked increase at the end of this period (Figure 2). Similar trends have been reported for dairy cattle (4-5).

For MY, the fractions of phenotypic variance due to permanent environmental effects were higher than the heritability estimates across lactation. The fractions of phenotypic variance for FY and PY were higher in the first month and between the fourth and seventh month of lactation, respectively. For MP, the fraction of phenotypic variance due to permanent environmental effects remained constant across lactation.

The estimates were higher between adjacent test days (early lactation) and decreased with increasing distance between test days, especially between tests performed at the beginning and end of lactation. Negative and near to zero genetic correlations between milk tests performed at the beginning and end of lactation have been reported for dairy cattle (13-14). Sesana et al (7) and Breda et al (15) reported similar trends in genetic correlations for buffaloes.

The correlations, variances and heritabilities obtained for the traits studied are probably due to the difficulty in modeling yields in early lactation, a phase during which the buffalo cow suffers from postpartum stress, as well as negative energy balance and nutritional deficiencies. In addition, genetic breeding programs of dairy buffaloes are only starting in Colombia and the herds and environment where the animals are raised, as well as the low selection rates and small number of observations may influence the early and late stages of lactation. This suggests the use of shorter lactation periods for buffalo populations in which selection programs are only started. De Roos et al (10) suggested the need to model fixed curves since the estimates of genetic variances at the beginning and end of lactation are increased if fixed herd curves are not modeled independently.

In conclusion, the RRM using LP are appropiate to describe the genetic variation in test-day MY, FY, PY and MP of buffaloes. Problems with the estimation of genetic parameters were observed at the beginning and at the end of the lactation curve, indicating that the traits studied may be affected by the number of available records at the extremes of the curve or by the high content of fat and protein, or may have been due to computing artifacts rather than biology (8). Residual variances varied across lactation and can be modeled by forming heterogeneous classes using low-order LP, explaining the changes in genetic variation according to days in milk.

Acknowledgments

This research had financial support from the Colombian Ministry of Agriculture and Rural Development, Colombian Cattle Breeders Federation (FEDEGAN), Colombian Buffaloes Breeders Association, Codi sostenibilidad E01808 and Fundação de Apoio à Pesquisa do Estado de São Paulo (FAPESP- Part of doctoral thesis of first author: Process N°2009/53773-1).

REFERENCES

1. Tonhati H, Lima A, Lanna D, Camargo G, Baldi F, Albuquerque L, Montrezor J. Milk fatty acid characterization and genetic parameter estimates for milk conjugated linoleic acid in buffaloes. J Dairy Res 2011; 78(2):178-183.         [ Links ]

2. Seno L O, Cardoso V, Tonhati H. Valores econômicos para as características de produção de leite de búfalas no Estado de São Paulo. R Bras Zootec 2007; 36(6):2016-2022.         [ Links ]

3. El Faro L, Cardoso Vera L, Albuquerque LG. Variance component estimates applying random regression models for test-day milk yield in Caracu heifers (Bos taurus Artiodactyla, Bovidae). Genet Mol Biol 2008; 31(3):665-673.         [ Links ]

4. Jamrozik J, Schaeffer LR. Estimates of Genetic Parameters for a Test Day Model with Random Regressions for Yield Traits of First Lactation Holsteins. J Dairy Sci 1997; 80(4):762-770.         [ Links ]

5. Silvestre A, Petim-Batista F, Colaço J. Genetic parameter estimates of Portuguese dairy cows for milk, fat and protein using a spline test-day model. J Dairy Sci 2005; 88(3):1225-1230.         [ Links ]

6. Aspilcueta RB, Sesana RC, Muñoz-Berrocal M, Seno LO, Bignardi AB, El Faro L, et al. Genetic parameters for milk, fat and protein yields in Murrah buffaloes (Bubalus bubalis Artiodactyla, Bovidae). Genet Mol Biol 2010; 33(1):71-77.         [ Links ]

7. Sesana RC, Bignardi AB, Aspilcueta RA, El Faro L, Baldi F, Albuquerque LG, Tonhati H. Random regression models to estimate genetic parameters for test-day milk yield in Brazilian Murrah buffaloes. J Anim Breed Genet 2010; 127(5):369-376.         [ Links ]

8. Martínez CA, Elzo M, Manrique C, Grajales LF, Jimenez A. Random Regression Models for Estimation of Covariance Functions, Genetic Parameters and Prediction of Breeding Values for Rib Eye Area in a Colombian Bos indicus-Bos taurus Multibreed Cattle Population. Rev Colomb Estad 2012; 35(2):309-330.         [ Links ]

9. Borghese A, Mazzi M. Buffalo Production And Research. Rome: Editorial Food and Agriculture Organization of The United Nations; 2005.         [ Links ]

10. Meyer K WOMBAT. A tool for mixed model analyses in quantitative genetics by REML, J Zhejiang Uni 2007; 8:815-821.         [ Links ]

11. De Roos APW, Harbers AGF, de Jong G. Random herd curves in a test-day model for milk, fat, and protein production of dairy cattle in The Netherlands. J Dairy Sci 2004; 87(8): 2693-2701.         [ Links ]

12. Rosati A, Van Vleck LD. Estimation of genetic parameters for milk, fat, protein and mozzarella cheese production for the Italian river buffalo Bubalus bubalis population. Livest Prod Sci 2002; 74 (2):203-217.         [ Links ]

13. Munera BOD, Herrera RAC, Cerón-Muñoz MF. Componentes de (co) varianza y parámetros genéticos para producción y valor económico de la leche a través de modelos de regresión aleatoria en hembras Holstein de primera lactancia. Livest Res Rural Dev 2013;25(12): art 208.         [ Links ]

14. Herrera RAC, Munera OD, Cerón-Muñoz MF. Variance components and genetic parameters for milk production of Holstein cattle in Antioquia (Colombia) using random regression models. Rev Col Cien Pec 2013; 26:90-97.         [ Links ]

15. Breda F , Albuquerque L, Euclydes RF, Bignardi A B. Baldi F, Torres RA, Barbosa L, Tonhati H. Estimation of genetic parameters for milk yield in Murrah buffaloes by bayesian inference. J Dairy Sci 2010; 93(2):784-791.         [ Links ]