Improving accuracy of genomic prediction in maize single-crosses through different kernels and reducing the marker dataset

Improving accuracy of genomic prediction in maize single-crosses through different kernels and reducing the marker dataset In plant breeding, genomic prediction (GP) may be an efficient tool to increase the accuracy of selecting genotypes, mainly, under multi-environments trials. This approach has the advantage to increase genetic gains of complex traits and reduce costs. However, strategies are needed to increase the accuracy and reduce the bias of genomic estimated breeding values. In this context, the objectives were: i) to compare two strategies to obtain markers subsets based on marker effect regarding their impact on the prediction accuracy of genome selection; and, ii) to compare the accuracy of four GP methods including genotype × environment interaction and two kernels (GBLUP and Gaussian). We used a rice diversity panel (RICE) and two maize datasets (HEL and USP). These were evaluated for grain yield and plant height. Overall, the prediction accuracy and relative efficiency of genomic selection were increased using markers subsets, which has the potential for build fixed arrays and reduce costs with genotyping. Furthermore, using Gaussian kernel and the including G×E effect, there is an increase in the accuracy of the genomic prediction models.


INTRODUCTION
Recent technologies for large-scale genotyping and the development of statistical tools have been used to increase selection in breeding programs. Molecular marker-assisted breeding, including marker-assisted selection (MAS) and genomic prediction (GP), in combination with high-throughput and precise phenotyping, can significantly accelerate the development of new varieties (Xu et al. 2012). MAS is a breeding approach based on significant associations of markers and target traits where only the significant markers are used for selection (Xu et al. 2012). However, quantitative traits are affected by many genes, and the limited of MAS is due to the genetic variance proportion explained by the quantitative trait loci (QTL). In GP, breeding values of progenies can be predicted by regressing phenotypic values on all available markers ).
According to , if a marker is in linkage disequilibrium with a QTL, same markers alleles will be correlated with positive effects on the quantitative trait, across all families, and can be used without the need to establish linkage phase in each family. In addition, GP is based on two phases, training and selection (Jonas and de Koning 2013). In training phase markers effects will be estimated based on a phenotypic data of a training population and then tested in a validation population, in order to verify the prediction accuracy of the model (Newell and Jannink 2014). In selection phase, genomic estimated breeding values (GEBV's) of the segregating population are predicted, in which only the individuals are genotyped. GP utilizes the advantage of high genome coverage based on molecular markers, simultaneously adjusting the genetic effects for all the markers, without using statistical tests. In addition, it allows to predict the performance of new breeding parents in early generations and generate crossing and selections based on the prediction model (Jannink et al. 2010).
Genomic prediction was first proposed in animal breeding , and then applied to plant breeding (Crossa et al. 2010(Crossa et al. , 2011de los Campos et al. 2010). Several genomic prediction studies were applied to simulated and real plant breeding data Crossa et al. 2010Crossa et al. , 2011de los Campos et al. 2013;Perez-Rodriguez et al. 2013;. In general, these studies showed good prediction accuracies for complex traits. In plant breeding, the GP provides opportunities to increase the genetic gain of complex traits per unit of time and reduce costs Heffner et al. 2009;de Oliveira et al. 2012). In a wheat program, genomic prediction schemes have similar costs to phenotypic selection, but the potential of increasing the gain per unit time should be the real driver .
Several factors influence the accuracy of genomic prediction, such as genetic architecture, non-additive effects, models, the presence of G×E interaction, and marker density Spindel et al. 2015;Cuevas et al. 2016;Ma et al. 2016). Many statistical methods have been proposed to estimate the GEBVs in the training population. The most common method used is the genomic best linear unbiased prediction (GBLUP) model, which uses the genetic marker information for computing associations between individuals by the Genomic Relationship Matrix (GRM, G) (Habier et al. 2007). The linear GBLUP model capture additive effects between markers. However, complex traits are affected by non-additive effects due to dominance and genetic interactions (epistasis) and their interaction with environment.
Thus, non-parametric and semi-parametric methods modeling the relationship between the phenotype and markers in a GS context were proposed by Gianola et al. (2006). Non-parametric methods are able to capture small epistatic interactions without explicitly modeling them. An example of a semi-parametric method is the Reproducing Kernel Hilbert Space (RKHS), that uses a kernel function to convert the marker matrix into a set of distances between pairs of individuals. Therefore, reducing the size of the parametric space and thus can capture small complex interactions between genes Gianola et al. 2014). Thus, RKHS methods had potential to increase accuracy prediction (Perez-Rodriguez et al. 2013;. Another factor that influences the GP is the presence of genotype × environment interaction (G×E) effect. Genomic-enabled prediction models were originally developed for a single trait in a single environment. However, multi-environment plant breeding trials are routinely conducted to estimate and take advantage of genotype × environment interaction (G×E). Therefore, to implement GS strategies in plant breeding, G×E needs to be estimated, modeled, and predicted. Resende Júnior et al. (2012) and Windhausen et al. (2012) observed that environment had a high influence on the prediction accuracy of the model when the effects were validated in different environments from where they were estimated. It shows the impact of G×E in the prediction of genotype performance. Thus, Burgueño et al. (2012) extended the GBLUP methodology to incorporate and model G×E effects. The Bayesian model of  is another GBLUP extension that introduces main and interaction markers effects and environmental covariables via covariance structures. Heslot et al. (2014) proposed using crop modeling for assessing genomic G×E. Recently, Lopez-Cruz et al. (2015) proposed a GBLUP prediction model that explicitly models G×E and marker × environment interaction (M×E) where marker effects and genomic values are partitioned into components that are stable across environments (main effects) and others that are environment-specific (interactions). In general, studies have shown that modeling G×E can give substantial gains in prediction accuracy (Burgueño et al., 2012;Heslot et al., 2014;Crossa et al., 2016). Cuevas et al. (2016) compared methods that applied the G×E interaction GS model of Lopez-Cruz et al.
(2015) using a linear kernel (GBLUP) and a nonlinear Gaussian kernel and observed that methods using multi-environment G×E interaction models with Gaussian kernel showed higher prediction ability.
Current applications of GP are typically based on single nucleotide polymorphism (SNP) genotypes called from SNP array data. Several studies had shown that low-density arrays can improve the accuracy of genomic prediction (Moser et al. 2010;Spindel et al. 2015;Hoffstetter et al. 2016;Ma et al. 2016). It was demonstrated that high genotyping density does not always increase accuracy and markers subset sometimes outperform the entire dataset (Zhang et al. 2010;Ma et al. 2016). Standard sets of selected SNPs are a tendency for low-cost of genotyping in plant breeding. Therefore, it can develop alternative genotyping platforms less expensive and more repeatable using selected marker subsets.
With the promise of accurate predictions, the genomic prediction has been quickly implemented by several breeding programs. However, strategies are still needed to make the implementation of genomic prediction reducing costs with genotyping and with more accurate predictions. In this context, the objectives were: i) to compare two strategies to obtain markers subset based on marker effect regarding their impact on the prediction accuracy of genome selection; and ii) compare the accuracy of four GP methods including genotype × environment interaction and two kernels (GBLUP and Gaussian).

COSTS OF GENOMIC SELECTION
ABSTRACT Different methodologies to obtain subsets of markers have been proposed as an alternative to developing a fixed array of SNP, reducing costs and time with genotyping. Thus, the objectives were: (1) to compare two strategies to obtain markers subsets by markers effects: original and re-estimated effects; (2) to compare the prediction accuracy, bias and relative efficiency of a main genotypic effect (MGE) model fitted with a linear kernel Genomic Best Linear Unbiased Predictor, GBLUP (GB), and a nonlinear kernel Gaussian kernel (GK) using markers subsets. The two model-kernel combinations for each marker subsets were applied in a public dataset of rice diversity panel (RICE) and two hybrids maize datasets (HEL and USP). These were evaluated for grain yield and plant height. Decreasing the number of markers, higher estimates of prediction accuracy were obtained by MGE model with GB and GK kernels. Overall, the relative efficiency of genomic to phenotypic selection was increased using markers subsets by the re-estimated effects method. We concluded that from a high-density panel, it is possible to select the most informative markers to improve accuracy and build a low-cost array to implement genomic selection in breeding programs. The best strategy to obtain markers subsets is the re-estimation of the marker effect, which increases the accuracy and reduces the bias. Keywords: Gaussian kernel; GBLUP; SNP arrays; Relative efficiency

