Random regression models for genetic evaluation of milk yield in the second lactation of tropical dairy goats

B C F Nogueira¹, H R Oliveira^2,3, N O Souza¹, V S Junqueira¹, M T Rodrigues¹, F F Silva¹ and L F Brito²

¹ Departament of Animal Sciences, Universidade Federal de Viçosa, Viçosa, 36570-900, MG, Brazil
hrojasde@purdue.edu
² Department of Animal Sciences, Purdue University, 270 Russell Street, West Lafayette, 47907, IN, USA
³ Department of Animal Biosciences, University of Guelph, 50 Stone Road, Guelph, N1G 2W1, ON, Canada

Abstract

Dairy goats breeding programs have mainly focused on improving milk yield in the first lactation. However, genetic analysis of second lactation is needed to increase the genetic gain over lactations. Given the lack of information about second lactation, test-day models such as random regression (RRM) are required to ensure better selection accuracies. In this context, we aimed to identify the most suitable RRM for genetic evaluation of second lactation of tropical dairy goats, as well as to estimate genetic parameters based on the chosen RRM. A total of 15,509 milk yield test-day records from the second lactation of 529 goats were used. The pedigree file was composed by 5,511 animals. Statistical analyses were performed using RRM considering different orders of Legendre polynomials to describe the curves for population mean, additive genetic and permanent environment effects. In addition, several classes of residual variances were tested. The best RRM were chosen based on the log likelihood, likelihood ratio test, number of parameters in the model, and modified Akaike and Bayesian Information criteria. The number of records, average, median, standard deviation, and the maximum and minimum values for the milk yield at test-day were 15,509, 2.04, 1.05, 2.0, 6.6 and 0.1, respectively. Among all 320 tested models, the RRM that used sixth order for the fixed and permanent environment effects, third order for the additive genetic effect, and four classes of residual variances was chosen as the most suitable one for genetic evaluation. Additive variances and heritabilities estimated using the best RRM were higher in the middle of lactation. Heritabilities obtained through the lactation curve were moderate to high (they ranged from 0.10 to 0.29), indicating that it is possible to promote genetic progress for milk yield in the second lactation of tropical dairy goats. Moreover, results reported in this paper suggest that it is possible to change the format of the lactation curve in the second lactation.

Keywords: Alpine, animal breeding, longitudinal, test-day

Introduction

The goats’ sector has become economically promising for Brazil. Currently, more than 9.78 million goats are being raised in the country, which shows the great importance of this sector in the national agriculture (Souza et al 2017). Due to the large number of animals, Brazil is one of the largest goat milk producers in America (FAO 2011). Among the several benefits of goat milk are its hypoallergenic properties, good digestibility, and high nutritional value due to biological proteins, essential fat acids, minerals and vitamins (Haenlein 2004). In order to increase the amount of milk produced per animal, investments have been made in the goats’ sector, such as the creation of national breeding programs to perform genetic evaluations (e.g., CAPRAGENE 2019). In summary, these programs have mainly focused on improving milk yield in the first lactation, as data from subsequent lactations are not included in the evaluations (CAPRAGENE 2019). However, selecting animals to become parents of the next generation based solely on the first lactation may not be ideal, as genetically superior animals in the first lactation may have reduced genetic merit in subsequent lactations. This fact is due to the genetic correlation between lactations, which is usually different from one (e.g., Oliveira et al 2016; Brito et al 2018). Genetic correlations between first and subsequent lactations might not be high because goats in first lactation are still developing their mammary gland (Zeng et al 1997). Therefore, the lack of information about genetic evaluations of dairy goats in second lactation is harmful, as this information is needed to increase the genetic gain over lactations.

The first step to include new traits (e.g., lactations) in a breeding program is to verify if they are heritable. Thus, the success of breeding programs relies on the heritability of the traits included in the selection indexes, and on the accuracy of estimated breeding values. In this context, random regression models (RRM) have been used to estimate more accurate genetic parameters and breeding values for longitudinal traits, such as milk yield (Meyer 2005; Oliveira et al 2019). Among the advantages of RRM is the fact that they allow the estimation of genetic parameters and breeding values along the lactation (Oliveira et al 2019), because they use statistical functions to model the trait over time. The most popular RRM functions used to describe the traits over time are the Legendre orthogonal polynomials (Kirkpatrick et al 1990). Therefore, we aimed to: 1) identify the most suitable order of Legendre orthogonal polynomials to model the fixed and random curves in the RRM used for genetic evaluation of milk yield in the second lactation of tropical dairy goats, and 2) estimate the covariance components and genetic parameters using the chosen RRM.

