In our previous study, a doubled haploid population of 202 lines (TNDH) was developed by microspore culture from the F₁ cross between Tapidor (European winter-type rapeseed cultivar) and Ningyou7 (Chinese semi-winter type rapeseed cultivar; Qiu et al., 2006). The 202 DH lines were used to generate an immortalized F₂ hybrid population (TNRC-F₂) with 404 single crosses. Each DH line served as a parent for single-cross hybrids (Shi et al., 2011). The lines of TNDH and hybrids of TNRC-F₂ were evaluated with Tapidor and Ningyou7 at three different environments, i.e., year-location combinations, in China (Supplementary Table 1). Details of the field evaluation were published elsewhere (Shi et al., 2011). Briefly, the plot size was 3.0 m² with a distance of 40 cm between rows and 25 cm between individuals. The average dry weight of the seeds was determined as Mg ha⁻¹.

Furthermore, 37 genotypes were sampled as an independent validation population. The validation population included 37 single-cross hybrids not considered in the TNRC-F₂ population and were derived from the crosses among two lines of the 180 DH lines and 18 new DH lines derived from the TN (Tapidor × Ningyou7) F₁ cross. The validation population was grown in one environment (2015–2016) in Wuhan in a trial with three replicates. Every plot comprised three rows with a total plot size of 1.8 m². The average dry weight of the seeds was determined as Mg ha⁻¹.

We genotyped 180 out of the 202 DH lines, as well as the two parents, Tapidor and Ningyou7, using a 60,000 SNP array (Zhang et al., 2016). After quality control, 13,753 SNP markers remained, which were polymorphic, with missing values <5% and a minor allele frequency (MAF) >5%. We further removed SNPs in a perfect linkage disequilibrium that resulted in 1,527 unique SNPs (Zou et al., 2016). The genotypes of the 318 TNRC-F₂ hybrids (Supplementary Figure 1) were inferred based on the genotypes of their respective parents. For the independent validation population with 37 new single crosses, we also genotyped the 18 new parental lines using the same SNP array and SNP marker filter parameters to get the genotypes of the 37 crosses based on the genotypes of the parents.

All quantitative genetic parameters were estimated based on the performance of 318 hybrids and the 180 DH parents. After outlier tests (Anscombe and Tukey, 1963), the adjusted means of the genotypes within each environment were estimated based on the following mixed model: yin=μ+gi+εin, where y_in was the seed dry weight of the nth observation for the ith genotype, μ was the mean value, g_i was the genetic value of the ith genotype, and ε_in was the corresponding residual. The mean value and genetic value were modeled as fixed effects, and the residual was treated as a random effect.

The adjusted means of seed dry weight for each genotype across environments were estimated with the following model: yik=μ+gi+lk+εik, where y_ik was the adjusted mean of the ith genotype in the kth environment, μ was the overall mean, g_i was the genetic value of the ith genotype, l_k was the effect of the kth environment confounded with replicate effects, and ε_ik was the residual. The mean and genetic values were modeled as fixed effects, environment effects, and residuals were treated as random effects.

In addition, we estimated the genetic variance components of hybrids and parental lines, as well as the variance of genotype × environment interactions using a one-step model: yijkn=μ+a+pij+(pl)ijk+mi+fj+sij+(ml)ik+(fl)jk+              (sl)ijk+lk+eijkn, where y_ijkn was the phenotypic performance of the nth observation for the ijth entry (line i = j, or hybrid i ≠ j) in the kth environment, μ was the overall mean, a was the group effect for lines and hybrids, p_ij was the genetic effect of the parental lines, (pl)_ijk was the interaction effect of the ijth parental line with the kth environment, m_i was the GCA effect of the ith male line, f_j was the GCA effect of the jth female line, s_ij was the SCA effect of crosses between lines i and j, (ml)_ik and (fl)_jk were GCA × environment effects of female and male lines, and (sl)_ijk was SCA × environment interaction effects, l_k was the effect of the kth environment, and e_ijkn was the residual. All effects were treated as random effects, except for the mean value and group effects.

