Parameter estimation and assessment of bias in genetic evaluation of carcass traits in Hanwoo cattle using real and simulated data

A preliminary analyses of variance was performed using the linear model in R [23] to identify the most appropriate fixed effects (contemporary groups) to be used in the subsequent analyses for all studied traits. The most significant factor was the interaction between herd (3,646 herds), birth-year (7-levels), and birth-season (4-levels), which were combined to form the contemporary group. The final dataset consisted of 469,002 animals from 31,403 contemporary groups (Table 1). A contemporary group of herd-year-season (HYS) with less than five records was omitted from the analysis. The distribution of animals across various ages and herds, as well as animals born in each year, and season, are shown in Fig. 1. The description of population structure, as well as the distribution of animals across age, herd, birth year, and season are summarized in Table 1 and Fig. 1.

Fig. 1. Distribution of animals across age, birth year, birth season, and herd top panel, and distribution of animals across contemporary group (HYS), lower panel. HYS, herd-year-season.

Download Original Figure

All individuals were slaughtered at various ages (28–35 months), and phenotypic measurements were taken on the chilled carcass. Phenotypic data on carcass traits included CWT, eye muscle area (EMA), backfat thickness (BFT), bodyweight (BW) at slaughter, and meat-index (MI). EMA and BFT were measured at the 12th and 13th rib junction after a 24-h chill, and cold CWT measurement was taken at that time. The fasted live BW in kilogram was taken prior to transport to the abattoirs. Only 52% of animals in the dataset had BW records at slaughter. The MI represents the retail cut percent, which is predicted from a linear index of CWT (kg), EMA (cm²), and BFT (mm). MI=68.184 − [0.625 × BFT] + [0.130 × EMA] − [0.024 × CWT] + 3.23.

Meat quality traits included MS, meat colour (MC), FC, and meat texture (MT) which were each graded at the time when carcass traits were measured. MS was recorded based on visual inspection by trained technicians using the Beef Marbling Standard (BMS) for grading the carcass. The BMS system classifies the meat from 1 to 9 scale of marbling, with a 2% intramuscular fat content difference between each point score [24]. Similarly, trained technicians graded the other three traits (MC, FC, and MT) manually. MC was assessed and graded from very light red (grade 1) to dark red (grade 7). Similarly, the FC was assessed and graded from polar white (grade 1) to creamy yellow (grade 6). The texture of the meat was evaluated on a scale from very fine (grade 1) to coarse (grade 3) (Table 2).

Variance components and heritabilities were estimated with a series of univariate animal models. For each trait, the model included the contemporary group (as defined above) fitted as a fixed effect and age fitted as a linear covariate in the model. Genetic and phenotypic correlations were estimated using a bivariate model with similar fixed and random effects as the univariate model. ASReml version 4.1 software [25] was used for the entire data analysis. Prior to the main analysis, we performed a model comparison between a model with and without maternal (permanent environmental) effect. Based on the Akaike information criterion (AIC) value, the model without maternal effect best fitted the data [26]. The model was represented as:

where X is an incidence matrix for observations y. Age as a linear covariate and contemporary groups (HYS) as a fixed effect in vector b, Z is an incidence matrix for random animal additive genetic values in u, and e represents random environmental effects. u and e were assumed to be distributed as $u|{\sigma }_a^2\sim N\left(0,A{\sigma }_a^2\right)$, $e|{\sigma }_{є}^2\sim N\left(0,I{\sigma }_{є}^2\right)$, respectively where ${\sigma }_a^2$ is the additive genetic variance and ${\sigma }_e^2$ is the residual. Where A is a numerator relationship matrix for all animals using 13 generations and I is an identity matrix.

Linear mixed models were used to estimate the breeding values and variance components [16] implemented in ASReml version 4.1 [25] using the model described in equation 1. The fixed effect (contemporary group) was fitted only to the LM trait in the sorted scenario. The bias in the variance of the estimated breeding values (EBV) was measured through the regression coefficient (slope) of the true breeding values of 100 base sires on their EBV in each replication. The estimated bias in EBV was the mean of the 10 replications.

As expected, the random allocation of animals to contemporary group was unbiased; however, bias was observed when animals were sorted into contemporary group based on EM traits. The degree of bias varied depending on the heritability and the correlations between EM and LM traits.