Materials and methods

The dataset used in this study was obtained from the goats’ sector of the Animal Science department of the Federal University of Viçosa, Viçosa, Minas Gerais, Brazil. The animals were considered healthy and kept in collective padlocks under free housing system, receiving a diet based on corn silage and a concentrated mixture offered according to their nutritional needs. The climate in this region is classified as mesothermic, in which summers are usually hot and rainy, and winters are cold and dry.

A total of 15,509 milk yield test-day records from 529 goats in second lactation were used, and the pedigree file contained 5,511 animals, from 15 different generations. The litter size per dam ranged from 1 to 3, and the total number of kids per dam ranged from 1 to 14. The total number of kids per father ranged from 1 to 298. The maximum inbreeding coefficient estimated for an animal was 0.3750, and the average inbreeding coefficient in the herd was 0.0109. The phenotypic recording was performed every Saturday twice a day (at 6am and 2pm), using mechanical milking. The milk yield at test-day (TDMY) represented the sum of amount of milk in these two records. As RRM do not require the use of complete lactations (Oliveira et al., 2019), all animals that had at least 5 TDMY collected between 5 and 290 days after goat parturition were used in the analysis.

The influence of environmental effects that could affect TDMY was evaluated using the Generalized Linear Model procedure (Proc GLM) of the Statistic Analysis System (SAS 9.0) software. The significative effects were: number of kids (i.e., litter size; ranging from 1 to 3 kids), year of birth (ranging from 2001 to 2014), season (taking 1 for the period from March to September; and 2 to the other months), and breed composition between Alpine and Saanen (Alpine breed as base proportion). Age of the goat at parturition (AGP), considering the linear and quadratic effects, was also significative and included as covariate in the RRM. Therefore, contemporary groups were formed by the number of kids, year and season of birth, and breed composition. The number of records by contemporary group ranged from 4 to 70, and the total number of contemporary groups was 1,514. Different orders of Legendre orthogonal polynomials were tested for days in lactation (DIL).

Analyses were performed using single-trait RRM, in which additive and permanent environment effects were included as random effects. The RRM used are described as:

in which y_ij corresponds to the TDMY j of the goat i; and EF_i corresponds to the group of fixed effects, i.e., year, season, number of kids, genetic group, and age of the goat at parturition. The are the m fixed regression coefficients related to TDMY for modelling the average curve of the population (i.e., fixed curve);α_im and Y_im are the m random regression coefficients of additive and permanent environment effects for the goat i, respectively; K_b, K_α, K_y are the orders of the Legendre orthogonal polynomials; ø_m (t_ij) is the Legendre polynomial for the DIL of the goat i in the test-day j (t_ij), standardized for the interval of -1 to 1 as suggested by Kirkpatrick et al (1990); and Ꜫ_ij is the residual effect.

A total of 320 different RRM were tested in this study, which had combinations of distinct orders for the polynomials of fixed, random genetic, and permanent environmental curves, besides a different number of classes for residual variance. In summary, four different orders (i.e., third, fourth, fifth and sixth orders) were tested for the fixed, additive and permanent environment curves. In addition, five different number of classes were tested for the residual variance (i.e., 1, 2, 3, 4 or 5 classes, equally spaced). After herein, the notation of the tested RRM follows the pattern: F_A_P_H_, with F_ referring to the order of the polynomial curve to model the population mean, the additive genetic (A_) and permanent environmental (P_) effects, and (H_) is the number of classes of residual variances. For instance, a model with regression of mean curve of order 3, additive genetic of order 4, and permanent environment of order 5, and that had two classes of heterogeneous residual variance was coded as F3A4P5H2.

