The training population included 888 accessions belonging to the Cassava Germplasm Bank of Embrapa Cassava and Fruits (Cruz das Almas, Bahia, Brazil). This germplasm comprised 835 landraces and 53 improved varieties. One hundred and eighty accessions were characterized as sweet cassava (< 50 ppm of cyanogenic compounds), 136 as containing intermediary cyanide content (50–100 ppm cyanogenic compounds), 560 as bitter cassava (> 100 ppm cyanogenic compounds), and 12 as unclassified. These accessions were collected from all 26 Brazilian states, with every state represented by at least one genotype. The genotypes were evaluated in the cities of Cruz das Almas and Laje in the state of Bahia, Brazil, in 21 trials over a six-year period (2011 to 2016).

For most experiments, 15–20 cm stem cuttings were planted in double lines during the rainy season in the region (May–July). The experimental plot consisted of two rows of eight plants per row. The rows were 0.9 m apart, while plants in the same row were 0.8 m apart, with 11.52 m² per plot. All recommended cassava cultural practices were employed (as in Souza et al., 2006). Trials were harvested 11–12 months after planting. The traits measured to estimate genomic selection efficiency were: 1) fresh root yield (FRY) at plot level (16 plants) and then adjusted to t.ha^-1, 2) dry matter content in the roots (DMC), according to Kawano et al. (1987), where approximately 5 kg of roots were weighed in a hanging scale (WA) and then, the same sample was weighed with the roots submerged in water (WW). DMC was estimated utilizing the following formula: D M C ( % ) = ( W A W A − W W x 158.3 ) − 142 and 3) dry root yield (DRY) in to t.ha^-1, estimated per plot by multiplying the FRY and DMC.

A joint analysis of 21 trials with complete randomized block design or augmented block design were used to obtain the phenotypic data. Three replicates were used in the complete randomized block design, while in the augmented block design, 10–16 replicates of the common checks were used, with equal distribution of accession number per block. Improved clones (9602-02, 9607-07, 9824-09, 9655-02) and improved varieties (BRS Dourada, BRS Gema de Ovo, and BRS Novo Horizonte) were used as checks in different field trials. More details from the phenotypic dataset could be seen in Table S1 and S2.

Due to unbalanced trials, we obtained the BLUP and deregressed BLUP (Garrick et al., 2009) for each clone. The BLUPs were obtained by the following mixed linear model: y_ijl =μ+c_i +β_j +r_l(j)+ϵ_ijl in which y_ijl is the vector of phenotypic observations; c_i is the clone random effect with c i ~ N ( 0 , I σ ^ c 2 ) β_j is the combination of location and year, assumed as fixed effect; r_l(j) s the replication nested within location and year, assumed as random effect with r j ( l ) ~ N ( 0 , I σ ^ r 2 ) and ϵijl is the residual with ϵ i j l ~ N ( 0 , I σ ^ e 2 ) The deregressed BLUPs were estimated by: d e r e g r e s s e d   B L U P = B L U P 1 − P E V σ ^ c 2 Garrick et al., 2009), where the PEV is the prediction error variance of each clone and σ ^ c 2 s the clonal variance component. The package lme4 (Bates et al., 2015) in R software version 3.5.2 (R Core Development Team, 2018) was used to obtain the BLUPs and deregressed BLUPs for each clone.

