Analysis of genomic data

All parental inbred lines were genotyped with the Illumina MaizeSNP50 BeadChip (Ganal et al. 2011). We removed all markers missing or heterozygous in >5% of the inbred lines. Remaining missing (0.2%) or heterozygous (0.3%) marker genotypes were replaced with the most frequent allele. A total of 35,478 markers were subsequently available for further analysis. The marker data are provided in File S1, File S2, and File S3.

Overall pairwise linkage disequilibrium (LD) between markers on the same chromosome was computed as r², separately for the Dent and Flint group, using only markers with a minor allele frequency (MAF) ≥ 0.025 in the respective group (24,242 markers for the Dent lines and 23,450 for the Flint lines). To diminish the confounding effect of varying marker density on regional LD patterns along the chromosomes, we reduced the marker density to ~4 markers per megabase (Mb), with a spacing of ~0.25 Mb, resulting in 4958 markers available for analysis within the Dent and 4929 within the Flint heterotic group. Nevertheless, some density differences could not be completely eliminated. This was because in some instances, no segregating markers could be found in the desired intervals. We then divided all chromosomes into bins of 5Mb width and computed the average pairwise LD, measured as r², between all markers in the bin. For each bin, we also determined the proportion of marker pairs with the same linkage phase, i.e., same sign of the r statistic in Dent and Flint (Technow et al. 2012), and the correlation between the r values of both groups. For this, 4397 markers with a MAF ≥ 0.025 in each group were used.

MAF patterns along the chromosomes were investigated using a similar approach. Again we used consecutive bins of 5Mb width and computed the average MAF for each bin in the sets of Dent and Flint lines as well as the average absolute difference between the reference allele frequencies in the two groups. These investigations were carried out using all 35,478 markers. The allele that had highest frequency across the combined set of Dent and Flint lines was defined as the reference allele.

Variance components and adjusted means

We used a two-stage analysis for estimation of variance components and adjusted entry means that closely followed Bernardo (1996) and Massman et al. (2013). Two-stage analysis is commonly used for analyzing plant breeding field trials and delivers in most cases results similar to those of considerably more complex one-stage approaches (Möhring and Piepho 2009). Its main advantage is the strongly reduced computational burden when numbers of genotypes and environments are large.

In the first stage, hybrid × environment means y were calculated with a standard α-lattice design analysis to adjust for the effects of the field design in these environments. In the second stage, we fitted the model

where vector y contained the phenotypic observations of the hybrids in the 131 environments obtained in stage one, β was the vector of fixed effects of environments, and X was the corresponding design matrix.

The design matrices Z_D and Z_F associated the random general combining ability (GCA) effects of the parental Dent lines (g_D) and Flint lines (g_F), respectively, to the observations of the hybrids in y. Z_S was the design matrix of the random specific combining ability (SCA) effects (s) for specific Dent × Flint hybrid combinations in y. The residuals were represented by vector e. The covariance matrix of g_D was $G_{\text{D}}{\sigma }_{\text{D}}^2,$ that of g_F was $G_{\text{F}}{\sigma }_{\text{F}}^2,$ and that of s was $S{\sigma }_{\text{s}}^2,$ where ${\sigma }_{\text{D}}^2,$ ${\sigma }_{\text{F}}^2,$ and ${\sigma }_{\text{s}}^2$ were the variance components pertaining to GCA and SCA effects. The covariance matrix of the residuals was $R{\sigma }_{\text{R}}^2,$ with ${\sigma }_{\text{R}}^2$ being the residual variance. The diagonal elements of R were the reciprocals of the number of replications in the environment of the corresponding data points. All other elements of R were zero. In the two-stage analysis applied in our study, the genotype × environment variance cannot be separated from the residual variance associated with the adjusted means in y (Möhring and Piepho 2009). Variance component ${\sigma }_{\text{R}}^2$ therefore contained the residual as well as the genotype × environment variance. This enabled also a direct comparison with the results of Massman et al. (2013), who used the same approach for computing variance components and entry means.

The genomic relationship matrix G_D was computed according to VanRaden (2008) as $G_{\text{D}}=W_{\text{D}}{{\mathbf{W}}^{\prime }}_{\text{D}}/m_{\text{D}},$ where m_D is the number of markers and $w_{uv}=\left(x_{uv}-2p_v\right)/\sqrt{4p_v\left(1-p_v\right)}$ (u being the index of the inbred line and v that of the marker), with x_uv coding the number of reference alleles, i.e., 0 or 2, and p_v being the allele frequency of the reference allele in the population of Dent lines. The genomic relationship matrix G_F was computed accordingly. For computing G_D and G_F, only markers were used that segregated in the respective heterotic group with MAF ≥ 0.025.

