## Abstract

A new mixed-model method was developed for mapping quantitative trait loci (QTL) by incorporating multiple polygenic covariance structures. First, we used genome-wide markers to calculate six different kinship matrices. We then partitioned the total genetic variance into six variance components, one corresponding to each kinship matrix, including the additive, dominance, additive × additive, dominance × dominance, additive × dominance, and dominance × additive variances. The six different kinship matrices along with the six estimated polygenic variances were used to control the genetic background of a QTL mapping model. Simulation studies showed that incorporating epistatic polygenic covariance structure can improve QTL mapping resolution. The method was applied to yield component traits of rice. We analyzed four traits (yield, tiller number, grain number, and grain weight) using 278 immortal F2 crosses (crosses between recombinant inbred lines) and 1619 markers. We found that the relative importance of each type of genetic variance varies across different traits. The total genetic variance of yield is contributed by additive × additive (18%), dominance × dominance (14%), additive × dominance (48%), and dominance × additive (15%) variances. Tiller number is contributed by additive (17%), additive × additive (22%), and dominance × additive (43%) variances. Grain number is mainly contributed by additive (42%), additive × additive (19%), and additive × dominance (31%) variances. Grain weight is almost exclusively contributed by the additive (73%) variance plus a small contribution from the additive × additive (10%) variance. Using the estimated genetic variance components to capture the polygenic covariance structure, we detected 39 effects for yield, 39 effects for tiller number, 24 for grain number, and 15 for grain weight. The new method can be directly applied to polygenic-effect-adjusted genome-wide association studies (GWAS) in human and other species.

EPISTATIC effects refer to interactions of allelic effects between different loci. Depending on the natures of populations, the variance of epistatic effects for a quantitative trait may be partitioned into several different types of variance components, *e.g.*, additive × additive, additive × dominance, dominance × additive, and dominance × dominance (Cockerham 1954). The relative importance of each variance component usually varies across different traits. Accurate partitioning of these epistatic variance components can help understand the genetic mechanisms of complex traits and develop more efficient breeding programs. Prior to the genome era, special designs of experiments were required to partition the epistatic variance into different variance components, *e.g.*, the NC design III (Cockerham and Zeng 1996). The well-known Cockerham’s epistatic model (Cockerham 1954) is the principle basis for data analysis. With molecular data, marker genotypes can be observed and thus complicated designs are no longer crucial. For example, the recombinant inbred line (RIL) design allows us to estimate the additive × additive effects between pairs of markers. The F2 design facilitates estimation of all the four types of epistatic effects. These simple designs of experiments are very popular in crops and laboratory animals.

Quantitative trait locus (QTL) mapping for epistatic effects has attracted much attention from geneticists (Cockerham and Zeng 1996; Kao and Zeng 2002; Yi and Xu 2002; Xu 2007; Xu and Jia 2007; Garcia *et al.* 2008). The main hurdle in epistatic effect QTL mapping is the large number of interaction effects to estimate. For low-density marker maps, all pairwise interactions can be fit simultaneously into a single model (Yi and Xu 2002; Xu 2007). With high-density marker maps, simultaneous fit of all genetic effects to a single model is prohibited. An *ad hoc* method is to fit one pair of loci at a time and the entire genome is then scanned in a two-dimensional approach (Bell *et al.* 2011). It is well known in genome-wide association studies (GWAS) that such a genome-wide scanning has ignored the polygenic effect and thus will inflate the residual error variance (Yu *et al.* 2006; Zhang *et al.* 2010; Zhou and Stephens 2012). As a consequence, the statistical power may be low. When the pedigree relationship of individuals is known, one can fit a polygenic effect to the model using the standard linear mixed-model approach (Henderson 1975), as done for the inbred lines of maize (Zhang *et al.* 2005). However, the pedigree relationship is rarely known in a randomly sampled population. Yu *et al.* (2006) proposed using marker-inferred kinship matrix to capture the polygenic variance. This approach has now become the standard method for GWAS. Efficient algorithms have been proposed afterward to improve the computational efficiency (Kang *et al.* 2008; Zhang *et al.* 2010; Lippert *et al.* 2011; Zhou and Stephens 2012). However, such a polygenic mixed linear model has not been investigated in GWAS for epistatic effects.