As indicated in Fig. 2A, the EBV was biased (overestimated) when the genetic and residual correlations between the simulated traits were −0.8 and 0.8 respectively. Slight underestimated EBV was observed when the genetic and residual correlation between the EM and LM traits was higher (0.8). These results indicated that the bias was higher for the sorting scenario when the correlation between the two studied traits was higher. But slightly biased EBV was observed when the genetic and residual correlation was lower (−0.2 and 0.2) (Fig. 2A). The result showed that changes in genetic and residual correlations affected the magnitude of bias in the EBV. The current results showed that the EBV was highly biased when the correlation between the sorting (EM) and the subsequent evaluation (LM) trait was higher, and sorting bias was lesser when the correlation between the two traits was weaker. At higher genetic (0.8) and residual (−0.8) correlations, considerably overestimated EBV with a bias of 0.54 was observed (Fig. 2D). However, less biased EBV was observed when the residual correlations were lower (Fig. 2D). Overall, overestimated, and underestimated breeding values were observed at the combination of high genetic and residual correlations (Fig. 2D). With similar genetic and residual correlations, a slight increment of bias was observed when the heritabilities of both simulated traits changed from high (0.5) to low (0.2).

Fig. 2. Change in bias due to sorting of animals and impact of the correlations and heritabilities of the EM and LM traits in genetic evaluation. Used genetic correlations on the X-axis and bias on the Y-axis. h², heritability; HH, high-high (0.5); HL, high-low (0.5; 0.2); LH, low-high (0.2, 0.5); LL, low-low (0.2, 0.2); re, residual correlations.

Download Original Figure

We assessed the bias generated by sorting animals in the genetic evaluation using traits with different heritabilities (high-low and low-high). In the randomly sampled scenario, the regression coefficients (slopes) in all alternatives were equal to one and no change of bias was observed with changing of input parameters including heritability of the traits (results not shown). However, we observed different magnitudes of bias with the changing of input parameters between the two studied traits when animals were sorted into various contemporary group depending on the EM trait (yearly weight). Highly biased EBV of 0.15 (±0.01) and 0.20 (±0.02) were observed when genetic and residual correlations were 0.8 and −0.8, respectively (Fig. 2B). These highly biased EBV (overestimated) were observed when the heritability of the LM trait (evaluation trait) was low (0.2). However, relatively less biased EBV were observed when the heritability of the LM trait was higher (Fig. 2C). The highly biased EBV with low heritability of LM trait indicated that at constant genetic and residual correlations, genetic evaluation is further affected by the heritability of the evaluation (LM) trait than the sorting (EM) trait. Relatively less biased EBV was observed at low correlations with both high and low heritability of the evaluation (LM) trait (Figs. 2B and 2C). The low genetic and residual correlations may be responsible for the less biased EBV even though the heritability of the LM trait varied. Overall, with the low heritability (0.2) of the LM trait, the EBV was highly biased (overestimated), however, when the heritability of the LM trait was higher (0.5) the EBV was less biased. In addition, the bias of EBV was less when the difference between the genetic and residual correlations was lower whereas highly biased EBV was observed when the difference between the two correlations was bigger. Overall, our results indicate that the magnitude of the sorting bias depends on the genetic architecture of the two traits.

The phenotypic data used in the current study was obtained from slaughterhouse and these carcass performance data were recorded on various age groups of animals ranged between 28 and 35 months. The estimated parameters in the current study are closer to previously reported estimates obtained from progeny testing experiment. This is indicating that commercial data obtained from slaughterhouse are useful and can be used to estimate genetic parameters in the evaluation of carcass traits. The commercial carcass performance data can be available in large amount from slaughterhouse compared to progeny testing experiment data, but we have observed some drawback in the commercial data. This is because, in most commercial farms, the desired carcass quantity and quality is specified by abattoirs based on the weight of animals, which is suitable for slaughter. As a result, most animals have arrived in abattoirs in homogenized weight which is called harvesting effect. Such selective formation of contemporary groups might lead to biased evaluation of genetic merit of animals for carcass trait. However, this type of bias in the genetic evaluation can be controlled using various models depending on the source of the bias. For instance, bias due to sorting of animals based on weight can be reduced or eliminated using multitrait genetic evaluation methods. The multitrait model can account for harvesting effect and thus produce more accurate and less biased parameter estimation in genetic evaluation.