The broad-sense heritability was calculated as the ratio of genotypic to phenotypic variance: H2= σG2σG2+σG×E2NE+σE2NE∗NR, where N_E referred to the number of environments, N_R was the average number of replications per environment, σG2 was the genotypic variance, σG×E2 was the variance of genotype multiplied by environment interaction, and σE2 referred to the error variance.

Design matrices for additive and dominance effects were specified for the hybrids and their parental lines according to the F_∞ metric (Falconer and Mackay, 1996). Data from each environment were used in QTL mapping. To correct for potential population stratification, the one minus the Rogers' distance matrix was used as a kinship matrix in the genome-wide QTL mapping scan (Zhao et al., 2013). Genome-wide scans for marker-trait associations were conducted to detect main-effect QTL, as well as all first-order epistasis effect QTL.

For main effect QTL, the model was defined as the following Yu et al. (2006): Y=Xβ+Ss+Zu+e.

Y stands for the adjusted entry means of the 498 genotypes, i.e., 180 DH lines and 318 TNRC-F₂ hybrids, of each environment, β was a vector of environment effects, s was a vector of SNP effects, u was a vector of polygene background effects, and e was a vector of residual effects. X, S, and Z were incidence matrices relating Y to β, s, and u. β and s were treated as fixed effects, while u and e were treated as random effects. A Quantile-Quantile plot was used to test for proper control of population stratification (Yu et al., 2006). The Bonferroni–Holm procedure (Holm, 1979) was applied to correct for multiple tests at a significance level of P < 0.1.

We performed a full two-dimensional scan to detect epistatic QTL using the following model: Y=Xβ+ZA1a1+ZA2a2+ZD1d1+ZD2d2+ZAAi11           +ZA1D2i12+ZA2D1i21+ZDDi22+Zu+e, where a₁, a₂, d₁, and d₂, are additive and dominance effects of the two loci and i₁₁, i₁₂, i₂₁, and i₂₂ correspond to all four different epistatic interaction effects. Z_A1, Z_A2, Z_D1, Z_D2, Z_AA, Z_A1D2, Z_A2D1, and Z_DD were incidence matrices for the effects defined above. In this model, all effects were treated as fixed effects except for u and e, which were treated as random effects. The epistasis model was implemented using the efficient mixed-model association (EMMA) approach, which significantly reduced the computation time (Kang et al., 2008). A permutation analysis with 1,000 repetitions was employed to correct for multiple testing of epistatic effects at a significance level P < 0.05 (Churchill and Doerge, 1994). The proportion of the phenotypic variance explained by single QTL was estimated using multiple regression with QTL ordered according to their P-values (Utz et al., 2000). The proportion of explained genotypic variance was determined as a proportion of explained phenotypic variance standardized by broad-sense heritability.

We applied five-fold cross validation to evaluate the accuracy to predict the genotypic values from the marker effects. The data set was randomly divided into an estimation set (100% of inbred lines and 80% of the hybrids) and a test set (20% of the remaining hybrids). QTL detection and estimation of marker effects were performed in the estimation set. The estimated marker effects were then used to predict the performance of the genotypes in the test set. The prediction accuracy was estimated as a Pearson correlation coefficient between predicted and observed values standardized with the square root of the heritability.

Based on the adjusted entry means of the 498 genotypes, we applied genomic best linear unbiased prediction (Vanraden, 2008; Zhao et al., 2015a) that considered additive, dominance, and epistatic effects (Zhao et al., 2015a). The model including only additive effects was: y=1nμ+ga+e.

The model including additive and dominance effects was: y=1nμ+ga+gd+e.