The variance components were estimated via restricted maximum likelihood (REML), using the WOMBAT software (Meyer 2007). The criteria used to choose the best RRM were the log likelihood (LogL); likelihood ratio test (LRT); number of parameters (NP); modified Akaike information criterion (AICm); and modified Bayes information criterion (BICm). The AIC and BIC criteria were modified by multiplying the traditional AIC and BIC by - 1/2, in order to make their scale similar to the LogL, facilitating the comparison between RRM.

Thus, the AICm and BICm used in this study are defined as:

where p is the number of parameters, N is the number of observations, and r is the rank of the fixed effects matrix. Regarding the LRT criteria, the comparison between complete and nested models were made at 5% probability, assuming the difference in the number of parameters between models as the number of degrees of freedom.

Once identified the best RRM, the heritability estimates throughout the lactation curve were obtained as:

where h²_j is the heritability estimate at test-day j, and are the additive genetic, permanent environmental, and residual variance estimates, respectively.

Results

The number of records, average, median, standard deviation, and the maximum and minimum values for the TDMY were 15,509, 2.04, 1.05, 2.0, 6.6 and 0.1, respectively. Among the 320 RRM evaluated, the complete model (F6A6P6H5) and the five best models are presented in Table 1. The LogL, AICm and BICm estimated for all RRM are shown in the Fig. 1, 2 and 3, respectively.

Figure 1. A 3D plot showing the log likelihood (LogL) of random regression models used for genetic evaluation of milk yield in the second lactation of dairy goats, using different orders of Legendre polynomials. Notation of the models follows the pattern: fixed (F), additive (A) and permanent environment (P) effects, and number of classes of residual variances (H)

Figure 2. A 3D plot showing the modified Akaike information criterion (AIC) of random regression models used for genetic evaluation of milk yield in the second lactation of dairy goats, using different orders of Legendre polynomials. Notation of the models follows the pattern: fixed (F), additive (A) and permanent environment (P) effects, and number of classes of residual variances (H)

Figure 3. A 3D plot showing the modified Bayes information criterion (BIC) of random regression models used for genetic evaluation of milk yield in the second lactation of dairy goats, using different orders of Legendre polynomials. Notation of the models follows the pattern: fixed (F), additive (A) and permanent environment (P) effects, and number of classes of residual variances (H)

According to the LTR, all five models were not statistically different from the complete model (Table 1). Increasing the order of genetic and permanent environment effects increased the LogL, AICm and BICm (Fig. 1, 2 and 3). The F6A3P6H4 was considered the most suitable RRM for genetic evaluations of second lactation of dairy goats, because it had smaller number of parameters (31), and/or higher LogL, AICm and BICm compared to the other alternative RRM. Variance components over DIL for milk yield in the second lactation, estimated using the most suitable RRM (i.e., F6A3P6H4), are shown in Fig. 4.

Figure 4. Additive, permanent environment, residual and phenotypic variances over days in lactation for milk yield in the second lactation of dairy goats, estimated using the best random regression model (i.e., F6A3P6H4; which used 6th order Legendre polynomials for the fixed and permanent environment curves, 3th order Legendre polynomials for the additive genetic curve, and four different classes of residual variances)

Additive genetic variances estimated in this study were higher in the middle of lactation. During all DIL, the permanent environment variances were higher than the additive genetic variances (Fig. 4). Heritabilities estimated using the referred variance components are shown in Fig. 5. The heritability estimates ranged from 0.10 to 0.29. The highest heritabilities were observed between the 220th and 240th DIL, and the lowest heritability estimates were observed during the beginning and end of lactation.

Figure 5. Heritabilities (h2) over days in lactation for milk yield in the second lactation of dairy goats, estimated using the best random regression model (i.e., F6A3P6H4; which used 6th order Legendre polynomials for the fixed and permanent environment curves, 3th order Legendre polynomials for the additive genetic curve, and four different classes of residual variances)

Discussion

