Milk yield is one of the main criteria for selecting dairy buffaloes. The recent application of genomic tools for the genetic evaluation of this trait has enabled genetic improvement of the species. In this scenario, the use of single-step genomic best linear unbiased prediction (ssGBLUP) methods in studies of multiple traits of buffaloes has also been found to give a substantial increase in prediction accuracy compared to the results of traditional methods based only on pedigree information (Aspilcueta-Borquis et al., Reference Aspilcueta-Borquis, Araujo Neto, Santos, Hurtado-Lugo, Silva and Tonhati2015; Cesarani et al., Reference Cesarani, Biffani, Garcia, Lourenco, Bertolini, Neglia, Misztal and Macciotta2021). This approach is used globally because it combines phenotypic data from genotyped and non-genotyped animals simultaneously in one step (Aguilar et al., Reference Aguilar, Misztal, Johnson, Legarra, Tsuruta and Lawlor2010).
As an alternative option to the single-step method applied in GBLUP, Fernando et al. (Reference Fernando, Dekkers and Garrick2014) proposed an adaptation of Bayesian regression models to include data from non-genotyped animals in the analysis as well, which might enable new perspectives for genomic prediction studies. Therefore, since there are few published studies exploring the single-step Bayesian regression (ssBR: Zhou et al., Reference Watanabe, Nascimento, Bernardes, Costa, Lôbo, Baldi and Munari2018; Santos et al., Reference Santos, Araujo Neto, Fernandez, Santos, Cunha, Aspilcueta-Borquis and Tonhati2024), we assessed the predictive ability of ssBR methods applied to milk yield of Murrah buffaloes.
Material and methods
This study was conducted using previously collected data, so no approval from animal-use ethics committees was required. We used test-day milk records of the first lactation of Murrah buffaloes born between 1995 and 2017, from farms located in the Brazilian states of São Paulo, Ceará and Rio Grande do Norte. Using these test day records, we calculated the 270-day cumulative milk yield (CMY270) and then formed contemporary group of cows that had calved on the same farm, year and season (dry or rainy). Animals with CMY270 records outside the range between ± 3 standard deviations from their contemporary group (CG) averages and animals belonging to CGs with fewer than five individuals were removed. The structure of the database used in the analyses after the consistency stage is described in Table 1.
Summary of the database used for the ability of prediction analysis performing several methods of single-step Bayesian regression models for milk yield of Murrah buffaloes cows
From a pedigree file containing 3,320 buffaloes, 553 animals were genotyped with the Axiom Buffalo® Genotyping Array 90 K. Only SNPs present on autosomal chromosomes according to the UOA_WB assembly, as well as those with call rate > 95%, MAF > 3% and significance level by the Hardy-Weinberg equilibrium test (p-value < 10−5) remained in the analysis. All samples had a call rate > 90%. After the consistency stage, 553 animals and 44,574 SNP markers were retained. The estimation of variance components and prediction of genomic values were carried out using the single-step Bayesian regression (ssBR) method proposed byFernando et al. (Reference Fernando, Dekkers and Garrick2014). The ssBR models combine all available data from genotyped and non-genotyped animals and use imputed marker covariates for animals that were not genotyped (Fernando, Dekkers, Garrick, Reference Fernando, Dekkers and Garrick2014). The general model can be described as:
\begin{equation*}\left[ {\begin{array}{*{20}{c}}
{{y_n}}\\
{{y_g}}
\end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
{{X_n}}\\
{{X_g}}
\end{array}} \right]\beta + \left[ {\begin{array}{*{20}{c}}
{{Z_n}}&0\\
0&{{Z_g}}
\end{array}} \right]\left[ {\begin{array}{*{20}{c}}
{{M_n}\alpha + \epsilon }\\
{{M_g}\alpha {\textrm{}}}
\end{array}} \right] + e\end{equation*} where vectors and matrices for non-genotyped animals are denoted with the subscript 
$n$ and those of genotyped animals have the subscript 
$g$. Thus, 
${y_n}$ and 
${y_g}$ represent the vectors of phenotypic observations; 
$\beta $ is the vector of systematic effects of GC and age at calving; 
${X_n}$ and 
${X_g}$ are design matrices for the fixed effects; 
${Z_n}$ and 
${Z_g}$ denote incidence matrices associated with the genomic values of the animals; 
$\alpha $, 
$\beta $ and 
$e$, represent, in this order, the vectors of marker effects, imputation residual and model residual; 
${M_g}$ is the marker matrix for the genotyped animals and 
${M_n}$ is the imputed marker matrix for the non-genotyped animals. The matrix 
${M_n}$ can be obtained by computing:
\begin{equation*}{M_n} = {A_{ng}}A_{22}^{ - 1}{M_g}\end{equation*} where 
$\textit{A}_{ng}$ represents the pedigree-based relationship matrix between genotyped and non-genotyped animals and; 
$\textit{A}_{22}^{-1}$ represents the inverse of the pedigree-based relationship matrix between genotyped animals. The following ssBR methods were used: BayesA (ssBA), BayesB (ssBB), BayesC (ssBC), Bayesian Lasso (ssBL) and Bayesian ridge regression (ssBRR). These are described in online Supplementary Table S1. For models based on mixture distributions (ssBB and ssBC), the proportions of markers having effect (π) were assumed both as fixed, with values of 99% and 90%, and as unknown, where two approaches to estimate π were applied (ssBayesBπ and ssBayesCπ). The “hibayes” package (Yin, Zhang, Li, Zhao, Liu, Reference Zhou, Mrode, Zhang, Zhang, Li and Liu2022), available in the R program, was used to perform the ssBR analyses. For the variance component estimation, a total of 40,000 samples and a burn-in interval of 5,000 samples were generated. The convergence was checked using graphical analysis and the Geweke (1992) and Heidelberger (Reference Heidelberger and Welch1983) tests, available in the Coda package (Plummer et al., Reference Plummer, Best, Cowles and Vines2006).
The predictive abilities of the methods were assessed by analyzing the full dataset and a partial dataset, obtained by omitting phenotypic records of genotyped animals born after 2010 (validation group), in a total of 80 animals. The predicted genetically estimated breeding values (GEBVs) from the genotyped animals from the validation group were used to calculate the following statistics:
\begin{equation*}\begin{array}{*{20}{c}}
{\textrm{Accuracy}}& {{\textrm{Dispersion}}} \\
{{r_{A\hat{A}}} = \frac{{cor\left( {{y_{adj}},{\hat{{u}}_p}} \right)}}{{\sqrt {{h^2}} }}}& {{b_{wp}} = \frac{{cov\left( {{\hat{{u}}_w},{\hat{{u}}_p}} \right)}}{{var\left( {{\hat{{u}}_p}} \right)}}}
\end{array}\end{equation*} where û w and û p denote the GEBVs of the animals from the validation group obtained by analyzing the full and partial datasets, respectively; and 
$\textit{Y}_{adj}$ is the CMY270 phenotype adjusted for the fixed effects of the models. Model details are given in the online Supplementary File.
Results
The heritability estimates obtained in this study, represented by the average of the posterior distributions and their 95% high-density interval (minimum - maximum), ranged from 0.207 (0.148 - 0.259) to 0.284 (0.246 - 0.318) for ssBayesC99 and ssBayesA, respectively. The accuracy (Table 2) found ranged from 0.55 (ssBayesRR) to 0.584 (ssBayesCπ). The ssBayesRR method presented the lowest GEBV dispersion (Table 2), followed by the ssBayesL method. These two were the only models with values greater than 0.90.
Estimates of accuracy and dispersion of the predicted GEBVs by different methods of single-step Bayesian regression for milk yield in Murrah buffaloes
Numbers in bold highlight the best values for each criterion.
Discussion
The heritability estimates had moderate magnitude, and were similar to those found in other studies with the species. Kumar et al. (Reference Kumar, Vohra, Rawan and Lathwal2023) described heritability estimates of 0.178 for 305-day cumulative milk yield of Murrah buffaloes. Cesarani et al. (Reference Cesarani, Biffani, Garcia, Lourenco, Bertolini, Neglia, Misztal and Macciotta2021) report values of 0.25 for Mediterranean buffaloes. In agreement, Aspilcueta-Borquis et al. (Reference Aspilcueta-Borquis, Araujo Neto, Santos, Hurtado-Lugo, Silva and Tonhati2015), performing ssGBLUP, estimated 0.25 for heritability of 305-day cumulative milk yield of Murrah buffaloes from Brazil. The similarity between the heritability estimates obtained in this study and those reported in the literature indicates that the ssBR method does not greatly affect the estimation of this parameter. This is possibly due to the polygenic structure of the trait, which also corroborates the observations of other studies, with small differences found when using traditional (pedigree) and genomic information for the estimation of genetic parameters of cumulative milk yield (Aspilcueta-Borquis et al., Reference Aspilcueta-Borquis, Araujo Neto, Santos, Hurtado-Lugo, Silva and Tonhati2015; Watanabe et al., Reference Yin, Zhang, Li, Zhao and Liu2021).
The values of accuracy (Table 2) found in our study ranged from 0.55 (ssBayesRR) to 0.584 (ssBayesCπ), magnitudes lower than those reported by other authors who have used ssGBLUP to predict milk production. Cesarani et al. (Reference Cesarani, Biffani, Garcia, Lourenco, Bertolini, Neglia, Misztal and Macciotta2021) described accuracy values for milk production of 0.60 for BLUP and values between 0.60 and 0.77 for ssGBLUP using different validation strategies. Aspilcueta-Borquis et al. (Reference Aspilcueta-Borquis, Araujo Neto, Santos, Hurtado-Lugo, Silva and Tonhati2015) reported values of 0.66 and 0.72 for milk yield by BLUP and ssGBLUP, respectively, when analyzing accuracy based on the average of the prediction error variance (PEV). With regard to the methods based on mixture distributions (ssBayeB and ssBayesC), a decrease in the proportion of markers with no effect led to an increase in accuracy (π = 90% to π = 99% - Table 2). Possibly, this trend is related to the quantitative nature of the trait in question. In this study, the estimates for the parameter π in models assuming an unknown value were 0.958 and 0.919 for BayesB and BayesC.
When π was considered as a parameter to be estimated, there was an increase in the accuracy of the ssBayesCπ method, opposite to the trend observed with the ssBayesBπ method. According to Lee et al. (Reference Lee, Cheng, Garrick, Golden, Dekkers, Park, Lee and Fernando2017), models that estimate π can be beneficial, but the quality of the estimation of this parameter requires large datasets. Determining an appropriate π may be less critical for ssBR methods, since these models allow the use of a larger volume of data. However, our results suggest that the ssBayesBπ method might have been more sensitive to our data structure than the ssBayesCπ method.
The ssBayesRR method presented the lowest GEBV dispersion, followed by the ssBayesL method. These were the only models with values greater than 0.90. Problems with GEBVs dispersion in prediction studies with milk yield were verified in both Holstein (Lee et al., Reference Lee, Dang, Choy, Do, Cho, Kim, Kim and Lee2019) and buffalo cows (Cesarani et al., Reference Cesarani, Biffani, Garcia, Lourenco, Bertolini, Neglia, Misztal and Macciotta2021; Herrera et al., Reference Herrera, Flores, Duijvesteijn, Moghaddar and van der Werf2021). Lower dispersion values indicate complications when selecting young animals (which do not yet have data or recorded daughters in the case of young bulls), since their GEBVs tend to be overestimated, causing distortions in the classification rank and reducing the accuracy of the final genetic gain.
In conclusion, the results of the ssBR performance indicate that the Bayesian Lasso method is the best choice for the genomic evaluation of milk yield in buffaloes. This model predicted GEBVs with one of the highest accuracy values, and at the same time with low dispersion compared to most of the assessed models.
Supplementary material
The supplementary material for this article can be found at https://doi.org/10.1017/S0022029925000317.
Acknowledgements
We thank the São Paulo State Research Foundation (FAPESP) for the resources provided for the development of this study.
