Genomic Prediction with Genotype by Environment Interaction Analysis for Kernel Zinc Concentration in Tropical Maize Germplasm

Zinc (Zn) deficiency is a major risk factor for human health, affecting about 30% of the world’s population. To study the potential of genomic selection (GS) for maize with increased Zn concentration, an association panel and two doubled haploid (DH) populations were evaluated in three environments. Three genomic prediction models, M (M1: Environment + Line, M2: Environment + Line + Genomic, and M3: Environment + Line + Genomic + Genomic x Environment) incorporating main effects (lines and genomic) and the interaction between genomic and environment (G x E) were assessed to estimate the prediction ability (rMP) for each model. Two distinct cross-validation (CV) schemes simulating two genomic prediction breeding scenarios were used. CV1 predicts the performance of newly developed lines, whereas CV2 predicts the performance of lines tested in sparse multi-location trials. Predictions for Zn in CV1 ranged from -0.01 to 0.56 for DH1, 0.04 to 0.50 for DH2 and -0.001 to 0.47 for the association panel. For CV2, rMP values ranged from 0.67 to 0.71 for DH1, 0.40 to 0.56 for DH2 and 0.64 to 0.72 for the association panel. The genomic prediction model which included G x E had the highest average rMP for both CV1 (0.39 and 0.44) and CV2 (0.71 and 0.51) for the association panel and DH2 population, respectively. These results suggest that GS has potential to accelerate breeding for enhanced kernel Zn concentration by facilitating selection of superior genotypes.

