Gaussian Process Modeling with

Genotype × Environment Kernels for Wheat

Performance Prediction
Lea Friedli∗, Tim Steinert†, Nathalie Wuyts‡, Fabian Guignard§,
arXiv:2508.17730v1 [[Link]] 25 Aug 2025

Lilia Levy Häner‡ , Didier Pellet, Juan M. Herrera‡ and David Ginsbourger†

Abstract

Optimizing wheat variety selection for high performance in different environmen-

tal conditions is critical for reliable food production and stable incomes for growers.
We employ a statistical machine learning framework utilizing Gaussian Process (GP)
models to capture the effects of genetic and environmental factors on wheat yield
and protein content. In doing so, selecting suitable covariance kernels to account
for the distinct characteristics of the information is essential. The GP approach is
closely related to linear mixed-effect models for genotype × environment predictions,
where random additive and interaction effects are modeled with covariance structures.
However, while commonly used linear mixed effect models in plant breeding rely on
Euclidean-based kernels, we also test kernels specifically designed for strings and time
series. The resulting GP models are capable of competitively predicting outcomes for
(1) new environmental conditions, and (2) new varieties, even in scenarios with little
to no previous data for the new conditions or variety. While we focus on a wheat test
case using a novel dataset collected in Switzerland, the GP approach presented here
can be applied and extended to a wide range of agricultural applications and beyond,
paving the way for improved decision-making and data acquisition strategies.

Keywords covariance kernels, genotype-environment predictions, leakage, multi-environment

variety trials, uncertainty quantification
∗
Engineering Risk Analysis Group, Technical University of Munich, Germany
†
Institute of Mathematical Statistics and Actuarial Science, University of Bern, Switzerland
‡
Agroscope, Plant-Production Systems, Switzerland
§
METAS, Federal Institute of Metrology, Switzerland

1
1 Introduction
Wheat is a staple crop essential for global food security, and its growth is influenced by an
interplay of genetic, environmental, and management factors. Understanding how different
wheat varieties respond to varying environmental conditions is crucial, especially in the face
of climate change. When the goal is to identify genotypes optimal for specific environments,
the patterns of genotype × environment (G×E; Kang & Gorman 1989) interactions must
be considered. Since not all combinations can be tested in the field, statistical and machine
learning models serve as efficient surrogates (e.g., Fernandes et al. 2024). In this paper,
we consider Gaussian Process (GP) modeling for G×E prediction, exploring various kernels
defined on the G×E product space.

We suppose that we observe some output values (potentially noisy) at certain input points,
and our goal is to predict the output at new input points. In the context of G×E predic-
tions for plant selection, the output usually represents a performance trait, such as yield
or grain protein content, and the inputs correspond to combinations of environments and
varieties/genotypes (treated from here on as synonyms). We consider a setting in which we
have additional information about the input combinations, as defined by the covariate tuple
x=(xG , xE ). For the environmental variables xE ∈ XE , we use meteorological variables that
describe temperature, precipitation, and solar radiation, which are measured over six time
periods within the wheat crop cycle (Costa-Neto et al. 2021). For the genetic information
vectors xG ∈ XG , we consider single nucleotide polymorphisms (SNPs; e.g., Rafalski 2002).

One modeling approach which has been widely used in G×E prediction is linear mixed-effect
models (LMMs; e.g., Piepho 1997; Crossa et al. 2006; Herrera et al. 2018; Buntaran et al.
2021). A baseline LMM for G×E prediction models the response as the sum of an overall
mean plus fixed and random effects due to the environment and/or the variety, plus an error
term. Thereby, the random effects and errors are assumed to be Gaussian and independent,
so no information is shared between different varieties or environments (Jarquı́n et al. 2014).
However, a key advantage of a LMM is its ability to account for covariance structures, en-
abling strength to be borrowed across groups (Crossa et al. 2006; Burgueño et al. 2007). The
numerous studies employing LMMs differ in their specification of fixed and random effects,
as well as in how they handle interactions. Building on the overview provided by Crossa
et al. (2022), we summarize the main modeling approaches used in the literature.

The basic single-environment genomic model analyzes genotypes within one environment,

2
fitting a separate model for each. When genetic markers are available for the genotypes,
the genetic effect is typically modeled using a parametric linear regression on the molecular
markers. Assuming that the vector of marker effects follows an independent and identically
distributed (i.i.d.) normal distribution leads to the ridge regression best linear unbiased
prediction (rrBLUP) model. Alternatively, the genetic effect can be modeled using a cen-
tered multivariate Gaussian distribution, where the covariance matrix is modeled with a
genomic relationship matrix derived from a kernel applied to the marker data. A specific
choice of linear kernel yields the genomic best linear unbiased prediction (GBLUP) model
(VanRaden 2007). Other kernels, such as the Gaussian kernel, have been proposed to better
capture non-linear marker effects and have shown improved performance in several studies
(Cuevas et al. 2016, 2017; Bandeira e Sousa et al. 2017). More recently, kernels informed by
deep learning techniques have also been explored (Cuevas et al. 2019; Costa-Neto et al. 2021).

Multi-environment trials allow to borrow information across environments. One initial ap-
proach allows marker effects or genetic values to vary across environments (Schulz-Streeck
et al. 2013; Burgueño et al. 2012). However, such models typically enable prediction only
within already observed environments. To enable prediction in new environments, environ-
mental covariates can be incorporated to model dependencies among the random environ-
mental effects, analogous to how genetic relationships are modeled using molecular markers.
Models which specify genotype-specific variations due to key environmental factors were af-
terwards called reaction norm models (Crossa et al. 2022). Capturing interactions between
molecular markers and environmental covariates is challenging due to structural differences
and potential high dimensionality. To address this, Jarquı́n et al. (2014) propose a modeling
approach assuming that the interaction effect follows a normal distribution, with the covari-
ance structure given by the element-wise (Hadamard) product of the two separate covariance
matrices. Cuevas et al. (2016, 2017) and Bandeira e Sousa et al. (2017) demonstrated that in-
corporating a Gaussian kernel into multi-environment genomic models significantly enhances
predictive performance compared to models employing a standard linear kernel. Building on
this idea, Costa-Neto et al. (2021) extended the approach of Jarquı́n et al. (2014) by using
both a Gaussian and a Deep kernel, within a Bayesian framework.

LMMs are typically fitted using (restricted) maximum likelihood and least square estima-
tion. However, LMM methodology can also be reformulated in a Bayesian framework by
treating all unknown quantities as random variables with assigned prior distributions (e.g.,
Meuwissen et al. 2001). Bayesian inference then proceeds via the posterior distribution,
allowing for joint estimation and full uncertainty quantification. Employing a full Bayesian

3
framework, Costa-Neto et al. (2021) use a Gibbs sampler to generate the posterior distribu-
tions of the hyperparameters, effects and predictions. Recently, Liu et al. (2025) employed a
similar kernel-based model, where predictions are obtained through analytical formulas for
Gaussian processes, following hyperparameter optimization based on Gibbs. In the work of
both Costa-Neto et al. (2021) and Liu et al. (2025), SNP markers are numerically encoded,
such that for both environmental and genetic information, Euclidean-based kernels are used.
However, genomic data (e.g., SNPs) and environmental variables (e.g., time series from me-
teorological data) have distinct structures—strings and temporal signals, respectively. This
motivates the use of more specialized kernels, here introduced within the framework of GP
modeling.