INTRODUCTION
With the advent of higher-density single-nucleotide polymorphisms (SNP) covering the whole genome of many plant species, breeding values can be predicted by regressing phenotypic values on all available markers Crossa et al. 2013). Afterward, genomic estimated breeding values (GEBV) are calculated based on these marker effects, which is the main output of genomic prediction (GP). Therefore, many statistical methods have been proposed to estimate the marker effects in the training population.
The first method was the Ridge Regression Best Linear Unbiased Prediction (RR-BLUP) ). This method is based on the assumption in that all marker effects are normally distributed and have equal variance . The most common method used, based on the same assumptions, and considered as equivalent is the genomic best linear unbiased prediction (GBLUP) model. This model uses the genetic marker information for computing associations between individuals through the Genomic Relationship Matrix (GRM, G) (Habier et al. 2007). Although it has shown satisfactory results, some authors proposed modifications in order to increase prediction accuracy. Besides the linear GBLUP kernel, semi-parametric methods to model the relationship between the phenotype and markers in a GP context were proposed by Gianola et al. (2006). One of them is the reproducing kernel Hilbert space (RKHS), a semi-parametric method that uses a kernel function to convert the marker matrix into a set of distances between pairs of individuals . Several studies have shown that using RKHS methods to obtain genetic relatedness between individuals has potential to increase accuracy prediction, mainly when genotype by environment interaction effect is added in the models (Cuevas et al. 2016;e Sousa et al. 2017).
Currently, high-density SNP chips provide the best coverage of genomes, but they tend to be expensive. On the other hand, genotyping by sequencing (GBS) is an affordable alternative.
However, it has shown relatively high rates of missing data and the quality of genotypes tends to be lower since it depends on the genome-wide sequence read depth (Gorjanc et al. 2015b). Thus, another approach is to develop low density and cost fixed arrays, which is highly efficient in smaller breeding programs (Spindel et al. 2015). The advantages of fixed arrays include robust allele calling, cost-effectiveness per data point and speed of genotyping turn-around (Thomson 2014a). Therefore, marker density has been investigated in some studies, which has shown divergent results regarding the use of low or higher density SNPs datasets (Habier et al. 2009;Weigel et al. 2009;Moser et al. 2010;Crossa et al. 2013;Perez-Rodriguez et al. 2013;Tayeh et al. 2015). Furthermore, higher genotyping density does not always increase accuracy and markers subset sometimes outperform the entire dataset (Zhang et al. 2010;Ma et al. 2016).
In this context, different methodologies to obtain subset of markers have been proposed, such as random (evenly distributed across the genome), marker effect, and markers significantly associated with the quantitative trait loci (QTL) (Resende et al. 2012;Spindel et al. 2015;Hoffstetter et al. 2016a). In general, the use of selected markers subset by their effects or positions has provided an efficient strategy to increase accuracy (Vazquez et al. 2010;Resende et al. 2012;Zhang et al. 2015). For instance, Tayeh et al. (2015) decreased the number of markers from 9824 to 2945 by retaining only a single marker per unique map position. They did not observe a reduction in prediction accuracies in any of the traits evaluated. Moreover, Ma et al. (2016) showed that the marker preselection based on haplotype block analyses is an interesting option to reduce the implementation cost of genomic selection.
In this sense, standard sets of selected SNPs are a tendency for low-cost of genotyping in plant breeding. However, for each targeted trait, it is necessary a specific SNP set, because the selected SNPs may vary with traits. In addition, the efficiency of a trait-specific low-density SNP chip depends critically on the linkage disequilibrium between the SNPs with large estimated effects and the true causative loci that affect the trait of interest (Wu et al. 2016). In this context, we compared two strategies to obtain markers subsets. The former is based on the original marker effect, obtained in the first estimate, using the whole dataset (original effect -ORI strategy). The latter is based on the re-estimated markers effects, that are re-estimated over the dataset reductions (reestimated effect -REE strategy). Thus, the objectives of our study were: (1) to compare two strategies to obtain markers subsets by markers effects: original and re-estimated effects; (2) to compare the prediction accuracy, bias and relative efficiency of a main genotypic effect (MGE) model fitted with a linear kernel Genomic Best Linear Unbiased Predictor, GBLUP (GB), and a nonlinear kernel Gaussian kernel (GK) using markers subsets.

Phenotypic data
For this study, we considered a rice and two maize datasets, using GY (grain yield, ton ha -1 ) and PH (plant height, cm): i) the first data is a public dataset available at Rice Diversity platform (http://www.ricediversity.org) (RICE). We used 270 elite breeding lines (F6-F7) from the International Rice Research Institute (IRRI) irrigated rice breeding program evaluated in a single location in Los Baños, Philippines during three years (2009 to 2011) in dry season (Spindel et al. 2015); ii) the second dataset is from Helix Sementes ® Company, São Paulo, Brazil (HEL).
HEL consisted of 452 maize hybrids obtained by crossing 111 inbred lines in a partial diallel. The experimental trial was carried out at five sites located in Southern, Southeast, and Midwest regions of Brazil, during the first growing season of 2014/15. The experimental design was a randomized block with two replicates per genotype and environment; iii) the third maize dataset is from University of São Paulo (USP). Data consists of 739 maize hybrids obtained by crossing 49 inbred lines in a partial diallel. The hybrids were evaluated at Piracicaba and Anhumas, São Paulo, Brazil, in 2016. The hybrids were evaluated using an augmented block design, with two commercial hybrids as checks. In both sites, the hybrids were evaluated under an ideal nitrogen (N) level, with 100 kg N ha -1 . Each plot had a range of 7 m, with 0.50 m spacing between rows and 0.33 m between plants; there was a phenotypic imbalance in both maize datasets.

Genotypic data
The RICE inbred lines were genotyped using about 73 K SNPs with GBS. The HEL and USP parent inbred lines were genotyped with an Affymetrix  Axiom  Maize Genotyping Array of 616 K SNPs (Unterseer et al. 2014). Standard quality control (QC) were applied to the data, removing markers with a Call Rate ≥ 0.95. The remaining missing data in lines were imputed with Synbreed package (Wimmer et al. 2015) using the algorithms from software Beagle 4.0 where is the response vector, is the replicate effect considered as fixed, is the vector of random effect of genotypes, where ~(0, 2 ), is the vector of fixed effects of environment; is the vector of random effects of Genotype by Environment interaction, where ~(0, 2 ) and is the vector of error, where ~(0, 2 ). The , , B and are incidence matrices. Environment effect for HEL and USP dataset is represented by site and for RICE dataset by year.
The significance of random effects was estimated using the Likelihood Ratio Test (LRT) (Gilmour et al. 2009). The variance components estimated for each model effect were used to calculate entry-mean based heritability ( 2 ), using means from each environment (considering each year or site as an environment).

Obtaining markers subsets
We used RR-BLUP (Random Regression -Best Linear Unbiased Predictor) method (Endelman 2011) to obtain markers subsets based on SNP effects: where is the vector of BLUPs, is intercept, is a matrix of markers, is a vector of marker effects, and is a residual effect. The following distribution were assumed ~(0, 2 ) and ~(0, 2 ).
Then, we applied two strategies to obtain the markers subsets, the reduction of the number of markers by i) their original effect (ORI strategy) and ii) re-estimation of markers effects (REE strategy). For each strategy, we first obtained a subset of 10000 markers, and from these markers, we obtained and subsequently evaluated five subsets (5000, 2500, 1000, 500, and 100 SNPs) for each strategy. For ORI strategy, the vector of markers effect was first estimated using all the markers and then ranked in decreasing order by their absolute values and selected 10000 SNPs with highest absolute value. Afterward, 5000 with highest absolute value was selected, this step was repeated until the subset of 100 SNPs was obtained. For REE strategy, the subset containing 10000 SNPs was obtained in the same way of ORI strategy. However, the markers effects were re-estimated, ranked in decreasing order by their absolute values and then, we selected 5000 SNPs with highest absolute value. This procedure was repeated until the subset of 100 SNP was obtained.