Let D and D* denote any two Dent lines and F and F* any two Flint lines. For a given pair of single crosses (D × F) and (D* × F*), the element of S was the product $g_{D{D}^{\mathrm{\ast }}}{g}_{F{F}^{\mathrm{\ast }}},$ where $g_{D{D}^{\mathrm{\ast }}}$ and $g_{F{F}^{\mathrm{\ast }}}$ are the corresponding elements of G_D and G_F, pertaining to D and D* and F and F*, respectively (Stuber and Cockerham 1966).

The variance components were estimated for the whole data set, using the EM algorithm for restricted maximum likelihood described by Henderson (1985) and adapted for variance component estimation in factorials by Bernardo (1996). The entry-mean heritability was computed as $H^2=\left({\sigma }_{\text{D}}^2+{\sigma }_{\text{F}}^2+{\sigma }_{\text{s}}^2\right)/\left({\sigma }_{\text{D}}^2+{\sigma }_{\text{F}}^2+{\sigma }_{\text{s}}^2+{\sigma }_{\text{R}}^2/e_{\text{H}}\right),$ where e_H was the harmonic mean of the diagonal elements of ${{\mathbf{Z}}^{\prime }}_{\text{s}}{R}^{-1}{Z}_{\text{s}},$ i.e., of the total number of replications per hybrid. Finally, environment-adjusted entry means of all hybrids (y*) were computed as $y^{\text{*}}={\left({{\mathbf{Z}}^{\prime }}_{\text{s}}{R}^{-1}{Z}_{\text{s}}\right)}^{-1}{{\mathbf{Z}}^{\prime }}_{\text{s}}{R}^{-1}\left(y-X\beta \right),$ following Bernardo (1996). The adjusted entry means are provided in File S4.

GBLUP

The performance of untested hybrids was predicted by GBLUP with the formula $C_{\text{UT}}{V}_{\text{TT}}^{-1}{y}_{\text{T}}^{*}$ (Henderson 1973). Here, C_UT is the genetic covariance matrix of untested and tested hybrids, V_TT is the phenotypic covariance matrix of the tested hybrids, and $y_{\text{T}}^{*}$ are the observed phenotypic values of the tested hybrids (a subset of y*). The elements of C_UT and V_TT were computed according to Bernardo (1996), using our estimates of $g_{D{D}^{\mathrm{\ast }}}$ and $g_{F{F}^{\text{*}}}.$

BayesB

Our BayesB-type model for the performance of the ith hybrid corresponded to model S₂ of Technow et al. (2012):

Here, the linear predictor of the performance of the ith hybrid is denoted as μ_i and β₀ is a common intercept. The row vectors $M_{{\text{D}}_i},$ $M_{{\text{F}}_i},$ and D_i are known marker genotype incidence vectors for the additive marker effects of the Dent parent lines in u_D and Flint parent lines in u_F and the dominance effects in d_DF. The likelihood of a single data point was a Gaussian density with mean parameter equal to μ_i and variance ${\sigma }_e^2.$

The elements of the matrices M_D and M_F code the presence or absence of the reference allele in the gametes produced by the parental Dent and Flint lines as 1/2 and −1/2, respectively. In contrast to M_D and M_F, which code the genotypes of parental gametes, matrix D directly reflects the genotypes of the single-cross hybrids, coding heterozygous genotypes as 1 and homozygous genotypes as 0. For example, if the allele contributed by the Dent parent was “C” and that by the Flint parent “T”, and T had the higher allele frequency, then the corresponding elements of M_D, M_F, and D were −1/2, 1/2, and, 1, respectively.

Additive effects were estimated only for markers with a MAF ≥ 0.025 within the set of tested inbred line parents of the respective heterotic group and dominance effects only for markers with a MAF ≥ 0.025 in at least one of the groups. We reduced the marker density to ~10 markers per megabase to facilitate computations. Using higher marker densities did not improve prediction accuracies, as far as we could see. In total, additive effects were estimated for m_D = 7500 markers of the Dent parental lines and for m_F = 6500 markers of Flint parental lines, on average. The average number of markers for which dominance effects were estimated was m_DF = 8900.

