Soybean [Glycine max (L.) Merr.] delivers the highest amount of protein per hectare than any crop and accounts for over 60% of total global oilseed production (United States Department of Agriculture, 2022a). It is the largest and most concentrated segment of global agricultural trade and one of the most essential crops to the world’s food security (Gale et al., 2019). Soybean production has nearly doubled over the last two decades (182,830 to 363,860 MT) (United States Department of Agriculture, 2022b). The genetic improvement of soybean cultivars, as well as advancements in farming technology and agronomic practices, have significantly contributed to this substantial increase (Specht et al., 1999; Rowntree et al., 2013; Koester et al., 2014; Rincker et al., 2014; Canella Vieira and Chen, 2021). A typical soybean breeding pipeline consists of several years of multi-environment field trials to select and advance high-yielding breeding lines (Canella Vieira and Chen, 2021; Yoosefzadeh-Najafabadi and Rajcan, 2022). The availability of high-dimensional genomic data (Song et al., 2013; 2020) and the advancements in genome-based prediction models (GP) have revolutionized and contributed to accelerated genetic gains as well as higher testing efficiency in soybean breeding programs (Jarquin et al., 2014a; Jarquin et al., 2014b; Persa et al., 2020; Widener et al., 2021; Canella Vieira et al., 2022).

The concept of GP revolves around using the information of all molecular markers—regardless of estimated effect size of significance—to develop prediction models of the genetic merit for the phenotype of interest in unobserved genotypes (Meuwissen et al., 2001). Thus, GP allows the identification and selection of desirable genotypes earlier in the breeding pipeline, which not only reduce cost, time, and space but enhance genetic gain by shortening the length of the breeding cycle and increasing selection intensity (Jarquin et al., 2014a; Crossa et al., 2017; Canella Vieira and Chen, 2021; Wartha and Lorenz, 2021; Canella Vieira et al., 2022). However, the presence of the genotype-by-environment G×E interaction, a change in the response patterns from one environment to another, complicates the selection of improved cultivars (Crossa et al., 2011). For this reason, it is necessary to establish multi-environment trials (METs) to evaluate the performance of genotypes under a wide range of weather conditions (environmental stimuli) allowing the selection of stable materials or materials with local adaptation only (Jarquin et al., 2020a). As expected, the high phenotyping cost does not permit the evaluation of all candidate genotypes in all of the environments of interest but a fraction of these combinations of genotypes-in-environments (Jarquin et al., 2020b). To overcome these disadvantages (i.e., G×E and the high phenotyping costs), the implementation of the reaction norm model (Jarquin et al., 2014a) leverages the borrowing of information of genotypes across environments helping to increase the predictability of unobserved combinations of genotypes-in-environments.

In addition, GP including G×E model parameters can substantially improve field testing design and efficiency, as well as resource allocation (Jarquin et al., 2020b; Montesinos Lopez et al., 2023a). For instance, GP can reduce the costs and space associated with field testing by using sparse testing designs in which only a subset of the genotypes are tested at each location (Jarquin et al., 2020b). Sparse testing allows the prediction of non-observed genotype-in-environment combinations reducing the costs at a fixed evaluation capacity (less expensive to make accurate inferences on the original set of genotypes-in-environment combinations) or increasing the overall evaluation capacity at fixed costs (inferences on more genotypes and environment combinations based on the original budget) (Jarquin et al., 2020b).