Main genotypic effect (MGE) model
The six marker subsets described above were used to perform genomic prediction, using the main genotypic effect model. This model fits the data separately from each subset and considers the main effect of genotypes. In matrix notation, the model is written as: where = ( 1 , … , ) ′ is the response vector, and represents the observations in the i th genotype ( = 1, … , ) in each markers subset; is the general mean; is the incidence matrix that connects the random genetic effects with phenotypes; is the random genetic effects, and the residual random effects. The MGE model (5) assumes that the distribution of the vector is multivariate normal with mean zero and a covariance matrix , with ~(0, 2 ), where 2 is the genetic variance component of and is a symmetric, positive semi-definite matrix, denoting the variance-covariance of the genetic values constructed from the genomic molecular markers by subsets. The error was assumed to be independent of each other and normally distributed ~( , 2 ). Therefore, is an approximation of the true unknown genetic values, and captures the residual genetic effects that were not explained by plus other non-genetic effects that approximate the errors.
The main genotypic effect model was used with two kernel regression methods for each markers subset.

Variance components, prediction accuracy, and bias parameters
The MGE model fitted with GB and GK methods were used on the whole RICE, HEL, and USP datasets for all the traits. Before fitting the models, the phenotypic data were centered and standardized (i.e., each phenotypic data point was centered by subtracting the overall mean and standardized by dividing by the sample standard deviation). These analyses were performed to derive estimates of variance components and genomic heritability. Since the data was standardized, the summation of the variance components of one specific model will approximate 1.
The prediction accuracy (PA) was performed separately for each trait and marker subset (all markers and the five markers subsets) for MGE-GB and MGE-GK models, and was assessed using thirty random partitions (replications), with 80% of the hybrids/lines comprising the training set (TRN), and the remaining 20% comprising the testing set (TST). For this validation procedure, all the parameters of the models (including variance components resulting from residual effects and genetic effects) were re-estimated from TRN data in each of the TRN-TST partitions.
For each TRN-TST partition, models were fitted to the TRN data set, and prediction accuracy was assessed by computing Pearson's correlation between predictions and phenotypes in the TST dataset. The same TRN-TST partitions were used to evaluate the prediction accuracy of each of the models. Thus, thirty correlations were computed for each model and trait.
We used two approaches to measuring the bias of GEBVs for ORI and REE strategies,

Relative efficiency of genomic prediction
The relative efficiency of genomic prediction for the all proposed scenarios was compared to the traditional phenotypic selection. For this, two parameters were used: i) the relative efficiency per cycle (REc), which was calculated by √ 2 ⁄ , were is the accuracy of model and is the phenotypic heritability of trait (Hoffstetter et al. 2016b); ii) the coincidence of selection, which was estimated based on the top 10% and 20% individuals identified by TRN-TST and the individuals with the best phenotypic performance for each trait.

Software
We used the rrBLUP-R package (Endelman, 2011) to run the RR-BLUP model. The BGLR package (de los Campos and Pérez-Rodríguez, 2016) was used to carry out the MGE models, using a total of 30,000 MCMC iterations, 5,000 for burn-in, and five for thinning.
Appendix B has the R codes used to obtained markers subsets by ORI and REE strategies, and how we fitted the MGE model with GBLUP and Gaussian kernels.

Phenotypic data
We observed for both traits, plant height (PH) and grain yield (GY), in the RICE, HEL, and USP dataset, genetic variability according to the likelihood ratio test based on the joint analysis (Appendix A, Table A1). Similarly, random effects of Genotype × Year (RICE) and interaction Genotype × Site (HEL and USP) were significant, except for PH in USP dataset. This significant effect suggests the differential performance of genotypes between sites or year for the tested traits. Entry-mean based heritability was 0.89 (PH) and 0.78 (GY) for RICE, 0.87 (PH) and 0.77 (GY) for HEL, and 0.83 (PH) and 0.59 (GY) for USP dataset (Additional Table S1), reflecting good accuracy of phenotypic evaluation.

Variance components of markers subset
Variance components (residual and genetic) and genomic heritability were obtained by MGE model through GBLUP and GK kernels (MGE-GB and MGE-GK).

RICE dataset
Using markers subsets obtained by REE strategy, for PH and GY trait, the estimated residual variance components were smaller than those obtained by ORI strategy (Figures 1A and   1D). The results taken through of the MGE-GK and MGE-GB models indicated that the use of small markers subset induced a larger reduction in the estimated residual variance. The variance component of genetic effects for each subset increased when MGE-GK model used rather than MGE-GB model ( Figures 1B and 1E). Furthermore, by the REE strategy, the genetic effects were slightly higher than ORI strategy, except for 500 markers dataset. The genomic heritability ranged from 0.44 (all markers MGE-GB model) to 0.90 (markers subset by REE strategy). When only 500 to 5000 markers were used to compute the GB or GK-matrix, the values ranged from 0.75 to 0.90. For PH, the results showed similar tendency with GY ( Figure 1E).

HEL dataset
The results of residual variance showed, for PH and GY, that MGE-GK tends to fit better the data than MGE-GB. It is because no consistent difference was found using all or subset of markers through the former method (Figures 2A and 2D). For GY, the genetic variance for MGE-GB (range from 0.20 to 0.30) was always smaller than MGE-GK (range from 1.20 to 0.60) ( Figure 2E). Genetic variance was higher, around 1.30 to MGE-GK and 0.80 to MGE-GB, using all markers for PH ( Figure 2B). However, using markers subsets by MGE-GK models, the genetic variances were always higher (around 0.70) than those obtained by MGE-GB (range from 0.50 to 0.30). Hence, genomic heritabilities by the former were greater than by the latter for all subsets of markers and both traits ( Figures 2C and 2F).

USP dataset
For GY, we observed higher values for residual variance using MGE-GB than by MGE-GK models ( Figure 3D). Moreover, the genetic variance was smaller, mainly through MGE-GB than by MGE-GK ( Figure 3E). Using all markers, we can see slightly higher values than using markers subset. The genomic heritability was lower using the former (range from 0.20 to 0.10) than the latter (range from 0.30 to 0.45) ( Figure 3F). The residual and genetic variance and heritability for PH followed the same trend of GY ( Figure 3A, 3B, and 3C).

RICE dataset
Results for GY showed that the prediction accuracy (PA) using markers subsets was better than using all markers ( Figure 4). The PA using all markers was 0.297 for MGE-GB and 0.291 for MGE-GK. However, when the markers subsets were used to predict, we observed a substantial increase in PA. For instance, for MGE-GB and MGE-GK models, the increase using 500 markers was 222% and 213%, respectively. This subset had the better PA and was obtained by REE strategy. Overall, comparing the two approaches, we verified REE was better than ORI.
The PA by MGE-GB ranged from 0.622 (100 markers) to 0.759 (1000 markers) via ORI strategy and ranged from 0.819 (100 markers) to 0.959 (500 markers) via REE strategy. In general, the MGE model using GB kernel is slightly better than GK kernel. For PH, the results are similar with GY ( Figure 4). The PA using all markers was low, where for MGE-GB model was 0.382 and for MGE-GK was 0.406. The increase in PA using markers subset was 156% and 132% using MGE-GB and MGE-GK models, respectively. The marker subset that showed the best PA (0.978 by MGE-GB and 0.945 by MGE-GK) was with 500 markers obtained by REE strategy.
The coefficient of regression (slope) of phenotype on GEBVs was calculated as a measurement of the bias of each strategy-model-subset (Appendix A, Table A2). Ideally, a value of β = 1 indicates no bias in the prediction. For GY and PH, comparing ORI or REE approaches and MGE-GB or MGE-GK models the slopes for MGE-GB model with ORI strategy were close to one, indicating no consistent bias in the prediction. In general, the fewer the markers, the smaller the bias. Furthermore, the value of β derived from MGE-GK model with REE strategy was higher for all traits (average across traits equal to 1.49). Moreover, the use of subsets with fewer markers resulted in slightly higher reliabilities for both traits (Appendix A, Table A2). Overall, for PH the REE strategy produced slightly better values of reliability than ORI approach.