QTL mapping (genome-wide linkage studies) is slightly different from GWAS because the target populations for the two approaches of genetic mapping are often different. However, if analyses are focused on markers only (ignoring pseudo-markers between observed markers), the statistical methods are much the same for the two approaches except that the random population GWAS requires an extra term in the model to control for population structures. The linear mixed-model GWAS incorporating marker-inferred kinship matrix as the covariance structure can be naturally adopted here for QTL mapping to capture the background genetic information. For a well-controlled mapping population, *e.g.*, the F2 population, the marker inferred kinship matrix is equivalent to the identical-by-descent (IBD) matrix of polygenes. Using this matrix to model the covariance structure plays a similar role to the cofactors added to the interval mapping, the so-called composite interval mapping as proposed by Jansen (1994) and Zeng (1994). A more appealing feature of this kinship-corrected QTL mapping approach is that it is easy to control the genetic background information caused by epistatic effects. The original composite interval mapping, in principle, can take into account pairwise interaction terms as cofactors, but no attempt has been made because of the difficulty of choosing the pairwise cofactors. If the covariance structure includes kinship matrices for various epistatic effects, the genetic background due to epistatic effects can be controlled.

Motivated by the above arguments, we propose a mixed-model GWAS for epistatsis using ultra-high-density SNP markers. Xie *et al.* (2010) and Yu *et al.* (2011) released a high-density SNP map containing about 270,000 SNPs. From the SNPs, they inferred recombination breakpoints for 240 RILs of a rice population derived from two inbred lines. The authors further combined consecutive SNP markers of cosegregation into bins. Within a bin, there are no breakpoints and thus all SNPs within the same bin have exactly the same genotypes. The >270,000 SNPs were eventually converted into 1619 bins. Each bin was then treated as a synthetic marker. Genetic analysis was then performed on the bin genotypes. Using the binned genotypes, Yu *et al.* (2011) conducted QTL mapping for seven yield component traits in rice and identified numerous QTL. Using the 240 RILs, Hua *et al.* (2002) and Hua *et al.* (2003) created 360 crosses to dissect epistatic effects and to investigate the genetic basis of heterosis for yield component traits. Most recently, Zhou *et al.* (2012) reinvestigated the genetic basis of heterosis using 278 of the 360 crosses. They called the RIL-generated crosses immortal F2 (IMF2) because the genotypes of the crosses mimic the genotypes of the regular F2 population. The IMF2 crosses are also called recombinant inbred intercrosses (RIX) in laboratory mice derived from crosses of multiple RILs (Zou *et al.* 2005). Zhou *et al.* (2012) used the 1619 bins as synthetic markers to test all pairwise interaction effects using a two-dimensional scanning approach by fitting one pair of bins at a time. The purpose of their study was to evaluate the relative importance of dominance and epistasis to heterosis.

In this study, we analyzed the same IMF2 population using genotypes of the 1619 bins. We first evaluated the relative importance of each type of epistatic effects on four yield component traits using variance component analysis and then mapped epistatic effects using marker-generated kinship matrices to control the epistatic genetic background.

## Material and Methods

### Plant materials and SNP bin map

An IMF2 population of rice was investigated in this study. The IMF2 population was generated by Dr. Qifa Zhang’s laboratory and previously described in detail by Hua *et al.* (2002) and Hua *et al.* (2003). It consisted of 360 crosses made by random matches of 240 RILs derived by single-seed descent from a cross between Zhenshan 97 and Minghui 63. Field data of yield (YIELD), number of tillers per plant (TILLER), number of grains per panicle (GRAIN), and thousand grain weight (KGW) were collected in the 1998 and 1999 rice growing seasons from replicated field trials on the experimental farm of Huazhong Agricultural University, Wuhan, China. Over a quarter million high-density SNP markers were used to infer recombination breakpoints (crossovers), which were then used to construct bins. The bins were treated as “new markers” for association studies. The number of bins in this rice population was 1619. The bin map was constructed by genotyping the RILs with population sequencing (Xie *et al.* 2010; Yu *et al.* 2011). Among the 360 crosses, only 278 of them were available in both phenotypes and bin genotypes. Therefore, the bin genotype data were stored in a matrix. Each genotype was coded by A for the Zhenshan 97 genotype, B for the Minghui 63 genotype, and H for the heterozygote. The genotypes and phenotypes were downloaded from the journal website of the study by Zhou *et al.* (2012). For each trait, there were two temporal replications (years 1998 and 1999). The phenotypic values of the two replicates were pooled for each cross after removing the year effects using (1)where and are the mean values of the trait measured in 1998 and 1999, respectively. This pooled trait value was treated as the actual phenotypic value for analysis. Apparently, we ignored the genotype by year interaction (G × E) effects, if there is any.

### Polygenic variance component analysis