Using two maize (Zea mays L.) data sets from the International Maize and Wheat Improvement Center (CIMMYT)’s breeding program in eastern Africa, Jarquin et al. (2020b) were the first to demonstrate that GP models could substantially reduce the testing footprint of breeding programs using sparse testing designs. Additional studies of GP-based sparse testing designs have been reported in wheat (Triticum L.) (Crespo-Herrera et al., 2021; He et al., 2021; Atanda et al., 2022; Montesinos Lopez et al., 2023a; Montesinos Lopez et al., 2023b), maize (Montesinos Lopez et al., 2023b), groundnut (Arachis hypogaea L.) (Montesinos Lopez et al., 2023a), and rice (Oryza sativa L.) (He et al., 2021; Montesinos Lopez et al., 2023b). To date, no applications of GP-based sparse testing have been reported in soybean. Therefore, the objective of this study is to investigate the potential of reducing the field testing footprint (less natural resources such as land and water associated with the in fields evaluation of RILs) in soybean breeding programs based on sparse testing designs, as well as the prediction accuracy derived from different GP models including the main effect of the molecular markers via covariance structures (M1), a multiplicative reaction norm model (M2) to account for the genotype-by-environment G×E interaction, and an extended reaction norm model also including the family structure (M3) in interaction with environmental stimuli (Persa et al., 2020). Two different methods for model calibration were considered. The first one, is initially based on RILs randomly selected then varying sample sizes and training composition (between non-overlapping [NO-RILs] and overlapping [O-RILs] genotypes across environments) for a fixed testing set size. While the second one, only varies the training set size since it is based on common sets of genotypes observed across environments and selected under a genetic criteria/algorithm, where the goal is to select a core sample of RILs that maximizes/minimizes the genomic diversity on a sample of fixed size. The objectives of implementing this second selection method were to assess the impacts in predictive ability using different levels of genomic diversity of the RILs when calibrating models, and evaluate the stability of these selected RILs across environments. The impacts in predictive ability using these selection methods were evaluated using a soybean population of 1,755 genotypes evaluated in nine environments (all genotypes in all environments).

The implementation of GP has revolutionized commercial and public soybean breeding programs by allowing plant breeders to predict the phenotype of interest in unobserved genotypes (Jarquin et al., 2014a; Jarquin et al., 2014b; Persa et al., 2020; Widener et al., 2021; Canella Vieira et al., 2022). The first report on GP in soybean was based on a standard G-BLUP model including only additive effects and an extended version of the G-BLUP model including additive-by-additive effects (Jarquin et al., 2014b). Ma et al. (2016) used ridge regression best linear unbiased prediction (rrBLUP) (Endelman, 2011) with fivefold cross-validations to explore strategies of marker preselection. Prediction accuracy based on pre-selected markers slightly increased compared with random or equidistant marker sampling (Ma et al., 2016). Xavier et al. (2016) showed that the training population size was the most impactful factor in the prediction accuracy when investigating the impacts of training population size, genotyping density, and different prediction models. Persa et al. (2020) expanded the reaction norm model proposed by Jarquin et al. (2014a) by incorporating the interaction between the family’s membership of the genotypes and the environment under the premise that the differential responses of families to environmental stimuli could be used for enhancing the selection process in target environments. These authors showed significant improvements in predictive ability by incorporating the family-by-environment interaction in the models compared to the conventional reaction norm model (Persa et al., 2020). However, they also pointed out that this model requires to observe, at least partially, the individuals of the family to predict otherwise no improvements in predictive ability are expected when predicting a yet to observe family. Thus, this model is benefited only when some individuals of the target family are observed. In this case, the family term here has a closely related to pedigree information. In addition, Habier et al. (2007) and (2013), showed that the addition of a model term modeling the differences between groups of indiviuals (families, pedigree, etc.) besides the genomic information, does a better job capturing relationship information rather than just relationships alone. In this way, the family term enables the genomic information to focusing on only capturing LD information rather than trying to capture or explain at the same time LD and the relationships.

Canella Vieira et al. (2022) investigated the potential of incorporating soil texture information and its interaction with molecular markers via covariance structures for enhancing predictive ability across breeding scenarios. The obtained results were more stable when predicting trait performance in novel environments compared to the conventional reaction norm model.

The expression of a phenotype is a function of the genotype, the environment, and the interaction between the genotype and environment (G×E) across different environments (de Leon et al., 2016). Grain yield is a highly complex and quantitative trait regulated by numerous large and small-effect genes, of which its expression largely depent on the genotype interaction with various components of the environment. Genomic prediction-based sparse testing is a promising concept to substantially increase the number of tested genotypes and environments while maintaining fixed requirements for land, seed availability, and costs (Jarquin et al., 2020b; Montesinos Lopez et al., 2023a). Hence, GP-based sparse testing enables a more comprehensive assessment of genotype stability and adaptation in the early stages of the breeding pipeline, which likely could result in a more accurate selection of superior genotypes widely adapted to various environments.

