Genomic prediction for targeted populations of environments in oat (Avena sativa)

Plant material, multi-environment trial data and genomic information

This study used grain yield from MET historical data of 5220 oat genotypes that were grown from 1997 to 2021 in 59 locations across nine states (Illinois, Indiana, Iowa, Michigan, Minnesota, New York, North Dakota, South Dakota, Wisconsin) in the Midwest United States (Supplementary Table S1). This dataset was compiled from different types of plant breeding trials including preliminary yield trials (PYT, F7), or stage 1, advanced yield trials (AYT, F8), or stage 2, elite yield trials (EYT, F9+), or stage 3, early maturity trials (ET) and late maturity trials (LT); several collaborative nurseries including the uniform early oat performance nursery (UEOPN), uniform mid-season oat performance nursery (UMPON), and Midwest collaborative nursery (MCON); and a few large experiments including the mega-environmental design (MED) and Metabolomics (META) studies. Genotypes were evaluated on partially replicated or augmented experimental designs (MCON, META), or with different number of replications in alpha latices (AYT, PYT, MED), or with a randomised complete block design (EYT, ET, LT, UEOPN, UMOPN), depending on the experiment type (available at T3/Oat database, https://oat.triticeaetoolbox.org). This dataset, typical of plant breeding evaluations, is highly unbalanced because not all genotypes were evaluated across all combinations of locations and years.

The study also used the genotypic data from the ‘GBS POGI’ genotype protocol (available at https://oat.triticeaetoolbox.org/breeders_toolbox/protocol/8) from T3/Oat database (available at https://oat.triticeaetoolbox.org), consisting of 1189 genotypes with genotyping-by-sequencing derived single-nucleotide polymorphisms (SNPs). After filtering for a minor allele frequency >0.05 and missing data <20%, there were 3034 high-quality SNPs.

Phenotypic data analyses for GGE biplot and ME delineation

Given the high level of sparsity in the dataset and the fact that proper two-step models tend to perform similarly to single-step models but are less computationally intensive and do not have as many convergence problems (Piepho et al. 2012), a two-step approach was used for obtaining best linear unbiased estimates (BLUEs) of grain yield for each genotype at each location for GEI characterisation. In the first step, genotypic BLUEs for grain yield were obtained after controlling for multiple experiments within a location–year combination using linear mixed models, with experiment within location–year combination modelled as random effects. In the second step, genotypic BLUEs for grain yield in each location–year combination were modelled using the following linear mixed model in the asreml-r (Butler et al. 2020) package of R statistical software (R Development Core Team 2017):

where y_ij is the genotypic BLUE for grain yield of the ith genotype in the jth environment (i.e. location–year combination), μ is the overall mean or intercept, g_i is the effect of the ith genotype, e_j is the effect of the jth environment, ge_ij is the effect of the interaction of the ith genotype evaluated in the jth environment with ge ~ N(0,∑), and ε_ij is the residual error with ε ~ N(0,R). Here, ∑ is the variance–covariance matrix of the GEI, with Σ( [ gxe ]x [ gxe ])=ΣG⊗ΣE where ∑_G(g×g) is the realised additive relationship matrix estimated with the molecular markers using the VanRaden (2008) method in the package AGHmatrix (Amadeu et al. 2016) in R statistical software, and ∑_E(e×e) is the genotypic variance–covariance matrix among environments modelled as a factor analytic of order 1 (FA1); R_{([g×e]×[g×e])} is the variance–covariance matrix of residuals and was modelled as a diagonal matrix with the reciprocals of the standard error of the BLUEs for each genotype in each environment to account for the heterogeneity in the precision of mean estimates in the two-step model following Piepho et al. (2012). However, not all genotypes were used in this second step because of the high level of sparsity. Therefore, a subset of 2000 genotypes that had been evaluated in at least four environments was included in the second step. This ended up accounting for ~82% of all data information. Of this subset, 401 genotypes had phenotypic and genotypic information. Variance components were estimated from a modification of Model 1, where genotype × environment was modelled fully by genotype × location × year interaction (see Table S2 for a full description of the model).

Mega-environments were then delineated by GGE (genotype + GEI) biplot analysis from the grain yield BLUEs for each location using the package GGEplot (Dumble et al. 2017) in R statistical software. Basically, a singular value decomposition of the genotypic main effect and GEI matrix was conducted (Yan et al. 2000; Yan and Tinker 2006) and then it was graphically represented in a biplot (Gabriel 1971; Yan and Tinker 2006). A GGE2 biplot was used for this representation. Winning genotypes were then used to delineate MEs from the GGE2 model (Yan and Tinker 2006). For this analysis a modification of Model 1 was used (see Table S2 for details).

Genomic prediction models

Genomic predictions were accomplished using small modifications from the general GBLUP model described in Model 1 with a full description of all of the prediction strategies in Table S2. When no explicit GEI was modelled, the general model was denoted GBLUP_M following the nomenclature of Lado et al. (2016). On the other hand, when GEI (ge_ij) was modelled with Σ( [ gxe ]x [ gxe ])=ΣG⊗ΣE, with some type of variance–covariance structure for ∑_E, the model was denoted as GBLUP_G×E following the Burgueño et al. (2012) model and the nomenclature of Lado et al. (2016). GBLUP_G×E assumes known environments where the variance–covariance relationship matrix can be estimated.

Genomic prediction cross-validation schemes

Several cross-validation schemes were implemented to address different prediction goals and scenarios. In all the cross-validation scenarios, a portion of the phenotypic dataset is masked, the remaining data are used to predict their performance, and the predictive ability, or correlation between the predicted and observed values, is estimated. One of the interesting prediction scenarios is the prediction of the performance of genotypes that have not yet been evaluated in any environment. Burgueño et al. (2012) called this strategy a CV1. We used the CV1 strategy with a 10-fold cross-validation scheme. Each iteration of the 10-fold cross-validation consists of (1) randomly splitting the dataset into 10 subsets; (2) using nine of the subsets to predict the remaining subset; (3) repeating the previous step multiple times for predicting the remaining nine subsets; and (4) estimating the correlation between observed and predicted values for the full set. The CV1 strategy was used for an overall prediction of the dataset and for predicting the performance in each ME using 100 iterations.

Another prediction scenario consisted of the leave-one-environment-out prediction strategy that Jarquín et al. (2017) called CV0. In this case, all of the phenotypic data for an environment at a time are masked and data from the remaining environments used for the prediction of the new environment. We evaluated six different scenarios with the CV0 approach. First, within MEs the GBLUP_M model was used to predict the new or masked location (CV0_LOC|M). With this strategy, we assume that GEI within MEs is very low and that locations can be properly grouped in MEs even before new performance data are available. Second, within MEs the GBLUP_G×E model was used to predict the new or masked locations (CV0_LOC|GEI). This situation assumes that the GEI within MEs, although small, is known and can be used to model the variance–covariance structure. Our next scenarios consisted of the leave-one-year-out strategy within each ME. For these scenarios, we used the GBLUP_M model assuming unknown genotype × year relationships. We then predicted either each year within a location (CV0_Y|LOC), or each year within an ME (CV0_Y|ME). The predictions by location were performed only for the three largest locations in each ME (i.e. Minot ND, Wayland NY and West Lafayette IN for ME1; Crookston MN, Ithaca NY and Morris MN for ME2; and Brookings SD, Madison WI and Urbana IL for ME3; see Table S1). Using only the largest locations may create an upward bias in predictive ability due to the population size effect. On the other hand, in the predictions for the whole ME, both genotype × year and genotype × location will affect the predictive ability. Although year-to-year variability is random and years cannot be reproduced, an understanding of the stability of the predictions from year to year can be valuable in the decision process of choosing genotypes. Finally, we used overall predictions (instead of by ME) for the CV0_LOC|M and CV0_Y|ME scenarios. All of these models and strategies are fully described in Table S2, and in all cases, predictive ability was evaluated as the correlation between predicted and observed genotypic values.

A third prediction scenario was the use of different numbers of previous years to predict the current year. We use a sliding-window approach where we consider periods of 3, 5, 10, 15, and 20 years to predict the genotypic performance, moving one year at a time in the range between 1997 and 2021 (i.e. we use 1997, 1998 and 1999 to predict the year 2000). We used a modification of Model 1 to account for the year, location and interactions. A full description of this model is presented in Table S2. The predictive ability was evaluated as the correlation between predicted and observed genotypic values for each year.

For all of the genomic predictions with cross-validations, a subset of 401 oat genotypes of the 5220 with 3034 high-quality GBS-derived SNP markers (minor allele frequency >0.05 and missing data <20%) was used. Additionally, predictions for specific MEs were only performed for ME1–ME3 because ME4 consists of a single location and therefore no cross-validation schemes were evaluated for this ME.

Genomic prediction model for variety decision tool

We also predicted the performance of all 1180 genotypes for which marker information was available for the variety decision tool. Genotypes had a large range of phenotypic information available in the dataset, ranging from no phenotypic evaluation to evaluation in 53 environments. The GBLUP_M model was used for this purpose in the asreml-r package of R statistical software. Prediction error variance (PEV), accuracy and reliability were estimated for each genotype following Gezan (2020). PEV was estimated as the squared standard error of the genotypic effects (g^i) as follows:

Accuracy (h) was estimated as the correlation between true and predicted (COR(g,g^)) effects:

Reliability (r) was estimated as the squared accuracy:

The variety decision tool has two goals: to serve farmers and to serve oat breeders in the Midwest US. The identification of superior oat genotypes for each TPE can aid oat breeders in the identification of release candidates for a specific region or in the identification of parents for future crosses. On the other hand, by filtering the dataset to genotypes that have been released (i.e. varieties) and that are available for sale in the Midwest (surveyed in November 2020), the decision tool can aid farmers in the identification of superior varieties for their local environment. Therefore, both the top 10 genotypes and the top 10 varieties were identified for each TPE. The variety decision tool is hosted on a web server by the Practical Farmers of Iowa (https://decide.practicalfarmers.org/).

Phenotypic data summary

This study evaluated 5220 genotypes for grain yield across 59 locations in 25 years resulting in 665 location–year (environment) combinations. BLUEs for grain yield ranged from 74.9 to 798 g m⁻² with a median of 416 g m⁻² across locations (Table S1), and substantial GEI was found (Fig. 1).

Fig. 1.

Genotype × environment interaction characterisation and mega-environment (ME) delineation for a historical dataset for grain yield of oat genotypes evaluated in 59 locations in nine states in the Midwest United States from 1997 to 2021: (a) GGE biplot for grain yield showing grouping of the 59 locations into four MEs; (b) smoothed map of the four oat MEs for grain yield; (c) variance decomposition of the phenotypic variation for grain yield; (d) pie chart of the variance decomposition of the phenotypic variation for grain yield.

GEI characterisation

Mega-environment delineation using GGE biplot analysis of grain yield BLUEs for each of the 59 locations explained 58.7% of the observed phenotypic variation (Fig. 1a). Four MEs were identified based on the winning genotypes on the genotypic polygon (Fig. 1a, b). ME3a and ME3b had different winning genotypes but similar overall ranking of the genotypes and were merged into a single ME3 for further analysis. The percentage of the variance explained by the GEI was 1.85 times higher than the genotypic variance (Fig. 1c, d).

Genomic prediction with different cross-validation strategies

The predictive ability for grain yield for un-phenotyped genotypes in the 10-fold CV1 strategy within MEs ranged from 0.45 (ME1) to 0.65 (ME3) with an overall (all MEs combined) predictive ability of 0.57 (Fig. 2, Table S3). When using CV0 for leave-one-location-out within ME with the GBLUP_M modelling approach (CV0_LOC|ME), ME2 and ME3 had the highest predictive ability with an average of 0.55 (range 0.20–0.89 for ME2 and 0.16–0.82 for ME3, Fig. 2), followed by ME1 with an average of 0.28 (range 0.15–0.51). When using the GBLUP_G×E modelling approach in CV0 for leave-one-location-out (CV0_LOC|GEI), ME3 had the highest predictive ability with an average of 0.55 (range 0.14–0.87, Fig. 2), followed by ME2 with an average of 0.47 (range 0.23–0.85) and ME1 with an average of 0.36 (range 0.19–0.32).

Fig. 2.

Predictive ability within and across mega-environments (MEs) for different cross-validation strategies for an oat dataset consisting of genotypes evaluated in 59 locations in the Midwest United States over 25 years. CV1 used a 10-fold cross-validation strategy for predicting the performance of genotypes that have not been evaluated in any environment within an ME (ME1–ME3) or overall. CV0_LOC|M used leave-one-location-out with an average GBLUP modelling approach (GBLUP_M) within each ME. CV0_LOC|GEI used leave-one-location-out with a GEI GBLUP modelling approach (GBLUP_G×E) within each ME. CV0_Y|LOC used leave-one-year-out with a GBLUP_M modelling approach within each location. CV0_Y|ME used leave-one-year-out with a GBLUP_M modelling approach within each ME. The predictive ability was estimated as the correlation between predicted and observed values.

The predictive ability for CV0 for the leave-one-year-out model was 0.29 for ME1 (range 0.20–0.55, Fig. 2), 0.29 for ME2 (range 0.10–0.62) and 0.25 for ME3 (range 0.12–0.50) when predicting years within ME (CV0_Y|ME, Fig. 2). When predicting for locations within MEs (CV0_Y|LOC) we obtained higher predictive values with a wide range of variation: 0.52 (0.10–0.83) for ME1, 0.51 (0.05–0.95) for ME2, and 0.53 (0.20–0.96) for ME3.

Non-repeatable GEI and years

When using the previous years to predict the next year, predictive abilities using 3, 5, 10, 15 and 20 years were similar in most cases. The overall predictive ability using the most recent 3 years was highest in 45% of the years and within 20% of the highest predictive ability in 20% of the years (Fig. 3). The three most recent years had highest or within the highest predictive ability in 75% (65% highest, 10% not different) of the years for ME1, in 75% (50% highest, 25% not different) of the years for ME2, and in 80% (40% highest, 40% not different) of the years for ME3. The predictive ability using >3 years tended to be better when the general predictive ability was low, such as in 2005, 2013 and 2021, possibly associated with an unusual weather or stress year.

Fig. 3.

Predictive ability using previous 3, 5, 10, 15 or 20 years to predict the current year across (overall) or within mega-environments (ME1, ME2, ME3) for an oat dataset where genotypes were evaluated in 59 locations in the Midwest United States over 25 years: (a) overall predictions, and predictions for (b) ME1, (c) ME2 and (d) ME3.

Genomic prediction for the variety decision tool

Although the accuracy and reliability of released varieties were higher than those of the advanced breeding lines or general genotypes, the actual predicted values were lower (Tables 1 and 2). None of the 10 genotypes differed >2% from the top-yielding genotype in the overall and ME2 predictions of genotypes (Table 1); meanwhile, the first six or four genotypes were within 2% of the top-performing genotype in ME1 and ME2, respectively. On the other hand, there were six commercial varieties with <2% difference from the top-performing variety in the overall predictions, four in ME1, and two in ME2 and ME3 (Table 2).

Genomic prediction success and strategies

In oat, genomic prediction models have been used to assess prediction accuracy showing potential to improve disease resistance (Haikka et al. 2020a), seed quality (Campbell et al. 2021; Hu et al. 2021), grain yield (Haikka et al. 2020b), nutritional compounds and other metabolites (Brzozowski et al. 2022a, 2022b, 2023), and other quantitative traits (Asoro et al. 2011), as well as testing optimisation (González-Barrios et al. 2019). Our study has looked at other aspects of genomic prediction including the prediction for specific TPEs, and we have achieved relatively high predictive abilities of ~0.45–0.65 for new genotypes and 0.15–0.89 for new environments for a trait as complex as grain yield.

Deploying big data in genomic selection for predicting for the TPE

Genotype × environment interactions create many challenges for breeding programs, some of which are easier to deal with than others (Piepho et al. 2008). Under GEI, the efficiency of the prediction relies on the representativity of the tested environments and how they are used for the construction of the training population (Haile et al. 2020; Atanda et al. 2022; Jarquin 2022). Genomic selection can take advantage of large historical datasets by connecting the phenotypic information from different genotypes through borrowing of information from relatives and thus controlling the sparsity (Bernardo 2010; Dawson et al. 2013; Lado et al. 2016; Monteverde et al. 2018) and producing higher predictive ability (Endelman et al. 2014; González-Barrios et al. 2019; Hoefler et al. 2020; Atanda et al. 2021a, 2021b). Therefore, the use of big data consisting of large historical datasets where individual genotypes, or their relatives, have been evaluated can substantially improve the identification of superior individuals that will perform better in the TPE.

Large MET historical datasets have their own constraints (i.e. sparsity and breeding structure) that need to be addressed in order to make efficient use of the available databases. How best to handle the level of sparsity is an ongoing area of research given that genotypes are not missing at random in the dataset because poorer performing genotypes are dropped every year (Aguate et al. 2019; Hartung and Piepho 2021; Wang et al. 2023). Other relevant limitations for the use of large datasets are due to the complexity of the data, the statistical models, and the requirement of great computational power (Atanda et al. 2021b). This study utilised a highly unbalanced dataset from several different METs across 25 years to investigate the use of a mixed model approach in characterising GEI, identifying MEs associated with the TPEs, and determining the potential of genomic selection of grain yield. Overall, the results obtained in the study show the high predictive abilities of large-scale historical data in genomic prediction.

The management of GEI and information leverage is determined by the size and overlap of the ME and TPE areas. MEs can be defined based on the breeding program necessities, adopting different spatial scales (i.e. at a regional scale, González-Barrios et al. 2019; country-wide scale, Crespo-Herrera et al. 2021; or worldwide, Braun et al. 1997). We used defined MEs by providing a broad idea of the GEI in the Midwest region of the United States. Although the TPE might be smaller or larger than the MEs, using the ME definition can be beneficial. For example, the use of information from other locations can still aid in the prediction for TPEs that are smaller than the MEs by increasing the number of evaluations in environments with low GEI. On the other hand, in the case where the TPEs are larger than the MEs or go beyond the limits of a single ME, the use ME units will allow, for example, tailoring of the MET evaluations to account for the variation that exists between MEs in the TPE, and the breeder will need to decide on the release strategies for this case. Furthermore, this work is focused on the use of the groups of locations defined as a TPE instead of predicting locations per se. The shift from a specific location to a TPE is grounded in the fact that performance of farmers’ fields will rarely coincide with the evaluations that breeders perform in research stations; therefore, a focus on predicting individual locations might not be useful for the goal of the breeding program unless it is used as an intermediate step in the prediction process. We were able to predict all MEs with a high level of predictive ability. Furthermore, we tested a situation where GEI within MEs was not modelled (CV0_LOC|M) and a situation where GEI within ME was modelled using a variance–covariance matrix that accounts for genetic correlations between locations (CV0_LOC|GEI). The use of the known correlation structure did not improve the predictive ability reported in other studies (Lado et al. 2016). There could be multiple reasons for this, including the level of GEI within the ME, an artifact due to the population size, or the level of unbalance in the population (Bassi et al. 2016; Lado et al. 2018; Berro et al. 2019), or just the fact that specific location-to-location GEI becomes less relevant when the number of locations is very large. Although the mean predictive ability did not improve with modelling the GEI, the variability in the predictive ability for each location decreased with the GBLUP_G×E model. Arguments can be made that the correlation between environments cannot be known for new locations or for farmers’ fields, and therefore, they should not be used in the prediction models. However, if the focus remains on the overall predictions within MEs, modelling the correlation is just another modelling tool that can be deployed, especially when large historical datasets are used where the estimation of those correlations might be more relevant. Another opportunity to improve predictions when no previous information is available is by incorporating environmental covariates explaining the genotypic performance in the ME (Bustos-Korts et al. 2016; van Eeuwijk et al. 2016; Monteverde et al. 2019; Neyhart et al. 2022; Rebollo et al. 2023). An advantage of the use of environmental covariates is that they may help to fine-tune the limits and boundaries of the ME through the definition of growing zones and modelling the crop growth (Heslot et al. 2013; Bustos-Korts et al. 2016, 2021; Monteverde et al. 2019). Although we did not use environmental covariates in this study, high–medium predictive abilities were obtained. However, we believe that the future use of environmental covariates could improve the predictive ability of these models and help in the delineation of MEs.

Although we obtained high levels of predictive ability for all MEs, the performance of genotypes was easier to predict in some MEs than in others. For example, ME1 had consistently low predictive ability (i.e. <0.4); this was expected because the number of locations within the ME was only five. This ME was the only one that benefited from modelling the genetic correlations among locations using the GBLUP_G×E model, where predictive ability increased from 0.27 to 0.36. This is also an artifact of the low number of locations, where the specific pairwise correlations have a larger impact than in the case of the other MEs with dozens of locations, where specific location-to-location interactions are overall less relevant.

Non-repeatable GEI and years

Year-to-year variation is known to explain an important part of the GEI variation in METs and, more importantly, for being nonreproducible. Furthermore, because many plant breeding programs rely on the evaluation of individuals in a small sample of years, genotypes with positive genotype × year interactions in testing years can be advanced when they would not necessarily have a superior performance in average years for the region. Therefore, datasets with more years may reduce the bias associated with the unbalance in METs, improve the variance component estimation (Aguate et al. 2019; Hartung and Piepho 2021), and improve the actual genotypic predictions (Oakey et al. 2016). Our results suggest that increasing the number of years beyond the last 3 years will not improve the predictive ability of the upcoming unknown year. This could be partially explained by the higher connectivity in the genotypic pool between consecutive years. Some years had a low predictive ability regardless of the set of years used in the prediction model. These are years with unusual weather patterns and some abiotic stresses. The odd performance in those years will be very difficult to predict regardless of the amount of previous information we use or have available. ME1 had 20% of years with a predictive ability <0.1, whereas ME3 had only 5% of years with low predictive ability. This can be partially explained by the population size and the frequency of climatic stresses present in the ME.

We also evaluated the predictive ability for specific years within the ME (CV0_Y|ME) and for specific years within locations (CV0_Y|LOC). The predictions for specific locations were higher (0.20–0.90) than the predictions for the whole ME (0.10–0.60). Multiple factors can explain these findings. The locations used to test the year-to-year variability had more data points and specifically more years than the average location in any given ME. Additionally, predictions for specific locations might not be as relevant as predictions for the ME, as previously discussed. Finally, predictions for the ME include both genotype × year and genotype × location interactions, and although in many cases genotype × year interactions are larger than genotype × location, in our case, they were of the same magnitude. Although in a breeding program we are not able to replicate a year or predict the type of year that will come next, evaluation of the predictions is relevant for understanding the stability in the predictions that can be expected from year to year in real future evaluations.

Variety decision tools and prospects

To the best of our knowledge, ours is the first study to deliver a variety decision tool for farmers in the Midwestern United States based on the farm’s location. There is a limited list of genotypes that are relevant to farmers, and those are genotypes that have been released as a cultivar and are currently commercially available. These genotypes are usually connected in METs because they were tested as advanced breeding lines and, later, most probably continued to be tested in regional trials. However, because there is no unified variety testing system in the United States, the list of varieties to be tested in each state is tailored and therefore differs from the lists to be tested in other states. The reports are then usually provided on a location-by-location basis, or in some cases, a state-wide recommendation is provided. How years are dealt in those reports is also state-specific but most of the states do not report any information beyond 3 years of data. This reporting system creates difficulties at several levels. For example, it is hard for farmers to decide which data to rely on for selecting the best performing varieties for their specific farm. This is especially challenging for farmers close to state borders or far from the specific, sparse testing locations for their state. Furthermore, there is no reason for limiting the flow of information across state borders. The use of MEs to delineate TPEs and the grouping of locations improves the accuracy of predictions and allows the distribution of resources in a more efficient way than by location-state evaluation sheets. This helps to create more efficient delivery of higher quality information to farmers. Finally, being able to include the information from advanced breeding lines is also key to improving the prediction of the future released varieties. This can be seen, for example, in the higher accuracy and reliability values for varieties than for advanced breeding genotypes, which may not have as many relatives in the dataset. This is also reflected in the mean differences between experimental genotypes and varieties, where varieties are >2% apart. This decision tool is being created as a first exploratory strategy. It represents an improvement over state reports, but further development would improve the predictions. For example, a very simple strategy was used to delineate the borders of the MEs, but environmental covariates could be used to refine these delineations further following the work of Rebollo et al. (2023) and Bustos-Korts et al. (2021). Additionally, environmental covariates could be used to improve the actual predictions. In summary, the strategy used in this study will assist in developing a tool for selecting the top-ranking genotypes in specific locations, using a zip-code-based system in the Midwest region, via the deployment of the genomic selection model in each TPE. Similar strategies could be expanded to other crops in any region around the world.