We used an experimental data set of winter barley (Osthushenrich et al., 2018) for obtaining estimates of genomic marker effects, genotypic values, and variance components. In our simulations, we assumed that these estimates were the true values of the respective parameters.

The experimental data set consisted of 250 doubled-haploid (DH) lines which were derived from 25 crosses of 10 parental lines. Five of the parental lines were elite lines and five were resistant donors. The derived lines and the parental lines were evaluated for yield at five locations in two years. The field trial was layed out as an augmented design with five blocks. In each block, 50 of the derived lines were tested together with the 10 parental lines that served as checks. In the first year, each of the derived lines was tested unreplicated in four locations, in the second year each of the derived lines was grown in two replications at five locations.

To estimate the adjusted treatment means of the of the parental lines and the derived lines we used the linear model

where l is the effect of the line, e is the effect of the environment, l:e is the genotype-by-environment interaction, r:e is the replication within environment effect, b:r:e is the block effect nested within replication and environment, and ε is the residual. The genotype was analyzed as a fixed factor, the remaining factors of the model were random. The adjusted treatment means ranged from 69.4 to 99.3 dt/ha (adjusted to 15% moisture), the standard errors of the means were 2.7 for the 10 check genotypes in the augmented design, and 3.1 for the 250 genotypes with one replication per block. The heritability for unbalanced trials, according to Equation (19) of Piepho and Möhring (2007), was h² = 0.83.

To estimate the variance components from the derived lines we used the linear model

where c is the effect of a cross, l:c is the effect of the line within the cross, e is the environmental effect, c:e the cross-by-environment interaction, l:c:e the line-by-environment interaction, r:e is the replication within environment effect, b:r:e is the block effect nested within replication and environment, and ε is the residual. Assuming all effects as random, the cross-by-environment variance σ^ce2=10.3, the line-by-environment variance σ^le2=11.1, and the residual variance σ^ε2=30.1 were estimated with software ASReml-R (Butler et al., 2017).

The lines were genotyped with a 50 K single-nucleotide polymorphism (SNP) chip (Trait Genetics, Gatersleben). SNPs with more than two recorded alleles, more than 10% missing values and a gene diversity smaller than 0.1 were excluded from the analysis, as well as genotypes with more than 15% missing information. After preprocessing the marker data, 9,597 SNP markers and 259 genotypes (249 DH lines and 10 parental lines) remained for the analysis.

The pairwise modified Roger's distances (cf Reif et al., 2005) between the lines were used in a principal coordinate analysis and a heatmap to illustrate the degree of relatedness of the lines (Figure 1).

From the adjusted treatment means, we estimated the genome-wide effects for yield with the RR-BLUP method (Meuwissen et al., 2001).

Starting with the marker genotypes of the parental lines, we simulated the crosses and the development of DH lines from the crosses. The simulations were genetic simulations of recombination along chromosomes using the code of the software Plabsim (Maurer et al., 2008). The result of a simulation run are the marker genotypes of a set of simulated DH lines. The genotypic values g of the simulated DH lines were determined from their marker genotypes and the corresponding genomic marker effects. For this purpose, we used the genomic marker effects obtained from the experimental data set and assumed that these were the true values. This procedure results in simulated DH lines of which the marker genotypes as well as the genotypic values g are known. The genotypic value of the l-th line of the c-th cross is denoted by g_cl.

The evaluation of the DH lines in field trials at E environments with R replications per environment was simulated using normally distributed random numbers. The masking effect for the l-th line of the c-th cross was m_cl = u_c+v_l+w_ε where u_c, v_l, and w_ε were realizations of independent random variables with distribution u~N(0,σce2/E), v~N(0,σle2/E), and w~N[0,σε2/(ER)]. The variance components σce2, σle2, and σε2 were determined on basis of the variance components estimated from the experimental barley data as described in the subsequent section on the simulated scenarios. The effects u_c were determined for each cross and the effects v_l for each DH line. With the simulated genotypic values g_cl and the masking effects m_cl, the phenotypic values were determined as p_cl = g_cl+m_cl.

The phenotypic values p_cl of the simulated DH lines were in turn used to estimate simulated genomic marker effects with the RR-BLUP method. The simulated genomic marker effects were then used to calculate genomic predictions y_cl of the true genotypic values g_cl.

In order to illustrate the principle of using genomic prediction for selection candidates evaluated in unreplicated field trials, we simulated two experiments with a mating scheme similar to that of Osthushenrich et al. (2018) as a starting point. These two experiments were simulated with only one simulation run. After a factorial cross of the 5 × 5 parental lines, 10 DH lines were generated from each cross, resulting in 250 DH lines.