In soybean breeding, this is particularly important in the early stages of yield trials (i.e., preliminary yield trials) where seed availability and the excessive number of testing genotypes are major constraints (Canella Vieira and Chen, 2021). In the public sector, preliminary yield trials often consist of thousands (∼1,000 to 3,000) of genotypes tested across several replicated environments (three to seven environments, two to three replications) resulting on average 20,000 preliminary yield plots (2,000 genotypes, five environments, two replications). As shown in this research (1,755 genotypes, nine environments, a total of 15,795 phenotypic records), GP-based sparse testing can sustain high accuracies at various training compositions which could reduce the testing footprint by as much as 90% (Figure 1). Given the cost of approximately $15 for a yield phenotype (plot), this methodology has the potential to reduce the total cost of preliminary trials (depending on the size of the program) by as much as $210,600 ($236,925–26,325; 1,755 vs 15,795) with the largest training set size (195) while about $224,100 ($236,925-$12,825; 855 vs 15,795) for the smallest one (95) while maintaining the number of genotypes and environments fixed. On the other hand, the number of tested RILs and/or environments can increase by 5- to 10-fold while maintaining the costs fixed. Considering the premise of same phenotyping cost across environments (fixed cost of $15 per phenotypic record), with the initial budget dedicated to test all the RILs in all of the environments ($236,925) the number of RILs or the number of environments can be increased by 9 folds or considering combinations of these.

In this study, different prediction models, as well as multiple training set sizes and compositions (ratio of NO-RILs and O-RILs across environments, SSD.max and SSD.min), were investigated. The most comprehensive model, including the interaction between families and the environment (M3), yielded the highest prediction accuracies independent of training set sizes and composition. The ability to borrow information between genotypes derived from the same families likely contributed to the superior performance of M3 because half of the cross was already observed a considerable number of times in other crosses and environments (Jarquin et al., 2020a). In GP-based sparse testing where most genotypes are untested within and across environments, borrowing information from tested individuals within the same family but observed in different environments has been shown to improve predictive ability by as much as 48% compared to models including the interaction of molecular markers and environments only or reaction norm model (Persa et al., 2020). Across all three models, reducing the training set size negatively impacted prediction accuracies. However, prediction accuracy was rapidly recovered by the addition of O-RILs. This indicates that the training set composition is critical to a successful GP-based sparse testing implementation when the training set is reduced. This observation is also important in advanced yield trials where the number of testing genotypes is significantly reduced as compared to preliminary trials (roughly 90% fewer genotypes), and the number of testing environments is significantly increased (Canella Vieira and Chen, 2021). Specifically, the advanced yield trials stage offers more flexibility in maximizing training set sizes and composition (increasing O-RILs) since seed availability should not be a major constraint.

The results and trends in predictive ability obtained with models M1 (main effects only) and M2 (reaction norm considering the interaction between markers and environments) were similar to those obtained by Jarquin et al., 2020b; Crespo-Herrera et al., 2021 analyzing maize and wheat data, respectively. However, these authors did not consider models including the interaction between family and environments. On the other hand, the results of model M3 (also including the interaction between family and environment) were similar to those obtained by Persa et al., 2020 analyzing information of the SoyNAM experiment comprising 1,358 RILs observed in 18 environments (not all RILs observed in all environments) but considering a conventional fivefold cross-validation predicting tested genotypes in observed environments (CV2) and untested genotypes in observed environments (CV1). In addition, Xavier and Habier (2022), using the SoyNAM population conducted an study to evaluate the effects in predictive ability of models similar to M1 and M2 when considering simulated data for the case where RILs are observed only once across environments for different heritabilities. The results obtained for these authors were higher compared with the results here presented probably due to the fact they considered simulated data as response as opposed to real data and a different cross-validation scheme.

