Genomic evaluation of carcass traits of Korean beef cattle Hanwoo using a single-step marker effect model

Abstract Hanwoo beef cattle are well known for the flavor and tenderness of their meat. Genetic improvement programs have been extremely successful over the last 40 yr. Recently, genomic selection was initiated in Hanwoo to enhance genetic progress. Routine genomic evaluation based on the single-step breeding value model was implemented in 2020 for all economically important traits. In this study, we tested a single-step marker effect model for the genomic evaluation of four carcass traits, namely, carcass weight (CW), eye muscle area, backfat thickness, and marbling score. In total, 8,023,666 animals with carcass records were jointly evaluated, including 29,965 genotyped animals. To assess the prediction stability of the single-step model, carcass data from the last 4 yr were removed in a forward validation study. The estimated genomic breeding values (GEBV) of the validation animals and other animals were compared between the truncated and full evaluations. A parallel conventional best linear unbiased prediction (BLUP) evaluation with either the full or the truncated dataset was also conducted for comparison with the single-step model. The estimates of the marker effect from the truncated evaluation were highly correlated with those from the full evaluation, ranging from 0.88 to 0.92. The regression coefficients of the estimates of the marker effect for the full and truncated evaluations were close to their expected value of 1, indicating unbiased estimates for all carcass traits. Estimates of the marker effect revealed three chromosomal regions (chromosomes 4, 6, and 14) harboring the major genes for CW in Hanwoo. For validation of cows or steers, the single-step model had a much higher R2 value for the linear regression model than the conventional BLUP model. Based on the regression intercept and slope of the validation, the single-step evaluation was neither inflated nor deflated. For genotyped animals, the estimated GEBV from the full and truncated evaluations were more correlated than the estimated breeding values from the two conventional BLUP evaluations. The single-step model provided a more accurate and stable evaluation over time.


Introduction
Hanwoo is an indigenous beef breed in the Republic of Korea and its meat is recognized for its tenderness, flavor, and taste (Nwogwugwu et al. 2020). Beginning in the 1980s, an extensive genetic improvement program was implemented in Hanwoo cattle, using routine performance recordings, and progeny tests. Rapid genetic progress has been achieved by applying the best linear unbiased prediction (BLUP) model for the routine genetic evaluation of economically important traits (Nwogwugwu et al. 2020). To hasten genetic progress, a routine genomic evaluation system was set up for Hanwoo in 2020 by the Korean Animal Improvement Association (KAIA).
Since the invention (Meuwissen et al. 2001) and introduction (VanRaden 2008) of genomic selection and evaluation, research projects on the routine genomic evaluation of Hanwoo  have been based on the single-step genomic BLUP model (ssGBLUP, Aguilar et al. 2010;Christensen and Lund 2010). Because the single-step model considers all phenotypes, genotypes, and pedigrees simultaneously, it has become the model of choice for routine genetic and genomic evaluation in diverse livestock species (Misztal et al. 2020). In this study, we chose an alternative to the single-step model that allows direct estimation of the effects of single nucleotide polymorphism (SNP): the single-step SNP BLUP model (ssSNPBLUP, Liu et al. 2014). Vandenplas et al. (2019) demonstrated the computational efficiency of the ssSNPBLUP model over other single-step models for large-scale genomic evaluations. Alkhoder et al. (2022) showed the advantages of ssSNPBLUP over the multistep SNP BLUP model for the genomic evaluation of conformation traits in the German Holstein population.
Over the past decade, numerous studies have been conducted on the genomic selection of Hanwoo (Lee et al. (2017) Mehrban et al. (2019)). These studies primarily used genotype data of bulls with progeny and a limited number of steers with their own performance data. To the best of our knowledge, no study has been published on the genomic evaluation of Hanwoo based on a large number of genotyped cows, steers, and young animals. Recently, Mehrban et al. (2021) demonstrated the enhanced accuracy of single-step genomic prediction using a multi-trait model compared to that using a single-trait model.
The objectives of this investigation were as follows: 1) To apply the ssSNPBLUP model for the multi-trait genomic evaluation of carcass traits in Hanwoo using genotype data from cows, steers, and young animals., 2) To conduct a genomic validation by simulating a previous genomic evaluation conducted 4 yr ago., and 3) To identify chromosomal regions important for carcass traits in Hanwoo.

