Wagyu Feeder Check: A genomic-based tool to identify performance differences of Australian Wagyu and Wagyu crossed cattle

Wagyu fullblood (FB) and crossbred (XB) population details

The Wagyu Feeder Check is based on a reference population of 8316 genotyped and phenotyped Australian fullblood Wagyu (FB, N = 2120) and Wagyu-crossed (XB, N = 6196) cattle, principally Wagyu × Angus F1 animals. Feedlot and carcase records were supplied by independent commercial supply chains from eight populations including three FB populations (POP1, N = 1455; POP2, N = 477; and POP3, N = 188) and five XB populations (POP4, N = 1049; POP5, N = 1285; POP6, N = 2456; POP7, N = 654; and POP8, N = 752). Animals were slaughtered from 2013 to 2022, with XB animals slaughtered only in 2021 (N = 2098) and 2022 (N = 4098). Phenotypes included daily gain at feedlot finishing (FADG), carcase weight (CWT), carcase eye muscle area (CEMA), carcase AUS-MEAT marbling score (MARB) and carcase rump fat at the P8 site (CP8). To accommodate high marbling content, MARB was measured using a modified AUS-MEAT scoring system (AUS-MEAT 2005), which ranges from 1 (nil) to 12 (abundant) in increments of 1.

For the analysis of phenotypes, a contemporary group (CG) was defined as a combination of population of origin (eight levels combining feedlot and abattoir), sex (two levels), and kill date. There were 71 and 13 kill dates for the FB and XB populations, respectively. For the FB populations all months were represented, whereas for the XB population all months were represented except for December and January. Initial edits aimed at removing animals without genotypes or from CG with less than three individuals. For the FB population, there were 82 CG with an average of 25.8 cattle and ranging from 3 to 145. For the XB population, there were 29 CG with an average of 213.6 cattle and ranging from 8 to 576.

Genotypes and genomic relationships

Genotypes for 82 504 autosomal single nucleotide polymorphisms (SNP) were available for all 8316 animals included in this study and used to compute the genomic relationship matrix (GRM) following Method 1 of VanRaden (2008) with the modification of Karoui et al. (2012) to make it invertible:

where S is the centred matrix relating SNP genotypes (recoded as 0, 1 or 2) in columns with animals in rows, and p_i is the frequency of the second allele of the i-th SNP, and I is an identity matrix included to make GRM invertible by enlarging the diagonal elements.

To obtain a measure of the genomic similarity between the two populations, FB and XB, we explored the SNP allele frequencies, the values of the GRM and performed a principal components analysis (PCA) based on a singular value decomposition of the GRM (Misztal and Legarra 2017). Additionally, using a smaller panel of 27 883 SNPs and numerical approaches outlined in Reverter et al. (2020) we estimated the genomic breed composition of FB and XB cattle across 10 breeds including: Wagyu, Angus, Brahman, Charolais, Hereford, Holstein, Limousin, Santa Gertrudis, Shorthorn and Simmental.

Finally, the genomic relationships among all individuals were processed following the network-based Pedigromics pipeline (Reverter et al. 2019) by establishing network connections after considering genomic relationships ≥0.125 corresponding to the equivalent of a great-grandparent to great-grand offspring relationship.

Genomic predictions and cross-validation accuracy

Variance components, heritability (h²), genetic (r_g) and residual (r_e) correlations were estimated based on GBLUP methodology (genomic best linear unbiased prediction) using the Qxpak5 software (Pérez-Enciso and Misztal 2011). For the genomic prediction models, we carried out GBLUP analyses including one pentavariate analysis with the entire dataset to produce the most accurate genomic predictions to use as reference; and a series of 20 univariate analyses from the cross-validation datasets (with 20 from five traits by four validation groups). In all cases, the GBLUP models contained the fixed effects of CG and the linear regression covariates of slaughter age (in days), and the first three principal components of the GRM. Fitting the principal components aims at accounting for hidden population structures, likely with FB and XB populations, that could have been missed if only fitting CG. Additionally, the random additive polygenic and residual effects were fitted in the GBLUP models with assumed distributions N(0, G⊗V) and N(0, I⊗R), respectively, where G represents the GRM described earlier, V is the genetic covariance matrix, I is an identity matrix, R is the residual variance–covariance matrix and ⊗ represents the Kronecker product.

