Wheat lines used in this study are from the breeding program of Washington State University (WSU) and were grown at various locations in the state of Washington on grower-cooperator fields using common agricultural practices for the region. Grain yield (GY) was collected using a Zürn 150 Combine (Zürn Harvesting GmbH & Co.) and was used for each of the four data sets.

• Dataset 1, Wheat_1 (Year 2019) contains 1,397 unique lines and three environments (Kincaid, Lind, Pullman) and contains 1,869 total observations since some lines are repeated in various environments.

Phenotypic data was collected using the Sentera Quad Multispectral Sensor (Sentera, St Paul, MN), which covered target bands of interest for winter wheat evaluation. The camera has four sensors that cover eight broad spectral bands between 450 and 970 nm. An unmanned aircraft system (UAS) mounted with the Sentera camera flew a programmed route at an elevation of 45 m capturing overlapping georeferenced images. Collected UAS images were stitched and prepped for data extraction in Pix4Dmapper (Pix4D Inc., Denver, CO), creating a single orthomosaic image for each sensor per location. Orthomosaic images were transferred to Geographic Information System (QGIS) for plot identification and then further processed with a custom R code for calibration, index calculation, and single plot mean data extraction. In 2019, a single reflectance panel (85% reflectance) was used for radiometric calibration on red, blue, green, (RBG) and red edge bands (RE1 and RE2). Quantum efficiency coefficients were used to calculate near infrared (NIR) using: N I R = 2.921 × B l u e − 0.754 × R e d

The NIR band was then normalized with a coefficient of 3.07 during the calculation of SRIs. In 2020 through 2022, a set of calibration panels were used (five panels ranging from 2%–85% reflectance, MosaicMill Oy, Vantaa, Finland). All raw band layers were adjusted based on the relationship: S R = D N × S l o p e ± i n t e r c e p t

Where the slope and intercept are based on the regression of the observed reflectance in calibration panels, digital numbers (DN) are the raw observed pixel values, and surface reflectance (SR) is the true reflectance value (Iqbal et al., 2018). All datasets used adjusted multispectral band values to calculate indices for further model analysis.

All the lines were genotyped using genotyping-by-sequencing (GBS; Poland et al., 2012). The original SNPs totaled 6,075,743, but after filtering for SNPs with homozygosity >80%, for less than 50% missing data, greater than a 0.05 minor allele frequency, and less than 5% heterozygosity, we end up with 19,645 SNPs. Markers with missing data were imputed using the “expectation-maximization” algorithm in the “R” package rrBLUP (Endelman, 2011). In each data set, the best linear unbiased estimates (BLUEs) were computed under two experimental designs:

The BLUEs for GY within each environment were calculated using the lmer function of the lme4 package (Bates et al., 2015) of the R statistical software with the following mixed linear model: y i j k l = μ + g i + c h e c k i + t j + r k j + b l j k + ε i j k l where y i j k l is the GY of the ith genotype in the jth trial, kth replicate and lth block, μ is the general mean, g i is the fixed effect of the genotype i, c h e c k i is the fixed effect of the check-genotype i, t j is the random effect of the trial, t j ∼ N I I D 0 , σ t 2 ; where NIID stands for normal, independent and identically distributed, r k j is the random effect of the replicate within the trial, r k j ∼ N I I D 0 , σ r 2 ; b l j k is the random effect of the incomplete block within the trial and the replicate, b l j k ∼ N I I D 0 , σ b 2 ; and ε i j k l is the residual ε i j k l ∼ N I I D 0 , σ 2 .

In this experimental design, the BLUEs for GY within each environment were calculated using the lmer function of the lme4 package of the R statistical software with the following mixed linear model: y i j = μ + g i + c h e c k i + b j + ε i j where y i j is the GY of the ith genotype in the jth block, μ is the general mean, g i is the fixed effect of the genotype i, c h e c k i is the fixed effect of the check-genotype i, b j is the random effect of the jth block, b j ∼ N I I D 0 , σ b 2 ; and ε i j is the residual ε i j ∼ N I I D 0 , σ 2 .