Data and materials
Our genomic evaluation of the Hanwoo animals followed the publicly available animal care and use standards in Korea and Germany, equivalent to the institutional animal care and use committee (IACUC) in the United States. Four carcass traits were routinely collected from steers and cows of Hanwoo cattle: carcass weight (CW), eye muscle area (EMA), backfat thickness (BF), and marbling score (MS). These carcass traits are among the economically important traits in Hanwoo (Kim et al. 2021;Bedhane et al. 2019;Mehrban et al. 2019;Nwogwugwu et al. 2020;Son et al. 2020;Seo et al. 2022). An updated scheme for data recording were introduced for Hanwoo in 2008. In recent years, the Korean breeding organization KAIA has begun genotyping male and female animals for genomic selection. Table 1 shows the number of steers and cows with recorded data on all four carcass traits, and the number of genotyped animals by year of birth. In total, carcass trait data were available for 8,023,666 cows and steers.
Over the last 10 yr, the number of genotyped female animals has increased continuously. The total number of genotyped animals was 29,965 using three SNP chips: Affymetrix Ancestors of the genotyped or phenotyped animals presented in Table 1 were traced back to the pedigree as far as possible. The total number of animals included in the final pedigree of the single-step evaluation was 13,534,295. Thirty-three phantom parent groups were defined for the missing parents of animals based on the four selection paths and birth years of the animals.

Multi-trait single-step SNP BLUP model
For an animal with phenotypic records of carcass traits, the following multi-trait model was used, as described by Son et al. (2020): where y is a vector of observations of the four carcass traits of the animal, f is a vector of fixed effects of herd-year and year-season on slaughter and sex of the animal, b is a vector of linear regression on the age of the animal at slaughter, u is a vector of breeding values of the four traits of the animal, and e is a vector of random error effects associated with the observations. X f is the incidence matrix for the classification fixed effects f, X b is the design matrix for the regression effect b.
If an animal is genotyped, the genomic breeding value (GEBV) of a carcass trait i can be decomposed into two components according to the ssSNPBLUP model (Liu et al. 2014). Therefore, model 1 can also be expressed for genotyped animals as follows: where g is a vector of all SNP marker effects, a is a vector of the residual polygenic effects (RPG; Liu et al. 2011) of the genotyped animals, and Z is the design matrix containing SNP genotype data of the animals. It was assumed that RPG effects explained 20% of the additive genetic variance, k = 0.2. This RPG parameter suggests that the SNP marker genotype information accounted for 80% of the additive genetic variance. Based on our validation study, the parameter k = 0.2 seemed to be an optimal value for the evaluated traits of Hanwoo. Let σ 2 ui represent the total additive genetic variance of the i-th carcass trait; then, the average genetic variance of SNP markers was calculated using the ssSNPBLUP model as follows: where m = 47,200 is the number of SNP markers, and p j is the allele frequency of the jth SNP marker.
The variance and covariance components for additive genetic and random error effects were obtained from a recent study (Son et al. 2020) on carcass traits of Hanwoo for single-step evaluation. The heritability values ranged from 0.28 to 0.48. Moderate positive genetic correlations were found between all trait pairs except for low negative correlations between BF and EMA or MS.
The software package MiX99 (Strandén and Lidauer, 1999) was used for the single-step evaluation. The ssSNP-BLUP model was implemented in MiX99 in a special manner with SNP markers treated as if they were animals with neither parents nor progeny known (Alkhoder et al. 2022).