DNA was extracted from cassava leaves following the CTAB (cetyltrimethylammonium bromide) protocol described by Doyle and Doyle (1987). To evaluate DNA integrity and standardize its concentration, 1.0% (w/v) agarose gels were stained with ethidium bromide (1.0 mg L^-1) for visual comparison of a series of DNA phage Lambda (Invitrogen) concentrations. The DNA samples were sent to the Genomic Diversity Facility at Cornell University (http://www.biotech.cornell.edu/brc/genomic-diversity-facility) for genotyping-by-sequencing (GBS) (Hamblin and Rabbi, 2014). Genotypic data were selected using a minimum call rate of 0.90 and the missing markers were imputed by Beagle 4.1 software (Browning and Browning, 2016). Finally, SNPs with minor allele frequency (MAF) > 0.05 were retained. After applying marker quality control, 48,655 SNPs were selected for genomic prediction.

The genomic best linear unbiased prediction (G-BLUP), Reproducing Kernel Hilbert Spaces (RKHS), and Bayes B prediction methods were evaluated, considering the additive (A) and additive-dominant (A+D) genetic models, except RKHS, which predicts genetic effects based on non-parametric—and thus neither additive nor dominance—covariances. The additive-dominant genetic model of G-BLUP is expressed: y_d =Jμ+Za+Hd+ϵ where yd is the deregressed BLUP vector; µ is the general mean; a is the additive effect vector, random a ~ N ( 0 , G σ ^ a 2 ) d is the dominant deviation effect vector, random d ~ N ( 0 , D σ ^ d 2 ) ϵ is the residual effect vector, ϵ ~ N ( 0 , I σ ^ e 2 ) J, Z and H are the incidence matrices for µ, a and d, respectively, as COV(a,d)=0 The additive relationship matrix G was: G = Z Z ' 2 ∑ ​ p i ( 1 − p i ) in which Z is the marker matrix (-1, 0 and 1) and p_i is the major allele frequency of i marker. Two additive-dominant genetic models were tested for the G-BLUP method, the Classical (Vitezica et al., 2013) and the Genotypic (Su et al., 2012), differing in the parameterization of the genomic relationship matrix due to dominance. The Classical dominant relationship matrix was parameterized by the following Vitezica et al. (2013):

The Genotypic dominant relationship matrix was estimated by the following equation (Su et al., 2012):

For the Bayes B method, the complete conditional prior distribution was used: y d i | a j , d j , Z i × j , H i × j ~ N ( μ + ∑ j Z i × j a j + ∑ j H i × j d j , σ ^ e 2 ) in which yd is the deregressed BLUP vector; µ is the general mean; a_j nd d_j re the additive and dominant marker effects, both random a j | σ ^ a j 2 ~ N ( 0 , I σ ^ a j 2 ) d j | σ ^ d j 2 ~ N ( 0 , I σ ^ d j 2 ) and COV(a_i ,d_i )=0 Z and H are the incidence matrix of a_j and d_j respectively.

The model of the RKHS method was: y_d =Jμ+Xg+ϵ where yd is the deregressed BLUP vector; µ is the general mean; g is the genotypic effect vector, random g ~ N ( 0 , K σ ^ g 2 ) ϵ is the residual effect vector, ϵ ~ N ( 0 , I σ ^ e 2 ) J and X are the incidence matrix of µ and g, respectively. K is a gaussian matrix estimated by: K = e x p ( − h D m e d i a n ( D ) ) h is the reduction coefficient to K values (in this work h was equal to 1), and D is the Euclidian distance of Z codified marker matrix (Gianola et al., 2006; Crossa et al., 2010).

The 5-fold cross-validation with three repetitions was performed to estimate the following parameters: 1) predictive ability ( r ^ y ^ y = C O R ^ ( P r e d ^ V a l , B L U P V a l ) ) in which P r e d ^ V a l are the genomic estimated breeding values (GEBVs) for additive genetic models, or genomic estimated genotypic values (GEGVs) for additive-dominant and RKHS models, and BLUP_Val are the BLUPs from the validation population; 2) bias ( b ^ = C O V ^ ( P r e d ^ T r a i n , B L U P T r a i n ) / σ ^ P r e d T r a i n 2 ) in which P r e d ^ T r a i n are the genomic estimated breeding values (GEBVs) for additive genetic models, or genomic estimated genotypic values (GEGVs) for additive-dominant genetic models, of the training population, BLUP_Train are the BLUPs from the training population, σ ^ P r e d T r a i n 2 s the variance of the GEBVs for additive genetic models, or genomic estimated genotypic values (GEGVs) for additive-dominant genetic models of the training population; 3) broad-sense genomic heritability ( H ^ 2 = σ ^ g 2 / ( σ ^ g 2 + σ ^ e 2 ) ) in which σ ^ g 2 s the genomic variance, σ ^ e 2 s the residual variance; 4) narrow-sense genomic heritability ( h ^ 2 = σ ^ a 2 / ( σ ^ g 2 + σ ^ e 2 ) ) which σ ^ a 2 s the additive genomic variance, σ ^ g 2 s the genomic variance, σ ^ e 2 s the residual variance. For each replicate of the cross-validation process, the population was split into five equal folds. Five genomic predictions were performed per fold used as test set (no phenotypes) each fold was predicted by the remaining four-folds training set (with phenotypes).