Firstly, the resulting GEBV from the pentavariate analysis were termed u^w to indicate that they are based on the whole dataset. Secondly, for the cross-validation of genomic predictions, we created four cross-validation datasets each with the phenotypes from a random 1549 (or 25%) XB cattle set as missing. In each cross-validation schema, the resulting GEBV were termed u^p to indicate that they are based on partial data.

Finally, traditional (Bolormaa et al. 2013) and LR method (Legarra and Reverter 2018) approaches were used to estimate accuracy, bias, and dispersion of GEBV. The following four metrics were employed:

For bias and dispersion, we constructed 95% confidence intervals based on ± 1.96 s.e. around the observed means across the 20 scenarios, i.e. 5 traits × 4 validation datasets.

Phenotypes

Table 1 provides summary statistics for all raw and unadjusted phenotypes, and the slaughter age covariate used in the analyses. Based on the fixed effects model used to adjust phenotypes, the CG effects and covariates combined accounted for 33.5%, 37.1%, 30.7%, 30.8% and 45.2% of the variation in FADG, CWT, CEMA, MARB and CP8, respectively, and with all effects being highly significant (P < 0.001) for all traits, except for slaughter age (P > 0.1) for CEMA, which was likely captured by the effect of CG. On average and compared to FB, XB cattle were 15 days older, with 6.5% heavier carcases and 27% less marble, but 46% more subcutaneous fat. For the XB cattle, the values in Table 1 are very similar to those reported by Connolly et al. (2019) with a Wagyu crossbred population with a Wagyu content estimated at 73.0 ± 0.7%. In that work, and after 163 days on feed and an average age at slaughter of 1147 days, the authors reported means ± s.e. for CWT, CEMA, MARB and CP8 of 434 ± 5.2 kg, 41.4 ± 0.85 cm², 5.9 ± 0.52 scores and 23.7 ± 1.01 mm, respectively.

With respect to FADG, our averages of 0.83 kg/day for FB and of 0.92 kg/day for XB are comparable to those from Vázquez-Mosquera et al. (2022), who reported daily weight gains of 0.916 kg/day for Wagyu purebred steers (N = 262) and 1.046 kg/day for Wagyu × Angus crossbred steers (N = 103) during the weaning to growing period, and 0.628 kg/day for purebred and 0.640 kg/day for crossbred during the growing to fattening phase. Similarly, the study of Alexandre et al. (2021) with 3408 Australian Angus steers slaughtered at an average of 734.53 days reported means ± s.d. for FADG, CWT, and CEMA of 1.59 ± 0.33 kg/day, 432.99 ± 65.60 kg and 90.06 ± 10.86 cm², respectively. More recently, Reverter et al. (2023) with a population of 2120 Angus-based steers from four commercial feedlots feed for an average of 222.46 day, reported means ± s.d. for CWT and MARB of 425.45 ± 43.94 kg and 3.64 ± 1.33, respectively.

Genotypes and genomic relationships

On average, the estimated percentage of genomic Wagyu content in the FB and XB cattle was 99.12% and 59.55%, respectively. Also on average, the sum of Wagyu and Angus content in the XB cattle was 94.12% and with the remaining 5.88% estimated to originate from mostly Holstein, Santa Gertrudis, and Shorthorn with more than 1% (Table 2). Given the large content of Angus in the XB population, caution will be needed when exploring the applicability of the Wagyu Feeder Check tool in Wagyu × Bos indicus crossbred animals. In particular, among the XB cattle, there were 129 (or 2.1%) for which the summed Wagyu and Angus content was <50%. For these animals, the average Wagyu and Angus content was 42.63% and 2.33%, respectively. Therefore, these 129 XB cattle could be considered F1 Wagyu × non-Angus. Other significant breed percentages represented among these 129 XB cattle were Shorthorn (17.80%), Santa Gertrudis (12.90%), Charolais (7.07%), Holstein (6.21%) and Hereford (5.46%).

