Multienvironment Models Increase Prediction Accuracy of Complex Traits in Advanced Breeding Lines of Rice

Genotype  ́ environment interaction (G  ́ E) is the differential response of genotypes in different environments and represents a major challenge for breeders. Genotype  ́ year-interaction (G  ́ Y) is a relevant component of G  ́ E, and accounting for it is an important strategy for identifying lines with stable and superior performance across years. In this study, we compared the prediction accuracy of modeling G ́ Y using covariance structures that differ in their ability to accommodate correlation among environments. We present the use of these approaches in two different rice (Oryza sativa L.) breeding populations (indica and tropical japonica) for predicting grain yield, plant height, and three milling quality traits—milling yield, head rice percentage, and grain chalkiness—under different cross-validation (CV) scenarios. We also compared model performance in the context of global predictions (i.e., predictions across years). Most of the benefits of multienvironment models come from modeling genetic correlations between environments when predicting performance of lines that have been tested in some environments but not others (CV2). For predicting the performance of newly developed lines (CV1), modeling between environment correlations has no effect compared with considering environments independently. Response to selection of multienvironment models when modeling covariance structures that accommodate covariances between environments was always beneficial when predicting the performance of lines across years. We also show that, for some traits, high prediction accuracies can be obtained in untested years, which is important for resource allocation in small breeding programs. E. Monteverde and S. McCouch, Plant Breeding and Genetics Section, School of Integrative Plant Science, Cornell Univ., Ithaca NY 14853; J.E. Rosas, Dep. of Statistics, College of Agriculture, Univ. de la República, Garzón 780, Montevideo, Uruguay; J.E. Rosas, P. Blanco, and F. Pérez de Vida, National Rice Research Program, National Institute of Agricultural Research (INIA), INIA Treinta y Tres, Villa Sara 33000, Uruguay; V. Bonnecarrère, Biotechnology Unit, National Institute of Agricultural Research (INIA), Estación Experimental Wilson Ferreira Aldunate, Rincón del Colorado 90200, Uruguay; G. Quero, Dep. of Plant Biology, College of Agriculture, Univ. de la República, Garzón 809, Montevideo, Uruguay; L. Gutiérrez, Dep. of Agronomy, Univ. of Wisconsin–Madison, 1575 Linden Dr., Madison, WI 53706. Received 19 Sept. 2017. Accepted 8 May 2018. *Corresponding author (srm4@cornell.edu). Assigned to Associate Editor Aaron Lorenz. Abbreviations: BLUE, best linear unbiased estimator; BLUP, best linear unbiased prediction; CV, cross-validation; GBLUP, genomic best linear unbiased prediction; GBS, genotyping by sequencing; GC, grain chalkiness; GK, Gaussian kernel; GS, genomic selection; GY, grain yield; G  ́ E, genotype  ́ environment interaction; G  ́ Y, genotype ́ year interaction; MTM, multiple-trait model; MY, milling yield; PH, plant height; PHR, head rice percentage; SNP, Single-nucleotide polymorphism. Published in Crop Sci. 58:1519–1530 (2018). doi: 10.2135/cropsci2017.09.0564 © Crop Science Society of America | 5585 Guilford Rd., Madison, WI 53711 USA This is an open access article distributed under the CC BY license (https:// creativecommons.org/licenses/by/4.0/). Published June 21, 2018

production has to be concomitant with an increase in grain quality.
Common to many traits of interest to breeders, yield and grain quality traits are controlled by many genes.The most promising way to improve yield and quality traits is through the application of new genomic tools and statistical approaches such as genomic selection (GS) (Meuwissen et al., 2001), which is designed to enhance the rate of genetic gain for complex traits (Heffner et al., 2009;Jannink et al., 2010).Initially, GS methods were focused on single-trait, single-environment analyses, but one major challenge for breeders is the differential response of genotypes in different environments, known as genotype ´ environment interaction (G ´ E).The G ´ E can affect trait heritability and line ranking over environments, frequently affecting decision making.For this reason, GS approaches capable of modeling G ´ E have increasingly gained popularity.
Several models accounting for G ´ E have been proposed to characterize the mean response of genotypes across environments, with the objective of developing locally adapted genotypes.Regression analysis (Yates and Cochran, 1938;Finlay and Wilkinson, 1963) measures the relative genotypic performance and stability across multiple locations.Multiplicative models combine univariate and multivariate approaches for reducing the dimensions of data and facilitating interpretation and selection decisions (Gauch, 1988;Smith et al., 2005).The most commonly used multiplicative model is the additive main effect and multiplicative interaction model (AMMI; Gauch, 1992) that combines ANOVA for estimating genotype and environment main effects and singular value decomposition to partition the G×E component.The use of mixed model approaches was later introduced to improve flexibility and to allow different correlation structures among environments (Piepho, 1998;Piepho and Möring, 2005;Burgueño et al., 2008Burgueño et al., , 2012)).
Modeling covariance matrices to account for G×E allows borrowing information among correlated environments.Curnow (1988) showed that using correlated information about treatment effects was beneficial for selecting the best treatment.Atlin et al. (2000) took up this idea to study the response to selection in a subdivided target region.Piepho and Möring (2005) applied this approach on a best linear unbiased prediction (BLUP) framework and showed that when there are different target regions in a breeding program, models that fit a genetic correlation across subregions showed better performance than only using data from the particular target region or using a global average value for genotypic performance.
Mixed models that allow the incorporation of a genetic covariance matrix calculated either from pedigree or molecular marker data, rather than assuming independence among genotypes, also improve the estimation of genetic effects (VanRaden, 2008).The advantage of using genetic covariance matrices in G ´ E mixed models is that the model relates genotypes across locations even if lines are not present in all locations.Burgueño et al. (2012) were the first to combine genetic and environmental covariance matrices using different covariance structures to model the environmental component.They showed that modeling both genetic and environmental covariances improved predictions in the target environment for genotypes that have already been evaluated in correlated environments.
If we consider G ´ E in the context of genotype ´ year interactions (G ´ Y), identifying lines with stable and superior performance across years is an important challenge for breeders.In this context, breeders are not only interested in predicting performance in future seasons (local predictions), pointing to the importance of repeatable G ´ Y interactions, but also in predicting rankings across years, taking G ´ Y as part of the noise (global predictions).The ability of models to predict G ´ Y in this context is crucial for crop improvement in current and future climate change scenarios.This is likely to be especially beneficial for rice breeding programs that are focused on high yield and quality, given that both components are highly affected by environmental conditions (Cameron and Wang, 2005;Chen et al., 2012;Lyman et al., 2013).
The objective of our study was to compare the effect on prediction accuracy of different types of prediction models for G ´ Y interaction.The models compared in this study differ in whether they allow borrowing of information across years or not, for both local and global prediction scenarios.We also compared the effect of using two different kernel regression approaches to construct the genetic covariance matrix: a linear kernel (genomic BLUP [GBLUP]), which assumes that allelic effects are additive, and a nonlinear reproducing kernel Hilbert space approach, which takes into account nonadditive as well as additive effects between genes.We present the use of these approaches in two different rice-breeding populations (indica and tropical japonica), which have been evaluated for five different traits across 3 yr in a single location.