Y were the adjusted entry means of the 498 genotypes across the three environments, 1_n was a vector of ones and n was the number of genotypes, μ referred to the overall mean across all three environments, and g_a and g_d represented the additive and dominance effects, respectively. In the model, μ was a fixed effect, and the remaining effects were all random effects following normal distributions ga~N(0,Gaσa2), gd~N(0,Gdσd2), and e~N(0,Iσe2), where G_a and G_d, were the relationship matrices corresponding to additive and dominance genetic effects, and σa2, σd2, and σe2 were the variance of additive effects, dominance effects, and the residuals. Details on the implementation of these relationship matrices can be found in Zhao et al. (2015a). We used a shrinkage method to calculate the relationship matrices (Ledoit and Wolf, 2004; Endelman and Jannink, 2012). All of the above GBLUP models were implemented using the R package BGLR (Perez and Campos, 2014). In addition, we employed a G-BLUP approach using a Gaussian kernel to evaluate the relevance of epistasis for prediction accuracies. The Gaussian kernel method has the benefit of considering all higher-order epistatic interact effects. The G-BLUP model using the Gaussian kernel was implemented by using the R package rrBLUP (Endelman, 2011). The prediction accuracies were evaluated using the cross-validation scenarios outlined above. We also evaluated the prediction accuracy for the independent validation set using the same training set.

The Best Linear Unbiased Estimations of seed yield for the single environments were significantly (P < 0.01) correlated, with Pearson moment correlation coefficients ranging from 0.48 to 0.59 (Supplementary Figure 2). The analyses across environments revealed genetic variance components, which were significantly (P < 0.05) larger than zero (Table 1). The broad-sense heritability estimates were 0.47 for the TNRC-F₂ hybrid population and 0.48 for the TNDH parental population. We further decomposed the genetic variance for seed yield of the TNRC-F₂ hybrid population into variances due to general (σGCA2) and specific combining ability effects (σSCA2). The σSCA2 was 1.4 times larger than σGCA2, and the variance of interaction effects between environments and SCA effects σEnvironment×SCA2 was 2.3 times larger than σEnvironment×GCA2.

Substantial transgressive variation was observed when comparing the seed yield of the TNDH parental population with the performance of the two founder lines, Tapidor and Ningyou7 (Figure 1). The average seed yield of the 318 single-cross hybrids of the TNRC-F₂ population was 2.09 Mg ha⁻¹ and was 1.27 times larger than the average seed yield of the 180 lines of the TNDH parental population.

Genome-wide QTL mapping was performed in the combined TNRC-F₂ hybrid and TNDH parental population. The Quantile-Quantile plots for the additive and dominance effects revealed that population stratification was properly controlled with a model including the kinship matrix (Supplementary Figure 3). In the genome-wide QTL mapping scan, there were 2 SNPs located on chromosomes A03 and A07, which significantly contributed to the additive genetic variation for seed yield (Figure 2). The SNP located on chromosome A03 explained 7% and the SNP on chromosome A07 explained 4% of the total genotypic variance (Supplementary Table 2). Moreover, 32 SNPs exhibited significant dominance effects and were mainly located on chromosomes A01, A02, A04, and C02. The SNPs with significant dominance effects explained 44% of the genotypic variance.

A full two-dimensional scan for epistatic effects was performed, which revealed 18 significant additive by additive effects, 129 significant additive by dominant effects, and 104 significant dominance by dominance effects (Figure 2). The epistatic interactions mainly involved SNPs, which were located on chromosome A03, A09, and C04. All 251 epistatic effects explained 31% of genetic total genetic variance (Supplementary Table 3).

Five-fold cross-validation was performed for the genome-wide QTL mapping study to obtain unbiased estimates of the total genetic variation explained by additive, dominance, and epistatic effects. The prediction accuracy was measured as a Pearson moment correlation between predicted and observed values of the test population standardized with the square root of the broad-sense heritability. The prediction accuracy of additive effects amounted to 0.18 and increased to 0.26 when modeling additional dominance effects. Adding epistatic interaction effects did not alter the prediction accuracy of the hybrid performance.

