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 ; VanRaden and Sullivan, 2010; Saatchi , 2012, 2013; Rius-Vilarrasa ). The landmark paper by Meuwissen 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 reference (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 ; Ibáñez-Escriche ; Habier , 2013; Kizilkaya ; Saatchi ; Toosi ).

Relatedness among animals and between breeds (in crossbred or composite populations) can directly impact prediction accuracy (de Roos ; Ibáñez-Escriche ), at least in part due to inconsistent LD patterns among cattle breeds (Gautier ; Bovine HapMap Consortium, 2009). Few studies have investigated the use of genomic prediction models using field or research data from crossbred populations (Hayes ; Snelling ). 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.

Materials and Methods

Cattle, genotypes, and traits

The cattle population, genotypes, and traits used in this study were previously described in Hulsman Hanna . Briefly, the cattle population consisted of Nellore-Angus F2 (embryo transfer and natural service reciprocal crosses), F3 and natural service half-siblings to embryo transfer F2 animals. The population was characterized by 3 distinct cycles (designated as Cycle 1, 2, and 3), where Cycle 1 included embryo transfer Nellore-Angus F2 and natural service half-sibling calves, Cycle 2 included natural service F2 calves, and Cycle 3 included natural service F3 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 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 ; Janss ). Previous cattle studies typically include polygenic effects to account for family (e.g., Michal ; Alexander ). In addition, studies with the current cattle population (Riley ; Hulsman Hanna ) 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 , and TEMP and WBSF were not found to be correlated in either analysis (p > 0.05).

Bayesian methods employed

BayesB (Meuwissen ) and BayesC (Kizilkaya ; Habier ) methods were utilized in this study as in Hulsman Hanna . Both methods require that π, a proportion of markers not contributing to the trait of interest, is known. The best fitting π (designated as π̃) was determined in Hulsman Hanna as the value that accounted for the majority of genetic variation with the fewest number of markers. For TEMP and WBSF, π̃ was determined to be 0.997 and 0.995, respectively. An additional analysis was run using BayesC procedures with π = 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 variance-covariance 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; i.e., BV based on genomic prediction) can be used in simulations to assess accuracy (Meuwissen ). 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 ). Accuracy of EBV (or DGV) predicted in this study was calculated following the standardization method proposed by Saatchi to adjust for underestimation, and takes the form of: where σp,EBV is the covariance of phenotype or trait with either EBV or DGV, σĝ 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 σ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 h2 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 is the phenotypic variance calculated from animals with genotypic data. The simple linear regression coefficient (β̂y,x), which can be used to give an indication of biasedness (Saatchi ) was calculated. The simple linear regression coefficient (β̂y,x) is the regression of phenotype on EBV (or DGV) and takes the form of: where σEBV,p is the covariance of EBV (or DGV) and phenotype and 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 , 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.

Results

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).

Table 1

Genetic parameter estimates (ĥ2, , and ), prediction accuracies (ρ̂gĝ), and linear regression coefficients (β̂y,x) for Warner-Bratzler shear force and overall temperament at weaning for 4 models.

	Family structure fixed effect included2			Family structure fixed effect excluded2
	BayesC (π = 0)	BayesC (π̃)	BayesB (π̃)	BayesC (π = 0)	BayesC (π̃)	BayesB (π̃)	Animal model
Warner-Bratzler shear force
ĥ2	0.155	0.136	0.234	0.138	0.148	0.236	0.055
σg2	0.064	0.056	0.097	0.057	0.061	0.097	0.023
σg^2	0.054	0.047	0.084	0.049	0.052	0.085	0.020
σEBV2	0.011	0.009	0.023	0.010	0.012	0.026	0.002
ρ̂gĝ	0.787	0.795	0.775	0.841	0.843	0.823	0.717
β̂y,x	4.157	4.610	3.030	4.359	4.115	2.886	9.326
Overall temperament at weaning
ĥ2	0.210	0.206	0.227	0.278	0.274	0.271	0.350
σg2	0.893	0.876	0.965	1.181	1.164	1.154	1.487
σg^2	0.740	0.721	0.802	1.025	1.010	1.004	1.313
σEBV2	0.182	0.193	0.232	0.558	0.562	0.563	0.629
ρ̂gĝ	0.671	0.702	0.698	0.918	0.923	0.918	0.862
β̂y,x	2.996	2.885	2.646	1.810	1.782	1.755	1.916

ĥ2 is the estimate of heritability, where ĥ2 for the Bayesian models is the posterior mean of the genomic heritability, is the additive genetic variance calculated as , which provides an estimate of what the additive genetic variance is expected to be in that population for the given trait. or 4.253 for Warner-Bratzler shear force or overall temperament at weaning, respectively, and is the estimate of phenotypic variance from the data. is the posterior mean or variance estimate of additive genetic variance from the training potulation calculated in that respective analysis. is the variance of the estimated breeding values estimated using either Bayesian or animal model procedures, where π̃ = 0.995 or 0.997 for Warner-Bratzler shear force or overall temperament at weaning, respectively, using genomic information in Bayesian models.