The sommer R package (Covarrubias-Pazaran, 2016) was used to fit the G-BLUP and RKHS models, while the BGLR R package (Perez and De Los Campos, 2014) was used to fit the Bayes B model. All methods were performed using R software version 3.5.2 (R Core Development, 2018). For Bayes B method, we ran 20,000 Markov Chain Monte Carlo (MCMC) iterations with the burn-in of the initial 4,000 iterations and thinning of 10, we applied different priori for π for each trait and genetic model, these values were previously estimated by Bayes Cπ (Table S3).

The training-validation partitions of the population used in cross-validation were set up to be identical across prediction models, using the set.seed() function of R software version 3.5.2 (R Core Development, 2018). The residual variances of Markov Chain Monte Carlo (MCMC) of the Bayes B method were used to evaluated the MCMC convergency by the Raftery and Lewis’s convergence diagnostic (Raftery and Lewis, 1992) applied in coda R package (Plummer et al., 2006).

Analysis of variance was performed to estimate the effects of the genomic selection methods for predictive ability and bias estimates for DMC, FRY, and DRY. These analyses were performed using the lme4 R package (Bates et al., 2015).

The following mixed model was used to estimate the efficiency of the genomic selection methods: y_ijk =m_i +s_jk +e_ijk which y is the dependent variable, as predictive ability and bias; mi is the mean of the genomic selection method I, assumed as fixed effect; s_jk is the effect of cross validation of the replication j and fold k, assumed as random effect s ~ ( 0 , σ ^ c v 2 ) and eijk is the residual effect of the i genomic selection method of the j replication and k fold, e ~ ( 0 , σ ^ e 2 ) The genomic prediction means were submitted to the Tukey multiple comparison test implemented in the emmeans R package (Russel, 2018).

The Cohen’s Kappa coefficient (Cohen, 1960) was used to analyze the coincidence of clone selection by the different genomic selection methods, considering a selection proportion (SP) amplitude ranging from 5–30%. The coincidence selection was performed using a binary code and the selected and unselected individuals received code “1” and “0”, respectively. The Kappa coefficient and coincidences selection index were calculated using R.

In general, the inclusion of the dominant genetic effects increased the genomic variance explained by the markers (Table 1), and reduced the genomic additive variance and residuals (Table 1 and Figure S1). Smaller changes in the broad-sense genomic heritability were observed for DMC, except for the Bayes B method, which demonstrated the highest broad-sense genomic heritability among the prediction methods with an additive-dominant genetic model.

Insert Table 1

A predominance of additive effects for DMC was identified with the G-BLUP method (Table 1 and Figure S1), while for FRY and DRY the dominant effects prevail. The Bayes B method showed the highest estimates of broad-sense genomic heritability and genomic variance components. However, the variation of the broad-sense genomic heritabilities between traits was smaller, suggesting a relatively large proportion of dominance variance. Even with the highest broad-sense genomic heritability, the Bayes B A+D method exhibited smaller narrow-sense genomic heritability than the G-BLUP A+D method, regardless of the dominant relationship matrix used (Table 1). However, all the additive-dominant genetic models overestimated the broad-sense genomic heritability because it was higher than the phenotypic heritability (0.337, 0.351, and 0.545 for FRY, DRY, and DMC, respectively).

