Genomic Prediction of Genotypic Effects with Epistasis and Environment Interactions for Yield-Related Traits of Rapeseed (Brassica napus L.)

Oilseed rape (Brassica napus) is an economically important oil crop, yet the genetic architecture of its complex traits remain largely unknown. Here, genome-wide association study was conducted for eight yield-related traits to dissect the genetic architecture of additive, dominance, epistasis, and their environment interaction. Additionally, the optimal genotype combination and the breeding value of superior line, superior hybrid and existing best line in mapping population were predicted for each trait in two environments based on the predicted genotypic effects. As a result, 17 quantitative trait SNPs (QTSs) were identified significantly for target traits with total heritability varied from 58.47 to 87.98%, most of which were contributed by dominance, epistasis, and environment-specific effects. The results indicated that non-additive effects were large contributions to heritability and epistasis, and also noted that environment interactions were important variants for oilseed breeding. Our study facilitates the understanding of genetic basis of rapeseed yield trait, helps to accelerate rapeseed breading, and also offers a roadmap for precision plant breeding via marker-assisted selection.

Introduction

Rapeseed (Brassica napus L.; AACC, 2n = 38) is one of the most important sources of vegetable oil in the world. Owing largely to the development of high yielding varieties of rapeseed in breeding programs, now hybrid rape has the largest market share due to its superior performance for multiple agronomic traits (Liu et al., 2002). While the genetic architectures of complex traits is still vague.

Genome-wide association studies (GWAS) are now becoming one of the standard methodology for revealing the genetic variation of complex traits for human, plants, and animals (Clark, 2010; Cockram et al., 2010; Zhang et al., 2010; Flint and Eskin, 2012; Segura et al., 2012; Li et al., 2013). Several GWAS studies in B. napus have been reported for seed quality traits (Li F. et al., 2014; Liu S. et al., 2016), branch angle (Liu J. et al., 2016), flower time (Wang et al., 2016; Xu et al., 2016), and yield-related traits mainly containing seed weight (Li F. et al., 2014), plant height (Li et al., 2016), primary branch number (Li et al., 2016) in recent years. However, these potential loci identified with GWAS methods only explain a small amount of phenotypic variance (4.3~25.17% for seed quality traits, 16.47~20.51% for branch angle, 5.62~15.75% for flower time and 3.72~13.87% for yield-related traits). A probable reason for the low heritability detected with GWAS is that only additive models are applied but not considering dominance, epistasis (Zuk et al., 2012; Hemani et al., 2013). In a study of Drosophila melanogaster populations, about 50% of phenotypic variation in adult olfactory behavior was assigned to genotype-by-environment (G×E) interaction (Sambandan et al., 2008). Similarly, G×E or epistatic interactions could explain considerable proportion of variance of flowering traits in rice (Uwatoko et al., 2008) and in Arabidopsis (Caicedo et al., 2004; El-Lithy et al., 2006). Thus, trying to find the so-called “missing heritability” could help to efficiently dissect the genetic mechanism of complex traits (Manolio et al., 2009; Ingvarsson and Street, 2011; Resende et al., 2012).

In the earlier studies of QTL mapping, epistasis was observed for the resistance to Sclerotinia sclerotiorum in rapeseed, and the additive by additive interactions were the predominant type of epistasis (Zhao and Meng, 2003). Additionally, total 190 plant architecture (PA)-related candidate genes for 91 unique PA QTLs and 2,350 plant yield (PY) epistatic interaction loci-pairs were identified, which explain 2.8–51.8 and 5.2–23.6% of phenotypic variation, respectively (Cai et al., 2016). However, the published studies relied on kinds of low density traditional markers. Recently, the revolution of inexpensive sequencing has been used for the identification of high-throughput genomic SNPs (Ganal et al., 2009), and a new mixed linear model approach has been developed to dissect the genetic architecture of multiple loci by partitioning into additive, dominance, epistasis, and environmental interactions (Zhang et al., 2015).

Rapeseed seed yield is usually determined by several factors including plant height (PH), main inflorescence length (IL), branch number (BN), number of seeds per silique (SS), effective silique number on main inflorescence (ISN), thousand seed weight (TSW), biomass yield per plant (BY), and seed yield per plant (SY). In the present study, the eight yield-related traits of 367 experimental materials were detected in two environments. GWAS performed with high-throughout SNPs was firstly used to find SNPs with additive, dominance, epistasis and environmental interactions for these traits. Additionally, prediction of genetic effects was conducted based on the detected quantitative trait SNPs (QTSs) of eight agronomic traits to predict high-yield genotype combination. The study will provide a better understanding of the genetic mechanism of yield-related traits, and the QTSs are beneficial for efficiently selecting ideal parental lines for breeding new hybrid rapeseed.

Materials and Methods

Plant Materials and Field Experiments

To generate genetic population used for the association studies, we adopted Yu5 (L155) and Zheyou18 (L157) as female parents to pollinate with another 149 male parents (L1-L149) and finally harvested 216 F₁ hybrids including 97 hybrids taken L155 as female and 119 hybrids taken L157 as female. All the lines came from China and adapted to the middle reaches of the Yangtze River agro-climatic conditions, China. The 367 experimental materials, including 151 inbred lines and 216 F₁ hybrids were grown in farm fields in Wuhan (29.58° N, 113.53° E) and Xiangyang (32.01° N, 112.08° E) in Hubei province, China, under normal conditions for crop production during the rapeseed growing season of 2012–2013. The annual mean temperature of Wuhan was 1~2° higher than Xiangyang, and the average rain days of Wuhan was 3 days more than Xiangyang during the rapeseed growing season of 2012–2013.