Sorting of animals based on yearling weight potentially affected the genetic evaluation of the latter measured traits. The current study has verified that the non-random evaluation of animals led to a bias in genetic evaluation.

The sorting bias observed in the current simulation study has a similar effect to the selection bias described in the literature [14,15,33]. The theory of selection bias was established clearly by Pollak et al. [15]. They demonstrated that potential bias was found in the genetic evaluation of the second trait, which was selected based on the first trait. Similarly, in the current simulation study, when animals were sorted based on the early measured trait, sorting bias was detected in the subsequent evaluation of the latter measured trait. Also, Pollak and Quaas [33] observed a selection bias in the predictors in single-trait analysis and the magnitude of the bias of the EBV depends on the correlation among studied traits which is in agreement with the current findings.

A recent study by Macedo et al. [34] has studied a selection bias in genetic evaluation models. They found a selection bias in genetic evaluation using data with environmental trends compared to randomly sampled data into different contemporary groups. Similarly, in the current study, we found biased EBV when animals were non-randomly assigned into various contemporary group. Selection bias has been studied in sheep, Eiríksson and Sigurdsson [35] demonstrated that a selected group of the lambs was kept for replacement and therefore not measured for carcass traits, led to bias in the genetic evaluations. The authors demonstrated that the bias was higher for the selected ram group, indicating that this group consisted of rams that had higher genetic merit for carcass conformation trait than the rest groups. Similarly, we found biased EBV when animals with higher yearling weight were assigned into same group.

In the current simulation study, the structures of input parameters including heritability, genetic and residual correlations, and the difference between the two correlations influenced the genetic evaluations. The genetic and phenotypic correlation of traits potentially affects the genetic evaluation of animals in breeding program [15,33,36,37]. Our results showed that highly biased EBV was observed when the genetic correlations were higher which coincided with a report by Pollak and Quaas [33]. In addition, our result showed that highly biased EBV was observed when the difference between genetic and residual correlation was higher. Conversely, when the difference between genetic and residual correlation was smaller, the EBV was less biased and the regression coefficient was closer to one.

In the current simulation study, the bias of the EBV was assessed by allocating unequal heritabilities (high-low and the vice versa) for the two simulated traits. In this case, the observed bias of the EBV was not the same as a similar heritability (high-high and low-low) was assigned to the studied traits. Highly biased EBV was observed when the heritability of the LM (evaluation) trait was lower. This could be explained as genetic evaluation is more affected by the heritability of the LM trait than the EM (sorting) trait with the constant genetic and residual correlations. A recent study showed that allocating wrong heritability for the trait in a genetic evaluation led to a biased EBV in simulation study [34]. These authors established that in pedigree-based predictions, the use of incorrect heritability generates a strong bias in genetic evaluation using simulation data. Our result showed that different magnitude of bias was observed with two different heritabilities assigned to the LM (evaluation) trait. However, we have not proven that the observed bias linked to incorrect heritability that was allocated to the evaluation trait during data simulation. Knowing the source of bias and then using appropriate models for the genetic evaluation of animals is a key component of the process to estimate accurate breeding values. For instance, to overcome selection bias in genetic evaluation, multivariate evaluations methods had been proposed [16]. The decision to use a multitrait model versus a single trait model depends on correlations among the studied traits. In the first scenario of the current study, the sorting bias that was observed in the single-trait model was reduced by the bivariate model (results not shown). Similarly, Pollak et al. [15] showed that the observed bias associated with selection in the univariate model was reduced or eliminated by multiple traits evaluation methods.

A recent study by Macedo et al. [34] has studied a selection or sorting bias in genetic evaluation models. They found a selection or sorting bias in genetic evaluation using data with environmental trend compared to a randomly sampled data into different contemporary groups. Similarly, in the current study, we found biased EBV when animals were non-randomly assigned into various contemporary group. Selection bias has been studied in sheep, Eiríksson and Sigurdsson [35] demonstrated that a selected group of the lambs was kept for replacement and therefore not measured for carcass traits, led to bias in the genetic evaluations. The authors demonstrated that the bias was higher for the selected ram group, indicating that this group consisted of rams that had higher genetic merit for carcass conformation trait than the rest groups. Similarly, we found biased EBV when animals with higher yearling weight were assigned into same group.