Family structure fixed effect refers to family nested within sire or type of cross (an effect based upon the combination of sire and dam breeds, n = 12) for overall temperament at weaning or Warner-Bratzler shear force, respectively.

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 with π̃ 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 (β̂ 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.

Table 2

Comparison of the number and percentage of individuals with estimated breeding values for Warner-Bratzler shear force that changed n quartiles between any two analyses.

Models compared2	Number of individuals that changed n quartiles1
	With family structure fixed effect3			Without family structure fixed effect3
	1	2	3	1	2	3
BayesC (π = 0) vs. BayesC (π̃)	28 (7.24%)	0	0	28 (7.24%)	0	0
BayesC (π = 0) vs. BayesB (π̃)	46 (11.89%)	0	0	66 (17.05%)	0	0
BayesC (π = 0) vs. Animal model	140 (36.18%)	40 (10.34%)	4 (1.03%)	136 (35.14%)	22 (5.68%)	2 (0.52%)
BayesC (π̃) vs. BayesB (π̃)	36 (9.30%)	0	0	44 (11.37%)	0	0
BayesC (π̃) vs. Animal model	141 (36.43%)	36 (9.30%)	5 (1.29%)	141 (36.43%)	22 (5.68%)	3 (0.78%)
BayesB (π̃) vs. Animal model	144 (37.21%)	38 (9.82%)	2 (0.52%)	139 (35.92%)	22 (5.68%)	3 (0.78%)

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).

π̃ = 0.995

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.

Table 4

Spearman rank correlations (p < 0.001) of breeding values predicted from different analyses for overall temperament at weaning and Warner-Bratzler shear force

	BayesC (π = 0)	BayesC (π̃)	BayesB (π̃)	Animal model
With family structure fixed effect2
BayesC (π = 0)		0.964	0.942	0.449
BayesC (π̃)	0.996		0.996	0.460
BayesB (π̃)	0.984	0.991		0.475
Animal model	0.652	0.655	0.658
Without family structure fixed effect2
BayesC (π = 0)		0.988	0.982	0.911
BayesC (π̃)	0.995		0.999	0.899
BayesB (π̃)	0.982	0.993		0.895
Animal model	0.810	0.805	0.796

Analyses listed are the Bayesian or traditional mixed model procedure used to estimate breeding values, where π̃ = 0.995 or 0.997 for Warner-Bratzler shear force (bottom diagonal) or overall temperament at weaning (top diagonal), respectively.

Family structure fixed effect refers to family nested within sire or type of cross (an effect based upon the combination of sire and dam breeds, n = 12) for overall temperament at weaning or Warner-Bratzler shear force, respectively.

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.

Table 3

Comparison of the number and percentage of individuals with estimated breeding values for overall temperament at weaning that changed n quartiles between any two analyses.

Models compared2	Number of individuals that changed n quartiles1
	With family structure fixed effect3			Without family structure fixed effect3
	1	2	3	1	2	3
BayesC (π = 0) vs. BayesC (π̃)	174 (22.63%)	1 (0.13%)	0	88 (11.44%)	0	0
BayesC (π = 0) vs. BayesB (π̃)	218 (28.35%)	4 (0.52%)	0	108 (14.04%)	0	0
BayesC (π = 0) vs. Animal model	309 (40.18%)	83 (10.79%)	7 (0.91%)	220 (28.61%)	7 (0.91%)	0
BayesC (π̃) vs. BayesB (π̃)	56 (7.28%)	0	0	30 (3.90%)	0	0
BayesC (π̃) vs. Animal model	309 (40.18%)	78 (10.14%)	7 (0.91%)	232 (30.17%)	10 (1.30%)	0
BayesB (π̃) vs. Animal model	317 (41.22%)	68 (8.84%)	9 (1.17%)	240 (31.21%)	12 (1.56%)	0

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 = 769).

π̃ = 0.997

Family structure fixed effect refers to family nested within sire for overall temperament at weaning.

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 quartile 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 ; reviewed by Minick ). 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 ; Smith ). 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, 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 bias in the estimates, it was reduced when the sample size increased (e.g., TEMP vs. WBSF), which was expected.

Habier 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 ; Piyasatian ). Habier 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 π̂ 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 ( ; either DGV or EBV, which are technically estimates of additive genetic effects and should be equivalent to ) to the additive genetic variance calculated through the analysis ( ). When family structure fixed effects were included in Bayesian analyses, the ratio of to 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 ( ) was always lower than what was expected ( ). Ultimately, the differences found between the three variance parameters of , , and provide an assessment of the estimation procedure, as is the value and 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 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 ( ) compared to the expected variance ( ) 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.