In 2012, all trials were designed as randomized complete blocks with three replications in each environment. Each plot consisted of three rows of 3.5 m length with 0.25 m distance between rows. During the last 5 days of September and the first 5 days of October, seeds were sown with a distance of 0.15 cm between plants in each row. The management of the field experiments was performed in accordance with local standard production practices. At maturity, 12 plants in the middle row were harvested (randomly within this row) from each plot for evaluation of the following eight quantitative traits: plant height (PH), main inflorescence length (IL), branch number (BN), number of seeds per silique (SS), effective silique number on the main inflorescence (ISN), thousand seed weight (TSW), biomass yield per plant (BY), and seed yield per plant (SY). Measurements for PH, IL, BN, SS, ISN, TSW, BY, and SY were completed as described by Shi et al. (2011).

Genotyping and Quality Control

Genotyping of the entire inbred line population was performed commercially by Emei Tongde Co. (Beijing) using the Brassica 60 K Illumina^® Infinium SNP array (http://www.illumina.com/technology/infinium_hd_assay.ilmn) according to the manufacturer's protocol. The SNP data were clustered and called automatically using Illumina BeadStudio genotyping software. SNPs with QQ or qq frequencies equal to zero, call frequency <0.9, or minor frequency <0.05 were excluded. The remaining SNPs were scrutinized visually, and those SNPs not showing three clearly defined clusters representing the three possible genotypes (QQ, Qq, qq) were also excluded. A total of 33,689 SNP markers among the 367 lines and hybrids that met the QC criteria were used for the association analyses.

The phenotypes of 8 traits were filtered based on the following two criteria: (1) the frequency of missing data for each subject <0.1; (2) outliers were discarded according to a distribution-based outlier detection of residuals (|ε − μ_ε|/ σ_ε > 3) by QQ plot (Quantile-Quantile plot) for each trait (Figure S1).

Statistical Analysis

The genetic model for the phenotypic value of the k-th line or hybrid in the h-th environment (y_hk) can be expressed by the following linear mixed model:

$\begin{array}{cc}\begin{array}{l}{y}_{hk}=\mu +{\displaystyle \sum _i{a}_i{x}_{A_{ik}}}+{\displaystyle \sum _i{d}_i{x}_{D_{ik}}}+{\displaystyle \sum _{i<j}a{a}_{ij}{x}_{A{A}_{ijk}}}\\ +{\displaystyle \sum _{i<j}a{d}_{ij}{x}_{A{D}_{ijk}}}+{\displaystyle \sum _{i<j}d{a}_{ij}{x}_{D{A}_{ijk}}}+{\displaystyle \sum _{i<j}d{d}_{ij}{x}_{D{D}_{ijk}}}+e_h\\ +{\displaystyle \sum _{i}a{e}_{ih}{u}_{A{E}_{ihk}}}+{\displaystyle \sum _{i}d{e}_{ih}{u}_{D{E}_{ihk}}}+{\displaystyle \sum _{i<j}aa{e}_{ijh}{u}_{AA{E}_{ijhk}}}\\ +{\displaystyle \sum _{i<j}ad{e}_{ijh}{u}_{AD{E}_{ijhk}}}+{\displaystyle \sum _{i<j}da{e}_{ijh}{u}_{DA{E}_{ijhk}}}\\ +{\displaystyle \sum _{i<j}dd{e}_{ijh}{u}_{DD{E}_{ijhk}}}+{\epsilon }_{hk}\end{array}& \left(1\right)\end{array}$

where μ is the population mean; a_i is the additive effect of the i-th locus with coefficient u_Aik (1 for homozygote major alleles QQ and −1 for homozygote minor alleles qq); d_i is the dominance effect of the i-th locus with coefficient u_Dik (1 for heterozygote Qq, 0 for homozygotes QQ and qq); aa_ij, ad_ij, da_ij, and dd_ij are the digenic epistasis effects with coefficients of random variables u_AAijk (1 for QQ × QQ and qq × qq, −1 for QQ × qq and qq × QQ), u_ADijk (1 for QQ × Qq, −1 for qq × Qq), u_DAijk (1 for Qq × QQ, −1 for Qq × qq), and u_DDijk (1 for Qq × Qq); e_h is the random effect of the h-th environment; ae_ih is the additive × environment interaction effect of the i-th locus in the h-th environment with coefficient u_AEihk; de_ih is the dominance × environment interaction effect of the i-th locus in the h-th environment with coefficient u_DEihk; aae_ijh, ade_ijh, dae_ijh, and dde_ijh are the digenic epistasis × race interaction effects in the h-th ethnic population with coefficients of random variables (u_AAEijhk, u_ADEijhk, u_DAEijhk, and u_DDEijhk); and ε_hk is the residual effect of the k-th line or hybrid in the h-th environment.

Heritability of individual genetic effects were estimated by $h_g^2=\alpha {\sigma }_g^2/V_P$ (α = 2 for additive effect, α = 1 for dominant effect, α = 4 for additive × additive, α = 2 for additive × dominant or dominant × additive, α = 1 for dominant × dominant), where phenotypic variance (V_P) is the sum of genetic variance (V_G), genetic by environment interaction variance (V_GE), and residual variance (V_ε),

$\begin{array}{cc}\begin{array}{l}{V}_P=V_G+V_{GE}+V_{\epsilon }\\ =(V_A+V_D+V_I)+(V_{AE}+V_{DE}+V_{IE})+V_{\epsilon }\\ =(V_A+V_D+V_{AA}+V_{AD}+V_{DA}+V_{DD})+\\ (V_{AE}+V_{DE}+V_{AAE}+V_{ADE}+V_{DAE}+V_{DDE})+V_{\epsilon }\end{array}& \left(2\right)\end{array}$

The total heritability was estimated by

$\begin{array}{l}{h}_{G+GE}^2=(h_A^2+h_D^2+h_I^2)+(h_{AE}^2+h_{DE}^2+h_{IE}^2)\\ =(h_A^2+h_D^2+h_{AA}^2+h_{AD}^2+h_{DA}^2+h_{DD}^2)\\ +(h_{AE}^2+h_{DE}^2+h_{AAE}^2+h_{ADE}^2+h_{DAE}^2+h_{DDE}^2)\\ ={\displaystyle \sum _i{h}_a^2}+{\displaystyle \sum _i{h}_d^2}+{\displaystyle \sum _{i<j}{h}_{aa}^2}+{\displaystyle \sum _{i<j}{h}_{ad}^2}+{\displaystyle \sum _{i<j}{h}_{da}^2}\\ +{\displaystyle \sum _{i<j}{h}_{dd}^2}+{\displaystyle \sum _i{h}_{ae}^2}+{\displaystyle \sum _i{h}_{de}^2}+{\displaystyle \sum _{i<j}{h}_{aae}^2}+{\displaystyle \sum _{i<j}{h}_{ade}^2}\\ +{\displaystyle \sum _{i<j}{h}_{dae}^2}+{\displaystyle \sum _{i<j}{h}_{dde}^2}\end{array}$

We used the GMDR module (Generalized Multifactor Dimensionality Reduction; Qi et al., 2013) in the QTXNetwork software (http://ibi.zju.edu.cn/software/QTXNetwork/) to scan 33,689 SNP markers in 367 subjects for 1D~3D significant candidate SNP markers, and obtained 539 candidate SNPs (260 in the A genome, 262 in the C genome, and 17 in the Scaffold group). The QTS mapping module in the QTXNetwork was then used to dissect the genetic architecture of the eight agronomic traits of oilseed rape (Zhang et al., 2015). Significant SNPs associated with phenotypic variants were analyzed by setting a total of 2,000 permutation tests to calculate the critical P-value for controlling the experiment-wise type I error. The effects were predicted by using a Markov Chain Monte Carlo method with 20,000 Gibbs sampler iterations (Yang et al., 2007). The correlation coefficient (R_Ŷ) between predicted breeding values and phenotypic values was estimated for each trait. Based on the predicted genetics effects of SNP loci for eight traits, we predicted total genotypic effects for the best lines and the best hybrids of the mapping population, and also predicted genotypic effects of superior lines and superior hybrids to inform further selection decisions (Yang and Zhu, 2005).

Results

Estimated Heritability and Predicted Genetic Effects

We conducted a genome wide association study for eight yield-related traits of B. napus on a population including 151 inbred lines and 216 F₁ hybrids obtained from mating two female parents (L155 and L157) to other 149 inbred lines as male parents. Eight yield-related traits were investigated in our study by using a full genetic model with genetic effects of additive, dominance, epistasis, and their environment interactions. A total of 17 QTSs controlling eight yield traits were detected: 4 QTSs for PH, 2 QTSs for BN, 3 QTSs for IL, SS, ISN, TSW, BY, and SY, respectively (Figure 1 and Table S1). Some loci exhibited pleiotropic effects, including C09_M34850_G/A for six traits (BN, BY, ISN, PH, SS, and SY), A07_M11103_A/G for both BY and IL, A08_M12337_C/A for both PH and TSW, and C04_M26614_A/G for both BY and SY.

FIGURE 1

Figure 1. G×G plot of detected significant QTSs (P_EW < 0.05) for eight traits. Circle, QTS with additive effect; Square, QTS with dominant effect; Line between two QTSs, epistasis effect; Red color, QTS with general effects for two environments; Green color, QTS with environment-specific effects; Blue color, QTS with both general and environment-specific effects; Black color, QTS with significant epistasis effects but without detected individual effects.

Estimated heritability and correlation coefficients ($R_{\widehat{Y}}$) between total genotypic values of detected QTSs and phenotypic values for the eight traits are listed in Table 1. The total heritability ranged from 58.47 to 87.98%, and was contributed by various types of genetic variance effects. With the exception of three yield traits (SS, BY, and SY), which were sensitive to the environment ($h_{GE}^2\stackrel{\wedge }{=}$ 60.24% for SY, 59.75% for BY, and 49.57% for SS), the other five traits were quite stable across the two environments ($h_{GE}^2\stackrel{\wedge }{=}$ 4.14~27.40%). With the exception of the SS trait, epistasis effects contributed a large portion of total heritability ($h_{I+IE}^2\stackrel{\wedge }{=}$ 26.07~62.14%). The correlation coefficient of genetic prediction with phenotype ($R_{\widehat{Y}}$ = 0.501~0.899) for each trait was very close to the estimated heritability ($h_T^2\stackrel{\wedge }{=}$ 58.47~87.98%), indicating that this statistic approach would be quite efficient for predicting the best lines/hybrids, and selecting the superior lines/hybrids by using the predicted genetic effects.

TABLE 1

Table 1. Estimates of heritability and correlation coefficient of detected QTSs for eight traits.

Highly significant (experiment-wise P_EW–value < 10⁻⁵) predicted genetic effects of 16 QTSs are presented in Table 2. Among the eight traits studied, plant height (PH) had the highest heritability ($h_T^2\stackrel{\wedge }{=}$ 85.06%) mainly contributed by epistasis ($h_I^2\stackrel{\wedge }{=}$ 50.30%). There was one locus (G/A of C09_M34850) detected with large and positive dominance effects. Homozygotes of this locus (G/G) could also increase plant height. The additive and dominance effects ($a\stackrel{\wedge }{=}$ −7.107 of G/G, $d\stackrel{\wedge }{=}$ −3.565 of G/A) were negative for locus Scaffold_M33906_G/A, which could be used in selecting for decreased plant height. Heterozygote C/A of A08_M12337 and homozygote G/G of C09_M34850 had a negative dominance × additive epistasis effect ($da\stackrel{\wedge }{=}$ −11.003), which could dramatically decrease plant height; however this large epistasis effect was counteracted by their main effects. Instead, the combination of homozygotes for major-allele C/C of A08_M12337 and minor-allele A/A of C09_M34850 implicated a lower plant height (Figure S2A). Epistasis effects were also the most important genetic effects on main inflorescence length (IL) ($h_I^2\stackrel{\wedge }{=}$ 53.01%) and branch number (BN) ($h_I^2\stackrel{\wedge }{=}$ 44.33%). Epistasis of heterozygote A/G of A03_M4640 × minor-allele homozygote A/A of A08_M12338, and heterozygote A/G of A07_M11103 could significantly increase IL (Figures S2B,C). Heterozygote G/A of C09_M34850 was the major locus for increasing BN. And due to its epistasis effects, combination of heterozygote C/A of A01_M1364 × minor-allele homozygote A/A of C09_M34850 could also significantly increase BN (Figure S2D).

TABLE 2

Table 2. Predicted genetic effects of highly significant QTSs for eight traits.

Number of seeds per silique (SS) had very strong environment-specific additive and dominance effects ($h_{AE}^2\stackrel{\wedge }{=}$ 35.91% and $h_{DE}^2\stackrel{\wedge }{=}$ 13.66). Heterozygote G/A of A08_M12212 ($h_d^2\stackrel{\wedge }{=}$ 6.43%), G/A of C09_M34850 ($h_d^2\stackrel{\wedge }{=}$ 1.63%, $h_{d{e}_2}^2\stackrel{\wedge }{=}$ 5.95%) could increase SS in different environments.

Effective silique number on main inflorescence (ISN) had high heritability ($h_T^2\stackrel{\wedge }{=}$ 72.57%) mostly due to dominance effects ($h_D^2\stackrel{\wedge }{=}$ 35.52%) of three loci (A/G of A05_M8315, C/A of C04_M31883, and G/A of C09_M34850). Based on their large main effects, these three loci could be used in selection for increasing ISN, despite the negative dominance epistasis effect between heterozygotes of C04_M31883 and C09_M34850 (Figure S2E).

Thousand seed weight (TSW) had relatively large additive and epistasis variances ($h_A^2\stackrel{\wedge }{=}$ 17.23%, $h_{I+IE}^2\stackrel{\wedge }{=}$ 26.07%). Increasing TSW could be expected by homozygote minor-alleles A/A of A08_M12337 ($h_a^2\stackrel{\wedge }{=}$ 15.03%), and also by heterozygote C/A of C01_M30818 ($h_d^2\stackrel{\wedge }{=}$ 7.29%), but the heterozygote of both loci should be avoided due to negative epistasis effects ($dd\stackrel{\wedge }{=}$ −0.212) (Figure S2F).

For biomass yield per plant (BY), the largest contributions of genetic variance were environment-specific dominance and epistasis ($h_{DE}^2\stackrel{\wedge }{=}$ 28.88%, $h_{IE}^2\stackrel{\wedge }{=}$ 27.88%). Heterozygote A/G of A07_M11103 ($d{e}_2\stackrel{\wedge }{=}$ 23.745, $h_{d{e}_2}^2\stackrel{\wedge }{=}$ 27.84%) could significantly increase BY (in environment Xiangyang). For the most important yield trait (seed yield per plant, SY), dominance and environment-specific epistasis ($h_D^2\stackrel{\wedge }{=}$ 11.11%, $h_{IE}^2\stackrel{\wedge }{=}$ 48.18%) were the major genetic recourses for increasing yield (Figure S2G). Heterozygote G/A of C09_M34850 could increase SY across environments ($h_d^2\stackrel{\wedge }{=}$ 8.94), and add extra selection response ($d{e}_2\stackrel{\wedge }{=}$ 4.743) in environment Xiangyang ($h_{d{e}_2}^2\stackrel{\wedge }{=}$ 5.96%). Although dominance of heterozygote A/G of C04_M26614 could only slightly increase SY ($d\stackrel{\wedge }{=}$ 2.096, $h_d^2\stackrel{\wedge }{=}$ 2.16%), its epistasis interaction with dominance of another heterozygote G/A of locus C09_M34850 could also have large increase for SY in environment Xiangyang ($dd{e}_2\stackrel{\wedge }{=}$ 4.365) (Figure S2H).

Predicted Genetic Effects for Different Genotypes

The Genotype-Phenotype (G-P) maps of epistasis SNPs for each trait, in both environments were presented in Figure S2. G-P maps exhibited various patterns as effects differed among different traits. High concordance was observed between phenotypic G-P maps and genotypic G-P maps. Based on predicted genetic effects of QTSs for each trait, we further predicted the maximum and minimum genotypic effects of the superior lines and superior hybrids in two environments on eight traits. We also predicted the genotypic effects of homozygotes (QQ, qq), and heterozygote (Qq) for eight traits in two environments, respectively (Table 3). All the predicted genotypic effects of eight traits were negative for major-allele homozygote QQ, but positive for minor-allele homozygote qq. The predicted genotypic effects of heterozygote Qq were positive for all the eight traits studied. Among the eight traits for all-locus heterozygote (Qq), the predicted genotypic values were much larger than minor-allele homozygote (qq) for eight traits but not for TSW in one environment. It was implied that for this rapeseed population, hybrid breeding could potentially increase breeding values of seven yield traits but not for shortening plant height.

TABLE 3

Table 3. Predicted genetic effects in two environments for genotype of QQ, qq, Qq, best and superior lines, best and superior hybrids of eight traits.

There was no difference between the best lines of mapping population and the predicted superior lines for four traits (BN, SS, TSW, and SY). It was suggested that pure-line variety breeding might have only limited potential for improving these traits based on the QTSs detected for this rapeseed population. For trait IL, breeding value of the best lines was smaller than the predicted superior lines, which was due to one locus in the best line L155 (minor-allele homozygote C/C of A01_M15289). The superior line could be obtained by replacing minor-alleles to major-alleles (A/A) of A01_M15289 for increasing $a\stackrel{\wedge }{=}$ 1.63 in two environments and extra $a{e}_2\stackrel{\wedge }{=}$ 1.27 in E₂. For another trait ISN, the breeding value of the best line ($G+GE\stackrel{\wedge }{=}$ 9.64 in E₁, and 7.88 in E₂) was larger than the superior line ($G+GE\stackrel{\wedge }{=}$ 4.15 in E₁, and 5.24 in E₂), because the best lines (L41 in E₁ and L45 in E₂) still had heterozygote locus (C/A of C04_M31883).

Among the eight traits studied, there were six traits (BN, SS, ISN, TSW, BY, and SY) with no difference of breeding values between the best hybrids and the superior hybrids in at least one environment. It was indicated that no selection advantage could be expected based on this mapping population for improving these six traits of hybrids. Expected gain of hybrid breeding could be obtained for two traits based on the best hybrids of this mapping population. For trait IL, the superior hybrid genotypes could be selected as A/G of A03_M4640 in two environments, and A/G of A07_M11103 in E₁ based on the best hybrids (L155 × L76 in E₁, L155 × L26 in E₂) and maintained homozygote A/A for other two loci (A08_M12338, A01_M15289). There could have dramatic decrease for plant height of hybrid by just selecting all detected four QTSs of PH as major-allele homozygotes (C/C of A08_M12337, G/G of A01_M114, C09_M34850, and Scaffold_M33906).

The genotype of receptor lines and donor lines was listed in Table S2. Superior hybrid for six traits (BN, SS, ISN, TSW, BY, and SY) could be obtained via hybridization of acceptor lines (L155 and L157) to certain donor lines (L1~L154). There was only limited number of donor lines (3–13) that could contribute to target traits, except for trait ISN, which could be improved by 103 donor lines. We also found L102 was a competitive donor line of high potential. L102 could simultaneously improve four traits (SS, ISN, and SY) mainly due to its ability to donor G allele to A allele of C09_M34850, whose heterozygote could largely improve trait SS, ISN and SY with its pleiotropic effects. L102 also carried A/A of A07_M11580 and C/C of A08_M12337, which enabled L102 to improve TSW via introducing heterozygote into acceptor lines.

Discussion

There are evidences that rare variants have large impacts on common human diseases (Cirulli and Goldstein, 2010; Yang et al., 2010; Zuk et al., 2014). Experimental evidence from association and linkage populations demonstrated that the rare genetic variation at β-carotene hydroxylase 1 (crtRB1) was associated with β-carotene concentration in maize kernels (Yan et al., 2010). In our study, full genetic model including epistasis and environment-specific effects was firstly used to excavate the missing heritability and dissect genetic architecture of important agronomic traits in B. napus. The average minor-allele frequency (MAF) in the 151 cultivars of this study was 11.07% (4.67~16.67%) for the detected 17 QTSs (10 A-alleles, 6 G-alleles, and 1 C-allele), and minor alleles had impacts on various genetic effects (−a of qq, d of Qq, aa of qq × qq, ad of qq × Qq, ± da of Qq × qq, and dd of Qq × Qq) according to the genetic model. These minor alleles (qq) increased breeding values for all the eight traits and made contributions to total heritability increasing by different genetic effects.

There were total 17 main additive effects loci and 14 environment-specific additive effects loci were identified, and most of them were negative (15 negative main additive effects and 8 negative environment-specific additive effects). However, the heritability of each additive locus was quite low ($h_{ā}^2\stackrel{\wedge }{=}$ 2.27%, $h_a^2\stackrel{\wedge }{=}$ 0.12~15.03%; $h_{\overline{ae}}^2\stackrel{\wedge }{=}$ 5.03%, $h_{ae}^2\stackrel{\wedge }{=}$ 1.66~23.13%). It suggested that the contribution to phenotype due to negative alleles could not be neglected and the additive variances were not the major genetic contribution for most traits studied. While previous study indicated that the additional effects of the alleles originated from both parents were detected to be important for yield components traits in rapeseed (Wang and Guan, 2010). The inconsistent of results reiterated the complex genetic mechanism of yield and yield-related traits, especially for allopolyploid species such as B. napus. All the dominance variants were also contributed due to minor-allele heterozygotes (Qq). Dominance had much larger impacts on eight traits studied due to 21 main dominance loci with 19 having positive effects ($h_{\bar{d}}^2\stackrel{\wedge }{=}$ 6.14%, $h_d^2\stackrel{\wedge }{=}$ 0.11~25.94%) and 20 environment-specific dominance loci with 13 having positive effects ($h_{\overline{de}}^2$ $\widehat{=}$ 4.76%, $h_{de}^2\stackrel{\wedge }{=}$ 0.30~27.84%). Which indicated that positive dominance effects were conducive to yield-related traits plasticity in B. napus. To data, epistatic effects were considered as important for complex traits in crops, such as plant height (Cao et al., 2001), yield (Huang et al., 2014), and salt tolerance (Wang et al., 2012) in rice seedlings, plant height in cultivated wheat (Zhang et al., 2008) and seed protein concentration in soybean (Qi et al., 2016). There were 10 pairs of loci identified with main dominance-related epistasis (dd, ad, and da) ($h_{ī}^2\stackrel{\wedge }{=}$ 19.05%, $h_i^2\stackrel{\wedge }{=}$ 1.35~50.30%) and eight pairs of loci detected with environment-specific dominance-related epistasis (dde, ade, and dae) ($h_{iē}^2\stackrel{\wedge }{=}$ 14.11%, $h_{ie}^2\stackrel{\wedge }{=}$ 0.74~38.80%). The most contribution of heritability for yield traits was due to 12 pairs of loci identified with epistasis between dominance effects and additive effects (ad and da, ade, and dae; $h_{AD+DA+ADE+DAE}^2\stackrel{\wedge }{=}$ 10.96~53.01%). It was revealed the major role of epistasis influencing rapeseed yield. None of the detected QTSs was common with the SNPs detected in the previous association studies (Li F. et al., 2014; Li et al., 2016; Liu J. et al., 2016; Liu S. et al., 2016; Wang et al., 2016; Xu et al., 2016). It confirmed again that yield-related traits are complex polygenic phenomenon in rapeseed. Four pleiotropic QTSs were found to be associated with more than one trait (Table 2). It indicated that these traits might share part of genetic basis. Thus, the pleiotropic loci should be a priority for further research, and multi-traits should be taken into account together in genetic breeding practice.

Fully characterizing the genetic mechanism mediating heterosis is helpful for increasing crop yield. While none of the current genetic models can completely explain the heterosis phenomenon. Previous study indicated that epistasis together with all levels of dominance from partial to overdominance is responsible for the expression of heterosis in rapeseed (Radoev et al., 2008). In the study, large and positive heteroses for eight yield traits were mostly due to minor-alleles in heterozygotes (Qq), and minor-allele homozygotes epistasis (qq × Qq of AD, Qq × qq of DA) was contributed for seven yield traits. It was concluded that epistasis together with heterozygotes loci play an important role in yield heterosis in B. napus.

Selection of best genotype combination is difficult due to the complexity of the genetic architecture. For example, to increase ISN, C/A of C04_M31883 × G/A of C09_M34850 should be avoided due to their negative additive effects, but their large positive dominant effects could overturn the epistasis effects and make heterozygotes the optimal choice for these two loci. Genotype-Phenotype maps of epistasis SNPs based on prediction were adopted to visually demonstrate the accumulated genetic effects of the epistasis SNP pairs. Due to the high heritability of epistasis effects, the G-P maps based on population mean exhibited a similar pattern with corresponding G-P maps based on genetic prediction. But G-P maps based on population mean may be biased away from prediction due to confounding with effects of other loci and residual error. So G-P maps based on genetic prediction could be a better choice for selection by visualizing the true effects of epistasis effects. For loci involved in multiple epistases, selection needs more caution. For IL, A08_M12338 interacted with both A03_M4640 and A07_M11103. Although G/G of A07_M11103 combined with G/A of A08_M12338 could increase IL by $ad\stackrel{\wedge }{=}$ 4.204, A/A of A08_M12338 should be chosen to obtain higher IL based on the total effects of the three loci. In this case, prediction function of superior line and superior hybrid is conductive to select the optimal genotype combination, thus efficiently utilize the heterosis in B. napus.

In the study, there were six traits (BN, SS, ISN, BY, TSW, and SY) having no difference of predicted breeding values between the best hybrids and the superior hybrids in at least one environment (Table 3). It was suggested that these six traits were already under strong breeding selection and conserved positive effects of minor-allele homozygotes (qq) and heterozygotes (Qq) for additive, dominance and epistasis effects. All superior hybrids for these six traits could be obtained via introducing heterozygote into L155 or L157, suggesting that L155 and L157 already offered suitable genetic background. Only limited amount (3~13) of donor lines could be crossed to acceptor line to gain desirable properties, except for trait ISN. It is indicated that donor lines should be chosen with caution, and that marker-assisted selection could help to reduce blindness in hybrid breeding based on the predicted breeding values.

Author Contributions

CM conceived and designed the experiments. XL, YY, and CM performed the experiments. JZ, YD and XL analyzed the data. CM and JZ contributed materials/analysis tools. JZ, CM, LZ and XL wrote the paper. JS polished the paper.

Conflict of Interest Statement

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.

Acknowledgments

This research was supported in part by grants from the National Basic Research Program of China (973) (2011CB109306), the National Key Research and Development Program of China (No. 2016YFD0100803-04), the National High Technology Research and Development Program of China (863) (2012AA101107), and the National Science Foundation of China (31301248, 31371250). The funders had no role in study design, data collection and analysis, decision to publish, or preparation of the manuscript.

Supplementary Material

The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fgene.2017.00015/full#supplementary-material

References

Cai, G., Yang, Q., Chen, H., Yang, Q., Zhang, C., Fan, C., et al. (2016). Genetic dissection of plant architecture and yield-related traits in Brassica napus. Sci. Rep. 6:21625. doi: 10.1038/srep21625

PubMed Abstract | CrossRef Full Text | Google Scholar

Caicedo, A. L., Stinchcombe, J. R., Olsen, K. M., Schmitt, J., and Purugganan, M. D. (2004). Epistatic interaction between Arabidopsis FRI and FLC flowering time genes generates a latitudinal cline in a life history trait. Proc. Natl. Acad. Sci. U.S.A. 101, 15670–15675. doi: 10.1073/pnas.0406232101

PubMed Abstract | CrossRef Full Text | Google Scholar

Cao, G., Zhu, J., He, C., Gao, Y., Yan, J., and Wu, P. (2001). Impact of epistasis and QTL × environment interaction on the developmental behavior of plant height in rice (Oryza sativa L.). Theor. Appl. Genet. 103, 153–160. doi: 10.1007/s001220100536

CrossRef Full Text | Google Scholar

Cirulli, E. T., and Goldstein, D. B. (2010). Uncovering the roles of rare variants in common disease through whole-genome sequencing. Nat. Rev. Genet. 11, 415–425. doi: 10.1038/Nrg2779

PubMed Abstract | CrossRef Full Text | Google Scholar

Clark, R. M. (2010). Genome-wide association studies coming of age in rice. Nat. Genet. 42, 926–927. doi: 10.1038/Ng1110-926

PubMed Abstract | CrossRef Full Text | Google Scholar

Cockram, J., White, J., Zuluaga, D. L., Smith, D., Comadran, J., Macaulay, M., et al. (2010). Genome-wide association mapping to candidate polymorphism resolution in the unsequenced barley genome. Proc. Natl. Acad. Sci. U.S.A. 107, 21611–21616. doi: 10.1073/Pnas.1010179107

PubMed Abstract | CrossRef Full Text | Google Scholar

El-Lithy, M. E., Bentsink, L., Hanhart, C. J., Ruys, G. J., Rovito, D., Broekhof, J. L. M., et al. (2006). New Arabidopsis recombinant inbred line populations genotyped using SNPWave and their use for mapping flowering-time quantitative trait loci. Genetics 172, 1867–1876. doi: 10.1534/genetics.105.050617

PubMed Abstract | CrossRef Full Text | Google Scholar

Flint, J., and Eskin, E. (2012). Genome-wide association studies in mice. Nat. Rev. Genet. 13, 807–817. doi: 10.1038/Nrg3335

PubMed Abstract | CrossRef Full Text | Google Scholar

Ganal, M. W., Altmann, T., and Roder, M. S. (2009). SNP identification in crop plants. Curr. Opin. Plant Biol. 12, 211–217. doi: 10.1016/j.pbi.2008.12.009

PubMed Abstract | CrossRef Full Text | Google Scholar

Hemani, G., Knott, S., and Haley, C. (2013). An evolutionary perspective on epistasis and the missing heritability. PLoS Genet. 9:e1003295. doi: 10.1371/Journal.Pgen.1003295

PubMed Abstract | CrossRef Full Text | Google Scholar

Huang, A., Xu, S., and Cai, X. (2014). Whole-genome quantitative trait locus mapping reveals major role of epistasis on yield of rice. PLoS ONE 9:e87330. doi: 10.1371/journal.pone.0087330

PubMed Abstract | CrossRef Full Text | Google Scholar

Ingvarsson, P. K., and Street, N. R. (2011). Association genetics of complex traits in plants. New Phytol. 189, 909–922. doi: 10.1111/J.1469-8137.2010.03593.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Chen, B., Xu, K., Gao, G., Yan, G., Qiao, J., et al. (2016). A genome-wide association study of plant height and primary branch number in rapeseed (Brassica napus). Plant Sci. 242, 169–177. doi: 10.1016/j.plantsci.2015.05.012

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, F., Chen, B., Xu, K., Wu, J., Song, W., Bancroft, I., et al. (2014). Genome-wide association study dissects the genetic architecture of seed weight and seed quality in rapeseed (Brassica napus L.). DNA Res. 21, 355–367. doi: 10.1093/dnares/dsu002

PubMed Abstract | CrossRef Full Text | Google Scholar

Li, H., Peng, Z. Y., Yang, X. H., Wang, W. D., Fu, J. J., Wang, J. H., et al. (2013). Genome-wide association study dissects the genetic architecture of oil biosynthesis in maize kernels. Nat. Genet. 45, 43–50. doi: 10.1038/Ng.2484

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, J., Wang, W., Mei, D., Wang, H., Fu, L., Liu, D., et al. (2016). Characterizing Variation of Branch Angle and Genome-Wide Association Mapping in Rapeseed (Brassica napus L.). Front. Plant Sci. 7:21. doi: 10.3389/fpls.2016.00021

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, R., Qian, W., and Meng, J. (2002). Association of RFLP markers and biomass heterosis in trigenomic hybrids of oilseed rape (Brassica napus x B-campestris). Theor. Appl. Genet. 105, 1050–1057. doi: 10.1007/S00122-002-1050-X

PubMed Abstract | CrossRef Full Text | Google Scholar

Liu, S., Fan, C., Li, J., Cai, G., Yang, Q., Wu, J., et al. (2016). A genome-wide association study reveals novel elite allelic variations in seed oil content of Brassica napus. Theor. Appl. Genet. 129, 1203–1215. doi: 10.1007/s00122-016-2697-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Manolio, T. A., Collins, F. S., Cox, N. J., Goldstein, D. B., Hindorff, L. A., Hunter, D. J., et al. (2009). Finding the missing heritability of complex diseases. Nature 461, 747–753. doi: 10.1038/Nature08494

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, H., Huang, J., Zheng, Q., Huang, Y., Shao, R., Zhu, L., et al. (2013). Identification of combining ability loci for five yield-related traits in maize using a set of testcrosses with introgression lines. Theor. Appl. Genet. 126, 369–377. doi: 10.1007/s00122-012-1985-5

PubMed Abstract | CrossRef Full Text | Google Scholar

Qi, Z., Pan, J., Han, X., Qi, H., Xin, D., Li, W., et al. (2016). Identification of major QTLs and epistatic interactions for seed protein concentration in soybean under multiple environments based on a high-density map. Mol. Breed. 36:55. doi: 10.1007/s11032-016-0475-x

CrossRef Full Text | Google Scholar

Radoev, M., Becker, H. C., and Ecke, W. (2008). Genetic analysis of heterosis for yield and yield components in rapeseed (Brassica napus L.) by quantitative trait locus mapping. Genetics 179, 1547–1558. doi: 10.1534/genetics.108.089680

PubMed Abstract | CrossRef Full Text | Google Scholar

Resende, M. D. V., Resende, M. F. R., Sansaloni, C. P., Petroli, C. D., Missiaggia, A. A., Aguiar, A. M., et al. (2012). Genomic selection for growth and wood quality in Eucalyptus: capturing the missing heritability and accelerating breeding for complex traits in forest trees. New Phytol. 194, 116–128. doi: 10.1111/J.1469-8137.2011.04038.X

PubMed Abstract | CrossRef Full Text | Google Scholar

Sambandan, D., Carbone, M. A., Anholt, R. R. H., and Mackay, T. F. C. (2008). Phenotypic plasticity and genotype by environment interaction for olfactory behavior in Drosophila melanogaster. Genetics 179, 1079–1088. doi: 10.1534/genetics.108.086769

PubMed Abstract | CrossRef Full Text | Google Scholar

Segura, V., Vilhjalmsson, B. J., Platt, A., Korte, A., Seren, U., Long, Q., et al. (2012). An efficient multi-locus mixed-model approach for genome-wide association studies in structured populations. Nat. Genet. 44, 825–830. doi: 10.1038/Ng.2314

PubMed Abstract | CrossRef Full Text | Google Scholar

Shi, J. Q., Li, R. Y., Zou, J., Long, Y., and Meng, J. L. (2011). A dynamic and complex network regulates the heterosis of yield-correlated traits in rapeseed (Brassica napus L.). PLoS ONE 6:e21645. doi: 10.1371/journal.pone.0021645

PubMed Abstract | CrossRef Full Text | Google Scholar

Uwatoko, N., Onishi, A., Ikeda, Y., Kontani, M., Sasaki, A., Matsubara, K., et al. (2008). Epistasis among the three major flowering time genes in rice: coordinate changes of photoperiod sensitivity, basic vegetative growth and optimum photoperiod. Euphytica 163, 167–175. doi: 10.1007/s10681-007-9584-2

CrossRef Full Text | Google Scholar

Wang, F., and Guan, C.-Y., (2010). Molecular mapping and identification of quantitative trait loci for yield components in rapeseed (Brasscia napus L.). Yi Chuan 32, 271–277. doi: 10.3724/SP.J.1005.2010.00271

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, N., Chen, B., Xu, K., Gao, G., Li, F., Qiao, J., et al. (2016). Association mapping of flowering time QTLs and insight into their contributions to rapeseed growth habits. Front. Plant Sci. 7:338. doi: 10.3389/fpls.2016.00338

PubMed Abstract | CrossRef Full Text | Google Scholar

Wang, Z., Cheng, J., Chen, Z., Huang, J., Bao, Y., Wang, J., et al. (2012). Identification of QTLs with main, epistatic and QTL x environment interaction effects for salt tolerance in rice seedlings under different salinity conditions. Theor. Appl. Genet. 125, 807–815. doi: 10.1007/s00122-012-1873-z

PubMed Abstract | CrossRef Full Text | Google Scholar

Xu, L., Hu, K., Zhang, Z., Guan, C., Chen, S., Hua, W., et al. (2016). Genome-wide association study reveals the genetic architecture of flowering time in rapeseed (Brassica napus L.). DNA Res. 23, 43–52. doi: 10.1093/dnares/dsv035

PubMed Abstract | CrossRef Full Text | Google Scholar

Yan, J., Kandianis, C. B., Harjes, C. E., Bai, L., Kim, E. H., Yang, X., et al. (2010). Rare genetic variation at Zea mays crtRB1 increases beta-carotene in maize grain. Nat. Genet. 42, 322–327. doi: 10.1038/ng.551

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., and Zhu, J. (2005). Methods for predicting superior genotypes under multiple environments based on QTL effects. Theor. Appl. Genet. 110, 1268–1274. doi: 10.1007/S00122-005-1963-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J., Zhu, J., and Williams, R. W. (2007). Mapping the genetic architecture of complex traits in experimental populations. Bioinformatics 23, 1527–1536. doi: 10.1093/Bioinformatics/Btm143

PubMed Abstract | CrossRef Full Text | Google Scholar

Yang, J. A., Benyamin, B., McEvoy, B. P., Gordon, S., Henders, A. K., Nyholt, D. R., et al. (2010). Common SNPs explain a large proportion of the heritability for human height. Nat. Genet. 42, 565–569. doi: 10.1038/Ng.608

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, F.-T., Zhu, Z.-H., Tong, X.-R., Zhu, Z.-X., Qi, T., and Zhu, J. (2015). Mixed linear model approaches of association mapping for complex traits based on omics variants. Sci. Rep. 5:10298. doi: 10.1038/srep10298

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, K., Tian, J., Zhao, L., and Wang, S. (2008). Mapping QTLs with epistatic effects and QTL × environment interactions for plant height using a doubled haploid population in cultivated wheat. J. Genet. Genomics 35, 119–127. doi: 10.1016/s1673-8527(08)60017-x

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhang, Z. W., Ersoz, E., Lai, C. Q., Todhunter, R. J., Tiwari, H. K., Gore, M. A., et al. (2010). Mixed linear model approach adapted for genome-wide association studies. Nat. Genet. 42, 355–360. doi: 10.1038/Ng.546

PubMed Abstract | CrossRef Full Text | Google Scholar

Zhao, J., and Meng, J. (2003). Genetic analysis of loci associated with partial resistance to Sclerotinia sclerotiorum in rapeseed (Brassica napus L.). Theor. Appl. Genet. 106, 759–764. doi: 10.1007/s00122-002-1171-2

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuk, O., Hechter, E., Sunyaev, S. R., and Lander, E. S. (2012). The mystery of missing heritability: genetic interactions create phantom heritability. Proc. Natl. Acad. Sci. U.S.A. 109, 1193–1198. doi: 10.1073/Pnas.1119675109

PubMed Abstract | CrossRef Full Text | Google Scholar

Zuk, O., Schaffner, S. F., Samocha, K., Do, R., Hechter, E., Kathiresan, S., et al. (2014). Searching for missing heritability: designing rare variant association studies. Proc. Natl. Acad. Sci. U.S.A. 111, E455–E464. doi: 10.1073/Pnas.1322563111

PubMed Abstract | CrossRef Full Text | Google Scholar