We first numerically coded the genotype of individual *j* at bin *k* into two variables, (2)where and represent the codes for additive and dominance effects, respectively. Let *y* be an column vector for the phenotypic values of all the IMF2 crosses. The full epistatic model for bins is described by(3)where represents some nongenetic effects, *e.g.*, intercept in this case, and is a vector of residual errors with a normal distribution of unknown variance . In addition, represents element-wise vector multiplication, and and are the additive and dominance effects, respectively, for bin *k*. The four terms, , , and are the additive × additive, dominance × dominance, additive × dominance and dominance × additive effects, respectively, between bins *k* and . These four terms are called the epistatic effects. The total number of genetic effects in the above model is (4)which is beyond the capability of any models currently available in genome-wide association studies in a simultaneous manner. By treating each genetic effect as a randomly distributed normal variable with mean zero and a common variance across all bins or pairs of bins, the model becomes a mixed model. Let , ,, , and be the distributions. The fact that the variance of each type of genetic effects does not depend on the bins or bin pairs implies that all effects of the same kind share a common genetic variance. These common variances are called polygenic variances. The expectation of *y* is and the variance–covariance matrix of *y* is (5)where the *K*’s are marker-generated kinship matrices; *i.e.*, they depend on and . Formulas to calculate these kinship matrices are given in Table 1. Note that the kinship matrices for epistatic variances presented here are different from the one given by Ober *et al.* (2012), who used the Hadamard square of the additive kinship matrix as the additive × additive kinship matrix. They did not detect any epistatic variance in their *Drosophila melanogaster* population, likely due to the use of wrong additive × additive kinship matrix. Given these marker-generated kinship matrices, the variance components were estimated using standard mixed-model procedures (Henderson 1975).

One may rewrite the variance–covariance matrix in Equation 5 by (6)where is the ratio of the additive genetic variance to the residual variance and other λ’s are defined similarly. A more useful measurement of the relative importance of a type of genetic effect, *e.g.*, the additive effect, is (7)where the denominator (the sum of all variance components) is called the phenotypic variance. The above ratio is called the narrow-sense heritability. The broad-sense heritability is define by (8)Note that the broad-sense heritability using marker-generated kinship is often close to unity when the marker density is sufficiently high. Therefore, it is not an informative piece of information. However, the proportion of each variance component relative to the total phenotypic variance is informative and is the focus of discussion in this study.

### Model for genome scanning

Individual bin effects and bin pair interaction (epistatic) effects can be estimated after fitting the polygenic effects. The idea of the mixed-model GWAS (Yu *et al.* 2006) can be adopted here. The polygenic term remains in the model except that we used six kinship matrices to describe the covariance structure. The genetic effects of each bin and the epistatic effects of each pair of bins are modeled like the usual GWAS (Yu *et al.* 2006) with the polygenic structure included. The epistatic model appears like (9)for loci *k* and , where ξ is the polygenic effect (the sum of all genetic effects across the entire genome). There are two additive effects, two dominance effects, and four epistatic effects (denoted by double symbols with double subscripts) for each pair of bins. A complete scan under the epistatic model requires a two-dimensional search for all pairs of bins. We noted that the additive and dominance effects would be estimated multiple times. The inclusion of polygenic effects allows us to search for bin effects and epistatic effects separately. The model that includes the additive and dominance effects is (10)A complete scan for this model (10) requires *m* calls of the above model and two genetic effects are estimated under each call. The epistatic effect model excluding the additive and dominance effects is (11)A complete scan for this model (11) requires calls of the above model and four epistatic effects are estimated under each call. The mixed model can be written generically as (12)where is a general error term with and (13)Matrix *K* is a weighted sum of all effect-specific kinship matrices, it is not the kinship matrix given by Yu *et al.* (2006) in the original GWAS study, (14)and the weights are the ratios of the variance components to the residual variance. For the additive and dominance effect model, and . For the epistatic effect model, (15)where the special notation *a*||*b* indicates column concatenation of matrices *a* and *b*. The vector of genetic effects γ has been treated as fixed effects in the classical GWAS (Yu *et al.* 2006; Zhou and Stephens 2012). In this study, we treat γ as random effects with a multivariate distribution. The shared common variance is estimated from the data. When is set to ∞, the genetic effects γ become fixed effects.

### A two-step approach to genome scanning

#### Polygenic analysis:

The mixed-model analysis would include seven genetic variance components (six polygenic variances and one genetic variance for each bin or bin pair). The six polygenic variances change very little across bins or bin pairs. Therefore, we adopted the two-step approach for parameter estimation to save computational time (Zhang *et al.* 2010; Lippert *et al.* 2011). In the first step, we estimated only the polygenic variances using the model (16)where (17)This model allows us to estimate the six variance ratios, , which are needed to construct the *K* matrix. The *K* matrix estimated from this model is then treated as a constant (known value) in the second step. The polygenic model was treated as the null model for the next step genome scan. The likelihood function evaluated under the polygenic model is denoted by in the subsequent likelihood ratio tests.

