Genomic Selection in Rubber Tree Breeding: A Comparison of Models and Methods for dealing with G × E

Several genomic prediction models incorporating genotype × environment (G×E) interactions have recently been developed and used in genomic selection (GS) in plant breeding programs. G×E interactions decrease selection accuracy and limit genetic gains in plant breeding. Two genomic data sets were used to compare the prediction ability of multi-environment G×E genomic models and two kernel methods (a linear kernel (genomic best linear unbiased predictor, GBLUP) (GB) and a nonlinear kernel (Gaussian kernel, GK)) and prediction accuracy (PA) of five genomic prediction models: (1) one without environmental data (BSG); (2) a single-environment, main genotypic effect model (SM); (3) a multi-environment, main genotypic effect model (MM); (4) a multi-environment, single variance GxE deviation model (MDs); and (5) a multi-environment, environment-specific variance GxE deviation model (MDe). We evaluated the utility of GS with 435 rubber tree individuals in two sites and genotyped the individuals with genotyping-by-sequencing (GBS) of single-nucleotide polymorphisms (SNPs). Prediction models were estimated for diameter (DAP) and height (AP) at different ages, with a heritability ranging from 0.59 to 0.75 for both traits. Applying the model (BSG, SM, MM, MDs, and MDe) and kernel method (GBLUP and GK) combinations to rubber tree data showed that models with the nonlinear GK and linear GBLUP kernel had similar PAs. Multi-environment models were superior to single-environment genomic models regardless the kernel (GBLUP or GK), suggesting that introducing interactions between markers and environmental conditions increases the proportion of variance explained by the model and, more importantly, the PA. In the best scenario (well-watered (WW / GK), an increase of 6.7 and 8.7 fold of genetic gain can be obtained for AP and DAP, respectively, with multi-environment GS (MM, MDe and MDS) than by conventional genetic breeding model (CBM). Furthermore, GS resulted in a more balanced selection response in DAP and AP and if used in conjunction with traditional genetic breeding programs will contribute to a reduction in selection time. With the rapid advances in and declining costs of genotyping methods, balanced against the overall costs of managing large progeny trials and potential increased gains per unit time, we are hopeful that GS can be implemented in rubber tree breeding programs.


Several genomic prediction models incorporating genotype x environment (G×E) interactions 20
have recently been developed and used in genomic selection (GS) in plant breeding programs. 21 G×E interactions decrease selection accuracy and limit genetic gains in plant breeding. Two 22 genomic data sets were used to compare the prediction ability of multi-environment G×E 23 genomic models and two kernel methods (a linear kernel (genomic best linear unbiased 24 predictor, GBLUP) (GB) and a nonlinear kernel (Gaussian kernel, GK)) and prediction accuracy 25 (PA) of five genomic prediction models: (1) one without environmental data (BSG); (2) a single-26 environment, main genotypic effect model (SM); (3) a multi-environment, main genotypic effect 27 model (MM); (4) a multi-environment, single variance GxE deviation model (MDs); and (5) a 28 multi-environment, environment-specific variance GxE deviation model (MDe). We evaluated 29 the utility of GS with 435 rubber tree individuals in two sites and genotyped the individuals with 30 genotyping-by-sequencing (GBS) of single-nucleotide polymorphisms (SNPs). Prediction 31 models were estimated for diameter (DAP) and height (AP) at different ages, with a heritability 32 1 Introduction 47 Generally, the rubber tree breeding program is characterized by breeding cycles of 25-30 years 48 and includes the production of crosses, evaluation, and selection of field progeny, and 49 propagation of selected superior material (Gonçalves et al., 2006). Compared to animal and 50 annual crop breeding, forest tree breeding is still in its infancy, and the most advanced programs 51 are in their third or fourth cycle of breeding, with very little differentiation of the bred 52 populations from natural populations (Isik, 2014). Rubber tree breeding programs are complex 53 and costly because the large size of trees requires experiments over large tracts of land to test 54 progeny, and the progeny tests are expensive to establish, manage over many years, and evaluate 55 via measurement. 56 The main objective of rubber tree breeding is the development of early selection methods that 57 support the accurate prediction of mature phenotypes at a younger stage and are therefore 58 important for shortening breeding cycles and, in the end, improving the cost efficiency of such 59 breeding programs. Hevea breeding needs to significantly reduce the time taken to derive a 60 clone, Priyadarshan (2017) proposed two strategies: (1) to cut short the breeding steps being 61 followed by conventional means and (2) to inculcate genomics into breeding programmes 62 specially to identify high-yielding genotypes in half-sibs, full-sibs and poly-cross seedlings 63 during juvenile stage that can minimize both space and time. 64 Traditional plant breeding programs depend mainly on phenotypes being evaluated in various 65 environments; selection and recombination are based solely on the resulting data plus pedigree 66 information, when available. Genomic selection (GS), a new approach using whole-genome 67 molecular markers, has the potential to quickly improve complex traits with low heritability, 68 significantly reduce the cost of the line and hybrid development and increase grain production in 69 less time to improve quantitative traits in large plant breeding populations (Meuwissen et al. 70 2001). 71 Genomic prediction combines marker data with phenotypic and pedigree data in an attempt to 72 increase the accuracy of the prediction of breeding and genotypic values. The method depends on dense genome-wide marker coverage to produce genomic estimated breeding values 74 (GEBVs) from an ensemble analysis of all markers. 75 According to Lorenz et al. (2011) , the accuracy of GS, which is measured as the correlation  76  between the GEBVs and true breeding values, is affected by the relationship between the training  77 and test sets, the number of individuals in the training set, linkage disequilibrium (LD) between 78 markers and quantitative trait loci (QTLs), the distribution of underlying QTL effects, the 79 statistical method used to estimate the GEBVs, and the trait heritability. 80 GS was proposed by Meuwissen et al. (2001)  In the rubber tree breeding program, pedigree-based analysis has been widely used to evaluate 87 field experiments, estimate genetic parameters, and predict breeding values (Furlani et al., 2005). 88 However, due to the decreasing costs of genotyping thousands or millions of markers and the 89 increasing costs of phenotyping (Krchov and Bernardo, 2015), GS is arising as an alternative 90 genome-wide marker-based method to predict future genetic responses. 91 Appropriate GS methods provide accurate predictions even for untested genotypes, allowing 92 considerable progress in breeding programs by reducing the number of field-tested genotypes 93 and, consequently, the costs of phenotyping (Krchov and Bernardo, 2015). The benefits of GS 94 are more evident when traits are difficult, time-consuming, expensive to measure, and several 95 environments need to be evaluated. 96 The objective of this paper was to evaluate the predictive capacity of GS implementation in 97 rubber trees using linear and nonlinear kernel methods and the performance of such prediction 98 when including GxE interactions in each of the four models described by Bandeira et al. (2017). 99 Thus, for all data sets, we fitted models with a linear kernel using the genomic best linear 100 unbiased predictor (GBLUP) (GB) or nonlinear Gaussian kernel (GK) with a bandwidth 101 parameter estimated according to (Pérez-Elizalde et al., 2015). We also compared the prediction 102 accuracy (PA) of the two kernel regression methods for the four models. The raw data were processed, and SNP calling was performed using TASSEL 5.0 (Glaubitz et  133 al., 2014). Initially, the FASTQ files were demultiplexed according to the assigned barcode. The 134 reads from each sample were trimmed, and the tags were identified using the following 135 parameters: a kmer length of 64 bp, minimum quality score within the barcode and read length of 136 20, minimum kmer length of 20 and a minimum count of reads for a tag of 6. All sequence tags 137 from each sample were aligned to the reference rubber tree genome (Tang et

GS analysis 146
For each character, the phenotypic analysis was carried out jointly for all years of evaluation 147 using the mixed model approach. across all environments and the interaction for each specific environment. 159 The GS models implemented with arrays of GK and GB pedigrees for the AP and DAP traits 160 were implemented in breedR software (Munõz and Sanchez, 2017); using frequentist statistics 161 with the function remlf90, em method, 5 folds and 5 repetitions, the training population (TRN)  162 was created with 4 folds, whereas the test population (TST) was created with one fold. The PA 163 was obtained from the correlation between the predicted BLUPs and the observed BLUPs. 164 For AP and DAP, five statistical prediction models were fitted to all data sets to study their PA 165 using random cross-validation (CV) schemes. The main objective was to compare the prediction 166 ability of the two proposed multi-environment G×E genomic models. 167 The PA of the two kernel regression methods was also compared for single environments and sets for all the traits, and the phenotypic data were centered and standardized. These analyses 174 were performed to derive estimates of variance components. The following variance components 175 resulting from the residual effects, main genetic effect, and genetic environment-specific effects 176 of the four models described above for the trait (LW -low water, and WW -well-watered) data 177 sets were computed. All models were fitted with GxE interactions using the software BGGE 178 . 179

Assessing PA by random cross-validation (CV) GxE 180
The PA of the SM model-method combinations was evaluated with 80% of the hybrids 181 comprising the TRN set, the remaining 20% of the individuals comprising the TST set and none 182 of the lines to be predicted in the TST set in the TRN set using 5 random partitions arranged in 5 183 folds with 100 random partitions each. This procedure was performed separately in each 184 environment, namely, LW and WW, and the SM models were fitted separately for each 185 environment. 186 In the multi-environment models, the PA of the model-method combinations was generated 187 using two different CV designs (Burgueño et al., 2012

Estimates of genetic parameters using SNP genotyping 230
With the genotyped SNPs, the population structure was assessed using a principal component 231 analysis (PCA), and the plots indicated that the 411 genotypes fell into two major clusters 232 (Supplementary Figure 1), which mainly contained hybrids derived from PR255 × PB217 and 233 hybrids derived from GT1 x RRIM701 and of GT1 x PB235 crosses. The first two PCs explained 234 19.51% and 2.18% of the total variance, respectively, clearly splitting the groups along the x and 235 y-axes. 236

Descriptive statistics 237
The genetic correlations between DAP and AP in both environments, namely, LW and WW, 238 were positive and significant, ranging from 0.72 (AP-LW x DAP-WW) to 0.99 (DAP x DAP-239 LW) 240 To assess how much the phenotypic variation is genetically controlled and thus efficient for GS, 241 we first estimated the broad-sense heritability (H) of DAP and AP. The heritability ranged from 242 0.60 to 0.75 for both traits, and when we analyzed the variation separately in each environment, 243 namely, LW and WW, the heritability varied between 0.33 and 0.34 for DAP and 0.41 and 0.42 244 for AP (Table 1). 245 For all the evaluated characters, environmental variation, progeny and the interaction between 246 these two factors were highly significant, indicating that the environments used were contrasting, 247 there was genetic variability among genotypes, and these genotypes presented differential 248 performance according to the environment, respectively. The coefficients of experimental 249 variation obtained in the joint analysis were following those reported in the literature for the 250 same characters, indicating that the experiments were performed with good precision. All 251 components of the variance were nonzero, and estimates were obtained with reasonable 252 accuracy, which was verified by the confidence intervals of these estimates (Table 1). 253

Estimates of variance components 254
Estimates of variance components for each of the GS models derived from the full data analysis 255 are presented in Table 2. 256 For the SM model of the two traits in each environment, the estimated residual variance 257 components for the GK method were smaller than those for the GB method (Table 2), and in 258 compensation, the variance in genetic effects in each environment was slightly more substantial 259 for the SM-GB model-method combination than for the SM-GK combination. Both the genetic 260 variance and phenotypic variance were more significant in LW than in WW. 261 The inclusion of the interaction term (G×E) when using the MM, MDe and MDs induced a more 262 significant reduction in the estimated residual variance for both traits (DAP and AP) ( Table 2). In 263 the six model-method combinations, the residuals were even lower for DAP with the GK 264 method, whereas for AP, the values were similar. 265 For the variance component associated with the genetic interaction effect, the values were, on 266 average, 65% lower for the GK method than for the GB method. 267

Assessment of PA 268
The PA without environmental data (BSG) is shown in Figure 1. The estimated correlations 269 between correlated phenotypes and predictions obtained in the CV test are shown in Figure 2 for 270 the single-environment model (SM) and Figure 3 for the multi-environment models (MM, MDs  271 and MDe). 272

PA without environmental data (BSG) 273
For DAP, the PA was not significantly different between the GB and GK models without multi-274 environmental data. The same result was obtained for GB and GK (0.18), whereas for DAP, 275 there was a small difference between GK (0.26) and GB (0.27) (Figure 1). 276

Single environment (SM) 277
For all traits, the random CV (CV1) applied to only one environment (LW or WW). The results 278 for AP and DAP showed that the PA of the SM-GK model-method combination was higher than 279 that for the SM-GB combination in both LW and WW. The PAs for AP in the LW conditions 280 were 0.17 for SM-GB and 0.18 for SM-GK and in the WW environment were 0.19 for SM-GB 281 and 0.20 for SM-GK. The results for DAP were 0.17 in LW for the SM-GB, and 0.19 for the 282 SM-GK and in WW environment were 0.27 for SM-GB and 0.28 for SM-GK. 283

Multi-environment (MM, MDe, and MDs) 284
The PA varied considerably between the CV1 and CV2 conditions, with an average PA 72.6% 285 higher for CV2 than for CV1 (Supplementary Figure 3). Considering only random CV2, the PA 286 was slightly higher in WW conditions for both characters (AP and DAP), ranging from 0.87 287 (MDs-GK-WW) to 0.80 (GK-MDe-LW). 288 The results obtained with the model-method combinations were very similar; generally, for the 289 LW environment, the best model was GK, without any difference between the methods. In DAP- for MM to 0.87 for MDe and MDs, respectively, for both the GK and GB (Figure 3). 296

EGG 297
For rubber tree breeding, the investigated alternative breeding strategies differed considerably in 298 the number of years required to finish one breeding cycle. For classic improvement, we consider 299 a minimum time of ten years for the beginning of the selection of the best genotypes. In the case 300 of GS, we consider three years for the first selection. 301 To determine the best EGG, we used the best scenario found in the analyses, so we chose to 302 work with the GK matrix and the phenotypic data from the WW environment; however, we 303 present the other results in Supplementary Table 2 In recent years, many statistical models have been proposed for applying GS in plant and animal 313 breeding programs. However, to the best of our knowledge, these models have not yet been 314 applied to rubber tree breeding programs. 315 The efficiency of early selection depends mainly on the early-mature correlation and the 316 heritability of juvenile traits. In this study, all pairwise correlations between environments were 317 strong and positive. This finding is important because the G×E model has the limitation of better 318 and more efficient prediction when applied to subsets of environments that have positive and 319 similar correlations (Crossa et al., 2016). 320 In a study by Moreti et al. (1994), with estimates of genetic parameters and expected gains with 321 the selection of juvenile characters in rubber tree progenies using classical breeding, some 322 parameters stood out positively (rubber production, bark thickness, and stem circumference). 323 Gonçalves et al. (1996) observed the same behavior of the results obtained by Moreti et al. 324 (1994), showing a correlation and its applicability in the selection process. Strong phenotypic 325 and genetic correlations were observed between yield and stem diameter, indicating the 326 possibility of obtaining young clones of good productive capacity and great vigor (Gonçalves et 327 al., 1984). Plants with rapid development of trunk circumference may be more productive, and 328 this may be a useful feature with which to predict more productive hybrids via GS, together with 329 the fact that latex production has a good heritability, better than growth in diameter, because the 330 influence of the rootstock is lower in the production, this will be very important in future studies 331 with this population. 332 The multi-environment genomic prediction was successfully implemented using a GBLUP 333 model; however, depending on the genetic architecture of the trait and germplasm, nonlinear semiparametric approaches, such the GK, could produce better accuracies (Cuevas et al., 2016). 335 We compared two models that used the GB and GK methods, and similar values were found for 336 the PA in all tested conditions. 337 Using genomic prediction without multi-environment data (BSG) for both GB and GK produced 338 the same results (0.18) for AP and a small change from 0.26 to 0.27 between GK and GB, 339 respectively. When a single-environment model was applied to only one environment (LW or 340 WW), the best PA values were in the environment of higher water availability. For AP, we 341 obtained the best values (0.20) using the GK matrix. On the other hand, for DAP, the best PA 342 was obtained via GK matrix (0.28) (Figure 2). 343 In this study, the PA of multi-environment models was assessed by applying the CV strategy. 344 The CV2 validation strategy performed better than the CV1 strategy when applied to the multi- Interactions in field trials affect mature selection and early selection; thus, when evaluating the 363 effectiveness of early selection, it is imperative to determine whether the GxE interaction among 364 sites has a meaningful impact on the early-mature genetic correlation. 365 The application of eight combinations of four models (SM, MM, MDs, and MDe) and two kernel 366 methods (GBLUP and GK) to rubber tree data sets showed that models with the nonlinear GK 367 had slightly higher PAs than the models with the linear GBLUP kernel. According to Gianola et 368 al. (2014), the GK has a better predictive ability and a more flexible structure than the GBLUP, 369 and the GK can capture nonadditive effects between markers. and production of the rubber tree and contribute to a large amount of variability in the behavior 410 of cultivars (Ortolani et al., 1996). 411 In rubber trees, the time required to complete a breeding cycle and recommend a clone for 412 commercial production can span multiple decades and is mainly divided into three selection 413 stages. First, the aim is to obtain progenies by controlled or open pollination to establish 414 nurseries. At two and a half years, based on yield evaluations performed by early testing of yield, 415 vigor, and tolerance to diseases, the breeding plants are selected and cloned for testing at a small 416 scale. In this second stage of the selection cycle, after the first two years of tapping, promising 417 clones are multiplied and subsequently evaluated in large-scale or regional trials . This last stage  418  usually takes some 12 to 15 years, until it is possible to recommend a clone for large-scale  419 cropping. It, therefore, takes approximately 30 years to complete the breeding cycle, from 420 controlled pollination to final cultivar recommendation (Gonçalves and Fontes, 2012). The use of 421 GS could dramatically reduce the time required for completion of a cycle of genetic 422 improvement by eliminating progeny phenotypic testing aimed at selecting the best individuals 423 (replaced by GS), significantly accelerating the genetic gain relative to that obtained by classical 424 breeding. Another advantage of GS is that more candidate genotypes are generated; therefore, 425 the population size for selection as well. All of these candidates are genotyped, and those with 426 the best-predicted test cross values are evaluated in the field, which can be regarded as an 427 indirect selection. 428 With the rapid advances in and declining costs of genotyping methods, balanced against the 429 overall costs of managing large progeny trials and the potential for increased gains per unit time, 430 our cautiously optimistic expectation is that GS has excellent potential to be implemented in 431 rubber tree breeding programs. However, further studies examining populations with a different 432 structure (which were not assessed in this initial work) are necessary before recommending GS 433 for operational implementation in tree breeding programs. 434

Conflict of interest 435
The authors have no conflicts of interest to declare. 436    Tables 596 Table 1. Phenotypic variation: heritability (H) of height (AP) and diameter (DAP), genotype x 597 environment interaction (GxE), residual (R) and genetic main effect (G) in the low-water (LW) 598

Author contributions
and well-watered (WW) environments considered together and alone, with p<.01 indicated by ** 599