Genomic Heritability: A Ragged Diagonal Between Bias and Variance

Many important traits in plants, animals, and microbes are polygenic and are therefore difficult to improve through traditional marker-assisted selection. Genomic prediction addresses this by enabling the inclusion of all genetic data in a mixed model framework. The main method for predicting breeding values is genomic best linear unbiased prediction (GBLUP), which uses the realized genomic relationship or kinship matrix (K) to connect genotype to phenotype. The use of relationship matrices allows information to be shared for estimating the genetic values for observed entries and predicting genetic values for unobserved entries. One of the key parameters of such models is genomic heritability , or the variance of a trait associated with a genome-wide sample of DNA polymorphisms. Here we discuss the relationship between several common methods for calculating the genomic relationship matrix and propose a new matrix based on the average semivariance that yields accurate estimates of genomic variance in the observed population regardless of the focal population quality as well as accurate breeding value predictions in unobserved samples. Notably, our proposed method is highly similar to the approach presented by Legarra (2016) despite different mathematical derivations and statistical perspectives and only deviates from the classic approach presented in VanRaden (2008) by a scaling factor. With current approaches, we found that the genomic heritability tends to be either over- or underestimated depending on the scaling and centering applied to the marker matrix (Z), the value of the average diagonal element of K, and the assortment of alleles and heterozygosity (H) in the observed population and that, unlike its predecessors, our newly proposed kinship matrix KASV yields accurate estimates of in the observed population, generalizes to larger populations, and produces BLUPs equivalent to common methods in plants and animals.

The LMM for this analysis is: 105 y = 1 n µ + I n g + where y is the vector of phenotypic LSMs of for n entries, n is 106 the number of entries, 1 n is an n-element vector of ones, µ is 107 the population mean, I n is the identity matrix of size n, g is 108 an n-element vector of random effect values for entries with 109 g ∼ N(0, Kσ 2 g ), and is the residual for the each entry with 110 ∼ N(0, I n σ 2 ). 111 In the sommer R package LMM (3) is expressed as: 112 mmer(fixed = Y~1, random =~vs(G, Gu = K), rcov =~units, data = data) where data is an n × 2 matrix with Y as a column of LSMs and 113 G is a column of factor coded entry IDs (Covarrubias-Pazaran 114 2016). 115 The ASV definition of total variance from LMM (3) is: where θ ASV y is the total phenotypic variance, V = Kσ 2 g + I n σ 2 117 is the variance covariance among observations, θ ASV g is the ge-118 nomic average semivariance, and θ ASV is the average semivari-119 ance of the residuals.

120
The ASV definition of the genomic variance is: = (n − 1) −1 σ 2 g tr(ZZ T ) = (n − 1) −1 tr(K)σ 2 g whereZ = PZ is the mean-centered marker matrix, and 122 K =ZZ T is the realized genomic relationship or kinship ma-123 trix described by VanRaden (2008) omitting the scaling constant 124 2 ∑ j p j (1 − p j ), where p j is the allele frequency of the j-th SNP, 125 which requires HWE to hold (de los Campos et al. 2015), and 126 tr(ZZ T P) = tr(ZZ T ). The trace ofZZ T is a function of het-127 erozygosity in the observed population ( Effect of heritability (h 2 g ), population size (n), and heterozygosity (H) on the accuracy of genomic heritability estimates. Phenotypic observations were simulated for 1,000 samples with n = 250, 500, and 1000 (left to right) genotyped for m = 5, 000 SNPs and the average heterozygosity H = 0, 25, 50, and 75%. The accuracy of genomic heritability estimates (ĥ 2 g ) from the six relationship matrices is shown for true genomic heritability (h 2 g ) = 0.2 (upper panel), 0.5 (middle panel), and 0.8 (lower panel). The upper and lower halves of each box correspond to the first and third quartiles (the 25th and 75th percentiles). The notch corresponds to the median (the 50th percentile). The upper whisker extends from the box to the highest value that is within 1.5 · IQR of the third quartile, where IQR is the inter-quartile range, or distance between the first and third quartiles. The lower whisker extends from the first quartile to the lowest value within 1.5 · IQR of the quartile. The dashed line in each plot is the true value from simulations. Legarra et al. 2018). When the observed population is in HWE, 129 n −1 tr(K) = 1, and when the population is not in HWE due to

132
In the general case, θ ASV is any form of the genomic relationship matrix calculated from Z, 134 without centering, orZ, with centering, because tr(K) = tr(KP).

135
If we assume G = Kσ 2 g , where G is the variance-covariance of 136 GEBVs, it can observed the magnitude tr(K) is directly inverse 137 to the magnitude of σ 2 g , assuming that V = Kσ 2 g + I n σ 2 and 138 that Kσ 2 g and I n σ 2 are uncorrelated. Thus, as different scales are 139 applied to K, the inverse should be applied to σ 2 g .

140
The ASV definition of the residual variance is: 141 θ ASV = (n − 1) −1 σ 2 tr(I n I T n P) = σ 2 Importantly, the genomic variance θ ASV g is on the same scale 142 as the residual variance θ ASV and both are defined such that 143 (4) is true. Note that REML estimates are directly equivalent to 144 ASV estimates when BLUEs or LSMs are used as the response 145 variable y. In 210 many plant quantitative genetic studies, the population sizes 211 are n ≈ 500, which may pose a general problem for variance 212 component and ratio estimation as those variance components 213 tended to have high variability between replicated experiments 214 (Fig 1). For very large populations, common in human and do-215 mesticated animal studies, it is possible to imagine precise (low 216 variance) but inaccurate (high bias) estimates of σ 2 g and h 2 g result-217 ing from different relationship matrices, unless the assumptions 218 of HWE happen to be perfectly met in the study population.