#### Genome scanning:

We proposed the following eigen decomposition algorithm that does not require repeated inverse of matrix . Eigen decomposition for GWAS was originally proposed by Lippert *et al.* (2011) and later extended by Zhou and Stephens (2012). The eigen decomposition for the kinship matrix is , where *U* is the eigenvector (an matrix) and is the eigenvalue (a diagonal matrix). One property of the eigen decomposition is , which has avoided matrix inversion because is simply the inverse of a diagonal matrix. The error variance is rewritten by (18)Before calling the mixed-model procedure, we performed data transformations so that the mixed model is rewritten by (19)Let be the transformed phenotypic values. The model for the transformed data is (20)The expectation and variance matrix of the transformed data are and (21)respectively. Let (22)be a weight matrix. The residual variance structure for the transformed data are (23)The mixed-model equations for BLUE of β and BLUP of γ are (24)The null hypothesis may be tested in one of two ways. One way is the likelihood ratio test under the null model . The other way is the Wald test under the null model . The mixed-model equations provide both the BLUP of γ and the variance–covariance matrix of the BLUP estimate. The Wald test statistic is (25)When is true, the Wald test statistic follows a chi-square distribution with 2 degrees of freedom for the additive-dominance model and 4 degrees of freedom for the epistatic effect model.

### Software implementation

All data analyses were implemented using the restricted maximum-likelihood (REML) method (Patterson and Thompson 1971) implemented via the MIXED and HPMIXED procedures in SAS (SAS Institute 2009b). PROC IML (SAS Institute 2009a) was also used to calculate eigenvalues and eigenvectors and to perform data transformations. The polygenic variance component analysis was conducted using PROC MIXED. Genome scans were implemented using PROC HPMIXED, which is a specialized high-performance version of the MIXED procedure. The SAS codes are provided in Supporting Information, File S11 for PROC MIXED and File S12 for PROC HPMIXED. The six types of kinship matrices are given in File S1. The phenotypes of four quantitative traits are given in File S2. Genotypes of two bins (bins 729 and 1064) are given in File S3.

## Results

### Results of a simulation experiment

Before analyzing the rice data, we wanted to show the difference of genomic scans between the model including the polygenic variances and the model without the polygenic variances using simulated data. Similar comparison has been performed in the traditional mixed-model GWAS that ignored the epistatic polygenic variance (Yu *et al.* 2006). We used the genotypes of 217 bins on the first chromosome for the 278 IMF2 crosses of rice as the simulated genotypes and added some effects on bins and bin pairs to generate genetic values of all crosses. We then added a normally distributed residual error to the genetic value of each cross to form a phenotypic value for the cross. To simplify the model, we simulated only the additive effects and the additive × additive effects. For 217 bins, the model included 217 additive effects and additive × additive effects. Therefore, the polygenic model contained two genetic variance components.

In the first simulation experiment, we added only one additive effect on bin 62 and one additive × additive effect on bin pair (112, 182). The true effects and estimated effects are illustrated in Figure 1, red and blue, respectively. Figure 1, A and B, shows the main effects and the epistatic effects, respectively, under the model that ignored the polygenic covariance structure. Figure 1, C and D, shows the main effects and the epistatic effects, respectively, under the model that incorporated only the additive covariance structure. The improvement can be seen for the main effects but not the epistatic effects. Figure 1, E and F, shows the main effects and the epistatic effects, respectively, under the model that incorporated both the additive and the additive × additive covariance structures. Further improvement can be seen, especially for the epistatic effects where the peak of the estimated effects is much sharper than the previous two models. For the simple experiment, the model ignoring polygenic covariance structure already performed well.

In the second simulation experiment, we added 9 additive effects and 13 additive × additive interaction effects. The true effects and the estimated effects are illustrated in Figure 2, red and blue, respectively. The improvement of the models including the polygenic effects over the model ignoring them is more obvious in the complicated experiment. For the model ignoring the polygenic covariance structure (Figure 2, A and B), peaks of the estimated main effects do not correspond to the true values at all. The estimated epistatic effects correspond to the true effects when the true effects are large. When the additive covariance structure is incorporated in the model, improvement for the main effects is obvious (Figure 2C) but the epistatic effects have very little improvement (Figure 2D). Figure 2, E and F, show the main and epistatic effects, respectively, for the model incorporated both the additive and additive × additive covariance structures. The improvement is much clearer compared with the previous two models. The simulation experiments by no means were exhausted but they provided visual evidence that the polygenic-effect-adjusted model performed better than the model without the polygenic effects. Incorporating epistatic covariance structure led to further improvement.