We considered two different field trial designs for the evaluation of the 250 simulated DH lines in order to compare unreplicated and replicated trials. The first investigated design was an unreplicated trial with one environment (E = 1) and one plot per genotype (R = 1). The second investigated design was a replicated trial with three environments (E = 3) and two replications per environment (R = 2). The variance components σce2, σle2, and σε2 that were estimated from the experimental data set were used for calculating the masking effects m_cl and thus the simulated phenotypic values p_cl for the two different field trial designs.

In order to investigate the effect of varying the variance components σce2, σle2, and σε2 of the masking variance for unreplicated (E = 1, R = 1) and replicated (E = 3, R = 2) field trials, we ran a set of replicated simulations with 2,000 simulation runs for each parameter setting. In addition to the original factorial crossing scheme of 5 × 5 lines, we also simulated a diallel mating scheme. In the diallel mating scheme, all 10 pairwise crosses between the 5 elite lines were carried out.

In the factorial mating scheme we simulated either family sizes of 10 DH lines per cross, resulting in a population of 250 lines, or family sizes of six lines per cross, resulting in a population size of 150 lines. For the diallel mating scheme, we generated 25 DH lines for each of the 10 crosses, resulting in a population size of 250 lines. Thus, we investigated three different mating schemes, all of which were evaluated in both an unreplicated (E = 1, R = 1) and a replicated (E = 3, R = 2) field trial.

For each of these six scenarios, we considered five different sets of variance components σce2, σle2, and σε2 of the masking variance in order to assess the effect on the correlations r(p, g), r(y, p), and r(y, g).

For the first set of variance components, we used a setting that is close to the variance components of the experimental design: σce2=10,σle2=10,σε2=30. For the second set of variance components, we used a smaller cross-by-environment variance and a greater line-by-environment variance than in the experimental data set: σce2=2,σle2=18,σε2=30. For the third set of variance components, we used a greater cross-by-environment variance and a smaller line-by-environment variance than in the experimental data set: σce2=18,σle2=2,σε2=30. In addition, we investigated a scenario with a greater error variance: σce2=10, σle2=10, σε2=60, and one in which the interaction variances as well as the error variance were considerably greater than in the experimental data set: σce2=20, σle2=20, σε2=120.

In phenotypic selection, the phenotypic values p assessed in field trials are used as estimates for the genotypic values g of the selection candidates. The correlation r(p, g) depends on the heritability and determines the response to selection. The response to genomic selection is determined by the correlation r(y, g) between the genomic predictions y and the genotypic values g. The correlation r(y, g) is referred to as prediction accuracy (Legarra et al., 2008). We used the mean values of the correlations across the replicated simulations to compare the efficiency of phenotypic selection with genomic selection.