Prior specifications as well as the Gibbs-sampling strategy were identical to those in Technow et al. (2012). The same uninformative prior distribution, a Gamma distribution with α = β = 0.1, was used for the scale parameter S². However, the hyperparameters ν and π were set to constant values, as in the original BayesB implementation of Meuwissen et al. (2001). Parameter ν was set to 4.001 for all types of marker effects, and π was chosen such that the number of markers fitted was 500, on average; e.g., for dominance effects (1 − π_DF)m_DF = 500.

Three independent Gibbs-sampling chains were run for 75,000 iterations, of which the first 74,000 iterations were discarded as burn-in. Using a higher number of iterations and chains did not improve prediction accuracy. The posterior means of marker effects were used to predict the performance of untested hybrids according to model (2).

For investigating the genetic architecture, namely the distribution and properties of marker effects, we fitted model (2), using all 1254 hybrids. The marker density was further reduced to ~1 marker per megabase, or 1617 markers used in total. This was done mainly to counter potential problems with likelihood identifiability that can occur when the number of effects is much larger than the sample size (Gianola 2013). All markers used segregated in each set of parental inbred lines with MAF ≥ 0.025. Thus, all three types of marker effects (additive effects for Dent and Flint and dominance effects) were estimated for each marker. For each trait, we ran 24 independent Gibbs-sampling chains for 1,000,000 iterations. We discarded the first 990,000 iterations as burn-in and afterward stored only samples from every 10th iteration. The posterior means of the marker effects were used as their point estimates.

Estimation of population-specific marker effects for genomic prediction of crossbreds and single-cross hybrids

There are many parallels between hybrid breeding in crops like maize and crossbreeding in livestock production. In simulation studies on genomic prediction with training sets consisting of crossbred individuals (Ibánez-Escriche et al. 2009; Zeng et al. 2013) or single-cross hybrids (Technow et al. 2012), it was found that genomic prediction models that fitted specific marker effects for the parental populations (i.e., purebred breed or heterotic group) had little or no advantage over simpler models that assumed marker effects to be the same across parental populations. One explanation the authors gave for this was that the linkage phase consistency across populations was sufficiently high at high marker densities. In addition, the authors argued that the strongly increased dimensionality of those models prevented them from efficiently capturing remaining across-population differences in marker effects. Knowing in which genomic regions marker–QTL linkage phases are consistent or not could also be used for developing models that estimate population-specific marker effects only where necessary. This would reduce the dimensionality of these models and might mitigate some of the problems associated with it.

We also observed that fitting marker effects to be the same across heterotic groups delivered virtually the same prediction accuracy as model (2) in which specific marker effects were estimated for Dent and Flint (results not shown). This seems to contradict our observation that marker effects are consistent only in pericentromeric regions. However, as is discussed later in detail, prediction of hybrid performance is mostly driven by the presence of close relatives in the training set, in particular for T2 and T1 hybrids. As shown by Habier et al. (2007), BayesB can capture such pedigree relationships, particularly when many markers are fitted. Capturing pedigree relationships with markers does not require physical linkage between them and the QTL (Habier et al. 2013). Consistency of marker–QTL linkage phase might therefore not be mandatory for accurate predictions when close relatives are present in the training set.

Hill–Robertson effect and heterosis

It is known that recombination is suppressed in the pericentromeric regions of maize chromosomes (Gore et al. 2009; Schnable et al. 2009; Ganal et al. 2011; Bauer et al. 2013) and while gene density is comparably low in these regions (Schnable et al. 2009), they still contain a considerable portion of genes (Gore et al. 2009). The Hill–Robertson effect (Hill and Robertson 1966; Felsenstein 1974) describes the influence of recombination on selection efficiency. This effect predicts a buildup of repulsion-phase linkage between QTL alleles when recombination is suppressed (McVean and Charlesworth 2000). One consequence of repulsion-phase linkage is pseudo-overdominance, because additive QTL effects cancel out. Based on the Hill–Robertson effect, McMullen et al. (2009) hypothesized that the strongly suppressed recombination in pericentromeric regions of maize results in pseudo-overdominance and is therefore a major cause of heterosis. Larièpe et al. (2012) mapped dominance QTL in Dent × Flint crosses for important agronomic traits, using a North Carolina III design, and found a large proportion of QTL with (pseudo-)overdominance in pericentromeric regions. Schön et al. (2010) observed the same for Stiff-Stalk Synthetic × Non-Stiff-Stalk crosses. They concluded that pseudo-overdominance in pericentromeric regions led to differential fixation of QTL alleles in each heterotic group. We found the largest allele frequency differences between Dent and Flint in pericentromeric regions and therefore conclude that also for the Dent × Flint heterotic pattern differential fixation in pericentromeric regions takes place.