The percentage of variation in genomic relationships accounted for the first three principal components was 4.97%, 1.02% and 0.90%, respectively. Fig. 1 displays the scatter plot of the PCA for the first two principal components (PC1 and PC2) with the two populations (FB and XB) distinctly highlighted. Also highlighted in Fig. 1 are the 129 XB cattle with low Angus content. Explaining 4.97% of the variation in genomic relationships, PC1 shows a clear separation between FB and XB. On closer examination, we found a very strong Pearson correlation coefficient (r) between the genomic estimate of Wagyu content and PC1 (r = 0.976 ± 0.002). This correlation remained strong when examined within the FB (r = 0.687 ± 0.016) and within the XB (r = 0.959 ± 0.004) populations. In multibreed beef populations, it is not uncommon for PC1 to distinguish between Bos taurus and Bos indicus cattle. For recent examples, see for instance Yonesaka et al. (2016) and Porto-Neto et al. (2023). The work by Yonesaka et al. (2016) is of particular relevance because although their PC1 explained 17.6% of the variation, indeed separating Bos taurus from Bos indicus population, PC2 explained 4.5% (very close to our 4.96%) and distinguished Japanese Black cattle in one extreme and Angus on the other.

Fig. 1.

Scatter plot of the first two principal component, PC1 and PC2, of the genomic relationship among individuals in the Wagyu fullblood (FB; black dots), Wagyu crossbreds (XB; red dots) populations, and in the 129 XB cattle with low Angus content (blue inverted triangles).

Fig. 2a shows the structure of the Wagyu Feeder Check population using the Pedigromics network approach (Reverter et al. 2019). The centre of the network is dominated by the three FB populations whereas the five XB populations, with smaller node sizes, are scattered around the periphery. To better illustrate the interconnectivity within and across populations, the insert in Fig. 2b shows a subset of the whole Pedigromics network where only the 10 most connected animals within each of the eight populations are displayed. As expected, due to the larger size of POP1 among the FB, most of the connectivity from XB to FB happens through POP1. This interconnectivity within and across populations is further exemplified in Fig. 2c where a heatmap of genomic relationships >12.5% is tabled. The diagonals are showing the connectivity within a population and populations are more highly correlated within themselves than across populations. Although there’s little connectivity among the three FB populations (POP1, POP2 and POP3) including no connectivity between POP2 and POP3, all five XB populations (POP4 to POP8) have connections with the FB populations. Also, the XB population with the highest Wagyu content (POP7 and 84.4% Wagyu) has the highest connectivity with the three FB populations. Importantly, the Wagyu content was not estimated based on connections to the FB populations.

Fig. 2.

(a) Structure of the Wagyu Feeder Check population using the Pedigromics network approach (Reverter et al. 2019), where animals from the eight populations (POP1 to POP8) are nodes, genomic relationships >12.5% are edges and the percent Wagyu content is mapped to the size of each node. (b) A subset of the whole Pedigromics network where only the 10 most connected animals within each of the eight populations are displayed. (c) Interconnectivity within and across the eight populations (POP1 to POP8) including number of records (N Rec), percentage of Wagyu content (% Wag) and based on genomic relationships >12.5%.

Upon closer inspection of the GRM, the proportion of genomic relationships >25% (indicating potential half-sibs) was 5.30%, 11.82%, 6.79%, 1.30%, 0.73%, 0.64%, 0.61%, and 0.27% for POP1 to POP8, respectively. Therefore, the proportion of half-sibs within each XB population (POP4 to POP8) does not seem to be higher than that observed within the FB populations (POP1 to POP3).

Genetic parameter estimates

