Comparison of breeding value prediction for two traits in a Nellore-Angus crossbred population using different Bayesian modeling methodologies

The objectives of this study were to 1) compare four models for breeding value prediction using genomic or pedigree information and 2) evaluate the impact of fixed effects that account for family structure. Comparisons were made in a Nellore-Angus population comprising F2, F3 and half-siblings to embryo transfer F2 calves with records for overall temperament at weaning (TEMP; n = 769) and Warner-Bratzler shear force (WBSF; n = 387). After quality control, there were 34,913 whole genome SNP markers remaining. Bayesian methods employed were BayesB (π̃ = 0.995 or 0.997 for WBSF or TEMP, respectively) and BayesC (π = 0 and π̃), where π̃ is the ideal proportion of markers not included. Direct genomic values (DGV) from single trait Bayesian analyses were compared to conventional pedigree-based animal model breeding values. Numerically, BayesC procedures (using π̃) had the highest accuracy of all models for WBSF and TEMP (ρ̂gĝ = 0.843 and 0.923, respectively), but BayesB had the least bias (regression of performance on prediction closest to 1, β̂y,x = 2.886 and 1.755, respectively). Accounting for family structure decreased accuracy and increased bias in prediction of DGV indicating a detrimental impact when used in these prediction methods that simultaneously fit many markers.


Introduction
Expected progeny differences (EPD) are widely used for selection in the US beef cattle industry. Current international research efforts are focused on incorporating genomic prediction methods into national evaluations that produce EPD (Hayes et al., 2009b;VanRaden and Sullivan, 2010;Saatchi et al., 2011Saatchi et al., , 2012Saatchi et al., , 2013Rius-Vilarrasa et al., 2012). The landmark paper by Meuwissen et al. (2001) described genomic predictions using Bayesian inference that they believed capitalized on linkage disequilibrium (LD) between observed markers and unobserved quantitative trait loci (QTL).
For these methods, prediction procedures can be a two-part process using the Bayesian framework of 1) training to calculate substitution effects of markers in historical data using both phenotypes and genotypes; and 2) using markers genotyped in individuals outside the training data to calculate direct genomic values (DGV) based on substitution effects estimated in training. The size of the refer-ence (i.e., training) population and its relationship to the prediction population are important considerations to achieve higher precision (accuracy) of genomic breeding values in animals with only genomic information available (de Roos et al., 2009;Ibáñez-Escriche et al., 2009;Habier et al., 2010Habier et al., , 2013Kizilkaya et al., 2010;Saatchi et al., 2010;Toosi et al., 2010).
Relatedness among animals and between breeds (in crossbred or composite populations) can directly impact prediction accuracy (de Roos et al., 2009;Ibáñez-Escriche et al., 2009), at least in part due to inconsistent LD patterns among cattle breeds (Gautier et al., 2007;Bovine HapMap Consortium, 2009). Few studies have investigated the use of genomic prediction models using field or research data from crossbred populations (Hayes et al., 2009a;Snelling et al., 2012). The objectives of this study were to 1) evaluate and compare four models (three Bayesian using genomic information and a pedigree-based animal model) in predicting genetic merit using all animals in the training population; and 2) evaluate the effects on genetic merit prediction of including or excluding fixed effects to account for family structure. These objectives were investigated using single trait analyses of two traits measured in a Nellore-Angus crossbred population.

Cattle, genotypes, and traits
The cattle population, genotypes, and traits used in this study were previously described in Hulsman Hanna et al. (2014). Briefly, the cattle population consisted of Nellore-Angus F 2 (embryo transfer and natural service reciprocal crosses), F 3 and natural service half-siblings to embryo transfer F 2 animals. The population was characterized by 3 distinct cycles (designated as Cycle 1, 2, and 3), where Cycle 1 included embryo transfer Nellore-Angus F 2 and natural service half-sibling calves, Cycle 2 included natural service F 2 calves, and Cycle 3 included natural service F 3 calves. All procedures involving animals were approved by the Texas A&M Institutional Care and Use Committee. Animals were genotyped on the Bovine SNP50 Version 1 & 2 assays (Illumina Inc., San Diego, CA) and, after quality edits, had 34,913 SNP available for genomic prediction analyses.
Overall temperament at weaning (TEMP), a subjective measurement on a 1 to 9 scale, where 1 is docile and 9 is wild or unruly (n = 772) was used in this study. Four evaluators scored each animal in this study for TEMP, but average TEMP score across all evaluators was used in subsequent analyses. Warner-Bratzler shear force (WBSF), a quantitative measurement using a United Testing machine (United 5STM-500, Huntington Beach, CA) with a 11.3 kg load cell and a v-notch WBSF attachment, was used to measure the peak force required to shear each core and reflect meat tenderness. Only electrically-stimulated data was available for all steers in Cycles 1, 2, and 3 and were therefore used. The average shear value of the electrically stimulated side of each steer, averaged across all shears per steak, for 14-d aged steaks (n = 390) were used for analyses in this study.