HEL dataset
The results obtained for GY indicated that using markers subsets showed an increase in prediction accuracy ( Figure 5). The PA using all markers was intermediate for MGE-GB (0.548) and MGE-GK (0.674). We observed an increase in PA using markers subsets. For instance, using just 500 markers, the PA was 0.676 (MGE-GB), and 5000 markers were 0.721 (MGE-GK), which represents an improvement of 23% and 7% for MGE-GB and MGE-GK models, respectively. Considering the comparison between the two strategies, the REE was better when the prediction used MGE-GK, where the PA ranged from 0.676 (100 markers) to 0.721 (5000 markers). For PH, the PA was high (>0.81) for both models using all and subsets of markers ( Figure 5). The increase using markers subset was 8% and 4% for MGE-GB and MGE-GK models, respectively.
The coefficient of regression (slope) for all scenarios was close to one, indicating no consistent bias in the prediction. In general, the reliabilities of the models for GY showed that MGE-GK (0.794 to 0.860) showed higher reliabilities than MGE-GB model (0.664 to 0.809) (Appendix A, Table A2). On the other hand, for PH, all scenarios had similar reliabilities (> 0.809).

USP dataset
For both traits, all scenarios presented similar accuracies ( Figure 6). The slope for GY and PH traits, for both models and subsets were close to one, indicating no consistent bias in the prediction (Appendix A, Table A2).

RICE dataset
Considering 10% of an intensity of selection (IS), we observed a higher increase in the coincidence of selection (CS) when markers subsets were used (Table 1 and Figure 7A-7B). For instance, the CS values for GY by MGE-GB model ranged from 26% (all markers) to 81% (500 markers) and by MGE-GK range from 29% (all markers) to 70% (500 markers). The REE strategy selects markers that provide greater coincidence of genomic with phenotypic selection.
Under 20% as IS, the CS showed the same tendency of that in 10% IS. However, the values were slightly higher, when MGE-GB model was used.
Regarding the relative efficiency (REc) per cycle, MGE-GB using 2500, 1000, and 500 markers, and MGE-GK using 500 markers subset, obtained by REE strategy, have values equal or greater than 1, showing the efficiency of genomic prediction over the phenotypic selection.

HEL dataset
Following the trend described above, when markers subsets are used, we observed an increase in the CS via MGE-GB (Table 1 and Figure 7C-7D). However, this increase is similar with those obtained using MGE-GK. For instance, for PH, the values through MGE-GB ranged from 42% (all markers) to 48% (500 markers) and by MGE-GK from 44% (all markers) to 47% (100 markers). The values considering 20% as IS were slightly higher than by selecting 10%.
Regarding the relative efficiency (REc) per cycle for PH and GY, the values were less than 1, showing that, in these cases, the efficiency of genomic prediction did not exceed phenotypic selection (Table 1 and Figure 7C-7D).

USP dataset
Under 10% as IS, for GY, there was a slight increase (4%) in the CS, considering markers subset rather than all markers by MGE-GB model (Table 1 and Figure 7E-7F). MGE-GK prone to produce lower values of CS using markers subset. In general, the CS values were similar for both models, except when all markers are used. Using 20% as IS, there was a higher coincidence via MGE-GK. On the other hand, for PH, the results reveal that the coincidences of selection are similar among the models, methods and intensities of selection.
Considering the relative efficiency (REc) per cycle for PH (0.74 to 0.77) and GY (0.65 to 0.71), the values were less than 1, showing that, in these cases, the efficiency of genomic prediction did not exceed phenotypic selection (Table 1 and Figure 7E-7F).

DISCUSSION
Several studies have documented the benefits of using markers subset to improve the prediction accuracy in genomic selection analyses (Moser et al. 2010;Resende et al. 2012;Spindel et al. 2015;Tayeh et al. 2015;Hoffstetter et al. 2016a). In this research, two strategies were used to obtain markers subsets where the prediction accuracy, using each subset or all markers, was compared by a main genetic effect (MGE) model using GB (GBLUP) and GK (Gaussian) kernels. We found that, by selecting subsets of SNPs, accuracy is considerably increased.
Furthermore, the REE strategy is the best to obtain markers subsets, due to the re-estimation of the marker effect, increasing accuracy and reducing bias of GEBVs. According to Porto-Neto et al. (2015), SNP selection would take into account different marker allele frequencies between individuals, and the missing heritability can be largely recovered, leading to improvements of accuracies. In general, many sequence variants associated with complex traits have small effects and low repeatability (Bian and Holland 2017). The ORI and REE strategies allow eliminating many redundant markers with small effect, increasing accuracy. This procedure allows reducing the multicollinearity problem between markers, which occurs especially because markers in close positions are expected to be highly correlated. Similarly, Resende et al. (2012) also used an approach of re-estimating markers effect and found higher performance compared to Bayesian methods for rust disease resistance and wood density traits in Loblolly pine.
Several studies showed that selecting markers from a high density genotyping, can be an effective alternative to genomic prediction (Spindel et al. 2015(Spindel et al. , 2016Tayeh et al. 2015;Ma et al. 2016). Another way would be to weight the genomic relationship matrix with previously selected an increase in accuracy add markers, selected with highest significance test, largest absolute effect, or largest genetic variance contributed by the locus, in high-density SNP chips for animal breeding. However, we used only markers with higher effects to build the GRM and concluded that this approach is an efficient and easy way to be implemented in breeding programs.
Recently for a rice breeding program, Yu et al. (2014) developed accurate fixed arrays with 6 K SNPs from a high density genotyping, selecting representative SNPs obtained by genetic diversity studies and markers located inside genes of important traits, which are evenly distributed on the 12 chromosomes. Furthermore, Spindel et al. (2016) using genomic prediction models that incorporate de novo GWAS, with the same rice dataset used in our work, found that approximately 5000 SNPs were effective for prediction, suggesting the possible use of smaller fixed SNP arrays.
In our results, eliminating markers with small effect, from 5000 to 500 SNPs, is a good strategy and has the potential for build a trait-specific low-density SNP chip, reducing costs of genotyping in a breeding program.
Considering a subset with 500 markers obtained by ORI and REE strategies, we observed a wide distribution of them on all chromosomes of maize and rice (Appendix A, Figure   A1, A2, and A3). In the RICE dataset, the re-estimation of markers effects provided a wider distribution of markers across the genome, than HEL and USP datasets, allowing a more pronounced increase in PA. Similar methods based on selected markers evenly distributed across the genome was reported in the literature (Spindel et al. 2015;Tayeh et al. 2015;Moser et al. 2010). For instance, Spindel et al. (2015) observed that subsets either distributed evenly throughout the genome was the most important contributor to PA and had less accuracy variance than chosen at random. Moreover, Moser et al. (2010) showed that a chip containing 3000 evenly spaced markers could provide approximately 90% of the accuracy achieved with a high-density.
Based on the maize hybrids datasets, we conclude that using Gaussian kernel we could obtain higher prediction accuracies compared to the linear kernel ( Figure 5 and 6). The Gaussian kernel can capture non-additive effects while GBLUP only additive effects Cuevas et al. 2016b). Maize hybrids present high heterosis due to non-additive effects (dominance and epistasis), explaining the results achieved. However, in the RICE dataset, the GBLUP model showed slightly higher prediction accuracies than GK in most scenarios.
Furthermore, the smaller number of markers used in genomic prediction, the smaller the difference of PA observed between GBLUP and GK. Therefore, a small number of markers is not able to capture non-additive effects between markers, especially in the complex trait as grain yield. Moreover, the heterosis present in maize hybrids may have contributed to the best result of using GK.
Regarding the bias parameters, we found slight differences among models and strategies (Appendix A, Table A2). According to Porto-Neto et al. (2015), while accuracy affects the ranking of individuals, bias affects the range of estimated breeding values. The authors worked with animal data and found lower bias using SNPs selected from related individuals. We found similar results, where the differences between strategies and kernels were larger in GEBVs bias and prediction accuracy when applied to a population with more genetic diversity. In addition, Ma et al. (2015) analyzed same strategies to reduce the bias of predicted genetic trend and affirmed that this parameter should be taken into consideration in the validation of genomic predictions. Nonetheless, for RICE dataset the markers subsets obtained by ORI strategy evaluated by GBLUP kernel, had performed better than GK concerning reducing bias of GEBVs.
Comparing the coincidence of the genomic with the phenotypic selection and relative efficiency per cycle (REc), we verified more coincidence using markers subsets (Table 1).
Although the coincidence was not 100%, the genomic selection is still a promising alternative, mainly because it has the advantage of reducing the time of a breeding cycle .
Overall, our results showed that from with a high-density panel, it is possible to select the most informative markers to improve accuracy and build a low-cost array to implement genomic selection in breeding programs. The best strategy to obtain markers subsets is the re-estimation of the marker effect, which increases the accuracy and reduces the bias.