The Bayesian models implemented only differ in the predictor they used. For this reason, the general model is given: Y i j = μ + E T A i j k + ϵ i j (1)

Where Y i j denotes the response variable in the jth line in the ith environment, μ denotes the general mean (intercept), and ϵ i j are random error components assumed to be independent normal random variables with mean 0 and variance σ 2 .   E T A i j k is the predictor for the jth line in the ith environment. The different E T A i j k ´ s (ETA) used are provided in Supplementary Table SA1.

The vector of length eleven H i j 1 , … , H i j 11 of the information include both the raw multispectral data and calculated indices denoted as: Blue, Green, Red, NIR, RE1, RE2, 900, 975 nm, NDRE1, NDVI, and Canopy Cover. While when the vector of length three ( I i j 1 , I i j 2 , I i j 3 ) was used as independent variables including only the indices: Can_Cover, NDRE1, and NDVI were used. The implementation of these models was carried out in the R statistical software (R Core Team, 2022) using the BGLR library of Pérez and de los Campos (2014). Equations used in the calculation of indices can be found in Table 1.

This study utilized the univariate PLS model, a statistical machine learning method introduced by Wold (2001) in econometrics and chemometrics for regression analysis. PLS is very useful for prediction problems where the number of independent variables ( p ) is larger than the number of observations ( n ) and when predictors are highly correlated. Under the univariate PLS framework, the response variable Y is a vector instead of a matrix of order n × 1 that is linked to a set of explanatory variables ( X ) of order n × p (Wold, 2001; Boulesteix and Strimmer 2006). In PLS, instead of regressing Y on X , we regressed Y on T , where T are the latent variables (LVs), also called X -scores or latent vectors; these LVs are related to the original X and Y matrices. The goal of PLS regression is to maximize the covariance between Y and T ; however, an iterative procedure is implemented for its computation. The main steps to compute the LVs under a univariate framework using the kernel algorithm for PLS are:

Step 1

Initialization of E = X and F = Y . Center each column of E and F ; scaling is optional.

Step 2

Compute S = X T Y (Cross product matrix) and then S S T = X T Y Y T X and S T S = Y T X X T Y .

Step 3

Compute the singular value decomposition (SVD) of S S T and S T S .

Step 4

Obtain w and q , the eigenvectors to the largest eigenvalue of S S T and S T S , respectively.

Step 5

Compute scores t and u as t = X w = E w and u = Y q = F q .

Step 6

Normalize the t and u scores as t = t / t T t and u = u / u T u .

Step 7

Next, compute X and Y loadings as p = E T t and q = F T t .

Step 8

Deflate matrices E and F as E n + 1 = E n − t p T and F n + 1 = F n − t q T .

Step 9

Use as input E n + 1 and F n + 1 , of Step 8, in Step 2, and repeat steps 2 to 9 until the deflated matrices are empty or the necessary number of components have been extracted.

With the outputs of w , t , p and q vectors, the matrices W, T, P, and Q, respectively, are built. Finally, after having all the columns of W , we compute R as: R = W P T W − 1

Next, with R we can compute the LVs, which are related to the original X matrix as: T = X R

Next, since we regressed Y on T , the resulting beta coefficients are b = T T T − 1 T T Y . However, to convert these back to the realm of the original variables ( X ) , we pre-multiplied with matrix R the beta coefficients ( b ); since T = X R , B = R   b

To reach optimal performance of the PLS method, only the first a components are used. Since regression and dimension reduction are performed simultaneously, all B , T , W , P and Q are part of the output. Both X and Y are considered when calculating the LVs in T . Thereafter, predictions for new data ( X n e w ) should be done with: Y ^ n e w = X n e w B = X n e w R b = T n e w b where T n e w = X n e w R . In this study, the optimal number of components was determined by cross-validation. We used the NRMSE, with an inner 10-fold cross-validation for selecting the optimal number of hyperparameters.

