Marker-Based Estimates Reveal Significant Nonadditive Effects in Clonally Propagated Cassava (Manihot esculenta): Implications for the Prediction of Total Genetic Value and the Selection of Varieties

In clonally propagated crops, nonadditive genetic effects can be effectively exploited by the identification of superior genetic individuals as varieties. Cassava (Manihot esculenta Crantz) is a clonally propagated staple food crop that feeds hundreds of millions. We quantified the amount and nature of nonadditive genetic variation for three key traits in a breeding population of cassava from sub-Saharan Africa using additive and nonadditive genome-wide marker-based relationship matrices. We then assessed the accuracy of genomic prediction for total (additive plus nonadditive) genetic value. We confirmed previous findings based on diallel crosses that nonadditive genetic variation is significant for key cassava traits. Specifically, we found that dominance is particularly important for root yield and epistasis contributes strongly to variation in cassava mosaic disease (CMD) resistance. Further, we showed that total genetic value predicted observed phenotypes more accurately than additive only models for root yield but not for dry matter content, which is mostly additive or for CMD resistance, which has high narrow-sense heritability. We address the implication of these results for cassava breeding and put our work in the context of previous results in cassava, and other plant and animal species.


INTRODUCTION 44 45
Understanding genetic architecture requires the decomposition of genetic effects 46 into additive, dominance, and epistatic components (Fisher 1918;Cockerham 1954;47 Kempthorne 1954). However, partitioning genetic variance components is notoriously 48 difficult, requiring specialized breeding designs (e.g. diallel crosses) and pedigree 49 information (Lynch and Walsh 1998) often limiting the genetic diversity that can be 50 sampled in any one given study. Genome-wide molecular marker data now enable the 51 accurate measurement of relatedness in the form of genomic realized relationship 52 matrices (GRMs) (VanRaden 2008;Heffner et al. 2009;Lorenz et al. 2011b). GRMs 53 provide more accurate relatedness information than pedigrees because they directly 54 measure Mendelian sampling (causing variation in relatedness within relatedness classes 55 such as full-siblings) (Heffner et al. 2009). Further, GRMs can measure relationships 56 even in diverse, nominally unrelated samples expanding the potential for studying 57 inheritance in natural and breeding populations (Lorenz et al. 2011a). 58 Estimation of narrow-sense heritability and prediction of breeding values in 59 genomic selection programs is becoming increasingly common using additive 60 formulations of GRMs (Visscher et al. 2008). Several recent studies have described 61 dominance and epistatic GRMs for the partitioning of non-additive genetic variance using 62 genome-wide SNP markers (Su et al. 2012;Vitezica et al. 2013;Muñoz et al. 2014;63 Here, y represents raw phenotypic observations. In our data, the only fixed 198 effect (β) was an intercept for all traits except RTWT, which contained a covariate 199 accounting for variation in the number of plants harvested per plot.  a  s  s  o  c  i  a  t  e  d  v  a  r  i  a  n  c  e  c  o  m  p  o  n  e  n  t  a  s  w  e  l  l  a  s  a  r  e  p  l  i  c  a  t  i  o  n  e  f  f  e  c  t  ,  n  e  s  t  e  d  i  n  l  o  c  a  t  i  o The 212 terms X, Z loc.year , Z rep , Z add , Z dom and Z epi are incidence matrices relating observations to 213 the levels of each factor. We list the different models fit in Table 1, each of which are  214 variations on the full model described above. 215 The formulation described above was used to fit the GG historical data in a single 216 model. For the C1 progenies only a single trial was available and therefore we fit all data 217 together. Since the C1 trial was conducted across three locations but with no replications 218 we fit the same model for C1 as GG excluding the replication term. The models described 219 above were fit using the regress package in R (R Core Team 2015). The regress function 220 finds REML solutions to mixed models using the Newton-Raphson algorithm. 221 For each trait, in both the C1 and GG we identified a "best fit" model among the 222 models listed in Table 1 based on the lowest Akaike Information Criterion (AIC; 2*k -223 2*log(likelihood), where k = number parameters estimated). We also examined the log-224 likelihood of each model and the total genetic variance explained (H 2 ). The precision of variance component estimates and the dependency among estimates was examined using 226 the asymptotic variance-covariance matrix of estimated parameters (V). Specifically, we 227 report standard errors for each variance component,defined as the square root of the 228 diagonal of V. We also converted V into a correlation matrix (F,as in (Muñoz et al. 229 2014)), where F is defined as L -1/2 VL -1/2 and L is a diagonal matrix containing one over 230 the square root of the diagonal of V. We use F to assess the dependency of variance 231 components estimates, especially for comparing results among traits and populations. 232 233

Within-trial analyses 234
We used only a subset of the GG trials to estimate variance components in a 235 single multi-environment model. In addition, we leveraged the entire historical GG data 236 by analyzing each trial (N=47, unique location-year combinations) separately. This 237 provided us with 47 estimates of additive, dominance and epistatic variance and we 238 examine the distribution of variance components estimates. As in the multi-environment 239 models, within-trial models were fit with regress in R. 240 241

Genomic prediction and cross-validation 242
We assessed the influence that modeling non-additive genetic variance 243 components have on genomic prediction using a cross-validation strategy. Because 244 single-step multi-environment models are computationally intensive, we used a two-step 245 approach here. In the first step, we combined data from all available GG and C1 trials 246 using the following mixed model: , where the covariance 252 structure was the identity matrix (I). The BLUP (ĝ) for the clone effect therefore 253 represents an estimate of the total genetic value for each individual. The mixed model 254 above was solved using the lmer function of the lme4 R package. 255 In our data, the number of observations per clone ranges from one to 131 with 256 median of two and mean of 5.97 excluding the checks TMEB1 and I30572 which had 257 941 and 902 observations respectively. Pooling information from multiple years and 258 locations, especially when there is so much variation in numbers of observations can 259 introduce bias. Much theoretical development, particularly in animal breeding has been 260 done to address this issue, and we followed the approach recommended by Garrick et al. 261 (Garrick et al. 2009 We implemented a 5-fold cross-validation scheme replicated 25 times to test the 268 accuracy of genomic prediction using the genomic relationship matrices and models 269 described above (Table 1). In this scheme, for each replication, we randomly divided the dataset into five equally sized parts (i.e. folds). We used each fold in turn for validation 271 by removing its phenotypes from the training population and then predicting them. We 272 calculated accuracy as the Pearson correlation between the genomic prediction and the 273 BLUP (ĝ, not-deregressed) from the first step. For each model, we calculated accuracy 274 both of the prediction from the additive kernel (where present) and the total genetic value 275 prediction, defined as the sum of the predictions from all available kernels (e.g. additive 276 Our first assessments of non-additive genetic effects in cassava were single-step 285 multi-environment models implemented in both the Genetic Gain (GG) and the Cycle 1 286 (C1) populations. The five models (Table 1) were initially compared using the AIC. 287 Tables 2 and 3 show the model results including AIC and variance components for the 288 GG and C1, respectively. For HI, MCMDS, and SPROUT the best models were A + D, A 289 + D + AD and A + D + AA, respectively. For these three traits, the best model according 290 to AIC was the same between the GG and C1. For DM the additive only model was best 291 in the GG but an A + D model was selected in C1. For RTWT, in the GG an A + D + AA 292 model was selected but in the C1 the dominance only model was selected. Finally, for MCBBS the additive only model was best in the GG but A + D + AD was selected in the 294

C1. 295
For every trait, when comparing the model achieving the highest broad-sense 296 heritability (H 2 ), we saw that the H 2 increased from GG to C1. This can be seen most 297 easily in Figure 1, which shows how total explainable genetic variance (H 2 ) is partitioned 298 among variance components in the C1 and GG (also see Tables  For HI, the additive only model was very similar between populations (h 2 = 0.34 306 in GG and 0.32 in C1). But for the best model, A + D, more of the total genetic variance 307 is partitioned to C1 (0.28) than in the GG (0.07). Like for DM, the three component 308 models (A + D + AA and A + D + AD) explained the most variance for the C1 (H 2 = 309 0.49), but the epistatic components had large standard errors (26.4±18.7 and 24.4±17.6 310 respectively). 311 RTWT was strongly non-additive in both populations. In the GG the A + D + AA 312 model, genetics explained 25% of the phenotypic variance with non-additive variances 313 (D + AA) explained over half of that amount (16% collectively) and the epistatic term 314 was significantly different from zero (0.033 ± 0.014). In the C1, 31% of the phenotypic 315 variance was explainable by dominance alone. 316 model, and 0.25 in the C1 for the A + D + AD model). While the limited apparent genetic 318 variance in the GG was best explained additively, in the C1 92% of total genetic variance 319 was non-additive. 320 For MCMDS, the additive model had the highest H 2 (0.79) in the GG but the AIC 321 best model, A + D + AD explained most (0.89) in the C1. The AIC best model for both 322 C1 and GG was the same. In both cases, significant AxD epistasis was detected, 323 explaining 44% of the variance in the GG and 18% in the C1. The additive variance 324 increased from GG (0.25) to C1 (0.64). 325 SPROUT was best modeled with A + D + AA in both populations. Broad-sense 326 heritability increased from 0.22 in GG to 0.39 in C1 for this trait. The dominance 327 component was not significant (5.8 ± 8.6) in the GG but was in the C1 (31.7 ± 26.7). 328 We examined the asymptotic correlation matrices of parameter estimates (F) to 329 ascertain the dependency of variance component estimation. Correlation matrices for 330 every trait + model combination are provided for GG in Tables S3-S6 and for the C1 in 331 Tables S7-S9. The correlation between genetic variance components was always negative 332 and was, in general, higher in the GG compared to the C1. Correlations between additive 333 and dominance components were highest in the A+D models (range -0.81 to -0.83 in the 334 GG and -0.5 to -0.61 in the C1). Correlations between A and D components dropped in 335 models with epistasis (range -0.42 to -0.63, GG and -0.26 to -0.58, C1). Correlations 336 between additive and AxA epistatic variances (range -0.09 to -0.29) and AxD epistasis 337 (range -0.07 to -0.22) were low. Correlations between dominance components and epistasis were higher ranging from -0.28 to -0.64 with AxA epistasis and -0.36 to -0.69 339 with AxD epistasis. 340 341

Within-trial analyses 342
We also examined variance partitioning within each of 47 GG trials for the 5 343 models described in Table 1. The mean and variability of model parameters (variance 344 components, heritability, etc.) across these trials are summarized in Table S10. Figure 2  345 provides a visual summary of the proportion of phenotypic variability explained by each 346 genetic variance component on average across the trials. We also compared the mean 347 AIC across trials and found them to agree overall with the results of the one-step multi-348 environment models (Table 2, Table S10). Specifically, the models that fit best in the 349 one-step models were best on average in the within trial analyses for the following traits: 350 DM (additive), HI (A+D), and MCMDS (A+D+AD). However, for RTWT (A+D instead 351 of A+D+AD) and MCBBS (D instead of A) the within trial AIC-best models were 352 different on average from the one-step multi-environment models. 353

Genomic Prediction of Additive and Total Genetic Value 355
We used cross-validation to assess the accuracy of genomic prediction of additive and 356 total genetic value for the five models (Table 1) Table S11-S12). The prediction accuracy for the additive kernel 359 was higher on average in the C1 for DM, MCBBS, MCMDS, and SPROUT compared to 360 the GG but lower for HI and RTWT. Additive accuracies by trait (across models) were highest for DM (range: 0.54-0.60), followed by MCMDS (range: 0.44-0.54) and HI 362 (range: 0.32-0.45). The rest of the traits had similar additive accuracies dependent on the 363 population and model. RTWT had the overall lowest accuracy (range: 0.09-0.36). 364 Across all multi-kernel analyses conducted, prediction of total genetic value was on 365 average 42% more accurate compared to the prediction of the additive kernel in the same 366 model. Compared to the single-kernel additive prediction however, total genetic value 367 predictions were an average of only 6% better (maximum of 26% improvement; Figure 4, 368 Tables S11-S12). By model, improvements in the correlation between total value and 369 phenotype over the additive only model were 5%, 7% and 8% for A+D, A+D+AA and 370 A+D+AD respectively. The additive only model predictions were on average 11% less 371 accuracy in the C1 than in the GG. Additive kernel predictions from models with non-372 additive components were 7% lower in C1 compared to GG. Total genetic value 373 predictions were also less accurate by 12% in the C1 relative to GG. 374 The models we fit for multi-kernel genomic prediction involved the estimation of 375 weight parameters corresponding to the partitioning of genetic variance among the 376 kernels. Overall, there was a tendency for non-additive kernels (mean of 0.40 for 377 dominance, 0.40 AxA epistasis and 0.46 AxD epistasis) to get more weight than additive 378 (mean of 0.32). Additive prediction accuracies were positively correlated (r = 0.81) with 379 weight placed on the additive kernel. The accuracy of total genetic value prediction was 380 negatively correlated (r = -0.64) with weight placed on non-additive kernels (Tables S11-381 S12). 382

DISCUSSION 385
In clonally propagated crops, non-additive genetic effects can be effectively 386 exploited by the identification of superior genetic individuals as varieties. For this reason, 387 we quantified the amount and nature of non-additive genetic variation for key traits in a 388 genomic selection breeding population of cassava from sub-Saharan African. Then we 389 assessed the accuracy of genomic prediction of additive compared to total (additive plus 390 non-additive) genetic value. Using several approaches and datasets based on genome-391 wide marker data, we confirmed previous findings in cassava based on diallel populations, 392 that non-additive genetic variation is significant, especially for yield traits (Cach et al.
. Results from the present study suggest 436 that the combination of genomic selection and hybrid breeding strategies should increase 437 the rate of genetic gain for complex traits such as yield. However, initial investment in 438 identification of complementary heterotic groups with good specific combining ability is 439 required. 440 One of the more interesting aspects of our study relative to previous ones is the 441 comparison between a parental generation (the Genetic Gain) and their offspring (Cycle 442 1), a collection of full-and half-sib families. From GG to C1, the H 2 generally increased. 443 For RTWT, MCBBS and HI this is largely attributable to increased non-additive variance 444 and CMD for which non-additive variance dropped in C1 relative to GG. In contrast to 445 our result, theory suggests that reduction (or fixation) of allele frequencies at some loci Based on the mean diagonal of the kinship matrix, C1 (0.66) does not appear notably 451 more inbred than GG (0.64). We also calculated mean pairwise LD (GG = 0.27, C1 = 452 0.29) and mean LD block size (21.7 kb in GG and 23.1 kb in C1) using PLINK (version 453 1.9, https://www.cog-genomics.org/plink2) and found the two generations to be similar. 454 Probably the strongest explanation for the difference in genetic variance 455 components between GG and C1 is the family structure (137 full-sib families from 83 456 outbred parents). In a population of half-sibs ¾ of the dominance variance is expressed 457 within families and all of it for full-sib populations (Hallauer et al. 2010;Ceballos et al. 458 2015). Indeed, increasing the number of full-sib relationships is known to increase the 459 non-additive genetic variance detectable in a population (Varona et al. 1998;Van Tassel 460 et al. 2010). 461 It is also conceivable that maternal plant effects could increase apparent non-462 additive effects in C1. The C1 clones in contrast to the GG clones are new, and were 463 derived from stem cuttings of seedling plants germinated in the previous field season 464 (2012)(2013). The suggestion is therefore that the quality and vigor of the seedling plant, 465 giving rise to the C1 clones may influence their performance in the 2013-2014 trial. We 466 further caution that comparison of GG and C1 may be biased by the disproportionate 467 amount of data available for the GG. 468 In our study, when additive and non-additive kernels were used together, the 469 ability of the additive kernel to predict the phenotype decreased, suggesting that the 470 additive kernel by itself was absorbing non-additive variance. Estimates of additive 471 genetic variance have previously been shown to capture some non-additive effects (Lu et 472 al. 1999;Su et al. 2012;Zuk et al. 2012;Muñoz et al. 2014). Predictions models that do 473 not explicitly incorporate non-additive components may therefore achieve gains in the 474 short-term that break down over the long-term (Cockerham and Tachida 1988; Walsh breeding) values may therefore provide a less biased, more accurate selection of parents 477 for genomic selection (Costa E Silva et al. 2004;Palucci et al. 2007). We note that our 478 predictions of total genetic value were focused on parametric models based in 479 quantitative genetic theory. However, many non-parametric and non-linear approaches 480 are available (e.g. RKHS and random forests) that may capture even more non-additive 481 variation than found in our study.    proportion of the total phenotypic variance explained by dominance (d 2 ), additive-by-695 additive epistasis (i 2 A#A ), additive-by-dominance epistasis (i 2 A#D ) and broad-sense 696 heritability (H 2 ) are provided. Model log-likelihoods and Akaike Information Criterion 697 (AIC) are also given. The best model for each trait (lowest AIC) is highlighted in grey. 698 699 Table 3. Single-step multi-environment model results for the Cycle 1 population. 700 Variance components (± standard errors), narrow-sense heritabilities (h 2 ), proportion of 701 the total phenotypic variance explained by dominance (d 2 ), additive-by-additive epistasis 702 (i 2 A#A ), additive-by-dominance epistasis (i 2 A#D ) and broad-sense heritability (H 2 ) are 703 provided. Model log-likelihoods and Akaike Information Criterion (AIC) are also given. 704 The best model for each trait (lowest AIC) is highlighted in grey.