The additive-dominant genetic models showed higher predictive ability than additive models and RKHS method for yield traits (FRY and DRY, Figure 1). The highest predictive ability was demonstrated by the G-BLUP A+D classical method (average of 0.484 for FRY and 0.492 for DRY), followed by Bayes B A+D (average of 0.479 for FRY and 0.488 for DRY). In addition, the predictions of dominant effects in genetic models for yield traits reduced the bias estimate, with the smaller bias at Bayes B method (Figure 1). The RKHS method showed the highest bias estimates for all traits.

Significant differences between the genomic selection methods with different genetic models were identified for predictive ability and bias for all agronomic traits except the predictive ability of DMC (Table 2). Although there were no significant differences in the predictive ability between the genomic selection methods with additive-dominant models, the G-BLUP A+D classical method showed the highest predictive ability for FRY (0.483) and DRY (0.492) (Table 2). Bayes B A+D and RKHS methods did not show significant differences for predictive ability in comparison with G-BLUP A+D classical method for DRY. On the other hand, for FRY only Bayes B A+D and G-BLUP A+D genotypic methods did not show significant differences with the G-BLUP A+D classical method.

Among the methods with non-additive effects, the G-BLUP A+D classical was significantly different from the RKHS method for DRY but not for FRY. As the RKHS method can predict additive and partial epistatic effects (Gianola et al., 2006; Crossa et al., 2010), it is possible that the epistatic effects were more important for FRY than DRY, as the RKHS method did not show a significant difference with the additive genetic models G-BLUP A and Bayes B A (Table 2).

DMC showed the highest phenotypic heritability and predictive ability of traits. However, there was no improvement in predictive ability when the additive-dominant genetic models were used to predict this trait, which reinforced the theory that DMC in cassava has a high influence from additive effects. On the other hand, for FRY and DRY, the additive-dominant models demonstrated increased predictive ability, suggesting a greater importance of dominant effects for these traits in cassava.

Although significant differences were detected between genomic prediction methods with different genetic models by ANOVA and Tukey’s mean test (Table 2), the expected genetic gains for genomic prediction were still smaller than those obtained by phenotypic selection, with expected selection gains equivalent to 67.5%, 67.1%, and 69.4% of the phenotypic selection for FRY, DRY, and DMC, respectively (Figure 2). Although selection gains with genomic predictions were similar for all traits, the non-additive genetic models, such as Bayes B A+D, RKHS, and G-BLUP A+D classical and genotypic, increased the gain by an average of 0.69 t/ha for FRY and 0.24 t/ha for DRY in comparison with the additive genetic models. For DMC, the differences between the selection gains of genomic prediction methods were lower (average of 0.04%), because there was no significant difference between the clone prediction methods for this trait (Table 2). Moreover, the selection differential for DMC in the roots was lower than for yield traits due to the smaller trait amplitude (17–38%).

There was a great uniformity in the differences between the selection gains of the phenotypic BLUP and the predicted gains in the different selection proportions, with a mean difference of selection gain of 6.18% and 7.79% of the Bayes B A+D model for FRY and DRY, respectively (Figure 2 and Table S2). For DMC, there were lower gains differences between the phenotypic selection and genomic prediction, with the largest difference observed in the G-BLUP A method (average of 1.40% of genetic gain) compared to others (Figure 2 and Table S4).

The genomic expected selection gain and its relative efficiency to phenotypic expected selection gain were calculated. According to Oliveira et al. (2012), the conventional breeding cycle of cassava is at least four years due to the need to include phenotypic information from a minimum of four breeding phases (clonal evaluation trial, preliminary yield trial, advanced yield trial, and uniform yield trial). The efficiency and the selection gains per time unit to simulate early selection assisted by genomic selection were calculated. The efficiency was determined by comparing time required to recombine the selected clones as parents in a conventional breeding program vs. one assisted by genomic selection.

Genomic selection based on the G-BLUP A+D classical method for FRY and G-BLUP A for DMC was more efficient than phenotypic selection when the breeding cycle was ≤ 0.60 of the conventional breeding cycle (Table 3 and Figure S2). However, for DRY the genomic selection was more efficient than phenotypic selection only with a selection proportion of 5–15%.