The accuracies of three different genome-wide prediction models were compared by applying five-fold cross validation (Figure 3). We observed an increase in the prediction accuracy from the model, which considered only additive effects (0.49) compared with the model with additive plus dominance effects (0.65) and with the highest accuracy of 0.72 observed when modeling additive, dominance, and epistatic effects. To summarize, the benefits were twice as much when including additionally dominance effects compared with epistatic effects.

The accuracy of the prediction model developed based on the TNDH and TNRC-F₂ populations was further evaluated in an independent validation. Therefore, an additional set of 37 RC-F₂ single-cross hybrids were generated and evaluated at an independent environment. The prediction accuracy was assessed again by standardizing the Pearson moment correlation between predicted and observed values by the square root of the broad-sense heritability. The latter was estimated based on the variance components observed for the TNRC-F₂ population. The prediction accuracy was high and amounted to 0.49 for the genome-wide prediction approach while also considering additive effects. Interestingly, the prediction accuracy dropped severely when additionally dominance (0.15) or epistatic effects (0.08) were included.

Distinguishing the variance of the hybrid performance into σGCA 2and σSCA2 provides the first insights into the role of additive and dominance effects. Assuming there is an absence of epistasis, σGCA 2is determined for immortalized F₂ populations through additive effects and σSCA2 was controlled by dominance effects (Lynch and Walsh, 1998). Thus, the predominance of σSCA2 was 1.4 times larger than σGCA2 (Table 1) and clearly points to the relevance of dominance effects in the TNRC-F₂ hybrid rapeseed population. Our findings agreed with previous results on grain yield reported for F₂ populations in maize (Wolf et al., 2000) and rice (Li et al., 2010) generated using the North Carolina Design III (Comstock and Robinson, 1952). By contrast, studies using factorial mating designs based on rapeseed (Brandle and McVetty, 1989), maize (Parisseaux and Bernardo, 2004), wheat (Zhao et al., 2015a), or barley diversity panels (Philipp et al., 2016) observed that σGCA2 is the main component of the variance of the hybrid performance. These discrepancies can be explained by allele frequencies substantially deviating from one-half in factorial crosses among diverse inbred lines, which severely impacts the ratio of σGCA2 vs. σSCA2 (Falconer and Mackay, 1996).

The findings of the genome-wide QTL mapping study (Supplementary Table 2) complement the picture of the relevance of dominance effects at the molecular level. We detected 16 times as many marker-trait associations that contributed to the dominance more than to the additive variation (Figure 2). The findings were in line with earlier studies reporting that dominance effects are key factors of heterosis in rapeseed (Radoev et al., 2008; Basunanda et al., 2010; Shi et al., 2011; Bu et al., 2015; Wen et al., 2015). Nevertheless, it is important to note that these results were not cross-validated, which was recommended to obtain an unbiased picture of the contribution of genetic effects to phenotypic variation (Utz et al., 2000). When we applied cross validation, the variance explained by all main effect QTL decreased from 55 to 7%, but the advantage for prediction accuracy when fitting data besides the additive, as well as the dominance effects, remains substantial (44%).

The small proportion of genetic variance explained by all QTL indicates a complex genetic architecture of seed yield in rapeseed, with the presence of many QTL each contributing only a little to the phenotypic variation. Genome-wide predictions are more suitable to tackle such complex traits (Riedelsheimer et al., 2012). Similarly, we observed a 2.5 times higher accuracy for the genome-wide predictions based on main effects (Figure 3) compared with the approach using QTL only. Nevertheless, the benefits of modeling the effects of the additive along with the dominance effects was comparable (33%, Figure 3) to the results observed for QTL mapping (44%). In summary, the phenotypic data analyses, the QTL mapping, and the genome-wide prediction study suggested that dominance effects were substantially contributing to the phenotypic variation of seed yield for the Tapidor × Ningyou7 rapeseed hybrid.