Statistical analysis
In a previous study with these animals, Hulsman Hanna et al. (2014) fitted a model for TEMP that included fixed effects of sex (n = 2), family nested within sire (n = 31), birth year-season combinations (n = 10), and temperament scoring pen nested within birth year-season combinations (n = 42). The corresponding model for WBSF included fixed effects of type of cross (an effect based upon the combination of sire and dam breeds, n = 12) and date of shear (n = 16).
Fixed effects that accounted for family structure such as family nested within sire or type of cross may complement or be confounded with genomic information in association analyses; the effect of inclusion of such effects on accuracy of prediction of breeding values is unknown. Genotype frequency differences between families could result in false associations while training markers for prediction of breeding values using Bayesian methods (Lander and Schork, 1994;Marchini et al., 2004;Janss et al., 2012).
Previous cattle studies typically include polygenic effects to account for family (e.g., Michal et al., 2006;Alexander et al., 2009). In addition, studies with the current cattle population (Riley et al., 2013;Hulsman Hanna et al., 2014) utilized fixed effects to account for family structure in association analyses, but as the effect on genomic predictions of including those effects is unknown, analyses in this study were evaluated with or without significant family structure fixed effects. A traditional animal model that accounted for family in the conventional manner using pedigree was also fitted as a control.
For both traits, only those animals with phenotypes, genotypes and all relevant fixed effects required in the model were used in the study (n = 769 or 387 for TEMP or WBSF, respectively). Correlations between TEMP and WBSF phenotype records were previously calculated using either Pearson's or Spearman's Rank correlation coefficients in Hulsman Hanna et al. (2014), and TEMP and WBSF were not found to be correlated in either analysis (p > 0.05).

Bayesian methods employed
BayesB (Meuwissen et al., 2001) and BayesC (Kizilkaya et al., 2010;Habier et al., 2011) methods were utilized in this study as in Hulsman Hanna et al. (2014). Both methods require that p, a proportion of markers not contributing to the trait of interest, is known. The best fitting p (designated asp) was determined in Hulsman Hanna et al. (2014) as the value that accounted for the majority of genetic variation with the fewest number of markers. For TEMP and WBSF,p was determined to be 0.997 and 0.995, respectively. An additional analysis was run using BayesC procedures with p = 0 (i.e., all markers were included in the model). Analyses were performed using GenSel software (see Internet References section) to generate DGV for animals included in the training group by summing across posterior means of random marker effects multiplied by the centered number of copies of the B allele for each individual. All animals with data available were used in the training population for comparison purposes, where cross-validation procedures were conducted following this study and will be reported separately. For each trait, comparisons were made between DGV and estimated breeding values (EBV) produced using a traditional animal model (i.e., traditional pedigree-based Best Linear Unbiased Predictions). Estimated breeding values from an animal model were generated using ASReml software (see Internet References section), fitting the same fixed effects as used in genomic analyses, but with random animal effects with variancecovariance determined from pedigree information, and random uncorrelated residuals.