time-consuming. However, with the recent advances in genomics, new methods for plant breeding such as genomic selection (GS) can be used to identify genotypes with enhanced Zn concentration more efficiently and rapidly.
Kernel Zn concentration is determined at the end of a plant's life cycle, so GS can enable selection of promising genotypes earlier in the life cycle. This reduces the time and cost of phenotypic evaluation and may increase the genetic gain per unit time and cost (Heslot et al. 2015;Manickavelu et al. 2017;Arojju et al. 2019). The utility and effectiveness of GS has been examined for many different crop species, marker densities, traits and statistical models and varying levels of prediction accuracy have been achieved (de los Campos et al. 2009Crossa et al. 2010Crossa et al. , 2013Crossa et al. , 2014Jarquín et al. 2014;Pérez-Rodríguez et al. 2015;Zhang et al. 2015;Velu et al. 2016). Although the number of markers needed for accurate prediction of genotypic values depends on the extent of linkage disequilibrium between markers and QTL (Meuwissen et al. 2001), a higher marker density can improve the proportion of genetic variation explained by markers and thus result in higher prediction accuracy (Albrecht et al. 2011;Zhao et al. 2012;Combs and Bernardo 2013;Liu et al. 2018). Importantly, higher prediction accuracies have been obtained when genotypes of a population are closely related than when genetically unrelated (Pszczola et al. 2012;Combs and Bernardo 2013;Spindel and McCouch 2016).
Initially, GS models and methods were developed for single-environment analyses and they did not consider correlated environmental structures due to genotype by environment (G x E) interactions . The differential response of genotypes in different environments is a major challenge for breeders and can affect heritability and genotype ranking over environments (Monteverde et al. 2018). Multi-environment analysis can model G x E using genetic and residual covariance functions (Burgueño et al. 2012), markers and environmental covariates (Jarquín et al. 2014), or marker by environment (M x E) interactions (Lopez-Cruz et al. 2015). This approach to GS can successfully be used for biofortification breeding of maize because multi-environment testing is routinely used in the development and release of varieties.
Modeling covariance matrices to account for G x E allows the use of information from correlated environments (Burgueño et al. 2012). Mixed models that allow the incorporation of a genetic covariance matrix calculated from marker data, rather than assuming independence among genotypes improves the estimation of genetic effects (VanRaden 2008). The benefit of using genetic covariance matrices in G x E mixed models is that the model relates genotypes across locations even when the lines are not present in all locations (Monteverde et al. 2018). GS models capable of accounting for multi-environment data have extensively been studied in different crops (Zhang et al. 2015;Cuevas et al. 2016Cuevas et al. , 2017Velu et al. 2016;Jarquín et al. 2017;Sukumaran et al. 2017a;Monteverde et al. 2018;Roorkiwal et al. 2018). In those studies, incorporating G x E demonstrated a substantial increase in prediction accuracy relative to single-environment analyses.
Kernel Zn has been investigated in several quantitative trait loci (QTL) analyses in maize and each study has reported that Zn concentration is under the control of several loci. The phenotypic variation explained by those loci ranges from 5.9 to 48.8% (Zhou et al. 2010;Qin et al. 2012;Simić et al. 2012;Baxter et al. 2013;Jin et al. 2013;Zhang et al. 2017a;Hindu et al. 2018). A Meta-QTL analysis across several of those studies identified regions on chromosome 2 that might be important for kernel Zn concentration (Jin et al. 2013). Additionally, genomic regions associated with Zn concentration were recently reported in a genome-wide association study of maize inbreds adapted to the tropics (Hindu et al. 2018). Whereas some of the regions were novel, four of the twenty identified were located in previously reported QTL intervals.
A wide array of maize genetic studies has reported considerable effects of G x E interactions for kernel Zn concentration (Oikeh et al. 2003(Oikeh et al. , 2004Long et al. 2004;Chakraborti et al. 2009;Prasanna et al. 2011;Agrawal et al. 2012;Guleria et al. 2013). However, genotypes with high-Zn concentration have been identified in both tropical and temperate germplasm (Ahmadi et al. 1993;Bänziger and Long 2000;Brkic et al. 2004;Menkir 2008;Chakraborti et al. 2011;Prasanna et al. 2011;Hindu et al. 2018). Additionally, evaluation procedures for kernel Zn are labor-intensive, expensive and time-consuming (Palacios-Rojas 2018). To the best of our knowledge, no study has examined the predictive ability of GS methods that incorporate G x E for Zn concentration in maize. Within the framework of the reaction norm model (Jarquín et al. 2014), the potential of GS for Zn using maize inbreds adapted to tropical environments were assessed. The objectives of this study were; (i) to evaluate the prediction ability for Zn using an association mapping panel and two bi-parental populations evaluated in three tropical environments, (ii) to assess and compare the predictive ability of different GS models, and (iii) to examine the effects of incorporating G x E on prediction accuracy for Zn.

Zinc association mapping (ZAM) panel
The ZAM panel consists of 923 inbreds from maize breeding programs of the International Maize and Wheat Improvement Center (CIMMYT). The panel represents wide genetic diversity for kernel Zn concentration (Hindu et al. 2018).

Bi-parental DH populations
From the ZAM panel, four inbreds with contrasting Zn concentration were selected and used to form two bi-parental (doubled haploid [DH]) populations. DH1 was derived from the F1 generation of a mating between CML503, a high-Zn inbred (31.21 mg/g) with CLWN201, a low-Zn inbred (22.62 mg/g). DH2 was derived from the F1 generation of a mating between CML465, another high-Zn inbred (31.55 mg/g) with CML451, a moderate-Zn inbred (27.88 mg/g). DH1 and DH2 were comprised of 112 and 143 inbreds, respectively.  (Hindu et al. 2018). In 2017, a randomized complete block design (RCBD) with two replications was used. The rows were 2.5 m long and 75 cm apart and each genotype was grown in a single row plot. All plots were managed according to the recommended agronomic practices for each environment.

Experimental design and phenotypic evaluation
From the ZAM panel and each DH population, four to six plants in each plot were self-pollinated, hand-harvested at physiological maturity, hand-shelled and dried to a moisture content of 12.5%. The bulked kernels from each plot are considered a representative sample and were used in subsequent Zn analyses as described (Hindu et al. 2018).

Genotypic data
Genomic DNA was extracted from leaf tissues of all inbred lines (ZAM panel and DH populations) using the standard CIMMYT laboratory protocol (CIMMYT, 2005). The samples were genotyped using the genotyping by sequencing (GBS) method at the Institute for Genomic Diversity, Cornell University, USA (Elshire et al. 2011;Crossa et al. 2013). The restriction enzyme ApeK1 was used to digest DNA, GBS libraries were constructed in 96-plex and sequenced on a single lane of Illumina HISeq2000 flow cell (Elshire et al. 2011). To increase the genome coverage and read depth for SNP discovery, raw read data from the sequencing samples were analyzed together with an additional 30, 000 global maize collections (Zhang et al. 2015).
SNP identification was performed using TASSEL 5.0 GBS Discovery Pipeline with B73 (RefGen_v2) as the reference genome (Elshire et al. 2011;Glaubitz et al. 2014). The source code and the TASSEL GBS discovery pipeline are available at https://www.maizegenetics.net and the SourceForge Tassel project https://sourceforge.net/projects/tassel. For each inbred, the pipeline yielded 955, 690 SNPs which were distributed on the 10 maize chromosomes. After filtering using a minor allele frequency of 0.05 and removing SNPs with more than 10% missing data, 181,889 (ZAM panel) and 170, 798 (bi-parental) SNPs were used for genomic prediction.

Phenotypic data analysis
For the ZAM panel, broad-sense heritability (H 2 ) across environments was estimated as: where s 2 G is the variance due to genotype, s 2 GE is variance due to genotype x environment, s 2 e is the error variance, l is the number of environments and r is the number of replications using multienvironment trial analysis with R (META-R) (Alvarado et al. 2016). For the DH populations, variance components based on the genomic relationship matrix were computed using BGLR package as implemented in GBLUP (Pérez and de los Campos 2014). An estimate of narrow-sense heritability (ĥ 2 ) for each DH population was calculated as:ĥ 2 ¼ŝ 2 ĝ s 2 g þŝ 2 g whereŝ 2 g is an estimate of the additive genetic variance andŝ 2 e is an estimate of the residual variance.
Correlation coefficients between Zn and environments, descriptive statistics and phenotypic data distribution using boxplots were generated in R (core Team 2018). Line means (genotypic values) for the ZAM panel were estimated as Best Linear Unbiased Estimators (BLUEs) with a random effect for replications nested within each environment. Raw data (values) were used for the DH populations.

Statistical models
Genomic models used in this study were based on the reaction norm model which models the markers (genomic) by environment interaction (Jarquín et al. 2014). This model is an extension of the Genomic Best Linear Unbiased Predictor (GBLUP) random effect model, where the main effects of lines (genotypes), genomic, environments and their interactions are modeled using covariance structures that are functions of marker genotypes and environmental covariates.
In this study, environment is the combination of site and year (site-by-year) and the adjusted means (BLUES) to be used in the genomic prediction models are obtained by fitting the phenotypes y ij as: this linear model represents the response of the j th (j = 1,. . .,J) genotype/line tested in the i th (i = 1,. . .,I) environment and fy ij g as the sum of an overall mean m plus random environmental main effect ½E i iid Nð0; s 2 E Þ; the random genotype effect½L j iid Nð0; s 2 L Þ, the random interaction between the j th genotype and the i th environment ½LE ji iid Nð0; s 2 LE Þ and a random error term ½e ij iid Nð0; s 2 e Þ. From this linear model, N(.,.) denotes a normal random variable, iid stands for independent and identically distributed responses and s 2 E , s 2 L , s 2 LE , s 2 e are the variances for environment, genotype, genotype by environment and residual error, respectively. The model above does not allow borrowing of information among genotypes because the genotypes were treated as independent outcomes.
Thus, models used in this study were derived from the baseline model above by subtracting terms or modifying assumptions and/or incorporating genomics/marker information. A brief description of the genomic models used in this study are given below.

M1. Environment + Line
This model is obtained by retaining the first three components from the baseline model (overall mean, random environment main effect and random line main effect) while their underlying assumptions remain unchanged. [1] Here environments were considered as site-by-year combinations.

M2. Environment + Line + Genomic
Another representation of the random main effect of line L j in the previous model is considering a linear combination between markers and their correspondent marker effects, g j ¼ P p m¼1 x jm b m , such that where b m iid Nð0; s 2 b Þ represents the random effect of the m th (m = 1,. . .,p) marker, x jm is the genotype of the j th line at the m th marker and s 2 b its correspondent variance component.
, is the vector of genetic effects, and follows a normal density with mean zero, and a co-variance matrix CovðgÞ ¼ Gs 2 g with G ¼ XX ' p being the genomic relationship matrix (Lopez-Cruz et al. 2015) that describes genetic similarities among pairs of individuals. In this model, the line effect L j is retained to account for imperfect information and model mis-specification because of potential imperfect linkage disequilibrium between markers and quantitative trait loci (QTL).
M3. Environment + Line + Genomic + Genomic 3 Environment This model accounts for the effects of lines L j ; of markers (genomic) g j , of environments (E i ) and the interaction between markers (genomic) and the environment ðgE ji Þ. The model includes the interaction between markers (genomics) and the environment via co-variance structure (Jarquín et al. 2014). The model is as follows: Where gE ji is the interaction between the genetic value of the j th genotype in the i th environment and gE ¼ fgE where Z g and Z E are the correspondent incidence matrices for the effects of genetic values of genotypes and environments, respectively, s 2 gE is the variance component of gE and # denotes the Hadamard product (element-to-element product) between two matrices.

Model assessment
Models were first fitted to the entire data set to estimate variance components using the R-package BGLR (de los Campos et al. 2010; de los Campos and Perez-Rodriguez 2016). The information generated from the full data analyses was not used as prior information for the cross-validation schemes (CV1 and CV2) used for assessing the prediction accuracy of the different models.
Prediction accuracy assessment using cross-validation Two distinct cross-validation schemes that mimic prediction problems that breeders may face when performing genomic prediction were used (Burgueño et al. 2012). One random cross-validation (CV1) evaluates the prediction ability of models when a set of lines have not been evaluated in any environment (prediction of newly developed lines). In CV1, predictions are entirely based on phenotypic records of genetically related lines. The second cross-validation (CV2) is related to incomplete field trials also known as sparse testing, in which some lines are observed in some environments but not in others. In CV2, the goal is to predict the performance of lines in environments where they have not yet been observed. In this study, CV2 mimics a situation where lines are evaluated in two environments but missing in the third environment. Thus, information from related lines and the correlated environments is used, and prediction assessment can benefit from borrowing information between lines within an environment, between lines across environments and among correlated environments.
In CV1 and CV2, a fivefold cross-validation scheme was used to generate the training and validation sets to assess the prediction ability for Zn within the ZAM panel and each DH population. The data were randomly divided into five subsets, with 80% of the lines assigned to the training set and 20% assigned to the validation set. Four subsets were combined to form the training set, and the remaining subset was used as the validation set. Permutation of five subsets taken one at a time led to five training and validation data sets. The procedure was repeated 20 times and a total of 100 runs were performed in each population. The average value of the correlations between the phenotype and the genomic estimated breeding values (GEBVs) from 100 runs was calculated for the ZAM panel, and each DH population for Zn in each environment and was defined as the prediction ability (r MP ).

Data availability
All models were fitted in R (core Team 2018) using the BGLR package (Pérez and de los Campos 2014). All phenotypic and genomic data can be downloaded from the link: http://hdl.handle.net/11529/ 10548331

Descriptive statistics
Mean values of kernel Zn concentration were estimated for each environment and across environments (Tables 1 and 2). For the ZAM panel, kernel Zn ranged from 14.76 to 39.80 mg/g in Celaya 2012, 15.16 to 42.52 mg/g and 17.05 to 46.52 mg/g in Agua Fria 2012 and 2013, respectively (Figure 1). The highest mean (29.53 mg/g) for Zn was observed in Agua Fria 2013. DH1 had Zn values ranging from 16.00 to 48.00 mg/g in Celaya 2012, 16.00 to 35.00 mg/g in Tlaltizapan 2015 and 15.50 to 39.00 mg/g in Tlaltizapan 2017, while the respective values for DH 2 were 17.70 to 43.14 mg/g, 15.60 to 37.80 mg/g and 14.70 to 37.60 mg/g (Figures 2A and 2B). The highest means for Zn were observed in Celaya 2014 (25.38 mg/g) and 2017 (27.96 mg/g) for DH1 and DH2, respectively (Table 2). Across environments, H 2 for n■  (Table 2). There were significant positive correlations between environments for Zn (Table 3), accounting for the moderate to high heritability estimates. Principal component analysis for the ZAM panel suggested presence of a relatively diverse set of lines, and 452 principal components (PCs) were needed to explain 80% of the genotypes' variance ( Figures  3A and 3B). The first two principal components explained 3.85% of the total variance. For the DH populations first two eigenvectors separated the two groups (DH1 and DH2) and 56 principal components were needed to explain 80% of the genotypes' variance ( Figures 3C and 3D). The first two principal components explained 27.50% of the total variation for the DH populations.

Estimates of variance components
Variance component estimates for all models were derived from the full-data analysis (Table 4). In the ZAM panel, variance components for M1 ranged from 16.18 to 7.01 with the main effect of environments explaining the largest proportion (46%). When marker information was incorporated (i.e., M2 and M3), the estimated variance due to environments was reduced to 9.46 for M2 and 8.11 for M3. Inclusion of the interaction term (genomic x environment) reduced the residual variance component by 30%, from an estimated residual variance of 7.01 in M1 to 4.92 in M3 suggesting that some components of differences among genotypes cannot be fully captured by the main effects of markers and environments. Similar trends were observed in DH1 and DH2 except for M1 of DH1 where the main effect of lines accounted for a slightly higher proportion of the total variance (38%) than the main effects of environments (35%).
The residual variance component values for all models of DH populations were slightly higher than those for the ZAM panel (Table 4), possibly because the populations were evaluated in singlereplicated experiments at Celaya and Tlaltizapan (2014 and 2015, respectively). However, estimates from M3 of the DH populations suggest that $ 60 of the within-environment variability can be explained by main effects of lines, genomic (markers) and their interaction term. The proportions of within-environment variation explained by the interaction term were $ 10% for the ZAM panel and DH populations suggesting the importance of considering such interactions in genomic prediction models.

Prediction ability in different populations
Cross-validated r MP values for kernel Zn were estimated for the ZAM panel and DH populations (Tables 5, 6 and 7). The average r MP values in CV1 were consistently lower than those in CV2, suggesting the importance of using information from correlated environments when predicting performance of inbred lines. The mean r MP values in CV1 and CV2 for the ZAM panel were 0.39 and 0.71, respectively (Table 5). For the DH populations, average r MP values were 0.53 for DH1-CV1, 0.44 for DH2-CV1 (Table 6), 0.70 for DH1-CV2 and 0.51 for DH2-CV2 (Table 7).
In the ZAM panel, the highest values in CV1 (0.47) and CV2 (0.72) were obtained in Celaya and Agua Fria 2012 (Table 5). For the bi-parental populations, both under CV1 and CV2, higher r MP values were observed for DH1 compared to DH2. The highest values in CV1 (0.56) and CV2 (0.71) were observed in Tlaltizapan 2017 and 2015, all for DH1 (Tables 6 and 7). The consistently higher r MP values in CV1 and CV2 of DH1 could be attributed to the higher (0.58 to 0.62) correlation values between environments (Table 3).

Prediction ability of different models
Comparing the r MP values obtained from each model, M1 had the lowest (-0.001, -0.03 and 0.04) accuracies in CV1 for the ZAM panel and DH populations (Tables 5 and 6). Those values were improved in CV2 because the predictions benefited from previous records  (collected from other environments) of lines whose Zn values were being predicted. When M1 was expanded to M2 by adding the main effects of markers, the r MP values at each environment and across environments were increased. For example, in CV1, M2, .100-fold increase in r MP values were observed for the ZAM panel and DH populations, and in CV2, M2, average r MP values increased by 2.98%, 2.94% and 11.11% for the ZAM panel, DH1 and DH2, respectively (Tables 5, 6 and 7).
The multi-environment model (M3), which includes the interaction between markers (genomic) and the environment ðgE ji Þ gave higher prediction accuracy than single-environment models (M1 and M2). In CV1, mean r MP values increased from 0.37 (M2) to 0.39 (M3) for the ZAM panel and from 0.43 (M2) to 0.44 for DH2 (Tables 5 and 6). Similar trends were observed in CV2 for the ZAM panel and DH2 (Tables 5 and 7). However, in both CV1 and CV2 of DH1, incorporating gE ji did not improve r MP values for Zn (Tables 6 and 7). For CV1, M3, r MP values for Zn in individual environments ranged from 0.34 to 0.47 for the ZAM panel (Table  5), 0.51 to 0.55 for DH1 and 0.35 to 0.50 for DH2 (Table 6). For CV2, M3, those values ranged from 0.69 to 0.72 for the ZAM panel, 0.68 to 0.70 for DH1 and 0.43 to 0.56 for DH2 (Tables 5,  6 and 7).

DISCUSSION
Overall, moderate to high prediction ability values for kernel Zn were observed for the ZAM panel and DH populations. This could be attributed to the heritabilities observed for kernel Zn (Tables 1 and 2). Similar observations were reported for Zn concentration in wheat (Velu et al. 2016;Manickavelu et al. 2017). Higher predicted values with high accuracy for GS programs are expected for traits with moderate to higher heritability estimates (Combs and Bernardo 2013;Lian et al. 2014;Muranty et al. 2015;Saint Pierre et al. 2016;Manickavelu et al. 2017;Zhang et al. 2017bZhang et al. 2019Arojju et al. 2019). Consistent with a study on Zn and iron (Fe) concentration in spring wheat, the prediction accuracies in this study are sufficient to discard at least 50% of the inbreds with low-Zn concentration (Velu et al. 2016).
Additionally, the moderate to high prediction accuracies reported in this study shows that GS can be used in maize breeding to improve kernel Zn concentration. Assuming two possible seasons of Zn evaluation per year, the predicted genetic gains can be estimated from prediction accuracies and genetic variances of the training populations. The genetic variances for the ZAM panel, DH1 and DH2 were 12.38,12.20 and 14.88, and prediction accuracies were 0.71, 0.70 and 0.51, respectively. If the inbreds in each predicted population are ranked based on their predicted Zn values and the top 10% selected, then their expected average Zn values can be estimated from the proportion of inbreds selected, their respective training population genetic variances, prediction accuracies and the time interval for evaluating the lines. With reference to this, the expected average values of Zn are approximately 31 mg/g for the ZAM panel, 30 mg/g for DH1 and 27 mg/g for DH2. These averages are higher n■  than the averages of the respective training populations (27 mg/g for the ZAM panel, 25 mg/g for DH1 and 26 mg/g for DH2) suggesting that the prediction accuracies achieved are sufficient to select at least 10% of the predicted inbreds with higher Zn concentration.
Data from both bi-parental populations and diverse collection of inbreds have been used for GS and based on cross-validation (CV), it has been established that prediction accuracies could also be affected by the relatedness between training and prediction sets (Habier et al. 2007;de Roos et al. 2009;Asoro et al. 2011;Daetwyler et al. 2013;Cericola et al. 2017;Crossa et al. 2017). In this study, average predicted accuracies were higher for CV1 of the bi-parental populations (0.53 for DH1 and 0.44 for DH2) compared to the ZAM panel (0.39). Higher predicted values in CV1 of the DH populations could be attributed to the closer relationship among DH lines in the training and prediction sets, maximum linkage disequilibrium (LD) between a marker and a QTL, and controlled population structure (Bernardo and Yu 2007;Albrecht et al. 2011;Zhang et al. 2015). In collections of diverse inbreds, prediction accuracy may depend on the ancestral relationships between the lines. So, in experiments using such collections of lines, prediction accuracies have been more variable than accuracies achieved using bi-parental populations (Spindel and McCouch 2016).
Cross-validation (CV) schemes are used in genomic prediction to estimate the accuracy with which predictions for different traits and environments can be made (Burgueño et al. 2012;Zhang et al. 2015;Saint Pierre et al. 2016;Velu et al. 2016;Sukumaran et al. 2017aSukumaran et al. , 2017bMonteverde et al. 2018;Roorkiwal et al. 2018). In this study, two CV schemes (CV1-predicting the performance of newly developed lines, and CV2-predicting the performance of lines that have been evaluated in some environments, but not in others) were used. The utility of these schemes indicated that prediction values for newly developed lines (CV1) were generally lower (0.39 for the ZAM panel, 0.53 for DH1 and 0.44 for DH2) than the values for lines which have been evaluated in different but correlated environments (CV2; 0.71, 0.70 and 0.51 for the ZAM panel, DH1 and DH2, respectively). Such observations indicate the importance of using information from correlated environments when predicting the performance of inbred lines. However, selection of new lines without direct field testing, as simulated in CV1, may enhance the breeding process by replacing the time and labor-intensive field testing for Zn with genomic-estimated breeding values. But, the prediction accuracy values obtained may be n■ E =Environment, L = Line, G = genomic (marker information), G x E = genomic x environment. a Relative to the total variance minus the variance due to main effect of the environment. The percentages of within-environment variance were computed without taking into account the variance of the environment.
n■ lower such that the annual rate of genetic progress in a GS program is compromised (Burgueño et al. 2012). So, the ultimate decision of how a breeding scheme should be structured could depend on the compromise between the desired prediction accuracy and the generation interval (Burgueño et al. 2012). Genotype by environment interaction is an important factor affecting kernel Zn concentration in maize and genomic prediction models that incorporate G x E may enhance the potential of GS for biofortification breeding. For different crop species and traits, genomic prediction models which incorporated G x E achieved higher prediction accuracies in both CV1 and CV2 schemes relative to models which did not include G x E (Burgueño et al. 2012;Guo et al. 2013;Jarquín et al. 2014;Lopez-Cruz et al. 2015;Zhang et al. 2015;Monteverde et al. 2018). In this study, the impact of modeling G x E variance structures for multi-environment trials was investigated and results indicated that the average predicted values from M3 (G x E model) were higher (0.39 and 0.44 for CV1 and 0.71 and 0.51 for CV2) than the values from M2 (non-G x E; 0.37 and 0.43 for CV1-M2, 0.69 and 0.50 for CV2-M2) for the ZAM panel and DH2. These findings agree with those reported on Zn concentration in wheat (Velu et al. 2016), providing evidence that incorporating G x E in GS models can enhance their power and suitability for improving maize for kernel Zn concentration. Conversely, the average predicted values for CV1 and CV2 of DH1 were higher in M2 (0.53 and 0.70) than in M3 (0.53 and 0.69). Except for differences in population size (112 lines vs. 143 lines), this was unexpected since DH1 and DH2 were grown in the same environments.
The gains in prediction accuracies for the GS model that accounted for G x E were dependent on the correlation between environments and CV method used. In this study, the phenotypic correlations between environments were all positive (ranging from 0.58 to 0.62 for DH1, 0.29 to 0.46 for DH2 and 0.61 to 0.66 for the ZAM panel). Such correlations can be exploited using multienvironment models to derive predictions that use information from across both the lines and environments (Burgueño et al. 2012). For instance, although the phenotypic correlations between environments for DH2 were positive (0.29 to 0.46), the lowest average prediction value (0.51) for CV2 was observed for this population. This was expected because CV2 uses phenotypic information from genotypes which have already been tested; hence, effectively exploiting the correlations between environments (Burgueño et al. 2012;Jarquín et al. 2014;Crossa et al. 2015;Pérez-Rodríguez et al. 2015;Saint Pierre et al. 2016;Monteverde et al. 2018). However, for CV1, the information between environments could only be accounted for through the genomic relationship matrix (Monteverde et al. 2018). Hence, the gains in CV1 may likely attribute to more accurate estimate of environment-specific marker effects (Guo et al. 2013). In contrast, when multiple environments are weakly correlated, prediction accuracies from across environment analyses can be negatively affected relative to prediction accuracies within environments (Bentley et al. 2014;Wang et al. 2014;Spindel and McCouch 2016). Thus, before designing a GS experiment, identifying correlated environments where environments can differ in terms of site, year or season in which data were collected is of great interest (Spindel and McCouch 2016).
However, the prediction accuracy values were of lower quality when genomic predictions were conducted across populations. For instance, when the ZAM panel was used as the training population, prediction accuracies for DH1, DH2 and DH1+DH2 were 0.15, -0.10 and 0.09, respectively. When DH1 and DH2 were used as a training and prediction set for each other, prediction accuracies were 0.08 and 0.16 (Unpublished data). These prediction accuracies are n■ Table 6 Correlations (mean 6 SD) between observed and genomic estimated breeding values for Zn in the three environments for three GBLUP models for cross-validation CV1 of DH populations n■ Table 7 Correlations (mean 6 SD) between observed and genomic estimated breeding values for Zn in the three environments for three GBLUP models for cross-validation CV2 of DH populations

Population Environment
Prediction accuracy in CV2 considerably lower than those reported in this study and the differences may be attributed to: (i) weak genetic relationships between the training and prediction population sets and (ii) different methods of analysis because the prediction accuracies reported in this study were partly achieved by modeling the random-effects environment structure to account for G x E while for the unpublished data, the randomeffects environment structure of G x E was not included. The ability to predict kernel Zn concentration using highthroughput SNP markers including G x E interactions creates an opportunity for efficiently enhancing Zn concentration in maize breeding programs. For instance, during early generations of a breeding program, GS can be utilized to identify genotypes with favorable alleles when numbers of progenies and families are large. This could potentially reduce the resource-intensive evaluation process and advancement of false-positive progenies (Velu et al. 2016). Coupled with advances in technologies for assessing Zn, plant scientists can more rapidly measure Zn concentration in maize kernels using the energy dispersive x-ray fluorescence (XRF) assays (Guild et al. 2017). Thus, with more validations and model refinements, GS can potentially accelerate the breeding process to enhance Zn concentration in maize for a wider range of environments.