### Polygenic variance analysis for yield component traits in rice

Results of the variance component analyses are summarized in Table 2. The variance ratios () are also listed in the table and they were used to calculate the weighted sum of different kinship matrices. The most important information from the table is the proportion of the phenotypic variance contributed by each type of genetic variance (). Several general conclusions were drawn from the observations: (1) None of the four traits were contributed by the dominance variance; (2) the additive variances contributed the most in two of the four traits (GRAIN and KGW); (3) the additive × dominance variance is important for YIELD and dominance × additive variance is important for TILLER; (4) the KGW trait was almost exclusively controlled by the additive variance; and (5) all the four traits had high broad sense heritability with KGW being the highest (). The goodness of fit (squared correlation between predicted and observed phenotypic values) was almost perfect for all traits. The likelihood ratio test (LRT) statistics for the polygenic effect model (including all six genetic variance components) were all high (37.9, 58.2, 82.2, and 248.9) with extremely small *P*-values (1.17*E*–06, 1.04 *E*–10, 1.25 *E*–15, and 7.04 *E*–51) for the four traits. The significant polygenic variance for each of the four traits implied that polygenic-effect-adjusted model should perform better than the model ignoring the polygenic effects. Note that the goodness of fit is not exactly the same as the broad sense heritability because the two are calculated using different formulas, in which was calculated from the data and was calculated from the sum of all variance components. The polygenic analysis provides only a general picture of the relative importance of each type of variance components. When a variance component is zero, it means that this type of effect is not as important as some other types of effects in general. It does not exclude the possibility that some bins may have significant effects of this type, which have been “diluted” in the polygenic analysis (see more discussion in the last section). File S11 gives the SAS code of PROC MIXED for the polygenic analysis.

### Genome scans for yield component traits in rice

#### Main effect analysis:

We used two models to perform the main effect analysis. One model contains the additive and dominance effects of each bin plus the polygenic effect captured by the polygenic covariance structure (all six types of polygenic variances were considered). The method is called weighted mixed-model analysis. The second model contains the additive and dominance effects for each bin but ignores the polygenic covariance structure. This method is called unweighted mixed-model analysis. The comparison provides a real life demonstration of the differences between the two models. The LRT statistic was used to declare statistical significance for each bin. The critical values for the LRT test statistics were determined via multiple permutations to adjust for multiple tests (Churchill and Doerge 1994). The polygenic covariance structures were reserved. The number of reshuffled samples in the permutation tests was 1000 per trait within each method. The maximum LRT was recorded for each permuted sample. For a total of 1000 permuted samples, we obtained 1000 such maximum LRT values. These maximum LRT values form an empirical null distribution. The LRT value of each bin from the original data analysis (not from the reshuffled data) was then compared to the null distribution. The *P*-value is the proportion of the 1000 permuted samples whose LRT values are greater than the LRT of the original data analysis. First, let us evaluate the critical values of the test statistics from the null distributions. The 90, 95, 99, and 100 percentiles of the null distributions are listed in Table 3 (the top and middle portions of the table). The critical values for the unweighted method are always higher than the weighted method. For example, the critical value for YIELD at type I error of 0.05 is 3.6711 for the weighted method, but the corresponding value for the unweighted method is 10.0605. The critical values also vary across different traits, indicating that permutation tests are necessary.

Figure 3 gives the LRT profiles for each trait under the two models. The weighted method usually has low LRT but with sharper peaks than the unweighted method. Due to lower empirical critical values for the weighted method, the powers are not necessarily low. In fact, many peaks co-occur between the two methods. For example, the peaks on chromosomes 6, 7, and 8 for YIELD are very similar between the two methods. Differences between the two methods are also obvious. Taking trait YIELD, for example, the weighted mixed model shows a sharp peak in the beginning of chromosome 12 (*P* < 0.05) but the unweighted method shows a much wider peak in the middle of the chromosome and this wide peak is not significant (*P* > 0.05). Because the peaks for the unweighted method are usually wider, more bins are significant but these significant bins are clustered due to linkage disequilibrium. For example, about 50 bins covering one-third of chromosome 3 for trait GRAIN are significant for the unweighted method due to high linkage disequilibrium, but the weighted method has narrowed down the peak into a small region consisting of only 10 significant bins. The estimated effects and the test statistics along with the *P*-values of the 1619 bins for all the four traits are provided in File S4 and File S5 for the weighted and unweighted methods, respectively.