Comparison criteria for breeding value prediction
Pearson correlation coefficients (r) between true breeding values (BV) and direct genomic values (DGV; 632 Hulsman Hanna et al. i.e., BV based on genomic prediction) can be used in simulations to assess accuracy (Meuwissen et al., 2001). Simple correlation coefficients between de-regressed EBV and DGV can be used to assess accuracy with field data, but that correlation underestimates true accuracy (Saatchi et al., 2012). Accuracy of EBV (or DGV) predicted in this study was calculated following the standardization method proposed by Saatchi et al. (2011) to adjust for underestimation, and takes the form of: where s p,EBV is the covariance of phenotype or trait with either EBV or DGV, s $ g is the standard deviation of the additive genetic effects estimated from the training population reported as the posterior mean or variance estimate (i.e., the mean estimate over the iterations or chains run) of the analysis conducted, and s g is the standard deviation of the additive genetic effects from the population (i.e., animals with phenotypic and genotypic information), which is the expected additive genetic effects of that population for the given trait and is calculated as: where h 2 is the heritability of the trait calculated during the analyses (e.g., the genomic heritability for Bayesian methods and based on variance components for the animal model), and s p 2 is the phenotypic variance calculated from animals with genotypic data. The simple linear regression coefficient ( $ b y,x ), which can be used to give an indication of biasedness (Saatchi et al., 2013) was calculated. The simple linear regression coefficient ( $ b y,x ) is the regression of phenotype on EBV (or DGV) and takes the form of: where s EBV,p is the covariance of EBV (or DGV) and phenotype and s EBV 2 is the variance of the EBV (or DGV) calculated from the DGV or EBV produced in that respective analysis using the sample variance function in Microsoft Excel 2013. This regression parameter should be 1 if unbiased (Saatchi et al., 2011(Saatchi et al., , 2013. The relative ranking of EBV (or DGV) is critical in determining candidates that will be selected and this was assessed in two ways: 1) using Spearman Rank correlation coefficients and 2) coarser assessments of rank changes calculated as the percentage and number of individuals whose EBV (or DGV) were in different quartiles from distinct analyses.

Breeding value prediction
The estimates of heritability for WBSF and TEMP using an animal model (single trait analysis, additive genetic component only) were 0.055 ± 0.091 and 0.35 ± 0.115, respectively. Originally, the type of cross effect was fitted for WBSF analyses using the animal model. This resulted in a heritability estimate of zero, however, and was subsequently excluded from runs using the animal model. Posterior means of genomic heritability estimated through Bayesian methods were higher for WSBF, but lower for TEMP (Table 1).
EBV from animal models for WBSF and TEMP were compared to DGV from three Bayesian methods with and without inclusion of fixed effects to account for significant family structure (Tables 2-4). Numerically, BayesC models withp had the highest accuracy for either trait, regardless whether fixed effects for family structure were included or not ( Table 1). The simple linear regression coefficient ( $ b y,x ) of DGV on phenotype was closest to 1 for BayesB for both traits and either statistical model, indicating less bias for that method. All analyses had significant Spearman rank correlation coefficients (Table 4) for trait breeding values predicted from the fitted models. Inclusion of fixed effects to account for family structure resulted in lower Spearman Rank correlation coefficients for DGV from any Bayesian model compared with EBV from a traditional animal model excluding such fixed effects. Spearman Rank correlation coefficients for EBV with DGV were always lower than those among DGV from the fitted Bayesian models. In general, correlation coefficients for DGV estimated from models that included fixed effects for family structure were lower, indicating substantial influence of that effect on the relative ranking of animals.
No more than 2% of individuals had EBV in different analyses that were more than 2 quartiles apart for WBSF or TEMP (Tables 2 and 3). When fixed effects accounting for family structure were included in the model, WBSF and TEMP were similar in the percentage of animals that changed n quartiles, although that percentage for WBSF was typically numerically lower than TEMP. When the statistical model did not include fixed effects for family structure, TEMP had a larger reduction in the number of animals (and, therefore, lower percentage) that changed n quartiles. When the statistical model included fixed effects accounting for family structure, between Bayesian analyses more than 71% of the DGV remained in the same quartile for either trait, but less than 55% were in the same quartile when comparing DGV from any Bayesian analysis to EBV from the animal model for that respective trait (Tables 2  and 3). Exclusion of fixed effects for family resulted in more stable (relative to quartile ranks) across-analysis estimation, where over 80% were estimated in the same quar-tile with the Bayesian analyses, and approximately 71% for any Bayesian estimation with those from animal models (Tables 2 and 3).