GP models originate from the geostatistical interpolation method known as kriging (Krige
1951), but have broadened their application far beyond geostatistics (Rasmussen & Williams
2006). GP modeling, aiming (in it its most common form) to learn a function f : X → R,
can be applied to a wide range of input spaces X , handling the diverse inputs through tai-
lored covariance kernels. Apart from its versatility, GP modeling is popular for enabling
uncertainty quantification, handling small training datasets, allowing the incorporation of
prior knowledge through the covariance kernel, and providing interpretable results. In recent
years, GPs have been increasingly applied across a variety of domains, including chemistry
(e.g., Griffiths et al. 2023), engineering (e.g., Su et al. 2017) and environmental modeling
(e.g., Ray 2021). In the considered GP approaches for G × E modeling, there are three
major degrees of freedom illustrated in Figure 1: the choice of the kernel for i) the genetic
effects, ii) the environmental effects, and iii) the strategy used to construct a joint kernel
over the product space of genotypes and environments. Note: Here, “GP” refers to Gaussian
process, not genomic-enabled prediction, which is also commonly abbreviated as “GP” in
agricultural literature.

G × E kernels and associated GP models offer promising potential for decision-making and
data acquisition strategies. Once a model is fitted, it provides probabilistic predictions at
any candidate input, enabling the search for inputs that yield optimal responses. As pre-
dictions are probabilistic, optimality would have to be defined in terms of the predictive
distribution, be it in from a multi-objective perspective (e.g., in the spirit of Markowitz
portfolio optimization) or via a single criterion such as a predictive quantile at a specified
level (e.g., Picheny et al. 2013). GP models have been used for global optimization for several
decades and the resulting field now referred to as Bayesian Optimization (BO) has become
very active from engineering to machine learning and beyond (Močkus 1974; Jones et al.

4
1998; Garnett 2023). GP-based modeling often separates controllable from non-controllable
variables, aiming to optimize performance over the controllable ones under average or worst-
case conditions (Williams et al. 2000; Janusevskis & Le Riche 2013; Ginsbourger et al. 2014).

In optimizing agronomic traits, one immediate application of GPs is using predictive dis-
tributions to compare and select varieties for specific environments. Identifying the best
genotype among a candidate population has also been addressed in the realm of LMMs in
Tanaka & Iwata (2018) and Tsai et al. (2021). However, environmental conditions have
been ignored, and since the exact conditions during the target period may be unknown at
the time of variety selection, GP methods are particularly promising for handling the asso-
ciated uncertainties. Apart from genotype selection, GP models may be leveraged to design
novel agronomical experiments. Experimental design in the G × E context has also been
targeted using LMMs, for instance to optimize the allocation of trials to sub-regions in Prus
& Piepho (2024). Here, the selection strategy aims to optimally predict the pairwise linear
contrasts of the genotypes among some pre-defined sub-regions. In comparison, GPs provide
a very flexible framework for any input variable space. Furthermore, GP models allow to
work out acquisition functions dedicated to targeted experimental design, e.g. focusing on
specific ranges of the response (Chevalier et al. 2014).

GP regression is closely related to LMM approaches; however, GPs typically rely on max-
imum likelihood estimation for hyperparameters and use analytical expressions for predic-
tions, without requiring sampling-based inference (see the analytical Gaussian predictions in
Fig. 1). Given the potential of GPs in the context of decision-making and data acquisition
strategies, our aim is to bridge the gap between the methods and pave the way for promising
future applications. Furthermore, while the state-of-the-art LMMs for G×E prediction rely
on Euclidean-based kernels for both environmental and genetic inputs, we also test kernels
specifically designed for strings and time series. While Costa-Neto et al. (2021) investigate a
test case involving maize hybrids, this manuscript analyzes a novel dataset targeting wheat
yield and grain protein content, collected between 1990 and 2023 in multi-environment tri-
als across the Swiss wheat production zone (Levy Häner et al. 2025). The manuscript is
organized as follows: Section 2 reviews GP modeling, G×E kernels, and links to LMMs.
Section 3 describes the dataset, implementation, metrics, and baselines. Finally, Section 4
presents the wheat test case, and Section 5 concludes with discussion and summary.

5
Figure 1: Illustration of GP modeling for G×E prediction: Given training data from specific
G×E configurations, the goal is to predict outcomes (here the wheat grain protein content)
for a new scenario. Kernels on G and E are used to quantify similarity between the training
data and the target scenario, effectively weighting the training points (indicated by trans-
parency of the colorful dots). The kernels are then combined to construct a joint kernel over
the product space, and predictions are made based on the weighted training data according
to this combined kernel. In GP modeling, the predictive distributions are Gaussian, as illus-
trated by the bell-shaped probability density functions depicted in prediction.

6
2 Gaussian Process Regression with combinations of
Genotype and Environment Kernels
2.1 Gaussian Process basics
Gaussian Processes (GPs) offer a flexible probabilistic approach to modeling unknown func-
tions. A GP, denoted by ξ = (ξ(x))x∈X , is a collection of real-valued random variables
defined over the same probability space. The defining property of a GP is that for any finite
set of inputs {x1 , . . . , xn } ⊂ X with n ≥ 1, the outputs ξ(x1 ), . . . , ξ(xn ) form a multivariate
Gaussian vector. The GP is fully specified by a mean function m : X → R, m(x) := E[ξ(x)],
and a covariance function (kernel) k : X × X → R, k(x, x′ ) := Cov(ξ(x), ξ(x′ )).

In the considered GP modeling context, the target function f : X → R is treated as a sample

from a GP ξ, which, in a Bayesian context, can be seen as a prior model. We assume that we
observe noisy measurements of f at n input points x1 , . . . , xn , resulting in observations of the
form Zi = f (xi ) + εi , where the noise terms are i.i.d. with εi ∼ N (0, τ 2 ). By conditioning
the GP on the observed data An = {Z n = z n }, we obtain a posterior GP characterized
by a predictive mean function mn (x) and a posterior covariance function kn (x, x′ ). We
consider the case of ordinary kriging (assuming a constant but unknown mean m), for which
the formulas of the posterior mean and covariance can be found in Roustant et al. (2012).
Selecting an appropriate positive definite kernel k(x, x′ ), x, x′ ∈ X is far from trivial and
requires careful consideration, particularly when dealing with non-Euclidean input spaces.
In the following section, we discuss ways to design kernels for the G×E product spaces.

2.2 Kernel combination approaches on product spaces

Let us first assume that valid kernels kG and kE are given on XG and XE , respectively. We
now discuss ways to combine them into valid kernels on XG × XE . Let us denote here by xG
a generic vector of genetic variables in XG , by xE a generic vector of environmental variables
in XE , and by (xG , xE ) a generic (concatenated) vector in the product space XG × XE . We
are thus interested in combining kernels kG and kE respectively defined on XG and XE (that
is, with respective input spaces XG2 = XG × XG and XE2 = XE × XE ) into valid kernels kGE
on XG × XE (thus taking their inputs in (XG × XE )2 = (XG × XE ) × (XG × XE )).

