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 μg 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) μg when marker covariates were not centered; and with (JC) or without (C) μg when all observed and imputed marker covariates were centered. Centering did not influence accuracy of genomic prediction, but fitting μg 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 μg in the model improved accuracy, but had little impact when the panel contained only markers. In populations undergoing selection, fitting μg 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 et al. 2014). 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 et al. 2014). 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 et al. (2014) proposed a solution for the marker effects model that involves fitting an additional fixed covariate that estimates the mean m g of the linear component of the genotypic value, which is denoted by a i 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 et al. 1989; genotypes, using their sample means can be done in addition to the Fernando et al. (2014) approach of fitting m g : 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 et al. (2011). 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 et al. 2011). 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 m g 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 m g changes from the mean of unselected individuals to become the mean breeding value in selected individuals, as claimed by Vitezica et al. (2011).

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 m i 9 denote the row vector of genotypes for individual i, which is often coded as 21,0,1. Then, under additive gene action, the genotypic value, g i ; which is the expected phenotypic value of an individual with genotypes m i 9, can be written as where b is the value of g i when m i 9 ¼ 09 and a is the vector of substitution effects. The scalar b and the vector a are constants, but g i will be a random variable because of randomness in a i ¼ m i 9a; owing to the randomness in the genotypes for a randomly sampled individual. Note that the expected value of the linear component a i of the genotypic value in Equation 1 is ; 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 ðb þ m g Þ is a constant, representing Eðg i Þ; and u i ¼ ðm i 9 2 k9Þa 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, k9; from the observed values, m i 9: In both Equations 1 and 2, a j is the same substitution effect for locus j. The intercepts in these models, however, are different. In Equation 1, the intercept is b, and it is the value of g i when m i 9 ¼ 09: In Equation 2, on the other hand, the intercept is ðb þ m g Þ; and it is the value of g i when m i 9 ¼ k9: More generally, k9 is not known, so genotypes are coded by subtracting a different vector v9 from the observed genotypes as m i 9 2 v9: Still, a j is the substitution effect for locus j, but the intercept will change to become ðb þ v9aÞ; which is the value of g i when m i 9 ¼ v9: Thus, as more rigorously shown in (Strandén and Christensen 2011), inference about a 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 et al. 2009) or explicitly (Fernando et al. 2014) using best linear prediction. Let M g denote the matrix of genotypes for individuals that were genotyped. Then, the genotypes of the individuals with missing genotypes are imputed as where A ng is the matrix of pedigree-based additive relationships between the nongenotyped and genotyped individual, and A gg is the matrix of additive relationships among genotyped individuals. Now, the model for the genotypic values, when genotypes are coded as in Equation 2, becomes where e is that part of g n that cannot be imputed from knowledge of the breeding values of genotyped relatives. In practice, the true value of k9 is not known, and data for its estimation may not be available. Rearranging these equations in terms of the uncentered M g rather than the centered matrix of genotype covariates ðM g 2 1k9Þ; results in and substituting m g ¼ k9a; as previously defined, and lettinĝ M n ¼ A ng A 21 gg M g ; results in which suggests that m g ¼ k9a could be treated as an unknown constant and estimated as a fixed effect from the data (Fernando et al. 2014). The covariate vector for m g is denoted by J n ¼ 2 A ng A 21 gg 1 for nongenotyped individuals, and by J g ¼ 2 1 for genotyped individuals. So, Equation 3 becomes which can be combined as ; and the 0 matrix in the last term of Equation 4 is required because e does not appear in the model for genotyped individuals. When the vector a represents the substitution effects of a large number of loci containing positive and negative effects, m g ¼ k9a will tend to have a value close to zero. Accordingly, we have simulated some scenarios with positive m a ¼ Eða i Þ so that the entire a vector is positive to exacerbate the impact of m g ¼ k9a: 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 m a ¼ Eða i Þ ¼ 0; in a population undergoing selection, it is expected that Eða i Þ ¼ Eðm i 9Þa 6 ¼ 0 in nonfounders. Suppose v9 is the mean of the observed genotype covariates in such a population undergoing selection, and these means are used to center the matrix M g of observed genotypes. Then, the model for the genotypic values can be written in terms of the matrix M Ã g ¼ M g 2 1v9 of centered covariates as and using J for the covariate corresponding to m Ã g; Equation 5 can be written as which can be combined as Note that the regression coefficients for J; m g in Equation 4 and m Ã g þ v9a in Equation 6, must be equal. This implies that m Ã g ¼m g 2 v9a: Similarly, the intercepts of these two models must be equal too, and this implies b Ã ¼ b þ v9a:

Simulations
Phenotypic and genotypic data were simulated using XSim (Cheng et al. 2015), 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 m a = 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 (h 2 ) of 0.1, 0.3, or 0.5. An additional scenario with m a = 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 m g , 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(Henderson , 1984: where y is a vector of phenotypic observations, 1 is a vector of 1s, m ¼ b þ m g is a general mean, u is a vector of random additive genetic effects, e is a vector of random residual effects, and Z is a known incidence matrix relating observations to u: In this model, EðuÞ ¼ 0; EðeÞ ¼ 0; so that EðyÞ ¼ 1m: Further, VarðuÞ ¼ As 2 g ; with A being the numerator relationship matrix, and VarðeÞ ¼ Is 2 e ; so that VarðyÞ ¼ ZAZ9s 2 g þ Is 2 e : Single-step Bayesian regression model: All genomic EBVs (GEBV) were obtained by single-step Bayesian regression (Fernando et al. 2014), using BayesC priors for marker effects with p = 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 m g (model J) was: where y n and y g are vectors of phenotypes for nongenotyped and genotyped individuals, 1 is a vector of 1s, m is a general mean, m g is the expected value of the linear component a i of the genotypic value if selection was absent, a is a vector of random substitution effects of markers, e is a vector of imputation residuals, Z n and Z g are incidence matrices relating the breeding values of nongenotyped and genotyped individuals to their phenotypes, J g ; which is defined for genotyped individuals, is a vector of 21s, J n ; which is defined for nongenotyped individuals, is a vector computed as A ng A 21 gg J g ;M n is the matrix of imputed marker covariates, M g is the matrix of observed marker covariates, and e n and e g are vectors of random residual effects for nongenotyped and genotyped individuals. This model can be represented as where y ¼ y n y g ; and e 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 m g was in the model, and whether or not the columns in the marker covariate matrix M 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).
Accuracy of genomic prediction was quantified using the correlation of TBV and GEBV (r g;ĝ ), 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 (b g;ĝ ). In models JC and J, GEBVs are obtained using Equation 24 in (Fernando et al. 2014): whereĝ is the GEBV,m g 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 M g contains the uncentered number of copies of the A allele at each locus, and the uncentered version is used to imputeM n : In model JC, the entire matrix of imputed,M n ; and observed, M g ; genotype covariates is centered. In models C and N, the GEBVs are computed in a corresponding manner, except that the covariate J and its coefficient m g 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 m a ¼ 0:2 to construct phenotypes with h 2 ¼ 0:5: The four analyses JC, J, C, and N were repeated using only genotype panels comprising QTL plus markers for three other scenarios: 200 loci, h 2 ¼ 0:1; m a ¼ 0:2; 200 loci, h 2 ¼ 0:3; m a ¼ 0:2; and 200 loci, h 2 ¼ 0:5; m a ¼ 0:2: Fernando et al. (2014) had observed that including J g in the model was necessary only when m a 6 ¼ 0 and the number of observations exceeds the number of markers. So, to examine the impact of m a ; the scenario with 200 loci and h 2 ¼ 0:5 was repeated with m a ¼ 0:0; and to examine the impact of the number of loci, the scenario with h 2 ¼ 0:5 and m a ¼ 0:2 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 et al. 2009;Aguilar et al. 2010;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, m g ¼ k9a will be close to zero, especially when the panel consists of only marker loci. However, in a population undergoing selection, m Ã g ¼ m g 2 v9a; which is the mean breeding value of genotyped individuals, is not expected to be zero even when the panel only includes marker loci and m a ¼ 0: 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, h 2 ¼ 0:5; and m a ¼ 0; when the matrix of observed genotype covariates were centered using their means as: M Ã g ¼ M g 2 1v9; where v9 ¼ 1=n g 19M g ; and n g 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 M g of uncentered covariates as shown by Equation 6. We a 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 m a = 0.2 and scaled to simulate a trait with a heritability 0.5.
b J: includes a covariate for m g 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. c The analyses were based on fitting covariates for only 50 QTL, only 150 markers, or both 50 QTL and 150 markers.
will compare the estimates of m Ã g in this model with those of m g 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 m a = 0.2 and h 2 = 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 1 m9a from the model equation, and this has no effect on the mixed model solutions for a: 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 m a = 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 m g in the model estimates the relevant mean from the data.
In panels that included causal variants (QTL), fitting m g 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 m g in the model had little impact. This is because in the population that was simulated here the founders were not selected, and thus m g 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 J 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 m was $9.90 for the analyses using J and N. The estimate of m g was 4.98 for the analysis using J. For the genotyped individuals, the covariate values in J g are all 21, so 1m + J gmg is a vector of values equal to 9.9024.98 = 4.92. For nongenotyped individuals, the covariate values in J n can vary widely, but many are close to 21 while others are close to 0. This means that 1m þ J nmg will include values that range from 9.90 to 4.92, accounting for variation in accuracy of imputation. When m g 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 n a 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 m a = 0.2 and scaled to simulate a trait with heritabilities 0.1, 0.3 or 0.5.
b J: includes a covariate for m g 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.
n a 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 m a = 0.2 and scaled to simulate a trait with a trait with a heritability 0.5. b The analyses were based on fitting covariates for only 50 QTL, only 150 markers, or both 50 QTL and 150 markers. c J: includes a covariate for m g 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.
which casem g will be close to zero as confirmed in the above mentioned replicate where the estimate was 21.18.
Here, results from the analysis with only QTL on the panel are used to show that m g is the mean of a i in the founder population and not the mean of the breeding value u i of the genotyped individuals. Recall that the QTL model was used with an intercept of b ¼ 0:0 to simulate the data. Thus, when only QTL are on the panel, the true value of b is zero. In analysis J because m ¼ ðb þ m g Þ; bothm andm g are estimates of m g ; and could be pooled, which for the replicate above would be ð9:91 þ 8:19Þ=2 ¼ 9:05: In that replicate, the actual mean of a i 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 u i in the 9000 genotyped individuals was 2.4, which is clearly not near the pooled estimate of 9.05 for m g : 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 u i 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 m g parameter corresponding to the covariate J is the mean of a i in the founder population and not the mean breeding value u i in the selected individuals. In analysis JC with the covariates centered, the intercept b is the value of g i when ðm9 i 2 v9 i Þa ¼ 0; which is the case when m i 9 ¼ v i 9: The estimatem was 16.11 in this analysis, butm g remained about the same value, namely 7.80. This shows thatm g 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 m g 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 m g 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.
Bias: Table 4 shows the regression coefficients of TBV on (G)EBV for h 2 ¼ 0:5 and m a = 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.
Effect of mean QTL effect (m a ¼ 0 vs. m a ¼ 0:2) We had hypothesized that the impact of omitting m g from the model would be greatest when m g departs significantly from 0, which is more likely to occur when m a departs from 0. For that reason, our base simulation used m a = 0.2. Results are shown in Table 5 for the panel including QTL and markers with h 2 = 0.5 for m a = 0.2 compared with m a = 0. These results confirmed that the benefit of fitting m g was greatest when m a = 0.2, but there was still an advantage to fitting m g when m a = 0. That advantage is likely to erode as the number of QTL increases.
n a 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. b Bias was quantified using the average regression coefficients of true breeding value on (genomic) estimated breeding values from 10 replications. c The true QTL effects were sampled from normal distributions with mean m a ¼ 0 or m a ¼ 0:2 and scaled to simulate a trait with a heritability 0.5. d J: includes a covariate for m g 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.
n a 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. b Bias was quantified using the average regression coefficients of true breeding value on (genomic) estimated breeding values from 10 replications. c The true effects for 50 or 500 QTL were sampled from a Normal distribution with mean m a = 0.2 and scaled to simulate a trait with a heritability 0.5. Changing the mean QTL effect had no impact on bias, except for a slight influence on PBLUP.
Effect of more QTL and markers (200 SNP vs.

SNP)
We had hypothesized that the improvement of accuracy from adding an extra covariate for m g would reduce as the number of QTL increases, because m g 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 m g 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.
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, h 2 ¼ 0:5, and m a = 0. The analyses were performed after centering: the entire matrix of imputed and observed genotype covariates (M Ã ¼ M 2 ð1=ðn g þ n n ÞÞ119M); only observed genotype covariates (M Ã g ¼ M g 2 ð1=n g Þ119M g ), which is the type of centering done in single-step GBLUP; or not centering the covariates (M Ã ¼ M). The accuracy of the genomic analysis with covariates centered as M Ã g 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 M Ã g ; the accuracy of prediction was markedly improved.
As explained previously, m g ¼ k9a; where k9 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 m a 6 ¼ 0: However, even if m a ¼ 0; in a population undergoing selection when selected individuals are genotyped, m Ã g ¼ m g 2 v9a 6 ¼ 0; where v9 is the expected value of the observed genotype covariates. In this study, selection was used to increase the mean of the trait. Thus, m Ã g is expected to be negative because most of the genotyped individuals were from G5, whereas m g is expected to be zero. The negative estimate ofm Ã g from 10 replicates of the JC Ã analysis, 22.69, confirms that m Ã g , 0: On the other hand, the mean ofm g from 10 replicates of the JC analysis was 20.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 et al. (2014) found that centering using M Ã g 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 M Ã g 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 et al. (2010) does not account for J. This was recognized by Vitezica et al. (2011), 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 = 2J and treating 2m Ã g as a random effect. In addition to this modification, Christensen et al. (2012) 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 et al. (2011) claimed that 2m Ã g 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.

ACKNOWLEDGMENTS
The authors are grateful to Bruce L. Golden and Hao Cheng for assistance in the implementation of SSBR models, to Jack C. M. Dekkers for his constructive comments on design of the simulation, and to the very helpful anonymous reviewers. R.L.F. acknowledges useful discussions with Andres Legarra and Zulma Vitezica.  a 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 m a = 0 and scaled to simulate a trait with a heritability 0.5. b Bias was quantified using the average regression coefficients of true breeding value on (genomic) estimated breeding values from 10 replications. c J: includes a covariate for m g 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.