Table 3 shows the estimates of heritability (h², ± standard error (s.e.)), genetic (r_g) and residual correlation (r_e) obtained from the pentavariate GBLUP model with the entire dataset. Heritability estimates were generally moderate to high, ranging from 0.347 ± 0.021 for CEMA to 0.526 ± 0.018 for FADG. Our h² estimates are very similar to those recently reported by Rostamzadeh Mahdabi et al. (2023) using a population of 9850 Wagyu steers and heifers. In that work, the authors reported h² estimates of 0.510 ± 0.014 for CWT, 0.430 ± 0.015 for CEMA, 0.473 ± 0.015 for CP8 and 0.486 ± 0.014 for MARB. In Japanese Black cattle, Onogi et al. (2014) reported h² for CWT (0.56), CEMA (0.43) and MARB (0.66), which are higher than those from our study. Similarly, h² estimates for carcase traits in Korean Hanwoo cattle have been reported recently. Estimates from Mehrban et al. (2019) include CWT (0.39), CEMA (0.45), MARB (0.64) and CP8 (0.51); and estimates from Naserkheil et al. (2021) were 0.42, 0.50, 0.59, and 0.56 for CWT, CEMA, MARB and CP8, respectively.

The review of Ríos Utrera and Van Vleck (2004) reported average h² estimates for CWT, CP8, CEMA and MARB of 0.40, 0.36, 0.40, and 0.37, respectively. Therefore, at 0.347 ± 0.021 our h² estimate for CEMA is somewhat lower than what has been generally reported. On the other extreme, at 0.526 ± 0.018 our h² estimate for FADG is higher than published values including the 0.33 of Torres-Vázquez et al. (2018) and the 0.30 of Alexandre et al. (2021) both with Angus cattle, and the 0.31 of Somavilla et al. (2017) with Nellore cattle.

Estimates of r_g were strong and positive between FADG and CWT (0.832 ± 0.012), and moderate and positive between CWT and CEMA (0.391 ± 0.045) and between CEMA and MARB (0.333 ± 0.048). For the same three pairs of traits, the study of Alexandre et al. (2021) with Angus cattle reported r_g estimates of 0.65, 0.37 and 0.14, respectively; with Hanwoo cattle, Naserkheil et al. (2021) estimated r_g of 0.56 between CWT and CEMA and 0.35 between CEMA and MARB. In general, all other r_g estimates were within two s.e. and likely not significantly different from zero. Also, in general, the estimates of r_e were lower and closer in magnitude to zero than the r_g estimates.

Heritabilities, genetic and residual correlations are trait- and population-specific parameters; therefore, a diversity of estimates can be observed in different studies even with the same breed and trait. Factors such as breed, environment, age at measurement, accuracy of measurement and data source (i.e. pedigree or genomic), as well as the analytical methodology may affect the estimates. Because of this diversity, it is always worthwhile to explore the quality of the resulting genomic predictions in terms of their accuracy, bias, and dispersion. This is particularly the case when, as in the current study, the main objective is to develop a genomic tool to be deployed at scale in commercial scenarios.

Genomic predictions and cross-validation results

Shown in Fig. 3 are the genomic prediction accuracies for all traits averaged across the four cross-validation samples and based on both accuracy metrics, ACC_T and ACC_LR. Also in Fig. 3 are the ACC_T for the 129 XB cattle with low Angus content. Across the 20 accuracy estimates obtained from five traits and four validation datasets, the correlation between ACC_T and ACC_LR was 0.922, an estimate higher than the 0.73 reported by Alexandre et al. (2021) across 49 estimates, and the 0.831 reported by Reverter et al. (2021) across 15 estimates and both studies with feedlot and carcase traits in Angus cattle.

Fig. 3.

Genomic prediction accuracy based on traditional accuracy (ACC_T) and LR method accuracy (ACC_LR) across the entire population and ACC_T in 129 XB cattle with low Angus content for average daily gain during feedlot finishing (FADG), hot carcase weight (CWT), carcase eye muscle area (CEMA), AUS-MEAT marbling score (MARB), and carcase subcutaneous fat depth at the rump or P8 site (CP8).

