The multi-site experiments included 169 potato breeding clones and cultivars in Helgegården, and 256 breeding clones and cultivars in Mosslunda and nearby Umeå (Supplementary Table 1, see link https://hdl.handle.net/11529/10548617). The breeding clones are in at least the fourth generation (T₄) of selection by Svensk potatisförädling of the Swedish University of Agricultural Sciences (Ortiz et al., 2020), while the cultivars are a sample of those released and grown in Europe during the last 200 years. Helgegården and Mosslunda are rural sites near the city of Kristianstad (56°01′46′′N 14°09′24′′E, Skåne, southern Sweden), while Umeå (63°49′30′′N 20°15′50′′E) is a city in the north of Sweden. The potato cropping season lasts about 3.5–4 months in Skåne (end of May–early September), while only 90 days in Umeå (early June–end of August). The average daily temperatures during the potato growing season vary between 12 and 18°C in the southern sites, and by 12.5–16°C in Umeå, while the average monthly precipitation amounts to 42 to 64 mm in sites near Kristianstad, and 48 to 75 mm in Umeå. The daylengths range between 11.5 h (toward harvest) and 17.5 h (about mid-growing season) in Skåne, and between 14.5 (around harvest) and ca. 21 h (at the beginning of the cropping season) in Umeå.

An incomplete block design, with two replications of 10 plants each, was the field layout for the field trials in Helgegården (13 × 13 simple lattice), Mosslunda (16 × 16 simple lattice), and Umeå (16 × 16 simple lattice). Fungicides were only used in Helgegården to avoid pests such as late blight (caused by the oomycete Phytophthora infestans) throughout the growing season, thus, allowing to estimate tuber yield potential at this site. Crop husbandry was used for potato farming at each site.

Total tuber yield per plot (kg), tuber weight by size (<40 mm, 40–50 mm, 50–60 mm, >60 mm; kg), and tuber flesh starch (measured as percentage based on specific gravity after harvest) were evaluated across all sites. Host plant resistance to late blight was evaluated using the area under the disease progress curve (AUDPC, Fry, 1978) in Mosslunda, while reducing sugars in the tuber flesh after harvest was determined using potato glucose strip tests (Mann et al., 1991) in Umeå.

Leaf samples --using 4 punches for each of the 256 breeding clones and cultivars included in the experiments-- were sent to Diversity Array Technology Pty Ltd (ACT, Australia) through AgriTech---Intertek ScanBi Diagnostics (Alnarp, Sweden) for further targeted genotyping following the genotype-by-sequencing approach.¹ The 2,000 single nucleotide polymorphisms (SNPs) used for genotyping were mostly derived from SolCAP SNPs based on chromosome positions and MAF >1 in germplasm from the Centro Internacional de la Papa (CIP, Lima, Perú) and the United States. According to Selga et al. (2021a), such a number of SNP already suffices for GEBVs without losing information. Although there were very few missing genotyping data (0.1%), one breeding clone (97) and two cultivars (‘Leyla’ and ‘Red Lady’) were not included further in the analysis because they were lacking enough SNP data.

We first briefly described the three different cases for codifying the molecular X matrix proposed by Slater et al. (2016) to be used in the genomic-enabled prediction models. Then, we defined the Bayesian linear single environment model and the multi-environment model, including the G × E using the GB and GK kernel methods.

Based on Slater et al. (2016), there are three cases for codifying the X matrix and thus, the type of genomic relationship matrices (Table 1). According to the authors, “there are at least two possible assumptions regarding the effect of marker allele dosage on phenotype for genomic selection. One assumption would be a pseudodiploid model, where all heterozygous genotypes have an equal effect on the genotype, and that the effects of the heterozygotes is at the midpoint of the two homozygotes.”

In this case, the marker matrix X is constructed as indicated in Table 1 (Slater et al., 2016) for the pseudo-diploid model with a column for each of the M markers with 0 for AAAA, 2 for BBBB, and 1 for any other form. Hence, the linear relationship between lines jandk for the GBLUP method (GB), can be constructed as:

The diagonal of matrix GB can be constructed as in Slater et al. (2016):

where M is the number of markers and p_i for the i^th marker is computed as:

where N is the total number of individuals and n_bbbb,n_abbb,n_aabb,n_aaab is the number of individuals of genotypes BBBB, ABBB, AABB, and AAAB, respectively.

To model a more complex relationship between the lines, the Gaussian kernel, defined as GK=exp⁡(-h⁢di⁢i′2/q), could be used where h is the bandwidth parameter that controls the rate of decay of the covariance between genotypes, and q is the median of the square of the Euclidean distance d_ii′=∑_k(x_ik−x_i′k)², which is a measure of the genetic distance between individuals based on molecular markers. The bandwidth parameter h was estimated based on the empirical Bayes proposed by Pérez-Elizalde et al. (2015).

Following Slater et al. (2016), Additive tetrasomic is adapted for estimating the additive marker effect by accounting for the tetraploid allele dosage. In this case, X has the dimensions N × M, but the new coding is now 0, 1, 2, 3, and 4, for AAAA, AAAB, AABB, ABBB, and BBBB, respectively.

In this study, matrix X is standardized by column (mean equals to zero and variance equals to 1) as such and according to Lopez-Cruz et al. (2015):

As in the previous case, the Gaussian kernel can be constructed as GK=e⁢x⁢p⁢(-h⁢di⁢i′2/q).

An alternative option for coding matrix X according to Slater et al. (2016) is considering additive and non-additive effects in a full tetrasomic, assuming each genotype has its own effect. In this case, there are five possible effects per SNP marker (Table 1). Then the genomic relationship between individuals j,k is computed as:

where M is the number of markers × 5. To compute the diagonal of this matrix, we can use:

where p_i is the frequency of each genotype, i.e., the frequency in each column. The Gaussian kernel can be calculated as in the previous cases.

Here we consider the single-environment and multi-environment models, each combined with two methods, linear kernel GBLUP (GB), and non-linear Gaussian (GK). In addition, each of these combinations of model/method were tested with the three single-kernel methods derived from the X with marker dosage A, B, and C and the two-kernel methods combining X with marker dosage B and C.

The two-kernel method attempts to exploit the additive effects of the genomic matrix (B) and the non-additive case (C), as explained below for a single environment and multi-environments under the two kernel methods (GB and G). Thus, each of the single-environment and multiple-environment models under GB and GK had three different single-kernel methods (for A, B, and C), and one two-kernel methods (B and C).

For the single site (models 1 and 2, Eqs. 1 and 2), we extracted 30 random samples to form groups; 70% for the training set (TRN) and 30% to be predicted (testing, TST set). For the multi-environment (models 3 and 4, Eqs. 3 and 4), we extracted four random folds each with 10 samples using a random cross-validation called CV2 that consists of predicting one line in one site, knowing the value of that line in at least one of the rest of the sites. We used Monte Carlo Markov Chain (MCMC), using the BGGE software (Granato et al., 2018), to fit the four models and to predict the individuals in the TST sets. For each of the samples, we computed the genome-based predictions and correlated them with the observed values. We reported the mean of the correlation between the predicted and the observed values and its standard deviations.

Markov Chain Monte Carlo (MCMC) diagnostics are tools that are used to investigate whether the quality of a sample generated with an MCMC algorithm is sufficient for providing an accurate approximation of the target distribution. The MCMC has diagnostic tools for testing (1) whether a large portion of the MCMC sample has been drawn from distributions that are significantly different from the target distribution or for observing (2) whether the size of the generated sample is too small. In this study in order to minimize the random error from the MCMC, we performed 30,000 iterations, with a burn-in of the first 5,000 and a thinning of 2. Supplementary Figure 1 displayed the MCMC results from the model 1 analyses of one trait; the genetic variance components are correctly mixed, achieving a correct convergence and thus, generating the posterior probability distribution.

The marker data as well the phenotype data for each of the three environments are stored at the link https://hdl.handle.net/11529/10548617.

Tables 2–4 and Figures 1–3 give the 2020 results of Helgegården, Mosslunda, and Umeå, respectively, for several traits, of the average correlation (between observed and predictive values) and their standard deviation for the single trial analysis for each of methods GBLUP (GB) and GK, considering (model 1) constructing matrix X as pseudo-diploid (A), additive tetrasomic polyploid (B), full tetraploid (C), and the combination of B-C (model 2).

The results (Tables 2–4) show unclear trends in the genome-based prediction accuracy comparing the different structures of the marker matrices and the methods (GBLUP vs. Gaussian kernel) for the different traits. At Helgegården, the two-kernel combinations (model 2) show an increase in prediction accuracy for most of the traits under the GK method as compared with those obtained under A, B, and C marker structures (Figure 1). However, at Mosslunda, model 2 increased the prediction accuracy of traits under GB, as well as GK (Figure 2). For Umeå, an unclear trend of marker forms and methods were found (Figure 3); however, model 2 and model 1 under B are always the best for starch % under GB and GK.

The trait with the highest prediction accuracy for all sites under A, B, C (model 1), model 2, and for the GB and GK methods was the highly heritable tuber starch percentage. Another trait with relatively high prediction accuracy was the total tuber weight. Concerning the single-kernel (model 1) vs. the two-kernel method (model 2), evidences show an increase in prediction accuracy of the combination of two kernels (model 2) over model 1. Non-linear Gaussian kernel (GK) does not show any clear advantage over the linear kernel GBLUP (GB).

Tables 5–7 and Figures 1–3 give the prediction results of Helgegården, Mosslunda, and Umeå, respectively, for several traits and their standard deviation for the multi-environment analyses for each of the methods (GB and GK), considering the single-kernel (model 3) for pseudo-diploid (A), additive tetrasomic polyploid (B), full tetraploid (C), and the two-kernel combination of B-C (model 4).

In general, results including G × E interaction in the multi-environment analyses exploit the information on the relationship between the location-year combinations and had prediction accuracy estimates higher than those obtained from the single-environment analyses (results from models 1 and 2 vs. models 3 and 4) (Figures 1–3). For the three sites, two-kernel model 4 with GB seems to be the best method in combination with the additive tetrasomic polyploidy structure B for predicting most of the tuber traits. Most of the traits gave relatively high prediction accuracy under this combination of marker structure (A, B, C, and B-C) and the combined methods of GB and GK, including the multi-environment with the G × E model.

In general, genome-based prediction accuracy obtained in this potato study, using different marker similarity matrices accounting for additive and non-additive marker relationship under single-environment and multi-environment models, show prediction accuracy patterns similar to those found in other studies using other species with different levels of ploidy. The process of borrowing information from multi-environment trial analyses modeling G × E provides a very useful increase of genomic-enabled prediction accuracy over the evidence obtained from the single-environment analyses. This increase in prediction accuracy of G × E models has been clearly and extensively documented, among others, in Burgueño et al. (2012), Jarquín et al. (2014), Crossa et al. (2017, 2019), and Cuevas et al. (2017, 2018, 2019) where the genomic similarity between cultivars is increased when modeling the phenomenon of G × E. That is, the appropriate statistical modeling of G × E allows borrowing information from correlated environments to the predictions of unobserved phenotypes in environments. For all the agronomy traits included in this study, the important increase in prediction accuracy when genomic prediction models include G × E models is clear. Supplementary Table 2 shows the relatively high and positive phenotypic correlations between the three sites and the six traits included in this study that explain part of the increase in genomic based prediction achieved by models including G × E as compared with the single trait models.

Two type of random cross-validation are usually employed for comparing different models and methods. Burgueño et al. (2012) and Lopez-Cruz et al. (2015) distinguished a random cross-validation 1 (CV1) when predicting lines that were never evaluated in any environment, and random cross-validation 2 (CV2) that consists of some lines tested in only some environments but not in others. Extensive results from Burgueño et al. (2012), Jarquín et al. (2014), and Lopez-Cruz et al. (2015) demonstrated that a genome-based prediction accuracy obtained from CV1 are similar to those obtained when using a single environment (site) genome-based prediction model. Lopez-Cruz et al. (2015) mentioned that “This feature of the M × E model can be exploited in prediction problems such as CV2; however, such borrowing of information within line is not possible in CV1 and, consequently, the M × E model performs similarly to the stratified analysis for prediction of performance of lines that have no phenotypic records.”

To investigate the results of CV1 in the this study we have assessed the genomic-enabled prediction accuracy of the trait tuber weight, using the full tetrapoid (C) GBLUP kernel with model 1 (single-site), and compared it with the G × E model 3 (multi-site). The results were similar for the three sites. For site Helgegarden, model 1 gave an average genomic prediction accuracy of 0.418, whereas model 3 gave an average prediction correlation of 0.433. Similar for site Mosslunda, model 1 gave an average prediction accuracy of 0.563 vs. 0.568 as mean prediction accuracy for model 3, whereas for site Umea, the mean genomic-enabled prediction for model 1 was 0.455 and for G × E model 3 model 3 was 0.450.

This is explained by the exchange of information (borrowing of information) that is achieved in the main effects component, that is, Z_2g, where g had the effect of each line that are predicted throughout the environments. For G × E model 3 (multi-site), borrowing (exchanging) information between lines occurs only if the lines have genetic and environmental similarity. These results showing a similar genome-based prediction accuracy between model 1 and CV1 vs. G × E model 3 (multi-site) are similar (with small differences between models) to those shown by Lopez-Cruz et al. (2015; Tables 5–7).

The use of Gaussian kernels has been extensively documented in genome-based studies as a non-linear kernel that increases prediction accuracy over the linear kernel given by the linear additive GBLUP. The non-linear kernels included in multi-environment G × E models have been shown to increase the genome-based prediction accuracy by around 5–10% in several studies (Cuevas et al., 2016, 2017, 2018, 2019; Crossa et al., 2019).

As noted in our study, including the GK in the multi-environment G × E model did not overcome the genome-based prediction accuracy over the GB method. One of the reasons could be that the some of the marker structures employed, like the full tetraploid, already account for the additive and non-additive structure of the markers; thus, no extra benefit is obtained by including models that will exploit these cryptic and small epistatic inter-locus interactions between markers. Another possibility is that the GK is not able to capture the residual epistatic interactions that may exist in the tetrasomic polyploid potato even after the use of a similarity marker structure that considers the linear additive kernel and the non-linear kernel of the full tetraploid. The two kernels (B and C), with the multi-environment G × E model using the linear GBLUP (GB) kernel, seem to capture most of the potential marker epistatic interactions without the need to add the non-linear GK kernel. More research is required in this area.

Table 8 provides an up-to-date summary information on all available journal articles regarding prediction accuracy estimates of GEBV for selection in potato. This table only includes traits that were evaluated in our research: i.e., tuber weight (total and by size), host plant resistance to late blight, tuber starch percentage, and crisp quality. The prediction accuracy estimates for tuber weight and tuber starch percentage are equal or above those available in the literature, or within the known ranges for both the host plant resistance to late blight (with a bias toward high correlations) and the crisp quality as measured by reducing sugars. These are encouraging results because they show that the multi-trait, multi-environment modeling of the GEBV increased the prediction accuracy estimates, which may also vary according to training population size and type, trial data quality, and method or model use.