In this application, we used the concatenation of the different sources of information for each predictor (ETA1 to ETA9 given in Table 1) as the matrix of independent variables X. For this reason, we first computed the design matrices of environments ( X E ) , the design matrix of genotypes ( X g ) and the design matrix of the Genotype × Environments term ( X g E ). But, since the PLS method does not allow direct inclusion, like the Bayesian GBLUP model, the genomic relationship matrix of lines G = M M T r , where M denotes the matrix of markers (coded as 0, 1 and 2) of order J × r ; J denotes the number of lines; r the total number of markers. The design matrices of lines and genotype × environments were post-multiplied by their corresponding square root matrices of their corresponding relationship matrices to incorporate into the design matrix this relationship information. That is, instead of using only X g ; X g E as input, we used X g L g (with L g = G 0.5 )   a n d   X g E L g E (with L g E = G G E 0.5 ) . For this reason, the final input matrix used for ETA1 to ETA9 under the PLS model was; X = X E , X g L g ; X = X E , X g L g , H ; X = X E , X g L g , H , X g E L g E ; X = X E , X g L g , H , X g E L g E , X g L g : H ; X = X E , X g L g , I ; X = X E , X g L g , I , X g E L g E ; X = X E , X g L g , I , X g E L g E , X g L g : I ; X = X E , H ; X = X E , I respectively, where H = H 11 , … , H 1 J , … , H I J T ; I = I 11 , … , I 1 J , … , I I J T ; X g L g : H represents the interaction term between the genomic information and the multispectral information; X g L g : I represents the interaction term between the genomic information and three indices built from the multispectral information. We did not post-multiply the design matrix of environments ( X E ) since we did not compute an environmental relationship matrix with environmental covariates, only with the dummy values of the position of environments. For this reason, under the PLS model were used as input the vector of response variables ( Y ) and the input matrix X , just defined above. The implementation of the PLS models was performed with the R statistical software (R Core Team 2022) using the PLS library (Mevik and Wehrens, 2007).

In each of the four datasets (corresponding to years 2019–2022), for implementing the type of cross-validation partially tested lines in tested environments, we used seven-fold cross validation (7FCV) (Montesinos-López et al., 2022). For this reason, 7 − 1 folds (85.71% of the data) were assigned to the outer-training set and the remaining fold (14.29% of the data) was assigned to the outer-testing set, until each of the 7 folds were tested once. Under the PLS model for tuning the number of principal components required ten nested cross-validations, that is, the outer-training was divided into ten groups where nine were used for inner training set (90% of the training) and one for the validation (inner-testing) set (10% of the outer training). This means that under the PLS model, the data set was divided in outer-testing (14.29% of data), inner-training (77.14% of data), and validation (8.57% of data). Using the validation set, the optimal number of principal components was selected. Also, it is important to point out that the sum of the inner-training plus the validation equals the outer-training. Next, the average of the ten validation folds was reported as the metric of prediction accuracy to select the optimal hyperparameter (number of principal components). Then, using this optimal hyperparameter, the PLS model was refitted with the whole outer-training set (the 7 − 1 folds), and finally, the prediction of each outer-testing set was obtained. For the selection of the hyperparameters under the inner-training, the average mean square error was computed and used as the metric of accuracy, but for the outer-training to evaluate the prediction accuracy under the partially tested lines in the test environments cross-validation, the average Pearson´s correlation was computed. It is important to note that under the GBLUP, not tuning was required and only the outer 7FCV was implemented and the average of the 7 folds was reported as prediction accuracy for each environment using the Pearson´s correlation. But the computation of Pearson´s correlation across environments (Global) under the outer 7FCV was done between averages of true and predicted phenotypes of lines over environments per year (data set). On the other hand, to implement the cross-validation partially tested lines in untested environments we used a leave one environment out (LOEO) approach where the training set was composed of the total number of environments ( n E ) minus one, and the remaining environment was used as testing set, meaning that each environment was used as testing set exactly one time. For this reason, only the average prediction accuracy for each environment was reported since only one-fold (testing set) was obtained for each environment. However, across environments, in addition to the average Pearson correlation, it was also possible to estimate the standard error. The Pearson´s correlation across environments (Global) in LOEO cross-validation was computed averaging the predictions resulting in each of the environments under study in each year. Under this approach the tuning process for the PLS was done exactly as was done under the 7FCV strategy.