Genomic validation of the single-step SNP BLUP model
To assess the stability of genomic prediction in Hanwoo, genomic validation was performed by simulating a previous genomic evaluation conducted 4 yr ago. All phenotypic records of carcass traits were removed for cows or steers born after 2016 (Table 1). Genotyped cows and steers with their own carcass traits, born in 2017 and thereafter, were chosen as validation animals. There were 480 and 419 validation cows and steers, respectively. Owing to the lack of genotype data, the total number of genotyped bulls with progeny was 355 in the full dataset. The number of genotyped bulls, defined as validation bulls, was 83. Therefore, we decided not to use bulls as validation animals for genomic validation.
To compare the predictive ability of the single-step model 1, we conducted a conventional BLUP evaluation using the same animals as those used for the single-step evaluation, with young animals included in the full and truncated datasets.

Results and Discussion
The pedigree error rate in Hanwoo was investigated using a genotyped subpopulation. Among the 29,965 genotyped animals, 25,501 and 6,604 were sire-and dam-genotyped, respectively. After parentage verification and discovery in the genotyped subpopulation, we obtained a pedigree error rate of 6.9% on the sire side and 6.1% on the dam side. Pedigree error rates in the genotyped population of Hanwoo seem to be within the range reported for beef or dairy cattle in literature.
All single-step and conventional BLUP evaluations with the full and truncated datasets were run on a Linux computer with 128 GB random access memory and 12 cores. For the single-step evaluation using the full dataset, the total number of effects to be estimated was 57,260,624. The peak memory usage was approximately 30.8 Gb and the total clock time was approximately 1 h to reach a predefined convergence criterion after 549 iterations. With the same dataset, the ssGBLUP model required a similar clock time to reach the same predefined convergence criterion, but with a smaller peak memory usage of 10.3 Gb. The ssSNPBLUP model required more memory than the ssGBLUP model, because the number of SNP markers was greater than the number of genotyped animals. With more genotyped animals added later, the ssSNPBLUP model will be more efficient in terms of memory usage (Alkhoder et al. 2022) compared to the ssGBLUP model.
Four evaluations were conducted: single-step evaluations with the full and truncated datasets, and conventional BLUP evaluations with the full and truncated datasets. The solutions of the mixed model equations for additive genetic effects were adjusted for the four evaluations for a common base population of animals (born in 2015) with carcass records. For a direct comparison of the evaluation results across the four traits, the adjusted GEBV were divided by the respective genetic standard deviation of each trait.

