BACKGROUND
genomic prediction  was first introduced in  <dig>  <cit>  as a method that allows the prediction of genomic estimated breeding values  for plants and animals by using information of genetic markers.
in plant breeding, gp has been adopted as another stage of the breeding scheme  <cit> , not diminishing the importance of the phenotypic analysis usually carried out in several environments.
merging the phenotype and the genotype analyses has been addressed through the so-called stage-wise analysis  <cit> .
in the first stage environments are analysed separately and genotype means are computed and then submitted in the gp stage to predict gebv based on dense genetic markers such as single nucleotide polymorphisms .
in plant breeding, assessing genotypic adaptability and stability, and predicting breeding values of the genotypes in other environments and other years, makes use of multi-environment trials , which aim to evaluate as many genotypes as possible in as many as possible locations .
these mets are typically laid out as generalised lattice designs testing a large number of different genotypes per trial.
the number of tested genotypes is limited by factors such as seed production, production cycle length and availability of physical resources, e.g.
land and budget  <cit> .
within years, genotypes are tested in series of trials, which are connected by checks.
checks are lines grown in every trial as controls because their performance is known and/or they are already commercial material.
checks can be also used to connect years.
in the rye breeding program considered in this paper, a completely different set of genotypes is tested in each year, but these genotypes are from the same breeding population.
the accuracy of a genomic prediction model depends on the number of genotypes used for calibration.
so there is definitely a need to combine data across years.
low connectivity across years is a challenge when trying to combine data across years, and this is one main motivation for this paper.
furthermore, the unbalancedness due to the design layout and the different and large number of evaluated genotypes increases the heterogeneity introducing high complexity to the variance-covariance structure among adjusted genotype means  <cit> .
analysis of mets could be done as single-stage analysis, modelling the complete observed data at the level of individual plots, or using a stage-wise approach, where experiments are analysed first at the level of environments , obtaining adjusted means per genotype, which are then summarised across environments  in the next stage  <cit> .
a single-stage analysis accounts entirely for the variance-covariance structure of the recorded observations  <cit> , therefore it is regarded as the gold standard.
however, it has been shown that in a stage-wise analysis, a loss of information occurring in the transition through stages can be minimized by an appropriate weighting scheme  <cit> .
if feasible, a single-stage approach is preferable to a stage-wise analysis  <cit> .
nevertheless, the latter is acceptable for gp, since it is simple, computationally more efficient and also allows to easily account for any specifics of randomisation layout and error modelling for each environment  <cit> .
it should be stressed, however, that in a stage-wise analysis the weights are chosen to approximate the variance-covariance matrix of adjusted means from previous stages.
we used here a three-stage approach and compared different spatial correlation structures in the first stage to correct field heterogeneity at the trial level.
spatial error models may provide more accurate estimates of genotype effects than models not accounting for spatial adjustment  <cit>  but they are computationally more demanding and convergence may be difficult to reach.
any effort in terms of improving the genomic predictions would include checking if these improved estimates have an effect on the predictive ability when markers are added to the model.
the performance of alternative spatial models can be assessed by k-fold cross validation .
similarly, the merits of different spatial models used to compute adjusted means in the first stage can be compared by the same cv procedure, if the same gp procedure is used for each analysis.
this suggests that genomic prediction-cross validation  can be used to identify the best-fitting mixed model in stage one.
the common method of model selection makes use of information criteria based on the log likelihood, e.g.
the akaike information criterion  or the bayesian information criterion   <cit> .
when the restricted maximum likelihood  method is used, models can only be compared by information criteria if they have the same fixed effects; otherwise, the maximum likelihood  method should be used  <cit> .
cv is, in this sense, not used to tune parameters as in many penalization methods , machine learning methods) but only as a tool to compare models that use reml.
reml is considered the best available method of variance parameter estimation, preferable to ml  <cit> .
consequently, it is of interest to devise model selection procedures that can use reml and also can compare models with different fixed effects.
gp-cv has already been used to judge environments in order to optimise the accuracy in gp  <cit> .
we used this tool here as model selection method in comparison to the traditional use of aic.
the aims of this work were: i) to assess the advantage for the predictive ability when using a spatial model for phenotypic analysis, ii) to compare stage-wise approaches for gp when the data are weakly connected across years, and iii) to compare aic and gp-cv as methods of selection of models for phenotypic data analysis towards gp in rye.
methods
field layout and data set
a commercial rye breeding program by kws-lochow established in poland and germany aims to develop superior hybrid varieties for the seed market.
the implementation of gp within the breeding program makes use of the measurements of hybrid performance of the first cycles of phenotypic evaluation of the material .
selections made in cycle <dig> are intensively evaluated in further cycles, aiming to double-check the selection decisions.
for our purposes, these additional cycles do not add much useful information.
hence, we used only the first cycles of the program.
the populations tested in each year consist of s  <dig> genotypes, which display genetic relatedness and population stratification due to complex genealogical history  <cit> .
besides the phenotypic data, a 16k infinium iselect hd custom beadchip was used to characterise  <dig> individuals from cycle1- <dig> and cycle1- <dig> and  <dig> checks.
several traits were evaluated during this project: grain dry matter yield, plant height and thousand kernel weight, as well as ordinal scores of rust, mildew and lodging among others.
in this work we used grain dry matter yield measurements of the phases of selection cycle1- <dig>  cycle1- <dig> and cycle1- <dig>  and marker information for the genotypes of  <dig> and  <dig>  although no marker information of year  <dig> was available, it makes sense to use this dataset to observe the trend in one additional year and in this way, support the results of the phenotypic analysis of previous years.
a cycle <dig> experiment consists of subsets of  <dig> genotypes from the s  <dig> populations tested in several locations within each of the two countries involving two testers .
we define a trial as the physical unit within a location, where a subset of genotypes that were testcrossed to the same tester is evaluated.
trials at a location were laid out as α-designs with two replicates.
each trial was randomized independently from the others using the software cycdesign .
.
in our notation, trials of a cycle <dig> experiment are labelled as s <dig>  s <dig>  …, s <dig>  row and column coordinates of the plots to account for spatial variation are available.table  <dig>
general representation of the testers by locations  by years classification of cycle <dig> year  <dig> and  <dig> in germany  and poland
series of trials are represented with the labels s <dig>  s <dig>  ⋯ s <dig>
general representation of the testers by locations  classification of cycle <dig> year  <dig> in germany  and poland
series of trials are represented with the labels s <dig>  s <dig>  ⋯ s <dig>
normally throughout the program, only a single tester was used per location and year, but in some locations, some subsets of genotypes were testcrossed with the two available testers.
this is the case, for example, for location g-l <dig> in cycle1- <dig>  where the genotypes evaluated in the trials s <dig>  s <dig> and s <dig> were testcrossed with both tester <dig> and tester <dig>  and it is also the case of locations p-l <dig>  p-l <dig>  p-l <dig> and p-l <dig> evaluating genotypes of trials s <dig>  s <dig> and s <dig> with both testers.
in each year, four common checks were testcrossed with the testers and grown twice in each trial.
over the years  <dig> and  <dig> one check was in common and none was shared with  <dig> .table  <dig>
year x check classification in germany  and poland
the field layout of some trials was not perfectly rectangular.
some trials at a given location and year had fewer blocks but larger size, i.e., there were two different block sizes within a few trials.
blocks were nested within rows of the field layout.
in the genetic dataset, homozygous marker genotypes were coded as - <dig> and  <dig>  and the heterozygous type, missing values and technical failures were coded as  <dig>   <dig> % of the markers corresponded to homozygous alleles and  <dig> % were heterozygous.
only a  <dig> % of the markers were recorded as missing values or technical failures; therefore, an imputation method would not have a strong impact on the subsequent analyses.
monomorphic markers and markers with minor allele frequency  less than 1% or missing information of more than 10% per marker were dropped.
a total of  <dig> markers passed the quality test and were used for gp.
models
in this section we present the models used in the first stage of the analysis and the models of the approaches followed to adjust the year effect either in the second or the third stage.
figures  <dig> and  <dig> depict a general scheme that helps visualizing the methodology.figure  <dig>
general representation of stage-wise approaches to compare year-effect adjustment.
factors were genotype , tester , location , year , trial , replicate  and block .
grain dry matter yield  is the response variable in the first stage, m
is the adjusted mean of genotypes across locations used in the second stage, m
is the year effect-corrected genotype adjusted mean,  represents the simple mean of genotypes of the r-th year.
in the genomic prediction  stage, m
is the n× <dig> vector of adjusted means of genotypes by year for approach 1a and across years for approach  <dig>  m
is the n× <dig> vector of adjusted means of year effect-corrected genotypes in approach 1b, x and β are respectively the design matrix and parameter vector of fixed effects, z is the n×p marker matrix, u is the p-dimensional vector of snp effects and e the error vector.
y=g·t:s/r/b is the shorthand notation of the model eq.
in the text: y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv, m
=g×l×t stands for the model eq.
in the text: , and m
=×g×l represents the extended model eq.
in the text: .
the final predictive abilities  are presented in the ellipses.
general representation of model comparison through all the stages of the analysis.
datasets generated from  <dig> spatial and non-spatial models plus two mixed datasets generated from best models given the akaike information criterion  and the predictive abilities .
factors in second stage were genotype , location  and tester .
m
represents the adjusted mean of genotypes across locations and years.
m
=g×l×t is the shorthand notation for .
in the genomic prediction  stage m
is the adjusted mean of genotypes across locations, x and β are respectively the design matrix and parameter vector of fixed effects, z is the n×p marker matrix, u is the p-dimensional vector of snp effects and e the error vector.
sampling methods in cross validation  were across crosses  and within crosses .
the final predictive abilities  are presented in the ellipses.
first stage
in the first stage we computed adjusted genotype means by location and year.
the factors used for the analysis were genotypes , testers , trials , replicates  nested within trials and blocks  nested within replicates.
we defined a baseline model as
<dig>
where yhijkv is the observed grain dry matter yield of the h-th genotype testcrossed with the v-th tester in the k-th block within the j-th replicate of the i-th trial, hv is the effect of the h-th genotype testcrossed with the v-th tester, si is the effect of the i-th trial , rij is the effect of the j-th replicate nested within the i-th trial , bijk is the effect of the k-th block nested within the j-th replicate of the i-th trial  and ehijkv is the plot error associated with the yhijkv observation .
in model equation  we assumed genotypes crossed with testers as a fixed effect to be able to compute genotype adjusted means per tester, whereas the other effects were considered as random effects due to the nested design structure  <cit> .
spatial and non-spatial models used for the first stage
y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv
y
hijkoqv=hv+s
i+r
ij+b
ijk
ijo+v
ijq+e
hijkoqv
y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv
y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv
y
hijkoqv=hv+s
i+r
ij+b
ijk
ijo+v
ijq+e
hijkoqv
y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv
y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv
y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv
y
hijkoqv=hv+s
i+r
ij+b
ijk
ijo+v
ijq+e
hijkoqv
y
hijkv is the observed dry matter yield of the h-th genotype testcrossed with the v-th tester in the k-th block within the j-th replicate of the i-th trial, hv is the effect of the h-th genotype testcrossed with the v-th tester, s
i is the effect of the i-th trial , r
ij is the effect of the j-th replicate nested within the i-th trial , b
ijk is the effect of the k-th block nested within the j-th replicate of the i-th trial  and e
hijkv is the plot error associated with the y
hijkv observation .
in the models including row and column effects, w
ijo is the effect of the o-th row within the j-th replicate of the i-th trial  and v
ijq is the effect of the q-th column within the j-th replicate of the i-th trial .
spatial variance-covariance structure were independent , autoregressive in one direction , one-dimension linear variance  and two-dimension autoregressive .
note that we use hv as fixed effect, which is necessary to obtain the genotype by tester means.
the purpose is also to recover the information of the entries that are grown in the same locations but using different testers , so that we captured the effect of the tester in the shared locations.
second stage
in the second stage we computed genotype means across locations and testers.
this was done either separately for each year  or also averaging across years .
the years  <dig> and  <dig>  where molecular marker data were available, were connected through only one check.
the resulting fundamental question is then how to fit the year effect.
either the year effect is estimated by the mean of all tested entries  or we rely on the adjustment by the one single check .
we assume that genotypes tested in each year can be regarded as a random sample from the same parent population.
based on the structure of the breeding program, this is a realistic assumption that motivates the approaches described in the following.
both approaches were compared using the m resulting from the analysis of the baseline model in the first stage.
approach 1: year-wise analysis
each year was analysed in the second stage using a three-way interaction model of genotypes , locations  and testers  as factors to obtain adjusted genotypes means of each year.
the model was
<dig>
where  represents the adjusted mean of grain dry matter yield of the h-th genotype, testcrossed with the v-th tester in the s-th location, gh, ls and tv are the main effects of the h-th genotype, the s-th location and the v-th tester, respectively, hs, hv and sv are the two-way interaction effects, hsv is the effect of the three-way interaction and ehsv is the residual error associated with , with  the variance of the hsv-th adjusted mean  obtained in the first stage.
location was considered as random effect  and hence, all the interactions containing this factor are random  <cit> .
the crossed effect of genotypes and testers  could have been a fixed effect since genotypes and testers are taken as fixed factors in this stage.
however, the crossed effects that include g were taken as random here because the factor genotype was used as random in the gp stage.
but note that in the first and the second stage we needed to take genotype main effects as fixed in order to compute adjusted means  <cit> .
besides, since not every genotype was tested with every tester , we needed to take hv random to be able to estimate genotype means across levels of testers.
in this approach, the year effect was adjusted in two ways, hereafter referred as to approach 1a and approach 1b.
approach 1a used years as fixed factors in the gp stage and approach 1b used a manual adjustment after the second stage by simply calculating the mean of the genotypes by year  and subtracting it to each genotype adjusted mean of the corresponding year .
the rationale behind the latter approach is the assumption that the correction for the year effect is better represented by the simple mean of the complete sample of genotypes per year than by just a few checks.
the resulting year effect-corrected genotype means  are forwarded to the gp stage, and through cv are evaluated as predictors.
as in the transition from the first to the second stage, there is a loss of information in passing on from the second to the third stage because the hsv effect is confounded with the residual error term.
this loss can be minimized by weighting the adjusted means  <cit> .
we used the smith et al.
scheme  <cit> , where adjusted means are weighted by the diagonal elements of the inverse of their variance-covariance matrix computed in the first stage.
at this stage, we computed the heritability for each year using the ad hoc method described in piepho and möhring  <cit>  as
<dig>
where  is the genetic variance and  is the mean variance of a difference of two adjusted genotype means, corresponding to the best linear unbiased estimators .
even though this is not the best method to estimate heritability  <cit> , the square root of this heritability estimate gives a rough idea of an upper limit for the predictive abilities.
approach 2: across years analysis
the model to account for the year effect in the second stage through the shared check was
<dig>
where  represents the adjusted mean of grain dry matter yield of the h-th genotype, testcrossed with the v-th tester, in the s-th location and r-th year, gh is the main effect of the h-th genotype, ls is the main effect of the s-th location and drv the main effect of the v-th tester within the r-th year, which can be extended as drv=ar+rv, with ar the effect of the year and t denoting the tester  <cit> .
hrv, hs and rsv are the two-way interaction effects, hrsv is the effect of the three-way interaction and ehrsv is the residual error associated to , with  the variance of the hrsv-th adjusted mean  obtained in the first stage.
the effects containing drv can be extended as hrv=hr+hrv, rsv=rs+rsv and hrsv=hrs+hrsv.
we considered genotypes and testers as fixed factors and location and year as random factors  and .
all effects involving ar are random except rv because we do not want to recover inter-year information since there are only two years and the year by tester classification is very disconnected .
moreover, the rv term is analogous to a block factor in an incomplete block design because it is free of gh; therefore, due to the unbalancedness and the small number of years, we can use it as a fixed effect.
furthermore, the main year effect  can be dropped considering that the adjustment of the genotype means is the same for ar+rv as for only rv.
including all the effects, the final model  is
to minimise the loss of information in the transition to the gp stage, we weighted the adjusted means using the inverse of the squared standard errors, which is also appropriate since we are not fitting random block effects  <cit> .
third stage: genomic prediction
at the third stage, the dataset of p markers was merged with the n grain dry matter yield adjusted means by years of evaluated models.
gp was performed using ridge-regression best linear unbiased prediction , where the genotypic values are predicted using the marker information by regressing each snp on the phenotype  <cit> .
the model was
<dig>
where, m is the n× <dig> vector of phenotypic records, here, containing the adjusted means calculated from the second stage, x and β are, respectively, the design matrix and parameter vector of fixed effects, z is the n×p marker matrix, whose elements zhm represent the snp genotype of the m-th marker of the h-th genotype entry and take the values - <dig>   <dig>  or + <dig> for the aa, aa, and aa genotypes  <cit> , u is the p-dimensional vector of snp effects and e is the error vector.
the term zu is interpreted as the genetic effect and its estimate  as the gebv.
the gebv of the h-th genotype corresponds to , with m= <dig> ⋯,p the number of markers,  is the estimated effect of the m-th marker and zhm the snp genotype of the m-th marker for the h-th genotype entry.
the assumptions of the model are that the error is normally distributed with zero mean and variance r and that u has a normal distribution with zero mean and variance .
r is a diagonal matrix with diagonal elements equal to the inverses of the diagonal elements of the inverse of the original variance-covariance matrix of the adjusted means of the second stage  <cit> .
ip is the p-dimensional identity matrix and  represents the proportion of the genetic variance contributed by each individual snp.
under the model equation  the variance of the observed data is , in which Γ=zzt and zt denotes the transpose of z
<cit> .
to speed up the computation, Γ was rescaled by replacing z with , with p the number of markers  <cit> .
in the year-wise analysis , the genotype adjusted means by year are merged in the m vector, and vector β contains the intercept and the year effect.
in the across-years analysis , where year effect was already accounted for, m contains the genotype adjusted means and vector β contains only the intercept.
in the year-wise analysis correcting genotype adjusted means for year effects , the model used did not include a fixed year factor  but a common intercept, thus the model was the same as for across-years analysis.
to measure the influence of the relationship among the genotypes on the predictions, we used the adjusted means obtained in the second stage and the pedigree information of the entries in a mixed model testing genotypes and crosses as random effects, so that the variances of both effects would give us an estimation of how much the variation is attributed to the pedigree, e.g.
the crosses.
the model was
<dig>
where  is the adjusted mean of the h-th genotype obtained in the second stage, gh is the effect of the h-th genotype, ca is the effect of the a-th grand parent  cross, e.g.
× , and eah the associated error.
additionally, we plotted the relationship heat-map of estimated coefficients of relatedness for individuals based on marker data computed according to wimmer et al.
<cit> .
cross validation for model comparison
to evaluate model performance, k-fold cv was carried out.
in cv, the data is split into k subsets t times.
k- <dig> subsets are used as the training set  and the one other subset is the validation set .
the ts is used to estimate the parameters that then are used to predict the observations in the vs. the performance of the model was assessed by the pearson correlation coefficient between the predicted gebv and the corresponding observations of the vs. this correlation is referred to as predictive ability  <cit> .
as in the first stage, the predictive ability was not adjusted by the square root of the heritability.
although breeding programs are most of the time operating with closely related genotypes, breeders are also interested in knowing the results in a scenario with more distantly related genotypes, for example, using genotypes that share the same grandparents either in the ts or in the vs but not in both.
hence, we wanted to check if accounting for the effect of population structure in the randomisation of cv would make the spatial error models improve the predictive abilities.
we chose two scenarios given the relatedness level of the entries and followed the suggested sampling schemes from albrecht et al.
<cit> , which takes into consideration this fact in the cv procedure.
in the first sampling scheme, hereafter called “within crosses” , random sampling is done using all genotypes in the dataset; in the second scheme, hereafter referred to as “across crosses” , genotypes were clustered by cross, so that complete cross-groups were used randomly either in the vs or the ts.
there were  <dig> crosses of different sizes, sharing none, one or two grand parents.
the general overview of the methodology is depicted in figure  <dig>
model selection
two strategies for selecting the best phenotypic model were used in the first stage.
in strategy one the best model for all locations is selected, that is, there is no model selection per location but across locations.
in strategy two, model selection is location-specific .
for both strategies we computed the aic and performed genomic prediction-cross validation , both per location-year combination.
to accomplish the gp-cv approach, we used the adjusted means per location and year of all spatial and non-spatial models.
then, means of genotypes by year-location combination were joined with the molecular marker data to perform gp-cv, in which genetic values were regressed on markers and validation of the model was done using k-fold cv.
predictions of unobserved records and predictive abilities of each model were obtained for each year-location combination.
we assessed the predictive ability of the models using the pearson correlation coefficient  between the predicted gebv and the observed phenotypic value.
hereafter we denote this predictive ability as ρ-gp-cv.
predictive abilities were not adjusted with the square root of the heritability, as suggested by dekkers  <cit> , since this adds an extra error due to heritability computation  <cit> .figure  <dig>
general representation of strategies to compare model selection methods.
factors were genotype , tester , trial , replicate  and block .
grain dry matter yield  is the response variable in the first stage.
y=g·t:s/r/b is the shorthand notation for the model y
hijkv=hv+s
i+r
ij+b
ijk+e
hijkv.
datasets of  <dig> spatial and non spatial models plus one mixed dataset  generated from best models given the akaike information criterion  and another mixed dataset  generated from best models given the predictive abilities .
for strategy one , the number of locations with the best fits  was counted, so that the model with the best fits in the majority of locations was identified as the best model.
for strategy two , two datasets were built: “mix 1”, containing the adjusted means of the locations with the best fit according to the aic and “mix 2”, containing the adjusted means of the locations with the highest ρ-gp-cv.
thus, after the first stage we had in total eleven data sets of adjusted means, nine corresponding to each tested model from strategy one, plus two more datasets from strategy two: a mixed data set  with the best models per location-year according to the aic, and another mixed set  with best models per location-year according to the ρ-gp-cv.
softwares
all analyses were performed using sas.
stage  <dig> and  <dig> used the mixed procedure and stage  <dig> used proc hpmixed.
relationship matrix was calculated using the synbreed package  <cit>  for r  <dig> .
RESULTS
first stage - strategy 1: model selection across locations
in the first stage - strategy  <dig>  we did model selection across locations using aic and predictive abilities  per location-year combination.
according to the aic, the results favoured the two-dimensional models .
to do a fair comparison between selection methods using aic and ρ-gp-cv, we first describe aic for years  <dig> and  <dig>  for which ρ-gp-cv were also available and then, as additional information, for year  <dig>  for which ρ-gp-cv was not available since the marker information was missing.table  <dig>
akaike information criterion  of models at first stage  by year and location  for grain dry matter yield
0
0
0
0
0
0
0
0
0
0
0
0
0
- <dig>
0
0
0
0
0
0
0
0
0
12
55%
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
0
6
43%
table shows Δaic relative to the best model.
boldfaced entries in the table indicate best model  within location.
empty cells  correspond to locations where the model did not converge.
in italics, we report the models that converged but the hessian matrix was not positive definite.
for years  <dig> and  <dig>  m <dig> and m <dig> had the majority of best fits across locations.
m <dig>  × ar + nugget) resulted in  <dig> out of  <dig> cases as the best model.
m <dig>  × ar + nugget) was best in  <dig> out of  <dig> cases.
the baseline model + row + column  fitted the best 9% of the times and m <dig> 5% of the times.
a similar tendency was observed in  <dig>  where 43% of the times  m <dig> had the best fit and m <dig> was best 29% of the times.
for this year  <dig>  models m <dig>  m <dig> and m <dig> could not be fitted in some locations.
another third of the times , m <dig> had best fits.
interestingly, m <dig> had the best fits in the locations that had convergence problems for models m <dig> and m <dig>  m <dig>  m <dig>  m <dig>  m <dig> and m <dig> never had best fits in any of both groups of years.
the predictive abilities  per location-year combination showed a rather different pattern for best models within locations; however, the two-dimensional models were also more frequently selected than one-dimensional models .
m <dig>  × ar + nugget) showed in seven of  <dig> settings the highest ρ-gp-cv per location-year combination followed by m <dig>  × ar + nugget) with six out of  <dig> times.
the baseline model + row + column  was selected twice and models m <dig>  m <dig> and m <dig> had also one, three and three selections out of  <dig>  respectively.
m <dig>  m <dig> and m <dig> had no best fits at all.table  <dig>
predictive abilities of observed and predicted values of a
5
-fold-cv by year-location combination of models at first stage  for grain dry matter yield , and repeatability  of the trait by location
r
<dig>
<dig>
<dig>
§
<dig>
§
<dig>
<dig>
<dig>
<dig>
<dig>
§
<dig>
<dig>
<dig>
<dig>
<dig>
§
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
<dig>
7
6
boldfaced entries in the table indicate best model  within location.
empty cells correspond to locations where the model did not converge.
in italics, we report the models that converged but the hessian matrix was not positive definite.
§ better than second best model at forth decimal place
one location of  <dig>  produced a negative predictive ability for all models.
we did not consider this location in the counting of best fits, since a higher negative number is actually a worse fit in regard to predictions, but low or high negative are both interpreted as zero prediction.
despite the negative correlations, this location was included in the mixed datasets produced from the site-specific model selection.
we used the adjusted means produced from the baseline model.
another location  showed way lower predictive abilities than the rest of the locations.
to understand these two situations, we calculated the repeatability of the trait in each location for the baseline model.
the repeatability r is defined as the ratio of the between-individual component to the total phenotypic variance  <cit> , which in our case, and following the methodology described by nakagawa and schielzeth  <cit> , corresponds to
<dig>
where  is the between-groups variance and corresponds to the variance of the effect hv fitted as random effect, and in the denominator, the total phenotypic variance given by the sum of the between-groups variance  and the within-groups variances, i.e.
replicates within trials  and blocks within replicates  plus the residual variance .
the interpretation of this repeatability strictly refers to the expected within-group correlations among measurements, i.e.
the agreement among measurements; thus, the gist of the definition of repeatability is related to the reproducibility of the absolute values of measurements.
a slightly higher repeatability in cycle1- <dig> was observed for location g-l <dig> , which involved more trials, i.e.
more genotypes, in comparison with other locations in germany.
the trend in cycle1- <dig> was in favour of the polish locations, which overal l had more homogeneous and higher repeatabilities.
we discuss the relation between repeatabilities and predictive abilities in the next section.
second stage: fitting genotypes by year vs. across years
from a methodological point of view, fitting the year effect in the gp stage was easier and more direct than accounting for the year effect in the second stage, in the sense that the model for the latter approach became too complex and the variance covariance matrix of adjusted means was not possible to be produced using the procedure hpmixed of sas given the high computer power required.
instead, we computed the adjusted means with corresponding standard errors, which were then used to do the weighting to pass on from the second to the third stage.
the adjusted means obtained from the across-years analysis  were plotted against the year effect-corrected genotype adjusted means  to compare the difference of adjustments, in the former case based on one single check against the adjustment given the simple mean of the genotypes in each year .
below the two principal lines, an observation corresponding to the shared check across years stood out from the others, reflecting the year adjustment.
at first glance, it is clear that the check was the only observation pulled down implying that the year adjustment of this check was not strong enough to pull down the observations of the whole year.
both approaches were examined later using the predictive abilities obtained in the gp stage.figure  <dig>
comparison of approaches for year adjustment.
in the x-axis, the genotype adjusted means across-year analysis are plotted.
in the y-axis, the year-effect-corrected adjusted means from the year-wise analysis are depicted.
third stage: genomic prediction
the predictive abilities of the gp stage were taken as the definitive decision criterion for identifying the best strategy for model selection, the best model, and the most reliable approach to account for year effects, and to identify the consequences of population stratification in gp.
we start by presenting results of the comparison of the approaches used for fitting the year effect, since with these we only used the baseline model.
then we present the differences between sampling methods for cv together with the comparison of the models and the model selection strategies.
comparison of approaches to account for year effect in gp
the gp-cv for the approach using the year as a fixed term in the third stage  yielded a predictive ability of  <dig>  , whereas predictive ability for the approach accounting for a fixed year effect in the second stage  was  <dig> .
the predictive ability reached  <dig> , using the year-effect-corrected adjusted means in the gp-cv .
the scatter plots of gebv () against the observed phenotypic values  in the three cases are depicted in figure  <dig>  in approach 1a, we plotted the gebv against the corrected observed phenotypic values, calculated as , where m is the vector of genotype adjusted means obtained in the second stage and  the predicted year effect .
for approach  <dig>  the observed phenotypic values m against  are shown .
for approach 1b, m against  are plotted, with m the year-effect-corrected adjusted means of genotypes .table  <dig>
predictive abilities between observed and predicted values for  <dig> spatial and non-spatial models  and mixed datasets using the best locations given the aic  and the
ρ
-gp-cv per location-year combination
same letters within rows indicate no significant differences  according to a paired t-test.
sampling strategies were: within crosses  and across crosses .
comparison between approaches to fit the year effect.
the y-axis represents the genotype adjusted means  in , m
in  and m
in  ] and the x-axis represents the gebv ().
year-wise analysis , fitting year as fixed effect in the gp stage,  across-years analysis , using year in the second stage and  year-wise analysis using the year effect-corrected genotype means .
ρ
gp represents the predictive ability.
comparison of model selection strategies using different sampling methods in cross validation
fitting model  to measure the influence of the relationship among genotypes on predictions yielded variance components for genotypes, crosses and error for year  <dig> of  <dig> ,  <dig>  and  <dig> , respectively, and for year  <dig> of  <dig> ,  <dig>  and  <dig> , respectively.
thus, the cross effect in  <dig> is contributing in about 40% and in the next year more than 60% to the total variation explained by the data.
the marker-based relationship heat-map  shows some clusters among genotypes of the same cross indicating genetic relatedness.
the predictive abilities using five times 5-fold cv of datasets resulting from first stage analysis of all spatial and non-spatial models plus the mixed datasets were in general very similar within sampling strategies .
for the across-crosses  sampling scheme, the predictive abilities were lower than the ones obtained with the within-crosses  sampling scheme.
in the ac sampling, we fixed the initial seed of the random number generator used for randomization in the cv procedure at the same value for all models to be able to compare the models when the same crosses were used in the training set.figure  <dig>
marker-based relationship heat-map.
visualised are pairwise relationship coefficients estimated from the maker data for genotypes of years  <dig> and  <dig>  higher values represent a stronger relationship.
we compared the models and the sampling methods using a paired t-test  by resembling a randomized complete block design, where the predictive ability of each repetition of the cv was taken as a block, thus accounting for the dependence among observations from the same samples .
for the first sampling method , three groups were identified with some overlaps, but showing not much of a difference among models.
from the across-crosses sampling strategy , five groups were distinguished with some overlaps: m <dig> had the highest predictive ability and models m <dig>  m <dig> and m <dig> had the worst predictive abilities.
potential bias of gp is another important element that could be used to compare models.
we computed the bias as suggested by  <cit> .
the comparison of the biases of all models followed a rather similar trend as the predictive abilities showed in table  <dig>  we present the analysis of bias as supplementary material .
the heritability  for the baseline model was estimated as  <dig>   for year  <dig>   <dig>   for year  <dig> and  <dig>   for  <dig> using the equation .
in principle, the ad hoc method may approximate the true value of heritability but making the unrealistic assumption of uncorrelated genotypes  <cit> .
we computed the heritability to have a rough idea of how much could we expect from the predictive abilities.
the predictive ability divided by square root of heritability is an estimate of the accuracy of gp  <cit> , and the square root of the heritability provides the upper bound for the predictive ability  <cit> , thus one expects that the predictive abilities are not very far from the square root of heritability.
in this case, the square roots of the heritabilities are somewhat larger than the corresponding predictive abilities, indicating that the predictions are not sufficiently accurate due to limited data size, thus not exhausting completely the genetic variance.
to explore in which extend could have our models explained the variance not captured by the markers, we fitted an additional component accounting for the polygenic effect in the gp stage  <cit> .
the baseline model  yielded a genotypic variance of  <dig> ; when we incorporated the polygenic effect, the genotypic variance was  <dig>  and polygenic variance was  <dig> , indicating that about 88% of the total genetic variance was captured by the rr-blup model.
discussion
selecting the models at the first stage produced different results than assessing them in the third stage.
aic had better scores for the models that used row and column effects, e.g.
models m <dig>  m <dig> and m <dig>  or m <dig> that had a 2-dimensional variance-covariance error structure.
ρ-gp-cv also picked m <dig> and m <dig>  but the choices were more spread over the models covering even the baseline model.
in general, in the first stage, both aic and ρ-gp-cv produced better scores for the two-dimensional models, whereas in the third stage the baseline and one-dimensional models seemed to be better than the more complex models .
the explanation of this pattern may be related to the second stage, where the interaction genotype × location played a role.
the two-dimensional models performed very well in modelling heterogeneity within field, but when the means were integrated across the whole experiment, including all locations and years, the two-dimensional spatial error models seemed to over-adjust the means, yielding a poorer predictive ability in the gp stage.
the one-dimensional spatial error models and the two-dimensional model without spatial error structure were sufficient to estimate appropriately adjusted means.
this corroborates piepho and williams  <cit>  who concluded that for small portions of a field, a particular spatial model may hold well but if fitted all across the field it may fail.
in a wheat experiment, lado et al.
<cit>  found that using moving averages as covariable significantly improved the predictive abilities of gp.
they recognised strong heterogeneous patterns of irrigation in the field, that were not controlled with a single blocking system.
models m <dig>  m <dig> and m <dig> were never selected as having the best fits either by aic or ρ-gp-cv.
these models had in common that none of them used rows and columns as additional factors, strengthening the conclusion that row-column designs may have the potential to correctly control field heterogeneity and thus enhance predictive ability of genomic prediction.
fitting a location-specific error model did not have an advantage over fitting a common model across locations.
neither did the dataset composed of means computed using models have best aic fits  nor the second dataset containing the means computed using models with highest ρ-gp-cv  produce better predictive abilities in the gp stage.
the models with nugget had better fits than the corresponding baseline model without the nugget.
the drawback was that fitting those models was not straightforward, since almost every location required a separate coding specifying initial values and lower boundary constraints on the covariance parameters.
good statistical and biological reasons have been presented of why including a nugget to analysis of field experiment is beneficial  <cit> .
if we ignore the two-dimensional spatial models , the aic privileges m <dig> and ρ-gp-cv yields more diverse results with the majority of choices for m <dig> and m <dig>  in fact, when the spatial component of a resolvable row-column design based on linear variance  does not lead to an improved fit, returning to classical row-column design provides randomisation protection  <cit> .
williams and luckett  <cit>  performed studies aiming to find the optimal plot size, the optimal plot arrangements and the best spatial model  and showed that in cotton and barley row and column designs are well suited for variety testing in plant breeding trials.
moreover, recent simulation studies from möhring et al.
<cit>  showed that designs including rows and columns outperformed one-dimensional blocking.
in the same work, the authors mention that blocking in the direction of plots with common long sides is preferable, which is common in cereal breeding  <cit> .
we cannot affirm that ρ-gp-cv was better than aic for model selection or vice versa, nor that the results showed the same trend; but if we would have used either of these two strategies to select the best model, we would have selected the m <dig> with aic or m <dig> with ρ-gp-cv.
the gp predictive ability obtained by m <dig>  was slightly better than m <dig> and m <dig> ; however, this model  was not highlighted by either of the two selection criteria .
in practice, the fact that there were no large statistical differences is good news for the breeders because the baseline model , or even better, the simplest model with row-column adjustment , are appropriate for phenotypic analysis towards gp.
as a model selection method, gp-cv is of interest because it may allow to compare models with different fixed effects, even when reml is used for estimating the variance parameters.
no simple recommendation has been reported concerning the best model selection criterion in the case of spatial models  <cit> .
predictive abilities have been used between environments as similarity measure and then to join similar environments into clusters  <cit> .
thus, in a sense ρ-gp-cv allows giving an interpretation to the environment under scrutiny and the displayed trend do not depart far from the classical aic.
the repeatabilities  presented in parallel to the ρ-gp-cv  show a low correlation  with the predictive abilities from the baseline model.
in fact, we expected that for location p-l <dig> of  <dig>  which had a negative predictability, the r was very low almost zero, but this was not the case; hence we could not conclude that the low predictive ability is mainly due to environmental effects.
riedelsheimer et al.
<cit>  also reported negative predictive accuracies when testing unrelated crosses in the cv procedure and observed that using unrelated crosses could have provided a negative prediction signal due to opposite linkage phases with important qtl displayed in the ts, suggesting that the negative predictive accuracies are associated with the marker pattern.
in this study we explored three ways to adjust the year effect given the weak connectivity across years.
using the single check  to make the year adjustment was not a better choice than adjusting by the simple year mean  or accounting for the year effect in the gp stage , even though the estimated predictive ability was the highest.
the “year clouds” produced using approach  <dig>  did not overlap perfectly, from which we concluded that the correction was not appropriate and generated an over-fitting of the markers in the gp-cv procedure due to the fact that markers also predicted the year effect and not the snp-effects alone.
using the year-mean correction for adjusted means in the second stage  produced a lower ρ-gp-cv, that, given the overlay of the clouds of predicted vs. observed values, seems to be more realistic.
however, fitting the year effect manually, i.e.
using ordinary least squares estimation  vs. fitting it as a fixed effect in the gp stage, i.e.
using generalised least squares estimation  can definitively yield a more precise estimate.
indeed, the residual variance in approach 1b using year effect-corrected adjusted means was around  <dig>   and in approach 1a using the year fixed effect in the gp stage yielded residual variance of  <dig>  .
in approach 1a, where we fitted the year in the gp stage, we removed the year effect from the observed adjusted means derived from the second stage  to avoid bias of the predictive abilities; however, there would still be some bias because the subtracted year effect was not the true effect but an estimate of the year effect.
models were eventually assessed and compared using the ρ-gp-cv in the third stage.
the two sampling scenarios to perform the cv procedure aimed to recreate the cases where the material was genetically close, with some individuals coming from the same parental cross, and more distantly related to avoid individuals from the same parental cross in the randomisation procedure of cv.
this more distantly related material shows some identical-by-state  similarity, therefore it was not unrelated in the theoretical sense of population genetics.
this more distantly related scenario may be seen also as a case where one tries to predict a scenario whose linking information is weak or lacking, e.g.
different genotypes and/or locations in the ts and vs .
the predictive abilities obtained for gp using wc sampling were located in the middle-high range and using ac sampling, predictive abilities were placed in the middle range.
the predictive ability of the ac sampling was significantly lower than wc, as expected for gp of a dataset showing population structure.
riedelsheimer et al.
<cit>  drew similar conclusions using unrelated biparental maize families.
they concluded that predictive accuracy could be increased by adding crosses  sharing both parents to the ts.
in this respect, the use of pedigree and marker information to borrow information from both sources is suggested  <cit> .
CONCLUSIONS
the main conclusions of this study are:  fitting a traditional model including row and column factors across all locations was good enough to account for field heterogeneity in the first stage under gp frame.
this also suggests that row-column designs may be preferable to designs with a single blocking factor;  aic and ρ-gp-cv did not have the same trend in selecting across models, but both favoured in the end models m <dig> and m9; however, none of the methods picked the model with highest predictive ability.
fitting a location-specific error model did not produce an advantage over fitting a common model across locations;  the baseline model  and the simplest row-column adjustment  had in overall the best results, which is very good news since in routine analysis complex models may require much programming expertise and powerful computers;  in a dataset weakly connected across years, a more reasonable model-wise structure is to account for the year factor in the genomic prediction stage rather than in a previous stage, to ensure that the effect is not confounded with the markers adjustment, and  datasets of distantly related genotypes may have a poor performance for gp purposes; however, increasing the size of the crosses may be an opportunity to enhance predictive ability in these cases of disconnected datasets on related sets of genotypes.
electronic supplementary material
additional file 1:
sas codes  used to implement first stage of phenotypic analysis referred in table
4
.
additional file 2:
analysis of bias of genomic prediction.
abbreviations
gpgenomic prediction
gebvgenomic estimated breeding value
snpsingle nucleotide polymorphism
metmulti-environment trial
cvcross validation
aicakaike information criterion
bicbayesian information criterion
remlrestricted maximum likelihood
mlmaximum likelihood
mafminor allele frequency
ρ-gp-cvpredictive ability
bluebest linear unbiased estimator
rr-blupridge regression-best linear unbiased predictor
vsvalidation set
tstraining set
wcwithin-crosses
acacross-crosses
olseordinary least squares estimation
glsegeneralised least squares estimation
ibsidentical-by-state.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
ambv participated in the design of the study, conducting the analysis, interpreting the results, writing and editing the manuscript.
jm participated in the design of the study, conducting the analysis and drafting the manuscript.
mschm supervised the collection of the kws-lochow data set and implemented the models for the whole dataset.
mschö and ccs were responsible for data controlling and participated in the conception of the study.
hpp conceived the study, participated in its design, writing and editing the manuscript and oversaw the project.
all authors read and approved the final manuscript.

 genomic prediction was first introduced in as a method that allows the prediction of genomic estimated breeding values for plants and animals by using information of genetic markers. in plant breeding, gp has been adopted as another stage of the breeding scheme, not diminishing the importance of the phenotypic analysis usually carried out in several environments. merging the phenotype and the genotype analyses has been addressed through the so-called stage-wise analysis. in the first stage environments are analysed separately and genotype means are computed and then submitted in the gp stage to predict gebv based on dense genetic markers such as single nucleotide polymorphisms. in plant breeding, assessing genotypic adaptability and stability, and predicting breeding values of the genotypes in other environments and other years, makes use of multi-environment trials, which aim to evaluate as many genotypes as possible in as many as possible locations. these mets are typically laid out as generalised lattice designs testing a large number of different genotypes per trial. the number of tested genotypes is limited by factors such as seed production, production cycle length and availability of physical resources, e.g. land and budget. within years, genotypes are tested in series of trials, which are connected by checks. checks are lines grown in every trial as controls because their performance is known and/or they are already commercial material. checks can be also used to connect years. in the rye breeding program considered in this paper, a completely different set of genotypes is tested in each year, but these genotypes are from the same breeding population. the accuracy of a genomic prediction model depends on the number of genotypes used for calibration. so there is definitely a need to combine data across years. low connectivity across years is a challenge when trying to combine data across years, and this is one main motivation for this paper. furthermore, the unbalancedness due to the design layout and the different and large number of evaluated genotypes increases the heterogeneity introducing high complexity to the variance-covariance structure among adjusted genotype means. analysis of mets could be done as single-stage analysis, modelling the complete observed data at the level of individual plots, or using a stage-wise approach, where experiments are analysed first at the level of environments, obtaining adjusted means per genotype, which are then summarised across environments in the next stage. a single-stage analysis accounts entirely for the variance-covariance structure of the recorded observations, therefore it is regarded as the gold standard. however, it has been shown that in a stage-wise analysis, a loss of information occurring in the transition through stages can be minimized by an appropriate weighting scheme. if feasible, a single-stage approach is preferable to a stage-wise analysis. nevertheless, the latter is acceptable for gp, since it is simple, computationally more efficient and also allows to easily account for any specifics of randomisation layout and error modelling for each environment. it should be stressed, however, that in a stage-wise analysis the weights are chosen to approximate the variance-covariance matrix of adjusted means from previous stages. we used here a three-stage approach and compared different spatial correlation structures in the first stage to correct field heterogeneity at the trial level. spatial error models may provide more accurate estimates of genotype effects than models not accounting for spatial adjustment but they are computationally more demanding and convergence may be difficult to reach. any effort in terms of improving the genomic predictions would include checking if these improved estimates have an effect on the predictive ability when markers are added to the model. the performance of alternative spatial models can be assessed by k-fold cross validation. similarly, the merits of different spatial models used to compute adjusted means in the first stage can be compared by the same cv procedure, if the same gp procedure is used for each analysis. this suggests that genomic prediction-cross validation can be used to identify the best-fitting mixed model in stage one. the common method of model selection makes use of information criteria based on the log likelihood, e.g. the akaike information criterion or the bayesian information criterion. when the restricted maximum likelihood method is used, models can only be compared by information criteria if they have the same fixed effects; otherwise, the maximum likelihood method should be used. cv is, in this sense, not used to tune parameters as in many penalization methods, machine learning methods) but only as a tool to compare models that use reml. reml is considered the best available method of variance parameter estimation, preferable to ml. consequently, it is of interest to devise model selection procedures that can use reml and also can compare models with different fixed effects. gp-cv has already been used to judge environments in order to optimise the accuracy in gp. we used this tool here as model selection method in comparison to the traditional use of aic. the aims of this work were: i) to assess the advantage for the predictive ability when using a spatial model for phenotypic analysis, ii) to compare stage-wise approaches for gp when the data are weakly connected across years, and iii) to compare aic and gp-cv as methods of selection of models for phenotypic data analysis towards gp in rye. methods field layout and data set a commercial rye breeding program by kws-lochow established in poland and germany aims to develop superior hybrid varieties for the seed market. the implementation of gp within the breeding program makes use of the measurements of hybrid performance of the first cycles of phenotypic evaluation of the material. selections made in cycle are intensively evaluated in further cycles, aiming to double-check the selection decisions. for our purposes, these additional cycles do not add much useful information. hence, we used only the first cycles of the program. the populations tested in each year consist of s genotypes, which display genetic relatedness and population stratification due to complex genealogical history. besides the phenotypic data, a 16k infinium iselect hd custom beadchip was used to characterise individuals from cycle1- and cycle1- and checks. several traits were evaluated during this project: grain dry matter yield, plant height and thousand kernel weight, as well as ordinal scores of rust, mildew and lodging among others. in this work we used grain dry matter yield measurements of the phases of selection cycle1- cycle1- and cycle1- and marker information for the genotypes of and although no marker information of year was available, it makes sense to use this dataset to observe the trend in one additional year and in this way, support the results of the phenotypic analysis of previous years. a cycle experiment consists of subsets of genotypes from the s populations tested in several locations within each of the two countries involving two testers. we define a trial as the physical unit within a location, where a subset of genotypes that were testcrossed to the same tester is evaluated. trials at a location were laid out as α-designs with two replicates. each trial was randomized independently from the others using the software cycdesign.. in our notation, trials of a cycle experiment are labelled as s s …, s row and column coordinates of the plots to account for spatial variation are available.table general representation of the testers by locations by years classification of cycle year and in germany and poland series of trials are represented with the labels s s ⋯ s general representation of the testers by locations classification of cycle year in germany and poland series of trials are represented with the labels s s ⋯ s normally throughout the program, only a single tester was used per location and year, but in some locations, some subsets of genotypes were testcrossed with the two available testers. this is the case, for example, for location g-l in cycle1- where the genotypes evaluated in the trials s s and s were testcrossed with both tester and tester and it is also the case of locations p-l p-l p-l and p-l evaluating genotypes of trials s s and s with both testers. in each year, four common checks were testcrossed with the testers and grown twice in each trial. over the years and one check was in common and none was shared with.table year x check classification in germany and poland the field layout of some trials was not perfectly rectangular. some trials at a given location and year had fewer blocks but larger size, i.e., there were two different block sizes within a few trials. blocks were nested within rows of the field layout. in the genetic dataset, homozygous marker genotypes were coded as - and and the heterozygous type, missing values and technical failures were coded as % of the markers corresponded to homozygous alleles and % were heterozygous. only a % of the markers were recorded as missing values or technical failures; therefore, an imputation method would not have a strong impact on the subsequent analyses. monomorphic markers and markers with minor allele frequency less than 1% or missing information of more than 10% per marker were dropped. a total of markers passed the quality test and were used for gp. models in this section we present the models used in the first stage of the analysis and the models of the approaches followed to adjust the year effect either in the second or the third stage. figures and depict a general scheme that helps visualizing the methodology.figure general representation of stage-wise approaches to compare year-effect adjustment. factors were genotype, tester, location, year, trial, replicate and block. grain dry matter yield is the response variable in the first stage, m is the adjusted mean of genotypes across locations used in the second stage, m is the year effect-corrected genotype adjusted mean, represents the simple mean of genotypes of the r-th year. in the genomic prediction stage, m is the n× vector of adjusted means of genotypes by year for approach 1a and across years for approach m is the n× vector of adjusted means of year effect-corrected genotypes in approach 1b, x and β are respectively the design matrix and parameter vector of fixed effects, z is the n×p marker matrix, u is the p-dimensional vector of snp effects and e the error vector. y=g·t:s/r/b is the shorthand notation of the model eq. in the text: y hijkv=hv+s i+r ij+b ijk+e hijkv, m =g×l×t stands for the model eq. in the text:, and m =×g×l represents the extended model eq. in the text:. the final predictive abilities are presented in the ellipses. general representation of model comparison through all the stages of the analysis. datasets generated from spatial and non-spatial models plus two mixed datasets generated from best models given the akaike information criterion and the predictive abilities. factors in second stage were genotype, location and tester. m represents the adjusted mean of genotypes across locations and years. m =g×l×t is the shorthand notation for. in the genomic prediction stage m is the adjusted mean of genotypes across locations, x and β are respectively the design matrix and parameter vector of fixed effects, z is the n×p marker matrix, u is the p-dimensional vector of snp effects and e the error vector. sampling methods in cross validation were across crosses and within crosses. the final predictive abilities are presented in the ellipses. first stage in the first stage we computed adjusted genotype means by location and year. the factors used for the analysis were genotypes, testers, trials, replicates nested within trials and blocks nested within replicates. we defined a baseline model as where yhijkv is the observed grain dry matter yield of the h-th genotype testcrossed with the v-th tester in the k-th block within the j-th replicate of the i-th trial, hv is the effect of the h-th genotype testcrossed with the v-th tester, si is the effect of the i-th trial, rij is the effect of the j-th replicate nested within the i-th trial, bijk is the effect of the k-th block nested within the j-th replicate of the i-th trial and ehijkv is the plot error associated with the yhijkv observation. in model equation we assumed genotypes crossed with testers as a fixed effect to be able to compute genotype adjusted means per tester, whereas the other effects were considered as random effects due to the nested design structure. spatial and non-spatial models used for the first stage y hijkv=hv+s i+r ij+b ijk+e hijkv y hijkoqv=hv+s i+r ij+b ijk ijo+v ijq+e hijkoqv y hijkv=hv+s i+r ij+b ijk+e hijkv y hijkv=hv+s i+r ij+b ijk+e hijkv y hijkoqv=hv+s i+r ij+b ijk ijo+v ijq+e hijkoqv y hijkv=hv+s i+r ij+b ijk+e hijkv y hijkv=hv+s i+r ij+b ijk+e hijkv y hijkv=hv+s i+r ij+b ijk+e hijkv y hijkoqv=hv+s i+r ij+b ijk ijo+v ijq+e hijkoqv y hijkv is the observed dry matter yield of the h-th genotype testcrossed with the v-th tester in the k-th block within the j-th replicate of the i-th trial, hv is the effect of the h-th genotype testcrossed with the v-th tester, s i is the effect of the i-th trial, r ij is the effect of the j-th replicate nested within the i-th trial, b ijk is the effect of the k-th block nested within the j-th replicate of the i-th trial and e hijkv is the plot error associated with the y hijkv observation. in the models including row and column effects, w ijo is the effect of the o-th row within the j-th replicate of the i-th trial and v ijq is the effect of the q-th column within the j-th replicate of the i-th trial. spatial variance-covariance structure were independent, autoregressive in one direction, one-dimension linear variance and two-dimension autoregressive. note that we use hv as fixed effect, which is necessary to obtain the genotype by tester means. the purpose is also to recover the information of the entries that are grown in the same locations but using different testers, so that we captured the effect of the tester in the shared locations. second stage in the second stage we computed genotype means across locations and testers. this was done either separately for each year or also averaging across years. the years and where molecular marker data were available, were connected through only one check. the resulting fundamental question is then how to fit the year effect. either the year effect is estimated by the mean of all tested entries or we rely on the adjustment by the one single check. we assume that genotypes tested in each year can be regarded as a random sample from the same parent population. based on the structure of the breeding program, this is a realistic assumption that motivates the approaches described in the following. both approaches were compared using the m resulting from the analysis of the baseline model in the first stage. approach 1: year-wise analysis each year was analysed in the second stage using a three-way interaction model of genotypes, locations and testers as factors to obtain adjusted genotypes means of each year. the model was where represents the adjusted mean of grain dry matter yield of the h-th genotype, testcrossed with the v-th tester in the s-th location, gh, ls and tv are the main effects of the h-th genotype, the s-th location and the v-th tester, respectively, hs, hv and sv are the two-way interaction effects, hsv is the effect of the three-way interaction and ehsv is the residual error associated with, with the variance of the hsv-th adjusted mean obtained in the first stage. location was considered as random effect and hence, all the interactions containing this factor are random. the crossed effect of genotypes and testers could have been a fixed effect since genotypes and testers are taken as fixed factors in this stage. however, the crossed effects that include g were taken as random here because the factor genotype was used as random in the gp stage. but note that in the first and the second stage we needed to take genotype main effects as fixed in order to compute adjusted means. besides, since not every genotype was tested with every tester, we needed to take hv random to be able to estimate genotype means across levels of testers. in this approach, the year effect was adjusted in two ways, hereafter referred as to approach 1a and approach 1b. approach 1a used years as fixed factors in the gp stage and approach 1b used a manual adjustment after the second stage by simply calculating the mean of the genotypes by year and subtracting it to each genotype adjusted mean of the corresponding year. the rationale behind the latter approach is the assumption that the correction for the year effect is better represented by the simple mean of the complete sample of genotypes per year than by just a few checks. the resulting year effect-corrected genotype means are forwarded to the gp stage, and through cv are evaluated as predictors. as in the transition from the first to the second stage, there is a loss of information in passing on from the second to the third stage because the hsv effect is confounded with the residual error term. this loss can be minimized by weighting the adjusted means. we used the smith et al. scheme, where adjusted means are weighted by the diagonal elements of the inverse of their variance-covariance matrix computed in the first stage. at this stage, we computed the heritability for each year using the ad hoc method described in piepho and möhring as where is the genetic variance and is the mean variance of a difference of two adjusted genotype means, corresponding to the best linear unbiased estimators. even though this is not the best method to estimate heritability, the square root of this heritability estimate gives a rough idea of an upper limit for the predictive abilities. approach 2: across years analysis the model to account for the year effect in the second stage through the shared check was where represents the adjusted mean of grain dry matter yield of the h-th genotype, testcrossed with the v-th tester, in the s-th location and r-th year, gh is the main effect of the h-th genotype, ls is the main effect of the s-th location and drv the main effect of the v-th tester within the r-th year, which can be extended as drv=ar+rv, with ar the effect of the year and t denoting the tester. hrv, hs and rsv are the two-way interaction effects, hrsv is the effect of the three-way interaction and ehrsv is the residual error associated to, with the variance of the hrsv-th adjusted mean obtained in the first stage. the effects containing drv can be extended as hrv=hr+hrv, rsv=rs+rsv and hrsv=hrs+hrsv. we considered genotypes and testers as fixed factors and location and year as random factors and. all effects involving ar are random except rv because we do not want to recover inter-year information since there are only two years and the year by tester classification is very disconnected. moreover, the rv term is analogous to a block factor in an incomplete block design because it is free of gh; therefore, due to the unbalancedness and the small number of years, we can use it as a fixed effect. furthermore, the main year effect can be dropped considering that the adjustment of the genotype means is the same for ar+rv as for only rv. including all the effects, the final model is to minimise the loss of information in the transition to the gp stage, we weighted the adjusted means using the inverse of the squared standard errors, which is also appropriate since we are not fitting random block effects. third stage: genomic prediction at the third stage, the dataset of p markers was merged with the n grain dry matter yield adjusted means by years of evaluated models. gp was performed using ridge-regression best linear unbiased prediction, where the genotypic values are predicted using the marker information by regressing each snp on the phenotype. the model was where, m is the n× vector of phenotypic records, here, containing the adjusted means calculated from the second stage, x and β are, respectively, the design matrix and parameter vector of fixed effects, z is the n×p marker matrix, whose elements zhm represent the snp genotype of the m-th marker of the h-th genotype entry and take the values - or + for the aa, aa, and aa genotypes, u is the p-dimensional vector of snp effects and e is the error vector. the term zu is interpreted as the genetic effect and its estimate as the gebv. the gebv of the h-th genotype corresponds to, with m= ⋯,p the number of markers, is the estimated effect of the m-th marker and zhm the snp genotype of the m-th marker for the h-th genotype entry. the assumptions of the model are that the error is normally distributed with zero mean and variance r and that u has a normal distribution with zero mean and variance. r is a diagonal matrix with diagonal elements equal to the inverses of the diagonal elements of the inverse of the original variance-covariance matrix of the adjusted means of the second stage. ip is the p-dimensional identity matrix and represents the proportion of the genetic variance contributed by each individual snp. under the model equation the variance of the observed data is, in which Γ=zzt and zt denotes the transpose of z. to speed up the computation, Γ was rescaled by replacing z with, with p the number of markers. in the year-wise analysis, the genotype adjusted means by year are merged in the m vector, and vector β contains the intercept and the year effect. in the across-years analysis, where year effect was already accounted for, m contains the genotype adjusted means and vector β contains only the intercept. in the year-wise analysis correcting genotype adjusted means for year effects, the model used did not include a fixed year factor but a common intercept, thus the model was the same as for across-years analysis. to measure the influence of the relationship among the genotypes on the predictions, we used the adjusted means obtained in the second stage and the pedigree information of the entries in a mixed model testing genotypes and crosses as random effects, so that the variances of both effects would give us an estimation of how much the variation is attributed to the pedigree, e.g. the crosses. the model was where is the adjusted mean of the h-th genotype obtained in the second stage, gh is the effect of the h-th genotype, ca is the effect of the a-th grand parent cross, e.g. ×, and eah the associated error. additionally, we plotted the relationship heat-map of estimated coefficients of relatedness for individuals based on marker data computed according to wimmer et al.. cross validation for model comparison to evaluate model performance, k-fold cv was carried out. in cv, the data is split into k subsets t times. k- subsets are used as the training set and the one other subset is the validation set. the ts is used to estimate the parameters that then are used to predict the observations in the vs. the performance of the model was assessed by the pearson correlation coefficient between the predicted gebv and the corresponding observations of the vs. this correlation is referred to as predictive ability. as in the first stage, the predictive ability was not adjusted by the square root of the heritability. although breeding programs are most of the time operating with closely related genotypes, breeders are also interested in knowing the results in a scenario with more distantly related genotypes, for example, using genotypes that share the same grandparents either in the ts or in the vs but not in both. hence, we wanted to check if accounting for the effect of population structure in the randomisation of cv would make the spatial error models improve the predictive abilities. we chose two scenarios given the relatedness level of the entries and followed the suggested sampling schemes from albrecht et al., which takes into consideration this fact in the cv procedure. in the first sampling scheme, hereafter called “within crosses”, random sampling is done using all genotypes in the dataset; in the second scheme, hereafter referred to as “across crosses”, genotypes were clustered by cross, so that complete cross-groups were used randomly either in the vs or the ts. there were crosses of different sizes, sharing none, one or two grand parents. the general overview of the methodology is depicted in figure model selection two strategies for selecting the best phenotypic model were used in the first stage. in strategy one the best model for all locations is selected, that is, there is no model selection per location but across locations. in strategy two, model selection is location-specific. for both strategies we computed the aic and performed genomic prediction-cross validation, both per location-year combination. to accomplish the gp-cv approach, we used the adjusted means per location and year of all spatial and non-spatial models. then, means of genotypes by year-location combination were joined with the molecular marker data to perform gp-cv, in which genetic values were regressed on markers and validation of the model was done using k-fold cv. predictions of unobserved records and predictive abilities of each model were obtained for each year-location combination. we assessed the predictive ability of the models using the pearson correlation coefficient between the predicted gebv and the observed phenotypic value. hereafter we denote this predictive ability as ρ-gp-cv. predictive abilities were not adjusted with the square root of the heritability, as suggested by dekkers, since this adds an extra error due to heritability computation.figure general representation of strategies to compare model selection methods. factors were genotype, tester, trial, replicate and block. grain dry matter yield is the response variable in the first stage. y=g·t:s/r/b is the shorthand notation for the model y hijkv=hv+s i+r ij+b ijk+e hijkv. datasets of spatial and non spatial models plus one mixed dataset generated from best models given the akaike information criterion and another mixed dataset generated from best models given the predictive abilities. for strategy one, the number of locations with the best fits was counted, so that the model with the best fits in the majority of locations was identified as the best model. for strategy two, two datasets were built: “mix 1”, containing the adjusted means of the locations with the best fit according to the aic and “mix 2”, containing the adjusted means of the locations with the highest ρ-gp-cv. thus, after the first stage we had in total eleven data sets of adjusted means, nine corresponding to each tested model from strategy one, plus two more datasets from strategy two: a mixed data set with the best models per location-year according to the aic, and another mixed set with best models per location-year according to the ρ-gp-cv. softwares all analyses were performed using sas. stage and used the mixed procedure and stage used proc hpmixed. relationship matrix was calculated using the synbreed package for r.  first stage - strategy 1: model selection across locations in the first stage - strategy we did model selection across locations using aic and predictive abilities per location-year combination. according to the aic, the results favoured the two-dimensional models. to do a fair comparison between selection methods using aic and ρ-gp-cv, we first describe aic for years and for which ρ-gp-cv were also available and then, as additional information, for year for which ρ-gp-cv was not available since the marker information was missing.table akaike information criterion of models at first stage by year and location for grain dry matter yield 0 0 0 0 0 0 0 0 0 0 0 0 0 - 0 0 0 0 0 0 0 0 0 12 55% 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 6 43% table shows Δaic relative to the best model. boldfaced entries in the table indicate best model within location. empty cells correspond to locations where the model did not converge. in italics, we report the models that converged but the hessian matrix was not positive definite. for years and m and m had the majority of best fits across locations. m × ar + nugget) resulted in out of cases as the best model. m × ar + nugget) was best in out of cases. the baseline model + row + column fitted the best 9% of the times and m 5% of the times. a similar tendency was observed in where 43% of the times m had the best fit and m was best 29% of the times. for this year models m m and m could not be fitted in some locations. another third of the times, m had best fits. interestingly, m had the best fits in the locations that had convergence problems for models m and m m m m m and m never had best fits in any of both groups of years. the predictive abilities per location-year combination showed a rather different pattern for best models within locations; however, the two-dimensional models were also more frequently selected than one-dimensional models. m × ar + nugget) showed in seven of settings the highest ρ-gp-cv per location-year combination followed by m × ar + nugget) with six out of times. the baseline model + row + column was selected twice and models m m and m had also one, three and three selections out of respectively. m m and m had no best fits at all.table predictive abilities of observed and predicted values of a 5 -fold-cv by year-location combination of models at first stage for grain dry matter yield, and repeatability of the trait by location r § § § § 7 6 boldfaced entries in the table indicate best model within location. empty cells correspond to locations where the model did not converge. in italics, we report the models that converged but the hessian matrix was not positive definite. § better than second best model at forth decimal place one location of produced a negative predictive ability for all models. we did not consider this location in the counting of best fits, since a higher negative number is actually a worse fit in regard to predictions, but low or high negative are both interpreted as zero prediction. despite the negative correlations, this location was included in the mixed datasets produced from the site-specific model selection. we used the adjusted means produced from the baseline model. another location showed way lower predictive abilities than the rest of the locations. to understand these two situations, we calculated the repeatability of the trait in each location for the baseline model. the repeatability r is defined as the ratio of the between-individual component to the total phenotypic variance, which in our case, and following the methodology described by nakagawa and schielzeth, corresponds to where is the between-groups variance and corresponds to the variance of the effect hv fitted as random effect, and in the denominator, the total phenotypic variance given by the sum of the between-groups variance and the within-groups variances, i.e. replicates within trials and blocks within replicates plus the residual variance. the interpretation of this repeatability strictly refers to the expected within-group correlations among measurements, i.e. the agreement among measurements; thus, the gist of the definition of repeatability is related to the reproducibility of the absolute values of measurements. a slightly higher repeatability in cycle1- was observed for location g-l, which involved more trials, i.e. more genotypes, in comparison with other locations in germany. the trend in cycle1- was in favour of the polish locations, which overal l had more homogeneous and higher repeatabilities. we discuss the relation between repeatabilities and predictive abilities in the next section. second stage: fitting genotypes by year vs. across years from a methodological point of view, fitting the year effect in the gp stage was easier and more direct than accounting for the year effect in the second stage, in the sense that the model for the latter approach became too complex and the variance covariance matrix of adjusted means was not possible to be produced using the procedure hpmixed of sas given the high computer power required. instead, we computed the adjusted means with corresponding standard errors, which were then used to do the weighting to pass on from the second to the third stage. the adjusted means obtained from the across-years analysis were plotted against the year effect-corrected genotype adjusted means to compare the difference of adjustments, in the former case based on one single check against the adjustment given the simple mean of the genotypes in each year. below the two principal lines, an observation corresponding to the shared check across years stood out from the others, reflecting the year adjustment. at first glance, it is clear that the check was the only observation pulled down implying that the year adjustment of this check was not strong enough to pull down the observations of the whole year. both approaches were examined later using the predictive abilities obtained in the gp stage.figure comparison of approaches for year adjustment. in the x-axis, the genotype adjusted means across-year analysis are plotted. in the y-axis, the year-effect-corrected adjusted means from the year-wise analysis are depicted. third stage: genomic prediction the predictive abilities of the gp stage were taken as the definitive decision criterion for identifying the best strategy for model selection, the best model, and the most reliable approach to account for year effects, and to identify the consequences of population stratification in gp. we start by presenting results of the comparison of the approaches used for fitting the year effect, since with these we only used the baseline model. then we present the differences between sampling methods for cv together with the comparison of the models and the model selection strategies. comparison of approaches to account for year effect in gp the gp-cv for the approach using the year as a fixed term in the third stage yielded a predictive ability of, whereas predictive ability for the approach accounting for a fixed year effect in the second stage was. the predictive ability reached, using the year-effect-corrected adjusted means in the gp-cv. the scatter plots of gebv () against the observed phenotypic values in the three cases are depicted in figure in approach 1a, we plotted the gebv against the corrected observed phenotypic values, calculated as, where m is the vector of genotype adjusted means obtained in the second stage and the predicted year effect. for approach the observed phenotypic values m against are shown. for approach 1b, m against are plotted, with m the year-effect-corrected adjusted means of genotypes.table predictive abilities between observed and predicted values for spatial and non-spatial models and mixed datasets using the best locations given the aic and the ρ -gp-cv per location-year combination same letters within rows indicate no significant differences according to a paired t-test. sampling strategies were: within crosses and across crosses. comparison between approaches to fit the year effect. the y-axis represents the genotype adjusted means in, m in and m in ] and the x-axis represents the gebv (). year-wise analysis, fitting year as fixed effect in the gp stage, across-years analysis, using year in the second stage and year-wise analysis using the year effect-corrected genotype means. ρ gp represents the predictive ability. comparison of model selection strategies using different sampling methods in cross validation fitting model to measure the influence of the relationship among genotypes on predictions yielded variance components for genotypes, crosses and error for year of, and, respectively, and for year of, and, respectively. thus, the cross effect in is contributing in about 40% and in the next year more than 60% to the total variation explained by the data. the marker-based relationship heat-map shows some clusters among genotypes of the same cross indicating genetic relatedness. the predictive abilities using five times 5-fold cv of datasets resulting from first stage analysis of all spatial and non-spatial models plus the mixed datasets were in general very similar within sampling strategies. for the across-crosses sampling scheme, the predictive abilities were lower than the ones obtained with the within-crosses sampling scheme. in the ac sampling, we fixed the initial seed of the random number generator used for randomization in the cv procedure at the same value for all models to be able to compare the models when the same crosses were used in the training set.figure marker-based relationship heat-map. visualised are pairwise relationship coefficients estimated from the maker data for genotypes of years and higher values represent a stronger relationship. we compared the models and the sampling methods using a paired t-test by resembling a randomized complete block design, where the predictive ability of each repetition of the cv was taken as a block, thus accounting for the dependence among observations from the same samples. for the first sampling method, three groups were identified with some overlaps, but showing not much of a difference among models. from the across-crosses sampling strategy, five groups were distinguished with some overlaps: m had the highest predictive ability and models m m and m had the worst predictive abilities. potential bias of gp is another important element that could be used to compare models. we computed the bias as suggested by. the comparison of the biases of all models followed a rather similar trend as the predictive abilities showed in table we present the analysis of bias as supplementary material. the heritability for the baseline model was estimated as for year for year and for using the equation. in principle, the ad hoc method may approximate the true value of heritability but making the unrealistic assumption of uncorrelated genotypes. we computed the heritability to have a rough idea of how much could we expect from the predictive abilities. the predictive ability divided by square root of heritability is an estimate of the accuracy of gp, and the square root of the heritability provides the upper bound for the predictive ability, thus one expects that the predictive abilities are not very far from the square root of heritability. in this case, the square roots of the heritabilities are somewhat larger than the corresponding predictive abilities, indicating that the predictions are not sufficiently accurate due to limited data size, thus not exhausting completely the genetic variance. to explore in which extend could have our models explained the variance not captured by the markers, we fitted an additional component accounting for the polygenic effect in the gp stage. the baseline model yielded a genotypic variance of; when we incorporated the polygenic effect, the genotypic variance was and polygenic variance was, indicating that about 88% of the total genetic variance was captured by the rr-blup model. discussion selecting the models at the first stage produced different results than assessing them in the third stage. aic had better scores for the models that used row and column effects, e.g. models m m and m or m that had a 2-dimensional variance-covariance error structure. ρ-gp-cv also picked m and m but the choices were more spread over the models covering even the baseline model. in general, in the first stage, both aic and ρ-gp-cv produced better scores for the two-dimensional models, whereas in the third stage the baseline and one-dimensional models seemed to be better than the more complex models. the explanation of this pattern may be related to the second stage, where the interaction genotype × location played a role. the two-dimensional models performed very well in modelling heterogeneity within field, but when the means were integrated across the whole experiment, including all locations and years, the two-dimensional spatial error models seemed to over-adjust the means, yielding a poorer predictive ability in the gp stage. the one-dimensional spatial error models and the two-dimensional model without spatial error structure were sufficient to estimate appropriately adjusted means. this corroborates piepho and williams who concluded that for small portions of a field, a particular spatial model may hold well but if fitted all across the field it may fail. in a wheat experiment, lado et al. found that using moving averages as covariable significantly improved the predictive abilities of gp. they recognised strong heterogeneous patterns of irrigation in the field, that were not controlled with a single blocking system. models m m and m were never selected as having the best fits either by aic or ρ-gp-cv. these models had in common that none of them used rows and columns as additional factors, strengthening the conclusion that row-column designs may have the potential to correctly control field heterogeneity and thus enhance predictive ability of genomic prediction. fitting a location-specific error model did not have an advantage over fitting a common model across locations. neither did the dataset composed of means computed using models have best aic fits nor the second dataset containing the means computed using models with highest ρ-gp-cv produce better predictive abilities in the gp stage. the models with nugget had better fits than the corresponding baseline model without the nugget. the drawback was that fitting those models was not straightforward, since almost every location required a separate coding specifying initial values and lower boundary constraints on the covariance parameters. good statistical and biological reasons have been presented of why including a nugget to analysis of field experiment is beneficial. if we ignore the two-dimensional spatial models, the aic privileges m and ρ-gp-cv yields more diverse results with the majority of choices for m and m in fact, when the spatial component of a resolvable row-column design based on linear variance does not lead to an improved fit, returning to classical row-column design provides randomisation protection. williams and luckett performed studies aiming to find the optimal plot size, the optimal plot arrangements and the best spatial model and showed that in cotton and barley row and column designs are well suited for variety testing in plant breeding trials. moreover, recent simulation studies from möhring et al. showed that designs including rows and columns outperformed one-dimensional blocking. in the same work, the authors mention that blocking in the direction of plots with common long sides is preferable, which is common in cereal breeding. we cannot affirm that ρ-gp-cv was better than aic for model selection or vice versa, nor that the results showed the same trend; but if we would have used either of these two strategies to select the best model, we would have selected the m with aic or m with ρ-gp-cv. the gp predictive ability obtained by m was slightly better than m and m; however, this model was not highlighted by either of the two selection criteria. in practice, the fact that there were no large statistical differences is good news for the breeders because the baseline model, or even better, the simplest model with row-column adjustment, are appropriate for phenotypic analysis towards gp. as a model selection method, gp-cv is of interest because it may allow to compare models with different fixed effects, even when reml is used for estimating the variance parameters. no simple recommendation has been reported concerning the best model selection criterion in the case of spatial models. predictive abilities have been used between environments as similarity measure and then to join similar environments into clusters. thus, in a sense ρ-gp-cv allows giving an interpretation to the environment under scrutiny and the displayed trend do not depart far from the classical aic. the repeatabilities presented in parallel to the ρ-gp-cv show a low correlation with the predictive abilities from the baseline model. in fact, we expected that for location p-l of which had a negative predictability, the r was very low almost zero, but this was not the case; hence we could not conclude that the low predictive ability is mainly due to environmental effects. riedelsheimer et al. also reported negative predictive accuracies when testing unrelated crosses in the cv procedure and observed that using unrelated crosses could have provided a negative prediction signal due to opposite linkage phases with important qtl displayed in the ts, suggesting that the negative predictive accuracies are associated with the marker pattern. in this study we explored three ways to adjust the year effect given the weak connectivity across years. using the single check to make the year adjustment was not a better choice than adjusting by the simple year mean or accounting for the year effect in the gp stage, even though the estimated predictive ability was the highest. the “year clouds” produced using approach did not overlap perfectly, from which we concluded that the correction was not appropriate and generated an over-fitting of the markers in the gp-cv procedure due to the fact that markers also predicted the year effect and not the snp-effects alone. using the year-mean correction for adjusted means in the second stage produced a lower ρ-gp-cv, that, given the overlay of the clouds of predicted vs. observed values, seems to be more realistic. however, fitting the year effect manually, i.e. using ordinary least squares estimation vs. fitting it as a fixed effect in the gp stage, i.e. using generalised least squares estimation can definitively yield a more precise estimate. indeed, the residual variance in approach 1b using year effect-corrected adjusted means was around and in approach 1a using the year fixed effect in the gp stage yielded residual variance of. in approach 1a, where we fitted the year in the gp stage, we removed the year effect from the observed adjusted means derived from the second stage to avoid bias of the predictive abilities; however, there would still be some bias because the subtracted year effect was not the true effect but an estimate of the year effect. models were eventually assessed and compared using the ρ-gp-cv in the third stage. the two sampling scenarios to perform the cv procedure aimed to recreate the cases where the material was genetically close, with some individuals coming from the same parental cross, and more distantly related to avoid individuals from the same parental cross in the randomisation procedure of cv. this more distantly related material shows some identical-by-state similarity, therefore it was not unrelated in the theoretical sense of population genetics. this more distantly related scenario may be seen also as a case where one tries to predict a scenario whose linking information is weak or lacking, e.g. different genotypes and/or locations in the ts and vs. the predictive abilities obtained for gp using wc sampling were located in the middle-high range and using ac sampling, predictive abilities were placed in the middle range. the predictive ability of the ac sampling was significantly lower than wc, as expected for gp of a dataset showing population structure. riedelsheimer et al. drew similar conclusions using unrelated biparental maize families. they concluded that predictive accuracy could be increased by adding crosses sharing both parents to the ts. in this respect, the use of pedigree and marker information to borrow information from both sources is suggested.  the main conclusions of this study are: fitting a traditional model including row and column factors across all locations was good enough to account for field heterogeneity in the first stage under gp frame. this also suggests that row-column designs may be preferable to designs with a single blocking factor; aic and ρ-gp-cv did not have the same trend in selecting across models, but both favoured in the end models m and m9; however, none of the methods picked the model with highest predictive ability. fitting a location-specific error model did not produce an advantage over fitting a common model across locations; the baseline model and the simplest row-column adjustment had in overall the best results, which is very good news since in routine analysis complex models may require much programming expertise and powerful computers; in a dataset weakly connected across years, a more reasonable model-wise structure is to account for the year factor in the genomic prediction stage rather than in a previous stage, to ensure that the effect is not confounded with the markers adjustment, and datasets of distantly related genotypes may have a poor performance for gp purposes; however, increasing the size of the crosses may be an opportunity to enhance predictive ability in these cases of disconnected datasets on related sets of genotypes. electronic supplementary material additional file 1: sas codes used to implement first stage of phenotypic analysis referred in table 4. additional file 2: analysis of bias of genomic prediction. abbreviations gpgenomic prediction gebvgenomic estimated breeding value snpsingle nucleotide polymorphism metmulti-environment trial cvcross validation aicakaike information criterion bicbayesian information criterion remlrestricted maximum likelihood mlmaximum likelihood mafminor allele frequency ρ-gp-cvpredictive ability bluebest linear unbiased estimator rr-blupridge regression-best linear unbiased predictor vsvalidation set tstraining set wcwithin-crosses acacross-crosses olseordinary least squares estimation glsegeneralised least squares estimation ibsidentical-by-state. competing interests the authors declare that they have no competing interests. authors’ contributions ambv participated in the design of the study, conducting the analysis, interpreting the results, writing and editing the manuscript. jm participated in the design of the study, conducting the analysis and drafting the manuscript. mschm supervised the collection of the kws-lochow data set and implemented the models for the whole dataset. mschö and ccs were responsible for data controlling and participated in the conception of the study. hpp conceived the study, participated in its design, writing and editing the manuscript and oversaw the project. all authors read and approved the final manuscript. 