Accuracies were in the range of 55–65% for all traits except for CEMA, which were in the 40–50% range. This lower accuracy for CEMA was attributed to its lower h² of ~0.35 compared to h² > 0.45 for the other traits (Table 3). Our GEBV accuracies are consistent with the theoretical expectation of ~55% and ~65% respectively for traits of 0.3 and 0.5 heritability and a reference population of ~7000 individuals (Goddard and Hayes 2009), which is the equivalent of 8316 total animals in our case minus the ~1500 XB animals set aside for each validation. If a whole sub-population of XB had their phenotypes set to missing the expectation is that the ACC would be diminished, particularly if that sub-population is poorly connected. However, all XB populations have connections among themselves and among all three FB populations (Fig. 2c). In addition, because most sub-populations of XB have less than 1549 animals or 25% of XB cattle (the exception being POP6 with 2456 cattle), setting their phenotypes to missing would imply that the remaining ‘reference’ is larger which might result in higher genomic prediction ACC. This phenomenon has been recently reported by Reverter et al. (2023) where data from four commercial feedlots were removed one feedlot at a time and compared to removing 50% of all the commercial data.

The study of Onogi et al. (2014) with Japanese Black cattle reported a predicted ability (correlation between the GEBV and the adjusted phenotypes) of 0.44, 0.42 and 0.39 for CWT, CEMA and MARB, respectively. Meanwhile, the study of Mehrban et al. (2019) with Hanwoo cattle estimated GEBV accuracies of 0.56, 0.44, 0.36 and 0.33 for CWT, CEMA, MARB and CP8, respectively.

For the 129 XB cattle with low Angus content GEBV accuracies were not markedly affected for FADG and CP8, whereas they were somewhat lower for CWT, and higher for CEMA and MARB. This rather unexpected result was attributed to the low number of animals affecting the precision of the ACC estimates, and to the large proportion of Wagyu content in these cattle (average 42.6%, range = 0–49.9%). Nevertheless, we maintain that caution will be needed when exploring the applicability of the Wagyu Feeder Check tool in Wagyu × non-Angus, particularly Wagyu × Bos indicus crossbred animals.

Averaged across all validation datasets, and indicating no bias, the 95% confidence interval for GEBV bias contained zero in all traits (Table 4, upper part). Similarly, the 95% confidence intervals for the dispersion in GEBV contained zero for all traits except FADG for which the 95% confidence interval showed a tendency for under-dispersion (or deflation) of GEBV (Table 4, lower part). Nevertheless, this under-dispersion vanished at 99% confidence interval.

Unlike the recent work of Reverter et al. (2023) with 3766 purebred Angus steers plus 2124 Angus-based commercial steers, in the present study we did not observe a tendency for over-dispersion (or inflation) of GEBV for CWT or MARB. Similarly, the recent work by Koo et al. (2023) exploring the quality of the genomic evaluation of Korean Hanwoo cattle based on a large number of genotyped cows, steers and young animals applied the LR method to confirm a slight negative bias (GEBV overestimation) for all traits and animal groups, but a slight positive bias (GEBV underestimation) for CEMA and MARB in steers. However, their estimates of bias did not deviate significantly from zero. Similarly, their LR estimate of dispersion did not deviate markedly from the expected value of 1 suggesting that GEBV were neither inflated nor deflated.

Like the case for genetic parameters, estimates of GEBV quality metrics (i.e. bias, dispersion, and accuracy) can be affected by a myriad of reasons including the size of the reference population, heritability of the trait, relatedness between reference and validation population, and marker density. Again, because of this diversity, it is always worthwhile to explore the quality of the resulting GEBV in terms of their ability to reflect phenotypic differences in the highest and lowest GEBV quartile.

Table 5 presents phenotypic values for validation XB animals assigned to quartiles on the basis of GEBV and the difference (Q1mQ4) between the top (Q1) and the bottom (Q4) quartile. In agreement with the theoretical expectation of a linear relationship between ACC and phenotypic differences in predefined percentiles (Reverter et al. 2022), we observed a Q1mQ4 difference of 0.154 kg/day in FADG, 40.92 kg of CWT, 4.368 cm² of CEMA, 1.419 scores of MARB and 5.325 mm of CP8.