Powerful one-dimensional scan to detect heterotic QTL

Powerful one-dimensional scan to detect heterotic QTL | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results Powerful one-dimensional scan to detect heterotic QTL View ORCID Profile Guoliang Li , View ORCID Profile Renate H. Schmidt , View ORCID Profile Yusheng Zhao , View ORCID Profile Jochen C. Reif , View ORCID Profile Yong Jiang doi: https://doi.org/10.1101/2025.05.06.652159 Guoliang Li 1 Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) Gatersleben , Corrensstraße 3, D-06466 Seeland, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Guoliang Li Renate H. Schmidt 1 Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) Gatersleben , Corrensstraße 3, D-06466 Seeland, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Renate H. Schmidt Yusheng Zhao 1 Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) Gatersleben , Corrensstraße 3, D-06466 Seeland, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yusheng Zhao Jochen C. Reif 1 Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) Gatersleben , Corrensstraße 3, D-06466 Seeland, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Jochen C. Reif For correspondence: jiang{at}ipk-gatersleben.de reif{at}ipk-gatersleben.de Yong Jiang 1 Leibniz Institute of Plant Genetics and Crop Plant Research (IPK) Gatersleben , Corrensstraße 3, D-06466 Seeland, Germany Find this author on Google Scholar Find this author on PubMed Search for this author on this site ORCID record for Yong Jiang For correspondence: jiang{at}ipk-gatersleben.de reif{at}ipk-gatersleben.de Abstract Full Text Info/History Metrics Supplementary material Data/Code Preview PDF Abstract To meet the growing global demand for food, increasing yields through heterosis in agriculture is crucial. A deep understanding of the genetic basis of heterosis has led to the development of a quantitative genetic framework that incorporates both dominance and epistatic effects. However, incorporating all pairwise epistatic interactions is computationally challenging due to the large sequencing depth and population sizes needed to uncover the genes behind complex traits. In this study, we developed hQTL-ODS, a one-dimensional scanning method that directly assesses the net contribution of each locus to heterosis. Simulations show that hQTL-ODS reduces computational time while offering higher power and lower false-positive rate. We applied this method to a population of 5,234 wheat hybrids with whole-genome sequenced profile, revealing key epistatic hubs that play a critical role in determining heterosis. Our findings offered valuable insights for improving breeding strategies to enhance crop yields. Main Heterosis, the enhanced performance of offspring compared to their parents, is a phenomenon whose exploitation contributes significantly to food security 1 , 2 . It is systematically leveraged in both livestock 3 – 5 and crop plants 6 – 8 through hybrid breeding. Heterosis can be explained genetically either by the overdominant effect of individual genes 9 – 11 , by the complementary combination of several genes with (partially) dominant gene action 12 – 14 or by specific epistatic interactions between genes 15 – 17 . The explanatory hypotheses are not mutually exclusive. Interestingly, the causes of heterosis in most crop and livestock species are poorly understood, but a proper theoretical framework and declining sequencing costs pave the way to fill these gaps through genetic association studies. The heterotic effect of a gene or quantitative trait locus (QTL) can be attributed to its dominance effect and its epistatic interaction effects with other loci 18 , 19 . Based on the quantitative genetic description of the heterotic effect, a multi-step scanning procedure 19 was developed to detect heterotic QTL (hQTL) in a hybrid population created by crossing different parents: QTL scans are performed for individual components, i.e., dominance effects of all markers and di-genic epistatic (additive-by-additive, additive-by-dominance, and dominance-by-dominance) effects of all marker pairs, followed by a test that integrates all significant components of a given putative heterotic hQTL (hQTL-MSS). This is the only existing method that considers both dominance and epistatic effects to detect hQTL in a diverse hybrid population. However, the challenge with this approach is that as the number of markers increases to for instance a million, the computational burden becomes prohibitive with trillions of digenetic epistatic interactions to be tested. This is particularly challenging in large populations which are needed to ensure high power for QTL detection. One solution to this problem is to avoid screening all individual components and test directly for the net heterotic effect of a locus, reducing the number of tests to the number of markers. In our study, we developed a novel one-dimensional scanning (hQTL-ODS) method to test the net contribution of individual loci to heterosis ( Fig. 1 ). This allowed us to substantially reduce the computational time of genome-wide searches for heterosis loci in large populations. In simulation studies, we showed that the method has higher statistical power with fewer false positives. Capitalizing on large hybrid wheat populations genotyped with whole-genome resequencing (WGS) data, we tested the strategy and showed the importance of the contribution of cumulative small epistatic interactions to heterosis in wheat. Download figure Open in new tab Figure 1. An overview of hQTL-ODS in comparison to the existing approach hQTL-MSS 19 . Results Methods overview The hQTL-ODS model directly tests the heterotic effect of a marker, avoiding extensive association tests for pairwise interaction effects ( Fig. 1 ). For each marker, the following two linear mixed models (LMMs) are compared. where y MPH is the vector of mid-parent heterosis (MPH) values for all hybrids. Model (1) is the null model containing covariate effects ( Xα ), multiple genetic background effects ( g D , g AA , g AD , g DD ) and the residuals ( ε ), but excluding any marker effects. The alternative model (2) includes the heterotic effect ( h i ) of the tested marker in addition to all effects in the null model. The heterotic effect h i is defined as a complex linear combination of (4 p -3) effects, i.e. the dominance effect d i and the epistatic effects aa ij , ad ij , da ij , dd ij ( j = 1, …, p and j ≠ i , where p is the number of markers) 19 . In the model, h i is treated as a random vector following a multi-variate normal distribution: , where H i is a covariance matrix and is a variance component. A key step to implement the model is to derive an alternative expression of h i which is less complex than (but equivalent to) the original definition. It is also crucial to find an efficient method of calculating H i (see Methods for more details). After solving the model by the restricted maximum likelihood (REML) method, the significance of h i is assessed by testing the null hypothesis using the likelihood ratio test. Namely, the test statistic is of the form LR i = −2In( L 0 / L 1 ), where L 0 is the maximum likelihood of the null model, and L 1 is the maximum likelihood of the alternative model. Under the null hypothesis, this test statistic follows approximately a mixture of and distributions with equal weights 20 , 21 , where and refer to chi-squared distributions with zero and one degree of freedom. hQTL-ODS is computationally more efficient than hQTL-MSS The new model hQTL-ODS has a clear advantage in computational time complexity, compared with hQTL-MSS. In hQTL-MSS, GWAS has to be performed for dominance effects of all markers and for digenic epistatic effects of all marker pairs based on an LMM similar to model (2), replacing h i by these component effects. Thanks to the P3D approximation 22 it is sufficient to solve the null model (1) once. Since the model contains multiple random vectors ( g D , g AA , g AD , g DD ) in addition to the residuals, the computational load cannot be reduced by applying eigen-decomposition to the covariance matrices, a commonly implemented technique for LMMs with a single random vector. Thus, the time complexity of this step is O ( tn 3 ), where n is the number of individuals and t is the number of iterations, which depends on the numerical method used to solve the model and is less than 10 in many cases. Computing the test statistics for dominance and epistatic effects takes O ( n 2 p ) and O ( n 2 p 2 ) time, respectively. In case of WGS data, it is expected that p ≫ tn . Thus, n 2 p 2 ≫ n 2 p ≫ tn 3 , indicating that the time complexity of hQTL-MSS is dominated by the procedure of testing digenic epistatic effects, namely O ( n 2 p 2 ). In hQTL-ODS, it is necessary to solve p + 1 LMMs (the null model once and the alternative model for each marker), resulting in O ( tn 3 p ) time. The time for producing the likelihood ratio test statistics is merely O ( p ), which is negligible. Additional time is needed for hQTL-ODS to calculate the covariance matrix H i for each marker. By implementing an efficient algorithm for calculating genomic epistatic relationship matrices 23 , the time required for calculating H i for all markers can be reduced from O ( n 2 p 2 ) to O ( n 2 p ) (see Methods for details). Hence, the time complexity of hQTL-ODS is dominated by solving p + 1 LMMs, namely O ( tn 3 p ). Considering p ≫ tn , we have n 2 p 2 ≫ tn 3 p . Therefore, the hQTL-ODS model is computationally much more efficient than hQTL-MSS. More precisely, the time required by hQTL-MSS is approximately times as long as that required by hQTL-ODS. For example, if t ≈ 10, hQTL-ODS is about 100 times faster than hQTL-MSS for a data set with 1,000 individuals and one million markers. The theoretical advantage of hQTL-ODS mentioned above was validated by using a wheat data set consisting of 1,557 hybrids with 1.2 million single-nucleotide polymorphisms (SNPs). The data set was termed Exp I (Table S1, see Methods for details). From this data set, we sampled subsets with four marker numbers (5,000, 10,000, 20,000 and 50,000) in combination with four sample sizes (200, 500, 1,000 and 1,557). hQTL-ODS and hQTL-MSS were performed for each subset on a computing platform with 50 CPU cores (Intel(R) Xeon(R) Gold 6130 CPU @ 2.10GHz) and 24 GB memory size per CPU. We observed that the running time of hQTL-ODS ( t ODS ) was substantially shorter than that of hQTL-MSS ( t MSS ), and the relative advantage measured as the ratio ρ = t MSS / t ODS depended on the population size and the number of markers ( Fig. 2a ). When the population size is fixed, the ratio linearly increased relative to the number of markers, which is well in accordance with the theory that ρ approximately equals p / tn ( Fig. 2b-e ). In the entire population with 1,557 genotypes, the observed trend of ρ values almost completely overlapped with the theoretical line of t = 5 ( Fig. 2e ). This means that hQTL-ODS would be about 154 times faster than hQTL-MSS for the full data set (1,557 genotypes with all 1.2 million markers). While it took 180 hours for hQTL-ODS to finish the analysis on the full data set, it would take more than three years for hQTL-MSS to complete the task. Download figure Open in new tab Figure 2. Comparison of running time, power and false-positive rate between hQTL-ODS and hQTL-MSS on experimental and simulated data sets. (a) The running time of hQTL-ODS ( t ODS ) and hQTL-MSS ( t MSS ) on experimental data sets with different sample sizes and marker numbers. The numbers in parentheses represent the ratio ρ = t MSS / t ODS . (b-e) The relative advantage of hQTL-ODS over hQTL-MSS. Theoretically, hQTL-ODS is approximately p / tn times faster than hQTL-MSS, where p is the number of markers, n is the sample size, and t is the average number of iterations for solving the linear mixed model. The lines depicted the expected values p / tn and the blue dots indicated the observed ratio ρ = t MSS / t ODS on experimental data sets. (f) The statistical power and (g) The false-positive rate of hQTL-ODS and hQTL-MSS evaluated on simulated data sets for five different scenarios (s1-s5). The simulated hQTL effect was contributed by epistatic effects (s1 and s2), by both dominance and epistatic effects (s3 and s4), or by dominance effects (s5). More details on the simulation scenarios were described in Methods and Table S2). Statistical power and false positive rate – hQTL-ODS outperforms hQTL-MSS hQTL-ODS is also theoretically superior to hQTL-MSS in terms of properly modeling heterotic effects. In hQTL-ODS, the heterotic effects were modeled unbiasedly according to the original definition, i.e., the dominance effect of a locus and its epistatic interaction effect with all other loci were considered, independent of their effect sizes. In contrast, hQTL-MSS filtered the component effects by testing them and applying a threshold, hence only modeled large effects. When many small component effects cumulatively make a large contribution, the heterotic effect estimated by hQTL-MSS is likely to be biased. To verify the above hypothesis, we compared the performance of hQTL-ODS and hQTL-MSS in a simulation study. MPH values were simulated following model (2) based on the genomic data of a subset of Exp I, consisting of 1,000 hybrids with 5,000 SNPs. Five scenarios were considered in order to take different patterns of component genetic effects composing hQTL into account (Table S2, see Methods for details). In all scenarios, only one hQTL was simulated and the proportion of phenotypic variance explained (PVE) by the hQTL was fixed at 2.5%. Despite the small PVE, the simulated hQTL in all scenarios were adequately captured by hQTL-ODS, with a detection power ranging from 45% to 63%, and a false positive rate (FPR) of less than 0.07% ( Fig. 2f, g ). When a few epistatic effects cumulatively contributed to the effect of the simulated hQTL, no matter whether the dominance effect contributed (Scenarios 3 and 4) or not (Scenarios 1 and 2), hQTL-ODS consistently outperformed hQTL-MSS, with higher power and lower FPR ( Fig. 2f, g ). Only in Scenario 5, where the hQTL effect consisted of the dominance effect of the QTL alone, the power of both methods was the same. These results clearly validated the theoretical superiority of hQTL-ODS over hQTL-MSS. A large resequenced hybrid wheat panel In order to test the performance of hQTL-ODS on experimental data, three hybrid wheat populations, each from independent experimental series (Exp I, II and III) 24 , were assembled into an integrated panel consisting of 5,234 hybrids derived from 597 parental lines representing a comprehensive selection of the Central European elite bread-wheat breeding pool. In each Exp, hybrids and parental lines were evaluated in field trials for grain yield and heading date in at least 10 year-by-location combinations (Table S1, S3). The parental lines were resequenced using WGS technology, and sufficient read coverage was established for 588 genotypes (Table S4). Read mapping and SNP calling against the reference genome assembly of cv. Chinese Spring v2.1 (RefSeq CSv2.1) 25 provided 7,835,467 SNPs with minor allele frequency (MAF) larger than or equal to 0.05 and with an average density of 0.54 SNPs per kb. SNPs whose pairwise linkage disequilibrium (LD, estimated by r 2 ) was higher than 0.9 were pruned within 50kb sliding window and 907,534 high-quality SNPs remained for subsequent analysis. The genotypes of the hybrids were inferred from the corresponding parental lines. Combining the genomic and phenotypic data resulted in 4,885 hybrids derived from 545 parental lines. A principal coordinate analysis (PCoA) revealed that neither the parental lines from each population nor across populations grouped into clusters ( Extended Data Fig. 1a ). Different Exps were linked by up to 26 common genotypes showing high trait correlations between different Exps (Table S5), suggesting that integration of the data did not introduce a systematic bias. Analysis in the integrated panel showed that the LD decayed on average to at least half of its maximum value at a distance of 10 Mb ( Extended Data Fig. 1b ). High broad-sense heritability estimates were observed for grain yield and heading date performance of hybrids (0.77, 0.94) and parents (0.88, 0.98) ( Extended Data Fig. 1c, d , Table S3). The MPH of grain yield varied widely from -0.86 to 2.02 Mg/ha and had a mean value of 0.81 Mg/ha, whereas for heading date the MPH ranged between -3.90 and 3.05 days with a mean value of -0.93 days ( Extended Data Fig. 1e, f , Table S3). The heritability estimates of MPH were 0.68 and 0.78 for grain yield and heading date, respectively ( Extended Data Fig. 1c, d , Table S6). Thus, the data provide a reliable foundation to dissect the genetic architecture of heterosis for the traits of interest. Download figure Open in new tab Extended Data Figure 1. Genetic characteristics of whole-genome re-sequenced wheat hybrid populations. (a) The first two principal coordinates (PCs) revealed by a principal coordinate analysis (PCoA) for 578 parental lines and 10 checks based on the pairwise Rogers’ distance matrix. Percentages in parentheses referred to the proportion of genotypic variance explained by the corresponding PCs. (b) Decay of linkage disequilibrium (LD) within each subgenome and across the whole genome estimated in the integrated panel. Heritability estimates of (c) grain yield (GY, Mg/ha) and (d) heading date (HD, days) for the population of hybrids and parents as well as for midparent heterosis (MPH) in each experimental series and the integrated panel, respectively. The phenotypic distribution of MPH within each experimental series and in the integrated panel for GY (e) and HD (f) . Resequencing data uncover the advantages of hQTL-ODS Before exploring hQTL in the integrated data set, we compared the performance of hQTL-ODS with hQTL-MSS in Exp I using a low-density marker panel. In a previous study, hQTL-MSS had been applied to Exp I genotyped by a 90k Illumina Infinium SNP chip to study the MPH of grain yield 19 . Thus, we directly adopted the results therein ( Fig. 3a ) and applied hQTL-ODS to the same data set ( Fig. 3b ). hQTL-MSS and hQTL-ODS detected 37 and 19 hQTL, respectively, of which only nine and eight colocalized with each other ( Fig. 3d ). A certain degree of inconsistency was not unexpected, as it has been shown both in theory and in simulation that hQTL-MSS might result in a biased estimation of the heterotic effect due to ignoring small genetic component effects. Download figure Open in new tab Figure 3. A comparison of hQTL for grain yield heterosis detected by hQTL-ODS and hQTL-MSS using experimental series I data. Manhattan plots showing the results of (a) hQTL-MSS using 90k SNP chip data (17,372 SNPs), which were adopted from Jiang et al. (2017), (b) hQTL-ODS using the same 90k SNP chip, and (c) hQTL-ODS using whole-genome sequencing data (about 1.2 million SNPs). Significance thresholds (P < 0.05 after Bonferroni-Holm correction for multiple testing) were indicated as red dashed horizontal lines. (d) Venn diagram showing the number of hQTL detected in (a), (b), and (c), (indicated in blue, purple and black, respectively) as well as the overlapping. Note that it is possible that different numbers of hQTL detected in two or three data sets were colocalized with each other, since significant SNPs were merged into hQTL in each of the three data sets separately. To be consistent with the results from Jiang et al. (2017), we only considered hQTL located on known chromosomes. As mentioned in the previous subsection, the population was also genotyped by WGS. A comparison of the SNPs obtained via WGS with those of the 90k SNP chip revealed an average SNP concordance rate of 95.6%, indicating that results of subsequent analyses based on the two marker panels can be reliably compared. Then, hQTL-ODS was applied to this population with WGS data and 165 hQTL were identified, 124 of which mapped to known chromosomes ( Fig. 3c ). Comparing with the results obtained using the SNP chip, we found that 105 hQTL were only detected with the WGS data ( Fig. 3d ), which demonstrated the superiority of high-density marker panels. Notably, 15 of the 105 hQTL (13.9%) showed neither significant dominance nor epistatic interactions effects with other loci. These hQTL can only be detected using hQTL-ODS as it does not ignore small component effects contributing to the heterotic effect. Interestingly, we observed that nine hQTL identified by hQTL ODS using the SNP chip were not detected by the same model with WGS data. Furthermore, six hQTL detected by hQTL-MSS using the SNP chip colocalized with 11 hQTL detected by hQTL-ODS when adopting WGS, but missed by hQTL-ODS using the SNP chip ( Fig. 3d ). Thus, using a low-density marker panel, the estimation of heterotic effect by hQTL-ODS may be biased because not all genetic variants are considered. The most reliable results are consequently expected by applying hQTL-ODS with a high-density marker panel, such as WGS. Application of hQTL-ODS highlighted pervasive cumulative epistatic effects We applied hQTL-ODS to the integrated panel as well as to the three individual populations (Exp I, II and III) using WGS data to dissect the genetic architecture of heterosis for grain yield and heading date. Thus, results were obtained in four different populations ( Fig. 4a-d , Extended Data Fig. 2a-d ). For convenience, hQTL detected in at least two populations were called common hQTL, and those identified in only one population were termed unique hQTL. For grain yield, 3,174 significant SNPs (including those without assigned chromosome information) merged into 188 hQTL were detected in the integrated population (Table S7), with a high proportion (74.5%) being common hQTL ( Fig. 4f ). In contrast, proportions of common hQTL identified in the three individual populations only amounted to 37.6% (Exp I), 39.7% (Exp II) and 25% (Exp III). Similar results were obtained for the trait heading date ( Extended Data Fig. 2a-d, f ). Thus, the reliability of hQTL detection was clearly enhanced in the integrated population because of its large sample size. Download figure Open in new tab Extended Data Figure 2. Genetic architecture of midparent heterosis (MPH) for heading date in wheat. Results of hQTL-ODS revealed genomic regions harboring loci associated with heading date MPH based on whole-genome sequencing data in the experimental series I ( a ), II ( b ), III ( c ), and the integrated population ( d ). (e) Results of a genome-wide scan for dominance effects of heading date MPH in the integrated population. In (a-e), significance thresholds ( P < 0.05 after Bonferroni-Holm correction for multiple testing) were indicated as red dashed horizontal lines. (f) Venn diagram showing the number of hQTL detected in (a), (b), (c) and (d) (indicated in blue, purple orange and black, respectively). Note that since the procedure of merging significant SNPs into hQTL was performed in each of the four data sets separately, an hQTL detected in one data set may colocalize with multiple hQTL detected in another data set. Therefore, the number of overlapping hQTL was separately indicated for each data set. The significant digenic epistatic interactions between hQTL and all other SNPs were shown by the colored links in the centers of the circles for additive-by-additive ( g ), additive-by-dominance ( h ), dominance-by-additive ( i ), and dominance-by-dominance interactions ( j ). The significant digenic epistatic effects were detected by P < 0.01 after Bonferroni-Holm correction for multiple testing. Download figure Open in new tab Figure 4. Genetic architecture of midparent heterosis (MPH) for grain yield in wheat. Results of hQTL-ODS revealed genomic regions harboring loci associated with grain yield MPH based on whole-genome sequencing data in the experimental series I ( a ), II ( b ), III ( c ), and the integrated population ( d ). (e) Results of a genome-wide scan for dominance effects of grain yield MPH in the integrated population. In (a-e), significance thresholds ( P < 0.05 after Bonferroni-Holm correction for multiple testing) were indicated as red dashed horizontal lines. (f) Venn diagram showing the number of hQTL detected in (a), (b), (c) and (d) (indicated in blue, purple, orange and black, respectively). Note that since the procedure of merging significant SNPs into hQTL was performed in each of the four data sets separately, an hQTL detected in one data set may colocalize with multiple hQTL detected in another data set. Therefore, the number of overlapping hQTL was separately indicated for each data set. The significant digenic epistatic interactions between hQTL and all other SNPs were shown by the colored links in the centers of the circles for additive-by-additive ( g ), additive-by-dominance ( h ), dominance-by-additive ( i ), and dominance-by-dominance interactions ( j ). The significant digenic epistatic effects were detected by P < 0.01 after Bonferroni-Holm correction for multiple testing. From now on, we focused exclusively on the results obtained with the integrated panel. The PVE of the 188 hQTL for grain yield ranged from 1.67% to 13.59%. The hQTL with the highest PVE was detected on chromosome 6B (Integrated_GY_hQTL139), spanning a large region from 47.65 Mb to 124.41 Mb and encompassing 595 significant SNPs. The two most significant hQTL (i.e., containing SNPs with lowest p values) were detected on chromosome 7B (Integrated_GY_hQTL163) and 4A (Integrated_GY_hQTL128), with 11.18% and 12.62% PVE, respectively. GWAS for dominance effects contributing to heterosis revealed only two dominance QTL (dQTL) on chromosomes 5B and 7B ( Fig. 4e , Table S8), both colocalized with hQTL. On the other hand, many significant epistatic interactions between hQTL and other loci in the genetic background were detected ( Fig. 4g-j ). For heading date, 71 hQTL were found with PVE ranging from 2.42% to 17.22% (Table S9). Notably, Ppd-D1 (TraesCS2D03G0156800) 26 , which colocalized with a dQTL, represents a candidate gene for an hQTL with 8.86% PVE. In total, five dQTL contributing to heterosis were found and three colocalized with hQTL ( Extended Data Fig. 2e , Table S10). Thus, for both traits, a large number of hQTL contrasted with few QTL for dominance effects, demonstrating the overriding importance of epistatic interactions for hQTL in wheat. Next, we exemplarily investigated the peak SNPs of three hQTL for grain yield which did not colocalize with any dQTL in more detail to disentangle the composition of component effects contributing to these hQTL effects. Indeed, the dominance effects of their peak SNPs were not significant ( Fig. 5a-c ). Thus, the hQTL effects must be comprised of epistatic effects, as shown in the distribution of phenotypic values for two-locus genotype classes consisting of the peak SNP and the SNP which has the strongest interaction with it ( Fig. 5d-f ). Interestingly, the peak SNPs of the three hQTL had distinct patterns of cumulative epistatic effects with other SNPs. SNP_4A_748345003 in Integrated_GY_hQTL128 had many significant (P < 0.05 after Bonferroni correction) epistatic interactions with other SNPs ( Fig. 5g ), whereas only several significant epistatic effects were identified for SNP_6B_653527460 in Integrated_GY_hQTL142 ( Fig. 5h ). Download figure Open in new tab Figure 5. Distinct patterns of cumulative epistatic effects highlighted by three heterotic QTL (hQTL). The three hQTL were detected for grain yield MPH in the integrated data set: Integrated_GY_hQTL128 (a, d, g) , Integrated_GY_hQTL149 (b, e, h) , and Integrated_GY_hQTL142 (c, f, i). (a-c) The phenotypic distribution of grain yield MPH for hybrids harboring different alleles for the peak SNP within hQTL. (d-f) The phenotypic distribution of grain yield MPH for hybrids with distinct genotype combinations of the peak SNP within hQTL and the SNP showing the most significant epistatic interaction effect with the peak SNP. In each group, the mean MPH value was indicated by red dots and the number of hybrids was shown below the boxplot. (g-i) . Circled Manhattan plot of epistatic effects between the peak SNP of hQTL and all other SNPs. in which the 10 SNPs with most significant epistatic effect was highlighted by red circles and linking lines in the center. Different colors were used to indicate the types of epistatic effects: orange (additive-by-additive), blue (additive-by-dominance), green (dominance-by-additive) and purple (dominance-by-dominance). Significance thresholds ( P < 0.05 after Bonferroni-Holm correction for multiple testing) were indicated as red dashed lines. Importantly, SNP_6B_180998007 in Integrated_GY_hQTL149 had no significant interactions with any other SNP ( Fig. 5i ). These latter results indicated that significant epistatic effects are not necessary for a marker showing a significant heterotic effect. Instead, small epistatic effects can cumulatively compose an hQTL. SNP_6B_180998007 is a representative for 66 out of 188 hQTL (35.1%) that can only be detected using hQTL-ODS as it does not ignore small component effects contributing to the heterotic effect. Discussion In this study, we developed hQTL-ODS, a reliable and time-efficient association analysis tool that comprehensively incorporates both dominance and epistatic effects to unravel the genetic architecture of heterosis, making it feasible to conduct GWAS for heterotic effects using large-scale diverse hybrid populations with WGS data. Previous studies on heterosis in various species (e.g. Arabidopsis 27 , 28 , rice 1 , 7 , 29 and maize 8 ) were based on WGS data with diverse hybrid populations. However, epistatic effects were either not considered or restricted to the interaction between dominance QTL in these studies, possibly due to the high computational burden. It was shown that considering epistasis is important for studying heterosis in wheat 19 , 30 , maize 31 and pigeonpea 32 . These studies applied the hQTL-MSS approach to test epistatic effects between all marker pairs, either based on a low-density marker panel or based on WGS in a relatively small population. As demonstrated in Results, it would be infeasible to apply hQTL-MSS to large populations with WGS. Our new model hQTL-ODS successfully removed the computational bottleneck by modeling the heterotic effect as a whole, hence avoiding the two-dimensional scan for digenic epistatic effects between all pairs of markers. Importantly, the computational efficiency is not the only advantage of hQTL-ODS over existing approaches. Treating a heterotic effect in its entirety averts the possible bias of focusing only on large component effects, as the heterotic effect of a locus can be significant due to the cumulation of small interaction effects between this locus and many other loci ( Fig. 5i ). Thus, even if the computational load of hQTL-MSS could possibly be reduced by integrating more efficient GWAS algorithms for the two-dimensional epistasis scan (e.g. REMMAX 33 , NGG 34 ), hQTL-ODS would still be superior because of higher power and fewer false positives ( Fig. 2f, g ). From the general viewpoint of statistical genetics, the framework of hQTL-ODS is a kernel-based association test for a set of genetic variants 35 , 36 . Thus, although the model presented here is specific to the study of MPH, its framework can be easily expanded to various other applications. For instance, it can be used to study the better-parent heterosis (BPH) instead of MPH, provided that the matrix which transforms the vector of original trait values to MPH is appropriately modified. It can also be applied to study original traits, for which one only needs to omit the step of applying the transformation matrix (see Supplementary Notes for details). In this regard, our model is similar to the so-called marginal epistasis test implemented in MAPIT 37 and FAME 38 . The difference is that in these two models only the cumulative additive-by-additive effects were considered, whereas all three types of digenic epistatic effects were integrated in hQTL-ODS. Since the model transformation from the original trait to MPH is only implemented in hQTL-ODS but not in the marginal epistasis tests, only the former model can be applied to study heterosis. Although hQTL-ODS delivers a substantial improvement compared with existing approaches in detecting heterotic QTL, there are still limitations in its computational efficiency if one considers the ever-growing data size. The theoretical time complexity of hQTL-ODS is about O ( tn 3 p ), indicating that the computing time will increase by a factor of 27 if the sample size triples. As a result, it took about 30 days to complete the analysis for the integrated panel with 4,885 hybrids and 907,534 SNPs, while it only took hours for the three individual populations (each with about 1,600 hybrids). Thus, in future, it is necessary to explore possibilities to further increase the computational efficiency. One option might be to switch from the LR test to the score test, possibly at the cost of decreasing power 35 . An even more important direction is to increase the efficiency of solving LMMs, which is the most time-consuming step in hQTL-ODS. Currently, we utilized a fast REML approach implemented in the R package gaston 39 . Note that the model involves multiple random vectors of genetic effects with dense covariance matrices. For this type of LMM, only a few fast algorithms are available: the Min-Max algorithm implemented in MM4LMM 40 , the Method-of-Moments approach in GEMMA 41 and RHE-mc 42 , the Monte Carlo REML method in BOLT-REML 43 , and a recently developed fast REML algorithm called MPH 44 . It would be worthwhile to investigate whether the computational efficiency of hQTL-ODS can be further improved by exploiting and/or combining some of the techniques implemented in these algorithms, or whether it may even be necessary to develop new efficient algorithms for the future. Methods MPH and the heterotic effect In this study, we consider a population consisting of n hybrids derived from crossing r diverse parental lines. For a trait in consideration, the MPH of a hybrid is defined as the difference between the trait values of the hybrid and the average of its parents. Let y ori be the ( n + r )-dimensional vector of the observed original trait values of all hybrids and parents, and y MPH be the n -dimensional vector of the MPH values of all hybrids. Then, it is clear from the definition that y MPH = Ty ori , where T is an n × ( n + r ) matrix of linear transformation. The heterotic effect of a marker is the net contribution of this particular marker to MPH, taking its dominance effect and its digenic epistatic interaction with the entire genetic background into consideration. In particular, it depends not only on the genotype of the hybrid, but also on the genotype combinations of the parental lines. In the following, we briefly recall its precise definition 19 . Assume that all markers are biallelic and let p be the number of all markers. The genotype of an individual at a marker is coded as 0, 1, or 2, depending on the number of reference alleles. Let F be a hybrid individual, P 1 and P 2 be its parental lines. Let d i be the dominance effect of the i -th marker, aa ij , ad ij , da ij (or equivalently, ad ij ) and dd ij be the additive-by-additive, additive-by-dominance, dominance-by-additive and dominance-by-dominance epistatic effect between the i -th and the j -th marker, respectively. Let R kl ( k, l = 0 or 2) be the subset of markers for which P 1 has genotype code k and P 2 has code l . Then, for any i (1 ≤ i ≤ p ), the heterotic effect of the i -th marker for this particular hybrid individual F , denoted by h i,F , is defined as follows: With the above definition, the MPH value of the individual F is the sum of the heterotic effects of all markers, i.e., . Let h i be the n-dimensional vector consisting of the heterotic effects h i,F for all hybrids. Let M A be the ( n + r ) × p matrix of markers coded as 1 (homozygous for the reference allele), 0 (heterozygous) or -1 (homozygous for the alternative allele), and M D be the ( n + r ) × p matrix of markers coded as 0 (homozygous) or 1 (heterozygous). Then, h i can be written as follows: where m i and m j are the i -th and the j -th columns of the matrix M A , l i and l j are the i -th and the j -th columns of M D , and “∘” denotes entry-wise product of two vectors. Note that it is not straightforward to obtain Eq. (4) from Eq. (3) , as the expression of h i,F in terms of the genetic component effects ( d i , aa ij , ad ij , ad ij and dd ij ) depends on the combination of parental genotypes (i.e., whether i ∈ R 20 , R 02 , R 22 or R 00 ). A proof is provided in the Supplementary Note . The hQTL-ODS model Since MPH is a derived trait, it is necessary to model the genetic effects in a way that they are consistent to those contributing to the original trait. Thus, the baseline model is the following model for the original trait, taking the additive and dominance effects of all markers, and digenic epistatic effects of all marker pairs into account: where y ori , M A and M D were defined in the previous subsection (recall that the columns of M A and M D are denoted by m i and l i (1 ≤ i ≤ p ), respectively), 1 n + r is the ( n + r )-dimensional vector of ones, μ is the common intercept term, α is a k -dimensional vector of covariate effects (e.g., subpopulation effects), X c is the corresponding ( n + r ) × k design matrix, a is the p -dimensional vector of additive effects for all markers, d is the p -dimensional vector of dominance effects, aa is the p ( p − 1)/2-dimensional vector of additive-by-additive epistatic effects between all marker pairs, M AA is an ( n + r ) × p ( p − 1)/2 matrix whose columns consist of m i ∘ m j for all i, j such that 1 ≤ i < j ≤ p. ad is the p ( p − 1)-dimensional vector of additive-by-dominance epistatic effects, M AD is an ( n + r ) × p ( p − 1) matrix whose columns consist of m i ∘ l j for all i, j such that 1 ≤ i, j ≤ p and i ≠ j. dd is the p ( p − 1)/2-dimensional vector of dominance-by-dominance epistatic effects, M DD is an ( n + r ) × p ( p − 1)/2 matrix whose columns consist of l i · l j for all i, j such that 1 ≤ i < j ≤ p, e is an ( n + r )-dimensional vector of residuals. Since y MPH = Ty ori , the model for MPH is naturally obtained by left-multiplying the matrix T to both sides of Eq. (5) : The definition of T implies that T 1 n + r = 0 n . For each hybrid, the marker coding of M A is exactly the average of the coding of its parents, implying TM A = 0 n × p . Let X = TX c , g D = TM D d, g AA = TM AA aa, g AD = TM AD ad, g DD = TM DD dd and ε = Te . Then, Eq. (6) can be written as the following: which is the same as Eq. (1) in the main text. In the model, α is assumed to be fixed effects, and , where , and * denotes D, AA, AD or DD. In particular, the covariance matrix for ε is TT ′ instead of an identity matrix, because the residuals e of the original trait are usually assumed to be independent (i.e., ) and ε = T e . Eq. (7) is the null model of hQTL-ODS, the role of random genetic effects g D , g AA , g AD and g DD is for controlling the genetic relatedness in the population. If there is strong population structure, additional covariates can be included in Xα . To test the heterotic effect of a marker, the heterotic effect h i is added into the model, resulting in the following alternative model: where h i has been defined in Eq. (4) . Let γ i be the (4 p − 3)-dimensional vector consisting of d i , aa ij , ad ij , ad ij and dd ij , and Z i be the n × (4 p − 3) matrix whose columns consist of and , for all 1 ≤ j ≤ p and j ≠ i . Then, Eq. (4) can be rewritten as h i = Z i γ i . Assume that , where I n is the n × n identity matrix, . Then, we have , where . With the above consideration, a classical likelihood ratio test is used to assess the significance of h i . Namely, the null model ( Eq. 7 ) and the alternative model ( Eq. 8 ) are solved by the REML approach, and the null hypothesis is tested by the following test statistic: where L 0 and L 1 are the maximum values of the likelihood function of the null and the alternative model, respectively. Under the null hypothesis, λ LR follows approximately an equal-weighted mixture of two chi-squared distributions with zero and one degrees of freedom 20 , 21 , i.e., . The models (7) and (8) are implemented by using the R package gaston 39 . Other details on the implementation are provided in the Supplementary Note . Here, we just mention three important points: 1) Unlike the model of GWAS for the original traits, the residual term of hQTL-ODS model has a non-trivial covariance matrix ( TT ′ ). Thus, a further transformation of the model has to be made before applying standard algorithms to solve the LMM. 2) Calculating the matrices K AA , K AD , K DD and H i is time-consuming if they are directly calculated from definition. Taking H i = Z i Z i ′ as an example, it takes O ( n 2 p ) time because the dimension of Z i is n × (4 p − 3). Thus, the total time complexity is O ( n 2 p 2 ) for calculating H i for all p markers. However, using the techniques similar to an efficient algorithm for calculating epistatic genomic relationship matrices 23 , the time complexity was eventually reduced to O ( n 2 p ). 3) The computational efficiency is further improved by fixing the ratio of the variance components and i.e., estimating these ratios only once in the null model and assuming that they are invariant in the alternative model for each marker. This approach is similar to the so-called P3D approximation in GWAS 22 . The model for testing genetic component effects The hQTL-ODS model provides a direct test for the heterotic effect of a marker without testing any component effect (i.e., the dominance effect of the marker and the digenic epistatic interaction effects between the marker and the entire genetic background). However, if a marker has significant heterotic effect h i , it would be interesting to investigate whether these component effects, namely d i , aa ij , ad ij , ad ij and dd ij for all j ≠ i , 1 ≤ j ≤ p , are significant or not. The following models are used for testing these effects. where all notations have been defined in the previous two subsections. In these models, the effects d i , d j , aa ij , ad ij , ad ij and dd ij are assumed to be fixed instead of random. Note that the additive effects ( a i and a j ) were not included in models (11), (12) and (13), which seems to be a violation to the principle that the corresponding main effects should be included in the model when an interaction effect is tested. But actually, TM A = 0 n × p (see the previous subsection) implies Tm i a i = Tm j a j = 0 n . Therefore, it makes no difference whether to include these effects or not. Let x be any of the effects d i , aa ij , ad ij , ad ij and dd ij . A Wald test is used to assess its significance, and the test statistic has the following form: where is the estimated value of x and is the estimated variance of . It is known that W follows asymptotically a chi-squared distribution with one degree of freedom . Hybrid wheat populations and field trials In this study, we integrated experimental data of three large hybrid wheat populations (Exp I, II, and III) 24 into a large panel. In total, there were 597 elite parental lines and 5,243 F 1 hybrid progenies. The parental lines were chosen to reflect a wide range of diversity in Central Europe. A factorial or partial factorial design was used to generate single-cross hybrids in each population. Exp I was composed by 135 parental lines (120 females and 15 males) and their 1604 hybrid progenies evaluated for grain yield (Mg ha -1 ) in 11 environments (five locations in 2012 and six in 2013) and heading date (days from January 1 st ) in 10 environments (six locations in 2012 and four in 2013) in Germany. Exp II was comprised of 226 parental lines (185 females and 41 males) and their 1815 hybrid progenies evaluated for grain yield in 12 environments (six locations each in 2016 and 2017) and heading date in 20 environments (11 locations in 2016 and nine in 2017) in Germany. Exp III consisted of 236 parental (196 female and 40 male) lines and their 1824 hybrid progenies evaluated for grain yield in 12 environments (six sites each in 2018 and 2019) and heading date in 18 environments (10 sites in 2018 and eight in 2019) in Germany. In each environment, the experimental design consisted of three trials, in which a partially replicated (Exp I) or un-replicated (Exp II and III) alpha lattice design was used. Different genotypes were evaluated in different trials linked by 10 (Exp I) or 11 (Exp II and III) common checks within the environment. For all genotypes, harvesting was performed mechanically, and the grain yield was adjusted to a moisture content of 140 g H 2 O kg −1 . A summary of the three populations was provided in Table S1. Details of the field trials and part of the phenotypic data for grain yield (Exp I, Exp II and the 2018 data of Exp III) have been described in previous studies 24 . DNA Isolation and genome sequencing For DNA extraction and sequencing, the seeds of all parental lines in the three populations were collected and grown in the greenhouse and a single leaf per genotype was harvested from a 10-day-old seedling. Genomic DNA was isolated from young leaf tissue of each genotype using a silica-membrane technology (NucleoSpin 96 Plant II) as described by the manufacturer (Macherey-Nagel), which were used as input for sequencing library generation. Subsequently, the WGS libraries were prepared using the Nextera DNA Flex Library Prep according to the manufacturer’s (Illumina) instructions. Libraries were pooled in an equimolar manner. The multiplexed pool was quantified by quantitative PCR and sequenced (paired-end, 2 × 151 cycles and 10 bp for the index reads) on NovaSeq 6000 at threefold coverage. WGS read processing and SNP calling To avoid bias resulting from base-calling duplicates and adapter contamination, raw read sequences were processed to remove adaptors and low-quality reads with a minimum read quality (q) cutoff of 20. High-quality data with an average of about 17Gb per sample were retained for subsequent analysis. The paired-end reads (172.8 billion) of all genotypes were aligned against the wheat reference genome (Chinese Spring, RefSeq v2.1) using the MEM algorithm of BWA with default parameters. The output was converted to binary alignment map (BAM) format file using SAMtools (v.1.9) and further sorted with NovoSort (v.3.06.05). After merging BAM files for samples sequenced on the same flowcell with Picard (v.2.21.9), variant calling was done using the mpileup and call functions from SAMtools with parameter ‘-DV’. Subsequently, a raw population genotype VCF file was generated, including 181,421,099 SNPs with an average density of 12.45 SNPs per kb. Bi-allelic SNPs with a minimum QUAL (i.e., mapping quality) score of 40, read depth for homozygous ≥ 1 and heterozygous calls ≥ 2 were recalled based on read depth ratios calculated from the DP (total read depth) and DV (depth of the alternative allele) fields in VCF file using a custom AWK script. The SNPs with missing rate ≤ 30%, heterozygosity ≤ 1% and counts of both homozygous genotypes ≥ 10 were kept. Then, the missing values were imputed by Beagle (v 5.2). Consequently, a total of 7,835,467 SNPs with MAF ≥ 5% were retained for the integrated panel (including all parental lines from Exp I, II and III). Pairs of SNPs with squared correlation greater than 0.90 were greedily pruned within 50kb sliding window using plink (v 1.90b6.9) with steps = 5 SNPs until no such pairs remained. The same data cleaning and filtering process was performed separately within each population ( Fig. S1 ). The genotypes of the hybrids were inferred from the parental genotypes. Download figure Open in new tab Figure S1. Schematic illustration of the curation and analysis of whole-genome re-sequencing data. QUAL: the quality value of site; mac: minor allele count; min-alleles: minimum allele number; max-alleles: maximum allele number; MissRate: missing rate of SNPs; Hom1&2: the number of homozygous genotypes for each allele; HetRate: heterozygosity of SNPs; MAF: minor allele frequency; “plink --indep-pairwise 50kb 5 0.90” means pairs of SNPs with squared correlation greater than 0.90 were greedily pruned within 50kb sliding window with steps = 5 SNPs until no such pairs remained. Genotype recalling based on read depth ratios, calculated from DP (total read depth) and DV (depth of the alternative allele), was performed using an AWK script, which is available in a Bitbucket repository: https://bitbucket.org/ipk_dg_public/vcf_filtering . Population structure, LD decay and duplicate genotypes Principal coordinate analysis (PCoA) was performed based on pairwise Rogers’ distances among genotypes 45 to investigate the population structure. The decay of pairwise LD between SNPs (measured as r 2 ) was analyzed by PopLDdecay 46 with parameters “-MaxDist 50000 -OutType 2”. The average physical distance at which the LD between SNPs decayed to half maximum value was about 10 Mb. Duplicate genotypes might occur when different labels were used for the same genotype in different populations. Thus, we investigated whether duplicate parental lines existed across populations by applying a principle of genetic distance 24 . Groups of genotypes with pairwise Rogers’s distances below 0.03 were defined as duplicates and were considered as the same genotype in subsequent analysis. As a result, 60 parental lines were merged into 28 groups, and 330 hybrids were accordingly merged into 162 groups. Curation of phenotypic data Firstly, all data were screened for data-entry errors and outliers by using the Bonferroni-Holm test to judge the residuals standardized by rescaled median absolute deviation 47 . Then, a two-step approach was used to analyze the phenotypic data for the integrated panel. In the first step, the phenotypic data within each environment was analyzed by fitting a linear mixed model including the effect of genotypes, trials, replications nested within trials, and blocks nested within trials and replications. Best linear unbiased estimations (BLUEs) of the genotypes in each environment were obtained and served as the input of the second step, in which the following linear mixed model was fitted: where y ijk is the BLUE of the i -th genotype in the j -th population (Exp I, II or III) evaluated in the k -th environment, t i is a dummy variable which equals 1 if the genotype is a parental line, and 0 if it is a hybrid, μ L is the mean of all lines, μ H is the mean of all hybrids, the effect of the i -th genotype is denoted by G L,i (if it is a line) or G H,i (if it is a hybrid), s j is the effect of the j -th population, E k is the effect of the k -th environment, and e ijk is the residual term. All effects except μ L , μ H and sj were assumed to be random. For any i and j , we assume that , and , where and are the variance components of the lines, hybrids and residuals, respectively, and is the residual variance in the k -th environment (i.e., we assume variance heterogeneity for the residuals in different environments). Covariance between each pair of these variables was assumed to be zero. The significance of each variance component was tested by the likelihood ratio test. The broad-sense heritability for the lines and hybrids was calculated separately using the following formulas: where is the average residual variance across all environments, N E,L and N E,H are the average numbers of environments in which the lines and hybrids were evaluated, respectively. Then, model (16) was fitted once again with slightly different assumptions, namely G L,i and G H,i were assumed to be fixed effects instead of random, to obtain the cross-environment BLUEs of the genotypic values for all hybrids and parental lines. Similar analysis was performed separately within each population. In this case, the model (16) was fitted without the effect s j . Estimating the heritability of MPH For the integrated panel, we first calculated the MPH value for each hybrid within each environment based on the within-environment BLUEs of all hybrids and parental lines. Then, the following linear mixed model was fitted: where y MPH,ijk is the MPH value of the i -th hybrid of the j -th population in the k -th environment, μ MPH is the mean across populations, G MPH,i is the genotypic effect of the i -th hybrid, s MPH,j is the effect of the j -th population, E MPH,k is the effect of the k -th environment, and ε ijk is the residual. In the model, μ MPH and s MPH,j are assumed to be fixed and the other effects are random: , and , where and are the variance components of genotypes and environments, is the residual variance in the k -th environment. Then, the broad-sense heritability of MPH was estimated using the following formula: where is the average residual variance across all environments and N E,H is the average number of environments in which the hybrids were evaluated, as in Eq. (17) . Similar analysis was performed separately within each population. In this case, the model (18) was fitted without the effect s MPH,j to estimate the within-population heritability of MPH. All linear mixed models in this and the previous subsection were implemented using the package ASReml-R 4.1 48 . Applying hQTL-ODS to the hybrid wheat data set For each population (Exp I, II and III) as well as the integrated panel, the genomic and phenotypic data were combined for detecting hQTL by using the hQTL-ODS model. The final number of parental lines, hybrids and markers are shown in Fig. S1 . In particular, there were 4,885 hybrids and 545 parental lines with 907,534 high-quality SNPs in the integrated panel. Since no clear population stratification was observed ( Extended Data Fig. 1 ), we didn’t include any covariate in the model (i.e., setting Xα = 0 in Eq. 7 and 8 ). P values obtained in the LR test for heterotic effects were corrected for multiple testing with the Bonferroni-Holm method 49 . The genome-wide threshold was determined to be P < 0.05 after correction. After the genome-wide scan of heterotic effects, the following approach was used to merge significant SNPs into hQTL. First, an interval was determined for each significant SNP by investigating the LD between the significant SNP and the flanking SNPs. On each side of the significant SNP, the nearest genomic position where the LD dropped below 0.3 was defined as the boundary of the interval. Next, the defined intervals were merged into independent hQTL based on the following criteria: 1) Two intervals were merged if they overlapped with each other; 2) Two non-overlapping intervals were merged if the distance between the two peak SNPs in the intervals was less than 10 Mb (determined by the average decay of LD across genome) and the average LD between all SNPs in the two intervals was higher than 0.3. The procedure was repeated until no more intervals could be merged, and the resulting intervals were treated as independent hQTL. In each hQTL, the most significant SNP was referred as the lead SNP. The PVE of each hQTL was estimated by the following formula: where the variance components were estimated by fitting the lead SNP of the hQTL in model (8). In additional to the test for heterotic effects, the dominance effects of all markers were tested using model(10) in the integrated panel. The same procedure as described above was used to merge significant SNPs into dQTL. Moreover, for each hQTL in the integrated panel, the digenic epistatic interaction effects between the SNPs in the hQTL and all SNPs across genome were tested using models (11)-(14). Simulation study A simulation study was performed to compare the performance of hQTL-ODS and hQTL-MSS in terms of statistical power and FPR. The simulation was based on the genomic data of Exp I. To reduce the computational load, we randomly sampled a subset of individuals and markers. As a result, the data used for simulation consisted of 90 parental lines (75 females and 15 males) and their 1,000 hybrid progenies with 5,000 SNPs. Five different scenarios were considered (Table S2). In each scenario, one marker was randomly sampled as hQTL and its PVE was fixed as 2.5%. The heritability was fixed as 0.5. In Scenario 1, the heterotic effect was contributed by cumulative small epistatic effects. One marker from each of the remaining 20 chromosomes (excluding the one on which the hQTL was located) was randomly sampled. The 20 markers interacted with the hQTL, and the pattern of interaction (additive-by-additive, additive-by-dominance, dominance-by-additive, and dominance-by-dominance) was randomly assigned. The PVE for each of the 20 epistatic interaction effects was fixed as 0.5%. Scenario 2 is similar to Scenario 1, but only five markers interacted with the hQTL and the PVE for each of the five epistatic effects was 2%. In Scenario 3, the significant heterotic effect was due to cumulative small dominance and epistatic effects. Thus, the only difference between Scenario 3 and 1 is that one epistatic effect was replaced by the dominance effect of the hQTL itself. Thus, the PVE for the dominance effect of the hQTL and for each of the 19 epistatic interaction effects was 0.5%. Scenario 4 is similar to Scenario 3, but only four markers interacted with the hQTL. The PVE for the dominance effect of the hQTL and for each of the four epistatic effects was 2%. In Scenario 5, the significant heterotic effect was solely due to the dominance effect of the hQTL with 2.5% PVE. Therefore, no epistatic interaction effect was simulated in this scenario. The simulated phenotypic data (MPH values of 1,000 hybrids) were produced by the following steps: Simulating the genetic background effect. Assuming that , we randomly sampled vectors and , where K D , K AA , K AD and K DD were defined previously (after Eq. 7 ). Then, the vector of genetic background effects was . Since the heritability was fixed as 0.5 and the PVE of hQTL was 2.5%, we knew that the PVE of the genetic background effects should be 47.5%. Thus, the total phenotypic variance V t was calculated as , where S (·) denote the sample variance. Simulating the hQTL effect h i using the following formula (reformulated from Eq. 4 ): where s c = 0 for scenarios 1 and 2, and s c = 1 for scenarios 3, 4 and 5; Q aa , Q ad , Q da , and Q dd are the subsets of markers that were sampled to interact with the hQTL with the assigned pattern. Note that each of them could be the empty set, depending on the scenario and the result of random assignment. For example, all of them were empty sets in scenario 5 because no epistatic effects were simulated. To simplify the description, we denote by x a component effect (i.e., x is d i , aa ij , ad ij , da ij or dd ij ), its coefficient in Eq. (21) by , and its PVE by P x . Thus, P x = 0.5% in scenarios 1 and 3, P x = 2% in scenarios 2 and 4, and P x = 2.5% in scenario 5. Then, the variance explained by x is var( c x x ) = var( c x ) x 2 , which should be V t · P x . Thus, the effect size x was calculated as . After all component effects were computed, they were summed up to produce h i by Eq. (21) . Then, the PVE of h i is S ( h i )/ V t , but this value is not guaranteed to be the predetermined value (i.e., 2.5%), nor is it necessarily equal to the sum of all P x , because there exists covariance between the coefficients c x . Practically, we changed the covariance structure by repeating the random sampling of the markers that interact with the hQTL until S ( h i )/ V t equaled 2.5%. Simulating the residual effect. Since the heritability was fixed as 0.5, the residual variance is . Thus, we randomly sampled . Generating the simulated MPH values according to Eq. (8) , i.e., . Note that for simplicity, we didn’t simulate any covariate effects. We applied hQTL-ODS and hQTL-MSS to each simulated data set and compared the results. Since only one hQTL was simulated in each data set, the statistical power in each scenario was estimated as the proportion of data sets in which the marker simulated as hQTL was successfully detected. The FPR was estimated as the average proportion of null markers (i.e., 4,999 markers that were not simulated as hQTL) that were incorrectly identified as significant across 100 simulations. Note that we used simple criteria based on markers instead of regions, as LD between markers was in most cases not high, given that only 5,000 markers were used for simulation. Data availability Raw sequence data collected in this study have been deposited at the European Nucleotide Archive ( https://www.ebi.ac.uk/ena/ ) under the following Project IDs: PRJEB48738 (parental lines from Exp I, released in a previous study 50 ) and PRJEB82869 (parental lines from Exp II and III). VCF files for whole WGS data are available from European Variation Archive ( https://www.ebi.ac.uk/eva/ ) under project PRJEB87554. The phenotypic data of all parents and hybrids, as well as the simulated data sets are provided in GitHub ( https://github.com/Ligl0226/hQTL-ODS/ ). Code availability The hQTL-ODS model was implemented as easy-to-use R functions. All codes are available in GitHub ( https://github.com/Ligl0226/hQTL-ODS/ ). Acknowledgement This study was partly supported by the German Research Foundation (DFG) within the Project “Developing a statistical genomics toolbox to decipher the genetic architecture of heterosis using whole-genome sequencing data” (grant number: 540803247). Guoliang Li was financially and logistically supported by a Sino-German (CSC-DAAD) postdoc scholarship program with China Scholarship Council (CSC), Grant no. 202006350258. Funder Information Declared Deutsche Forschungsgemeinschaft, https://ror.org/018mejw64 , 540803247 China Scholarship Council , 202006350258 Footnotes https://www.ebi.ac.uk/ena/ https://www.ebi.ac.uk/eva/ References 1. ↵ Huang , X. et al. Genomic architecture of heterosis for yield traits in rice . Nature 537 , 629 – 633 ( 2016 ). OpenUrl CrossRef PubMed 2. ↵ Gornicki , P. & Faris , J. D. Rewiring the wheat reproductive system to harness heterosis for the next wave of yield improvement . Proc. Natl. Acad. Sci. U.S.A . 111 , 9024 – 9025 ( 2014 ). OpenUrl FREE Full Text 3. ↵ Ribeiro , A. M. F. et al. Genetic parameters, heterosis, and breed effects for body condition score and mature cow weight in beef cattle . Journal of Animal Science 100 , skac017 ( 2022 ). OpenUrl CrossRef PubMed 4. Ni , A. et al. Genetic parameters, reciprocal cross differences, and age-related heterosis of egg-laying performance in chickens . Genet Sel Evol 55 , 87 ( 2023 ). OpenUrl CrossRef PubMed 5. ↵ Cassady , J. P. , Young , L. D. & Leymaster , K. A. Heterosis and recombination effects on pig reproductive traits . Journal of Animal Science 80 , 2303 ( 2002 ). OpenUrl CrossRef PubMed 6. ↵ Paril , J. , Reif , J. , Fournier-Level , A. & Pourkheirandish , M. Heterosis in crop improvement . The Plant Journal 117 , 23 – 32 ( 2024 ). OpenUrl CrossRef PubMed 7. ↵ Gu , Z. et al. Structure and function of rice hybrid genomes reveal genetic basis and optimal performance of heterosis . Nat Genet 55 , 1745 – 1756 ( 2023 ). OpenUrl CrossRef PubMed 8. ↵ Xiao , Y. et al. The genetic mechanism of heterosis utilization in maize improvement . Genome Biol 22 , 148 ( 2021 ). OpenUrl CrossRef PubMed 9. ↵ East , E. M. Heterosis . Genetics 21 , 375 – 397 ( 1936 ). OpenUrl FREE Full Text 10. Hull , F. H. Recurrent selection for specific combining ability in corn . J. Am. Soc. Agron 37 , 134 – 145 ( 1945 ). OpenUrl CrossRef 11. ↵ Crow , J. F. Alternative hypotheses of hybrid vigor . Genetics 33 , 477 – 487 ( 1948 ). OpenUrl FREE Full Text 12. ↵ Bruce , A. B. The Mendelian Theory of Heredity and the Augmentation of Vigor . Science 32 , 627 – 628 ( 1910 ). OpenUrl FREE Full Text 13. Jones , D. F. Dominance of Linked Factors as a Means of Accounting for Heterosis . Proc. Natl. Acad. Sci. U.S.A . 3 , 310 – 312 ( 1917 ). OpenUrl FREE Full Text 14. ↵ Birchler , J. A. , Yao , H. & Chudalayandi , S. Unraveling the genetic basis of hybrid vigor . Proc. Natl. Acad. Sci. U.S.A . 103 , 12957 – 12958 ( 2006 ). OpenUrl FREE Full Text 15. ↵ Richey , F. D. Mock-Dominance and Hybrid Vigor . Science 96 , 280 – 281 ( 1942 ). OpenUrl FREE Full Text 16. Powers , L. An Expansion of Jones’s Theory for the Explanation of Heterosis . The American Naturalist 78 , 275 – 280 ( 1944 ). OpenUrl CrossRef 17. ↵ Jinks , J. L. & Jones , R. M. Estimation of the components of heterosis . Genetics 43 , 223 – 234 ( 1958 ). OpenUrl FREE Full Text 18. ↵ Melchinger , A. E. , Utz , H. F. , Piepho , H.-P. Zeng , Z.-B. & Schön , C. C. The Role of Epistasis in the Manifestation of Heterosis: A Systems-Oriented Approach . Genetics 177 , 1815 – 1825 ( 2007 ). OpenUrl Abstract / FREE Full Text 19. ↵ Jiang , Y. , Schmidt , R. H. , Zhao , Y. & Reif , J. C. A quantitative genetic framework highlights the role of epistatic effects for grain-yield heterosis in bread wheat . Nat Genet 49 , 1741 – 1746 ( 2017 ). OpenUrl CrossRef PubMed 20. ↵ Chernoff , H. On the Distribution of the Likelihood Ratio . The Annals of Mathematical Statistics 25 , 573 – 578 ( 1954 ). OpenUrl CrossRef 21. ↵ Visscher , P. M. A Note on the Asymptotic Distribution of Likelihood Ratio Tests to Test Variance Components . Twin Res Hum Genet 9 , 490 – 495 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 22. ↵ Zhang , Z. et al. Mixed linear model approach adapted for genome-wide association studies . Nat Genet 42 , 355 – 360 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 23. ↵ Jiang , Y. & Reif , J. C. Efficient Algorithms for Calculating Epistatic Genomic Relationship Matrices . Genetics 216 , 651 – 669 ( 2020 ). OpenUrl Abstract / FREE Full Text 24. ↵ Zhao , Y. et al. Unlocking big data doubled the accuracy in predicting the grain yield in hybrid wheat . Sci. Adv . 7 , eabf9106 ( 2021 ). OpenUrl FREE Full Text 25. ↵ Zhu , T. et al. Optical maps refine the bread wheat Triticum aestivum cv. Chinese Spring genome assembly . The Plant Journal 107 , 303 – 314 ( 2021 ). OpenUrl CrossRef PubMed 26. ↵ Beales , J. , Turner , A. , Griffiths , S. , Snape , J. W. & Laurie , D. A. A Pseudo-Response Regulator is misexpressed in the photoperiod insensitive Ppd-D1a mutant of wheat (Triticum aestivum L .). Theor Appl Genet 115 , 721 – 733 ( 2007 ). OpenUrl CrossRef PubMed Web of Science 27. ↵ Seymour , D. K. et al. Genetic architecture of nonadditive inheritance in Arabidopsis thaliana hybrids . Proc. Natl. Acad. Sci. U.S.A . 113 , ( 2016 ). 28. ↵ Yang , M. et al. Genomic architecture of biomass heterosis in Arabidopsis . Proc. Natl. Acad. Sci. U.S.A . 114 , 8101 – 8106 ( 2017 ). OpenUrl Abstract / FREE Full Text 29. ↵ Huang , X. et al. Genomic analysis of hybrid rice varieties reveals numerous superior alleles that contribute to heterosis . Nat Commun 6 , 6258 ( 2015 ). OpenUrl CrossRef PubMed 30. ↵ Boeven , P. H. G. et al. Negative dominance and dominance-by-dominance epistatic effects reduce grain-yield heterosis in wide crosses in wheat . Sci. Adv . 6 , eaay4897 ( 2020 ). OpenUrl FREE Full Text 31. ↵ Sang , Z. et al. Epistasis Activation Contributes Substantially to Heterosis in Temperate by Tropical Maize Hybrids . Front. Plant Sci . 13 , 921608 ( 2022 ). OpenUrl CrossRef PubMed 32. ↵ Saxena , R. K. et al. Characterization of heterosis and genomic prediction-based establishment of heterotic patterns for developing better hybrids in pigeonpea . The Plant Genome 14 , e20125 ( 2021 ). OpenUrl CrossRef PubMed 33. ↵ Wang , D. et al. Rapid epistatic mixed-model association studies by controlling multiple polygenic effects . Bioinformatics 36 , 4833 – 4837 ( 2020 ). OpenUrl CrossRef PubMed 34. ↵ Carré , C. et al. Next-Gen GWAS: full 2D epistatic interaction maps retrieve part of missing heritability and improve phenotypic prediction . Genome Biol 25 , 76 ( 2024 ). OpenUrl CrossRef PubMed 35. ↵ Lippert , C. et al. Greater power and computational efficiency for kernel-based association testing of sets of genetic variants . Bioinformatics 30 , 3206 – 3214 ( 2014 ). OpenUrl CrossRef PubMed 36. ↵ Fu , B. et al. Fast kernel-based association testing of non-linear genetic effects for biobank-scale data . Nat Commun 14 , 4936 ( 2023 ). OpenUrl CrossRef PubMed 37. ↵ Crawford , L. , Zeng , P. , Mukherjee , S. & Zhou , X. Detecting epistasis with the marginal epistasis test in genetic mapping studies of quantitative traits . PLoS Genet 13 , e1006869 ( 2017 ). OpenUrl CrossRef PubMed 38. ↵ Fu , B. et al. A biobank-scale test of marginal epistasis reveals genome-wide signals of polygenic epistasis . Preprint at doi: 10.1101/2023.09.10.557084 ( 2023 ). OpenUrl Abstract / FREE Full Text 39. ↵ Dandine-Roulland , C. & Perdry , H. Genome-wide data manipulation, association analysis and heritability estimates in R with Gaston 1.5 . in 46th European Mathematical Genetics Meeting (EMGM , 2018 ) ( Cagliari, Italy ). 40. ↵ Laporte , F. , Charcosset , A. & Mary-Huard , T. Efficient ReML inference in variance component mixed models using a Min-Max algorithm . PLoS Comput Biol 18 , e1009659 ( 2022 ). OpenUrl CrossRef PubMed 41. ↵ Zhou , X. A unified framework for variance component estimation with summary statistics in genome-wide association studies . Ann. Appl. Stat . 11 , ( 2017 ). 42. ↵ Pazokitoroudi , A. et al. Efficient variance components analysis across millions of genomes . Nat Commun 11 , 4020 ( 2020 ). OpenUrl CrossRef PubMed 43. ↵ Loh , P. R. et al. Contrasting genetic architectures of schizophrenia and other complex diseases using fast variance-components analysis . Nat Genet 47 , 1385 – 1392 ( 2015 ). OpenUrl CrossRef PubMed 44. ↵ Jiang , J. MPH: fast REML for large-scale genome partitioning of quantitative genetic variation . Bioinformatics 40 , btae298 ( 2024 ). OpenUrl CrossRef PubMed 45. ↵ Rogers , J. S. Measures of genetic similarity and genetic distance . in Studies in genetics VII 145 – 153 ( 1972 ). 46. ↵ Zhang , C. , Dong , S.-S. , Xu , J.-Y. He , W.-M. & Yang , T.-L. PopLDdecay: a fast and effective tool for linkage disequilibrium decay analysis based on variant call format files . Bioinformatics 35 , 1786 – 1788 ( 2019 ). OpenUrl CrossRef PubMed 47. ↵ Bernal-Vasquez , A.-M. Utz , H.-F. & Piepho , H.-P. Outlier detection methods for generalized lattices: a case study on the transition from ANOVA to REML . Theor Appl Genet 129 , 787 – 804 ( 2016 ). OpenUrl CrossRef PubMed 48. ↵ Butler , D. G. , Cullis , B. R. , Gilmour , A. R. , Gogel , B. G. & Thompson , R. ASReml-R Reference Manual Version 4 . ( VSN International Ltd., Hemel Hempstead, HP1 1ES, UK , 2017 ). 49. ↵ Holm , S. A simple sequentially rejective multiple test procedure . Scand. J. Stat . 65 – 70 ( 1979 ). 50. ↵ Schulthess , A. W. et al. Genomics-informed prebreeding unlocks the diversity in genebanks for wheat improvement . Nat Genet 54 , 1544 – 1552 ( 2022 ). OpenUrl CrossRef PubMed View the discussion thread. Back to top Previous Next Posted May 12, 2025. Download PDF Supplementary Material Data/Code Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following Powerful one-dimensional scan to detect heterotic QTL Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share Powerful one-dimensional scan to detect heterotic QTL Guoliang Li , Renate H. Schmidt , Yusheng Zhao , Jochen C. Reif , Yong Jiang bioRxiv 2025.05.06.652159; doi: https://doi.org/10.1101/2025.05.06.652159 Share This Article: Copy Citation Tools Powerful one-dimensional scan to detect heterotic QTL Guoliang Li , Renate H. Schmidt , Yusheng Zhao , Jochen C. Reif , Yong Jiang bioRxiv 2025.05.06.652159; doi: https://doi.org/10.1101/2025.05.06.652159 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Genetics Subject Areas All Articles Animal Behavior and Cognition (7637) Biochemistry (17705) Bioengineering (13899) Bioinformatics (41970) Biophysics (21463) Cancer Biology (18605) Cell Biology (25526) Clinical Trials (138) Developmental Biology (13385) Ecology (19911) Epidemiology (2067) Evolutionary Biology (24329) Genetics (15615) Genomics (22514) Immunology (17743) Microbiology (40424) Molecular Biology (17194) Neuroscience (88650) Paleontology (667) Pathology (2835) Pharmacology and Toxicology (4827) Physiology (7648) Plant Biology (15160) Scientific Communication and Education (2046) Synthetic Biology (4302) Systems Biology (9825) Zoology (2271)