We now focus only on the weighted mixed-model analysis and report all bins that passed a preliminary screen to eliminate all bins with LRT less than 4.61, equivalent to LOD score less than one. All bins with *P*-values <0.05 were declared as significant at the genome-wide type I error of 0.05. Using the *P*-value <0.05 criterion, the number of significant bins are 25, 4, 32, and 23 for traits YIELD, TILLER, GRAIN, and KGW, respectively. These numbers are not the numbers of detected QTL because consecutive bins tended to show cosegregation. Therefore, the actual numbers of QTL are less than the number of significant bins. For example, among the four significant bins (bins 534, 1004, 1005, and 1262) for trait TILLER, only three QTL can be declared because bins 1004 and 1005 are counted as one due to their high linkage. File S6 lists all the bins passing the preselection, including the significant bins drawn from the permutation tests for the four traits. The estimated numbers of QTL are presented later using another round of screening to remove redundant bins caused by close linkage.

#### Epistatic effect analysis:

Due to limitation of computing power, we used only the polygenic-effect-adjusted model (weighted method) to analyze the yield component traits. The LRT statistic was used to declare statistical significance for bin × bin interactions. The critical value for the LRT statistics was determined via a permutation analysis (Churchill and Doerge 1994) to adjust for multiple tests. For the epistatic effect model, we first eliminated all bin pairs with LRT <9.22, equivalent to LOD score <2. We then used the survived bin pairs to perform permutation analysis. The maximum LRT test statistics of 1000 permuted samples form an empirical null distribution. The 90, 95, 99, and 100 percentiles of the permutation-drawn null distributions for the four traits are given in Table 3 (bottom). These critical values also vary across traits. Trait KGW had the highest critical values and trait TILLER had the lowest critical values. The percentiles for the epistatic effect model are substantially higher than the percentiles for the main effect model. For example, the LRT critical value of KGW at 0.05 type I error is 5.47 for the main effect model but it is 15.77 for the epistatic effect model. The *P*-values were drawn for each pair of bins and the significant bin pairs along with their estimated effects are listed in File S7. The numbers of significant bin pairs are 1071, 205, 573, and 66 for traits YIELD, TILLER, GRAIN, and KGW, respectively. Again, these numbers are not the estimated numbers of epistatic effects due to tight linkage of consecutive bins. The estimated number of QTL interactions is addressed below using another cycle of screening.

For trait KGW in the epistatic effect model analysis, bin pair had an LRT of 18.23049 with a *P*-value of 0.016. We now provide a detailed analysis for this pair of bins. We used a model that includes the additive and dominance effects for both bins and all the four epistatic effects. The six polygenic covariance structures were also included in the model. This revised model now contains eight genetic effects and an intercept. Reanalysis of this bin pair produced eight estimated genetic effects, which are further converted into nine genotypic values using the transformation (26)where the vector in the right-hand side stores the estimated genetic effects in the following order: . The vector in the left-hand side gives the nine genotypic values. These genotypic values along with the marginal values are arranged in the following matrix:

For bin , the marginal effects are not significantly different among the three genotypes. However, the differences are very significant under the background. This detailed analysis shows the importance of the epistatic effects between the two bins for trait KGW. The two bins jointly contribute a genetic variance of 0.2991, which explains of the total trait variance. This is a significant contribution from only two loci for a quantitative trait. File S12 gives the SAS code of PROC HPMIXED for the analysis of the two bins. The genotypes of the two bins are given in File S3.

### Numbers of bins and bin × bin interactions

This step provides the last screen of the bins and bin × bin interactions that have survived all previous selections. We noted that due to close linkages between consecutive bins, multiple bins are detected in a wide range of chromosome regions. In the original study of the IMF2 population by Zhou *et al.* (2012), the authors set up a subjective criterion to eliminate the extra significant bins. For example, they pooled several bins in the neighborhood of a few centimorgans of the genome into a single cluster. Each cluster was treated as one QTL. The authors realized that this approach is subjective and may eliminate some true effects. Here, we adopted a different and more objective strategy to eliminate those superfluous bins. Since the numbers of significant bins and bin × bin interactions are substantially smaller than the total number of bins and bin × bin pairs available in the data, we can fit all these significant bins and bin pairs simultaneously in a single model. For example, there are 23 significant bins (*a* and *d*) and 66 significant bin pairs (*aa*, *dd*, *ad*, and *da*) for trait KGW, and simultaneous analysis of all effects with a population size of 278 is possible, provided that a penalized regression is used to overcome the overfitting problem. For trait YIELD, the total number of significant bins and bin × bin pairs is 25 + 1071 = 1096, leading to effects, which is even harder to deal with than the KGW trait. Therefore, we used a penalized regression approach to perform parameter estimation, the Lasso method (Tibshirani 1996) implemented via the GlmNet/R package (Friedman *et al.* 2010). The method selects effects with nonzero effects (variable selection). After a minor modification of the Lasso method, the program also provides a test for each effect included in the model (Xu 2013).