Discussion
The heritability estimate for WBSF is low in comparison to most values reported for Bos taurus cattle with electrically-stimulated carcasses, which have ranged from 0.11 to 0.53 depending on breed or cross and number of measured animals (Shackelford et al., 1994;reviewed by Minick et al., 2004). Although the estimate of heritability for WBSF did not differ from zero in the present study, others have reported low values using Bos indicus cattle due to an apparently low additive genetic component for this trait in these types of cattle (Riley et al., 2003;Smith et al., 2007). The effect of electrical stimulation on these carcasses could be reducing variation in WBSF resulting in lower estimates of heritability. Reported estimates of heritability for temperament measured as flight speed have been reported as high as 0.40 (Burrow, 2001) and the estimate of 0.35 using an animal model from the present study is consistent with those estimates.
Traditional animal models account for family structure using pedigree information (Henderson, 1973(Henderson, , 1975, but with the inclusion of genomic information, differences in genotype frequency due to family structure could lead to false associations, and therefore impact the accuracy of DGV. Results from the present study indicated that inclusion of fixed effects to account for differences between families decreased accuracy and increased bias associated with prediction of DGV in this crossbred population. When fixed effects for family structure were not included, Bayesian methods were more accurate, had less bias, and had less re-ranking between analyses compared to an animal model. As all animals with data available were included in the training population, accuracy was expected to be moderately high, even though the sample sizes were lower than 1,000 individuals (the typical standard minimum for running genetic evaluations) because the information used to estimate effects were from those animals. Without including family structure fixed effects, accuracy for Bayesian and animal model analyses were higher than may typically be expected given the sample sizes, but training population accuracies are not typically reported and, therefore, comparisons are not available. Although all analyses showed 634 Hulsman Hanna et al. s 2 , which provides an estimate of what the additive genetic variance is expected to be in that population for the given trait. s p 2 = 0.413 or 4.253 for Warner-Bratzler shear force or overall temperament at weaning, respectively, and is the estimate of phenotypic variance from the data. s $ g 2 is the posterior mean or variance estimate of additive genetic variance from the training potulation calculated in that respective analysis. s EBV 2 is the variance of the estimated breeding values estimated using either Bayesian or animal model procedures, wherep = 0.995 or 0.997 for Warner-Bratzler shear force or overall temperament at weaning, respectively, using genomic information in Bayesian models. bias in the estimates, it was reduced when the sample size increased (e.g., TEMP vs. WBSF), which was expected. Habier et al. (2013) investigated genomic BLUP and its ability to capture relationships, model LD and exploit cosegregation (the deviation from independent segregation of alleles if loci are linked). They found that genomic BLUP can exploit LD, cosegregation, and additive-genetic relationships captured by SNP, but suggest that Bayesian methods with t-distributed priors (e.g., BayesB or BayesC in this study) may be more beneficial to account for rapid decay in LD. This would explain why the Bayesian methods utilized in this study did not need fixed effects to account for family structure and further suggested that inclusion of those fixed effects detrimentally impacted the ability of Bayesian methods to utilize genomic information from the current population.
Simulation studies have reported superior accuracy values from Bayesian methods relative to traditional animal models (e.g., Meuwissen et al., 2001;Piyasatian et al., 2006). Habier et al. (2011) concluded that the optimal Bayesian modeling method to use was trait dependent when they compared six Bayesian models. In the case of TEMP and WBSF, there was little difference between Bayesian models, although BayesC procedures using $ p had numerically higher accuracies, but larger bias than BayesB for both traits.
As accuracy and the regression coefficient were calculated using genetic parameter estimates from the respective analyses (see Table 1), it is interesting to note the differences seen when comparing the variance of the breeding values estimated (s EBV 2 ; either DGV or EBV, which are technically estimates of additive genetic effects and should Compare crossbred BV prediction 635 The number of quartiles changed was calculated by first assigning an animal's quartile for any given analysis, then finding the difference of each animal's quartile between the two analyses compared. Percentage was calculated by dividing the number of individuals within that category by the total number of animals (n = 387). 2p = 0.995 3 Family structure fixed effect refers to type of cross (an effect based upon the combination of sire and dam breeds, n = 12) for Warner-Bratzler shear force. be equivalent to s $ g 2 ) to the additive genetic variance calculated through the analysis (s $ g 2 ). When family structure fixed effects were included in Bayesian analyses, the ratio of s EBV 2 to s $ g 2 ranged from 0.246 to 0.289, which was much lower than when excluded (0.544 to 0.561 for Bayesian models and 0.479 for the animal model) of TEMP. This difference was not observed for WBSF, most likely due to sample size.
Furthermore, the estimate (either the posterior mean or variance component) of additive genetic variance (s $ g 2 ) was always lower than what was expected (s g 2 ). Ultimately, the differences found between the three variance parameters of s EBV 2 , s $ g 2 , s g 2 and provide an assessment of the estimation procedure, as s g 2 is the value s EBV 2 and s $ g 2 are expected to be. This is not surprising or novel, as these models do not take into account interactions between or within causative loci for a given trait, thereby creating missing variance (e.g., see Zuk et al. (2012) discussion on phantom heritability). Differences may have been strongly influenced by sample size available for this study. The results from this study show, however, that the variance estimated from the analysis (s $ g 2 ) compared to the expected variance (s g 2 ) are relatively similar, the difference is consistent across the models within a trait, and the use of family structure fixed effects in the model should not be included for prediction of genetic merit.