We note in Table 2 that the heritabilities of years 2019, 2020 and 2021 are larger than 0.70; however, the heritability for year 2022 is low (0.099). Also, we can see in Table 2 that the GE interaction term is not relevant in years 2019, 2020 and 2022. We note that each of the four data sets is unbalanced since each wheat line was evaluated on average in 1.462, 1.252, 1.711, and 1.309 environments in 2019, 2020, 2021 and 2022, respectively. It is important to note that the environments evaluated in each year were 3, 6, 8, and 6, respectively. It can also be observed in Table 2 that the replications of each line in each environment were less than two since some lines were evaluated in replicated experiments while the remainder were examined in unreplicated experiments.

The results under 7FCV for each year are provided. It is important to note that for each model (GBLUP and PLS), nine predictors were evaluated to see how much each part contributes to improve the prediction accuracy and the results are reported for each environment and across environments (Global) for each year.

In Figure 1 and Supplementary Table SA2 we observe that under the GBLUP model, the best prediction performance was observed in environment Kincaid and the worst in Lind, while under the PLS model, the best predictions were also observed in the Kincaid environment, and the worst in Pullman. In Figure 1 and Supplementary Table SA2, the worst predictions were observed under the predictor (E + g). We observed when both types of information are integrated (genomic + multispectral information under its H or I versions) the best prediction performances were obtained and the predictions across environments and for some environments with Pearson’s correlation values close to one. However, Figure 1 shows that adding the interaction terms gE, gH, and gI to the predictors does not significantly increase prediction performance. Also, it is observed that simple predictors that contain only multispectral information, like E + H and E + I, produce similar performance to predictors that incorporate the genotypic information and interaction terms. However, the predictors that integrate both sources of information were more stable and consistent. Regarding the predictions using H or I, we observe that using I information is better since more consistent results than the H information and larger Pearson’s correlation are observed. In the models, we observe that both models are effective in this prediction problem and with this data, yet the predictions using GBLUP were better in most cases (Supplementary Table SA2).

Under the GBLUP and PLS models, the best predictions were observed in environment Ritzville and the worst in Kincaid; however, very competitive predictions were observed in most environments except for Kincaid (Figure 2; Supplementary Table SA3). Also, the worst predictions were observed under the predictor E + g (see Figure 2; Supplementary Table SA3). The best prediction performances were observed when both types of information are integrated (genomic + multispectral l information under its H or I versions) with Pearson’s correlation values close to one in some environments and across environments (Figure 2). Again, we observe in Figure 2; Supplementary Table SA3 that adding interaction terms gE, gH, and gI in the predictors does not provide a relevant increase in performance. In the 2020 data, we observed simple predictors with only multispectral information like E + H and E + I produce a similar performance to more complex predictors that incorporate the genotypic information and some interaction terms, but predictors with both sources of information are more stable and consistent. Regarding the predictions using E + H and E + I, we observed the predictor E + H produced results closer to one for Pearson’s correlation values. Both models were effective in this prediction scenario and with this data, but the predictions using GBLUP were better (Supplementary Table SA3).

In Figure 3 and Supplementary Table SA4 , we can note that under both models (GBLUP and PLS), the worst predictions were observed in Pullman. In the remaining environments, the predictions were effective since they were close to one in terms of Pearson’s correlation. The worst predictions were observed under the predictor E + g and the best when both types of information were integrated (genomic + multispectral information under its H or I versions) with Pearson’s correlation values also close to one. Again, it was observed in Figure 3; Supplementary Table SA4 that adding interaction terms gE, gH, and gI in the predictors did not improve the prediction performance. It was observed that simple predictors with only multispectral information like E + H and E + I, provided similar accuracies to predictors with both types of information (genotypic + multispectral information) and interaction terms. However, we observed more stable and consistent predictions in predictors that integrated both sources of information. The predictions using H or I did not produce relevant differences since, in some cases, using I information provided slightly better results than using H or vice versa. In the models, we observed both models were effective in this prediction scenario, and with this data, the predictions using GBLUP were slightly better (Supplementary Table SA4).