in the case of met poorly connected between years, analysing each year separately and fitting year as a fixed effect in the genomic prediction stage yielded the most realistic predictive abilities. we compared stage-wise approaches to analyse a real dataset of a multi-environment trial in rye, which was connected between years only through one check, and used different spatial models to obtain better estimates, and thus, improved predictive abilities for genomic prediction. the aims of this study were to assess the advantage of using spatial models for the predictive abilities of genomic prediction, to identify suitable procedures to analyse a met weakly connected across years using different stage-wise approaches, and to explore genomic prediction as a tool for selection of models for phenotypic data analysis. using complex spatial models did not significantly improve the predictive ability of genomic prediction, but using row and column effects yielded the highest predictive abilities of all models. predictive abilities can also be used to select models for phenotypic data analysis. the trend of the predictive abilities was not the same as the traditionally used akaike information criterion, but favoured in the end the same models. it makes use of genotypic information to make predictions used for selection decisions. genomic prediction is becoming a daily tool for plant breeders. a proper phenotypic analysis is a crucial prerequisite for accurate calibration of genomic prediction procedures. stage-wise analysisgenomic predictioncross validationspatial modelsmulti-environment trials restricted maximum likelihood issue-copyright-statement© the author 2014
 genomic prediction is becoming a daily tool for plant breeders. it makes use of genotypic information to make predictions used for selection decisions. the accuracy of the predictions depends on the number of genotypes used in the calibration; hence, there is a need of combining data across years. a proper phenotypic analysis is a crucial prerequisite for accurate calibration of genomic prediction procedures. we compared stage-wise approaches to analyse a real dataset of a multi-environment trial in rye, which was connected between years only through one check, and used different spatial models to obtain better estimates, and thus, improved predictive abilities for genomic prediction. the aims of this study were to assess the advantage of using spatial models for the predictive abilities of genomic prediction, to identify suitable procedures to analyse a met weakly connected across years using different stage-wise approaches, and to explore genomic prediction as a tool for selection of models for phenotypic data analysis.  using complex spatial models did not significantly improve the predictive ability of genomic prediction, but using row and column effects yielded the highest predictive abilities of all models. in the case of met poorly connected between years, analysing each year separately and fitting year as a fixed effect in the genomic prediction stage yielded the most realistic predictive abilities. predictive abilities can also be used to select models for phenotypic data analysis. the trend of the predictive abilities was not the same as the traditionally used akaike information criterion, but favoured in the end the same models.  making predictions using weakly linked datasets is of utmost interest for plant breeders. we provide an example with suggestions on how to handle such cases. rather than relying on checks we show how to use year means across all entries for integrating data across years. it is further shown that fitting of row and column effects captures most of the heterogeneity in the field trials analysed. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords stage-wise analysisgenomic predictioncross validationspatial modelsmulti-environment trials restricted maximum likelihood issue-copyright-statement© the author 2014 