The Accuracy and Bias of Single-Step Genomic Prediction for Populations Under Selection.

In single-step analyses, missing genotypes are explicitly or implicitly imputed, and this requires centering the observed genotypes using the means of the unselected founders. If genotypes are only available for selected individuals, centering on the unselected founder mean is not straightforward. Here, computer simulation is used to study an alternative analysis that does not require centering genotypes but fits the mean [Formula: see text] of unselected individuals as a fixed effect. Starting with observed diplotypes from 721 cattle, a five-generation population was simulated with sire selection to produce 40,000 individuals with phenotypes, of which the 1000 sires had genotypes. The next generation of 8000 genotyped individuals was used for validation. Evaluations were undertaken with (J) or without (N) [Formula: see text] when marker covariates were not centered; and with (JC) or without (C) [Formula: see text] when all observed and imputed marker covariates were centered. Centering did not influence accuracy of genomic prediction, but fitting [Formula: see text] did. Accuracies were improved when the panel comprised only quantitative trait loci (QTL); models JC and J had accuracies of 99.4%, whereas models C and N had accuracies of 90.2%. When only markers were in the panel, the 4 models had accuracies of 80.4%. In panels that included QTL, fitting [Formula: see text] in the model improved accuracy, but had little impact when the panel contained only markers. In populations undergoing selection, fitting [Formula: see text] in the model is recommended to avoid bias and reduction in prediction accuracy due to selection.

In pedigree-based analyses, the expected value of breeding values is zero. In order to achieve similar properties in whole-genome analyses, marker genotype covariates are often transformed. When all individuals are genotyped, it has been shown that inference on genotype effects does not depend on how the covariates are transformed (Strandén and Christensen 2011). However, when data includes genotyped and nongenotyped individuals, inference on marker effects from single-step analyses may depend on how the covariates are transformed (Fernando ). In single-step analyses using marker effects models, the breeding values of nongenotyped individuals are partitioned into components representing the prediction of nongenotyped individuals conditional on their genotyped relatives and an independent deviation (Fernando ). The prediction of nongenotyped individuals conditional on their genotyped relatives is done based on best linear prediction, which requires the first moments to be known without error. This is straightforward if the mean of the genomic breeding value is zero in the absence of selection. Centering the observed genotype covariates using what their means would be in the absence of selection would result in genomic breeding values with null means. However, such genotype covariate means are typically unavailable. Fernando proposed a solution for the marker effects model that involves fitting an additional fixed covariate that estimates the mean of the linear component of the genotypic value, which is denoted by in Equation 1 below, in a population where selection is absent. Using this approach, even when there is selection, the selection process can be ignored, because the analysis is conditional on the data used for selection (Goffinet 1983; Gianola and Fernando 1986; Im ; Sorensen ). In Markov chain Monte Carlo (MCMC) analyses, centering results in better mixing (Strandén and Christensen 2011), reducing the number of iterations required to obtain converged genomic predictions. In practice, centering the entire matrix of genotype covariates, including the observed and imputed genotypes, using their sample means can be done in addition to the Fernando approach of fitting This type of centering of the entire genotype matrix does not affect inference on marker effects.

The same issue with centering of the observed genotype covariates that we described above for the marker effects model is also implicit for the single-step breeding value model (single-step GBLUP), and a similar solution was proposed by Vitezica . In their proposed solution, the observed genotype covariates are centered using their means; in addition, the genomic covariance matrix is corrected for the change in the mean breeding value of the genotyped individuals (Vitezica ). It was shown in that paper that this is equivalent to fitting the change in breeding value due to selection as a random effect.

Here, we use simulated data to compare the accuracy and bias in genomic prediction applied to populations under selection with and without centering the entire matrix of genotype covariates, and with and without fitting as a fixed effect. Further, we show that when the observed genotype covariates are centered using means calculated from selected individuals rather than means from all individuals, the meaning of changes from the mean of unselected individuals to become the mean breeding value in selected individuals, as claimed by Vitezica .

Methods

Theory

