Genetic Dissection of Hybrid Performance and Heterosis for Yield-Related Traits in Maize

Heterosis contributes a big proportion to hybrid performance in maize, especially for grain yield. It is attractive to explore the underlying genetic architecture of hybrid performance and heterosis. Considering its complexity, different from former mapping method, we developed a series of linear mixed models incorporating multiple polygenic covariance structures to quantify the contribution of each genetic component (additive, dominance, additive-by-additive, additive-by-dominance, and dominance-by-dominance) to hybrid performance and midparent heterosis variation and to identify significant additive and non-additive (dominance and epistatic) quantitative trait loci (QTL). Here, we developed a North Carolina II population by crossing 339 recombinant inbred lines with two elite lines (Chang7-2 and Mo17), resulting in two populations of hybrids signed as Chang7-2 × recombinant inbred lines and Mo17 × recombinant inbred lines, respectively. The results of a path analysis showed that kernel number per row and hundred grain weight contributed the most to the variation of grain yield. The heritability of midparent heterosis for 10 investigated traits ranged from 0.27 to 0.81. For the 10 traits, 21 main (additive and dominance) QTL for hybrid performance and 17 dominance QTL for midparent heterosis were identified in the pooled hybrid populations with two overlapping QTL. Several of the identified QTL showed pleiotropic effects. Significant epistatic QTL were also identified and were shown to play an important role in ear height variation. Genomic selection was used to assess the influence of QTL on prediction accuracy and to explore the strategy of heterosis utilization in maize breeding. Results showed that treating significant single nucleotide polymorphisms as fixed effects in the linear mixed model could improve the prediction accuracy under prediction schemes 2 and 3. In conclusion, the different analyses all substantiated the different genetic architecture of hybrid performance and midparent heterosis in maize. Dominance contributes the highest proportion to heterosis, especially for grain yield, however, epistasis contributes the highest proportion to hybrid performance of grain yield.