Reducing breeding cycle time by 60% using genomic selection could result in gains of 65%, 69%, and 77% over those provided by phenotypic selection for FRY, DRY, and DMC, respectively, in a selection proportion of approximately 15% of the best clones (Table 3). If the breeding cycle was reduced to 20% of the conventional breeding cycle (four years to ten months), the genetic gains would be 163%, 170%, and 183% over those provided by phenotypic selection for FRY, DRY, and DMC, respectively.

The selection proportion affected significantly the relative efficient of the genomic prediction only in breeding cycle time reductions biggest then 35% of the conventional Cassava breeding cycle (Table 3).

According to Oliveira et al. (2015b), heritability estimates can assist selection strategies in increasing genetic gain, as well as defining the breeding method and experimental design. Given the broad- and narrow-sense genomic heritability, the G-BLUP A+D classical method showed that cassava yield traits demonstrate a predominance of dominant effects. In addition, the broad-sense genomic heritability of G-BLUP A+D was closer to the phenotypic heritability values (0.337 for FRY, 0.351 for DRY, and 0.545 for DMC [Table 1]). Stability of FRY and DRY are important agronomic attributes for any cassava variety to ensure high market competitiveness in the starch industry, especially as there is a minimum acceptable DMC threshold for processing the raw material. Roots with DMC index below this threshold are not processed by the starch industry due to the high industrial cost and low starch yield.

Knowledge about trait heritability and variation gained during field evaluation in different environments may assist in optimizing selection of cassava breeding programs, with the goal of developing new cassava varieties with higher starch yield stability. Optimizing the selection proportion and evaluated traits in each breeding phase can maximize the probability of selecting the best clone. This is because low heritability traits such as FRY and starch yield are generally evaluated in the final breeding phases due to greater stem cutting availability (more plants per plot across multiple locations).

Wolfe et al. (2016b) also related the predominance of additive and dominant deviation effects for DRY and FRY, respectively. They found similar estimates of broad- and narrow-sense heritability for the first genomic selection cycle of IITA population using the G-BLUP A+D method (0.12 and 0.35 for narrow- and broad-sense heritability, respectively, for FRY, and 0.47 and 0.52 for narrow and broad-sense heritability, respectively, for DMC). Wolfe et al. (2016b) found that the appropriate genetic model for DMC was the additive-dominant, while in the present study the additive-dominant models obtained similar results to the other models’. The reduction of the variance explained by the additive component was noted previously in cassava (Wolfe et al., 2016b) and other species such as Pinus taeda L. (Muñoz et al., 2014) and hybrids of Eucalyptus urophylla and E. grandis (Bouvet et al., 2016). Several authors reported that during prediction using additive genetic models, part of the dominant deviation was predicted along with the additive effects; however, when using additive-dominant genetic models, this dominant deviation predicted by the additive effects is then computed by the dominant variance (Zuk et al., 2012; Vitezica et al., 2013; Muñoz et al., 2014; Wolfe et al., 2016a). According to Vitezica et al. (2013) genetic models with assumptions of additive and dominant deviation effects result in better genomic predictions.

There were significant differences in predictive ability between the methods for FRY and DRY, mainly due to different genetic models (additive and non-additive). The additive-dominant genetic models showed higher predictive ability than additive genetic models for FRY and DRY. Among the genomic selection methods, the G-BLUP A+D classical (Vitezica et al., 2013) had high predictive ability and low bias, statistically similar to other additive-dominant genetic models. Therefore, the additive-dominant genetic models allow for exploration of part of the non-additive effects by increasing cassava clone selection accuracy. Other authors have evaluated relationship matrices with classical and genotypic dominant models and verified the lack of differences in the genomic predictions of these matrices, although the broad-sense heritability has been somewhat lower in the matrix (H*) of genotypic dominant (Vitezica et al., 2013; Wolfe et al., 2016a). However, the correlation between the additive and dominant parameters was higher in the G-BLUP genotypic method in comparison with the classical one (Wolfe et al., 2016a), which corroborates the correlations found in this work.