The selected (nonzero) effects and their *P*-values obtained from the Lasso method are given in File S8 for all four traits. In this data set, if bin1 and bin2 have the same number, the effect represents a main effect; otherwise, the effect is an epistatic effect. The variable named “type” provides information about the type of effect. For example, type = “*aa*” indicates the additive × additive effect. Significant effects are labeled 1 in the last column named “significant.” The numbers of significant effects sorted by effect type for the four traits are given in Table 4, where the top shows the number of (nonzero) effects selected by the Lasso method and the bottom gives the corresponding numbers that are significant at *P*-value <0.05. There are 39 significant effects apiece for YIELD and TILLER. Traits GRAIN and KGW have 24 and 15 significant effects, respectively. The number of significant epistatic effects is often higher than the number of main effects due to the substantially larger number of bin × bin pairs than the number of bins. The only exception occurs for trait KGW, which has 8 significant main effects and 7 significant epistatic effects. Table 5 shows the relative contributions of different types of effects to the total phenotypic variance. The significant effects collectively contributed 0.83, 0.69, 0.68, and 0.61 of the phenotypic variances for traits YIELD, TILLER, GRAIN, and KGW, respectively. These proportions are different from the broad sense heritability in the polygenic analysis (see Table 2). File S9 stores only the significant effects for the four traits (a subset of File S8). File S10 gives the independent variables converted from the genotypes for all the significant effects. For example, trait KGW has 15 significant effects, and thus the data set (sheet KWG in File S10) is stored in a matrix plus the bin and effect type information.

Detailed examination of Table 5 shows that the relative contributions of different types of effects vary across traits. Epistatic effects play more important roles than main effects for traits YIELD, TILLER, and GRAIN while the additive effect is more important for trait KGW. This general conclusion is consistent with that of the polygenic analysis.

## Discussion

Epistasis is often considered important to the variation of complex traits (Cockerham 1954). However, dissection of epistasis is difficult prior to the genome era because epistatic variances are often confounded with the additive and residual variances. To separate epistatic variances from additive and dominance variances, one needs pedigree information and the pedigrees must be very complicated to make the IBD matrices of epistasis sufficiently different from the kinship matrix for additive effects. In human populations, pedigree information is often incomplete and shallow (traced back to just one or two generations) and thus does not give us sufficient power to dissect epistasis. In the genome era, genome-wide marker information is available, which gives us an ample opportunity to model epistatic effects. Pedigree information is no longer crucial and a simple line cross experiment may be sufficient. In this study, we used binned genotypes of an IMF2 population to derive all the four types of kinship matrices to facilitate estimation of all types of epistatic variance components.

The result of the polygenic variance component analysis showed that adding epistatic variances can increase the model goodness of fit. We treated the additive model as the basis and progressively added more variance components to the model in the following sequence: *a*, *a* + *d*, *a* + *d* + *aa*, *a* + *d* + *aa* + *dd*, *a* + *d* + *aa* + *dd* + *ad*, and *a* + *d* + *aa* + *dd* + *ad* + *da*. Table 6 shows the increase of goodness of fit expressed as the squared correlation coefficient between the observed and predicted phenotypic values for the 278 IMF2 crosses. Depending on the characteristics of different traits, the goodness of fit has increased from 51.5 to 99.8% for YIELD and from 89.8 to 99.6% for KGW. The high goodness of fit implies an overfitting issue for the full model. Increasing sample size will partially correct the overfitting problem. Interestingly, the additive model alone already fits well for all the traits, especially for trait KGW. The reason may be explained for the high correlation between different estimated variance components, which may be caused by high similarities of various types of kinship matrices. The estimated variance components under different model sizes are listed in Table S1. Under the full model (all six genetic variance components are fit), the estimated dominance variances are zero for all the four traits. However, when the model size decreases, this estimated variance component can be nonzero, indicating that the actual dominance variance may not be zero but it is captured by other variance components due to the high correlation. The high correlation of estimated genetic variance components and high similarity between different kinship matrices suggest that benefit from fitting multiple polygenic covariance structures may not be significant for small sample sizes. Large sample sizes are necessary to separate the different genetic variances. How large of a sample size is sufficiently large? Theoretical investigation along with more simulation studies are required to have a clear answer for this question. The idea of incorporating multiple epistatic variance components in GWAS may be extended to incorporating multiple chromosome specific variance components in GWAS. How much improvement can we get by fitting multiple chromosome variances? Again, there is no clear answer at this moment. We may need a very large sample size to show the benefit.