In our knowledge, there are no studies in the literature that address exclusively the second lactation of dairy goats. However, studies have compared the use of different lactations in the genetic evaluations. In this context, Silva et al (2014), comparing the first four lactations of Alpine dairy goats, concluded that breeding programs should not evaluate solely the first lactation, because there is a large re-ranking of animals when different lactations are analyzed together. In a study using first and second lactations of Alpine and Saanen breeds, Irano et al (2012) concluded that there is a decrease in genetic variability from the first to the second lactation. In addition, Menéndez-Buxadera et al (2010) suggested that the genetic correlations between different DIL are larger and more stable in the first lactation than in the second one, even though the amount of milk yield is bigger in the second lactation. Thus, it is paramount to define the appropriate RRM to genetically evaluate the second lactation of dairy goats.

The average, maximum and minimal values for the TDMY in the second lactation found in this study were similar to those reported by Silva et al (2013). Regarding to the model comparison (Table 1, and Figures 1, 2 and 3), higher values for LogL are always expected for more parameterized models. However, it is important to use other comparison criteria that provide a choice of more parsimonious models (Brito et al 2018). The AIC and BIC allow the comparison between nonaligned models, and they penalize models with higher number of parameters (Nunez-Antón and Zimmerman 2000). In this context, the BIC criterion is even more severe than the AIC (Nunez-Antón and Zimmerman 2000). Combining all model selection criteria, the F6A3P6H4 model was chosen as the most suitable RRM for genetic evaluation of second lactation of tropical dairy goats.

The pattern of additive genetic and permanent environment variances estimated in this study (Fig. 4) were similar to those reported by Silva et al (2013), studying TDMY in Alpine goats. On the other hand, Menezes et al (2011) found higher additive genetic variances at the beginning of lactation, and higher permanent environment variances at the end of lactation for Saanen goats. Residual variances were heterogeneous in this study, and they decreased along the lactation period. This fact may be related to the scale effect, as TDMY are larger in the beginning of lactation (closer to the lactation peak) than in the end of lactation. More details about the scale effect in goats can be obtained from Silva et al (2014), and Oliveira et al (2016). Phenotypic variances estimated in this study were higher between 45th and 105th DIL, which may be related to the great differences in the lactation peak and persistency among animals in this herd. The slight increase in phenotypic variances observed in the end of lactation is probably due to the lower number of records available during this period.

The higher heritability estimates observed between the 220th and 240th DIL suggest that selection during this period may bring desirable genetic gain. In addition, the low heritability estimates were observed in the beginning and end of lactation (Fig. 5), which indicates great influence of the environment at these stages. The pattern of heritabilities estimated in this study is similar to that observed by Sarmento et al (2008), studying Alpine dairy goats. The levels of heritability estimates reported in this study show that it is possible to increase milk yield in the second lactation of tropical dairy goats. Moreover, the different heritabilities estimated along the lactation curve suggest that it is possible to select for changing the patter of the lactation curves.

Conclusion

• The random regression model that used sixth order for the fixed and permanent environment effects, third order for the additive genetic effect, and four classes of residual variances was chosen as the most suitable one for genetic evaluation of tropical dairy goats during second lactation.
  
  
• Heritabilities obtained through the lactation curve were moderate to high (they ranged from 0.10 to 0.29), indicating that it is possible to promote genetic progress for milk yield in the second lactation of tropical dairy goats.
  
  
• Results reported in this paper suggest that it is possible to change the format of the lactation curve in the second lactation.

References

Brito L F, Silva F G, Oliveira HR, Souza N O, Caetano G C, Costa E V, Menezes G R O, Melo A L P, Rodrigues M T and Torres R A 2018 Modelling lactation curves of dairy goats by fitting random regression models using Legendre polynomials or B-splines. Candian Journal of Animal Science, 98, 73-83, from https://www.nrcresearchpress.com/doi/full/10.1139/cjas-2017-0019#.XY4tAS2ZOX1

CAPRAGENE 2019 Available at: http://srvgen.cnpc.embrapa.br/pmgcl/. Accessed on July 1, 2019.

FAO 2011 Food and Agriculture Organization of United Nations. FAOSTAT – FAO Statistics Division/ProdSTAT: livestock. Available at: http://www.fao.org/faostat/en/#home. Accessed on July 12, 2019.

Haenlein G F W 2004 Goat milk in human nutrition. Small Ruminant Research, 51, 155-163, from https://www.sciencedirect.com/science/article/pii/S0921448803002724