CONCLUSION
From with a high-density panel, it is possible to select the most informative markers to improve accuracy and build a low-cost array to implement genomic selection in breeding programs. The best strategy to obtain markers subsets is the re-estimation of the marker effect, which increases the accuracy and reduces the bias.              Appendix B Box B1. The simplified scripts that can be used to obtain the markers by ORI and REE strategy through rrBLUP R package.

GENOTYPE × ENVIRONMENT INTERACTION
ABSTRACT Multi-environment trials are routinely conducted in plant breeding to select candidates for the next selection cycle. In this study, we compare the prediction accuracy of four developed genomic-enabled prediction models: (1) singleenvironment, main genotypic effect model (SM), (2) multi-environment, main genotypic effects model (MM), (3) multi-environment, single variance G×E deviation model (MDs), and (4) multi-environment, environment-specific variance G×E deviation model (MDe). Each of these four models were fitted using two kernel methods: a linear kernel Genomic Best Linear Unbiased Predictor, GBLUP (GB), and a non-linear kernel Gaussian kernel (GK). The eight model-method combinations were applied to two extensive Brazilian maize data sets (HEL and USP data sets), having different numbers of maize hybrids evaluated in different environments for grain yield (GY), plant height (PH) and ear height (EH). Results show that the MDe and the MDs models fitted with the Gaussian kernel (MDe-GK, and MDs-GK) had the highest prediction accuracy. For GY in the HEL data set, the increase in prediction accuracy of SM-GK over SM-GB ranged from 9% to 32%. For the MM, MDs, and MDe models, the increase in prediction accuracy of GK over GB ranged from 9% to 49%. For GY in the USP data set, the increase in prediction accuracy of SM-GK over SM-GB ranged from 0% to 7%. For the MM, MDs, and MDe models, the increase in prediction accuracy of GK over GB ranged from 34% to 70%. For traits PH and EH, gains in prediction accuracy of models with GK compared to models with GB were smaller than those achieved in GY. Also these gains in prediction accuracy decreased when a more difficult prediction problem was studied.
Keywords: Genomic-enabled prediction; Maize multi-environment trials; Genotype × environment interaction ( developed and applied to simulated and real plant breeding data de los Campos et al., 2009de los Campos et al., , 2013Crossa et al., 2010Crossa et al., , 2011Pérez-Rodríguez et al., 2013;. In general, these studies showed good prediction accuracies for complex traits evaluated by means of random cross-validation partitions of the data. Genomic-enabled prediction models were originally developed for a single trait in a single environment. However, multi-environment plant breeding trials are routinely conducted to estimate and take advantage of genotype × environment interaction (G×E). Therefore, to implement GS strategies in plant breeding, G×E needs to be estimated, modeled, and predicted.
When genetic marker information is used for computing associations between individuals through the Genomic Relationships Matrix (GRM) G (Habier et al., 2007), this model is also

Phenotypic experimental data
This study considered two hybrid maize data sets. The first data set is from the Helix There was an imbalance in both data sets because not all hybrids were evaluated in all locations (incomplete field trials). The main traits of the two data sets were GY (grain yield in ton ha -1 ), PH (plant height in cm) and EH (ear height in cm). Plant height was measured from the ground to the flag leaf and EH from the ground to the base of the ear. The empirical distribution of the three evaluated traits (GY, PH, and EH) were symmetric in most of the environments ( Fig. 1). For the HEL data set, PH and EH were evaluated in three environments. Broad sense heritability (repeatability) were computed using the standard formula based on plot means ℎ 2 /[ ℎ 2 + ℎ 2 / + 2 / where ℎ 2 is the variance of the hybrids, ℎ 2 is the variance of the hybrids × location interaction and 2 is the residual error variance and s and r are the number of environments and replicates in each environment, respectively.

Genotypic data
The USP and HEL parent lines were genotyped with an Affymetrix  Axiom  Maize Genotyping Array of 616 K SNPs (Unterseer et al., 2014). Standard quality controls (QC) were applied to the data, removing markers with a Call Rate ≥ 0.95. The remaining missing data in lines were imputed with the Synbreed package (Wimmer et al., 2015) using the algorithms from the software Beagle 4.0 . The hybrid genotypes were obtained by genomic information of the parent inbred lines. Markers with Minor Allele Frequency (MAF) of ≤ 0.05 were removed. After applying QC, 52,811 and 54,113 SNPs were available to make the predictions in the HEL and USP data sets, respectively.

Data availability
The phenotypic and genotypic data for the maize hybrids included in this study can be found at http://hdl.handle.net/11529/10887. Each HEL and USP data set contains the data corresponding to each GB and GK kernel, as well as phenotypic data for each trait (PH, EH, and GY) in each environment.

Statistical models
Statistical models for genomic predictions taking into account G×E were proposed by In this study, four statistical prediction models were fitted to both data sets to study their prediction accuracy using random cross-validation schemes (Table 1). We also compared the prediction accuracy of two kernel regression methods in the four models. The models were: a single-environment, main genotypic effect model (

The single-environment, main genotypic effect model (SM)
The SM model fits the data from each environment separately and takes into account the main effect of genotypes. In matrix notation the model is written as: where = ( 1 , … , ) ′ is the response vector, and represents the observations in the i th line ( = 1, … , ) in each environment; is the general mean; is the incidence matrix that connects the random genetic effects with phenotypes; is the random genetic effects for each environment, and the residual random effects for each environment. The SM model (1) assumes that the distribution of the vector is multivariate normal with mean zero and a covariance matrix 2 , that is, ~( , 2 ), where 2 is the genetic variance component of in the j th environment and is a symmetric, positive semi-definite matrix, denoting the variance-covariance of the genetic values constructed from the genomic molecular markers. It is also assumed that the errors in each environment are independent with homogeneous variance, 2 , and distributed as ~( , 2 ) (where I is the identity matrix, and 2 is the residual variance in the j th environment). Therefore, is an approximation of the true unknown genetic values, and captures the residual genetic effects that were not explained by plus other nongenetic effects that approximate the errors. Note that matrix K may have different dimensions depending on the number of lines evaluated in each environment.
When markers have a more complex function, the Gaussian kernel (GK) has proved to be more effective for capturing complex cryptic interactions between markers, thereby improving the prediction accuracy of the model . The entries of the Gaussian kernel are computed as (`) = exp (−ℎ2 ), where ` is the Euclidean distance between the individuals th and ′th ( = 1, … , ) given by the markers; ℎ > 0 is the bandwidth parameter that controls the rate of decay of values (Pérez-Rodríguez et al., 2012;. In this study, GK takes the form (`) = exp (−ℎ2 /median(2 )), where ℎ = 1 and the median of the distances is used as a scaling factor (Crossa et al., 2010). de los Campos et al.

The multi-environment, main genotypic effect model (MM)
One multi-environment model considers the main fixed effects of environments, as well as the random genetic effects across environments  (Table 1).

Multi-environment, single variance genotype × environment deviation model (MDs)
This model extends (2)  Therefore, under these conditions where represents the kernel constructed from the molecular markers of the lines in the j th environment. Just as in model (2), matrix is used in the variance-covariance for of model (3) and is also a component of the variance-covariance of . Then, kernel matrix can be constructed using the linear kernel (MDs-GB) or the Gaussian Kernel (GK) (MDs-GK) ( Table   1).

Multi-environment, environment-specific variance genotype × environment deviation model (MDe)
The model of López-Cruz et al. (2015) separates the genetic effects of markers into the main marker effects across all environments and the specific marker effects in each environment.
Thus, this model in matrix notation is For the MDe model (4) of the vector, represents the main effect of markers (across all environments) with a variance-covariance structure similar to those used in models (2) and (3) (3), that is, ~( , ), where is: Note that matrix can be expressed as the sum of matrices and the effects given by are specific for the th environment; thus, it has a normal distribution with mean and variance equal to 0, except for the th environment, which has a variance-covariance matrix 2 . The importance of these two terms ( and ) of MDe model (4) (Table 1).

Estimating variance components using full data analyses
The four models -SM, MM, MDs, and MDe fitted with GBLUP and GK methodswere used on the entire HEL and USP data sets for all the traits. To fit the models, the phenotypic data were centered and standardized (i.e., each phenotypic data point was centered by subtracting the overall mean of all environments and standardized by dividing by the sample standard deviation across all environments). These analyses were performed to derive estimates of variance components. The posterior variance components resulting from the residual effects, genetic main effect, and genetic environment-specific effects of the four models described above for three traits and for environments in HEL and USP data sets were computed and reported.
Since the data were standardized, the summation of the variance components approximates 1.

Assessing prediction accuracy by random cross-validation
The prediction accuracy of the single-environment model-methods was assessed using 50 random partitions organized in 5 folds with 10 random partitions each, with 80% of the hybrids comprising the training (TRN) set and the remaining 20% of the individuals comprising the testing (TST) set. This was performed separately in each environment as the singleenvironment models were fitted separately for each environment. For this validation procedure, all the parameters of the models (including variance components resulting from residual effects, and genetic effects) were re-estimated from TRN data in each of the TRN-TST partitions.
For the multi-environment models, the prediction accuracy of the model-methods for predicting the pattern of missing values in each environment was generated using two different cross-validation (CV) designs (Burgueño et al., 2012) Therefore, adjusting models for full data, making inference on the parameters, and assessing predictions in TRN-TST partitions were based on 30,000 samples collected from the posterior and predictive distributions after discarding 5000 samples as burn-in.

Software
The aforementioned models can be implemented using the R (R Core Team, 2016) package Bayesian generalized linear regression (BGLR) (de los Campos and Pérez-Rodríguez, 2016; Pérez-Rodríguez and de los Campos, 2014 ). Appendix C gives the R codes used in the models included in this study when the numbers of individuals in different environments vary. It should be pointed out that in previous studies using the  and López-Cruz et al. (2015) models, the number of entries in each environment was the same; thus, the model fitting and the random cross-validation partition scheme were simpler than in this case, where there are different numbers of entries in different environments.

Descriptive Statistics
Box-plots of GY, PH and EH in each environment for the two data sets are depicted in had better potential for all traits (Fig. 1A).
All pair-wise environments for all traits in the HEL and USP data sets had high and positive correlations (Tables 2 and 3). As expected, for PH and EH, the correlations were higher than for the complex trait GY. In all environments, GY, PH and EH showed high heritability in the HEL data set: GY ranged from 0.60 to 0.

Estimating variance components
The distribution of the residuals after fitting the eight model-method combinations to the two data sets was approximately normal (data not shown). Phenotypic (and marker) data were standardized; thus, summation of variance components approximated 1 for GB models.
Deviations may be due to the high degree of imbalance in the number of hybrids in different environments, especially for data set HEL.

Single-environment
For the SM model for the three traits in each environment, the estimated residual variance components for the GK method were smaller than those for the GB method (Table 4).
In contrast, the variance component of genetic effects for each environment increased for the SM-GK model-method as compared to the genetic effects of model-method SM-GB.

Multi-environment
For the MM model, the residual variance components for MM-GK were smaller than the residuals for MM-GB for all traits, whereas for the genetic effects, the opposite occurred: the variance component of MM-GK was higher than that of MM-GB. For the MDs models, the residual variance components for MDs-GK were always smaller than the residual variances for MDs-GB, whereas the opposite occurred for the genetic main effect and genetic interaction effect ( Table 4).
The results of the MM and MDs models indicated that the inclusion of the interaction term (G×E) induced a larger reduction in the estimated residual variance for GY; for traits PH and EH, these reductions in residual variance were smaller than those found for GY from the MM and MDs models. For the MDe model, the residual variance components of MDe-GK were smaller than those of the MDe-GB for all traits, whereas the variance components for the genetic main effect and genetic environment specific effect were higher for the GK models than for the GB model. The genetic main effect values are related to the correlations between environments, which are positive and range from medium to high. For the variance component associated with the genetic interaction effect, the values are lower because they decreased G×E in environments with a positive correlation. In general, the fit of MDs is slightly worse than the fit of MDe models because the genetic interaction effect component has a small effect on all traits, but mainly on EH and PH (Table 4).

Single-environment
For all traits, the variance components of genetic effects were higher when using the GK than the GB method in all environments and traits for the single-environment, main genotypic effect model (SM) ( Table 5). The residual variance for SM-GK was always smaller than the residual for SM-GB for all traits.

Multi-environment
The residual variance components of the MM and MDs models for all traits were very similar, as were most of the variance components for the genetic main effects. For all the environments and traits, the residuals from the GK were smaller than the residuals from the GB for the MM, MDs, and MDe models. The variance components of the genetic main effects of MDs-GB and MDe-GB were consistently smaller than the genetic main effect of model-methods MDs-GK and MDe-GK. Also, the variance components of the genetic environment specific effects of the GK models were higher than those of the GB models for all traits (Table 5).

Prediction accuracy of the models with GBLUP and GK methods
Results of the prediction accuracy of the eight model-method combinations for the random cross-validation CV2 of the two data sets are given in Figs. 2-5 and Tables 6 and 7.
Results of the prediction accuracy of the eight model-method combinations for random crossvalidation CV1 of the two data sets are given in Tables D1 and D2 in Appendix D. Note that for the single-environment models, prediction accuracy was assessed by 50 random partitions of the data into 80% TRN and 20% TST.

Single-environment
Results for CV2 for GY showed that the prediction accuracy of the SM-GK modelmethod was higher than for SM-GB in all environments (Fig. 2). For GY in environment IP, prediction accuracies were 0.64 for SM-GB and 0.73 for SM-GK, while in environment PM, prediction accuracies were 0.68 for SM-GB and 0.74 for SM-GK (Table 6). The percent change in accuracy for single-environment GB vs. single-environment GK for GY ranged from 9 to 32% for the PM and NM environments, respectively. For trait PH, prediction accuracies were high: >0.65 for SM-GB and >0.69 for SM-GK. For trait EH, prediction accuracies were >0.64 for SM-GB and >0.67 for SM-GK; the increase in prediction accuracy ranged from 3 to 6% for traits EH and PH (Table 6).

Multi-environment
For GY, the best model-methods for CV2 were MDs-GK and MDe-GK for the five environments, followed by MDs-GB in the IP environment, MM-GK in the NM, PM and SO environments and by MDs-GB and MDe-GB in the SE environment (Fig. 3). For trait GY, the percent change in accuracy for MM-GB vs MM-GK was high, ranging from 9% in PM to 48% in SE environment; for MDs-GB vs MDs-GK, the % change ranged from 13% in PM to 47% in the SO environment; and for MDe-GB vs MDe-GK, the % change ranged from 13% in PM to 49% in SE environment (Table 6). These results show that for trait GY, the increase in prediction accuracy was high using the GK method. For less complex traits PH and EH, the percent change of GK over GB for all models showed a small increase in prediction accuracy, with the percent change ranging from 1% to 6%; there was no difference for trait EH in the IP environment.
For GY, MDs-GK and MDe-GK were, for CV2, the models that had the best prediction accuracy in all environments, with values ranging from 0.58 in SE to 0.81 in the IP environment. These values were higher for a complex trait like GY (Fig. 3). For the PH trait, MDs-GK and MDe-GK were the best models, with a prediction accuracy of 0.81-0.82 for all environments; for trait EH, these values ranged from 0.77 in IP to 0.81 in PM. In general, MDs models had similar results across environments.
Prediction accuracy for random cross-validation CV1 decreased (Table D1, Appendix D) as compared with those computed for CV2 for all traits and models. For example, for environment IP, MM-GB in CV2 gave a correlation of 0.56, whereas for CV1, prediction accuracy was 0.48, that is, an 8% decrease in accuracy, whereas for environment NM, the decrease is about 10%. Similar decreases in accuracy from CV1 to CV2 were found for other traits and model-method combinations. However, for GY, models with GK showed a greater increase in prediction accuracy over models with GB under CV1 (Table D1, Appendix D) than the increases achieved by those models under CV2 (Table 6). The % change in prediction accuracy of CV1 vs CV2 for models with GK vs GB did not change for traits PH and EH.

Summary of results
Prediction accuracy of trait GY in the 5 testing environments increased from 9% to 48% for models with method GK compared to models with method GB. For GY, increases in prediction accuracy of environments from multi-environment models MM, MDs, and MDe with GK vs MM, MDs, and MDe with GB tended to be higher than those achieved by the singleenvironment model (SM-GK vs SM-GB). For traits PH and EH in 3 environments, the increase in prediction accuracy of method GK over method GB ranged only from 0% to 6% and increases in prediction accuracy from multi-environment models MM, MDs, and MDe with GK vs MM, MDs, and MDe with GB tended to be lower than those achieved by the singleenvironment model (SM-GK vs SM-GB).
Although prediction accuracy for all model-method combinations decreased in CV1 vs CV2, the % change in the prediction accuracy of models with GK vs models GB did not change and even increased for GY.

USP data set
Single-environment Table 7 shows the correlations between observed and predicted values for all models in the USP data set. In general, for traits GY, PH and EH, all prediction accuracies of the singleenvironment models were similar using the GB or GK methods (Fig. 4). In these cases, the percentage change in accuracy for SM-GB vs SM-GK was 7% (GY) in environment P-LN and 5% (PH) in environment A-LN. For GY, the prediction accuracy in the A-IN environment (0.40 for SM-GB and SM-GK) was greater than in other environments. Gains in prediction accuracy of SM-GK over SM-GB for trait GY in this data set are much lower than those achieved in the HEL data set.

Multi-environment
Between the GB and GK methods, there were large increases in prediction accuracy using the GK method, especially for traits GY and PH (Fig. 5). For trait GY, the MM model gave a high percentage change in accuracy for GB vs GK ranging from 34 to 59%. In the MDs model, the percentage change ranged from 35 to 65%, and for the MDe model, the percentage change ranged from 37 to 70% (Table 7).
For PH, the increase using the GK method was similar to those found for the GY trait. The prediction accuracy of random cross-validation CV1 decreased (Table D2, Appendix D) compared to the accuracy obtained with CV2 for all traits and models. In general, compared with other models, MDs and MDe decreased the accuracy of CV1 vs CV2 more; this was more pronounced for GY in some environments than in other environments. Also, the decrease in the accuracy of CV1 vs CV2 was higher for GY than for PH and EH.

Summary of results
The prediction accuracy of trait GY in 4 environments increased from 34% to 70% for multi-environment models MM, MDs, and MDe with GK vs MM, MDs, and MDe with GB.
These gains in accuracy were much higher than those achieved by single-environment SM-GK over SM-GB (from 0% to 7%). Similar patterns were found for traits PH and EH, where the increase in prediction accuracy of method GK over method GB ranged from -2% to 5% for the single-environment model (SM) and drastically increased for multi-environment models MM, MDs, and MDe with GK vs MM, MDs, and MDe with GB, ranging from 15% to 57%.
Prediction accuracy for all model-method combinations decreased in CV1 vs CV2 and the % change in accuracy of models with GK vs models with GB decreased for all three traits.
The increase in accuracy of models with GK vs models with GB did not occur under the CV1 design.

DISCUSSION
Several studies have documented the benefits of using a non-linear Gaussian kernel in multi-environment models for capturing small complex interactions among markers and increasing prediction accuracy (Gianola et al., 2006Crossa et al., 2010;Pérez-Rodríguez et al., 2013;Cuevas et al., , 2016b. In this research, genomeenabled prediction accuracy was studied in two different hybrid maize data sets using multienvironment models with the GBLUP and Gaussian kernel methods. In this context, we proposed including the Gaussian kernel in the model of Jarquin et al. (2014), which accounts for environment information without and with interaction terms. Prediction accuracy was obtained using traits with different genetic architecture and heritability values that ranged from low to high.

Prediction accuracy differences in datasets, methods, cross-validation designs and G×E
We compared the single-environment model with two G×E models. The first G×E model can accommodate environmental covariables ) that were not available in this study; the second G×E model (López-Cruz et al., 2015) accommodates environment main effects and environment-specific variances; both models used the GB and Gaussian kernel methods proposed by . These two G×E models have demonstrated good prediction accuracy when used in genomic-enabled prediction studies (López-Cruz et al., 2015;Zhang et al., 2015;Saint Pierre et al., 2016). We found that GK methods improved the prediction ability of all single-environment and multi-environment models for CV2 design for all traits in HEL and USP datasets. Prediction of all traits with multi-environment models incorporating the G×E term and the GK method gave better prediction accuracy (especially for complex traits such as GY) than the other model-methods. The increases in prediction accuracy using models with GK under random cross-validation 2 (CV2) are also reflected under random cross-validation 1 (CV1) for HEL data set, but to a lesser extent due to the difficulty of borrowing information of unobserved (unphenotyped) lines in all environments.
Similar decreases in prediction accuracy were found by López-Cruz et al. (2015) when attempting to predict wheat lines in untested (unobserved) environments under the CV1 random partition scheme.
Furthermore, for trait GY, differences in CV1 vs CV2 for GB and GK methods for the two data sets might also be due to other factors such: (1) differences in dataset repeatability, higher in HEL data sets (0.81, 0.60, 0.86, 0.69, and 0.78 for each of the 5 environments) than in USP (0.22,0.42,0.19,and 0.42 for each of the four environments) (Tables 2 and 3, respectively) due to better quality of the replicated trials in HEL than the unreplicated trials in USP; (2) more opportunity to borrow information from the close relatives hybrids in HEL (genetic diversity =0.175) than from the less related hybrids included in USP dataset (genetic diversity =0.372). For h=1 (bandwidth parameter) the Gaussian kernel is a direct function of the squared Euclidean distance between hybrids based on markers (`), thus for a set of unrelated hybrids `→ large and →I (identity matrix), that is, the Gaussian kernel weight more heavily the within (own) hybrid information compared to the between hybrid (from relatives) information than GB (the marker information is of no use) (de los Campos et al., 2010). On the other hand, for a set of related hybrids then `→ 0 and →1 (matrix of ones), that is, marker information given by the Gaussian kernel weight equally or less the own hybrid information compared to information from relatives given by GB. Therefore, we speculate that for related hybrids of HEL dataset GK gives heavier weights to the within hybrid as well as between hybrids than the GB method does, that is why for GY, CV1 and CV2 for the multi-environment models gave very good increase in prediction accuracy over the GB, (3) for CV2 G×E modulates both values (within and between hybrid information); hybrids in USP (doubled the genetic diversity) less related than those from HEL data set, then it is expected negligible values between hybrids, and similar values assign by GK and GB to within hybrid information; for CV2, predictions accuracy depends, to a great extent, on the correlations between environments.
In summary, for CV2, for the prediction accuracy of GY in the HEL dataset, the method GK weights the own and the between hybrids performance as well as the relationship between environment (G×E modeled by MDs y MDe), whereas for USP, GK weights only the own hybrid performance and the relationship between environment becomes more important that in the HEL data set.
In the studies of López-Cruz et al. (2015) and , the data were balanced in the sense that all the individuals were included at the same time in all the environments. On the contrary, in this study, we had a great deal of imbalance because different numbers of maize hybrids were included in different environments; this was more pronounced in the HEL data set than in the USP data set and different R scripts were necessary for implementing the model-method combinations.

Prediction accuracy using linear and non-linear kernel methods
According to Gianola et al. (2014), Gaussian kernel has better predictive ability and a more flexible structure than GBLUP. Another point is that GK can capture non-additive effects between markers. Jiang and Reif (2015) evaluated maize data sets and investigated whether the prediction accuracy across connected bi-parental families can be increased by modeling additive × additive epistasis; the authors found that the prediction accuracy of RKHS (including epistasis) was superior to that of GBLUP (ignoring epistasis). Our study clearly shows a large increase in predictive ability when the G×E model and the GK method are combined for all traits, but mainly for GY (both data sets) and PH (USP data set). This indicates that to increase predictive ability, it is important to consider the non-additive effect (i.e., epistasis) as the genetic relatedness across connected populations.
There are different choices for computing kernel functions: for example, linear kernel matrices incorporate only additive effects of the markers, polynomials kernels of different orders might incorporate different degrees of marker interactions, and the Gaussian kernel function uses complex epistatic marker interaction (Akdemir and Jannink, 2015). In genomic selection, additive kernels of the GBLUP type are employed when predicting breeding values, whereas when attempting to predict genetic values, Gaussian kernels would be more appropriate. Akdemir and Jannink (2015) demonstrated that epistatic marker effects in local regions of the chromosome with low recombination are stable through generations and offer the opportunity to exploit epistasis for improving genomic-enabled prediction accuracy; the authors defined local kernels in regions of the genome and calculated separate kernels for each region. The results shown in this study clearly showed the benefit of exploiting these local epistatic effects captured in the Gaussian Kernel and their interaction with environments.
The flexible structure of the MDs and MDe models is important, especially when combined with kernels that capture non-additive effects, as previously proven. In general, in the two maize data sets, the MDs-GK and MDe-GK models had better prediction accuracy than other models. The increase was not higher due to the occurrence of intermediate-to-high positive correlations between the analyzed environments, resulting in low G×E interaction. These results corroborate those obtained by López-Cruz et al. (2015) and , who observed that the intensity of the environmental correlation is related to the proportion of the genomic variance explained by the genetic main effects of markers across environments and geneticspecific effects of markers in environments.

Better fit of the G×E Gaussian kernel models
For GY of the HEL data set, the better fit of the model-method MDe-GB over the

Prediction accuracy using multi-environment models
In this study, all pairwise correlations between environments were high and positive.
This is important, because the G×E model of López-Cruz et al. (2015) has the limitation of better and more efficient prediction when applied to subsets of environments that have positive and similar correlations . (2016b) developed multi-environment Bayesian genomic models that allow an arbitrary genetic covariance structure between environments, because an unstructured covariance matrix was used and its parameters were estimated from the data.
In this study, the prediction of multi-environment models was assessed by applying the cross-validation strategy called CV2 in Burgueño et al. (2012) and , where some lines are represented in some environments but not in others. This CV2 validation strategy performed better than cross-validation strategy CV1 (where lines have not been evaluated in any field trials) proposed by Burgueño et al. (2012), when applied in multi-environment models López-Cruz et al., 2015;Crossa et al., 2016b;Saint Pierre et al., 2016).
When introducing interaction effects, models MDs and the MDe showed increases in prediction accuracy for GY in most environments for the HEL and USP data sets over the single-environment model. Furthermore, models with GK had superior prediction accuracy than GB models. For traits PH and EH in both data sets, not much increase in prediction accuracy of MDs and MDe models was achieved over the single-environment model. These results were expected because the genetic architecture of traits PH and EH is less complex than that of GY and less influenced by environmental factors. When main marker effects and interaction effects are introduced in the model using covariance structures this improves prediction accuracy.
The increase in accuracy when the interaction term for complex traits was included Mexico showed a consistent increase in the genomic prediction accuracy of the MDs-GB model over the MM-GB and SM-GB models (Sukumaran et al., 2017). One limitation of the data used in this study is that only one year was available for assessing the genomic-enabled accuracy of the various models.
Our results show that predictions with medium-to-high accuracy for genomic selection programs can be expected in environments with low-to-high heritability. Despite the low heritability of GY in P-LN and A-LN environments, the results were similar to results in other environments that have medium-to-high heritability. But in general, for the USP data set, prediction accuracy of all target traits under stress conditions was lower than under ideal nitrogen conditions, while in stress environments, heritability was lower than under ideal nitrogen conditions. These results agree with those obtained by Zhang et al. (2015), who found lower genomic prediction accuracy under water stress conditions than under ideal conditions, especially for complex traits such as GY.

CONCLUSION
Incorporating the Gaussian kernel method into the model of  increased prediction accuracy as compared to the model used with the linear kernel GBLUP; these results were found in both data sets with the models used in this study: MM across environments, MDs and MDe. Other results of the application of eight model-method combinations between four models (SM, MM, MDs, and MDE) and two kernel methods (GBLUP and Gaussian kernel) in two extensive maize data sets show that (1) genomic models incorporating G×E interaction had higher prediction accuracy than single-environment models; and (2) models with non-linear Gaussian kernel had higher prediction accuracy than models with linear kernel GBLUP. The model-method combinations with the highest prediction accuracy were MDe-GK and MDs-GK. Results of this study indicated that by employing appropriate statistical genomic-enabled prediction models the researchers and plant breeders can improve the prediction of hybrids that were not evaluated in several environments. Random cross-validation 2 (CV2) mimics sparse (incomplete) testing that allow to safe resources to the breeding program while improving prediction accuracy. Further research is required to compare the genomicenabled G×E kernel prediction models used in this study with the models recently developed by Cuevas et al. (2016b), who studied genomic-enabled G×E models by means of Kronecker products applied to unstructured covariance matrices. It is also necessary to develop efficient computing software for fitting the structures considered in the models of this study while maintaining the inferential advantages of the Markov Chain Monte Carlo.      Table 6. HEL data set. Mean correlation (for 50 random partitions CV2) for models singleenvironment, main genotypic effects model with GBLUP kernel method (SM-GB) and singleenvironment, main genotypic effects model with Gaussian kernel method (SM-GK), multienvironment, main genotypic effect model GBLUP kernel (MM-GB), multi-environment, main genotypic effect Gaussian kernel (MM-GK), multi-environment, single variance G×E deviation model GBLUP kernel (MDs-GB), multi-environment, single variance G×E deviation model Gaussian kernel (MDs-GK), multi-environment, environment-specific variance G×E deviation model GBLUP kernel (MDe-GB), multi-environment, environment-specific variance G×E deviation model Gaussian kernel (MDe-GK) for three traits, grain yield, plant height, and ear height (standard deviation in parentheses) Model

GENERAL CONCLUSION
From with a high-density panel, it is possible to select the most informative markers to improve accuracy and build a low-cost array to implement genomic selection in breeding programs.
The best strategy to obtain markers subsets is the re-estimation of the marker effect, which increases the accuracy and reduces the bias.
The use of Gaussian kernel and the including G×E effect, there is an increase in the accuracy of the genomic prediction models.