To simplify the presentation of the genetic model without loss of generality, we will assume that the unconditional expectation of the phenotypic value for all individuals is the same. Let denote the row vector of genotypes for individual i, which is often coded as −1,0,1. Then, under additive gene action, the genotypic value, which is the expected phenotypic value of an individual with genotypes , can be written aswhere β is the value of when and is the vector of substitution effects. The scalar β and the vector are constants, but will be a random variable because of randomness in owing to the randomness in the genotypes for a randomly sampled individual. Note that the expected value of the linear component of the genotypic value in Equation 1 is where which may not be equal to zero. Thus, it is customary to write the model for the genotypic value, as can be derived from Equation 1, as follows:where is a constant, representing and is a random variable that has null expectation, which is the breeding value predicted in a pedigree-based BLUP (PBLUP) evaluation. When genotypes are observed and used in a genomic analysis, they may be transformed or coded by subtracting their expectations, from the observed values, In both Equations 1 and 2, is the same substitution effect for locus j. The intercepts in these models, however, are different. In Equation 1, the intercept is β, and it is the value of when In Equation 2, on the other hand, the intercept is and it is the value of when More generally, is not known, so genotypes are coded by subtracting a different vector from the observed genotypes as Still, is the substitution effect for locus j, but the intercept will change to become which is the value of when Thus, as more rigorously shown in (Strandén and Christensen 2011), inference about does not depend on how the genotypes are coded. A simpler but rigorous proof is given in the Appendix of this paper.

In single-step analyses, where some individuals are not genotyped, the missing genotypes are imputed either implicitly (Legarra ) or explicitly (Fernando ) using best linear prediction. Let denote the matrix of genotypes for individuals that were genotyped. Then, the genotypes of the individuals with missing genotypes are imputed aswhere is the matrix of pedigree-based additive relationships between the nongenotyped and genotyped individual, and is the matrix of additive relationships among genotyped individuals. Now, the model for the genotypic values, when genotypes are coded as in Equation 2, becomeswhere is that part of that cannot be imputed from knowledge of the breeding values of genotyped relatives. In practice, the true value of is not known, and data for its estimation may not be available. Rearranging these equations in terms of the uncentered rather than the centered matrix of genotype covariates results inand substituting as previously defined, and letting results inwhich suggests that could be treated as an unknown constant and estimated as a fixed effect from the data (Fernando ). The covariate vector for is denoted by for nongenotyped individuals, and by for genotyped individuals. So, Equation 3 becomeswhich can be combined aswhere and the 0 matrix in the last term of Equation 4 is required because does not appear in the model for genotyped individuals.

When the vector represents the substitution effects of a large number of loci containing positive and negative effects, will tend to have a value close to zero. Accordingly, we have simulated some scenarios with positive so that the entire vector is positive to exacerbate the impact of Nevertheless, when marker rather than causal alleles are fitted in the model, the sign of the substitution effects depends on the phase relationship between marker and causal allele, which is equally likely to be positive or negative.

Even if in a population undergoing selection, it is expected that in nonfounders. Suppose is the mean of the observed genotype covariates in such a population undergoing selection, and these means are used to center the matrix of observed genotypes. Then, the model for the genotypic values can be written in terms of the matrix of centered covariates asand using for the covariate corresponding to Equation 5 can be written aswhich can be combined asNote that the regression coefficients for in Equation 4 and in Equation 6, must be equal. This implies that Similarly, the intercepts of these two models must be equal too, and this implies

Simulations

Phenotypic and genotypic data were simulated using XSim (Cheng ), based on haplotypes from 10 genomic regions of 721 US Hereford beef cattle that were genotyped with the Illumina 770K BovineHD BeadChip, and reported in terms of the number of copies of the A allele at each locus. The selected regions came from choosing every 10th single nucleotide polymorphism (SNP), starting at SNP 5001, to get a low-density panel of 200 SNPs from each of chromosomes 1–10 (BTA1–BTA10), after eliminating SNPs with MAF < 0.01. These 2000 SNPs represent 10 0.1 M chromosomes. Average linkage disequilibrium (LD) between adjacent SNPs was 0.300. 50 SNPs from each chromosome were randomly chosen to represent QTL. The QTL effects were sampled from a normal distribution with mean = 0.2 and multiplied by the number of copies of the A allele to produce the true breeding value (TBV). The TBVs were added to a normally distributed residual term scaled by the sample variance of the TBV to simulate a trait with a heritability of 0.5. The first 20 SNPs from each of the 10 chromosomes were also used to simulate a smaller panel, with 5 QTL and 15 markers per chromosome. The average LD between adjacent SNPs was 0.289. TBV were simulated in the same way as for the 2000 SNP scenario, then scaled to simulate traits with heritabilities () of 0.1, 0.3, or 0.5. An additional scenario with = 0 was used to simulate TBV for a trait with heritability 0.5.