Irano N, Bignardi A B, Rey F S B, Teixeira I A M A and Albuquerque L G 2012 Genetic parameters for milk yield in Saanen and Alpina breed goats. Revista Ciência Agronomica, 43, 376-381, from http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1806-66902012000200022

Kirkpatrick M, Lofsvold D and Bulmer M 1990 Analysis of the inheritance, selection and evolution of growth trajectories. Genetics, 124, 979-93, from https://pdfs.semanticscholar.org/6b29/c9f35766dd1a32825501820c719094cbebf6.pdf

Menéndez-Buxadera A, Molina A, Arrebola F, Gil M and Serradilla J 2010 Random regression analysis of milk yield and milk composition in the first and second lactations of Murciano-Granadina goats. Journal of Dairy Science, 93, 2718-2726, from https://www.sciencedirect.com/science/article/pii/S0022030210002766

Menezes G R O, Torres R A, Sarmento J L R, Rodrigues M T, Brito L F, Lopes P S and Silva F G 2011 Random regression models in the milk yield evaluation in Saanen goats. Revista Brasileira de Zootecnia, 40, 1526-1532, from http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1516-35982011000700018

Meyer K 2005 Advances in methodology for random regression analyses. Australian Journal of Experimental Agriculure, 45, 847-858.

Meyer K 2007 WOMBAT- A tool for mixed model analyses in quantitative genetics by restricted maximum likehood (REML), Journal Zhejiang University Science B, 8, 815-821, from https://www.ncbi.nlm.nih.gov/pubmed/17973343

Nunez-Antón V N and Zimmerman D L 2000 Modelling nonstationary longitudinal data. Biometrics, 56, 699-705, from https://www.ncbi.nlm.nih.gov/pubmed/10985205

Oliveira H R, Brito L F, Lourenco D A L, Silva F F, Jamrozik J, Schaeffer L R and Schenkel F F 2019 Invited review: Advances and applications of random regression models: From quantitative genetics to genomics. Journal of Dairy Science, 102, 1-20, from https://www.sciencedirect.com/science/article/pii/S0022030219305521

Oliveira H R, Silva F F, Siqueira O H G B D, Souza N O, Junqueira V S, Resende M D V, Borquis R R A and Rodrigues M T 2016 Combining different functions to describe milk, fat, and protein yield in goats using Bayesian multiple-trait random regression models. Journal of Animal Science, 94, 1865-1874, from https://www.ncbi.nlm.nih.gov/pubmed/27285684

Sarmento J R, Albuquerque L G and Torres R A 2008 Comparison of random regression models for the estimation of genetic parameters in dairy goats. Revista Brasileira de Zootecnia, 37, 1788-1796, from http://www.scielo.br/scielo.php?script=sci_arttext&pid=S1516-35982008001000011

Silva F G, Torres R A, Brito L F, Euclydes R F, Melo A L P, Souza N O, Ribeiro Junior J I and Rodrigues M T 2013 Random regression models using Legendre orthogonal polynomials to evaluate the milk production of Alpine goats. Genetics and Molecular Research, 12, 6502-6511, from https://www.ncbi.nlm.nih.gov/pubmed/24390996

Silva F G, Torres R A, Silva L P, Ventura H T, Silva F F, Carneiro A P S, Nascimento M and Rodrigues M T 2014 Genetic evaluation of milk yeld in Alpine goats for the first four lactations using random regression models. Genetics and Molecular Research, 13, 10943-10951, from https://www.ncbi.nlm.nih.gov/pubmed/25526215

Souza J D F, Magalhães K A, Filho Z F H and Souza J D F 2017 Evolução do rebanho caprino entre 2007 e 2016, Boletim do Centro de Inteligência e Mercado de Caprinos e Ovinos, Sobral, 1, 5-7, from https://www.infoteca.cnptia.embrapa.br/infoteca/handle/doc/1078494

Zeng S S, Escobar E N and Popham T 1997 Daily variations in somatic cell count, composition and production of Alpine goat milk. Small Ruminant Research, 26, 253-260, from https://naldc.nal.usda.gov/download/12199/PDF

Received 27 September 2019; Accepted 5 January 2020; Published 1 February 2020

Go to top