Dominance effects occurs due to the interaction between alleles at the same locus and its main benefits are expected in crossbreeding, since dominance has been suggested as one of the genetic mechanisms explaining heterosis (Shull, 1908). Indeed, hybrid vigor for yield components in cassava over better-parent values has been reported (Parkes, et al., 2020), indicating that heterosis should be explored in order to develop superior cassava genotypes. Therefore, genomic predictions for traits such as FRY and DRY must be based on the assumption that non-additive effects are an important component that should be considered in the predictions to optimize crossing designs, such as in mate-pair allocation (Almeida Filho et al., 2019). As a further step, the role of dominance effects on the genetic architecture of FRY and DRY should be evaluated in the breeding population generated in this study.

For DMC, there was no significant difference between genetic models for predictive ability, although the additive genetic models showed better prediction abilities than additive-dominant models. Similar results were reported in the first genomic selection cycle of IITA cassava population (Wolfe et al., 2016a). According to Denis and Bouvet (2013) the decision of which genetic model to use in genomic selection depends on the training population as well the traits under selection. Specifically, in cassava, the genomic prediction of FRY and DRY would be more efficient if applying additive-dominant genetic models, while for DMC the additive models are satisfactory. Among genetic models, additive gene action had highest response to selection. Therefore, in the case of DMC, the population improvement focused on genetic additive effects can achieve large medium‐to‐long term genetic gains.

Incorporating non-additive effects into the genetic model reduces the additive genomic variance and the bias of the GEBVs, as well increasing the accuracy for selecting the best parent (Vitezica et al., 2013; Wolfe et al., 2016a). On the other hand, some simulated studies did not find any differences in predictive ability of GEBVs between additive and additive-dominant genetic models (Almeida Filho et al., 2016; Heidaritabar et al., 2016). Therefore, it is expected that non-additive effects prediction may increase the genetic gains for yield traits of new cassava varieties in the breeding programs (Muñoz et al., 2014; Wolfe et al., 2016a).

The expected selection gains for FRY and DRY were high due to the training population being composed of germplasm accessions with high genetic variability (Figure S2) for several traits, including yield traits (Oliveira et al., 2015a; Oliveira et al., 2016). Within a group of clones that deviated from the FRY and DRY mean (Figure S2), some were from high-yield, improved varieties (FRY potential of > 30 t/ha) such as BRS Novo Horizonte, BRS Poti Branca, BRS Kiriris, and BRS Tapioqueira.

According to the Bayes B A+D method, if the cassava breeding cycle was reduced by 40%, the genomic selection gains would be on average 12.48% and 11.92% higher than phenotypic selection gains for FRY and DRY, respectively. A similar observation was made for DMC (22.16%). Oliveira et al. (2012) reported that with a 25% reduction in breeding cycle time, the relative efficiency of RR-BLUP genomic prediction was 4.6%, 15.96%, and -7.05% for FRY, starch yield, and DMC, respectively. According to these authors, higher selection gains may be achieved by reducing the time required to identify and recombine the parents in the breeding cycle. These results may assist in the planning of genomic selection implementation to increase the frequency of new cassava varieties with good agronomic traits and adaptations to new biotic and abiotic stress challenges.

Reducing the selection proportion is not feasible as reducing breeding cycle time to improve genetic gain, but it may be the next milestone to improve genetic gain in cassava breeding, by increasing the number of clones evaluated in earlier stages by genomic prediction.