Half the observed diplotypes from US Hereford cattle were assigned to represent males and the remainder to represent females. Those 360 males and 361 females were sampled in pairs, with replacement, to produce 4000 male and 4000 female offspring representing generation G4. There were no mutations. Four more nonoverlapping generations of random mating were carried out, with one male and one female offspring per dam mated to randomly chosen sires to produce the founder population (G0).

The G1 generation was produced by mass phenotypic selection of the top 200 G0 males, and this was repeated for five more generations. Each female was randomly mated twice to selected males to produce one offspring of each sex each generation. Across nonoverlapping generations G0–G5, a total of 48,000 individuals with phenotypes, genotypes, and TBVs were simulated for each scenario.

The training data included phenotypes from all individuals in G0–G4 (n = 40,000), and genotypes from all 1000 sires and all 8000 G5 animals. Fixed loci, if any, were filtered from the panel before genomic prediction analyses. The genetic and residual variances used in genomic prediction were the sample variance of the TBV in G0 and the corresponding residual variance used to define the desired heritability in the founder population.

Models

Five statistical models were compared for differences in accuracy and bias of prediction. These included models with or without , and with or without centering of observed and imputed marker covariates, and a model that used pedigree relationships but not marker covariates.

Mixed linear model:

Accuracy of PBLUP was quantified using the correlation of TBV and estimated breeding values (EBV), where TBV was as simulated and EBV were obtained by fitting the mixed linear model (Henderson 1973, 1984):where is a vector of phenotypic observations, 1 is a vector of 1s, is a general mean, is a vector of random additive genetic effects, is a vector of random residual effects, and is a known incidence matrix relating observations to In this model, so that Further, with being the numerator relationship matrix, and so that

Single-step Bayesian regression model:

