Abstract
In this article, I propose a mixedmodel method to detect QTL with significant mean effect across environments and to characterize the stability of effects across multiple environments. I demonstrate the method using the barley dataset by the North American Barley Genome Mapping Project. The analysis raises the need for mixed modeling in two different ways. First, it is reasonable to regard environments as a random sample from a population of target environments. Thus, environmental main effects and QTLbyenvironment interaction effects are regarded as random. Second, I expect a genetic correlation among pairs of environments caused by undetected QTL. I show how random QTLbyenvironment effects as well as genetic correlations are straightforwardly handled in a mixedmodel framework. The main advantage of this method is the ability to assess the stability of QTL effects. Moreover, the method allows valid statistical inferences regarding average QTL effects.
THERE are several different strategies to map quantitative trait loci (QTL; Kearsey and Farquhar 1998), e.g., singlemarker locus analysis (Liu 1998); simple interval mapping (IM; Lander and Botstein 1989); composite interval mapping (CIM; Zeng 1993, 1994), also called multiple QTL mapping (MQM; Jansen and Stam 1994); simplified CIM (Tinker and Mather 1995); marker regression (Kearsey and Hyne 1994; Wu and Li 1994); Bayesian methods (Sillanpää and Arjas 1998); and multiple interval mapping (MIM; Kaoet al. 1999; Zenget al. 1999). The latter methods have been shown to yield better power of QTL detection than IM and singlemarker locus analysis (Liu 1998; Lynch and Walsh 1998). In this article the main focus is on CIM for its simplicity and its importance in practice.
Within a frequentist framework, two different modelfitting approaches can be distinguished among procedures for QTL mapping by CIM: those based on maximumlikelihood (ML) estimation (Lander and Botstein 1989; Zeng 1994) and those using multiple linear regression (Haley and Knott 1992; Martinez and Curnow 1992, 1994). The latter approach usually gives results very similar to the former, though it must be seen as an approximate method. The main advantage of the multiple regression approach is its computational simplicity (Martinez and Curnow 1994). The computational tradeoff may become quite considerable in complex situations, i.e., when the model is extended to cover several random effects, as in this article. I therefore prefer the regression approach to ML and also to Bayesian methods using cofactors (Sillanpää and Arjas 1998).
In this article, I analyze the barley dataset by the North American Barley Genome Mapping Project (Han and Ullrich 1993). This dataset comprises 16 trials performed in different environments. The aim of the analysis is to detect QTL with significant mean effect across environments and to characterize the stability of effects across multiple environments. A need for mixed modeling arises in two ways. First, for an assessment of mean effects and of stability, it is necessary to regard environments as a random factor. Second, I expect a genetic correlation among pairs of environments caused by undetected QTL (Korolet al. 1998). Random QTLbyenvironment effects as well as genetic correlations are straightforwardly handled in a mixedmodel framework. The main advantage of my method is the ability to assess the stability of QTL effects. Moreover, the method allows valid statistical inferences regarding average QTL effects.
MATERIALS AND METHODS
The data: I use the sixrow Steptoe/Morex mapping population by the North American Barley Genome Mapping Project (Han and Ullrich 1993). The data comprise 150 doubled haploid (DH) lines, 223 markers, and 16 environments in the United States and Canada in 1991 and 1992. The same 150 genotypes were tested in all environments. The data are partially replicated with two replications. I use genotypebyenvironment means for different analyses.
The model for a single environment: The model will be adjusted to DH and backcross progeny data, but it is readily extended for F_{2} and other populations. Assume that the F_{1} cross is M_{1}QM_{2}/m_{1}qm_{2} with respect to the two flanking markers bordering the interval of interest (M_{1} and M_{2}) and the QTL (Q). Define a random variable g_{i} from the ith genotype taking value g_{i} = 1 if the F_{1} gamete carries the Q allele and g_{i} = 0 if the F_{1} gamete carries the q allele. In the regression approach to QTL mapping (Haley and Knott 1992; Martinez and Curnow 1992), the phenotype is regressed on the expected value of g_{i}, given the flanking markers. Table 1 gives explicit expressions for z_{i} = E(g_{i} flanking markers) assuming no interference (Martinez and Curnow 1992). Missing marker data can be handled using the method of Martinez and Curnow (1994). In the case of backcross progeny, the genetic effect for z_{i} constitutes a mixture of additive and dominance effects, whereas for doubled haploids, the genetic effects are confined to additive effects. For CIM using DH data, the following basic regression model can be used (Lynch and Walsh 1998, p. 465),
Model for data from multiple environments: The basic model (1) may be taken as a building block for more refined modeling to account for the design. Here, I am particularly interested in modeling data from multipleenvironment trials (METs). Specifically, I contend that a realistic model for genotypebyenvironment effects is needed to make unbiased inferences regarding QTL effects and positions. Moreover, the problem of genetic correlation among environments needs to be taken into account in case the same set of genotypes is tested in the different environments (Korolet al. 1998). Let y_{ij} be the observation of the ith genotype (i = 1,..., N) in the jth environment (j = 1,..., M). The proposed model is
All environmentspecific deviations (u_{j}, a_{j}, g_{kj}, and d_{ij}) and the residual genetic effect e_{i} are regarded as random normal deviates. To fully state the model, I need to specify the variances and covariances of all random terms. The variancecovariance structure should allow sufficient generality to realistically model real data. Before stating the full second moment assumptions, model (3) is modified to allow more generality. Model (3) contains a residual genetic effect e_{i} and a residual d_{ij}. This model corresponds to the usual factorial partitioning of main effects and interaction for genotypebyenvironment data in plant breeding and quantitative genetics (Lynch and Walsh 1998). To highlight the fact that d_{ij} contains both error and genotypebyenvironment interaction, I may write d_{ij} = f_{ij} + h_{ij}, where f_{ij} is the interaction and h_{ij} is error. The customary assumption in mixedmodel analyses, henceforth denoted as compound symmetry assumption, is that e_{i}, f_{ij}, and h_{ij} are identically distributed with zero mean and constant variances,
Genotypebyenvironment effects: I regard the genetic components in e_{ij} as random. When testing the same set of genotypes in the various environments, as has been assumed so far, genetic correlation among observations on the same genotype made in different environments needs to be allowed for. Many articles on the mapping of QTL based on multienvironment data corresponding to this design (Jansenet al. 1995; Beavis and Keim 1996; Romagosaet al. 1996; SariGorlaet al. 1997; Korolet al. 1998) have employed a multiple regression model with independent errors. This model implicitly assumes absence of genetic correlation, which is an unrealistic assumption. It is to be expected that in the presence of positive genetic correlation, the information on QTL position provided by data from a sample of environments is smaller than when genetic correlation is absent. This is best explained by a simple example (not specifically designed for QTL analysis). Consider two observations y_{1} and y_{2}. The mean (y_{1} + y_{2})/2 has variance
QTLbyenvironment effects and environmental main effect: In my analysis inferences are to be drawn with respect to a target population of environments. I regard the testing environments as a random sample from the target. The purpose of a mixedmodel analysis is to reveal mean QTL effects across environments. Here, random QTLbyenvironment interaction essentially plays the role of an error term. Moreover, the stability of QTL effects across environments is an important aspect. The larger the variance of QTLbyenvironment effects the lower the stability. Finally, I can make environmentspecific inferences, employing best linear unbiased predictions (BLUPs) of QTLbyenvironment effects.
It is necessary to allow for correlation among the effects in b_{j}. For example, a perfect correlation must be assumed for regression coefficients pertaining to adjacent markers, since both are linear in the additive genetic effect of the flanked QTL (Whittakeret al. 1996). Also, variances of two components in (b_{j}) corresponding to a pair of markers will have to be heterogeneous, considering the explicit expressions given in Whittaker et al. (1996, p. 25, right, bottom). Moreover, different QTL may be responding similarly to differential environments, giving rise to positive correlation among regression coefficients corresponding to markers adjacent to different QTL. For these reasons I make the general assumption that
The effect of the putative QTL in the jth environment is given by (α + a_{j}). The diagonal element in G corresponding to a_{j} can therefore be interpreted as a measure of stability of the effect of the putative QTL. The larger the variance of a_{j}, the more variable/less stable are the environmentspecific QTL effects (α + a_{j}). The breeder will seek a large absolute value for α and a small variance for a_{j}, i.e., a high stability. This interpretation of a variance component as a measure of QTL effect stability is akin to approaches for assessing yield stability in MET data (Linet al. 1986; Becker and Léon 1988; Piepho 1998b).
It is convenient to write the full model for y_{ij} in matrix form as
Why analysis of means across environments is problematic: If the same set of genotypes is tested in different environments, it is tempting to compute genotype means across environments and subject these to standard CIM. Such an analysis assumes that means of different genotypes,
Model selection: In what follows I assume random environments and consider CIM for mean QTL effects across environments. Model selection is necessary regarding three aspects: (1) the markers to be used as cofactors, (2) the variancecovariance structure for R, and (3) the model for QTLbyenvironment interaction, i.e., for G. These selection problems are briefly discussed. Unfortunately, the three problems cannot be tackled in an entirely independent manner. For example, in a joint analysis of the data, the chosen variancecovariance model will have an effect on the selected set of markers and vice versa. From a theoretical point of view it seems desirable to handle these three model components simultaneously. Due to the large number of candidate models (i.e., combinations of choices for 13), however, this simultaneous approach is not usually feasible in practice. Thus, some form of sequential approach is preferable. While this may entail the risk of missing some goodfitting models, it has the important advantage of reducing the total number of models to be considered. I suggest to first select the markers, then the variancecovariance structure for R, and finally the model for QTLbyenvironment interaction (G). At each step, I use the Schwarz Bayesian Criterion (SBC) to choose among options (Wolfinger 1993; McQuarrie and Tsai 1998). The criterion is given by SBC = log L  ^{1}/_{2}p log(n), where L is the maximized likelihood, p is the number of parameters, and n is the number of observations. Models with large values for SBC are preferred.
Cofactor selection: I select cofactors by multiple linear regression of means y on marker types using the marker pair selection (MPS) approach by H.P. Piepho and H. G. Gauch (unpublished results). This procedure has three distinctive features: (i) markers are selected in adjacent pairs to increase the chance of selecting flanking markers while reducing the risk of selecting markers not linked to QTL; (ii) an exhaustive search per chromosome is used in place of simple forward selection, which reduces the risk of missing the best fitting model; (iii) a model selection criterion such as SBC is employed to select the final model among a sequence of models. Among a selected pair, I use the marker that fits best as a cofactor for CIM. The procedure was developed for models with a single error term. Extension to the mixed model (4) for the replicate data is straightforward in principle but not generally feasible at present, mainly because of the prohibitive workload of having to fit a multitude of complex models by ML or restricted maximum likelihood (REML). This is the reason for my suggestion to apply MPS to means y. The approach is of an ad hoc nature, considering the fact that strictly speaking the means violate the independence assumption. As more efficient software becomes available, applying MPS to the replicate data using a mixedmodel framework is preferable.
Genotypebyenvironment interaction (R): I initially assume the most general model for QTLbyenvironment interaction on the basis of the genetic model corresponding to the selected cofactors (Wolfinger 1993). At this stage, the model does not yet contain the covariate z_{i} for the pair of markers flanking the putative QTL. Interactions are regarded as fixed, although at later steps of the analysis I take QTLbyenvironment interaction as random. No specific model is assumed for the interactions. The model I select for R is used subsequently when modeling QTLbyenvironment interaction and when scanning the genome for putative QTL.
QTLbyenvironment interaction and environmental main effect (G): I now take QTLbyenvironment interactions and the environmental main effect as random. Thus, I have the task of selecting an appropriate variancecovariance structure for G. The need for invariance under recoding of the markers dictates the unstructured model for G, while the parsimony principle suggests that simpler approximating structures may be worth considering. I propose to generally fit an unstructured model, except when the dimension of G is large and an unstructured model is difficult to fit.
Scanning the genome: The same model as selected for both R and G in the previous steps is used when scanning the genome for QTL. Of course, R and G are reestimated at each putative QTL position. Note that G is extended at the scanning stage by the covariate z_{i} for the two flanking markers. Assuming the same structure as selected for the case where x_{i} contains only the cofactors may not be optimal. It would not be practicable, however, to select a different model at each step during the genome scan. Also, the type of model finally chosen for R and G not only depends on model fit but also on ease of estimation. Some structures for R and G may be well fitting but difficult to estimate (convergence problems, difficulty in choosing good starting values, etc.), thus making them infeasible for automated QTL scans, where the same model has to be estimated a large number of times.
EXAMPLE
I used the barley data to exemplify the proposed methods. MPS picked the marker pair M81/M82. Between these two, M82 had a better fit than M81 if fitted alone. Subsequently, the cofactor M82 and the interactions with environments were included in the fixed part of the model and various structures were fitted to R by the REML method. At this stage, x_{i} did not yet contain a covariate z_{i} for the putative QTL. All mixedmodel analyses were done using ASREML (Gilmouret al. 1999). The same genotypes have been tested in all environments, so genetic correlation needs to be modeled in R. Since there are M = 16 environments and hence R is a M × M matrix, there are M(M + 1)/2 = 136 parameters for the unstructured model. The results for different models are shown in Table 2. On the basis of SBC the factoranalytic model
The Wald Fstatistic for cofactor(M82)byenvironment interaction in b was 10.84, which is significant at α = 5%. Thus, the interaction term was retained. In the next modelbuilding step, the environmental main effect and interaction with M82 was regarded as random and different structures were fitted to G, using the model
I computed parameter estimates and standard errors of the fixed effects based on the selected models for R and G. For comparison, I also computed parameter estimates and standard errors on the basis of a multiple linear regression with means y. The results are reported in Table 4. The parameter estimates do not differ much in the two analyses, but the standard errors are larger with the mixedmodel analysis. The analysis shows that it is not appropriate to ignore correlations among genotypes as is done in a simple analysis of means.
The next step is to scan the genome for a QTL by CIM. For this purpose, x_{i} needs to be augmented by the covariate z_{i} for the putative QTL. This again raises the question of how to model G. A priori, the unstructured model seems most appropriate, especially since there is only one cofactor so that parsimony is not a pressing issue. Thus, I used the unstructured model. The window size was 10 cM; i.e., for putative QTL within 10 cM of a cofactor, the cofactor was dropped from the model. The step size of the chromosome scan was 1 cM. During the chromosome scan, parameter estimates for G and R from the present putative QTL position were used as starting values at the next position. This resulted in convergence of the REML algorithm within a few iterations (typically 510). [One referee noted that this choice of starting values may cause convergence to a local, but not the global maximum of the (restricted) loglikelihood. In my experience the likelihood of this problem is small with a short step size. In the present example the problem was not observed. The same referee indicated that using fixed starting values gives more stable results.] At each position, I computed a Wald statistic (F) for the test of the null hypothesis of no QTL. Conditionally on the position, F asymptotically follows a χ^{2} distribution with 1 d.f. On chromosome segments with the same set of cofactors for CIM, the type I error rate was controlled using the quick method proposed by Davies (1977, 1987). An overall threshold controlling the genomewise type I error rate was obtained by the Bonferroni method. The critical value was 12.83. The F profile for the seven barley chromosomes is shown in Figure 1. There is a significant QTL on chromosome 3 near the 58.9cM position, where the Fstatistic reaches a maximum of 20.62. Chromosome 5 shows some indication of a second QTL, though the F values are not significant. For comparison, I also computed the F profile on the basis of a regression analysis using means y. The critical value was 13.01. The profiles are similar in shape, though there are some notable differences to the mixedmodel analysis. The peak on chromosome 3 is much more pronounced, reaching a maximum of ∼80, while for the mixedmodel analysis the maximum value was ∼20. At the estimated QTL position on chromosome 3 (58.9 cM), the estimated QTL effect was
DISCUSSION
A common feature of QTL analyses is that QTL effects depend on environment. Many researchers have dealt with this problem by analyzing each environment separately. This approach is quite useful, if one is interested in the particular test environments. As pointed out by Tinker and Mather (1995), “separate analysis by environment circumvents the problem of dealing with QTLbyenvironment interaction and avoids complications due to environmental heterogeneity. However, the results of separate analyses are difficult to interpret, and they do not take advantage of the builtin replication provided by multiple environments.” Quite frequently, test environments are just a sample from a target population, and the breeder is interested in making broad inferences not restricted to the particular test environments (Melchingeret al. 1998). This objective calls for a mixedmodel analysis with random environments (Beavis and Keim 1996). The present article has shown how to implement such an analysis for CIM. The approach is easily extended to cover multipleQTL and epistatic effects by appropriate modifications of β and b_{j} (MorenoGonzalez 1992; Wanget al. 1999). This is useful as an additional analysis step, when several QTL have been detected by CIM.
A simple alternative to a mixedmodel analysis of MET QTL data is to proceed in two steps as follows: first, the quantitative trait is analyzed by ANOVA techniques to obtain (adjusted) genotype means across environments. Second, the means together with the marker data are submitted to a routine for QTL analysis. My theoretical considerations and analysis of a real dataset led me to conclude that such analyses may lead to inappropriate inferences, mainly because standard error estimates are inappropriate. A mixedmodel framework allows more valid inferences to be obtained by incorporation of different random components of variance that appropriately account for the environmental and genetic structure of the data.
Wang et al. (1999) use a mixedmodel approach similar to the one presented here. There are some differences, the main one being that Wang et al. (1999) do not allow for genetic correlation. They exploit map information to model the equivalent of my G. This approach requires a specific genetic model, including epistasis and specification of multiple QTL. It assumes that all genetic effects are modeled by G and that all correlation among effects is solely due to the map. It is to be expected that the approach is susceptible to misspecification of the model. By contrast the model used in this article allows unexplained genetic effects and genetic correlation to be subsumed in the residual e_{ij}. As a result there is no need to assume a specific overall model. Moreover, using a general model for G allows for covariance among QTL that is due to correlated response to differential environmental condition. Such correlation can arise even for QTL on different chromosomes, which would be regarded as independent under the model used by Wang et al. (1999). My model could be easily extended to exploit the map in modeling the correlation among cofactors and the putative QTL in G, but I prefer to work with a general model for G both for ease of computation and to reduce the danger of model misspecification. Also, if the approach is extended to multiple QTL and epistatic effects, it is still preferable to work with general G for the same reasons. When the dimension of G becomes large, it is reasonable to consider more parsimonious models such as those used for R.
In this article I have demonstrated how to use mixed models for assessing mean and stability of QTL effects based on genotypebyenvironment means from MET data. My mixedmodel framework is easily extended to other settings, e.g., when spatial heterogeneity needs to be modeled at the plot level (nearest neighbor analyses; Moreauet al. 1999) and when the experimental designs give rise to random effects (incomplete blocks, etc.; Haley and Knott 1992). Error variance heterogeneity among different environments in MET can be accounted for by weighting (Culliset al. 1996), though in my experience weighting has little effect on final parameter estimates and standard errors.
My analysis has assumed random environments. In a model with fixed environments, the focus is on studying QTLbyenvironment interactions. It has been stressed by Korol et al. (1998) that the number of interaction parameters in models such as (4) increases linearly with the number of environments. These authors used the regression approach of Eberhart and Russell (1966) to model interactions with fewer parameters. The regression approach by Eberhart and Russell (1966) was originally proposed for the analysis of genotypebyenvironment interaction, but it can be applied in the same way to model QTLbyenvironment interaction, as demonstrated by Korol et al. (1998). In fact, there are a large number of different models for genotypeenvironment interaction (Linet al. 1986; Becker and Léon 1988; van Eeuwijket al. 1996; Piepho 1998b), all of which are potentially useful for modeling QTLbyenvironment interaction. This potential seems to have gone largely unnoticed in QTL work. For example, extensions of the EberhartRussell regression such as the additive main effects multiplicative interaction model (Gauch 1988) typically explain a much larger fraction of the total interaction and so promise improved performance (Romagosaet al. 1996). These approaches can be incorporated into my mixed model by taking b as fixed and imposing some structure such as in EberhartRussell regression.
Acknowledgments
Thanks are due to Hugh Gauch Jr. and Susan McCouch for inspiring discussions. H.F. Utz (University of Hohenheim, Germany) is thanked for helpful comments on an earlier draft. Support of the Heisenberg Programm of the Deutsche Forschungsgemeinschaft is gratefully acknowledged. Part of the research for this article was conducted while the author was visiting the Department of Biometrics and the Department of Plant Breeding, College of Agriculture and Life Sciences, Cornell University, Ithaca, NY.
Footnotes

Communicating editor: J. B. Walsh
 Received March 22, 2000.
 Accepted July 25, 2000.
 Copyright © 2000 by the Genetics Society of America