A previously recommended method for clone and parent selection in seedling trial phases was assessment of the harvest index (the ratio of FRY and the biomass yield), used for FRY indirect selection in seedling nursery trials and clonal evaluation trials (Kawano et al., 1998) due to its high correlation (0.730). However, when analyzing the historical data (2000–2013) of the cassava breeding program at the International Center for Tropical Agriculture (CIAT), Barandica et al. (2016) found very low correlation between FRY and harvest index (0.11). In addition, the harvest index assessment is more labor-intensive than the FRY evaluation alone. According to Barandica et al. (2016) the correlation of FRY between the clonal evaluation trials and the uniform yield trials was 0.29, while in the present study the correlation between the GEBVs and the uniform yield trials for the G-BLUP A+D classical method was 0.483. In future studies, the efficiency of genomic selection in the seedling trial phase for FRY could be better understood by determining realized genetic gains.

Ceballos et al. (2016a) stated that one issue in the selection of good parents is the high intra-family genetic variability due to high heterozygosity. Thus, new cassava varieties may derive from crosses between parents with low agronomic performance. Indeed, Kawano et al. (1998) evaluated almost 327,000 clones from 4,130 crosses during 14 years of research, and among all those evaluated clones only three were officially released as new varieties.

Commonly the standard methods used for parent selection are the per se performance and, less commonly, general combining ability (Ceballos et al., 2004). Unfortunately, there is no linear relationship between the per se performance and the progeny’s breeding values due to dominant deviation (Ceballos et al., 2004). In addition, the diallel analysis in cassava breeding programs is problematic because the crosses in cassava are laborious and usually imbalanced due to issues of flowering synchronization (Ceballos et al., 2017), the considerable unpredictability of the flowering season (between four to ten months after planting), and the time for seed maturity after harvest demanding at least one year to obtain the seeds of controlled crossings (Ceballos et al., 2004). Su et al. (2012) reported that genetic models with additive and non-additive effects prediction might allow for exploitation of specific combining ability. Therefore, applying genomic selection with genetic models that consider both genetic effects may be a faster alternative for selecting clones for advancement in the breeding pipeline, parents for crossings, inheritance studies, and variation of traits at the different stages of the cassava breeding program.

Another strategy for selecting promising parents is the pedigree-based best linear unbiased prediction (P-BLUP) method for predicting breeding values (Ceballos et al., 2016a). This strategy attempts to estimate breeding values after obtaining clone phenotypic data. According to Piepho et al. (2008) this method allows for dissection of the genotypic value in additive and non-additive effects, and it can be approached by identity-by-descent (P-BLUP) or identity-by-state (G-BLUP) information. However, this method requires a large amount of kinship information, which is not always available once several crosses have been carried out between germplasm accessions with no kinship data available. Nevertheless, the lack of kinship/pedigree information can be efficiently compensated for by identity-by-state (IBS) performed using an additive relationship matrix proposed by VanRaden (2007), as Hayes et al. (2009) considered the additive genomic relationship matrix as accurate as the kinship matrix. Bouvet et al. (2016) found that the G-BLUP prediction method had higher predictive ability than P-BLUP in several genetic models. According to Zhang et al. (2015) the GEBVs may be even more accurate when using a genetic architecture-enhanced relationship matrix for each trait, with the parametrization of relationship matrix composed by markers with high effect for the trait.

Using G-BLUP for breeding value estimation at preliminary, advanced and uniform yield trials, we assume that there is a genetic correlation between clones due to relationship-by-state (Piepho et al., 2008). Since cassava is vegetatively propagated, the additive genomic matrix may be used as a genetic covariance matrix for selecting promising parents by applying mixed models in the different breeding phases (such as the clonal evaluation trial, preliminary yield trial, advanced yield trial, and uniform yield trial). As the correlation between breeding values vs. genotypic values is not perfect (0.716 of selection coincidence at 13.3% selection proportion for FRY, and 0.690 of selection coincidence at 8.3% selection proportion for DRY), the coincidence in the selection of clones to be used as parents and for advancement in the breeding program tends to be low. Therefore, by using the genetic covariance matrix in mixed models, the selection of parents with high potential to generate promising clones would be performed based on their breeding values even if the clones had low genotypic value and/or low per se performance. This strategy can increase the parent selection accuracy, estimate the narrow-sense heritability, and predict the GEBVs and GEGVs across the field trials, assisting parent and clone selection, respectively.