Phenotypic Data: Experimental Design
The data used in this study belong to the National Institute of Agricultural Research (INIA-Uruguay) and consists of 3 yr of phenotypic data from field trials (2010,2011,2012) collected from 309 tropical japonica and 327 indica elite rice lines in one location, Paso de la Laguna Experimental Station (UEPL), Treinta y Tres, Uruguay (33°15¢ S, 54°25¢ W).These two breeding populations were originally developed for a genomewide association study project for grain quality traits, and results from that study have been reported by Quero et al. (2018).
Each year, the tropical japonica and indica populations were planted independently in replicated trials in six-row plots using an augmented randomized complete block design with two is the vector of the response variable of the BLUE values for the jth environment; m = (1 n1 ¢m 1 , …,1 nj ¢m j , …, 1 nm ¢m m )¢, where m represents the vector of intercepts for each environment and 1 nj is a vector of ones of order n j ; ( ) is the random vector of genetic values with u j ?N(0,S), where S is a genotypic covariance matrix between environments; and e = (e 1 , …, e j , …, e m )¢ is the vector of random residuals with ( ) . In this study, R takes a diagonal form, where each environment has its own residual variance.I n is an identity matrix of order n, and Ä is the Kronecker product.
This general mixed model can be used to predict genotypes that have not been evaluated in either a particular environment or in any environment, via the S covariance matrix, that allows borrowing information from evaluated lines and/or environments.In this study, we discuss two different models for the structure of S (diagonal and unstructured) and their implications for predicting genotype performance under different scenarios.