In Figure 4 and Supplementary Table SA5, we note that in the GBLUP model, the best and worst predictions were observed in Prescott and Farmington, respectively, while under the PLS model, the environment’s performance showed no difference. However, in both models, the predictions were not as high as those observed in the previous years. Also, the worst predictions were observed with the predictor E + g (see Figure 4; Supplementary Table SA5) and the best joining of the two types of information (genomic + multispectral information under its H or I versions). Again, we observed that adding interaction terms gE, gH, and gI, in the predictors did not improve the prediction performance regarding the additive integration of both types of information. Also, we observed that simple predictors with only multispectral information, like E + H and E + I, provided more competitive accuracies than predictors with both types of information (genotypic + multispectral information) and interaction terms. Yet we did not observe using the predictor E + H to give the best predictions. More stable and consistent predictions were observed in predictors that integrated both sources of information. Also, regarding the predictions using H or I, we observed non-relevant differences between them since sometimes using I information provides slightly better results than using H or vice versa. We observed that both models faced difficulties in predicting some environments but that the GBLUP model outperformed the PLS model (Supplementary Table SA5).

In this section, the results under LOEO for each year are reported. For each model (GBLUP and PLS), five predictors were evaluated to see how much each input contributed to the improvement of prediction accuracy in a complete environment. For each year, the results are reported for each environment and across environments (Global).

Under the LOEO cross-validation, in both types of models (GBLUP and PLS), we observed the best predictions under Kincaid (Cor>0.75) and the worst under the Lind and Pullman environment with Cor values less than 0.5 (Figure 5; Supplementary Table SA6). The worst predictions were observed only when the genomic information was used (g) and the best when both types of information were integrated (Figure 5; Supplementary Table SA6). It was observed that the predictors with only multispectral information (H or I) provided similar accuracies to predictors with both types of information. Of note, in most environments, the best predictions were observed with the predictors with only multispectral information (H and I) and within those observations, the best predictions were observed using only the H information. The observed models displayed respectable predictions, albeit with considerably lower accuracy than those observed under the 7FCV strategy.

Under both types of models (GBLUP and PLS), again under the LOEO cross-validation, we observed the best predictions found under the Farmington, Harrington, Ritzville and Walla Walla (with Cor>0.75 most times) and the worst under the Kincaid and Lind environments, with Cor values less than 0.3 (Figure 6; Supplementary Table SA7). In addition, the worst predictions were observed using only the genomic information (g) and the best when only the multispectral information was used, with better performance using H in place of I, but with no large difference between the two sources of multispectral information. However, very competitive predictions were observed when both types of information were used. We noted good predictions with both models (Cor>0.75) for some environments but modest for others (Cor<0.3).

Under the LOEO cross-validation and both models in this data set, for all environments, it was possible to obtain strong predictions (Cor>0.75 in most cases), with the exception of Pullman, where the predictions in terms of Pearson´s correlation were around 0. Using only the genomic information (g) provided the worst predictions and the best predictions were obtained when only the multispectral information was used (Figure 7; Supplementary Table SA8), with better performance using I than H, but with similar performance between the two sources of multispectral information. However, joining both types of information provided competitive predictions. We observed that both models had strong predictions (Cor>0.75) for many environments, with only one poor result (Cor ≈ 0.0).

Under both models and the LOEO cross-validation for all environments, it was not possible to obtain good predictions since, in many environments and some predictors, the predictions in terms of Pearson´s correlation were less than 0.2. But even in this year with lower predictions, the worst performance was obtained using only the genomic information (g) and the best when only the multispectral information was used (Figure 8; Supplementary Table SA9), with better performance using I than H, but with similar performance between the two types of multispectral information. However, joining both types of information provided very competitive predictions. With both models, it was not possible to obtain good predictions (with Cor<0.5) with some predictions lower than 0.2.