In this study, in addition to the ratio of NO-RILs and O-RILs, the selection of RILs to compose the training set based on maximizing (SSD.max) and minimizing (SSD.min) genetic diversity was also investigated. The Super Saturated Design (SSD) methodology was implemented to identify the set of 195 RILs that either decreased or increased E(S ²) (SSD.max and SSD.min, respectively). The method successfully created two fixed-size sample groups with contrasting genetic diversity. This was further observed in the population structure of each group, where SSD.max samples showed wider distribution and minimal clustering as compared to the SSD.min samples (right panel in Figure 2). Interestingly, the SSD.max samples' prediction accuracy outperformed the SSD.min samples in all scenarios (Figure 5). Superior prediction accuracy in the SSD.max samples could be explained by the higher availability of diverse alleles associated with stress resilience, and therefore, resulted in higher stability across environments. A narrower genetic diversity, although may prioritize high-yielding alleles, can be more susceptible to environmental stressors. This was also observed in the higher phenotypic correlation of the SSD.max samples across all environments compared to SSD.min samples (Supplementary Figure S2).

Also, in general the results of the SSD.max sample slightly under performed the results of the random samples for most of the sample sizes of common RILs across environments and across ancestry groups, except for the case when the sample size was fixed in 115 common RILs (Figure 5A). Considering the different ancestry groups, mixed results were obtained for G1, while for G2 better results were obtained with the SSD.max sample, and for G3 the random samples were superior. In all of the cases, the SSD.min returned the worse results. Seems like the G1 group of families (high yielding) is very susceptible to the presence of genotype-by-environment interactions thus reduced sample sizes with a high propotion of non-overlapped genotypes significantly reduced predictive ability. On the other hand, in G2, the reduction of the training set sizes is the only factor that reduces predictive ability, and the results are more stable across the different training sets compositions. Finally, for G3, neither the sample size nor the training set composition seems to affect the predictive ability since similar results were obtained across these.

Additionally, higher genetic diversity (here represented by G3) yielded stable prediction accuracies across various training set sizes and compositions. Like the SSD.max observation, this could be attributed to the higher availability of diverse stress resilient-alleles and therefore require a smaller number of samples, as well as reduced O-RILs to achieve maximum prediction accuracy. These observations are fundamental to the establishment of a successful GP-based sparse testing design and should be further explored across populations with various genetic backgrounds, including high-yielding bi-parental populations.

Similarly to the analysis of the different genetic backgrounds, a more detailed analysis can be done by considering the model’s predictive ability at the family level. As expected, a significant reduction in predictive ability was obtained (Supplementary Figures S4–S6); however, there was not a unique model always outperforming the other two. In general, model M3 slightly outperformed the other two models, specially when considering the largest training set size and with a lower number of overlapping genotypes. Also, it was superior when only considering the families from G3. Arguably the results obtained when computing the correlation across families might be inflated; however, the selection of the superior cultivars are not made within families but across these. In addition, considering the three models, there was not a unique model systematically outperforming the others. The lack of a unique model significanly outperforming the within families predictive ability is due to the fact that no relationship is factored in when ranking individuals within families (Habier et al., 2013). All the models have the same ability to model LD. In our case, at the family level, the most complex model M3 showed slight improvements in predictive ability by including the family term as main effect and in interaction with environments. This could help breeders to better screening RILs across families and environments while reducing phenotyping cost via sparse testing designs.

Genomic prediction-based sparse testing design is a promising approach to further maximize the applications of high-dimensional genomic data and predictive models toward improving cultivar development. The increase of testing environments in the early stages of the breeding pipeline can provide a more comprehensive assessment of genotype stability and adaptation which are fundamental for the precise selection of superior genotypes widely adapted to various environments. Various training set sizes and compositions, as well as prediction models, have been investigated. Overall, the training set size and the inclusion of O-RILs appear to be the main factors impacting prediction accuracy given a fixed training set size while the genetic diversity seems to be a secondary factor, except when it was minimized returning the worse results. Additional studies investigating the real-world effectiveness of prediction accuracy based on genotype ranking and advancement breeding decisions can help determine the ideal protocols for GP-based sparse testing in soybean. In summary, GP-based sparse testing can either improve/increase testing capacity (represented as the number of genotypes and environments) at a fixed cost or substantially decrease the cost of a breeding pipeline at a fixed testing capacity.