Rice is one of the most important staple foods worldwide. Yearly global rice consumption is estimated to reach 500 million tonnes by 2050 (Abdullah et al. 2008). Arable land area per capita has declined substantially in recent decades, from 0.37 ha in 1961 to 0.196 ha in 2014 (http://data.worldbank.org/indicator/%20AG.LND.ARBL.HA.PC), and supplying sufficient food to a growing population is a major challenge for future human development. Utilization of high yield potential traits will be needed to facilitate the development of high yield rice production systems to overcome the global rice deficit. One of the most effective approaches was used in China, where hybrid rice programs utilizing heterosis contributed to increased yields of up to 13.5 tonnes/ha in 2012 (FAORAP and APSA 2014).

Heterosis, also known as hybrid vigor, is observed in plants, animals, and fungi and is characterized by the superior phenotypic performance of a hybrid compared with its parents. Hybrid vigor is observed in a range of environmental conditions, and phenotypic performance is stably observed. Despite substantial research on heterosis, the genetic underpinnings of the phenomenon are not yet fully understood. Dominance, overdominance, and epistasis (or digenic interactions) are three major genetic components involved in heterosis (Li et al. 2001; Luo et al. 2001; Li et al. 2008; Shen et al. 2014). Dominance and overdominance are interactions of different alleles at the same locus, and epistasis describes interactions between different loci.

Currently available commercial rice hybrids utilizing heterosis optimize the yield potential of crosses between indica varieties. Several recent studies investigated the potential of indica-japonica hybrids (Chin et al. 2011; Xin et al. 2011; Chu et al. 2012; Dan et al. 2014); however, severe hybrid sterility was observed, creating a major barrier to developing super high-yielding hybrid rice (Chu et al. 2012; Dan et al. 2014). The relationship between hybrid sterility and heterosis is unclear, but yield decreases in F₁ plants of inter-subspecific crosses were attributed to severely low seed set. A method of yield adjustment, namely, presupposition yield, was suggested as a means of experimentally minimizing the effect of hybrid sterility (Chu et al. 2012).

Broad-sense heritability of yield-related traits was 21.47–79.85% for seven yield-related traits of testcross populations derived from indica-japonica chromosome segment substitution lines (Tao et al. 2016). The contribution of additive effects (narrow sense heritability) to heterosis was particularly large for grain yield and thousand grain weight (TGW). Dominant effects for spikelet and panicle number also had a major impact on heterosis.

In this study, two recombinant inbred lines (RILs) derived from indica (Korean Tongil type) / japonica crosses were used to investigate the genetic contributions of dominance, overdominance, and epistasis to yield and yield-related traits. Major quantitative trait loci (MQTLs) and EpQTLs (epistatic QTLs or digenic interactions) were identified for each RIL. Comparisons of QTLs/genes in two RILs allowed identification of putative heterosis QTL hotspots. MQTL analysis used average values for the three heterotic allele types (AiAi/AiAj, AiAj/AiAi, and AiAj/AiAj) present in a BC₁F₁ population. The genetic effects of each allele type were then analyzed separately to understand which heterotic allele types were the major contributors to yield and yield-related traits. Lastly, the direct effects of yield-related traits on yield heterosis QTLs were determined.

Plant materials

A total of 179 F₈ RILs derived using the single-seed descendant method from a cross between Dasanbyeo (Tongil type, indica) and TR22183 (temperate japonica) developed by Seoul National University in Korea (DT-RILs) were used in this study (Chin et al. 2011; Jiang et al. 2011). Each RIL was crossed with Dasanbyeo as the female parent, and 178 BC₁F₁ lines (except for line 40) were produced (BC₁F₁ line). F₁ plants of Dasanbyeo and TR22183 were generated and cultivated alongside the BC₁F₁ line to allow identification of specific heterosis traits.

Experiment design and phenotyping

All plant materials were grown at Seoul National University Experimental Farm. Field preparation and crop cultivation were as described previously (Chu et al. 2012). Nine traits were measured and scored: culm length (CL, cm), days to heading (DTH, days), panicle number per plant (PN), panicle length (PL, cm), spikelet number per panicle (SN), spikelet fertility (SF, %), TGW (g), grain yield per plant (GYP), and presupposition yield per plant (PYP). Of these, PN, SN, SF, and TGW were regarded as yield-component traits. To remove the effect of hybrid sterility, PYP was estimated and applied where substantial hybrid sterility was observed in BC₁F₁ lines (Chu et al. 2012). PYP was calculated if the specific BC₁F₁ line exhibited fertility <90%. As SF had a proportional effect on GYP, an adjustment for grain yield value was used (Supplementary Fig. S1). The equation used was PYP = GYP × 90 / SF. Phenotyping results of RILs and BC₁F₁s and calculated mid-parent heterosis (MPH) values are described as the ‘RIL’, ‘BCF1’, and ‘MPH’ datasets, respectively.

Linkage map construction and QTL mapping analysis

Genomic DNA was extracted from 2 cm fresh leaves from seedlings at the three-leaf stage using cetyl-trimethylammonium-bromide (CTAB) as previously described (Rogers and Bendich 1988). Across the 12 chromosomes, 162 SSR and subspecies-specific STS markers (Chin et al. 2007) that showed polymorphisms between Dasanbyeo and TR22183 were used to construct a linkage map with JoinMap version 3.0 software (Stam 1993). Only markers with LOD scores >3.0 were used.

MPH values were calculated for each RIL using the following equation, as given in Chu et al. (2012): (value of MPH) = (average value of the corresponding BC₁F₁) – (average value of RIL and Dasanbyeo). All datasets were analyzed using QTLmapper version 1.6 for the identification of MQTLs and QTLmapper version 2.0 for EpQTLs (Wang et al. 1999). For MQTL analysis, interval mapping using a mixed model approach was applied. For EpQTLs (digenic interactions), only significant two-locus interactions below P < 0.001 were used for further analysis.

Raw genotyping and phenotyping data from another set of indica-japonica RILs, MT-RILs (Milyang 23 × Tong 88-7) (Chu et al. 2012), were re-analyzed using the same software and options to allow comparative analysis with DT-RILs.

Determination and detection of heterosis QTLs

The genetic components of heterosis were suggested to be overdominance, dominance, and epistasis (Li et al. 2008). Different types of dominance (overdominance, partial-to-complete dominance, and underdominance) can be determined by estimating additive effect (a) and dominance effect (d) in RIL and MPH datasets, respectively (Supplementary Fig. S1). Theoretically, sum of dominance and additive effects (a + d) estimated from hetero alleles of linked loci can act as another indicator of heterosis QTLs in the BCF1 dataset. The QTLs detected in both the BCF1 and MPH datasets were regarded as candidates for heterosis QTLs. QTLs detected only in RILs and/or the BC₁F₁ population were classified as additive QTLs (‘A’ type). Dominant QTLs (‘D’ type) and overdominant QTLs (‘OD’ type) were determined using the values of d/a and a + d. QTLs with 0 < d/a < 1 or (a + d) ≥ 2d were designated as ‘D’ type QTLs. QTLs with d/a > 1, or 2d > a + d, or those detected only in MPH datasets, were designated as ‘OD’ type QTLs (Mei et al. 2005; Chu et al. 2012). QTLs detected in MPH datasets but showing negative ‘d’ values were defined as underdominance QTLs (‘UD’ type). Of the detected QTLs, only ‘D’, ‘OD’, and ‘UD’ type QTLs were used for subsequent comparative studies and direct effect analysis.

The allele types that contributed to heterosis in the three heterotic combinations (homo/hetero, hetero/homo, and hetero/hetero) were identified for MQTLs for GYP and PYP. The additive values of yield and yield-related traits (PN, SN, SF, and TGW) of the three heterotic allele types in BC₁F₁ populations were calculated and compared. All statistical analyses were conducted using Minitab version 14.

Comparative mapping study and direct effect analysis

Re-analyzed QTL results and the map of MT-RILs reported by Chu et al. (2012) were compared with those of DT-RILs. Primer sequences of linked markers for all the identified heterosis QTLs were searched using GRAMENE Blast tools (http://ensembl.gramene.org, December 15, 2015) to identify their physical positions in the Nipponbare genome. QTLs anchored in those regions were compared with the heterosis QTLs from this study using the QTARO database (http://qtaro.abr.affrc.go.jp/).

SPSS-AMOS version 21 was used to estimate the standardized direct and indirect effects of yield-related traits (PN, SN, SF, and TGW). The contribution fractions of MQTLs and EpQTLs for each trait were calculated and analyzed together with the standardized direct effects of the yield-related trait components. For each GYP and PYP QTL in two RILs, the traits with the strongest positive effects were proposed if the standardized direct effect value was >0.5.

Trait performance with respect to yield and yield-related traits

Significant differences were observed between Dasanbyeo (P1) and TR22183 (P2) for four traits (DTH, CL, PN, and SN), with no significant differences seen for the remaining traits (PL, SF, TGW, GYP, and PYP) (Table 1). TR22183 exhibited earlier heading (88 days after transplanting) than Dasanbyeo but maintained the same productivity. Values of CL, PL, SN, TGW, GYP, and PYP were significantly elevated in F₁ plants, but very low fertility was observed. In DT-RILs, transgressive segregation was observed in CL, PL, SN, SF, TGW, GYP, and PYP. With the exception of PN and SF, the trait value distributions differed between the BC₁F₁ and RIL populations (Fig. 1). The median panicle trait values were smaller for BC₁F₁ than for F₁ plants in SN, PL, TGW, and PYP. After fertility recovery by backcrossing, median SF and GYP values were larger in BC₁F₁ than in F₁ plants.

MQTL identification

A total of 66 QTLs (15 QTLs for yield traits and 51 QTLs for the other traits) were identified in the RIL, BCF1, or MPH datasets for the traits used in this study (Table 2). Of these, 19 and 14 QTLs were identified in the BCF1 and MPH datasets, respectively.

Culm length (CL)

Eight QTLs in RIL, four QTLs in BCF1, and one QTL in MPH explained 81.4%, 36.4%, and 21.6% of total variance, respectively. No common QTLs were identified across all datasets, but two QTLs of BCF1 (qCL1.5 and qCL4.3) were adjacent to RIL QTLs (qCL1.4 and qCL4.2). One overdominance (‘OD’ type) QTL, which had the effect of extending CL by 2.11 cm, was detected on chromosome 9 (qCL9.1). All other QTLs were classified as additive (‘A’ type) QTLs.

Days to heading (DTH)

Four QTLs in RIL and one QTL in BCF1 explained 42.3% and 12.4% of total variance, respectively. No QTLs showing the effects of heterotic alleles were identified.

Panicle number per plant (PN)

Three QTLs in RIL and one QTL in BCF1 explained 23.5% and 16.3% of total variance, respectively. No heterosis QTLs were identified for this trait.

Panicle length (PL)

Three QTLs in RIL, two QTLs in BCF1, and two QTLs in MPH were responsible for 59.1%, 16.4%, and 24.8% of total variance, respectively. Two heterosis QTLs were identified in MPH. One dominance (‘D’ type) QTL with relatively large additive (36.8% phenotypic variance explained; PVE) and dominance effects (15.4% PVE) was identified in RIL on chromosome 8 (qPL8.1). Another overdominance QTL (qPL4.1), which had a smaller dominance effect (9.4% PVE) than qPL8.1, was identified in MPH.

Spikelet number per panicle (SN)

Five QTLs in RIL, two QTLs in BCF1, and one QTL in MPH explained 43.9%, 17.9%, and 10.2% of total phenotypic variance, respectively. One overdominance heterosis QTL (qSN9.1), which resulted in a 6.2% increase in SN, was identified on chromosome 9 in MPH.

Spikelet fertility (SF)

Two QTLs in RIL and two QTLs in BCF1 were responsible for 13.7% and 27.6% of total phenotypic variance, respectively. For SF, examination of hybrid sterility in BCF1 is more meaningful than consideration of heterosis QTL. Two QTLs in BCF1 (qSF4.1 and qSF6.1) contributed to 7.86% and 19.69% reductions in SF, respectively.

Thousand grain weight (TGW)

Eight QTLs in RIL, two QTLs in BCF1, and six QTLs in MPH explained 50.4%, 22.4%, and 43.9% of total phenotypic variance, respectively. Three overdominance heterosis QTLs (qTGW3.3, qTGW4.1, and qTGW10.1) and two dominance heterosis QTLs (qTGW1.1 and qTGW7.1) were identified. One underdominance QTL (qTGW4.2) with a negative dominance effect (−0.41) was identified on chromosome 4, adjacent to qSF4.1. Excepting qTGW4.2, the total PVE of the five heterosis QTLs was 41.7%.

Grain yield per plant (GYP)

Three QTLs in RIL, three QTLs in BCF1, and one QTL in MPH were responsible for 30.7%, 34.4%, and 16.0% of total phenotypic variance, respectively. No heterosis QTLs were identified for GYP in DT-RIL. One identified QTL on chromosome 5 (qGYP5.1) showed a negative dominance effect (−2.72%).

Presupposition yield per plant (PYP)

Three QTLs in RIL, two QTLs in BCF1, and three QTLs in MPH explained 27.3%, 15.1%, and 23.4% of total phenotypic variance. Two QTLs (qPYP8.1 and qPYP9.1) were also QTLs for GYP (qGYP8.1 and qGYP9.1). One QTL for GYP on chromosome 11 (qGYP11.1), which showed a negative additive effect (−0.9), was not detected in the PYP dataset, as expected. One QTL on chromosome 2 (qPYP2.1) was not found for GYP. Three overdominance heterosis QTLs were identified on chromosomes 1, 7, and 8. The removal of hybrid sterility effects allowed possible QTLs for GYP to be identified.

The QTLs identified only in the BCF1 population were not classified as heterosis QTLs in this study. However, many of these QTLs showed positive values of combined additive and dominant effects: all four QTLs for CL; one QTL for DTH; one QTL for PN; one of two QTLs for PL (qPL9.1); one of two QTLs for SN (qSN4.2); one of two QTLs for TGW (qTGW1.2), and both QTLs for PYP.

Many of the QTLs were clustered within specific regions on chromosome 1 (RM297-S01143A), chromosome 4 (S04120-S04129B), and chromosome 8 (RM25-RM72) (Fig. 2, left chromosome shown).

EpQTL identification

When all traits were considered, 23 EpQTLs in RIL, 15 EpQTLs in BCF1, and two EpQTLs in MPH were identified in this study (Table 3). The total PVE of EpQTLs for each trait was 10.7–27.5% for RIL, 11.2–26.8% for BCF1, and 7.9–22.0% for MPH. No EpQTLs were identified for SN and PYP in RIL or for PN and TGW in BCF1.

For grain yield traits, different EpQTLs for GYP were identified in RIL and BCF1. Only one EpQTL, between chromosome 1 (S01143A-S01157B, qSF1.1 of RIL) and chromosome 8 (RM72-S08052B, qGYP8.1 of RIL), contributed positively to higher grain yield (1.75 g/plant). For PYP, one EpQTL between chromosome 3 (S03010B-S03020) and chromosome 10 (S10071-RM590, qSN10.1 of RIL) had positive additive effects.

One EpQTL for PL and one for SN were identified in the MPH dataset. No MQTLs were identified within the chromosomal regions linked with EpQTLs for PL. However, for SN, qPL12.1 in RIL interacted with the region on chromosome 4.

Trait contribution to heterosis QTLs for grain yield

First, to identify the yield-component traits contributing substantially to the three overdominance heterosis QTLs and one underdominance QTL for grain yield in DT-RIL, the additive effects of each yield-component trait on the QTLs were presented in radar charts (Fig. 3). Interval mapping through backcrossing was used to identify QTLs between markers. There were three types of heterotic allele combination for each marker (AiAi/AiAj, AiAj/AiAi, or AiAj/AiAj), and the additive effects of each allele type were calculated by dividing the value of the corresponding allele type with that of the AiAi/AiAi (homo) allele genotype. Additive effects identified in this study were in the −15–15% range (Table 4).

For qPYP1.1, the AiAi/AiAj (S01140/S01143A) genotypes did not contribute to heterosis. However, AiAj/AiAi genotypes had a significantly larger PN, and significantly or highly significantly smaller SN, SF, and GYP, than homoallele genotypes. Likewise, AiAj/AiAj genotypes had significantly higher PN and PYP values than homoallele genotypes. This suggested that the PN of AiAj/AiAj genotypes might be the major factor explaining the significant heterotic effect on PYP.

For qPYP7.1, AiAi/AiAj (RM234/S07101) genotypes showed negative additive effects on PYP. The AiAj/AiAi and AiAj/AiAj genotypes had significantly higher PYP values than homoallele genotypes. The additive effect of AiAj/AiAi on SN might explain the significantly positive additive effect of PYP. Although a significant positive additive effect of SF was found for AiAi/AiAj, the PYP values of the genotypes were negative, thus demonstrating that SF alone could not be a factor of heterosis.

For qPYP8.2, AiAi/AiAj (S08090/S08106) and AiAj/AiAj showed significant additive effects for PYP, which may have been caused by low hybrid fertility. AiAj/AiAj genotypes showed a positive effect of SF, but PYP did not show additive effects.

Unlike the PYP cases, no heterosis QTLs were identified in GYP. However, one underdominance GYP QTL (qGYP5.1) was analyzed to see which allele types were contributory. Only AiAj/AiAj (S05014B/RM267) showed a significant negative effect on SF and GYP. Two other heterotic allele genotypes of linked markers (AiAj/AiAi and AiAi/AiAj) contributed to significantly higher SF and GYP values. This result indicates the importance of studying different heterotic allele combinations of linked markers for underdominance QTLs.

Second, the direct effects of each trait on grain yield were analyzed using the path analysis module of SPSS-AMOS software (Table 5). Only DTH, PN, SN, SF, and TGW were considered as yield-related traits for this analysis, because CL was a vegetative trait and PL was highly correlated with SN. In each heterosis QTL for grain yield, the standardized direct effects of each trait >0.5 or <−0.5 were considered as positive effect QTLs. As expected from the radial chart analysis, PN and SN were the most important traits explaining heterosis of qPYP1.1. Likewise, PN, SN, and TGW variously contributed to qPYP7.1 for each allele combination. qPYP8.2 was related to PN and SN for the AiAj/AiAi type only. PN, SN, and SF showed direct effects on GYP via the contribution of the AiAj/AiAi and AiAi/AiAj genotypes at qGYP5.1. SF was a negatively effective trait in all heterosis QTLs for PYP.

Heterosis QTLs of intra- and inter-locus interactions

Heterosis is governed by various genetic effects including overdominance effects, dominance effects, additive effects, and their interactions. In this study, indica-japonica RILs were used to identify heterosis QTLs for grain yield and its related traits. Theoretically, heterotic effects (inter-locus heterosis) should be explained solely by the genetic interactions between different recombinant blocks originating from different parents. However, two-QTL interactions (or digenic interactions) are difficult for plant breeders to utilize due to their complexity and the difficulties in dealing with the large numbers of samples in breeding programs. Accessing heterosis through intralocus interactions between different alleles (dominance and overdominance) may thus be more amenable. Intra-and inter-locus heterosis can be assessed using BC₁F₁ populations, which were used to study overdominance and dominance effects in this study. Inter-locus heterogeneity is more prevalent in RIL rather than BC₁F₁ populations, and estimation of intra-locus heterosis may be improved using BCF1 and MPH datasets. The heterotic effects on inter-locus and intra-locus interactions governing grain yield may not be identified at the same loci in different populations. For example, two EpQTLs for GYP in RIL were unidentified in the BCF1 and MPH datasets (Tables 2 and 3). On the other hand, two EpQTLs for GYP and two EpQTLs for PYP in MT-RILs (RILs derived from Milyang 23 and Tong 88-7; Chu et al. 2012) were not identified in heterosis MQTL analysis using BCF1 and MPH datasets in DT-RILs (Supplementary Tables S1 and S2).

In DT-RIL, there were nine heterosis MQTLs for all traits and one underdominance heterosis MQTL for GYP (Table 2). EpQTLs for heterosis are rarely detected if the population size is small, because the effects of multi-locus allele combinations require a large number of samples and meaningful statistical analysis for detection. Thus, the subsequent heterosis QTL analyses in this study were conducted using only MQTL results. Nevertheless, there were clearly no EpQTLs for GYP or PYP found using MPH datasets in DT-RIL. Conversely, in MPH datasets of MT-RIL, three overdominance MQTLs for GYP and two for PYP were identified (Supplementary Table S1).

Common heterosis QTL clusters in two indica-japonica RILs indicating relatedness of heterosis for yield and for yield-related traits

MQTLs for different traits were co-located at some chromosomal regions (Fig. 2), and several QTL clusters were observed in MT-RILs, on chromosome 1 (S01143A-S01157B), chromosome 4 (RM348-S04128), chromosome 6 (S06015-S06019B, RM50-RM527, and S06087-S06100), and chromosome 7 (S07035-S07055B). In both RILs, QTL clusters included QTLs identified in RIL, BCF1, and MPH datasets for several yield-related traits. Thus, an additive effect of preferred alleles in RILs would be expected in marker-assisted breeding programs that utilized the heterosis QTLs in the QTL clusters. This may explain the plant growth vigor seen in recombinant plants generated through introgression of donor QTLs in some generations of breeding programs.

Heterosis QTLs for each trait were compared between the two RILs (Supplementary Table S3), and no common heterosis QTLs were found. However, the QTL clusters may represent common loci for heterosis (Fig. 2). QTL clusters in DT-RIL (RM297-S01143A) and MT-RIL (S01143A-S01157B) were at adjacent regions on chromosome 1, and might share common genetic elements for heterosis. This was also observed for S04120-S04129B (CL, SF, and TGW) in DT-RIL and RM348-S04128 (PL, GYP, and SN) in MT-RIL on chromosome 4. Three QTL regions linked to grain yield traits are close in comparative linkage maps and should be studied further to see whether common regions for heterosis are shared between the RILs (Fig. 2): RM584-RM204 in DT-RIL and S06015-S06018 in MT-RIL on chromosome 6, RM253-S06031 (SF in DT-RIL) and RM50-RM527 (PYP in MT-RIL) on chromosome 6, and RM21-S11117 (GYP in DT-RIL) and S11112-S11117B (PYP in MT-RIL) on chromosome 11.

Previously reported QTLs that are linked to the identified heterosis QTLs in DT and MT-RILs and are linked to a range of agronomic traits are listed in Supplementary Table S4. QTL locations were compared between the physical map locations of the markers in this study and the linked markers from various mapping populations. In DT- and MT-RILs, qTGW4.2, qGYP8.1, qSN8.1, and qCL9.1/qSN9.1 were not previously reported as yield or yield-related QTLs. With the exception of these, all the QTL regions in this study were previously linked to yield or yield-related traits. For example, qSN1.1 of MT-RIL was co-located with gp1a in F₂ (grains per panicle, Shen et al. 2003), qPN1 in BC₄F₄ (panicles per plant, Tian et al. 2006), and others. Thus, some previously identified QTLs linked to different traits may be co-located with or the same as QTLs identified in this study. This may explain the significant positive correlations observed between some yield traits. For example, with respect to the linked QTLs given above, grain number (gp1a) is correlated with spikelet number (qSN1.1) and panicle number (qPN1.1).

Previous studies of RILs rarely discussed whether these genes/QTLs exhibited heterosis. Some previously reported QTLs/genes were co-located with the heterosis QTLs identified in this study, and further analysis is needed to determine whether heterosis is seen in heterotic alleles. For DT-RIL, these were as follows: qTGW3.3 (gw3a, grain weight, Yu et al. 1997); qGYP5.1 (S31, hybrid sterility, Li et al. 2005; qSF5.2, hybrid sterility, Reflinur et al. 2012); qPYP7.1 (qPLLN7, pollen fertility, Yan et al. 2003; hd7a, heading date, Yu et al. 2002); qTGW7.1 (S23, hybrid sterility, Sobrizal et al. 2000); qPL8.1 (RH8, heading date, Li et al. 2016; qPL8, Hd5, heading date, Lin et al. 2003; PL, Miyata et al. 2007); and qPYP8.2 (WFP, panicle branching and yield, Miura et al. 2010). For MT-RIL, these were qCL1.2 (sd1, plant height, Li et al. 2003); qGYP6.1 (DPL2, sterility, Mizuta et al. 2010; qSPTF6, SF, Yan et al. 2003); qPYP6.1 (spf6, SF, Song et al. 2005; S5, SF, Yan et al. 2000); qCL7.2 (Hd4, DTH, Lin et al. 2003); and qPYP11.1 (FON3, floral organ number, Jiang et al. 2007). These QTLs/genes were reported in various population types as well as in backcrossed populations, implying that QTLs showing additive effects in homozygous alleles may give heterotic effects.

Comparative analysis indicated that the negative additive values of some QTLs of DT- and MT-RILs could be partially explained by hybrid sterility QTLs/genes. qGYP5.1, qGYP6.1, and qGYP10.1 were co-located with known sterility QTLs (qSF5.2 and S31) and QTLs for pollen fertility (Doi et al. 1998; Tan et al. 1998). However, there were no known hybrid sterility QTLs/genes for qTGW4.2, qSN8.1, qTGW9.1, and qGYP10.1. Conversely, qPYP6.1 in MT-RIL did not show negative additive values, indicating a lack of allele-induced hybrid sterility between Milyang 23 and Tong 88-7.

Tao et al. (2016) also reported some heterosis QTLs from testcrossed CSSLs (Chromosome Segment Substitution Lines). The QTL hotspots on chromosomes 1, 6, and 8 identified in this study were also identified by Tao et al. (2016).

Three heterotic allele genotypes on heterosis QTLs of MT-RILs

For QTL determination, the effects of three heterotic allele combinations (AiAi/AiAj, AiAj/AiAi, and AiAj/AjAj) were combined, hampering the identification of the allele types that had the most substantial impact on phenotypic values (Fig. 3 and Supplementary Fig. S2).

The range of additive values was larger for MT-RILs than for DT-RILs, indicating that the effect of MQTLs was larger for MT-RILs than for DT-RILs. Five QTLs for GYP and two QTLs for PYP were investigated in this study (Supplementary Fig. S2 and Supplementary Table S5). For qGYP1.2 (RM259/RM600), the heterotic allele of RM259 contributed to the additive effect of SF and heterosis of GYP. The heterotic alleles linked to both the markers contributed to the increase of PN. A QTL for panicle number (ppp1.1) was reported on qGYP1.2 (Cho et al. 2007). This was not significant (P = 0.051), but a significant positive effect was seen for heterotic allele type for SF contribution to qGYP5.1. On qPYP6.1 (S06015/S06018), all three heterotic allele types exhibited similar contributions to heterosis of PYP. SN was the most important contributory trait for PYP increase, though SF was significantly low for all three types. A QTL for SN was identified on qPYP6.1 (Zhuang et al. 2001). At qGYP8.1 (S08059/S08066), heterotic allele effects of S08059 and S08066 gave contrasting effects for SF and GYP. S08059 heterotic allele type positively impacted the additive effect of GYP. On qPYP11.1 (S11112/S11117B), the AiAi/AiAj and AiAj/AiAj types showed significantly low SF but high PYP. The presupposition yield calculation was influenced by S11117B heterotic allele type. FON3 (Jiang et al. 2007) reported a QTL for floral organ numbers in the F₂ population on the same region. However, it was unclear which trait contributed to the PYP increase of qPYP11.1 in this population. Low SF on qGYP6.1 (RM50/RM527) was impacted by RM527 heterotic allele type, giving significant positive additive effects of PYP and significant negative effects of GYP. Likewise, qGYP10.1 (S10019/S10026C) also showed an underdominance heterosis type.

Direct effect analysis for MQTLs and EpQTLs

It is desirable to identify the yield-related traits or yield-component traits that are most influential on yield. In this study, to simplify the complicated mediation networks, only the standardized direct effect of each yield-related trait on PYP and GYP was estimated by path analysis (Table 5).

The direct effect of each trait was estimated for each different allele type. One homoallele type (AiAi/AiAi) in RIL and three heterotic allele types (AiAi/AiAj, AiAj/AiAi, and AiAj/AiAj) in BCF1 populations were analyzed separately. PN and SN were the yield-related traits with most direct impact on heterosis QTLs in DT-RILs. For qPYP1.1, the strongest trait, with an estimated direct effect >0.5, was SN, particularly for AiAj/AiAi and AiAj/AiAj types (62 of 179 RILs; 34.6% (Table 4). However, for the same locus, PN was the most effective trait for qPYP1.1 with respect to the homoallele (AiAi/AiAi) and complete heterotic allele (AiAj/AiAj) types. The same approach could be applied to other heterosis yield QTLs. Different combinations of yield-related traits provided direct effects for each qPYP7.1 allele type. Only the AiAj/AiAi allele type of qPYP8.2 was linked with PN and SN in BCF1. Conversely, for qGYP5.1, only the AiAi/AiAj type was directly impacted by PN, without the negative effect of SF.

A similar approach was used for MT-RILs (Supplementary Table S6). Unlike with DT-RILs, the positive direct effect on SF was apparent across many QTLs. A direct negative effect of SF was found only for qPYP6.1 in all allele types. qGYP1.2, qGYP5.1, and qGYP10.1 all positively impacted SF. Conversely, SN had the most positive effect in AiAj/AiAj of qGYP1.2, AiAj/AiAj of qPYP6.1, AiAj/AiAi of qGYP6.1, and AiAj/AiAj of qGYP10.1. TGW was influenced by AiAi/AiAj of qPYP6.1 and AiAj/AiAi of qGYP6.1. PN and DTH directly affected qPYP6.1 in partial heterotic allele types (AiAi/AiAj or AiAj/AiAi).

In brief, the standardized direct effects of different heterotic allele types on yield-related traits differed between two indica-japonica RILs. Heterosis effects were identified using the accumulated additive effects of PN, SN, SF, and TGW (Dan et al. 2015). For each QTL, the yield-related trait with the most positive impact could be identified by this analysis. Followed by the development of near-isogenic lines (NILs) for each QTL, the most important yield-related traits could next be studied to estimate the consequent yield heterosis.

The contribution of MQTLs and EpQTLs to heterosis can be estimated by the sum of total PVEs for the identified QTLs. Fig. 4 shows the total effects of PVE of each yield-related trait (center), PYP (outer sides), and GYP (inner sides) connected by arrows representing the direct effects. The total PVEs of yield-related traits were larger in DT-RILs than in the BCF1 population. However, total PVEs of GYP and PYP were larger in the BCF1 population than in RILs. In the DT-BCF1 population, contributions of EpQTLs for GYP (61.2%) and PYP (36.4%) were large, and consequently, the heritability (total PVE) of yield heterosis was larger than GYP (45.8%) and PYP (27.3%) in RILs. This differed for MT-RILs (Supplementary Fig. S3), where heritability of yield heterosis was smaller than in RILs. This implies that the contribution of the dominance effect was larger than that of epistasis between Dasanbyeo (indica) and TR22183 (japonica).

Suggestions for plant breeding strategies utilizing heterosis QTLs

Heterosis QTLs identified by current software provide average values for different allele types of heterotic alleles. The additive values of different heterotic allele combinations could be calculated from the identified QTLs to identify the allele type that contributed to the preferable yield heterosis. Hybrid sterility, which is common in indica-japonica crosses in rice, hampers the interpretation of heterotic effects of yield traits in specific heterotic allele types at those loci. The breakdown of allele types may help resolve this.

The heterosis QTLs identified in the two RILs in this study mostly coincided with previously reported yield and yield-related traits. This similarity was also noted when specifically expressed genes of F₁ plants were compared to the linkage map of Nipponbare/93-11 (Venu et al. 2014). This also showed that heterosis could be explained by certain specific yield-related traits associated with specific loci. These loci might be detected as MQTLs in conventional QTL studies.

The identified EpQTLs explained the presence of epistatic interactions between loci for heterosis. To utilize EpQTLs for plant breeding, marker-assisted recurrent selection or genomic selection strategies would be more effective rather than the simple marker-assisted backcrossing strategies used for MQTLs.

Common heterosis QTLs identified across different populations and studies would be highly valuable for heterosis breeding programs; however, further quantitation of the complex genetic networks for each yield-related trait is needed to design breeding programs with high predictability.

This study was supported by a grant from the Next-Generation BioGreen 21 Program (no. PJ01102401) of the Rural Development Administration, Korea. This research was also supported by the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (no. 2017R1D1A1B04034862).