The R code used to carry out the simulations is available at Github (https://github.com/JT-Giessen/Terraillon_2022).

Our experimental data was collected in two years, in the first year one plot per genotype was planted in each of four locations (R = 1, E = 4), in the second year two replications were grown in each location (R = 2, E = 5). In the retrospective analysis we regarded the data from each of the four environments in year one as an unreplicated trial (R = 1, E = 1). This results in four sets of phenotypic values. We denote the sets of phenotypic values obtained from these unreplicated trials with p^*.

We then used the four sets of phenotypic values p^* to estimate four sets of genomic effects using the RR-BLUP method. From these we calculated four sets of genomic performance estimates y^*. For each set of values p^* we used the remaining three locations of the first year and the data from the second year to determine the adjusted treatment means of the lines according to Equation (1). These phenotypic values were then regarded as estimators for the genotypic values ĝ^*.

For each of the four sets of values p^*, y^*, and ĝ^* we determined the correlations r(p^*, ĝ^*) between the phenotypic values and the estimate of the genotypic value, and r(y^*, ĝ^*) between the genomic prediction and the estimate of the genotypic value.

A one-sided Pearson z-test with the alternative hypothesis H_A: r(p^*, ĝ^*) < r(y^*, ĝ^*) was carried out with the R package cocor (Diedenhofen and Musch, 2015) to test whether the genomic predictions are superior to the phenotypic values in predicting the estimated genotypic value.

In our study, we use a model of a breeding program for our simulations. Here we discuss briefly the rationales behind choosing its components.

In a field trial, the environment, the genotype-by-environment interaction, the design factors, such as blocks or replications, and the experimental error contribute to the phenotypic value of the tested genotypes. This is modeled in our approach by adding a masking effect to the genotypic value. If, for a given genotype, n replications with respect to a given variance component σ² are available, then averaging over the n replications results in a transformed random variable with variance σ²/n. Summing up realizations of transformed random variables for all variance components of a trial allows to model environmental conditions similar to those of the trial from which the variance components were estimated.

In applied breeding programs, it is a common approach to select the best lines available from one cycle of material development and to recombine these as parents to obtain the selection candidates for the next cycle. If the lines are from one material group, then one option for recombining the lines is to mate every line with every other line, this corresponds to a diallel mating scheme. If the lines originate from two pools, for example, a resistance pool and a pool with high yielding genotypes, then one option to recombine the lines is to mate every line from the first pool to every line of the second pool. This corresponds to a factorial mating scheme. This was the rationale for investigating both, the diallel and the factorial mating scheme in our study.

In experimental evaluations of genomic prediction, the correlation r(y, p) between genomic predictions y and phenotypic values p is often used as a measure for the precision of prediction (e.g., Hofheinz et al., 2012; Albrecht et al., 2014; Lorenz and Smith, 2015; Zenke-Philippi et al., 2017; Werner et al., 2018). The reason for this is that the correlation r(y, g) between the genomic predictions y and the true genotypic values g remains unknown. Some authors divide the correlation r(y, p) by the square root of the heritability (e.g., Albrecht et al., 2011; Zhao et al., 2013; Technow et al., 2014; Sallam et al., 2015). This is a linear transformation, and the ranking of selection candidates remains the same as in the untransformed data assuming constant heritability. The simulation methodology employed in this study provides, in contrast to experimental evaluations, a direct assessment of the correlations r(p, g) and r(y, g) between the phenotypic values p and genotypic values g and between the genomic predictions y and the genotypic values g, respectively. This allows the comparison of the selection criteria p and y with respect to their precision to predict the true unknown genotypic values g of selection candidates. Moreover, it facilitates the investigation of quantitative genetic factors that affect the two correlations.

The genomic marker effects and the components of the masking variance from our experimental barley data set were used as the basis for the simulations. This methodology implies simplifying assumptions, the most important of which we point out here briefly. The purely additive genetic model neglecting epistasis as well as the estimation method for genomic marker effects might have an effect on the results. It is neglected that a certain SNP variant might be in linkage disequilibrium with alleles having different genomic effects in breeding material that originates from genetically different sources. Moreover, we assume that the genomic effects of SNPs were estimated in the original barley data set without a residual error. In particular this simplification might result in an overestimation of the correlations of y and p with the true genotypic values g. In addition, the genetic structure of the parental lines can be assumed to have an effect on the size of the estimated genomic marker effects.

We used genome-wide effects estimated from an experimental data set as the basis for our simulations. An alternative approach is to use genetic effects that were drawn with a random number generator from a probability distribution. A prominent example for this methodology is Meuwissen et al. (2001). We have chosen to use effects estimated from experimental data, because we think that this might be closer to reality than effects from a random number.

In consequence, the results reported here can only illustrate the concept of using genomic predictions of tested material. Further research is needed on the basis of other experimental data sets or quantitative genetic scenarios to confirm the transferability to other situations of the results reported here.

In the simulation using the variance components and genomic marker effects of Osthushenrich et al. (2018), the correlation r(p, g) between the phenotypic values p and the genotypic values g in an unreplicated field trial was 0.53 (Figure 2). The strength of this correlation might just be sufficient for selection, but the selection gain that can be reached is low.

In a field trial in three environments with two replications (E = 3, R = 2), the correlation between the phenotypic and the genotypic values amounted to r(p, g) = 0.82, which should allow for efficient selection with reasonable selection gain. However, the increase in precision in comparison to the unreplicated field trial comes at the cost of a six-fold increase of the resources required for the field trial.

In contrast, the correlation r(y, g) between the genomic predictions and the true genotypic values amounted to r(y, g) = 0.87 in the unreplicated trial (Figure 2). Hence, the marker-based genomic predictions y of the genotypic value have a much higher correlation with the true genotypic value g than the phenotypic values p of the unreplicated trial. The correlation reaches values that are typically reached by replicated trials and even surpasses the correlation r(p, g) reached by the replicated field trial with E = 3, R = 2. This enables an efficient selection of genotypes with high performance which results in a high selection gain.

From an applied point of view, the costs of genotyping a line might equal roughly the costs of one field plot. In consequence, our results suggest for the investigated data set that using one field plot and in addition genotyping the lines can realize a selection gain roughly corresponding to that of a replicated trial, but requires only a third of the resources.

The correlations r(y, g) for unreplicated trials reached values that surpassed the correlations r(p, g) for the corresponding replicated trials for all investigated scenarios (Table 1). For example, an experimental design close to the experimental data set resulted in the correlation r(y, g) = 0.88 for the unreplicated trial (l. 1 of Table 1) and in the correlation r(p, g) = 0.80 for the replicated trial (l. 6 of Table 1). We conclude that using unreplicated trials combined with genomic selection might be a useful strategy for a wide range of genetic situations, which we address subsequently.

The lower limit for the population size that is required to apply the proposed approach of predicting tested lines is not yet reached with 150 individuals. This conclusion is supported by the r(y, g) = 0.86 of l. 11 of Table 1 and r(p, g) = 0.80 presented in l. 16 of Table 1.

In our analysis we divided the genotype-by-environment interaction variance into two components, the cross-by-environment variance σce2 and the line-by-environment variance σle2. With an increasing importance of the cross-by-environment component of σce2=2,10,18 (l. 2, 1, 3 of Table 1), the correlation between genomic predictions and genotypic values r(y, g) decreased while the correlation r(y, p) increased. This result suggests that a large genotype-by-cross variance might render the assessment of prediction accuracy using the correlation r(y, p) difficult as there might be situations in which changing parameters of a genomic selection program increases r(y, p) without actually increasing r(y, g).

With increasing error variances and increasing genotype-by-environment variances (e.g., l. 1, 4, 5 of Table 1), the correlations r(p, g) decreased considerably. The correlations r(y, g) remained above 0.75 nevertheless in the factorial mating schemes. This suggests that even when the correlations between the phenotypic and genotypic values of an unreplicated trial get low, using genomic predictions instead can still enable successful genomic selection of candidate genotypes.

In the diallel mating scheme of the elite lines, the genetic variance was only a fraction of the genetic variance in the complete factorial. This observation might be caused by the contrasting yield of the two parental groups in the factorial mating scheme. The small genetic variance in comparison to the components of the masking variance in the diallel mating scheme results in a low heritability in the unreplicated trial accompanied by a correlation r(p, g) of 0.19 or 0.20 (l. 21–23 of Table 1). Even in this extreme scenario, using the genomic predictions for the genotypic values for selection results in correlations r(y, g) of 0.37 to 0.52.

These comparisons suggest that using genomic predictions y of the genotypic values might be preferable over using the phenotypic values p directly for a wide range of variance components and populations sizes with our data set.

The high correlation between the true genotypic values g and their genomic predictions y observed in this study might be explained by the replicated evaluation of chromosome segments in different lines. Due to the replicated use of the parental lines, the parental chromosome segments are present and evaluated in the field in a large number of derived lines. The amount of this pseudo-replication of chromosome segments strongly depends on the degree of relatedness between the evaluated lines and the crossing scheme that was used to develop them.

From a quantitative genetics point of view, the approach suggested here can be related to older best linear unbiased prediction (BLUP) approaches. Assuming a large number of markers with small genomic effects of similar size, the genomic prediction model converges to a GBLUP model that employs the realized relationship matrix between the tested genotypes. Those type of prediction models showed a good performance in maize (Bernardo, 1994). A direct consequence of these hypotheses on the mechanism behind the observed results is that the approach suggested here can only be expected to work if the material under investigation is related. This, however, is typically the case in breeding programs for cultivar development, because often several selection candidates were derived from the same crosses, and the parents of the crosses are used in several crosses.

A thorough experimental validation of the results obtained from the presented simulation study requires an experimental data set that estimates the genotypic values of the tested lines with high precision. To achieve highly precise estimates, typically trials over several years (e.g., more than tree) at many locations (e.g., more than 10) are required. The precision of the field trial used to collect our experimental data can be considered as sufficient in applied breeding programs, but it does not provide highly precise estimates of the genotypic value. In consequence, this data set can not be used for a rigid experimental validation of our simulation results. Nevertheless, it can be used to demonstrate the potential usefulness of our approach.

In the retrospective re-analysis we estimated the genotypic values ĝ^* from three locations with one plot per genotype in the first year and five locations with two replications per genotype in the second year. The fourth location in the first year was considered as a “preliminary yield trial” (in terms of Endelman et al., 2014) to obtain the phenotypic values p^*. The research question was, whether the estimations of the genotypic values ĝ^* can be better approximated by the field data p^* of single-plot experiments or by genomic predictions y^* obtained from the field data of single-plot experiments.

For all four data sets, which correspond to the four single-plot experiments in the first year, we observed a greater correlation r(y^*, ĝ^*) than r(p^*, ĝ^*). Hence, the genotypic predictions from single-plot experiments approximated the genotypic values better than the original phenotypic values of the experiments.