Estimation of SNP effects
The single-step model ssSNPBLUP (Liu et al. 2014) directly estimated the effects of all the SNP markers, besides the GEBV of all the animals in the pedigree and all the fixed effects of model 1. We compared the estimates of SNP effects between the full and truncated evaluations of the single-step model. The stability of the GEBV over time is one of the most important quality control criteria for routine genomic evaluations. If the estimates of SNP effects are highly correlated between two evaluations at earlier and later times, the GEBV of the animals will also be highly correlated, particularly for young animals without phenotypic data. A high SNP-effect correlation between the two evaluations at different times is desirable for stable genomic prediction. Figure 1 shows the correlation coefficients of estimates of SNP effects between the full and truncated evaluations for each trait. The correlation coefficients of the SNP effects range from 0.88 to 0.92 for traits CW and BF, respectively. The regression coefficients of the estimates of SNP effects of the full and truncated evaluations are shown in Figure 2. The regression coefficients are close to their expected values of 1, varying from 0.97 for BF and 1.00 for EMA. The regression slopes indicate that the estimates of SNP effect were neither inflated nor underestimated. The relatively high correlations of estimates of the SNP effects in Figure 1 may be, to a certain degree, caused by the high proportion of common phenotypic data shared between the truncated and full evaluations. However, the high percentage of common phenotypic data between the full and truncated evaluations should have a negligible impact on the regression coefficients of the estimates of SNP effects, as shown in Figure 2.
In contrast to the usual single-marker genome-wide association study (GWAS), the genomic model [1] jointly estimated the effects of all 47,200 SNP markers. Figures  3 and 4 show the plots of SNP effects for the CW and MS traits, respectively. All the estimates of the SNP effects were divided by the genetic standard deviation of SNP markers, σ SNP , such that they can be compared directly across the four traits. In contrast to the regular GWAS, we did not include or exclude the SNP markers by explicitly calculating the probability of the t-test for the effect estimate of every SNP marker because the focus of our genomic prediction was not on the significance testing of the SNP markers and both significant and nonsignificant SNP were included. Aguilar et al. (2019) provided a statistical method for approximating the P values of the estimates of SNP effects from the ssGBLUP model. Further research is necessary to develop similar methods for approximating the P values or prediction error variances of the estimates of SNP effects achieved using our ssSNPBLUP model. As shown by Liu et al. (2011), the variance in the estimates of SNP effects increased with the size of the reference population. More SNP markers with larger effects would yield an even larger estimate of effect with an increasing number of reference animals. Our plots of SNP effects resembled those of regular Manhattan plots because the standard errors of the estimates of SNP effects of all the SNP markers did not differ significantly. Our plots of SNP effects provide sufficient information to detect chromosomal regions with large effects. Chromosomes 4, 6, and 14 appear to harbor large genes or mutations in the CW trait ( Figure 3). Compared to CW, the SNP effects of trait MS have a higher variation, with many SNP effects greater than 50% of σ SNP .
In their GWAS, Lee et al. (2013) also reported SNP markers with large effects on chromosome 14 in CW. Similar to our plot for CW, Kim et al. (2021) found peaks on chromosomes 4 and 6, in addition to the highest peak on chromosome 14. These two previous studies confirmed that our chromosomal regions harbored major genes or mutations in CW.
In contrast to CW, no dominant peak was observed in the plot for trait MS ( Figure 4); however, many SNP markers with medium-sized effects were identified for trait MS. Similar SNP effects were also revealed in studies by Bedhane et al. (2019) and Kim et al. (2021).
For trait BF (data not shown), we also identified a region on chromosome 19 as described by Kim et al. (2021). Our plots of SNP effects for MS resembled that of Kim et al. (2021). As shown in Figure 4, for the MS trait, no major loci appeared to exist for the EMA trait (data not shown). However, there were many SNP markers with medium-sized effects, with a standardized marker effect between 0.3 and 0.8. We did not see a strong indication of major genes for trait EMA in our study, which does not agree with the findings of Seo et al. (2022), who identified significant regions on chromosomes 3 and 14 via a single-marker GWAS. One  reason for the difference in significant chromosome regions for trait EMA may be that most genotyped animals were cows in our study versus bulls by Seo et al. (2022). Additional research with more genotyped animals is necessary to verify our results concerning the estimates of SNP effects for trait EMA in Hanwoo.

Genomic validation results
A linear regression test (LR test, Legarra and Reverter, 2018) was performed between the GEBV of the single-step full and truncated evaluations of validation cows or steers. For comparison with the single-step model, the LR test was applied to the conventional BLUP EBV with the full and truncated  datasets. Table 2 shows the results of genomic validation for the two groups of validation animals of Hanwoo for carcass traits.
It can be seen that the R 2 values of the linear regression model, being related to the accuracy of genomic prediction, are higher for the single-step than for the conventional BLUP model in either group of validation animals for each of the four carcass traits. The difference in R 2 values between the two models ranges from 0.13 for trait BF of validation steers to 0.23 for trait MS of validation cows. The single-step model had a markedly higher correlation of GEBV compared to the conventional BLUP model, suggesting a higher stability of GEBV in the single-step evaluation over time.
The regression intercept in Table 2, b 0 , is expressed in the unit of the genetic standard deviation of the corresponding trait. The lowest b 0 value, −0.03, is found for trait MS and CW of validation cows and steers, respectively, whereas the highest value was 0.10 for trait EMA of validation steers. A slight overestimation of the intercept, shown as a negative value, was observed in all combinations of traits and animal groups, except for validation steers in EMA and MS, which showed a slight underestimation with a positive b 0 value. Overall, the regression intercept values do not deviate significantly from the expected value of 0.
The regression slope, b 1 , indicates the inflation of genomic prediction if b 1 < 1, and an underestimation of genomic prediction if b 1 > 1. None of the b 1 values in Table  2 deviated markedly from the expected value of 1, suggesting that the single-step genomic prediction or conventional BLUP prediction was neither inflated nor deflated. Relative to the single-step evaluation, conventional BLUP evaluation showed a greater deviation in b 1 from the expected value of 1.