According to previous results (Radoev et al., 2008; Basunanda et al., 2010; Shi et al., 2011; Bu et al., 2015; Wen et al., 2015), a large number of 251 significant digenic epistatic effects were detected in the genome-wide QTL mapping study (Supplementary Table 3). Around 31.5% (79) of the digenic interaction were additive by dominance effects involving SNPs on chromosome A03 and C04 (Supplementary Table 3). These 79 epistatic pairs trace back to eight clusters of markers, five on chromosome A03, and three on chromosome C04. Some of the SNPs of the eight clusters were also associated with main effect QTL for flowering time, seed weight, seed number, and oil content (Supplementary Table 3). Marker Bn-scaff_16534_1-p2320270 located on C04 showed significant dominance effects (Supplementary Table 2) and was together with closely linked markers involved in 26 of the epistatic interactions between chromosomes A03 and C04. Thus, Bn-scaff165341-p2320270 is an interesting candidate for further fine mapping studies.

Despite the large number of 251 significant digenic epistatic effects, the cross-validated prediction accuracy did not increase when modeling epistatic effects beyond main effects. This can either point to the irrelevance of epistasis for seed yield in rapeseed or can be explained by the challenge to detect epistatic effects for the complex trait seed yield (Mackay, 2014). The results of the genome-wide prediction study support the latter explanation. The prediction accuracy increased by 11% compared with the prediction approach based on additive and dominance effects only (Figure 3). The observed increase is much higher than that reported in previous studies on genome-wide prediction of hybrid performance. Zhao et al. (2015a) observed a 2% increase in prediction ability when modeling epistatic effects in a large population of 1,604 wheat hybrids and their 135 parental inbred lines. The lower benefit observed for wheat vs. rapeseed can be explained by the more marginal allele frequencies in factorial mating designs compared with immortalized F₂ design. Surprisingly, Xu et al. (2014) observed no benefits when modeling epistatic effects in a large immortalized F₂ population comprised of 240 rice inbred lines and 360 F₂ genotypes. The discrepancy between our findings and that of Xu et al. (2014) might point to differences in the genetic architecture among crop species, as has been observed for the genetic basis of heterosis in the self-pollinating species rice vs. the outcrossing species of maize (Garcia et al., 2008) and deserves further research.

The total population was used to estimate the genetic components of variance with a Bayesian generalized linear regression based on the design matrices of additive, dominance, and all types of digenic epistatic effects. We observed a prominent role of the variance of dominance effects contributing to 49% of the genetic variance (Supplementary Figure 4). Additive (27%) and the sum of all digenic epistatic effects (24%) contributed nearly equally to the genetic variance of the hybrid performance. This suggests that the Tapidor × Ningyou7 rapeseed hybrid successfully exploits all types of genetic effects. Whether our findings are also valid for further rapeseed hybrids or not, deserve further research in other genetic backgrounds.

The phenotypic data analyses revealed that σGCA 2 was less affected by varying environmental conditions than σSCA2 (Table 1). This result contrasted with previous findings based on an F₂ population derived from the B73 × Mo17 maize hybrid (Wolf et al., 2000), which again points to differences in the genetic architecture among crop species. Validating the developed prediction model in an independent sample of genotypes and environments revealed that the observation at the phenotypic level was also reflected at the molecular level. The prediction accuracy is stable for additive effects (0.49) but collapsed when adding dominance (0.15) or epistatic effects (0.08) (Figure 3). Thus, the development of robust hybrid prediction models, including dominance and epistatic effects, requires a much deeper phenotyping analysis in multi-environmental trials. The flip side of our finding is that prediction models focusing exclusively on additive effects yields already stable and high prediction accuracies. Additive models can be easier implemented simplifying the application of hybrid prediction in plant breeding programs.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.