As allele frequencies in opposite heterotic groups drift apart during reciprocal recurrent selection (Labate et al. 1999), the ratio of SCA variance to GCA variance decreases (Reif et al. 2007) and dominance effects are increasingly absorbed into the population mean or become inseparable from additive effects (i.e., when QTL are fixed in one group but still segregate in the other). In particular, QTL with strongly positive (pseudo-)dominance or (pseudo-)overdominance effects are expected to be affected by differential fixation. The dominance effects of these QTL increase the “baseline” heterosis of the Dent × Flint heterotic pattern but are not detectable with statistical means in our set of Dent × Flint interpool hybrids. This can explain why dominance marker effects had positive and negative signs in almost equal proportions even though dominance effects for grain yield in maize are expected to be mostly positive (Schön et al. 2010). It also explains the absence of any noticeable accumulation of major dominance marker effect estimates in pericentromeric regions.

Comparison of prediction methods

GBLUP and BayesB achieved nearly identical prediction accuracies. Both GY and GM are considered to be highly polygenic traits, based on QTL mapping results (Schön et al. 2004; Huang et al. 2010). Several authors found in simulation studies and for real data sets that GBLUP models were superior to or equally well performing as Bayesian whole-genome regression methods for such traits (Zhong et al. 2009; Hayes et al. 2010; Clark et al. 2011; Kärkkäinen and Sillanpää 2012; Technow and Melchinger 2013; Wimmer et al. 2013). Daetwyler et al. (2010) arrived at the same conclusion based on theoretical results. Zhao et al. (2013) compared several methods for genomic prediction of grain yield of wheat hybrids and also found that GBLUP delivers the same or slightly higher prediction accuracy than BayesB. Thus, no substantial differences between both methods are expected for prediction of hybrid performance for traits like GY and GM.

If the effects of single QTL in polygenic traits vary considerably in size, adaptively shrinking Bayesian whole-genome regression methods could potentially outperform GBLUP. Furthermore, we hypothesize that BayesB could have an advantage for prediction of performance of T0 hybrids from lines distantly related to the parents of the training set hybrids, because then prediction accuracy would mainly come from short-range LD, which is not captured optimally by GBLUP (Habier et al. 2013).

Bayesian whole-genome regression methods can suffer from a lack of likelihood identifiability, when the number of markers is much larger than the size of the training set (Gianola 2013). This can lead to computational and convergence problems in Gibbs sampling (Gelfand and Sahu 1999). As reported by Technow and Melchinger (2013), nonidentifiability can impair prediction accuracy. Our BayesB model for prediction of hybrid performance fits up to three effects per marker, thereby exacerbating the problem. Consequently, Bayesian whole-genome regression methods require larger sizes of the training set for realizing a potential advantage.

Technow et al. (2012) confirmed in a simulation study that BayesB can achieve slightly higher prediction accuracy than GBLUP under a polygenic trait architecture, with a training set comprising 800 hybrids. Assembling large training sets is possible even for moderately sized breeding programs, like the one of the University of Hohenheim. With our data set, for example, a training set of 1254 hybrids could have been assembled, albeit without the possibility of performing a thorough cross-validation. Nonetheless, given the considerably greater computational demands of Bayesian whole-genome regression methods, GBLUP seems to be a very pragmatic and robust method for genomic prediction of hybrid performance for polygenic traits.

Prediction accuracy of T2, T1, and T0 hybrids

We confirmed the sizeable differences in prediction accuracy between T2, T1, and T0 hybrids found in the simulation study of Technow et al. (2012). The same was observed by Maenhout et al. (2010) and Schrag et al. (2010), who compared only T1 and T0 hybrids. These differences can be explained by the different numbers of parents of the hybrids that are also parents of training set hybrids (i.e., two, one, and zero for T2, T1, and T0 hybrids, respectively); the more that are shared, the higher the accuracy that can be expected (Technow et al. 2012). The paramount importance of pedigree relationships relative to other potential sources of accuracy like LD between markers and QTL was convincingly substantiated by Wientjes et al. (2012). In the human genetics context, De Los Campos et al. (2013) have derived an upper limit for the prediction accuracy that is a function of the accumulated relationship between individuals in the training and testing sets, respectively. The importance of close relatives for achieving highly accurate predictions was also observed in an animal breeding context (Legarra et al. 2008; Habier et al. 2010).