INTRODUCTION
Heterosis is the phenomenon that a hybrid outperforms its two parents (Birchler et al., 2006;Lippman and Zamir, 2007). Maize is the most successful example for the utilization of heterosis in crops to improve agricultural production, as single-cross varieties of maize have substantially contributed to the improvement of maize production in the past decades (Hochholdinger and Baldauf, 2018). There are three hypotheses to explain the genetic basis of heterosis: dominance (Bruce, 1910;Jones, 1917), overdominance (East, 1936) and epistasis (Powers, 1944). Many studies were performed to test these hypotheses, but the results often varied, depending on the populations and the traits studied, suggesting that heterosis is a complex genetic phenomenon. One commonly used design to study heterosis is the North Carolina Design III (NCIII) or Triple Testcross Design which allows to estimate the contribution of additive, dominance, and epistasis effects to heterosis (Melchinger et al., 2007b;Garcia et al., 2008). In a maize study, a total of 264 F 3 genotypes were generated by intercrossing B73 and Mo17, and the F 3 genotypes were then backcrossed to the two parents. The results showed that nearly all heterozygous individuals performed better than the homozygous individuals, supporting the overdominance (or pseudo-overdominance) hypothesis (Stuber et al., 1992). Conversely, the analysis of hybrid maize data from another NCIII design showed that dominance loci contributed the most to heterosis in maize, while the additiveby-additive effects contributed the most to the heterosis of rice (Garcia et al., 2008).
An alternative design is the North Carolina Design II (NCII) or factorial design, where a set of males is crossed with a set of females in a balanced or unbalanced way. In a partial NCII of maize, eight main effect (additive and dominance) QTL and 37 epistatic QTL pairs were identified (Bu et al., 2015). In addition to the NC mating designs, advanced maize populations were also developed and used for analysis of heterosis. For example, Wei et al. (2016) detected 36 heterotic loci from a series of singlesegment substitution lines. Using near-isogenic lines for QTL detection, many additive QTL and additive-by-additive QTL pairs were detected (Melchinger et al., 2007a;Reif et al., 2009). An immortalized F 2 population (IMF 2 ) was also a promising mating design for dissecting the genetic basis of heterosis and epistasis QTL (Hua et al., 2003;Xu, 2013;Guo et al., 2014;Yi et al., 2019).
Linear mixed models (LMM) are a powerful tool for the genetic dissection of complex traits and are widely used in plant and animal breeding (Yu et al., 2006;Xu et al., 2014;Cui et al., 2020). In a hybrid population of rice, a LMM incorporating multiple polygenic covariance structures to control the genetic background was developed (Xu, 2013). In wheat, a quantitative genetics approach was proposed to dissect the genetic basis of grain-yield heterosis, allowing QTL mapping of dominance, epistasis and heterotic loci for midparent heterosis (MPH) (Jiang et al., 2017). In addition to QTL mapping, genomic selection (GS) has become a new tool for plant breeding and the genetic dissection of complex traits (Meuwissen et al., 2001) and has been applied to hybrid wheat (Zhao et al., 2013(Zhao et al., , 2015bJiang et al., 2017), hybrid rice (Cui et al., 2020) and hybrid maize (Albrecht et al., 2014;Technow et al., 2014).
The general combining ability (GCA) is a measure for the average performance of a line in different hybrid combinations, while the specific combining ability (SCA) describes the deviation of a hybrid from the performance expected based on the GCA of its two parental lines. The additive and additive-by-additive variances contribute to the variation of GCA, while the nonadditive polygenic variances contribute to the variation of SCA (Reif et al., 2007). A two-step approach has been widely used to study the genetics underlying hybrid performance, where the first step consists of estimating the GCA, SCA and the MPH (Guo et al., 2014;Zhou et al., 2018;Yi et al., 2019) and the second step represents the QTL mapping step with the GCA, SCA and MPH treated as the traits of interest. In a previous genome-wide association study (GWAS) with an NCII population, different coding schemes for the genotypes were applied, namely the additive, dominance and recessive coding (Hyun et al., 2008;Liu et al., 2021). However, the additive model was usually not sufficient to explain hybrid performance and MPH. Thus, more elaborate models incorporating non-additive effects should be used to study heterosis.
In this study, we developed a NCII population of maize by crossing a set of 339 recombinant inbred lines (RILs) with two elite inbred lines, resulting in two populations of hybrids. A total of 10 traits were recorded in four to five environments and high-density genotypic data were obtained by genotyping-bysequencing of the RILs and resequencing of the parents. The aims of this study were to (1) evaluate the heritability of MPH and the relative contribution of various traits to grain yield, (2) perform QTL mapping for main (additive and dominance) and non-additive effect loci for hybrid performance and MPH, (3) identify QTL hotspots for yield-related traits, (4) explore the mechanisms of heterosis and hybrid performance, and (5) assess the accuracy of genomic prediction in various breeding schemes.

Plant Materials
A RIL population consisting of 365 F 11 lines was developed by crossing inbred lines Qi319 as the male parent and Ye478 as the female parent originating from two different heterotic groups of maize (Zhou et al., 2016). Two hybrid populations were developed by crossing the RILs to two female testers, Chang7-2 and Mo17, and the two populations Chang7-2 × RIL and Mo17 × RIL were named TC and TM, respectively (Zhou et al., 2018). Different numbers of offspring were obtained from the two hybrid populations. A total of 339 common lines from the RIL, TC, and TM populations were retained for further analysis.
We used a randomized incomplete block design with two replicates in each environment. To avoid competition, the RIL and the hybrid populations were planted separately. Each genotype was planted in two rows with a row interval of 0.6 m, a row length of 4 m and a plant interval of 0.25 m.

Phenotypic Data Analysis
The combination between year and location was considered as an environment (a total of 4 or 5 environments). The studentized residual razor method (Bernal-Vasquez et al., 2016) was used to remove outliers with a threshold of 2.8. The best linear unbiased estimations (BLUE) of the fixed effects and the variance components of the random effects were estimated using the following model: where y ijk was the phenotypic value of the kth replicate of genotype i from the jth environment, µ was the overall mean, G i (i = 1, 2,..., 339) was the effect of the ith genotype, E j (j = 1, 2,..., 5) was the effect of the jth environment, G * E ij was the genotype-by-environment interaction effect, R k (E j ) (k = 1, 2) was the effect of the kth replicate nested in the jth environment, ε ijk was the residual. For estimation of variance components, all random effects were assumed to be normally distributed with mean 0 and variances denoted by σ 2 G , σ 2 G×E and σ 2 ε for G i , G * E ij and ε ijk , respectively. The broad-sense heritability of a trait was defined as (Falconer and Mackay, 1996), where N E = 4 or 5 was the number of environments and N R = 2 was the number of replicates within each environment. Genetic analysis of MPH was conducted in two steps (Jiang et al., 2017). The first step was represented by BLUE of the trait value for each parent and each hybrid. The BLUE of the trait value obtained from the two replicates in one environment was calculated with the following formula: where y ik was the trait value for the kth replicate of genotype i, µ was the mean of the trait under the current environment, G i was the genetic value of the ith genotype and R k was the effect of the kth replicate assumed to follow a N(0, σ 2 R ) distribution, ε ik was assumed to follow a N(0, σ 2 ik ) distribution.
The MPH was defined as (Melchinger et al., 2007b): Where H was the BLUE value of hybrids, P 1 was the BLUE value of Chang7-2 or Mo17 (corresponding to female parent of hybrid), P 2 was the BLUE value of RIL (corresponding to the male parent of hybrid).
The second step in the MPH analysis required the following mixed model: where MPH ij was the MPH value calculated in the first step for hybrid (genotype) i in environment j, G i (i = 1, 2, . . ., 339) was the genetic effect of MPH for the ith hybrid, E j was effect of the jth environment and ε ij was the residual. Noted that G i was treated as a fixed effect in the BLUE calculation or a random effect following a N(0, σ 2 G ) distribution in variance estimation, E j was treated as a random effect following a N(0, σ 2 E ) distribution and ε ij was assumed to be N(0, σ 2 ε ) distributed. The variance components of the above linear mixed model were implemented using the ASReml 3.0 package in R (Gilmour et al., 2009).
In addition, the hybrid performance was decomposed into GCA, SCA and interaction with the environment using a twostep method. In the first step, the BLUEs in RIL, TC and TM populations were calculated within each environment following the same formula above. In the second step, the following formula was applied to the hybrid performance (Zhao et al., 2015a): Where y was the hybrid performance, µ was the mean, E was the environment effect, GCA RIL was the GCA of RILs, GCA Tester was the GCA of testers, the rest was the interaction between GCA, SCA, and environment, ε was the error. All effects were treated as random following normal distributions. The variances were estimated in ASReml 3.0 package in R (Gilmour et al., 2009).
Path analysis can be used to determine the relative contribution of independent variables to a response variable. Path analysis was implemented in the R package sem by taking GY as the response variable and the other traits as independent variables. Path coefficient p i of variable X i was obtained by where b i was the partial correlation, SS X i and SS Y were sum of square for X i and the response variable Y, respectively. Path diagrams were drawn with the semPlot package in R, where values above 0.14 (p = 0.01, n = 339) were displayed.

Genomic Data Analyses
The genotyping procedures for the RILs, the two parents of the RILs, and the two testers were described in a previous study (Zhou et al., 2016). In brief, for the four parents, the pairedend sequencing libraries were created with a fragment length of ∼500 bp and were sequenced on an Illumina HiSeq 2000 sequencer. The resequencing depth was ∼30×. For the RILs, a genotyping-by-sequencing (GBS) strategy was applied. A total of 137,699,000 reads were generated. On average, there were 357,376 reads per individual, which was approximately a 0.07fold coverage of the maize genome. The cleaned reads were obtained after quality control.
The filtered high-quality reads of the four parents and the RILs were mapped to the reference genome (B73_RefGen_v4) with BWA (Li and Durbin, 2009). SAMtools (Li, 2011) were used to call SNPs with quantity over 20 and a total of 41,791,163 SNPs were finally produced. Details regarding the parameters for the SNP calling process can be found in a previous report (Zhou et al., 2016). After filtering of all SNPs for minor allele frequency < 0.05, missing rate > 0.1 and unknown physical positions, a total of 36,095 SNPs remained in the data set for analysis. Missing genotypes of SNPs were imputed using the BEAGLE software package (version v5) with the default parameters (Browning and Browning, 2016).
The low-coverage high-throughput sequence technologies like GBS generate sequences that are often error-prone, which might lead to errors for detection of genetic variants (Ma et al., 2019). Therefore, the hmm.vitFUN.rils function in the R package MPR.genotyping was used to correct the genotyping errors using a Hidden Markov model with errorRate = 0.05 (Xie et al., 2010). The SNPs with high error probabilities were either corrected or set to missing values.
The bin function in the ICIMapping package was used to bin redundant markers with missing rate > 0.2 and a distortion p-value < 0.001, while missing values and anchor information were considered at the same time (Meng et al., 2015). After the above imputation and correction, there were still a little proportion of missing values left, then the argmax method in qtl/R was used to perform the final imputation additive-bydominance and (Broman et al., 2003). Finally, a total of 4,141 bins were discovered across the entire maize genome. The genetic map was constructed using the map function in the IciMapping package with the default parameter values (Meng et al., 2015).

Mapping Quantitative Trait Loci in Recombinant Inbred Line, TC, TM, and Pooled TC-TM Populations
To determine the contribution of each genetic component to hybrid performance and MPH variation and identify significant non-additive QTL, firstly, we combined the TC and TM populations to form a pooled population called TC-TM. For 341 lines (337RILs, 2 parents; 2 testers) lines, if the genotypes were the same as Ye478, it was coded as "1"; if the genotypes were the same as Qi319, it was coded as "−1", then the genotypes of the hybrids were inferred from their parents (the RILs and the testers). The additive and dominance coding matrices, Z and W, for individual j at marker k were coded as Z jk = 1 0 −1 for the additive effect and W jk = 0 1 0 for the dominance effect. The linear mixed model for variance component analysis was (Xu, 2013;Jiang et al., 2017): where y was an n × 1 vector of phenotypic values of the hybrids and Xβ captured the fixed effects of the model that were not relevant to genetic effects. The design matrix for the fixed effects was X = [X 0 , X 1 ], where X 0 was an n × 1 vector of unity (a vector with all elements being 1) and X 1 was an n × 1 vector indicating one of the two populations, X j1 = 0 for TC and X j1 = 1 for TM. The last term of model (1) was a vector of residuals. The remaining terms in model (1) were various polygenic effects (each polygenic effect was an n × 1 vector) and were defined below. ξ a = m k=1 Z k a k was the polygenic additive effect; The operator # represented element-wise product of matrices. In the formulas above, a k and d k were the additive and dominance effect for marker k, (aa) kk , (ad) kk , and (dd) kk were the additive-by-additive, additiveby-dominance and dominance-by-dominance effect between markers k and k' , respectively. The distributions for the polygenic and residual effects were where K a , K d , K aa , K ad , and K dd were the corresponding kinship matrices calculated using the method given by Xu (2013). The six variance components (five genetic variance components and the residual variance) were estimated using the BGLR package in R (Pérez and De Los Campos, 2014) with the number of iterations set at 15,000 and the number of burn-in set at 5,000.
Let λ x = σ 2 x /σ 2 , where σ 2 x was one of the five genetic variance components and σ 2 was the residual variance. The above variance could be rewritten as Model (1) could be rewritten as which was the null model for the GWAS of main effect and epistatic effect detection. On this null model, we added a specific marker or marker pair to the model to test the putative effect.
To test the additive effect of marker k, we added Z k a k to the null mode so that the linear mixed model became: Frontiers in Plant Science | www.frontiersin.org Let e = ξ + ε , so that the model was rewritten as: The expectation of model (4) was E(y) = Xβ + Z k a k and the variance was: Let us perform eigenvalue decomposition for matrix K, K = UDU T , where U was the eigenvector matrix and D was a diagonal matrix holding the eigenvalues. So, The above linear mixed model was The variance of the transformed residuals was The expectation and variance of y * were E(y * ) = X * β + Z * k a k and var(y * ) = Iσ 2 . Therefore, model (6) became a simple linear model with a homogeneous residual variance. The conventional least squares method could be used to estimate the parameters and test for the marker effect. Since the model of the transformed phenotypic values was very simple, the "lm" function in R was applied to estimate the marker effect and test the significance of the marker.
Considering the dominance and epistatic effects, we adopted a more general likelihood ratio test (LRT) for a particular effect. The likelihood ratio test for the additive effect of marker k was where L 0 (β) was the likelihood value evaluated from the null model given in equation (7) below, and L 1 (β,α k ) was the likelihood value evaluated from the full model given in equation (6). The LRT statistic was eventually converted into the log of odds (LOD) score using LOD = LRT/4.61. If the intervals of different QTL were overlapped or the genetic distance of peak SNP of two QTL was within 0.65 cM (the average density in the whole genome), we called such QTL as a pleiotropic QTL (a QTL affecting more than one trait).
Dominance effect of marker k was detected using the same model as the additive effect except that Z k was replaced by W k . In the following, we called the significant additive and dominance QTL as the main effect QTL.
The additive-by-additive effect was detected by the following likelihood ratio test, where the null model was and the full model was Similarly, the additive-by-dominance effect was detected using The null model and the full model were and respectively. Similarly, the dominance-by-additive effect was detected using The null model and the full model were and respectively. Finally, the dominance-by-dominance effect was tested using The corresponding null model and full model were and respectively. LOD scores were converted the same way as we did for the additive effect. An empirical threshold of 2.5 for the LOD score was used to determine significance of an additive or a dominance effect. A LOD threshold of 5.0 was used to determine the significance of an epistatic effect (Churchill and Doerge, 1994;Xu, 2013). A confidence interval in the genome was determined for each detected QTL with the following steps: (1) all significant SNPs passing the threshold were selected; (2) the most significant SNPs were kept within a 10 cM interval; (3) the QTL interval was formed using a 1.5-LOD drop-off method (Broman, 2001). The names of QTL referred to McCouch's method (McCouch et al., 1997), and a dash (-) was added to designate different datasets.
The estimated additive and dominance effects for each QTL were extracted from the estimated regression coefficients (a k and d k ) from the models presented above. The proportion of the phenotypic variance explained (PVE) contributed by each QTL was calculated using (Utz et al., 2000;Garin et al., 2017), where RSS Full was the residual sum of squares of the full model and RSS Null was the residual sum of squares of the null model. We also performed QTL mapping in the RIL, TC and TM population separately. The model was the same as described above except that only the additive and additive-by-additive polygenic effects were used to control genetic background. QTL mapping for MPH was conducted using a similar linear mixed model to the original traits. Details of the MPH analysis can be found in a previous study (Jiang et al., 2017).

Genomic Selection
The genetic effects of single-cross hybrids can be dissected into additive, dominance and epistatic polygenic effects as mentioned before. Here, we only considered the first two components in the genomic prediction model. The linear mixed model was (Su et al., 2012;Xu et al., 2014).
where y was the phenotype vector, Xβ represented the fixed effect, ξ a was the additive polygenic effect with an assumed distribution of ξ a ∼ N(0, K a σ 2 a ), ξ d was the dominance polygenic effect with a distribution of ξ d ∼ N(0, K d σ 2 d ), K a was the additive kinship matrix and K d was the dominance kinship matrix.
Three genomic prediction schemes were proposed to mimic the scenarios in practical genomic hybrid breeding. Scheme (1), abbreviated as CV1: to predict the trait values for the TM population from the phenotypes and genotypes of the TC population or vice versa. Scheme (2), abbreviated as CV2: to select the hybrids sharing the same RILs in TC and TM population as the training set to predict the rest of the population. Scheme (3), abbreviated as CV3: to select the hybrids having the different RILs in TC and TM population as the training set to predict the rest of the population. Scheme (1) and (2) belong to the so-called T1 case, and scheme (3) is in the category of T2 (Technow et al., 2014;Zhao et al., 2015a). The three scenarios are illustrated in Supplementary Figure 1.
The across population prediction in scheme (1) was conducted using a model that contained only the additive polygenic effect. For schemes (2) and (3), the prediction models contained both the additive and the dominant polygenic effects. The predictions were implemented with the BGLR software package in R (Pérez and De Los Campos, 2014). The prediction accuracy was assessed with a two-fold cross-validation scheme. In each run, 1/2 of the lines were removed from the training set and then the correlation between the predicted values and the observed values of the removed lines was calculated. The two-fold cross-validation was repeated 200 times. In addition, significant SNPs were treated as fixed effects in the prediction model, which has been termed wGS (Bernardo, 2014;Würschum et al., 2018). For example, when the TC population was used to predict the TM population, the QTL detected in the TC population were treated as fixed effects in the linear mixed model used to predict the TM population. For schemes (2) and (3), QTL detected from the pooled population of TC and TM were treated as fixed effects included in the models to predict the rest of the population. For comparison, the additive model with kinship matrix inferred from RILs population was used to yield the prediction accuracy of 10 traits using a two-fold cross-validation scheme in TC and TM population, respectively. This process was repeated 200 times. Data visualization was done with the ggplot2 and ggpubr packages in R (Wickham, 2016).

Phenotypic Variation and Heritability in the Recombinant Inbred Line, TC, and TM Populations
The RIL population showed a larger variation for the 10 investigated traits than the TC and TM populations ( Table 1 and Figure 1). The genetic variance components were significant (p < 0.01) for all traits in the three populations. Except for the trait VW in the TM population, the variance of the genotype-by-environment interaction was also significant (p < 0.01) for all traits. The estimated broad-sense heritability ranged from 0.68 for VW to 0.95 for PH in the RIL population, from 0.57 for VW to 0.91 for PH and KT in the TC population, and from 0.60 for GY to 0.89 for PH in the TM population. In general, PH had the highest heritability and VW or GY had the lowest heritability. The obtained moderate to high heritability implied that the experimental designs and phenotyping procedures were appropriate and accurate.
The MPH for traits KT and VW was negative on average in both the TC and TM populations, which means that the hybrids often had phenotypic values lower than the mean of the two parents (Figure 1). In both the TC and TM populations, GY had the highest heterosis, followed by PH and EH (Supplementary Table 1), and the genetic variances of MPH were statistically significant for all traits. The heritability of MPH in the TC population ranged from 0.36 for VW to 0.81 for PH, while the heritability of MPH ranged from 0.27 for VW to 0.78 for PH in the TM population. The moderate to high heritability of MPH lay the foundation to dissect the genetic architecture of heterosis. For PH, EH, RNPE, KT, and KL the GCA variance of testers (σ 2 GCA Tester ) had larger values than the GCA variance of RILs (σ 2 GCA RIL ), which indicated that the testers played an important role in hybrid performance. The SCA/GCA ratios indicating a relative contribution of additive and nonadditive (dominance and epistasis) effects to phenotypic variation ranged from 0.04 for KT and 0.77 for GY (Supplementary Table 2). And for GY, the variance of SCA (σ 2 SCA ) was higher than both σ 2 GCA Tester and σ 2 GCA RIL , which was consistent with the large MPH variation in phenotype ( Figure 1J).

Trait Correlation and Path Analysis in the Three Populations
Relatively high correlations between traits were observed in the three populations (Figure 2A). The correlation between traits KNPR and GY was r = 0.69 (p < 0.01) and the correlation between KL and GY was 0.57 (p < 0.01) in the RIL population, which were the highest among the correlations between GY and the other traits. In the TC population, the highest correlations of GY occured between GY and HGW (r = 0.50, p < 0.01) and between GY and KL (r = 0.45, p < 0.01). In the TM population, the highest correlations were between GY and KNPR (r = 0.54, p < 0.01) and between GY and KL (r = 0.41, p < 0.01).
It is difficult to determine which trait contributes the most to the variation of grain yield only through correlation analysis between GY and the other traits. We therefore next performed a path analysis of all traits with GY (Figures 2B-D). By taking GY as the response variable and all other traits as independent variables, we estimated the path coefficients for every trait. In the RIL population, the highest path coefficients occurred for KNPR (0.60) and for HGW (0.36). The trait HGW had the highest path coefficient (0.79), followed by KNPR (0.54) in the TC population. In the TM population, the highest path coefficient was 0.79 for KNPR, followed by 0.62 for HGW. In summary, KNPR and HGW contributed most to the variation of grain yield.

Main Effect Quantitative Trait Loci Mapping in the Recombinant Inbred Line, TC, and TM Populations
A high-density genetic map was constructed using 4,141 bins, covering 2669.49 cM of the maize genome (Supplementary Table 3 and Supplementary Figure 2). The average density of the marker map was 0.64 cM/bin in the whole genome, enabling a high resolution for QTL mapping.
To dissect the genetic architecture of the 10 traits, we first examined the additive model with the additive polygenic effect plus the additive-by-additive polygenic effect to control the genomic background (Supplementary Table 4). The additive (narrow-sense) heritability in the RIL population ranged from 0.25 for VW to 0.69 for PH. In the TC population, it ranged from 0.31 for VW to 0.70 for RNPE and in the TM population, it ranged from 0.38 for VW to 0.72 for EH. Generally, the proportion of phenotypic variance explained by the additive effects was greater than that explained by the additive-by-additive effects for all 10 traits. We also found that the proportion of variance explained by the additive-by-additive effects for the traits RNPE, KT, KW, KL, and GY was larger in the RIL population than the corresponding proportion in the TC and TM populations (Supplementary Table 4), illustrating that further studies are needed to understand the non-additive genetic architecture of these traits.
We also mapped QTL for the 10 traits in the RIL, TC, and TM populations, respectively ( Figure 3A and Supplementary  Table 5). In the RIL population, a total of 16 QTL were identified on eight chromosomes and five superior alleles were from the Ye478 parental line over the 10 traits. In the TC population, a total of 18 QTL were identified, among which 10 superior alleles came from the Ye478 parent. 19 QTL were identified in the TM population, among which eight superior alleles originated from the Ye478 parent. Three common QTL were jointly identified in the RIL and TM populations, two QTL were shared by the RIL and TC populations, and six QTL were jointly detected in the TC and TM populations ( Figure 3B). One QTL (MQTL15) located in the interval from 162.93 to 172.23 Mb on chromosome 3 was shared among all three populations. This QTL was associated with KW in the RIL and the TC population and with HGW in all three populations (Figure 3C). A few other QTL also showed pleiotropic effects. For example, MQTL9 located between 1.71 and 4.67 Mb on chromosome 2 was associated with RNPE and KL in the TC population and with KNPR and KL in the TM population.
In the above QTL mapping results, low overlapping ratios among TC, TM and RIL populations were observed. In addition, in phenotype, the top 10 lines in the TC population did not match those identified in the TM population or vice versa (Supplementary Figures 3A,B). The two genetic phenomenons suggested that non-additive effects were important for hybrid performance. In this case, it is interesting to further dissect the hybrid performance to mine dominance and epistatic QTL.

Multiple Variance Components Dissection and Main Effect Quantitative Trait Loci Mapping for Hybrid Performance and Midparent Heterosis in the TC-TM Population
We dissected the contribution of all five variance components (additive, dominance, and three epistatic polygenic variances) by Bayesian generalized linear regression (Pérez and De Los Campos, 2014) based on the hybrid performance and MPH in the TC-TM population. The results for the hybrid performance showed that additive-by-additive was the most important polygenic effect for the traits PH, EH, and KT, additive-by-dominance was predominant for VW and dominance was the most important polygenic effect for the remaining traits (Supplementary Table 6). For the analysis of MPH, the additive-by-dominance variance contributed the most for traits KT, KW, VW, and HGW, while the dominance variance contributed the most for the other six traits (Supplementary Table 7). Different proportions of dominant variances among 10 traits showed the complexity of heterosis.
We implemented a mixed model to test the main (additive and dominance) effects of a specific marker for both the hybrid performance and the MPH for all traits in the pooled TC-TM population. A total of 21 main effect QTL were identified for the 10 traits for hybrid performance (Supplementary Table 5 and Figure 4A). Among them, one had a significant dominance effect for KNPR and was located in the interval 210.29-211.57 Mb on chromosome 2 ( Figure 4B and Supplementary Figure 4). For the other 20 QTL, the additive and dominance effects were confounded due to the fact that there were only two genotypes per locus (Supplementary Figure 4). Moreover, a total of 17 dominance QTL were detected for MPH for the 10 traits ( Figure 4C and Supplementary Table 5). Interestingly, only two detected QTL were in common between MPH and hybrid performance (Figure 4A). The pleiotropic QTL MQTL43 located in the interval around 80.08-112.87 Mb on chromosome 10 was associated with EH and GY in the MPH dataset and with PH, EH and GY in the TC-TM-Main dataset (Supplementary Table 5 and Figure 4C). The lack of common QTL between MPH and hybrid performance implies that the two phenomena might have different genetic architectures, consistent with the results of the variance component analysis.

Epistasis Plays an Important Role in Hybrid Performance
For hybrid performance in the TC-TM population, we scanned the entire genome to identify significant epistasis loci for the 10 traits and 197, 176, 131, and 112 significant epistatic pairs of loci were identified for additive-by-additive, additive-by-dominance, dominance-by-additive and dominanceby-dominance effects, respectively (Supplementary Table 8).
The number of significant locus pairs varied across traits and the proportion of explained variance of an epistatic interaction ranged from 3.46 to 4.52%. For grain yield, only one significant additive-by-dominance QTL were detected. We observed the phenomenon of a continuous region interacting with another locus in the genome. For example, for additive-byadditive mapping, the interaction between a cluster of adjacent SNPs on chromosome 8 (Chr8_180048590, Chr8_180913576, Chr8_181023046, and Chr8_180032314) and a locus on chromosome 6 (Chr6_166754537) was significantly associated with PH (Supplementary Table 8).
EH had a more simple genetic architecture compared to GY and the variation of MPH for EH was also higher (Supplementary Table 1). We therefore used the trait EH as an example to investigate the epistatic effects in the RIL and the two hybrid populations. In the RIL population, no QTL was identified ( Figure 5A). However, in the pooled TC-TM population, two main effect QTL for hybrid performance were identified on chromosomes 1 (TC-TM-Main-qEH1 represented by the peak SNP Chr1_131115160) and 10 (TC-TM-Main-qEH10 represented by the peak SNP Chr10_91890676) ( Figure 5B). For MPH, however, only the MPH-qEH10 QTL had a significant dominance effect as well as several additional small-effect QTL ( Figure 5C). We further tested the additive-by-additive interactions between TC-TM-Main-qEH1 and all other SNPs (4,140 in total). None of the tested effects were significant in the RIL population ( Figure 5D). However, several significant interactions were identified in the pooled TC-TM population ( Figure 5E). Further analysis confirmed the interaction between the two loci TC-TM-Main-qEH1 and TC-TM-Main-qEH10 in the pooled TC-TM population (Figure 5F). The specific type of epistatic effect between the two loci in the TC-TM population could not be determined because there were only two different genotypes at each locus. However, as we observed that the additive-by-additive effect was not significant between these two loci in the RIL population, we concluded that it is likely the additive-by-dominance or dominance-by-dominance effects that led to the detection of this epistatic QTL in the hybrid population but not in the RIL population.

Correlation Between the Number of Favorable Quantitative Trait Loci and Hybrid Performance
We chose a slightly lower significance threshold of LOD = 2.0 to obtain more loci for this analysis, which yielded four and six significant QTL for GY in the TC and TM population, respectively. If the performance of heterozygous genotypes was better than that of homozygous genotypes at one QTL, it was called a heterozygous favorable QTL; otherwise, it was called a homozygous favorable QTL. The correlations between the number of favorable QTL and the hybrid performance were calculated for all 10 traits (Supplementary Table 9). The correlations between the hybrid performance and the number of favorable homozygous QTL (r 1 ), the number of favorable heterozygous QTL (r 2 ) and the total number of favorable QTL (r 3 ) varied across traits, but were significant for most of traits in both the TC and TM populations. In the TC population, three of the four QTL for GY were heterozygous favorable QTL and r 1 , r 2 , and r 3 were 0.16 (Supplementary Table 9), 0.43 ( Figure 6A) and 0.41 (Figure 6B), respectively. In the TM population, only two of the six detected QTL were heterozygous favorable QTL and r 1 , r 2 , and r 3 were 0.42 (Figure 6C), 0.25 (Figure 6D), and 0.47 (Figure 6E), respectively. These results illustrate that superior hybrids can be selected by combining favorable alleles at significant loci.

Genomic Selection Accuracy in Different Breeding Schemes
For genomic selection within populations, the prediction accuracy ranged from 0.70 for RNPE to 0.40 for VW in TC and ranged from 0.63 for KT to 0.51 for GY in TM population ( Figure 7A). Generally, traits with a low heritability usually had a low prediction accuracy, like VW in TC and GY in TM population.
The three cross-validation schemes were shown in Supplementary Figure 1. For strategy 1 (CV1), when using the TC population to predict the TM population (TC_TM), the prediction accuracy ranged from 0.305 for GY to 0.699 for EH ( Figure 7B). Conversely, the prediction accuracy ranged from 0.287 for GY to 0.73 for EH when the TM population was used to predict the TC population (TM_TC). When the significant QTL identified in the training population were included as fixed effects in the prediction model, this did not result in an improvement of the prediction accuracy for most of the traits. Only a few traits, e.g., EH and KL, showed a slight improvement. For some traits, e.g., KNPR and KT, the prediction accuracy even decreased.
For the second cross-validation strategy (CV2), the lowest prediction accuracy was 0.49 for VW and the highest prediction accuracy was 0.90 for RNPE. For the third cross-validation strategy (CV3), the lowest prediction accuracy was 0.53 for GY and the highest prediction accuracy was 0.92 for KT and RNPE ( Figure 7C). The results also showed that the prediction accuracy of CV3 was higher than within population scheme and CV1, CV2, regardless of whether GS or wGS was applied. The wGS taking potential QTL as fixed had a higher prediction accuracy than GS in both CV2 and CV3 ( Figure 7C).

Hundred Grain Weight and Kernel Number per Row Significantly Contribute to the Variation of Grain Yield
Grain yield is a complex trait, affected by many genetic and non-genetic factors. The three traits that were found to mainly contribute to GY are HGW, RNPE, and KNPR. In earlier studies, the focus has been placed on correlation analysis between traits. In general, moderate to high correlations were observed between GY and many other traits (Cui et al., 2016;Li et al., 2020Li et al., , 2021. However, it is difficult to determine which trait contributes the most to grain yield. Path analysis is an alternative approach that allows examining the relative importance of a component trait to the variation of the target trait. Our results revealed that regardless of whether the population was the RIL population or the hybrid population, HGW and KNPR had the highest path coefficients thus, contributed the most to the variation of GY (Figures 2B-D). Consequently, HGW and KNPR are promising indirect traits to improve GY in hybrid breeding in maize.

Identification of Quantitative Trait Loci Using an Additive Genetic Model
Previous studies in the underlying populations (RIL, TC, and TM) focused on additive genetic models, where the genotypes and phenotypes of the RIL population were used to map QTL by treating the GCA as traits (Zhou et al., 2018;Lu et al., 2020). In this study, the genotypes and phenotypes of two hybrid populations were used to detect significant main (additive and dominance) effects and epistatic QTL. Only two common QTL were identified between the RIL and the TC hybrid population, and three common QTL were identified between the RIL and the TM hybrid population (Figure 3B). These results indicated that the non-additive genetic effect played an important role, which meant that a line with a moderate value of GY can still yield a high GY when crossing with testers (Supplementary Figure 3). And an additive model was not enough to explain heterosis.

Non-additive Polygenic Effects Play an Important Role in Hybrid Performance
The proportion of phenotypic variance explained by the additiveby-additive effects in the RIL population was higher for most traits than that in the TC and TM populations (Supplementary Table 4). And the unparallel relationship between RILs and hybrid populations (Supplementary Figure 3) inspired further exploration of the genetic basis in the combined TC-TM hybrid population using a model integrating non-additive polygenic effects. Based on the estimated variance components, we conclude that the prominent gene action varies across traits (Supplementary Table 6). For example, the prominent variance was the additive-by-additive component for the trait PH, but for GY the prominent variance was the dominance component. This presents a considerable challenge for selecting elite singlecross hybrids and for uncovering the importance of non-additive genetic effects because additive variance is not a prominent factor to control the variation of any of the investigated traits.

Identification of Quantitative Trait Loci Using a Non-additive Model
For genetic dissection of the hybrid performance, Xu (2013) proposed a new mixed model method for QTL mapping by incorporating multiple polygenic covariance structures, which consist of the additive, dominance, and epistatic variance components. In theory, each particular effect could be tested in a model by controlling all other genetic effects as background. Similarly, a quantitative genetic framework was proposed for the genetic dissection of MPH (Jiang et al., 2017). The above two linear mixed models are very similar in the polygenic background control. The main difference is the response variable, in the former it is the performance of the hybrid (Xu, 2013) and in the latter the MPH (Jiang et al., 2017). It should be noted that the hybrid performance could not be simply replaced by MPH with just the removal of the additive effect from the linear mixed model (Jiang et al., 2017). In our study, the GBS technology only covered about 0.07-fold of the genome in our populations. The two tester lines had the same genotype at 95% of the loci, which resulted in the pooled TC-TM population having just two genotypes at 95% of loci. Two genotypes per locus mimick a backcross population so that the additive effects are confounded with the dominance effects, which explains why the TC-TM-Add model was nearly the same as the TC-TM-Dom model (with one exception for a dominance QTL) (Supplementary  Figure 4). Regardless of the high similarity between the hybrid populations and the hypothetical BC population, the hybrid FIGURE 7 | Genomic prediction accuracy of different cross-validation strategies. (A) Prediction accuracy within Chang7-2 × RIL (TC) and Mo17 × RIL (TM) population. RIL, the recombinant inbred line population developed by Ye478 × Qi319. (B) Prediction accuracy of 10 traits of cross-validation strategy 1 (CV1). TC_TM represents TM predicted by the TC population; TM_TC represents TC predicted by the TM population. (C) Prediction accuracy of 10 traits of cross-validation strategy 2 (CV2) and cross-validation strategy 3 (CV3). wGS is a weighted genomic selection that incorporates the peak single nucleotide polymorphisms (SNP) of potential target quantitative trait loci as fixed effects. PH, plant height; EH, ear height; RNPE, row number per ear; KNPR, kernel number per row; KT, kernel thickness; KW, kernel width; KL, kernel length; VW, volume weight; HGW, hundred grain weight; GY, grain yield per plant.
populations still had advantages. In a dominance test for MPH, 17 significant dominance loci were detected ( Figure 4A) and for hybrid performance, a set of significant epistatic loci was identified (Supplementary Table 8).

Improve Prediction Accuracy by Integrating Functional Markers
In genomic selection models like GBLUP or rrBLUP, all SNPs were treated equally or had the same distribution when treated as random. Actually, significant QTL contributed more to the variation of traits. In such a case, significant SNPs should be treated differently. In this study, those SNPs were included in the fixed effect in the GS model to explore whether prediction accuracy could be improved. For cross-validation scheme 1, just two of 10 traits, namely EH and KL showed slight improvement. We guess that TC and TM populations had different genetic backgrounds  and those QTL made different effect within two different populations. So the QTL effect was estimated biasedly when only one tester population as training population. But for CV2 and CV3, the population used for QTL mapping consisted of lines from both TC and TM, in this case, improvements were observed for all traits harboring QTL ( Figure 7C). By comparison, it showed that CV3 scheme had superiority over the within population scheme, CV1, and CV2, which inspired us that when there were some known functional QTL in a target population, a strategy treating known QTL as fixed effects with CV3 design was a better choice for genomic prediction.

Relationship Between Midparent Heterosis and Hybrid Performance
Hybrid performance is the phenotypic value of the hybrid, which is the sum of the midparent heterosis and the midparent value. Hybrid performance is controlled by the additive, the dominance and all four epistatic polygenic effects, whereas MPH is not affected by the additive effect because the additive effect does not contribute to heterosis (Jiang et al., 2017). In this study, we confirmed the different genetic architecture of hybrid performance and MPH as both had only two QTL in common (Figure 4A), which is consistent with a previous study (Hua et al., 2003). Furthermore, the variance component ratios were also different between hybrid performance and MPH (Supplementary Tables 6, 7). In a wheat study, the midparent value showed a negative correlation with MPH but was positively correlated with the hybrid performance (Boeven et al., 2020). In our study, we observed the correlation between the hybrid performance and MPH was 0.77 (p < 0.01) (Supplementary Figure 5A), and a positive correlation of 0.23 (p < 0.01) between the hybrid performance and midparent value (Supplementary Figure 5B), while a negative correlation between midparent value and MPH for grain yield of −0.45 (p < 0.01) (Supplementary Figure 5C). And the path analysis also highlighted the superior contribution of MPH to hybrid performance in hybrid population (Supplementary Figure 5D). In plant hybrid breeding, we aim to select a single-cross hybrid with both high MPH and midparent value, but these seem to be contradictory goals considering the negative correlation. Consequently, hybrid breeding must balance the two and target hybrid performance to achieve high performing hybrids.

Mechanisms of Midparent Heterosis and Hybrid Performance
Although dominant and additive effects couldn't be separated and dominant degree couldn't be estimated in this study, multiple variance components dissection provided possibility to assess the mechanism of heterosis and hybrid performance. Results showed the dominance contributed the highest proportion for MPH of most traits, especially for GY and KNPR (Supplementary Table 7). However, it was found that the epistasis (sum of additive-by-additive, additive-by-dominance, and dominanceby-dominance) contributes the highest proportion to hybrid performance of GY, PH, EH, and KNPR (Supplementary Table 6). The results were similar to a previous report in maize (Tang et al., 2010). A series of linear mixed models incorporating multiple polygenic covariance structures together with NCII population provide possibility to explore the genetic factors and mechanism of heterosis.

DATA AVAILABILITY STATEMENT
The original contributions presented in the study are included in the article/Supplementary Material, further inquiries can be directed to the corresponding authors.

AUTHOR CONTRIBUTIONS
WL, ML, and SX designed the study. ZZ collected the phenotypic data. DL, XiaL, GL, JL, and HW performed data analysis. DL drafted the manuscript. DL, XiaL, YJ, SC, TW, JR, XinL, SX, ML, and WL revised the manuscript. All authors read and approved the final manuscript.