219
Analysis of real-world data confirms that ASV does not impact 220 BLUPs 221 The differences in the prediction accuracy for unobserved indi-222 viduals in the examples we selected, when present, appeared 223 to be negligible and did not lend themselves clearly to "better" 224 or "worse" categorisations (Fig 2). Strandén  and Gianola 2014). Thus, for genomic selection applications, re-240 searchers can apply whichever form of the realized relationship 241 matrix they choose to obtain GEBVs (ĝ) and estimated random 242 marker effects (β) and to calculate prediction accuracy (r(ĝ, y)). 243 While the choice of K does not impact BLUP, the choice of 244 GRM impacts estimates of genomic varianceσ 2 g , genomic her-245 itabilityĥ 2 g , predictive abilityĥ −1 g r(ĝ, y), average prediction er-246 ror variance PEV, and selection reliability 1 − σ −2 g PEV, which 247 all rely onσ 2 g . The PEV of the BLUPs was similar across rela-248 tionships indicating that differences in selection reliability and 249 predictive ability will be driven by differences inσ 2 g that arise 250 from different forms of K. These results are also preceded by 251 Strandén and Christensen (2011) who show that PEV depends 252 on the centering and scaling and that different methods change 253 the amount of certainty, or reliability, in the BLUPs. 254 We observed the expected patterns that were exposed by our 255 simulations (Fig 1) in the real-world data sets (Table 1). When 256 the study populations are fully inbred as in wheat, Arabidop-257 sis, or inbred per se evaluations in hybrid crops, such as maize, 258 tomato, rice, the ASV corrections produce larger VCEs (Table 1) (Table 1) For the Arabidopsis data set (second row), K Y systematically produced singular systems in sommer::mmer() and prediction accuracy was not estimated for either FT10 µ or FT16 µ . The upper and lower halves of each box correspond to the first and third quartiles (the 25th and 75th percentiles). The notch corresponds to the median (the 50th percentile). The upper whisker extends from the box to the highest value that is within 1.5 · IQR of the third quartile, where IQR is the inter-quartile range, or distance between the first and third quartiles. The lower whisker extends from the first quartile to the lowest value within 1.5 · IQR of the quartile.
for LD.

271
The relationship between ASV and normalized genomic rela-272 tionship matrices 273 As our simulations confirmed, the scaling of the trace of K pro-274 foundly affects the estimation of h 2 g (Figure 1), consistent with 275 Speed and Balding (2015) and Legarra (2016), but has virtually 276 no effect on BLUP or genomic prediction accuracy (Fig 2).

299
Recall that the unbiased estimator of a sample variance among breeding values is: whereḡ is the average genetic value of individuals in the ob-300 served population and g i is the genetic value of the i-th observed 301 individual. The population variance, which is considered by 302 Legarra (2016), among breeding values is: We use S 2 g , the population variance, to indicate a different metric 304 relative to s 2 g , the sample variance.

305
Taking expectations of S 2 g with respect to the distribution of 306 g: 307 E(S 2 g ) = n −1 E(g T Pg) = n −1 tr(KP)σ 2 g (10) Thus, the genomic variance in the reference population (S 2 g ) will 308 be a function of the variance component associated with K in 309 the mixed model (Legarra 2016 eral studies. There are likely to be differences between these 315 K ASV and K GN when n is small. When n is reasonably large 316 K ASV ≈ K GN as n → ∞ (Fig 1). This work, Forni et al. traits. We simulated marker effects for all m = 5, 000 loci fol-364 lowing a normal distribution µ = 0 and σ = 1. When multi-365 plied by the marker genotypes and summed, the score is taken 366 as the true genetic value g of each individual. Residuals are 367 simulated with µ = 0 and σ 2 = (1 − h 2 )/(h 2 s 2 g ) to obtain a 368 trait with the desired genomic heritability (Endelman 2011) and 369 s 2  The form proposed by VanRaden (2008) is: whereZ is the marker matrix centered on column means (2p j ), 384 and p j is the minor allele frequency for the j-th SNP. This form 385 assumes HWE and that p j is obtained from a historical refer-386 ence population, and not the observed population. When p j 387 is estimated in the population of observed entries, the center-388 ing by 2p j is equivalent to column centering and K VR only 389 differs from K ASV by a scaling factor, such that: The form of the relationship matrix proposed by Endelman 393 and Jannink (2012) is : where δ ≈ (n/m)CV −2 is a shrinkage factor, CV 2 is the co-395 efficient of variation of the eigenvalues of S, S = m −1ZZT − 396 Z •k Z T •k , S ii is the mean of diagonal elements of S. Notably, 397 at high marker densities, when δ = 0, Endelman and Jannink 398 (2012) is equivalent to VanRaden (2008) when p j is estimated 399 from the observed entries.
where z ij is the j-th SNP in the i-th individuals, z jk is the j-th 403 SNP in the k-th individual when j = k, and m is the number 404 of markers. It is clear that the diagonals are treated differently, 405 relative to the off-diagonals.
where z j is the i-element vector of the j-th SNP.

407
The classical identity-by-state definition is given in Astle et al. 408 (2009) as: We analyzed four publicly available data sets using six methods  Endelman and Jannink (2012)