Because of the rapidly expanding arrays of genotyped lines, the number of T1 and T0 hybrids will eclipse the number of T2 hybrids. For example, if 1000 lines are available per heterotic group, of which 100 are parents of hybrids in the training set, the number of T1 hybrids reaches 180,000 and the number of T0 hybrids a staggering 810,000, while there are “only” 10,000 T2 hybrids (minus those in the training set). Thus, by sheer numbers, the best hybrids are most likely found among T1 and T0 hybrids. However, owing to the lower prediction accuracies, it will be more difficult to identify them, compared to identifying superior T2 hybrids. Breeders are unlikely to rely solely on genomic predictions when selecting potential hybrids for commercialization. Rather, genomic prediction will be employed as an initial stage in a multistage selection scheme, involving field testing of the most promising experimental hybrids. The number of experimental hybrids that can be tested in such a manner is limited by budget constraints. For practical application of genomic prediction, it is therefore important to investigate how the preselection of hybrids should be informed by the different prediction accuracies observed in the three groups.

In an earlier study on genomic prediction of hybrid performance for GY and GM, Massman et al. (2013) also found high prediction accuracies with training set sizes comparable to ours. The most likely explanation for the high prediction accuracies generally observed is that both H² and the realized relationships among parental lines tend to be very high in commercial maize breeding programs. For example, our high estimates of H² were for both traits in close agreement to those of Schrag et al. (2006) and Massman et al. (2013), the latter of which analyzed data from a U.S. corn-belt breeding program. Massman et al. (2013) also found similarly high pairwise realized relationships to those in our study (details not shown). For prediction of breeding values under an additive genetic model, a trait with high H² is expected to have higher r_A values than a trait with low H² (Daetwyler et al. 2010). In our study, however, the r_A values observed for GM, which had a considerably higher H² than GY, were higher than r_A for GY only for T2 hybrids, but lower for T1 and T0 hybrids. Interestingly, Massman et al. (2013) reported exactly the same findings, with GM having higher r_A values than GY for T2 hybrids but lower values for T1 hybrids (there were no T0 hybrids in their study). Regional differences in LD, as found in our study, are a possible explanation why the relationship between heritability and prediction accuracy differs strongly among traits (Habier et al. 2013). We hypothesize that the contribution of information from LD to the prediction accuracy differs between T2, T1, and T0 hybrids. Regional differences in LD, therefore, can also explain why the relationship between heritability and prediction accuracy is inconsistent not only among traits, but also between T2, T1, and T0 hybrids.

Composition of training set

Prediction accuracy increased with increasing training set size N_H, as expected (Table 2). However, the increase was relatively small, even when N_H was tripled. This is in contrast to studies on genomic prediction of additive breeding values in plant breeding, where tripling N_H could double the accuracy (Asoro et al. 2011; Technow et al. 2013). One explanation for this is the already rather high level of prediction accuracy reached. On the other hand, accuracy increased for T2 hybrids more than for T0 hybrids, even though their prediction accuracy was already higher for small N_H.

The key point is that increasing N_H for constant N_D and N_F does not eliminate the weakness of limited sampling of different GCA effects from each parental germplasm pool. It only increases the number of crosses ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}},$ in which a line is tested, i.e., the number of replicates per GCA effect. Thus, the precision of estimates of GCA effects of tested lines is increased, but under a high H², as in our study, this has only little impact on r_A. Another, more important consequence of increasing the number of hybrids per line is that separation of GCA and SCA effects becomes easier, improving the predictability of both. However, the contribution of SCA variance to total genetic variance was comparatively small in our data, which again limits the benefit of increasing N_H under constant N_D and N_F. Therefore, increasing N_H when N_D and N_F are constant, i.e., increasing ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}},$ might have a greater impact on r_A under low H² and in crops or breeding programs with less defined or no heterotic groups, where the relative contribution of SCA variance is expected to be larger (Reif et al. 2007).

