A Comparison of Genomically Enhanced Breeding Values Predicted by Different Single-Step Approaches

The analyzed set (Table 1) represented data from the Polish national genetic evaluation of Holstein cattle from December 2021 for stature, which represents a highly heritable trait (h²=0.54), and foot angle, representing a low heritability trait (h²=0.09). It comprised 1,098,611 cows with stature phenotypes and 1,098,766 cows with foot angle phenotypes born between 1992 and 2019 as well as 141,397 bulls (stature) and 117,482 bulls (foot angle) born between 1986 and 2017 with pseudo-phenotypes expressed by their deregressed proofs (DRP) from the multiple across country evaluation (MACE) carried out by Interbull (interbull.org). Genomic data in the form of 46,118 SNP markers were available for 134,960 individuals, including 70,134 cows with phenotypes born between 2009 and 2021, as well as 64,826 bulls born between 1985 and 2021, of which 26,471 represented young individuals without pseudo-phenotypes and 38,355 were bulls with MACE-DRP. In our study, the pedigree of animals with phenotype data was truncated after the fifth generation, resulting in 1,555,995 individuals and 33 genetic groups.

The following single-step single-trait models were considered for the prediction of GEBV.

The ssG-BLUP (Christensen and Lund, 2010): (1) y=Xb+WG a+e, {\bf{y}} = {\bf{Xb}} + {{\bf{W}}_{\bf{G}}}\;{\bf{a}} + {\bf{e}}, where: y is the vector of dependent variables represented by cows’ measured phenotypes and bulls’ pseudo-phenotypes expressed by their MACE DRPs, b represents a vector of fixed effects scored at cows’ first parity including age at calving, lactation phase, and herd, while for bulls the corresponding fixed effects were represented by assigning common classes, a represents a vector of breeding values, and e is the vector of residuals. The underlying covariance structure of the random effects is given by a∼N(0,H_G σ_a²) and e∼N(0,Rσ_e²). H_G is given by (Liu et al., 2016): A11A12A21A22 ++ A12A22−1G−A22A22−1A21A12A22−1G−A22G−A22A22−1A21G−A22, \matrix{ {\left[ {\matrix{ {{{\bf{A}}_{11}}} & {{{\bf{A}}_{12}}} \cr {{{\bf{A}}_{21}}} & {{{\bf{A}}_{22}}} \cr } } \right]\; + } \cr { + \;\left[ {\matrix{ {{{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}\left( {{\bf{G}} - {{\bf{A}}_{22}}} \right){\bf{A}}_{22}^{ - 1}{{\bf{A}}_{21}}} & {{{\bf{A}}_{12}}{\bf{A}}_{22}^{ - 1}\left( {{\bf{G}} - {{\bf{A}}_{22}}} \right)} \cr {\left( {{\bf{G}} - {{\bf{A}}_{22}}} \right){\bf{A}}_{22}^{ - 1}{{\bf{A}}_{21}}} & {{\bf{G}} - {{\bf{A}}_{22}}} \cr } } \right],} \cr } where: A₁₁, A₁₂/A₂₁, and A₂₂ are the components of the numerator relationship matrix constructed based on the pedigree corresponding to non-genotyped animals, the covariance between non-genotyped and genotyped animals, as well as between genotyped animals, respectively, while G represents the genomic relationship matrix between genotyped animals. R is a diagonal matrix containing 1.00 for cows with phenotypes or n_i for bulls with MACE DRPs, with n_i representing a difference in effective daughter contributions of i-th bull between the MACE and the national evaluation. X and W_G denote the corresponding design matrices.

The ssSNP-BLUP (Gao et al., 2018): (2) y=Xb+Wsa+e, {\bf{y}} = {\bf{Xb}} + {{\bf{W}}_{\bf{s}}}{\bf{a}} + {\bf{e}}, where y, b, a, and e are the same as defined above, but a is parameterized as Zg + u with g being the vector of random SNP effects and being a vector of random additive residual polygenic effects. The covariance structure imposed on the residual effect (e) is the same as defined above, while for the additive genetic effect a the distribution is defined as a∼N(0,H_S σ_a²) where the structure of H_s is expressed by TGT′+DTGTZBGT′GZBBZ′T′BZ′B \left[ {\matrix{ {{\bf{TG}}{{\bf{T}}^\prime} + {\bf{D}}} & {{\bf{TG}}} & {{\bf{TZB}}} \cr {{\bf{G}}{{\bf{T}}^\prime}} & {\bf{G}} & {{\bf{ZB}}} \cr {{\bf{B}}{{\bf{Z}}^\prime}{{\bf{T}}^\prime}} & {{\bf{B}}{{\bf{Z}}^\prime}} & {\bf{B}} \cr } } \right] with D = (A¹¹)⁻¹, T = A₁₂A₂₂⁻¹ and B representing a diagonal matrix of the form I1∑i=1N2pi1−pi {\bf{I}}{1 \over {\sum\nolimits_{i = 1}^N {2{p_i}\left( {1 - {p_i}} \right)} }} with p_i denoting the allele frequency of the i-th SNP, N is the number of SNPs, and A¹¹ is the upper left block corresponding to A11A12A21A22−1 {\left[ {\matrix{ {{{\bf{A}}_{11}}} & {{{\bf{A}}_{12}}} \cr {{{\bf{A}}_{21}}} & {{{\bf{A}}_{22}}} \cr } } \right]^{ - 1}} . X and W_s denote the corresponding design matrices (Gao et al., 2018). Note that the variance components and the proportion of residual additive polygenic variance were not estimated, instead their values corresponding to the parameters used in the Polish national genetic and genomic evaluation (Table 2).

Solutions for the effects fitted in the above models were obtained using the MiXBLUP 3.0 software (Liu et al., 2014; Vandenplass, 2022) that solves the following equation: D⁻¹ M⁻¹ Cx = D⁻¹ M⁻¹ p where C represents the coefficient matrix corresponding to the Mixed Model Equations (MME) for solving (1) or (2), x is the vector of model parameters, and p is the RHS of MME, while M and D respectively represent the first-level and the second-level preconditioning matrices. The computations were performed on a dedicated computing server running the Linux Red Hat operating system with 260GB of RAM, 16 Intel Xeon CPUs with 2.20GHz, and 600GB of hard disk space. For a large national dairy population, such as the one used in our study, the numerically severe challenge is to obtain the inverse of the H_G and H_S matrix, in particular of its component related to genotyped individuals – the dense sub-matrix G. In our study, the following approaches were considered.

The GT approach (Lourenco et al., 2020; Mäntysaari et al., 2020), in which the G⁻¹ matrix is represented by 1wA22−1−1wT′T {1 \over w}{\bf{A}}_{22}^{ - 1} - {1 \over w}{{\bf{T}}^\prime}{\bf{T}} , with T = L^–1 Z′ A₂₂⁻¹ where w denotes the proportion of a residual polygenic variance and L is defined by Z′ A₂₂⁻¹ Z + wI = LL′.

The APY approach (Legarra et al., 2014) that divides the genotyped individuals into a core and non-core sub-groups for which the inverse is handled differently in such a way that the exact inverse is computed only for the core animal sub-group (G_c) and the covariance between core and non-core individuals, while the part of the matrix corresponding to the non-core genotyped animals (G_n) is handled as a diagonal matrix: G−1≈Gc−1000+−Gc−1GcnIMn−1−GncGc−1I {{\bf{G}}^{ - 1}} \approx \left[ {\matrix{ {{\bf{G}}_c^{ - 1}} & {\bf{0}} \cr {\bf{0}} & {\bf{0}} \cr } } \right] + \left[ {\matrix{ { - {\bf{G}}_c^{ - 1}{{\bf{G}}_{cn}}} \cr {\bf{I}} \cr } } \right]{\bf{M}}_n^{ - 1}\left[ {\matrix{ { - {{\bf{G}}_{nc}}{\bf{G}}_c^{ - 1}} & {\bf{I}} \cr } } \right] where the subscripts n and c represent non-core and core individuals respectively, and M_n is a diagonal matrix. In our study, four approaches to the selection of core animals were considered: APY3000top – where 3,000 genotyped bulls with the highest effective daughter contributions (EDC) were selected as core individuals, APY3000random – where the 3,000 core individuals were selected randomly from the genotyped population and the corresponding versions termed APY15000top, APY10000random, and APY15000random implementing 10,000 and 15,000 core individuals respectively. APY10000random proposed by Misztal et al. (2014) formed a basis to assess differences between scenarios implementing a smaller (APY3000) as well as a larger (APY15000) selection of core animals.

The ssSNP-BLUP approach (Gao et al., 2018) does not meet the numerical burden of the models based on G, since in terms of the modelling of genomic information only the inverse of the diagonal matrix B is required: HS−1=A11A120A21A22+1w−1A22−1−1wA22−1Z0−1wZ′A22−111−wB−1+1wZ′A22−1Z, {\bf{H}}_S^{ - 1} = \left[ {\matrix{ {{{\bf{A}}^{11}}} & {{{\bf{A}}^{12}}} & 0 \cr {{{\bf{A}}^{21}}} & {{{\bf{A}}^{22}} + \left( {{1 \over w} - 1} \right){\bf{A}}_{22}^{ - 1}} & { - {1 \over w}{\bf{A}}_{22}^{ - 1}{\bf{Z}}} \cr 0 & { - {1 \over w}{{\bf{Z}}^\prime}{\bf{A}}_{22}^{ - 1}} & {{1 \over {1 - w}}{{\bf{B}}^{ - 1}} + {1 \over w}{{\bf{Z}}^\prime}{\bf{A}}_{22}^{ - 1}{\bf{Z}}} \cr } } \right], where A11A12A21A22 \left[ {\matrix{ {{{\bf{A}}^{11}}} & {{{\bf{A}}^{12}}} \cr {{{\bf{A}}^{21}}} & {{{\bf{A}}^{22}}} \cr } } \right] represent the blocks corresponding to A11A12A21A22−1 {\left[ {\matrix{ {{{\bf{A}}_{11}}} & {{{\bf{A}}_{12}}} \cr {{{\bf{A}}_{21}}} & {{{\bf{A}}_{22}}} \cr } } \right]^{ - 1}} .

For the validation of GEBV prediction, bulls born after 2013 (7,296 bulls for stature; 6,200 bulls for foot angle) and cows born after 2016 (96,772 cows for stature; 96,771 cows for foot angle) were removed from the phenotype vector (y) of equations (1) and (2). Based on such truncated data sets, the GEBVs of bulls with EDC greater than 20 were predicted by models (1) and (2). Prediction accuracy was assessed by a weighted linear regression (Mäntysaari et al., 2010): GEBV_f = b₀ + b₁GEBV_t + e, with GEBV_f representing the vector of GEBVs predicted based on the full data set with all available individuals while GEBV_t contains the GEBVs predicted based on the truncated data set. For i-th bull, weights were defined as EDCiEDCi+k {{{EDC_i}} \over {{EDC_i} + k}} with k=4−h2h2 k = {{4 - {h^2}} \over {{h^2}}} . The unbiased evaluation is then represented by b₁ equal to 1.00.

The above linear regression equation was fitted using the lm function and Pearson correlation coefficients, used to assess differences between GEBVs across scenarios were calculated using R software (www.rstudio.com).

The validation data set for stature and foot angle comprised 1,727 and 1,725 bulls with EDC > 20, respectively. For stature, generally, no marked differences in the estimated slope of the linear regression were observed between the models, varying between 0.94 (APY3000top) and 1.02 (APY10000random). Furthermore, except for APY3000, the differences between R² corresponding to each model were small. For stature, the GT approach achieved the highest R² of 0.83, while the lowest R² corresponded to 3,000 core animal scenarios with 0.57 for APY3000top and 0.60 for APY3000random. Similar results were observed for foot angle, albeit with overall lower R². GT and SNP-BLUP achieved the highest value of 0.76, while the lowest of 0.60 was attributed to APY3000random and 0.62 for APY3000top. Details of the validation results are summarized in Table 3. Pearson correlations (Figure 1) between GEBVs predicted for the validation bulls from the full and the truncated datasets were calculated for all pairs of models. For stature, they demonstrated a good agreement only between SNP-BLUP and GT, as expressed by correlations of 0.996. Interestingly, the lowest correlations of 0.810 were observed between both models implementing APY3000 (i.e. APY3000random and APY3000top), as well as between APY3000top and APY15000top. For foot angle correlations are generally higher, and the pattern is very similar, with the correlation between SNP-BLUP and GT reaching 0.999 and the lowest correlation of 0.870 observed between both APY3000 scenarios, between GT and APY3000top, as well as between GP and APY3000random. Notably, for both traits, the correlations within APY3000top and within APY3000random are markedly lower than for the other models

Figure 1.

Pearson correlation coefficients of GEBVs predicted by different model variants

Table 3.

Model variants	b^0 {{\hat b}_0}	SEb^0 {\rm{SE}}( {{{\hat b}_0}})	b^1 {{\hat b}_1}	SEb^1 {\rm{SE}}( {{{\hat b}_1}})	R2
Stature, 1,727 validation bulls, h2 = 0.54
SNP-BLUP	−3.90	0.32	1.01	0.01	0.77
GT	−4.40	0.27	1.01	0.01	0.83
APY3000top	−1.81	0.44	0.94	0.02	0.57
APY3000random	−2.20	0.44	0.99	0.02	0.60
APY10000random	−3.74	0.36	1.02	0.02	0.72
APY15000top	−2.59	0.32	0.96	0.01	0.75
APY15000random	−2.32	0.33	0.96	0.01	0.73
Foot angle, 1,725 validation bulls, h2 = 0.09
SNP-BLUP	−2.03	0.18	1.03	0.01	0.76
GT	−1.96	0.18	1.04	0.01	0.76
APY3000top	−0.68	0.21	0.97	0.02	0.62
APY3000random	−0.52	0.26	1.03	0.02	0.60
APY10000random	−2.04	0.21	1.02	0.02	0.69
APY15000top	−2.15	0.18	1.02	0.01	0.75
APY15000random	−2.10	0.19	1.03	0.02	0.73

Figure 2 shows the differences between the GEBVs predicted by ssSNP-BLUP compared to the GEBVs predicted by the six ssG-BLUP implementations. For the difference in GEBV between SNP-BLUP and GT, all genotyped animals were considered, while for all comparisons involving APY, only the core individuals were used. The original solutions were rescaled for each model by subtracting the mean of a cow base population to provide GEBVs comparable across all models. For stature, the differences were generally small across bull birth years (2A), however, for foot angle a different pattern emerged (2B) as for all comparisons, except GT and the most informative APY15000top, the differences were unstable across bull birth years. Furthermore, we compared correlations of GEBVs predicted by the full and the truncated models for individuals representing base cows (i.e. cows with phenotype records) and bulls (i.e. bulls, which have a minimum of one daughter). The dashed line divides the plot into animals defined as old and young in the validation process. For stature, regardless of the model, correlations between full and truncated data sets, for old animals, were very close to unity. For young candidates, we observed a declining trend with the lowest correlations calculated for APY3000random (Figure 3 A). For foot angle, similar results were obtained, except for APY3000random for females, where we observed a decreasing correlation trend as early as from 2004. From 2008, correlations dropped below 0.9 (Figure 3 B). Finally, we compared the overlap of bulls with the top 50 GEBV predictions resulting from different models (Figure 4). The number of bulls common across all models for stature (Figure 4 A) was 39, while the highest number of exclusive bulls (3) was observed in the top 50 list resulting from the APY15000top. Regarding foot angle (Figure 4 B), the number of common bulls was even lower, being 31, while 8 were exclusive for the APY3000random.

Figure 2.

Mean GEBV difference divided by the genetic standard deviation between SNP-BLUP and GT, as well as SNP-BLUP and APY models implementing different scenarios of core individual selection. (A) stature, (B) foot angle

Figure 3.

Correlation between full and truncated data set for cows with phenotype records and bulls with a minimum of one daughter. (A) stature, (B) foot angle

Figure 4.

An upset plot of bulls with the top 50 GEBVs. (A) stature, (B) foot angle

The wall clock time corresponding to setting up and solving models (1) or (2) using the MiXBLUP software and the peak memory consumption varied considerably between solvers. The exact values are specified in Table 4, but generalizing for both traits, SNP-BLUP and APY3000 were the fastest, closely followed by APY10000 random. The wall clock time of GT was the longest, twice the time of APY15000top. The peak memory consumption by SNP-BLUP was on the order of ten times lower than for the remaining solvers. In the case of iterations, SNP-BLUP required the most iterations (673 for stature and 1,027 for foot angle) and an average of 2.3 seconds per iteration. The lowest number of 335 iterations and an average of 0.18 seconds per iteration were used for stature by APY3000top and 499 and 0.18 seconds for foot angle by APY3000random.

Regarding the prediction of GEBV on the active population scope, no marked differences between solvers (except the APY with only 3,000 core individuals) were observed, still GT and SNP-BLUP solvers resulted in the highest prediction accuracy, while the effectiveness of the APY model depends on the size of the set of core animals: small cores lead to a noticeable decrease in accuracy, while cores with large numbers minimize these differences. Another factor influencing the choice of the solver is its computational efficiency expressed by the computation time and memory resources required, which was also indicated by Koivula et al. (2012). However, it should be kept in mind that the ranking of the top bulls is not identical between models, which has implications for the breeding industry in terms of semen pricing and selection since, typically, the top-ranking bulls are the most widely used.

Considering the abovementioned aspects, the results of our study indicate that GT and SNP-BLUP solvers offer more modelling and computational advantages.