Random regression models using B-splines functions provide more accurate genomic breeding values for milk yield and lactation persistence in Murrah buffaloes.

There is a great worldwide demand for cheese made with buffalo milk, due to its flavour and nutritional properties. In this context, there is a need for increasing the efficiency of buffalo milk production (including lactation persistence), which can be achieved through genomic selection. The most used methods for the genetic evaluation of longitudinal data, such as milk-related traits, are based on random regression models (RRM). The choice of the best covariance functions and polynomial order for modelling the random effects is an important step to properly fit RRM. To our best knowledge, there are no studies evaluating the impact of the order and covariance function (Legendre polynomials-LEG and B-splines-BSP) used to fit RRM for genomic prediction of breeding values in dairy buffaloes. Therefore, the main objectives of this study were to estimate variance components and evaluate the performance of LEG and BSP functions of different orders on the predictive ability of genomic breeding values for the first three lactations of milk yield (MY1, MY2, and MY3) and lactation persistence (LP1, LP2, and LP3) of Brazilian Murrah. Twenty-two models for each lactation were contrasted based on goodness of fit, genetic parameter estimates, and predictive ability. Overall, the models of higher orders of LEG or BSP had a better performance based on the deviance information criterion (DIC). The daily heritability estimates ranged from 0.01 to 0.30 for MY1, 0.08 to 0.42 for MY2, and from 0.05 to 0.47 for MY3. For lactation persistence (LP), the heritability estimates ranged from 0.09 to 0.32 for LP1, from 0.15 to 0.33 for LP2, and from 0.06 to 0.32 for LP3. In general, the curves plotted for variance components and heritability estimates based on BSP models presented lower oscillation along the lactation trajectory. Similar predictive ability was observed among the models. Considering a balance between the complexity of the model, goodness of fit, and credibility of the results, RRM using quadratic B-splines functions based on four or five segments to model the systematic, additive genetic, and permanent environment curves provide better fit with no significant differences between genetic variances estimates, heritabilities, and predictive ability for the genomic evaluation of dairy buffaloes.