Nonetheless, increasing the number of hybrids ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}}$ for lines serving as parents of hybrids in the training set will increase the prediction accuracy for GCA effects of these lines, which is especially beneficial under low contribution of SCA variance, because then the performance of a hybrid can be approximated by the sum of the parental GCA effects. This explains why prediction accuracy of untested hybrids profits most from increasing N_H if both parents were parents of other hybrids in the training set (T2 hybrids) and least if none of them were involved in any training set hybrid (T0 hybrids). As expected under this rationale, the r_A increase for T1 hybrids, of which only one parental line is a parent of hybrids in the training set, was between that of the T2 and T0 hybrids.

Increasing the number of tested lines N_D and N_F while keeping N_H constant decreased the number of hybrids per line ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}}$ in the same manner as decreasing N_H while keeping N_D and N_F constant. The same reasoning therefore applies here, which explains why the prediction accuracy of T2 hybrids decreased when N_D and N_F were increased. In contrast to the scenario of constant N_D and N_F and varying N_H, the decrease in ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}}$ did not lead to decreasing prediction accuracy for GY of T1 and T0 hybrids and for GM, the r_A values even increased. The reason is that increasing N_D and N_F not only decreased ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}}$ but also at the same time widened the array of germplasm covered in the training set. Thus, untested parent lines of T1 and T0 hybrids are represented better, which improved prediction accuracy of their GCA effects. Not all allele combinations encountered in T1 and T0 hybrids are present in narrow training sets and, consequently, their effects are not predictable. More diverse and larger training sets, therefore, might also improve predictability of SCA effects.

An increase of 40 lines and a 34% and 53% decrease in ${\overline{c}}_{\text{D}}$ and ${\overline{c}}_{\text{F}},$ respectively, between the high and low N_D and N_F scenarios might have been too small to observe major effects. Further research is warranted to design the training set in an optimum manner so that the prediction accuracies of T2, T1, and T0 hybrids are balanced in a way that achieves maximum selection gain across all three groups. In the short term, this is possible only with simulations, because the resources required for phenotyping larger factorials are prohibitive.

Habier et al. (2010) studied the role of pedigree relationships on accuracy of genomic prediction in German Holstein cattle. They found that prediction accuracy decreased only slightly when the training set size was halved as long as the number of close relatives per validation set individual remained constant. In a study on genomic prediction in maize, Albrecht et al. (2011) also observed that the drop in prediction accuracy was small when the training set size was halved. This was most likely a consequence of the presence of close relatives, too, because in their cross-validation scheme, an individual in the validation set had in most cases several full sibs in the training set. These studies demonstrate the disproportional importance of the closest relatives for prediction accuracy. Per definition, hybrids from the T2 and T1 groups always have very close relatives in the training set that share 50% of their genome. However, for hybrids from the T0 group, too, the maximum genomic relationships (i.e., estimated from marker data) to hybrids in the training set remained virtually unchanged by varying N_H or N_D and N_F (results not shown). This is a consequence of the high degree of relationship between the inbred lines in a closed, medium-sized breeding program. The rather small differences in prediction accuracy between the various scenarios investigated reflect the presence of close relatives that determined prediction accuracy. A comparison of the prediction accuracies of our genomic methods with those of pedigree-based methods (Bernardo 1996) could be used to quantify the contribution of pedigree relationships to the prediction accuracy. However, different from the situation in animal breeding, pedigrees of our lines were often incomplete and rarely extended more than two generations. A simulation study, in which pedigree relationships are known without error and maximum relationships can be varied, might help to clarify the role that close relatives have in determining the accuracy with which hybrid performance can be predicted.

Estimates of r_A showed a larger variation for the T0 group compared to the T2 and T1 groups. The lower variability for the T1 and T2 groups can be explained by the guaranteed presence of close relatives in the training set. The maximum relationship of T0 group hybrids with hybrids in the training set fluctuates between replications (even though it is high and similar between scenarios, on average). However, technical limitations, such as sampling constraints and different sizes of validation groups, most likely also contributed to the observed differences.

We acknowledge the excellent technical assistance of the staff members of the maize breeding program of the University of Hohenheim, especially Dietrich Klein, in conducting part of the field trials for this study. We are indebted to the group of Ruedi Fries, from Technische Universität München, for the SNP genotyping of the parental lines. We also thank the three anonymous reviewers of the initial version of the manuscript for their very helpful input and suggestions. This research was funded by the German Federal Ministry of Education and Research within the AgroClustEr “Synbreed—Synergistic plant and animal breeding” (FKZ:0315528d).