A zero-estimated variance component for a particular type of effect does not mean that there is no genetic effect of this type for a particular bin or bin pair. The polygenic variance component of each type is the average of all bin- or bin pair-specific variances. A small number of bins or bin pairs may have large effects, but their effects are “diluted” under the polygenic model, which causes a zero estimation of the variance component. Our data analysis shows a lot of significant dominance effects but the estimated dominance variances obtained from the full model are zero for all traits.

In the original study of the IMF2 population, Zhou *et al.* (2012) used a different model to identify epistatic effects. They fit a two-locus model using a 3 × 3 factorial design (two loci and three genotypes per locus) and used a two-dimensional genome scan. They did not fit the polygenic variances. How much improvement have we achieved by fitting the polygenic variances? We do not have the answer at the moment because that would involve reanalysis of all the entire genome using the unweighted mixed model. Analysis of one trait, including permutation tests, took about 10 days of computing time because the large number of model effects fit to the model under the epistatic effect model, which is , where the 4 in the expression represents the four types of epistatic effects. Advantages of fitting the polygenic variances (weighted mixed model) have been demonstrated under the main effect (*a* + *d*) model using both simulated data as well as the rice data. Our results are also not comparable with the results of Zhou *et al.* (2012) because we pooled data of the 2 years together and ignored the potential G × E interactions.

Our study is different from the conventional GWAS (Yu *et al.* 2006) because we fit epistatic effects. However, the statistical model in our study shares similar properties with the improved eigenvalue decomposition methods (Lippert *et al.* 2011; Zhou and Stephens 2012) due to the fit of polygenic variances. We fit six variance components while the improved methods fit only one. Another difference is that we treated the bin effects and bin × bin interactions as random effects while Lippert *et al.* (2011) and Zhou and Stephens (2012) treated them as fixed effects. The advantages of fitting the random effects come from the shared variances and the tests of variance components. For instance, all the four epistatic effects share the same epistatic variance and are tested as one unit using a single LRT statistics. This approach may share some of the natures in rare variant detection (Wu *et al.* 2011), in which a group of rare variants are detected together using a single score test. Compared with the additive effect, an epistatic effect is a rare variant because the coefficient, *e.g.*, , has a smaller variance across individuals than the variance of or . Testing all four epistatic effects together will increase the statistical power for the same reason as the rare variant score test.

Finally, existing SAS programs were used here and the SAS codes are extremely simple (see File S11 and File S12). Theoretically, for each bin or bin pair, the six genetic variance components and the bin or bin pair specific variance can be estimated simultaneously using the HPMIXED procedure in SAS. This would be called the exact method (Zhou and Stephens 2012), as opposed to the approximate method (Lippert *et al.* 2011), where the polygenic variance is estimated only once before the genome is scanned. Our method is considered the approximate method because we estimated the polygenic variances prior to the two-dimensional scan. We did not adopt the exact method for two reasons: (1) there is a computational concern because we would have to estimate seven genetic variances per bin pair (high computational cost), and (2) there is no empirical evidence that the exact method outperforms the approximate method, although the former is theoretically superior over the latter. Using simulated data for the main effect model, we did not see much difference between the exact and approximate methods (data not shown). Zhou and Stephens (2012) proposed the exact method not because the exact method is better than the approximate method, but because it will avoid any “potential risk” of the approximate method and the exact method does not increase the computing time due to a special algorithm (eigen decomposition) for likelihood function evaluation. Once the polygenic variances are treated as constants in the approximate method, we provided a weight variable (see Equation 22 for the weight) to the HPMIXED procedure without modifying the code. In addition, the input variables were transformed using the eigenvector matrix prior to calling the SAS procedure. For readers without easy access to the SAS software packages, an R code is available from the author on request.

## Acknowledgments

The author is grateful to two anonymous reviewers for their comments on early versions of the manuscript. The author also appreciates Dr. Qifa Zhang for sharing some additional data for the IMF2 population of rice beyond the data posted on the journal website. The project was supported by the U.S. Department of Agriculture National Institute of Food and Agriculture Grant 2007-02784.

## Footnotes

*Communicating editor: N. Yi*

- Received July 7, 2013.
- Accepted September 16, 2013.

- Copyright © 2013 by the Genetics Society of America