Models for S
The genetic covariance matrix S can be decomposed into a genomic related matrix K and an environment related matrix U E .When the number of individuals is the same in each environment (n j = n j¢ = n) as is the case in this study, , where Ä is the Kronecker product.
There are several structures for modeling matrix U E in mixed models (Burgueño et al., 2012;Malosetti et al., 2016).In this study, we will use two contrasting structures for the U E matrix that differ as to whether the information between environments can be borrowed or not.The diagonal (U DIAG ) structure fits a separate variance component per environment assuming no genetic correlations between environments.The diagonal values of the diagonal matrix are thus variance components for each environment, with off-diagonal elements equal to zero.
The most general way to model covariances between environments is the unstructured m ´ m covariance matrix (U UN ): where the jth diagonal element of the variance-covariance matrix is the additive genetic variance ( 2 u m s ) within the jth environment, and the off-diagonal element ( u u j j¢ s ) is the genetic covariance between environments j and j¢.

Genomic Kernels
The matrix K is a genetic relationship matrix that can be derived either from pedigree information or from molecular marker data.In this study, we compared the effect of two different kernels (K) on genotype predictions.Model 1 uses the linear kernel, also known as GBLUP, where K = (XX¢/p), where X represents the matrix of p centered and standardized molecular markers (originally coded as 0, 1, or 2) in the jth environment.
Model 2 uses the Gaussian kernel (GK) as covariance structure.The GK used in this work was ( ) or three replications and included two Uruguayan cultivars as checks for each population (El Paso 144 and INIA Olimar for the indica dataset, and EEA-404 and INIA Tacuari for the tropical japonica dataset, http://www.inase.org.uy/Sitio/RegistroNa-cionalCultivares/Default.aspx).Trials were conducted under irrigated conditions using appropriate pest and weed control.
The agronomic traits analyzed each year were grain yield (GY of paddy rice in kg ha −1 ) and plant height (PH measured in cm from the soil surface to the tip of the flag leaf ).The grain quality traits measured were yield after milling (MY measured in g, as the weight of grain recovered after milling divided by the weight of rough rice before milling, using a 100-g sample of rough rice), percentage of head rice recovery (PHR measured in g, as the weight of whole milled kernels, divided by the weight of rough rice, using a 100-g sample of rough rice), and the percentage of grain chalkiness (GC, measured as % of chalky kernels in a subsample of 50 g of total milled rice, where the area of chalk core-white back or white belly-was larger than half the kernel area based on visual inspection) (Supplemental Table S1).The grain quality traits were measured as described in Quero et al. (2018) and in Siebenmorgen et al. (2012).
Genotype data for both datasets were obtained using genotyping by sequencing (GBS).Single-nucleotide polymorphism (SNP) calling was performed using the TASSEL 3.0 GBS pipeline (Bradbury et al., 2007), and SNPs were aligned to the Michigan State University Nipponbare reference genome version 7.0 using Bowtie 2 (Langmead and Salzberg, 2012).Imputation of missing data was performed with the FILLIN algorithm implemented in TASSEL 5.0 (Swarts et al., 2014) for both datasets separately.The GBS datasets were filtered to retain markers with <50% missing data after imputation, and a minor allele frequency <0.05, as reported by Quero et al. (2018).The final indica dataset contained 92,430 markers and the tropical japonica dataset had 44,598 markers.The final datasets was transformed to numeric coding (−1, 0, and 1 for Class I homozygotes, heterozygotes, and Class II homozygotes respectively) to facilitate statistical analysis.

Statistical Models
The best linear unbiased estimators (BLUEs) were calculated on a per line basis for each year and trait separately.The model used to calculate BLUEs for each year is where y ijkl is the trait score, m is the overall mean, b i is the random block effect with where s 2 b is the block variance, g j is the genotypic effect, r k(i) and c l(i) are the random row and column effects nested within blocks with

Prediction Models
The general model for the G ´ E interaction is given by the following equation that fits the data for n lines (i = 1, …, n) in m environments ( j = 1, …, m) in the following way:

Model Implementation
All the models tested in this work were fitted using R (R Core Team, 2017).The multienvironment models were fitted using the multiple-trait model (MTM) software developed by de los Campos and Grüneberg (2016) that also uses a Bayesian approach and assumes that K is the same in all environments.This package assumes that the prior distribution for u is multivariate normal with mean zero and a variance-covariance matrix The intercepts are assigned a flat prior, and the prior distribution for U E is an inverse Wishart p(U E |S 0 , df 0 ) = W −1 (S 0 , df 0 ), where the parameters are a scale matrix (S 0 ) equal to an identity matrix of order m (number of environments), and degrees of freedom (df 0 = m).The error variances for each environment ( ) are the diagonal elements of R and are assigned a scaled inverse c 2 distribution with df and scaled factor equal to 1.
The MTM package uses Markov chain Monte Carlo (MCMC) and the Gibbs sampler to fit the models.In this work, we used 55,000 iterations with a burn-in of 5000 and a thinning of five for the Gibbs sampler.For the multienvironment models, the posterior variance-covariance matrices between environments were also calculated for each model using the full datasets.An example of the R code to run these models is provided in the supplemental material.

Assessing Prediction Accuracy
Prediction accuracies for each trait for both the U DIAG and U UN models were assessed using 50 training-validation random partitions.We used two different cross-validation (CV) designs (Burgueño et al., 2012).The CV1 mimics the situation when newly developed genotypes have not been tested in any environments; in this scenario, 30% of the lines were missing in all three environments and were treated as unknown.In this case, predictions for the missing lines are based on the phenotypic records of the other lines that were tested.The CV2 design mimics the situation where some lines were evaluated in some environments but were missing in others.For the two CV scenarios, we assigned individuals to the training and validation sets using the same procedure as Lopez-Cruz et al. (2015).
In each random partition, the Pearson's correlations (r) between the predicted and the observed values were computed.As in Malosetti et al. (2016), and to comply with the normality assumption, the Fisher's z transformation was used, where Simple factorial ANOVA and Tukey's tests were performed at a significance level a = 0.05 over the z-transformed values to determine significant differences in prediction accuracies between the two genetic kernels (GBLUP and GK) and the two covariance structures (U DIAG and U UN ) within each subpopulation (indica and tropical japonica) and each year.Mean prediction accuracies were presented using the original scale after back transformation: Another prediction problem studied here was that of predicting future seasons and was denoted as "leave-1-yr-out."This prediction was performed by using 2 yr in the training population to predict a third year where no phenotypic data were collected.

Response to Selection under Global Prediction Scenario
When dealing with G ´ Y, it is interesting to explore the line performance across years (global prediction).For the global prediction scenario, we consider the response to selection as a way of comparing the predictive ability of the two covariance matrices (U DIAG and U UN ) and the two genetic kernels (GBLUP and GK), as reported by Piepho and Möring (2005).
In the global scenario, the genotypic value of the ith genotype accounting for all the evaluated years is being some environment weighting parameter such as the relative growing areas in the m subregion, as discussed in Piepho and Möring (2005), or the inverse of the heritability, as evaluated in the present study.
As in Piepho and Möring (2005), we regard the estimator m to be an indirect trait, in the special case that w¢ = v¢W, where L i = v¢g i and assuming i g  is the BLUP of genetic effects with m and W = U E (U E + S e ) −1 .We formulate the selection response as a correlated response to selection as where r g is the genetic correlation between L i and i L  and can be calculated as i is the selection intensity, and h is the square root of the heritability (Falconer and Mackay, 1996) that can be calculated per Piepho and Möring (2005) as

w U w w U w
To evaluate different methods, we consider the ratio of R values, since both the selection intensity (i) and the ( ) terms cancel out.In this way, the ratio R UN /R DIAG is calculated as where r UN , r DIAG , h UN, and h DIAG are the r g and h values calculated by substituting U E with U UN and U DIAG , respectively.
Prediction Accuracies for U DIAG and U UN Models under Different Cross-Validation Scenarios For the indica dataset, among all the prediction models tested with the CV1 and CV2 methods, the relative advantage of using the U UN covariance matrix over the U DIAG was evident for most traits only when the CV2 method was used, although GY was an exception (Fig. 2).There were no significant differences between CV1 and CV2 when the U DIAG was used.As expected, mean prediction

Descriptive Statistics
Boxplots of the five traits (GY, MY, PHR, GC, and PH) in each of the years (2010-2012) for the indica and the tropical japonica datasets are shown in Fig. 1.
For the indica dataset, the trait distributions for each year are similar, although we see a slight increase for GY and decrease for PH in 2012 (Fig. 1A).For the tropical japonica dataset, 2011 was a bad year for GY, showing the lowest value for this trait.On the other hand, 2012 was the worst year in terms of grain quality and PH, showing the lowest values for MY, PHR, and PH and the highest for GC (Fig. 1B).The ANOVA analyses performed on each trait and each population on both datasets confirmed significant G ´ Y interaction (Supplemental Tables S2 and S3).
Trait heritabilities by year were intermediate to high for both datasets.For the indica dataset, the lowest values were MY in 2012 (0.42) and GY in 2010 (0.46), whereas for tropical japonica, the lowest values were GY in 2010 (0.43) and GC in 2011 (0.41) (Table 1).accuracies for CV1 were lower than those obtained with CV2 for U UN , since this scenario benefits from borrowing information across environments (Fig. 2).
Phenotypic correlations among environments for both datasets were all moderate to high and positive (Supplemental Table S4).In general, PH showed the highest correlation for both indica and tropical japonica, whereas GY for indica showed the lowest correlations.Estimated genetic variance within environments (diagonal) and covariances and correlations between environments (offdiagonals) were higher for GK than for GBLUP for all traits, except GC (Supplemental Table S5).Among all traits, GY showed the lowest gains from modeling the U UN covariance structure under CV2, compared with U DIAG (Fig. 2).Grain yield also showed generally low variance-covariance values and correlations compared with other traits (Supplemental Tables S4 and  S5).On the other hand, PH was the trait that showed the highest gain from U UN vs. U DIAG for CV2 (Fig. 2) and the highest genetic variance-covariance values and correlations among all traits (Supplemental Tables S4 and S5).
For the tropical japonica dataset, the CV2-GK method when using the U UN matrix was the best performing method, with the exception of GY (when predicting 2010 and 2011), MY (when predicting 2012), and GC (when predicting 2012) (Fig. 3).For this dataset, PH was again the trait that showed the highest gains from prediction when modeling covariance between environments under CV2.There were no significant differences when predicting with the GBLUP vs. GK, except for PHR where the CV2-GK method was the best when using the unstructured covariance matrix (Fig. 3).Supplemental Table S6 shows that the elements of U E for the GK are larger than those for the GBLUP method, and the opposite occurred with the diagonal elements of the R matrices.
The results of the performance of the models in untested environments (leave-1-yr-out) are shown in Table 2.The advantage of modeling U UN over U DIAG is more evident for those traits that show a higher correlation between environments, such as MY, PHR, and PH in indica, and PHR and PH in tropical japonica.In all cases prediction accuracies for U UN were either equal or higher than those obtained with U DIAG .

Estimation for Global Adaptation in the indica and Tropical japonica Datasets
To study the gain from borrowing information between environments in the global adaptation scenario, we calculated the ratio of the correlated response to selection obtained from modeling U UN and U DIAG .We assumed two different weighting scenarios: (i) by growing area, with subregions of equal size where v = (v 1 , v 2 , v 3 ) with v r = 1/3, and (ii) by using the inverse of the heritability, where . In both cases, w = v¢W, where W is either U UN or U DIAG and R UN and R DIAG are the responses calculated with both matrices, respectively.Ratios R UN /R DIAG are shown in Table 3.
The results show that, for both weighting scenarios, modeling a covariance structure that allows for borrowing information across environments is always beneficial, in both datasets (Table 3).The traits that showed highest R UN / R DIAG ratios in the indica dataset were GY, MY, and PH, with gains >50%, whereas in tropical japonica, the traits with highest R UN /R DIAG ratios were GY, MY, and PHR.
For both datasets, the relative gains from borrowing information across years are higher for the GBLUP kernel than when the GK kernel is used (Table 3).

DISCUSSION
One of the major challenges in plant breeding is the differential response of genotypes in different environments, known as G ´ E interaction.Finding lines with predictable G ´ E across seasons (G ´ Y) is becoming an increasingly pressing issue due to the highly variable and unpredictable weather conditions in current and future climate change scenarios.Several approaches for modeling G ´ E have been proposed in the literature (Burgueño et al., 2012;Jarquín et al., 2014;Heslot et al., 2014;Lopez-Cruz et al., 2015), and all confirm the relative benefits of using multienvironment models in comparison with single-environment analyses.
In this work, we used a multienvironment modeling approach for GS where the genetic effects of rice lines were calculated using the Kronecker product of genetic variance-covariance matrices between environments (U E ) and marker-genomic relationship matrices calculated with two different kernel methods, GBLUP and GK.We compared the effect of using two different covariance structures for the U E matrix: (i) a covariance structure that fits a separate variance component per environment and does not allow borrowing information across environments (U DIAG ), and (ii) a covariance structure that has a variance component for each environment and a separate covariance parameter for each pair of environments (U UN ).
These prediction models were implemented as reported in previous studies (Burgueño et al., 2012;Malosetti et al., 2016;Cuevas et al., 2017).Prediction accuracies from both multienvironment models were calculated for five traits  with different genetic architectures and medium to high heritabilities: GY, three different grain quality traits, and PH, all measured in three different years (2010, 2011, and 2012).The two datasets used in this study were an indica and a tropical japonica rice breeding population belonging to the Uruguayan National Rice Breeding Program.Population structure and genome-wide association study analyses on these same two populations and the three grain quality traits analyzed in this study were previously reported (Quero et al., 2018).In our dataset, the environments represent 3 yr of evaluation (2010, 2011, and 2012), so G ´ E is considered as G ´ Y.In this context, we performed both local predictions by predicting the traits for each of the tested years, and global predictions by predicting line performance across years, taking G ´ Y as part of the noise.For local prediction, we considered three different CV scenarios: (i) CV1, where we predicted the performance of untested lines on tested environments, (ii) CV2, where we predicted the performance of tested lines in tested environments, and (iii) leave-1-yr-out, where we predicted tested lines in untested environments.For global prediction, we compared the correlated response to selection calculated from models using both covariance structures (U UN and U DIAG ), based on the method reported by Piepho and Möring (2005).
Genomic selection models capable of accounting for multienvironment data have been demonstrated to increase prediction accuracies relative to single-environment analyses (Lopez-Cruz et al., 2015;Malosetti et al., 2016;Lado et al., 2016;Saint Pierre et al., 2016;Bandeira e Sousa et al., 2017;Cuevas et al., 2017).In our study, we found that multienvironment models that allow borrowing information across environments were superior to or not different from those that do not exploit genetic covariance between environments.We also found that the gains in prediction accuracy from models that exploit the genetic covariance between environments depend on the CV method used.This was expected, because when using phenotypic information from tested genotypes in already tested environments (CV2), we are better exploiting the correlation between environments.On the other hand, under CV1, the information between environments flows in a more indirect fashion through the kinship relatedness matrix only.These results were in accordance with previous studies that used similar CV designs (Burgueño et al., 2012;Jarquín et al., 2014;Zhang et al., 2014;Crossa et al., 2016;Malosetti et al., 2016;Saint Pierre et al., 2016;Bandeira e Sousa et al., 2017).For this reason, it is always beneficial to use as much data as possible to train our models, and to partition our resources by phenotyping as many different genotypes as possible, and spreading them across the tested environments, rather than fully phenotyping only a part of the population, as also suggested by Malosetti et al. (2016).This could be helpful for resource allocation in a small breeding program like the Uruguayan breeding program.
Previous studies in wheat (Triticum aestivum L.), maize (Zea mays L.), and rice have demonstrated that multienvironment GS models can give improved prediction accuracies by borrowing information across highly correlated environments (Crossa et al., 2014;Lado et al., 2016;Saint Pierre et al., 2016;Spindel et al., 2016).For these multienvironment models, the genetic correlations between environments depend in part on the off-diagonal values of matrix U UN .When these values are close to zero, the U UN matrix tends to U DIAG , producing prediction accuracies similar to those obtained by using U DIAG models.When the off-diagonal values of U UN are moderate to high (either positive or negative), the correlations between environments will be far from zero, allowing the borrowing of information across environments.In this case, the U UN models for those environments will perform better than the U DIAG models under CV2.
In our study, the phenotypic and genotypic correlations between environments were all positive and moderate to high.This led to higher prediction accuracies in all environments using the U UN models in comparison to U DIAG models when using CV method CV2.Plant height was the trait that showed the highest genetic and phenotypic correlations between environments for the two datasets.On the other hand, traits with lower genetic covariance between environments (off-diagonal values of U UN ) showed similar prediction accuracies between U DIAG and U UN models.This was the case for GY in indica, where genetic covariances were the lowest, leading to comparable prediction accuracies between U DIAG and U UN models.The accurate estimation of the genetic covariance structure between correlated environments is critical to increase prediction accuracy with multienvironment models.Unstructured covariance matrices have m(m + 1)/2 unique parameters, where m is the number of environments.When m is large, estimation of an unstructured covariance matrix may be unfeasible due to convergence problems.Also, the estimation of multiple covariance parameters poses the need for sufficiently large datasets to support accurate estimation of the numerous parameters involved (Denis et al., 1997;Meyer and Kirkpatrick, 2008;Meyer, 2009;Elias et al., 2016).For our two balanced datasets with only three environments, fitting an unstructured covariance matrix did not pose any convergence issues and worked better than a diagonal matrix model.However, as the number of environments increases, fitting a more parsimonious covariance structure (such as factor analytic) should be considered (Piepho, 1998;Smith et al., 2001).
Prediction of untested environments is a highly relevant topic for breeders because it allows predicting which lines will be the best performers in future testing.Another good reason to use CV designs such as leave-1-yr-out is to avoid the overlap between training and validating populations.Predicting for already tested environments, such as in CV1 and CV2 scenarios, assures the overlap between the testing and validation populations along the environmental and genotype dimensions.This could lead to overfitting, and thus inflation of the prediction accuracies.Here, we evaluated the predictions for different traits in untested years (leave-1-yrout), where the training set included the same genotypes as the testing set, but no phenotypic information for the year tested, thus avoiding the overlap of information along the environment dimension.Our analysis revealed that years where no genotypes have been evaluated could be predicted with good accuracy (>40% in most cases) when using the U UN matrix.Previous studies have remarked that predicting environments that have not yet been tested could be a difficult task, pointing to the need to establish a link between the tested and untested environments.A good way to establish this link could be using environmental covariates (Malosetti et al., 2016;Saint Pierre et al., 2016); however, positive correlation between environments is still an important factor for achieving good prediction accuracies in unobserved environments, as shown in this study.Jarquín et al. (2016) optimized training sets for genomic prediction of soybean [Glycine max (L.) Merr.] accessions using designs similar to our leave-one-site-out with no environmental covariates and also showed high prediction accuracy for percentage protein and GY.High correlation among sites could be also the result of stable climate conditions, which tend to be more common in temperate areas than in the tropics, which experience greater seasonal variation, as well as more disruptive weather events (e.g., typhoons).Spindel and McCouch (2016) stressed the importance of highly correlated phenotypic and environmental data to achieve higher prediction accuracies and improve breeding outcomes.As long as our environments remain highly correlated, with relatively stable weather conditions, and as long as the breeders do not introduce any novel genetic diversity into their populations, prediction accuracies will remain high, as shown in this study.Another factor that contributed to the prediction accuracies observed in this study were the moderate to high heritabilities found for our traits across years in both datasets.Good agronomy practices and experimental design, as long as the irrigated rice paddy system homogenizes field conditions, improve heritability in rice compared with other crops such as wheat or maize, where water availability is an important variable affecting trait expression.
Besides prediction for specific years and environments (local prediction), another relevant scenario is to predict genotype rankings across years (global prediction).This is particularly relevant given the unpredictable weather patterns related to years.In this work, we evaluated the response to selection of multienvironment models when information is borrowed across environments using the R UN /R DIAG ratio.To calculate the responses, we used a method previously reported by Piepho and Möring (2005) that involves a weighting scheme based on BLUP.Our conceptual framework for this global scenario assumes an extensive target region (i.e., the Uruguayan rice breeding program) that is subdivided into subregions (years), and the objective is to predict the best performing genotypes across years considering G ´ Y as part of the noise.We concluded that modeling covariance structures that accommodate covariances between environments is always beneficial (or at least has no penalty) when predicting the performance of lines across years, and for this reason, it is always useful to consider all data.We acknowledge that our dataset is small, as it includes data from one location and only 3 yr of testing, and that this may limit the conclusions we can draw.However, the increase in prediction accuracy observed under CV2, and the increased response to selection when using unstructured covariance matrices strongly suggests that the addition of more information should improve the prediction ability of our U UN models.A more comprehensive analysis with more years of testing and the incorporation of environmental covariates is the subject of current research.
Various studies have documented the benefits of using GK in genomic prediction in both single-and multi-environment settings in crops such as maize, wheat, and rice (Gianola et al., 2014;Iwata et al., 2015;Onogi et al., 2015;Pérez-Elizalde et al., 2015;Spindel et al., 2015;Cuevas et al., 2016;Bandeira e Sousa et al., 2017;Cuevas et al., 2017).However, in our study, prediction accuracies obtained with GK were in general not significantly different than those obtained with a linear kernel (GBLUP), with a few exceptions detailed in the results section.
Genomic selection can shorten breeding cycle length and increase gains per unit time by replacing time-intensive phenotypic evaluation of complex traits with genomic estimated breeding values.In temperate crops, GS can be exploited to accelerate gain from selection through the use of off-season nurseries, where phenotyping is difficult and costly.This could be beneficial for temperate rice growing countries, like Uruguay, where only one rice cycle per year is possible.There are several points where GS could be implemented in a standard pedigree selection-based rice breeding program to either avoid a cycle of phenotyping or to help to decide which lines should be advanced (or discarded) for the next generation.In general, a GS step before conducting multienvironment trails is suggested, and some studies propose the incorporation of phenotypic data from preliminary yield trials into the GS framework (Endelman et al., 2014;Michel et al., 2017).
Our results show that for the germplasm and environmental conditions used in this work, most of the benefits of multienvironment models come from modeling genetic correlations between environments under CV2.For predicting the performance of newly developed lines (CV1), modeling between-environment correlations has no effect compared with considering environments independently.However, our study showed that for some traits, high prediction accuracies can be obtained in untested years, which is important for resource allocation in small breeding programs.Additional research incorporating the use of environmental covariates to model future seasons and bringing historical data to model climate response patterns are promising steps to improve the prediction accuracies as a part of our effort to accelerate the rate of genetic gain for both yield and grain quality in rice breeding programs in temperate areas.
error variance.Broad-sense heritability for each environment was calculated on a per line basis as() variance among rice lines and r is the number of replicates.
the Euclidean distance between individuals i and i¢ based on molecular markers, and h is the bandwidth parameter that controls the rate of decay of the values of K. FollowingCrossa et al. (2010)

Table 2 .
Correlation between the observed and predicted values of the leave-1-yr-out prediction problem for grain yield (GY), milling yield (MY), head rice percentage (PHR), grain chalkiness (GC), and plant height (PH) for indica and tropical japonica Uruguayan breeding populations.

Table 3 .
Response to selection for global adaptation scenario for grain yield (GY), milling yield (MY), head rice percentage (PHR), grain chalkiness (GC), and plant height (PH) for indica and tropical japonica Uruguayan breeding populations.Response was measured as the ratio R UN /R DIAG where R UN is the response to selection when information between environments is borrowed (unstructured matrix), and R DIAG is the response to selection when information between environments is not borrowed (diagonal matrix).