Mapping Quantitative Trait Loci by Controlling Polygenic Background Effects

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 $n\times m=278\times 1619$ 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

where ${\overline{y}}_1$ and ${\overline{y}}_2$ 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,

where $Z_{jk}$ and $W_{jk}$ represent the codes for additive and dominance effects, respectively. Let y be an $n\times 1=278\times 1$ column vector for the phenotypic values of all the IMF2 crosses. The full epistatic model for $m=1619$ bins is described by

where $X\beta$ represents some nongenetic effects, e.g., intercept in this case, and $\epsilon \sim N(0,I{\sigma }^2)$ is a vector of residual errors with a normal distribution of unknown variance ${\sigma }^2$. In addition, $Z_k\#W_k$ represents element-wise vector multiplication, and $a_k$ and $d_k$ are the additive and dominance effects, respectively, for bin k. The four terms, ${(aa)}_{kk’}$, ${(dd)}_{k{k}^{\prime }}$, ${(ad)}_{k{k}^{\prime }}$ and ${(da)}_{k{k}^{\prime }}$ are the additive × additive, dominance × dominance, additive × dominance and dominance × additive effects, respectively, between bins k and $k^{\prime }$. These four terms are called the epistatic effects. The total number of genetic effects in the above model is

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 $a_k\sim N(0,{\sigma }_a^2)$, $d_k\sim N(0,{\sigma }_d^2)$,${(aa)}_{k{k}^{\prime }}\sim N(0,{\sigma }_{aa}^2)$, ${(dd)}_{k{k}^{\prime }}\sim N(0,{\sigma }_{dd}^2)$, ${(ad)}_{k{k}^{\prime }}\sim N(0,{\sigma }_{ad}^2),$ and ${(da)}_{k{k}^{\prime }}\sim N(0,{\sigma }_{da}^2)$ 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 $\text{E}(y)=X\beta$ and the variance–covariance matrix of y is

where the K’s are marker-generated kinship matrices; i.e., they depend on $Z_k$ and $W_k$. 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

where ${\lambda }_a={\sigma }_a^2/{\sigma }^2$ 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

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

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

for loci k and $k^{\prime }$, 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

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

A complete scan for this model (11) requires $m(m-1)/2$ calls of the above model and four epistatic effects are estimated under each call. The mixed model can be written generically as

where $e=\xi +\epsilon$ is a general error term with $E(e)=0$ and

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,

and the weights are the ratios of the variance components to the residual variance. For the additive and dominance effect model, $Z=Z_k\left|\right|W_k$ and $\gamma ={[\begin{array}{cc}{a}_k& d_k\end{array}]}^{\text{T}}$. For the epistatic effect model,

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 $N(0,I{\sigma }_{\gamma }^2)$ distribution. The shared common variance ${\sigma }_{\gamma }^2$ is estimated from the data. When ${\sigma }_{\gamma }^2$ is set to ∞, the genetic effects γ become fixed effects.

A two-step approach to genome scanning

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 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 $217(217-1)/2=23437$ 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 (${\lambda }_X={\sigma }_X^2/{\sigma }^2$) 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 ($h_X^2={\sigma }_X^2/{\sigma }_{\text{P}}^2$). 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 ($H=0.9636$). 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.17E–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 $\mathrm{var}(y)$ was calculated from the data and ${\sigma }_{\text{P}}^2$ 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

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 $23\times 2+66\times 4=310$ 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 $25\times 2+1071\times 4=4334$ 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 $15\times 278$ 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.