A first important remark is that kG and kE can directly be used as kernels on the product
space, in the sense that kernels (informally) defined by kGE G
((xG , xE ), (x′G , x′E ) = kG (xG , x′G )
and kGEE
((xG , xE ), (x′G , x′E ) = kE (xE , x′E ) are both valid on the product space. Either

7
option would amount to ignore part of the variables, which is expected to be detrimental
in prediction but would not alter the non-negative definiteness of resulting kernel matrices,
thereby delivering valid kernels on XG × XE . The combination approaches considered in this
work essentially revolve around (blockwise) tensor sums and products, in the following sense.
First, the tensor sum,
+
kGE ((xG , xE ), (x′G , x′E ) = kG (xG , x′G ) + kE (xE , x′E ),

defines a valid kernel XG × XE , following the same mechanism as for usual additive kernels.
This operation is immediately extended to non-negatively weighted sums. Beyond this,
×
kGE ((xG , xE ), (x′G , x′E ) = kG (xG , x′G ) × kE (xE , x′E ),

is also known to provide valid kernels on XG × XE , in virtue of the stability of the cone
of symmetric non-negative matrices by matrix product. Using again summation and non-
negative linear combinations, we arrive at a class of kernels synthesizing and extending the
latter approaches, whereby α, β, γ ≥ 0:
∼
kGE ((xG , xE ), (x′G , x′E ) = αkG (xG , x′G ) + βkE (xE , x′E ) + γkG (xG , x′G ) × kE (xE , x′E ). (1)

Let us also remark that taking a different (kG , kE ) pair (e.g., with different hyperparameter
values) in the sum and product parts of Equation (1) still leads to valid kernels on XG × XE .

2.3 Linear mixed-effect models as Gaussian Processes Regression

Before discussing specific kernels, we want to connect the GP modeling framework to the
LMM setting. Using LMM notation accounting for G×E interactions and covariates, we can
write the observation vector Z n as (Costa-Neto et al. 2021; Liu et al. 2025),

Z n = 1n m + U G + U E + U GE + ε,

where m denotes the general mean and ε ∼ N (0, τ 2 I n ) a n-dimensional noise vector.
Furthermore, the random effects of genotype U G ∼ N (0, σG 2
Z GK GZ ⊤ G ) and environment
2 ⊤ 2 2
U E ∼ N (0, σE Z E K E Z E ) are considered, where σG and σE denote the variances. The
incidence matrices Z G and Z E are based on a one-hot encoding of the inputs and K G
and K E denote the relationship matrices between the genotypes and the environments,
respectively. If the relationship matrices are based on the covariance kernels kG and kE ap-
plied to the covariates, we can write Z E K E Z ⊤ ⊤
E = [kE ((xE )i , (xE )j )]1≤i,j≤n and Z G K G Z G =
[kG ((xG )i , (xG )j )]1≤i,j≤n . Finally, the interaction term U GE ∼ N (0, σGE 2
Z GK GZ ⊤ ⊤
G ◦Z E K E Z E )
2
is modeled using the variance σGE , where ◦ denotes the Hadamard product of two matrices.

8
By Gaussianity of sums of independent Gaussian-distributed random variables,
2
Z GK GZ ⊤ 2 ⊤ 2 ⊤ ⊤ 2


Z n ∼ N m, σG G + σE Z E K E Z E + σGE Z G K G Z G ◦ Z E K E Z E + τ In ,

∼
what coincides with the GP assumption under the combination approach kGE (Eq. 1).

To predict at a new input x, it is assumed that the response at x follows a joint Gaussian
distribution with Z n . If x is then conditioned on Z n , we obtain the same posterior Gaussian
as with GP modeling (up to some variations due to trend parameters, see Roustant et al.
2012). In the same way, a LMM approach without G×E interaction could be linked to GP
+
modeling by using the additive kernel combination kGE and a LMM with fixed effects could
be modeled for instance by using a Universal Kriging approach (Matheron 1969; Handcock &
Stein 1993). However, while LMMs typically rely on frequentist or full Bayesian approaches,
GP models employ empirical Bayes and analytical formulas for prediction.

2.4 Kernels
Commonly used families of covariance kernels for Euclidean input spaces include the isotropic
Matérn kernels, which also cover the isotropic exponential and Gaussian (or square-exponential)
kernels as special cases (Stein 2012). In this context, isotropy implies that the kernel values
for pairs of locations depend solely on the Euclidean distance between them. For multidi-
mensional inputs, anisotropic versions of these kernels are frequently employed. Here, we
consider the isotropic exponential and Gaussian kernel for x, x′ ∈ Rν ,
||x − x′ || ||x − x′ ||2
   
′ ′
kEXP (x, x ) = exp − , kGAU (x, x ) = exp − ,
θ θ2
p
with ||x − x′ || = (x − x′ )⊤ (x − x′ ) and the correlation length θ > 0. Additionally, we are
considering kernels on non-Euclidean spaces, which we discuss in the next sections.

2.4.1 Considered kernels on G

The genetic information xG ∈ XG considered is given by sequences of strings (SNPs).

Gaussian-GBLUP Kernel (kG:GAU−GBLUP ) In this approach, the SNP data are numeri-
cally encoded based on the major allele frequency. A Gaussian kernel is then applied to the
resulting Euclidean distance matrix, following the methodology described by Cuevas et al.
(2016) and Costa-Neto et al. (2021).

Exponential-Hamming Kernel (kG:EXP−HAM ) The Hamming distance for string se-

quences xG , x′G ∈ XG is given by dH (xG , x′G ) = ν1G νi=1 1(xG [i] ̸= x′G [i]) and can be used
PG

9
in the exponential kernel kEXP instead of the Euclidean distance. The Hamming distance,
which counts the number of differing entries between two vectors, can be viewed as the L∞
distance. Hutter et al. (2014) introduced a generalized (weighted) version of this kernel and
demonstrated its positive definiteness. While this result applies to the exponential-Hamming
kernel, it does not hold for the Gaussian-Hamming kernel (following a foundational theorem
by Schoenberg, see Berg et al. 1984).

Spectrum Kernel (kG:SPE ) The k-spectrum kernel, introduced by Leslie et al. (2001),
captures similarity by comparing substrings of length k. The spectrum kernel relies on
Mercer’s theorem (Mercer 1909), defining the kernel as the inner product in a feature space.
Therefore, it considers all possible subsequences a of length k from the alphabet A and
k
defines a feature map from XG to R|A| as Ψk (x) = (ψa (x))a∈Ak , where ψa (x) represents
the frequency of subsequence a occurring in x. The k-spectrum kernel is then given by
kG:SPE (x, x′ ) = ⟨Ψk (x), Ψk (x′ )⟩.

2.4.2 Considered kernels on E

In the considered wheat test case, xE ∈ XE contains meteorological variables that describe
temperature, precipitation, and solar radiation over six time periods.

Exponential-Euclidean (kE:EXP−EUCL ) and Gaussian-Euclidean Kernel (kE:GAU−EUCL )

For these two kernels we apply the exponential and Gaussian kernel to the Euclidean dis-
tances of the environmental variables.

Global Alignment Kernel (kE:GAK ) Global Alignment Kernels (GAK; Cuturi et al. 2007)
compute similarities between time series by building on Dynamic Time Warping (DTW;
Sakoe & Chiba 1970). Since DTW does not satisfy the triangle inequality, it is strictly not
a metric and cannot be used directly as a kernel. As a way around this, Cuturi et al. (2007)
propose to rely on a soft-minimum, and get the following positive definite kernel,
|ℓ|
X X
′ − dx,x′ (ℓ)
φ xℓ1 (i) , x′ℓ2 (i) ,


kE:GA (x, x ) = e , where dx,x′ (ℓ) =
ℓ∈L(x,x′ ) i=1

using the set of all alignments L(x, x′ ) and with the function φ(·) typically chosen as the
squared Euclidean distance.

10
3 Dataset and Implementation
3.1 Wheat dataset
The wheat data (Levy Häner et al. 2025) were collected in multi-environment trials across
the Swiss wheat production zone (from Western to Eastern Switzerland, between 390 and
640 MASL), conducted by Agroscope (the Swiss Federal Center for agricultural research)
and partners as part of its official Swiss winter wheat variety trials for cultivation and use
(Strebel et al. 2025). The data includes the performance traits yield (15% moisture content,
average of three replicates, dt per ha) and grain protein content (measured on a mixture of
three replicates, %) per location and year, for a total of 98 varieties. These were grown at 16
locations over the period 1990 to 2023. The management conditions of the trials were low
input, i.e. a standard nitrogen fertilization rate and no chemical inputs in the form of growth
regulators, insecticides or fungicides. The term ‘environment’ is referring to a combination
of trial location and harvest year; we have 263 environments available.

Genetic information consists of 12106 SNPs (Single Nucleotide Polymorphisms) marker data.
The data were characterized by a SNP call frequency above 80%, a frequency of homozygous
calls above 90%, and a sample call rate above 90%. The SNPs are in IUPAC notation with
‘G’, ‘A’, ‘T’ and ‘C’ representing guanine, adenosine, thymine and cytosine bases, and ‘K’,
‘M’, ‘R’ and ‘Y’ representing combinations of these bases as specified in Cornish-Bowden
(1985). For the Gaussian-GBLUP kernel, bi-allelic marker data is employed with homoge-
neous ‘G = GG’, ‘A =AA’, ‘T=TT’, ‘C=CC’ and heterogeneous ‘K=GT’, ‘M=AC’, ‘R=AG’,
‘Y=CT’. Meteorological data are averaged over six periods covering the Swiss winter wheat
growing season: winter (October-February, sowing and vernalization), March (tillering),
April (stem elongation), May (heading-flowering), June (flowering and seed filling) and July
(harvest). Daily minimum, maximum, and mean temperatures, total and average daily
precipitation, and monthly dry days were extracted for each location from the MeteoSwiss
Spatial Climate Analysis products (1 km grid; MeteoSwiss 2024b). Daily solar radiation was
obtained from the nearest automated meteorological station (MeteoSwiss 2024a).

3.2 Model evaluation

For evaluation, data are split into training x1 , . . . , xn and test set xn+1 , . . . , xm and we
perform cross-validation on 80 % of the data using 30 splits. We evaluate predictive per-
formance in two scenarios: new environment and a new variety. Group-based splitting, by
environment or variety, prevents information from leaking into the training set. We also con-

11
sider controlled leakage, allowing one observation in the training set to mimic minimal prior
knowledge. For a more detailed discussion of data leakage and its implications in group-based
cross-validation, we refer to Guignard et al. (2024). In agricultural studies, leave-one-out ap-
proaches exclude and predict individual environments or varieties (e.g., Tadese et al. 2024).
Our approach extends this by predicting multiple targets per split while keeping the number
of training data roughly consistent.

To assess the predictive accuracy of the GP, we first use the Mean Squared Error (MSE).
This metric assesses the mean mn (x) as a predictor, ignoring full probabilistic performance.
To capture the probabilistic aspect, the Gaussian predictive distribution must be com-
pared to a single observed value. This can be achieved through the use of a scoring rule
(Gneiting & Raftery 2007). Scoring rules assess performance by assigning lower scores to
“better” probabilistic predictions, ideally capturing both calibration (agreement with ob-
served outcomes) and sharpness (concentration of the predictive distribution). We consider
the Continuous Ranked Probability Score (CRPS; Sanders 1963, Murphy 1973), defined as
CRPS(F, y) = −∞ [F (u) − 1{y ≤ u}]2 du, where F (·) is the cumulative distribution function
R∞

of the predictive distribution (Gneiting & Raftery 2007). For a deterministic point predictor,
where the predictive distribution reduces to a Dirac delta, the CRPS simplifies to the abso-
lute difference between the predicted and the true value. Finally, we compute the median
CRPS for the test set. As a second scoring rule, we consider the log-score (logS; Good
1952), logS(f, y) = − log p(y), where p(·) is the probability density function of the predictive
distribution. Again, we consider the median score over the test set.

3.3 Implementation
A key implementation step is estimating the hyperparameters. The covariance matrix of the
training data in the noisy case with noise variance τ 2 is given by K obs = σK2
K + τ 2 I n , where
2
K = [k(xi , xj )]1≤i,j≤n . We summarize ν = σK + τ 2 and denote the proportion explained
2
by the GP as ς = σK / (τ 2 + σK
2
) . This gives K obs = νK ς , with K ς = ςK + (1 − ς)I n .
The optimal ν, as a function of (θE , θG , ς), can be estimated by likelihood maximization
(Roustant et al. 2012). The advantage of this approach lies in including the observation
noise without increasing the number of hyperparameters, while restricting the search space
of ς to [0, 1]. Recall the kernel sum and product combinations on the product space (Section
2.2). To avoid over-parameterization, we constrain (α, β, γ) to sum to one. With one hy-
perparameter per kernel, denoted as θE for the environmental and θG for the genetic kernel,

12
 ∼ + G E
this gives five parameters for kGE , four for kGE (γ = 0), and three for kGE (β = γ = 0), kGE
×
(α = γ = 0), and kGE (α = β = 0).

Hyperparameters are tuned via maximum likelihood (empirical Bayes). A coarse grid search
provides first estimates to initialize the Adam optimizer (Kingma & Ba 2014), run with a
learning rate of 0.01 decayed by 0.8 every 5 steps, for up to 1000 iterations. In order to
significantly decrease computation time, we use batching of size 0.5, computing the mean
gradient from two randomly selected, disjoint batches, each the size of half the training
dataset. For kernels with discrete parameters, we fix them via grid search and furthermore,
we restrict the length scales θG and θE to the training set’s maximum point distance.

The main parts of the R code used in this study can be found in a GitHub (GitHub 2025).
We implemented the GP model manually, for the spectrum kernel, we used the kernlab
package (Hornik et al. 2024), and for the global alignment kernel, we used the code available
on Marco Cuturi’s personal website. In the Hamming distance, we treat missing letters
as matches but heterozygotes as separate instances. As the meteorological variables have
different orders of magnitude, we normalize them within each time period.

3.4 Comparison approaches

As stressed in Section 2.3, a Bayesian LMM with a covariance structure on environments and
genotypes is closely linked to a GP model. To underline this, we run the GP with the state-
of-the-art kernels using the R-package for Bayesian Generalized Linear Regression (BGLR;
Pérez & de Los Campos 2014), which performs a fully Bayesian analysis and samples the
hyperparameters, effects and predictions with Gibbs. For the hyperparameters of the kernels
(length scale), we use the ML fit of the GP model. Furthermore, we consider a traditional
(frequentist) LMM, with fixed effects on the environmental covariates and (independent)
group random effects on the variety (and environment). Finally, we compare our GP models
to baseline averages: the global average, the environmental average (over all points of the
same environment), and the variety average (over points of the same variety as the prediction
target). These methods are summarized together with the GP approaches in Table 1.

13
 Table 1: Considered prediction methods.

Method Fixed Random Estimation/prediction

GP m GPG : G
kGE Empirical Bayes:
GPE : E
kGE
+ ML for hyperparameters,
GP+ : kGE
× Analytical for prediction
GP× : kGE
(Sct. 2.1)
GP∼ : ∼
kGE
GP1 : kE:GAU−EUCL , kG:GAU−GBLUP
GP2 : kE:GAU−EUCL , kG:EXP−HAM
GP3 : kE:GAU−EUCL , kG:SPE
GP4 : kE:EXP−EUCL , kG:GAU−GBLUP
GP5 : kE:EXP−EUCL , kG:EXP−HAM
GP6 : kE:EXP−EUCL , kG:SPE
GP7 : kE:GAK , kG:GAU−GBLUP
GP8 : kE:GAK , kG:EXP−HAM
GP9 : kE:GAK , kG:SPE

BGLR m BGLR∼ : g G + g E + g GE Full Bayes with Gibbs

kE:GAU−EUCL , kG:GAU−GBLUP

LMM xE Independent group effects Frequentist:

LMM1 : variety
REML for hyperparam.,
LMM2 : variety and environment
BLUP for prediction
Average GLOA All training data
VARA Data from same variety
ENVA Data from same environment

14
4 Results
4.1 New environment
In this section, we are interested in predicting for a new environment. We begin by evaluat-
ing GP∼ using the full kernel combination kGE ∼
(Eq. 1) in conjunction with the nine different
kernel configurations introduced in Table 1. Table 2 summarizes the results for yield and
protein content by depicting the median MSE, CRPS and logS. First we consider the scenario
without leakage (left of vertical bars). For both traits and with respect to all three metrics,
the GPs using the exponential-Euclidean kernel for the environmental information (GP∼ 4,
∼ ∼ ∼ ∼ ∼
GP5 , GP6 ) outperform the ones employing the Gaussian-Euclidean (GP1 , GP2 , GP3 ) and
the GAK-based kernel (GP∼ ∼ ∼ ∼ ∼
7 , GP8 , GP9 ). Among the models GP4 , GP5 and GP6 , we
∼

observe minimal performance differences between the genetic kernels. When allowing one
controlled leakage point per environment (right of bars), the scores improve drastically. In
this scenario, there is generally less difference between the kernels, especially for yield pre-
diction. However, regarding protein, the GAK-based kernel for the environmental variables
(GP∼ ∼ ∼ ∼
7 , GP8 , GP9 ) and also the very first kernel (GP1 ) perform worse than the others.

Next, we compare the different kernel combination strategies within GPG , GPE , GP+ , GP× ,
and GP∼ . Thereby, we focus on GP5 , which demonstrated convincing performance in the
first comparison of the kernels within GP∼ . To simplify the interpretation, we complement
Table 2 with boxplots depicting the MSE values of selected methods for protein content
and yield prediction (Fig. 2a+b). We observe that GP+ , GP× , and GP∼ generally exhibit
comparable performance, while GPG and GPE perform worse than the others. While for
protein and without leakage (Fig. 2b, blue bars), the GP using the genetic kernel only (GPG5)
E
performs clearly better than the one using the environmental kernel only (GP5 ), for yield
(Fig. 2a, blue bars) the performances are similar. When considering a single controlled leak-
age observation of the target environment (red bars), the performance of the environmental
kernel only (GPE 5 ) improves and becomes clearly better than that of the genetic only, even
more so for yield than for protein content. Also the medians in CRPS and logS of Table 2
show very high scores for the GP using the genetic kernel only (GPG 5 ). The GP using the
E
combinations of kernels still performs better than GP5 and for the leakage case, we notice
that GP+ and GP∼ perform better than GP× .

15
Table 2: Model evaluation results for predictions in a new environment without | with one
controlled leakage point. Y represents yield (dt per ha) and P represents grain protein
content. The analysis is based on 30 different train/test splits and we show the median of
the obtained metrics. For the details regarding the methods we refer to Table 1.

Method MSE Y CRPS Y logS Y MSE P CRPS P logS P

GP∼
1 127.91 | 43.34 6.31 | 3.67 3.88 | 3.34 1.54 | 0.56 0.69 | 0.43 1.64 | 1.18
GP∼
2 128.22 | 41.82 6.35 | 3.58 3.88 | 3.31 1.54 | 0.57 0.69 | 0.43 1.63 | 1.19
GP∼
3 130.10 | 41.59 6.32 | 3.58 3.89 | 3.30 1.53 | 0.57 0.69 | 0.43 1.64 | 1.17
GP∼
4 119.72 | 41.14 6.09 | 3.55 3.83 | 3.30 1.43 | 0.55 0.67 | 0.42 1.60 | 1.17
GP∼
5 119.77 | 41.09 6.09 | 3.54 3.83 | 3.29 1.43 | 0.55 0.67 | 0.42 1.60 | 1.16
GP∼
6 119.09 | 41.16 6.10 | 3.54 3.83 | 3.29 1.43 | 0.55 0.67 | 0.42 1.60 | 1.15
GP∼
7 151.40 | 47.20 6.95 | 3.76 3.94 | 3.39 1.74 | 0.60 0.75 | 0.44 1.70 | 1.22
GP∼
8 151.43 | 47.27 6.95 | 3.76 3.94 | 3.39 1.75 | 0.60 0.75 | 0.44 1.70 | 1.20
GP∼
9 151.70 | 47.48 6.96 | 3.75 3.93 | 3.37 1.75 | 0.60 0.75 | 0.44 1.70 | 1.18
GP+
5 125.21 | 45.69 6.22 | 3.78 3.85 | 3.50 1.44 | 0.59 0.68 | 0.44 1.62 | 1.35
GP×
5 121.12 | 54.71 6.16 | 4.12 3.97 | 3.42 1.48 | 0.70 0.69 | 0.47 1.73 | 1.25
GPG
5 147.51 | 149.34 8.68 | 8.80 34.04 | 34.54 1.72 | 1.72 0.96 | 0.95 29.33 | 29.44
GPE
5 152.61 | 82.53 6.90 | 5.32 3.98 | 3.92 1.95 | 1.41 0.81 | 0.70 1.87 | 2.01
GLOA 177.91 | 177.75 10.48 | 10.45 - 2.34 | 2.35 1.21 | 1.21 -
VARA 148.19 | 148.28 9.62 | 9.63 - 1.70 | 1.71 1.06 | 1.06 -
ENVA —— | 96.59 —– | 9.63 - —– | 1.81 —– | 1.06 -
LMM1 156.51 | 148.56 9.66 | 9.40 - 1.70 | 1.63 1.03 | 1.01 -
LMM2 142.66 | 46.05 9.34 | 5.14 - 1.66 | 0.61 1.02 | 0.61 -
∼
BGLR 139.95 | 46.94 6.56 | 3.81 3.93 | 3.53 1.69 | 0.60 0.73 | 0.45 1.71 | 1.38

16
 (a) (b)

(c) (d)

Figure 2: Boxplots of the MSE values for prediction of (a,c) yield and (b,d) protein content.
The first line (a,b) is for the setting considering a new environment and the second line (c,d)
for a new variety. The blue bars depict the values when no leakage is considered, and the red
bars the ones where we use one controlled leakage observation per environment or variety.
The boxplots show the values across 30 different train/test.

To gain a deeper understanding of the inner workings of the GP approaches, Figure 3 presents
a boxplot of the estimated hyperparameters for the different kernel combinations. We only
show the hyperparameters for the scenario targeting yield without leakage, but the results
are very similar for protein content and do not drastically change when adding one leakage
observation. We immediately observe that θG is always estimated at the boundary value of
one, suggesting that the ML approach aims to flatten the effect of the genetic kernels. In
contrast, the length scale θE of the environmental kernel is consistently estimated around
0.25 for all GP+ , GP× , and GP∼ models with kernels 1 and 5. The coefficients for the com-
binations in GP+5 assign approximately α = 0.65 of the weight to the environmental kernel.
In the full model, both GP∼ ∼
1 and GP5 allocate about γ = 0.2 of the weight to the kernel
product. Specifically, GP∼ 1 assigns roughly α = 0.35 to the genetic kernel and β = 0.45
to the environmental kernel, while GP∼ 5 assigns about α = 0.25 to the genetic kernel and

17
Figure 3: Estimated hyperparameters for the different GP models, fitted using the Adam
optimizer. The boxplot shows the values across 30 different train/test splits for a new
environment, without leakage, and for the trait yield.

α = 0.55 to the environmental kernel. Finally, the proportion of variance explained by the
GP, denoted ς, is estimated similarly across all methods. For GP5 , it is lowest for GP+ ,
followed by GP∼ , and highest for GP× .

Figure 2 and Table 2 also show the scores of the considered state-of-the-art methods: BGLR∼ ,
LMM1 and LMM2 . As expected, GP∼ ∼
1 and BGLR yield similar results for both traits and
leakage scenarios. Also, LMM1 (without a random effect for the environment) generally
performs worse than most methods, particularly in the leakage scenario. In the controlled
leakage situation, the MSE performance of LMM2 (including a random effect for the envi-
ronment) improves substantially, though it still falls slightly short of the best GP models.
Note that the LMMs traditionally do not account for uncertainty, so the logS is missing, and
the CRPS values based on the MAE are worse than those of models that explicitly model
uncertainty. Finally, we consider the averaging approaches. The global average is clearly
outperformed by all methods. The variety and environmental averages perform better, but
still clearly worse than the best GP approaches. Interestingly, the variety average performs
about as well as the GPG 5 approaches, indicating a lack of predictive power in the genetic
kernels. Conversely, the GPE approach outperforms the environmental average, suggesting
more effective kernel choices.

4.2 New variety

In this section, we consider predicting the performance of a new variety. Table 3 and Fig-
ure 2c+d report the same assessment metrics as in the previous section. First, we observe
that performance in the no-leakage case is substantially better than in the setting targeting

18
a new environment. The difference between leakage and no leakage is less pronounced here,
although for the LMM models it still has the largest impact. The scores for GP∼ using the
nine kernels in Table 3 show that, in this setting, the spectrum kernel for genotype (GP∼ 3,
GP∼ ∼ ∼ ∼ ∼
6 , GP9 ) performs worse than both the Gaussian-GBLUP (GP1 , GP4 , GP7 ) and the
exponential-Hamming (GP∼ ∼ ∼
2 , GP5 , GP8 ). Here, the choice of environmental kernel makes
little difference. Once again, the differences between the kernels diminish under leakage.

We also compare again GPG , GPE , GP+ , GP× , and GP∼ , focusing on GP5 . Figure 2c+d
show little difference between the combinations GP+ , GP× , and GP∼ . For yield, GPE
performs remarkably well, and also for protein it still performs reasonably. GPG performs
very poorly in this setting for both traits, likely due to the lack of training data from the
same variety. As before and as expected, GP∼ ∼
1 and BGLR produce very similar results.
LMM1 again cannot compete, while LMM2 performs well in the leakage scenario. Finally,
considering the averaging approaches in Table 3, the global average performs poorly here
as well. The environmental average performs similarly to the GPE approaches when not
accounting for uncertainty quantification (MSE). In the leakage case, the variety average
relies on one single measurement of the same variety and performs very poorly.

5 Discussion and Conclusions

Optimizing multi-environment trials and recommending varieties based on the local envi-
ronment are two very common themes in agricultural research with significant potential for
application. This paper introduces a GP modeling approach for G×E prediction, and links
it to state-of-the-art Bayesian LMM methods. Thereby, we explore the impact of kernel
choice and kernel combination, and confirm that GP models and Bayesian LMMs operate in
essentially the same way. In doing so, we not only test the Gaussian kernels already used,
but also alternative kernels specifically designed for the structure of the data at hand.

There is considerable scope to discuss and improve the choice of kernels within the GP mod-
eling. We observed that the performance of the state-of-the-art approaches (BGLR∼ and
GP∼ 1 ) could be improved using alternative kernels, primarily when targeting new environ-
ments (Tab. 2). However, improvements did not come from the more complex kernels for
time series and strings, but rather from the simpler exponential kernel approaches (GP∼
5 ). In
particular, the time series kernel based on global alignment performed poorly when predict-
ing for a new environment (Tab. 2). This may be because, for wheat performance, not only

19
Table 3: Model evaluation results for predictions for a new variety without | with one
controlled leakage point. Y represents yield (dt per ha) and P represents grain protein
content. The analysis is based on 30 different train/test splits and we show the median of
the obtained metrics. The methods are detailed in Table 1.
Method MSE Y CRPS Y logS Y MSE P CRPS P logS P
GP∼
1 51.66 | 39.04 3.99 | 3.46 3.42 | 3.30 0.65 | 0.48 0.45 | 0.39 1.23 | 1.10
GP∼
2 52.42 | 39.48 4.01 | 3.47 3.42 | 3.28 0.65 | 0.48 0.46 | 0.39 1.22 | 1.08
GP∼
3 62.37 | 43.17 4.35 | 3.64 3.48 | 3.33 0.74 | 0.49 0.49 | 0.4 1.27 | 1.08
GP∼
4 52.53 | 38.99 4.02 | 3.45 3.42 | 3.3 0.66 | 0.49 0.45 | 0.39 1.22 | 1.10
GP∼
5 51.96 | 39.19 4.03 | 3.46 3.41 | 3.28 0.66 | 0.48 0.46 | 0.39 1.21 | 1.09
GP∼
6 61.48 | 41.24 4.39 | 3.59 3.49 | 3.33 0.73 | 0.5 0.48 | 0.4 1.26 | 1.09
GP∼
7 51.5 | 38.96 3.98 | 3.47 3.44 | 3.35 0.66 | 0.48 0.46 | 0.4 1.24 | 1.13
GP∼
8 53.21 | 39.6 4.03 | 3.48 3.44 | 3.34 0.65 | 0.48 0.46 | 0.39 1.22 | 1.10
GP∼
9 62.32 | 43.7 4.42 | 3.66 3.5 | 3.33 0.76 | 0.5 0.49 | 0.4 1.27 | 1.10
GP+
5 55.67 | 40.36 4.16 | 3.61 3.62 | 3.69 0.67 | 0.5 0.46 | 0.41 1.35 | 1.33
GP×
5 52.33 | 44.69 4.00 | 3.70 3.44 | 3.38 0.72 | 0.56 0.48 | 0.42 1.28 | 1.16
GPG
5 159.91 | 158.79 7.92 | 7.98 6.42 | 6.83 2.07 | 2 0.9 | 0.89 3.33 | 3.85
GPE
5 53.75 | 53.48 4.77 | 4.76 7.86 | 7.89 0.96 | 0.95 0.65 | 0.65 6.39 | 6.48
GLOA 177.71 | 177.77 10.48 | 10.49 - 2.40 | 2.40 1.25 | 1.25 -
VARA —— | 257.13 —— | 12.57 - —— | 3.07 —— | 1.41 -
ENVA 53.73 | 53.21 5.76 | 5.76 - 0.94 | 0.94 0.77 | 0.77 -
LMM1 144.07 | 129.05 9.33 | 8.81 - 2.04 | 1.68 1.14 | 1.02 -
LMM2 86.85 | 47.16 7.42 | 5.36 - 1.02 | 0.57 0.83 | 0.60 -
∼
BGLR 57.88 | 40.22 4.23 | 3.61 3.66 | 3.68 0.68 | 0.51 0.48 | 0.41 1.40 | 1.35

20
the sequence of weather events matters, but also their specific timing within the growing
cycle. Similar behavior was observed for the spectrum kernel for strings when predicting a
new variety (Tab. 3), suggesting the use of alternative kernels. The spectrum kernel con-
sidered here is also quite simplistic and has been extended, for instance, by Eskin et al. (2002).

We observe that the kernel combination approaches (GP+ , GP× , and GP∼ ) overall did not
lead to substantial differences. Only in the leakage case did the additive model GP+ perform
worse than the others. Across all settings, the full model GP∼ consistently performed slightly
better. In general, the genotype kernels perform poorly when predicting new varieties. As
shown in Fig. 2c and d, the environmental-only model GPE is not substantially improved
by adding genotype information in GP∼ , while the genotype-only model GPG performs al-
most as poorly as the global average (Tab. 3). Also in the setting concerned with a new
environment, the length scale θG is maximized to flatten out the genetic effect as much as
possible (Fig. 3), and the variety average performs about as well as GPG . This performance,
especially in the scenario predicting for a new variety, indicates that the considered genetic
kernels effectively only ‘work’ for the same variety, meaning that their measure of genetic
similarity is not very effective in-between varieties. This clearly highlights the potential of
exploring alternative kernels beyond the Gaussian-GBLUP and our proposed kernels.

The proposed environmental kernel kE:EXP−EUCL performs better, and in Figure 2a and b,
which target a new environment, using a kernel combination clearly outperforms using a
single kernel. However, including one leakage observation from the same environment sub-
stantially improves prediction accuracy, highlighting the importance of site- and year-specific
management and soil characteristics for wheat performance. This effect of local information
is also evident in the stronger prediction results for new varieties (Table 3), compared to pre-
diction for new environments (Table 2). In this setting, all available local environmental data
can be leveraged, providing insights into local management and soil characteristics, which
leads to good prediction performance even without variety leakage. To mitigate this effect,
incorporating additional environmental covariates describing soil and management practices,
as suggested by Buntaran et al. (2021), could be beneficial. Extending the GP with kernels
for these inputs enables a comprehensive model integrating genetic, environmental (weather
and soil) and management information.

We hope this paper serves as a useful proof of concept demonstrating that the GP frame-
work, with its flexibility for handling diverse input spaces, provides a promising foundation
for improved G×E predictions. We see various potential applications: Local variety rec-

21
ommendations at the farm level with the GP model enabling the inclusion of new varieties
from the official Swiss recommended list (Strebel et al. 2025), investigation of future weather
scenarios during the winter wheat growing season and improvement of multi-environment
trials. Moreover, GPs enable joint modeling of yield and protein content, allowing for local
analysis of their relationship. Finally, GP modeling opens the door to sequential design
strategies (Chevalier et al. 2014), and, as it inherently provides uncertainty quantification,
to risk-averse strategies and decision-aid tools.

Data and Code The anonymized wheat dataset used in this work is available on Zenodo
(Levy Häner et al. 2025). Key code is available on a Github repository (link).

Acknowledgements Lea Friedli acknowledges support by the Swiss National Science

Foundation (grant number: 225353). Tim Steinert and David Ginsbourger acknowledge the
support of the Digitization Commission (DigiK) of the University of Bern via the project
“Perception in Statistics, Econometrics and Probability”. Part of this research was per-
formed while Tim Steinert was visiting the Institute for Mathematical and Statistical Inno-
vation (IMSI), supported by the National Science Foundation (Grant No. DMS-1929348).
David Ginsbourger would like to thank the Isaac Newton Institute for Mathematical Sci-
ences, Cambridge, for support and hospitality during the program Representing, calibrating
& leveraging prediction uncertainty from statistics to machine learning, where work on this
paper was undertaken that was partially supported by EPSRC grant EP/Z000580/1 and
by a grant from the Simons Foundation. Nathalie Wuyts, Lilia Levy Häner, Didier Pellet
and Juan M. Herrera acknowledge support by Agroscope, swiss granum, the Swiss Federal
Office for Agriculture [project ’Wheat Advisor’, grant no. 19.03], the Schweizerischer Getrei-
deproduzentenverband, Prometerre, Fresh Food & Beverage Group and Timac Agro Swiss.
Computations were performed on UBELIX ([Link] the HPC cluster
at the University of Bern

References
Bandeira e Sousa, M., Cuevas, J., de Oliveira Couto, E. G., et al. (2017). Genomic-enabled
prediction in maize using kernel models with genotype×environment interaction. G3:
Genes Genomes Genet., 7(6), 1995–2014. [Link]
Berg, C., Christensen, J. P. R., & Ressel, P. (1984). Harmonic analysis on semigroups:
theory of positive definite and related functions, volume 100. Springer. [Link]
10.1007/978-1-4612-1128-0.

22
Buntaran, H., Forkman, J., & Piepho, H.-P. (2021). Projecting results of zoned multi-
environment trials to new locations using environmental covariates with random coefficient
models: accuracy and precision. Theor. Appl. Genet., 134, 1513–1530. [Link]
10.1007/s00122-021-03786-2.
Burgueño, J., Crossa, J., Cornelius, P. L., et al. (2007). Modeling additive×environment and
additive×additive×environment using genetic covariances of relatives of wheat genotypes.
Crop Sci., 47(1), 311–320. [Link]
Burgueño, J., de los Campos, G., Weigel, K., & Crossa, J. (2012). Genomic prediction
of breeding values when modeling genotype× environment interaction using pedigree
and dense molecular markers. Crop Sci., 52(2), 707–719. [Link]
cropsci2011.06.0299.
Chevalier, C., Bect, J., Ginsbourger, D., et al. (2014). Fast parallel kriging-based stepwise
uncertainty reduction with application to the identification of an excursion set. Techno-
metrics, 56(4), 455–465. [Link]
Cornish-Bowden, A. (1985). Nomenclature for incompletely specified bases in nucleic acid
sequences: rcommendations 1984. Nucleic Acids Research, 13(9), 3021–3030. https:
//[Link]/10.1093/nar/13.9.3021.
Costa-Neto, G., Fritsche-Neto, R., & Crossa, J. (2021). Nonlinear kernels, dominance, and
envirotyping data increase the accuracy of genome-based prediction in multi-environment
trials. Heredity, 126(1), 92–106. [Link]
Crossa, J., Burgueño, J., Cornelius, P. L., et al. (2006). Modeling genotype×environment
interaction using additive genetic covariances of relatives for predicting breeding values of
wheat genotypes. Crop Sci., 46(4), 1722–1733. [Link]
11-0427.
Crossa, J., Montesinos-Lopez, O. A., Pérez-Rodrı́guez, P., et al. (2022). Genome and envi-
ronment based prediction models and methods of complex traits incorporating genotype ×
environment interaction. Genomic Prediction of Complex Traits: Methods and Protocols,
(pp. 245–283). [Link]
Cuevas, J., Crossa, J., Montesinos-López, O. A., et al. (2017). Bayesian genomic prediction
with genotype×environment interaction kernel models. G3: Genes Genomes Genet., 7(1),
41–53. [Link]
Cuevas, J., Crossa, J., Soberanis, V., et al. (2016). Genomic prediction of genotype× en-
vironment interaction kernel regression models. Plant Genome, 9(3). [Link]
10.3835/plantgenome2016.03.0024.
Cuevas, J., Montesinos-López, O., Juliana, P., et al. (2019). Deep kernel for genomic and near
infrared predictions in multi-environment breeding trials. G3: Genes Genomes Genet.,
9(9), 2913–2924. [Link]
Cuturi, M., Vert, J.-P., Birkenes, O., & Matsui, T. (2007). A kernel for time series based
on global alignments. In 2007 IEEE International Conference on Acoustics, Speech and
Signal Processing-ICASSP’07, volume 2 (pp. II–413).: IEEE.

23
Eskin, E., Weston, J., Noble, W., & Leslie, C. (2002). Mismatch string kernels for svm
protein classification. Adv. Neural Inf. Process. Syst., 15.
Fernandes, I. K., Vieira, C. C., Dias, K. O., & Fernandes, S. B. (2024). Using machine
learning to combine genetic and environmental data for maize grain yield predictions
across multi-environment trials. Theor. Appl. Genet., 137(8), 189. [Link]
1007/s00122-024-04687-w.
Garnett, R. (2023). Bayesian optimization. Cambridge University Press. [Link]
org/10.1017/9781108348973.
Ginsbourger, D., Baccou, J., Chevalier, C., et al. (2014). Bayesian adaptive reconstruction of
profile optima and optimizers. SIAM/ASA J. Uncertain. Quantif., 2(1), 490–510. https:
//[Link]/10.1137/130949555GetAccessBibTeX.
GitHub (2025). Wheat advisor. [Link]
anonymous-E380/. GitHub.
Gneiting, T. & Raftery, A. E. (2007). Strictly proper scoring rules, prediction, and
estimation. J. Am. Stat. Assoc., 102(477), 359–378. [Link]
016214506000001437.
Good, I. J. (1952). Rational decisions. J. R. Stat. Soc. Ser. B, 14(1), 107–114.
Griffiths, R.-R., Klarner, L., Moss, H., et al. (2023). Gauche: a library for gaussian processes
in chemistry. Adv. Neural Inf. Process. Syst., 36, 76923–76946.
Guignard, F., Ginsbourger, D., Levy Häner, L., & Herrera, J. M. (2024). Some combinatorics
of data leakage induced by clusters. Stoch. Environ. Res. Risk Assess., (pp. 1–14). https:
//[Link]/10.1007/s00477-024-02715-1.
Handcock, M. S. & Stein, M. L. (1993). A Bayesian analysis of kriging. Technometrics,
35(4), 403–410. [Link]
Herrera, J. M., Häner, L. L., Holzkämper, A., & Pellet, D. (2018). Evaluation of ridge
regression for country-wide prediction of genotype-specific grain yields of wheat. Agr.
For. Meteorol., 252, 1–9. [Link]
Hornik, K., Zeileis, A., & Buchta, C. (2024). kernlab: Kernel-based Machine Learning Lab.
R package version 0.9-32.
Hutter, F., Xu, L., Hoos, H. H., et al. (2014). Algorithm runtime prediction: Methods &
evaluation. Artif. Intell., 206, 79–111. [Link]
003.
Janusevskis, J. & Le Riche, R. (2013). Simultaneous kriging-based estimation and optimiza-
tion of mean response. J. Glob. Optim., 55(2), 313–336. [Link]
s10898-011-9836-5.
Jarquı́n, D., Crossa, J., Lacaze, X., et al. (2014). A reaction norm model for genomic
selection using high-dimensional genomic and environmental data. Theor. Appl. Genet.,
127, 595–607. [Link]
Jones, D. R., Schonlau, M., & Welch, W. J. (1998). Efficient global optimization of expensive
black-box functions. J. Glob. Optim., 13(4), 455–492. [Link]

24
 1008306431147.
Kang, M. & Gorman, D. (1989). Genotype × environment interaction in maize. J. Agron.,
81(4), 662–664. [Link]
Kingma, D. P. & Ba, J. (2014). Adam: A method for stochastic optimization. CoRR,
abs/1412.6980.
Krige, D. G. (1951). A statistical approach to some basic mine valuation problems on the
witwatersrand. J. S. Afr. Inst. Min. Metall., 52(6), 119–139.
Leslie, C., Eskin, E., & Noble, W. S. (2001). The spectrum kernel: A string kernel for svm
protein classification. In Biocomputing 2002 (pp. 564–575). World Scientific.
Levy Häner, L., Strebel, S., Visse-Mansiaux, M., Wuyts, N., Herrera, J. M., & Pellet, D.
(2025). Winter wheat performance data for varieties grown under low input conditions in
switzerland over the period 1990–2023, together with trial site weather data and variety
genetic marker data (1.0.0) [data set]. [Link]
Liu, J., Gock, A., Ramm, K., et al. (2025). Incorporating gene expression and environment
for genomic prediction in wheat. Front. Plant Sci., 16, 1506434. [Link]
3389/fpls.2025.1506434.
Matheron, G. (1969). Le krigeage universel (Universal kriging), volume 1 of Cahiers du
Centre de Morphologie Mathématique. Fontainebleau: École des Mines de Paris. https:
//[Link]/10.1016/[Link].2003.08.018.
Mercer, J. (1909). Functions of positive and negative type, and their connection with the
theory of integral equations. Philos. Trans. R. Soc. A, 209, 415–446. [Link]
10.1098/rsta.1909.0016.
MeteoSwiss (2024a). Automatic measurement network. [Link]
[Link]/weather/measurement-systems/land-based-stations/
[Link].
MeteoSwiss (2024b). Spatial climate analyses. [Link]
climate/the-climate-of-switzerland/[Link].
Meuwissen, T. H., Hayes, B. J., & Goddard, M. (2001). Prediction of total genetic value
using genome-wide dense marker maps. Genetics, 157(4), 1819–1829. [Link]
10.1093/genetics/157.4.1819.
Močkus, J. (1974). On Bayesian methods for seeking the extremum. In IFIP Tech. Conf.
Optim. Techn. (pp. 400–404).: Springer. [Link]
Murphy, A. H. (1973). A new vector partition of the probability score. J. Appl. Meteo-
rol. Climatol., 12(4), 595–600. [Link]
ANVPOT>[Link];2.
Pérez, P. & de Los Campos, G. (2014). Genome-wide regression and prediction with
the BGLR statistical package. Genetics, 198(2), 483–495. [Link]
genetics.114.164442.
Picheny, V., Ginsbourger, D., Richet, Y., & Caplin, G. (2013). Quantile-based optimization
of noisy computer experiments with tunable precision. Technometrics, 55(1), 2–13. https:

25
 //[Link]/10.1080/00401706.2012.707580.
Piepho, H.-P. (1997). Analyzing genotype-environment data by mixed models with multi-
plicative terms. Biometrics, (pp. 761–766). [Link]
Prus, M. & Piepho, H. P. (2024). Optimizing the allocation of trials to sub-regions in
crop variety testing with multiple years and locations. J. Agric. Biol. Environ. Stat.
[Link]
Rafalski, A. (2002). Applications of single nucleotide polymorphisms in crop genetics. Curr.
Opin. Plant Biol., 5(2), 94–100. [Link]
Rasmussen, C. E. & Williams, C. K. I. (2006). Gaussian Processes for Machine Learning.
The MIT Press. [Link]
Ray, A. (2021). Bayesian inversion using nested trans-dimensional Gaussian processes. Geo-
phys. J. Int., 226(1), 302–326. [Link]
Roustant, O., Ginsbourger, D., & Deville, Y. (2012). DiceKriging, DiceOptim: Two R
packages for the analysis of computer experiments by kriging-based metamodeling and
optimization. J. Stat. Softw., 51(1), 1–55. [Link]
Sakoe, H. & Chiba, S. (1970). A similarity evaluation of speech patterns by dynamic pro-
gramming. In IEICE National Meeting (pp. 136).
Sanders, F. (1963). On subjective probability forecasting. J. Appl. Meteorol. Climatol., 2(2),
191–201. [Link]
Schulz-Streeck, T., Ogutu, J. O., Gordillo, A., et al. (2013). Genomic selection allowing for
marker-by-environment interaction. Plant Breeding, 132(6), 532–538. [Link]
10.1111/pbr.12105.
Stein, M. L. (2012). Interpolation of Spatial Data: Some Theory for Kriging. Springer
Science & Business Media. [Link]
Strebel, S., Levy Häner, L., Imhoff, Y., et al. (2025). Liste der empfohlenen Getreidesorten
für die Ernte 2026. Agroscope Transfer, 591. [Link]
publication/59585.
Su, G., Peng, L., & Hu, L. (2017). A Gaussian process-based dynamic surrogate model
for complex engineering structural reliability analysis. Struct. Saf., 68, 97–109. https:
//[Link]/10.1016/[Link].2017.06.003.
Tadese, D., Piepho, H.-P., & Hartung, J. (2024). Accuracy of prediction from multi-
environment trials for new locations using pedigree information and environmental co-
variates: the case of sorghum (Sorghum bicolor (L.) Moench) breeding. Theor. Appl.
Genet., 137(8), 181. [Link]
Tanaka, R. & Iwata, H. (2018). Bayesian optimization for genomic selection: a method for
discovering the best genotype among a large number of candidates. Theor. Appl. Genet.,
131(1), 93–105. [Link]
Tsai, S. F., Shen, C. C., & Liao, C. T. (2021). Bayesian optimization approaches for iden-
tifying the best genotype from a candidate population. J. Agric. Biol. Environ. Stat., 26,
519–537. [Link]

26
VanRaden, P. (2007). Genomic measures of relationship and inbreeding. Interbull Bull.,
(37), 33–33.
Williams, B. J., Santner, T. J., & Notz, W. I. (2000). Sequential design of computer
experiments to minimize integrated response functions. Stat. Sin., 10(4), 1133–1152.
[Link]

27