All genomic EBVs (GEBV) were obtained by single-step Bayesian regression (Fernando ), using BayesC priors for marker effects with π = 0, which gives predictions that are identical to those from single-step GBLUP. The model was implemented in Julia (http://julialang.org) based on the SSBR package (http://QTL.rocks) to construct an MCMC chain of 50,000 samples. Individuals were separated into two groups designated with subscripts g or n according to whether or not simulated genotypes were assumed to be observed or missing. The single-step Bayesian regression model including a covariate J for (model J) was:where and are vectors of phenotypes for nongenotyped and genotyped individuals, 1 is a vector of 1s, μ is a general mean, is the expected value of the linear component of the genotypic value if selection was absent, is a vector of random substitution effects of markers, is a vector of imputation residuals, and are incidence matrices relating the breeding values of nongenotyped and genotyped individuals to their phenotypes, which is defined for genotyped individuals, is a vector of −1s, which is defined for nongenotyped individuals, is a vector computed as is the matrix of imputed marker covariates, is the matrix of observed marker covariates, and and are vectors of random residual effects for nongenotyped and genotyped individuals. This model can be represented aswhere U = and is a vector of random residual effects.

This model was used for four variants of the single-step Bayesian regression analysis, depending on whether or not the covariate J corresponding to the mean was in the model, and whether or not the columns in the marker covariate matrix were centered using means of the observed and imputed genotype covariates. The analyses with J or without J are denoted as J or N when covariates were not centered, and as JC or C when the entire matrix of imputed and observed genotype covariates were centered, respectively (Table 1).

Table 1

Four combinations of the single-step Bayesian regression analyses

Models	Marker Covariates
	Centereda	Not Centeredb
With J and μ_g	JC	J
Without J	C	N

Centered, e.g., genotype values represented as −1, 0, 1 when the uncentered genotype covariate with values 0, 1, 2 has mean 1.

Not centered, e.g., genotype values represented as the number of copies of the A allele.

Accuracy of genomic prediction was quantified using the correlation of TBV and GEBV (), where GEBVs were obtained from each of the four analyses described above. Bias of genomic prediction was quantified using the deviation from unity of the coefficient of regression of TBV on GEBV (). In models JC and J, GEBVs are obtained using Equation 24 in (Fernando ):where is the GEBV, is the best linear unbiased estimate of the mean of breeding values, is the BLUP of the vector of random substitution effects of all markers, and is the BLUP of the imputation residual.

In model J, the matrix contains the uncentered number of copies of the A allele at each locus, and the uncentered version is used to impute In model JC, the entire matrix of imputed, and observed, genotype covariates is centered. In models C and N, the GEBVs are computed in a corresponding manner, except that the covariate and its coefficient are not included in the model.

The four analyses J, N, JC, and C were all applied to three different genotype panels, comprising the causal QTL plus markers, just the causal QTL, or just the markers. All 12 combinations of four analyses and three genotype panels were applied to data simulated with 200 loci comprising 50 QTL whose effects were sampled from a normal distribution with to construct phenotypes with The four analyses JC, J, C, and N were repeated using only genotype panels comprising QTL plus markers for three other scenarios: 200 loci, 200 loci, and 200 loci,

Fernando had observed that including in the model was necessary only when and the number of observations exceeds the number of markers. So, to examine the impact of the scenario with 200 loci and was repeated with and to examine the impact of the number of loci, the scenario with and was repeated with 2000 loci.

Every scenario was replicated 10 times with each replicate having been constructed starting from the sampling of G5 which represented simulated offspring from the haplotypes of real animals. Every phenotypic dataset was also fitted using the PBLUP model. All reported correlations and regression coefficients are the means of 10 replicates. These are presented along with the SEs of those means.

In single-step GBLUP (Legarra ; Aguilar ; Christensen and Lund 2010), the missing genotypes are not explicitly imputed, and only the observed genotype covariates are centered using their means. We argued earlier that when the number of loci is large, will be close to zero, especially when the panel consists of only marker loci. However, in a population undergoing selection, which is the mean breeding value of genotyped individuals, is not expected to be zero even when the panel only includes marker loci and So, in addition to the above, single-step Bayesian regression analyses, using the model given in Equation 7 with and without J (models JC* and C*), were applied to a marker panel with 200 loci, and when the matrix of observed genotype covariates were centered using their means as: where and is the number of individuals with genotypes. Recall that even when the means of the observed genotypes are used for centering, the model for the genotypic values can be written in terms of the matrix of uncentered covariates as shown by Equation 6. We will compare the estimates of in this model with those of from Equation 4 where the means of observed genotype covariates are not used for centering.

Data availability

The genotypes representing G0 from each replicate and scenario are available at: https://figshare.com/s/d7798b811a9a6a4172fc. These genotypes and the methodology described previously are sufficient to reproduce the simulations used in this study.

Results and Discussion

Effect of fitting a genotypic mean and centering of observed and imputed marker covariates

Accuracy:

The accuracies of genomic prediction as assessed by validation in G5 after training using G0–G4 for a trait with 50 QTL whose effects were sampled from a normal distribution with = 0.2 and = 0.5 are in Table 2. The accuracy of PBLUP in predicting breeding values for individuals without phenotypes was 41.5%, accounting for <20% of genetic variance. All analyses using genotypes resulted in more accurate predictions than using pedigree alone. Centering the entire matrix of imputed and observed genotype covariates had no effect on the accuracy of the genomic analyses regardless of the nature of the marker panel. This is because, as shown in the Appendix, this type of centering is equivalent to adding and subtracting from the model equation, and this has no effect on the mixed model solutions for

Table 2

Correlations (%, ) between TBV and (G)EBV for alternative analyses

Genotype Dataa	Analyses
	JCb	Jc	C	N	PBLUP
50 QTL + 150 markers	97.59 ± 0.00	97.63 ± 0.00	96.32 ± 0.00	96.31 ± 0.00	—
50 QTL only	99.44 ± 0.00	99.45 ± 0.00	90.20 ± 0.05	90.18 ± 0.05	—
150 markers only	80.47 ± 0.03	80.47 ± 0.03	80.44 ± 0.03	80.44 ± 0.03	—
No genotypes	—	—	—	—	41.56 ± 0.01

Average correlation between true breeding value (TBV) and (genomic) estimated breeding values from 10 replications validated in Generation 5, comprising 8,000 individuals with genotypes but no phenotypes. The true QTL effects were sampled from a Normal distribution with mean = 0.2 and scaled to simulate a trait with a heritability 0.5.

J: includes a covariate for μ in the model, C: entire matrix of imputed and observed genotype covariates centered, JC: both J and C, N: neither J or C, and PBLUP: pedigree-based BLUP.

The analyses were based on fitting covariates for only 50 QTL, only 150 markers, or both 50 QTL and 150 markers.

In this study, selection resulted in the successive advances in mean TBV from G0 to G5 being 10.14, 10.63, 11.18, 11.68, 12.15, and 12.57. The mean genotypic value was not zero in G0 because QTL genotypes were not centered and the mean QTL effect was = 0.2. Since the QTL effects do not change with selection, the advance in TBV reflects changes in the frequencies of the favorable alleles of the 50 QTL. So centering using the allele frequency means of the genotyped sires in G1–G4 and all individuals in G5 does not closely approximate the centering that would have occurred if the allele frequency means had been obtained from the unselected population. By contrast, fitting in the model estimates the relevant mean from the data.

In panels that included causal variants (QTL), fitting in the model substantially improved the accuracy to being near perfect. This is not surprising, given that there were only 50 QTL, the heritability was 0.5, and there were 40,000 phenotyped ancestors, including 200 genotyped sires per generation in the training. However, in the panel that contained only markers with no causal variants, fitting in the model had little impact. This is because in the population that was simulated here the founders were not selected, and thus was close to zero for the panel with only markers. This would not be the case in a population where even the founders are descendants of selected parents. Thus, for traits that have been subject to selection, fitting in the model is expected to improve accuracy.

Using one replicate as an example, for the panel including both QTL and markers, the estimate of μ was ∼9.90 for the analyses using J and N. The estimate of μ was 4.98 for the analysis using J. For the genotyped individuals, the covariate values in are all −1, so + is a vector of values equal to 9.90−4.98 = 4.92. For nongenotyped individuals, the covariate values in can vary widely, but many are close to −1 while others are close to 0. This means that will include values that range from 9.90 to 4.92, accounting for variation in accuracy of imputation. When is not included in the model, these effects are ignored, which can reduce the accuracy of predicting nongenotyped individuals. Failing to account for these effects will propagate errors in and the latter impacting the accuracy of predicting genotyped individuals. Collectively, these errors reduced accuracy from 97 to 96% for the panel including QTL and markers, and from 99 to 90% for the panel including only QTL. However, when the panel comprised only markers, the estimates will include both positive and negative values because the phase of markers and QTL are equally likely to take either sign, in which case will be close to zero as confirmed in the above mentioned replicate where the estimate was −1.18.

Here, results from the analysis with only QTL on the panel are used to show that is the mean of in the founder population and not the mean of the breeding value of the genotyped individuals. Recall that the QTL model was used with an intercept of to simulate the data. Thus, when only QTL are on the panel, the true value of β is zero. In analysis J because both and are estimates of and could be pooled, which for the replicate above would be In that replicate, the actual mean of in G0 was 10.16, which was estimated in the analysis to be 9.05. On the other hand, the mean of the breeding value in the 9000 genotyped individuals was 2.4, which is clearly not near the pooled estimate of 9.05 for These genotyped individuals included 1000 selected sires, of which 200 were genotyped in each generation from G0 to G4, and 8000 offspring from G5. The mean values of for the selected sires were 1.00, 1.59, 2.08, 2.50, and 2.88, respectively, for G0 through G4, and 2.45 for the offspring in G5. It is apparent that the parameter corresponding to the covariate is the mean of in the founder population and not the mean breeding value in the selected individuals. In analysis JC with the covariates centered, the intercept β is the value of when which is the case when The estimate was 16.11 in this analysis, but remained about the same value, namely 7.80. This shows that has the same interpretation whether the entire matrix of observed and imputed genotypes is centered or not. In neither case does it represent the mean breeding value of selected individuals.

Accuracy of PBLUP increased with heritability, as expected (Table 3). Further, genomic predictions using panels that include causal mutations were near perfect when was included in the model. These high accuracies are a reflection of these phenotypes being influenced by only 50 QTL and there being a large training dataset. Accuracy was reduced when was not fitted in the model. There was no advantage in terms of accuracy to centering the covariates, but MCMC mixing may have been improved, although this was not investigated.

Table 3

Correlations (%, ) between TBV and (G)EBV for alternative analyses for different heritabilities

Heritabilitiesa	Analysesb
	JC	J	C	N	PBLUP
h^2=0.1	93.68 ± 0.01	93.71 ± 0.01	92.42 ± 0.01	92.42 ± 0.01	31.33 ± 0.02
h^2=0.3	96.79 ± 0.01	96.81 ± 0.00	95.19 ± 0.01	95.19 ± 0.01	37.07 ± 0.02
h^2=0.5	97.59 ± 0.00	97.63 ± 0.00	96.32 ± 0.00	96.31 ± 0.00	41.56 ± 0.01

Average correlation between true breeding value (TBV) and (genomic) estimated breeding values from 10 replications validated in Generation 5, comprising 8000 individuals with genotypes but no phenotypes. The true quantitative trait loci (QTL) effects were sampled from a normal distribution with mean = 0.2 and scaled to simulate a trait with heritabilities 0.1, 0.3 or 0.5.

J: includes a covariate for μ in the model, C: entire matrix of imputed and observed genotype covariates centered, JC: both J and C, N: neither J or C, and PBLUP: pedigree-based BLUP. Covariates were fitted for both 50 QTL and 150 markers.

Bias:

Table 4 shows the regression coefficients of TBV on (G)EBV for and = 0.2, the same scenarios represented in Table 2. The regression coefficients of TBV on GEBV for each scenario were close to 1 with very low SE, which indicates that the genomic predictions exhibited almost no bias. The differences in regression coefficients between analyses were very small, but the marker panel comprising only markers was biased upward, whereas the marker panels that included causal mutations were biased slightly downward.

Table 4

Regression coefficients () of TBV on (G)EBV

Genotype Dataa	Analysesb
	JCc	J	C	N	PBLUP
50 QTL + 150 markers	1.06 ± 0.02	1.06 ± 0.02	1.06 ± 0.02	1.06 ± 0.02	—
50 QTL only	1.05 ± 0.02	1.05 ± 0.02	1.12 ± 0.04	1.12 ± 0.04	—
150 markers only	0.95 ± 0.03	0.95 ± 0.03	0.95 ± 0.03	0.95 ± 0.03	—
No genotypes	—	—	—	—	0.95 ± 0.04

Average Regression coefficients of true breeding value (TBV) on (genomic) estimated breeding values from 10 replications validated in Generation 5, comprising 8,000 individuals with genotypes but no phenotypes. The true QTL effects were sampled from a Normal distribution with mean μ= 0.2 and scaled to simulate a trait with a trait with a heritability 0.5.

The analyses were based on fitting covariates for only 50 QTL, only 150 markers, or both 50 QTL and 150 markers.

J: includes a covariate for μ in the model, C: entire matrix of imputed and observed genotype covariates centered, JC: both J and C, N: neither J or C, and PBLUP: pedigree-based BLUP.

Effect of mean QTL effect ( vs. )

We had hypothesized that the impact of omitting from the model would be greatest when departs significantly from 0, which is more likely to occur when departs from 0. For that reason, our base simulation used = 0.2. Results are shown in Table 5 for the panel including QTL and markers with = 0.5 for = 0.2 compared with = 0. These results confirmed that the benefit of fitting was greatest when = 0.2, but there was still an advantage to fitting when = 0. That advantage is likely to erode as the number of QTL increases. Changing the mean QTL effect had no impact on bias, except for a slight influence on PBLUP.

Table 5

Accuracy and bias of genomic prediction () for alternative QTL distributions and analyses

Substitutiona Effects	Analysesb
	JCc	Jd	C	N	PBLUP
Correlations (%)
μ_α=0	97.91 ± 0.00	97.91 ± 0.00	97.75 ± 0.00	97.48 ± 0.00	42.66 ± 0.01
μ_α=0.2	97.59 ± 0.00	97.63 ± 0.00	96.32 ± 0.00	96.31 ± 0.00	41.56 ± 0.01
Regression coefficient
μ_α=0	1.05 ± 0.04	1.05 ± 0.04	1.05 ± 0.04	1.05 ± 0.04	0.97 ± 0.07
μ_α=0.2	1.06 ± 0.02	1.06 ± 0.02	1.06 ± 0.02	1.06 ± 0.02	0.95 ± 0.04

Accuracy was quantified using the average correlation between true breeding value and (genomic) estimated breeding values from 10 replications validated in Generation 5, comprising 8000 individuals with genotypes but no phenotypes.

Bias was quantified using the average regression coefficients of true breeding value on (genomic) estimated breeding values from 10 replications.

The true QTL effects were sampled from normal distributions with mean or and scaled to simulate a trait with a heritability 0.5.

J: includes a covariate for in the model, C: entire matrix of imputed and observed genotype covariates centered, JC: both J and C, N: neither J or C, and PBLUP: pedigree-based BLUP. Covariates were fitted for both 50 QTL and 150 markers.

Effect of more QTL and markers (200 SNP vs. 2000 SNP)

We had hypothesized that the improvement of accuracy from adding an extra covariate for would reduce as the number of QTL increases, because is likely to be closer to zero for a trait that is more polygenic. Table 6 shows that PBLUP was largely unaffected by changes to genetic architecture, but the accuracy of genomic prediction declined slightly as the number of QTL increased. This reflects the fact that the precision of estimating QTL effects is greater when the effects are large, and polygenic traits with more QTL must have on average smaller effects when compared at the same genetic variance. The benefit of fitting in the model was virtually eliminated when the number of substitution effects to estimate increased from 150 to 1500. In contrast to the results for accuracy, there was no impact of QTL number on bias. Centering had no impact on accuracy or bias.

Table 6

Accuracy and bias of genomic prediction () for different numbers of QTL and alternative analyses

Numbers of QTLa	Analysesb
	JCc	Jd	C	N	PBLUP
Correlations (%)
50 QTL + 150 markers	97.59 ± 0.00	97.63 ± 0.00	96.32 ± 0.00	96.31 ± 0.00	41.56 ± 0.01
500 QTL + 1500 markers	90.45 ± 0.01	90.49 ± 0.01	89.99 ± 0.01	89.99 ± 0.01	41.62 ± 0.02
Regression coefficient
50 QTL + 150 markers	1.06 ± 0.02	1.06 ± 0.02	1.06 ± 0.02	1.06 ± 0.02	0.95 ± 0.04
500 QTL + 1500 markers	1.08 ± 0.03	1.08 ± 0.03	1.08 ± 0.03	1.08 ± 0.03	0.98 ± 0.05

Accuracy was quantified using the average correlation between true breeding value and (genomic) estimated breeding values from 10 replications validated in Generation 5, comprising 8,000 individuals with genotypes but no phenotypes.

Bias was quantified using the average regression coefficients of true breeding value on (genomic) estimated breeding values from 10 replications.

The true effects for 50 or 500 QTL were sampled from a Normal distribution with mean μ = 0.2 and scaled to simulate a trait with a heritability 0.5.

J: includes a covariate for μ in the model, C: entire matrix of imputed and observed genotype covariates centered, JC: both J and C, N: neither J or C, and PBLUP: pedigree-based BLUP. Covariates were fitted for either 50 QTL and 150 markers or 500 QTL and 1500 markers.

Centering using the entire matrix of genotype covariates or only the observed genotype covariates

Table 7 shows the accuracies and regression coefficients of TBV on (G)EBV for the genotype panel with 150 markers, , and = 0. The analyses were performed after centering: the entire matrix of imputed and observed genotype covariates (); only observed genotype covariates (), which is the type of centering done in single-step GBLUP; or not centering the covariates (). The accuracy of the genomic analysis with covariates centered as but without J (model C*) was ∼17% lower than the other genomic analyses. However, when J was included in the model with covariates centered as the accuracy of prediction was markedly improved.

Table 7

Accuracy and bias of genomic prediction () when centering for all genotypes or observed genotypes

Analysesa	Correlations (%)b	Regression Coefficientc
JC	82.08 ± 0.06	0.95 ± 0.07
J	82.08 ± 0.06	0.95 ± 0.07
C	82.05 ± 0.06	0.95 ± 0.07
N	82.05 ± 0.06	0.95 ± 0.07
JC*	82.09 ± 0.06	0.95 ± 0.07
C*	65.35 ± 0.06	1.16 ± 0.11
PBLUP	42.66 ± 0.01	0.97 ± 0.07

Accuracy was quantified using the average correlation between true breeding value and (genomic) estimated breeding values from 10 replications validated in Generation 5, comprising 8,000 individuals with genotypes but no phenotypes. The true QTL effects were sampled from Normal distributions with mean μ = 0 and scaled to simulate a trait with a heritability 0.5.

Bias was quantified using the average regression coefficients of true breeding value on (genomic) estimated breeding values from 10 replications.

J: includes a covariate for μ in the model, C: entire matrix of imputed and observed genotype covariates centered, JC: both J and C, N: neither J or C, C*: only observed genotype covariates centered, JC*: both J and C*, and PBLUP: pedigree-based BLUP. Covariates were fitted for 150 markers.

As explained previously, where is the expected value of the covariates in the founders, will tend to zero for the marker panel that does not include QTL, even with However, even if in a population undergoing selection when selected individuals are genotyped, where is the expected value of the observed genotype covariates. In this study, selection was used to increase the mean of the trait. Thus, is expected to be negative because most of the genotyped individuals were from G5, whereas is expected to be zero. The negative estimate of from 10 replicates of the JC* analysis, −2.69, confirms that On the other hand, the mean of from 10 replicates of the JC analysis was −0.75. This explains why fitting J in the model improved the accuracy of genomic prediction when covariates were centered as in single-step GBLUP.

Fernando found that centering using improved the accuracy without J in the model when the population was not under selection and the genotyped individuals were unselected. In that study, mating was random with no selection, so the allele frequency means of the genotyped individuals were a reasonable approximation of the allele frequency means in the founder population. By contrast, our simulation here shows that centering using can reduce the accuracy when the population is under selection, unless J is fitted in the model.

In single-step GBLUP, the observed genotypes are commonly centered by subtracting their mean and used to construct a genomic relationship matrix, for example, by using the first method proposed by VanRaden (2008). Using that genomic relationship matrix in the single-step GBLUP formula in Aguilar does not account for J. This was recognized by Vitezica , who proposed a modification for populations under selection that involved adding a constant to all elements of the genomic relationship matrix that they derived by equating the sum of the elements of the genomic relationship matrix to the sum of the elements of the numerator relationship matrix. In the appendix of that paper they showed that this modification is equivalent to fitting a covariate Q = −J and treating − as a random effect. In addition to this modification, Christensen proposed a multiplicative scaling to the genomic relationship matrix such that its diagonals have the same mean as the diagonals of the numerator relationship matrix. Vitezica claimed that − represents the mean breeding value of selected individuals, and we have confirmed here that this is true provided the observed genotype covariates are centered by their mean.

Most populations are under natural or artificial selection. In many cases, genotypes are only available on selected individuals. In single-step genomic analysis that combine genotyped and nongenotyped individuals in a joint analysis, the mean of observed genotypes are available for centering. If the observed genotypes include QTL, the accuracy of genomic prediction can be severely compromised, unless the J covariate is fitted in the model. If the observed genotypes are only markers, the accuracy of genomic prediction may not necessarily be improved by fitting J in the model, but it doesn’t do any harm. However, if centering is applied only to the observed genotypes, which is the type of centering used in single-step GBLUP, accuracy could be severely compromised, unless the J covariate is fitted in the model or an equivalent approach is adopted.