Genetic trends in steers
The average GEBV of all steers with their own carcass records, genotyped, and non-genotyped, was calculated by the year of birth for the single-step evaluations with full and truncated datasets. Figure 5 shows the genetic trends in steers of the single-step model for trait MS. Notable genetic progress has been made in trait MS of Hanwoo steers, amounting to 1.6 genetic standard deviations from birth year 2000 to 2020. The single-step evaluation with the truncated dataset exhibited a slightly lower trend than that of the full evaluation. The genetic trends of the conventional BLUP evaluation (data not shown here) seem to be highly similar to those of the single-step model, suggesting that genomic selection started only recently in Hanwoo, with no clear impact. Figure 6 shows the genetic trends of BF in Hanwoo steers. The decreasing trend in trait BF indicates favorable genetic  progress since the update of the data-recording scheme in 2008. A flat genetic trend was observed before the birth year 2007/2008. We also observed a nearly flat genetic trend curve in the recent birth years from 2016 and later, reflecting less emphasis on the selection objective imposed on trait BF. As shown in Figure 5 for trait MS, we observed a slightly lower genetic trend for the single-step evaluation with the truncated dataset than with the full dataset. In comparison to trait MS, less genetic progress has been made in trait BF, which can be explained by the lower heritability value of BF and lower weight of BF in the carcass index of Hanwoo. Similar genetic trends were found in steers for traits CW, EMA, and MS (data not shown). For the CW trait, a larger trend was observed in the single-step evaluation, notably with the truncated dataset, than with the full dataset. Compared to steers with carcass data, cows with their own carcass records seem to have less genetic progress in trait MS but similar progress in traits BF, EMA, and CW.

Correlations of predictions between the full and truncated evaluations
The GEBV correlation between the full and truncated evaluations of data from 4 yr ago indicates the predictive ability of the single-step genomic model. Figure 7 shows the correlations of the GEBV or conventional BLUP EBV of steers with the CW records between the truncated and full evaluations for the single-step model (solid blue line) or the conventional BLUP model (dotted red line). Steers born in 2016 or before have CW records in both the truncated and full datasets, and the correlations of GEBV or EBV between the two evaluations are unity for either model. Steers born in 2017 and later had CW records in the full evaluation; however, no data in the truncated evaluation showed a much lower GEBV or EBV correlation, as low as 0.70. Despite a smaller difference in the correlations, the GEBV of the truncated evaluation were consistently higher and more correlated with the full evaluation for the single-step model than those for the conventional BLUP model. The higher GEBV correlation indicated a higher prediction accuracy of the single-step model than that of the conventional BLUP model, despite the relatively small difference in the correlations owing to the limited number of genotyped animals with CW records. Similar correlations were observed for the other three traits (BF, EMA, and MS) for steers with carcass records.
As shown in Table 1, some of the 29,965 genotyped animals, including those born in the early years, had no carcass records. Figure 8 shows the GEBV and EBV correlations between the full and truncated evaluations for the single-step model (solid blue line) and conventional BLUP model (dotted red line), respectively. Across all birth years of the genotyped animals, the single-step model showed a markedly higher GEBV correlation between the two evaluations than the conventional BLUP model. Only in the birth years 2014 through 2016 did the two models give equal GEBV or EBV correlations between the two datasets because these three birth years have the most genotyped animals with their own carcass records in either the full or the truncated evaluation. The higher correlations between the full and truncated evaluations indicated a higher accuracy of genomic prediction and a more stable genomic prediction over time for the single-step model.