## Abstract

The effects of quantitative trait loci (QTL) on phenotypic development may depend on the environment (QTL × environment interaction), other QTL (genetic epistasis), or both. In this article, we present a new statistical model for characterizing specific QTL that display environment-dependent genetic expressions and genotype × environment interactions for developmental trajectories. Our model was derived within the maximum-likelihood-based mixture model framework, incorporated by biologically meaningful growth equations and environment-dependent genetic effects of QTL, and implemented with the EM algorithm. With this model, we can characterize the dynamic patterns of genetic effects of QTL governing growth curves and estimate the global effect of the underlying QTL during the course of growth and development. In a real example with rice, our model has successfully detected several QTL that produce differences in their genetic expression between two contrasting environments. These detected QTL cause significant genotype × environment interactions for some fundamental aspects of growth trajectories. The model provides the basis for deciphering the genetic architecture of trait expression adjusted to different biotic and abiotic environments and genetic relationships for growth rates and the timing of life-history events for any organism.

A number of statistical methods have been developed to map quantitative trait loci (QTL) affecting complex traits in well-structured pedigrees (Lander and Botstein 1989; Jansen 2000) or natural populations (Wu* et al.* 2002a). These methods have been instrumental in the identification of QTL responsible for various quantitative traits of agricultural, biological, or biomedical value (Wu* et al.* 2000; Mackay 2001; Zimdahl* et al.* 2002). However, the derivations of these statistical methods that were based on a simplified assumption that there is a direct relationship between the genotype and phenotype did not consider sequential developmental pathways that form the phenotype. It is likely, therefore, that these methods lead to biased estimation of QTL positions and effects and have limited power to detect QTL involved.

More recently, a host of new statistical models have been proposed to map a complex trait by taking into account its underlying developmental changes (Ma* et al.* 2002; Wu* et al.* 2002b, 2004a,b). These models integrate universal growth laws derived from fundamental biological principles and described by mathematical models into a QTL mapping framework through statistics as a bridge. They do not directly estimate the expected means of different QTL genotypes at different time points, but rather fit these means using growth curves as a function of time. Thus, the estimation of QTL effects on growth is equivalent to the estimation of model parameters describing the shape of growth curves. These models, called *functional mapping*, have successfully detected dynamic QTL that affect developmental changes during ontogenetic growth in a long-lived forest tree (Wu* et al.* 2003) and the animal model system—mouse (Zhao* et al.* 2004).

Current functional mapping models that have examined the effects of individual QTL on growth trajectories have focused on a particular phenotype measured at a finite set of time points in a single environment (Ma* et al.* 2002; Wu* et al.* 2002b). Other studies have examined genotype × environment interactions involving multiple environments (Piepho 2000), but have focused on a single measurement of a phenotype and hence cannot provide a comprehensive picture of how a QTL affects the dynamic process of growth under a range of environmental conditions. Wu* et al.* (2004a) proposed a model for characterizing the effects of epistatic QTL on growth trajectories but did not incorporate environmental influences into the model.

A great wealth of evidence suggests that different expressions of the same genotype across environmental contexts, referred to as phenotypic plasticity (Wu 1998), are controlled by specific QTL (Mackay 2001) through two different but not exclusive genetic mechanisms, allelic sensitivity and gene regulatory control (Wu 1998; Vieira* et al.* 2000). Yan* et al.* (1998) identified several QTL on the rice genome that display significant genotype × environment interaction effects on plant height measured at multiple time points of development, although their results were derived from interval mapping of single phenotypes. In her seminal review, Mackay (2001) claimed the existence of environment-specific QTL that direct organismic development toward the best utilization of resources in heterogeneous environments.

In this article, we extend our previous functional mapping model to map QTL that display genotype × environment interactions for developmental trajectories. Although this extended model is similar to the previous epistatic model (Wu* et al.* 2004a) in modeling the mean vector for time-dependent genotypic values, these two models are methodologically different in modeling the (co)variance matrix for longitudinal measurements of growth. The epistatic model needs to model only the structure of one (co)variance matrix, whereas the genotype × environment interaction model needs to model multiple (co)variance matrices each corresponding to an environment as well as to model the covariances between different environments.

Our extended model for characterizing genotype × environment interactions is constructed within the maximum-likelihood-based mixture model framework and unifies the genetic effects of QTL expressed in different environments through biological principles and statistical models. We incorporate the EM algorithm to solve the mixture model and provide estimates of growth curves specified by different QTL genotypes and of the parameters for modeling the structure of (co)variance matrices. Our model allows for the tests of a number of biologically meaningful hypotheses regarding QTL × environment interactions on growth processes. We derive a general model for converting the effects of genotypic curves to additive, dominant, and/or epistatic variance components. A real example with rice (Yan* et al.* 1998) validates our model that detects five environment-specific QTL responsible for plant height growth. When synthesized into evolutionary developmental biology (*evo-devo*; Raff 2000; Arthur 2002), our model helps to address an old but still unsolved question of how genetic factors regulate developmental processes in an organism from embryo to adult.

## THE GROWTH LAW

According to von Bertalanffy (1957), growth in weight, length, or size will occur whenever the anabolic or metabolic rate exceeds the rate of catabolism. Thus, the ratio of these two processes indicates the occurrence and change of growth. When the ratio approaches and then drops below unity, growth will decrease and eventually cease. von Bertalanffy also noted that the metabolic rate of an animal [in fact, any organism, as observed by succeeding researchers (West* et al.* 2001)], scales as the *k*th power of its weight (Niklas 1994; West* et al.* 1997) but the catabolic rate is proportional to the weight itself. Therefore, the growth rate, *i.e.*, the difference between these two rates, becomes 1where *g*(*t*) represents the growth at time *t* and η and λ are the constants of metabolism and catabolism, respectively. For small values of *k*, integration of Equation 1 leads to the growth equation, 2where *g*_{0} is the growth at *t* = 0. This growth function is sigmoidal (S-shaped), approaching asymptotically the value (η/λ)^{1/(1−}^{k}^{)} as *t* → ∞. Such an S-shaped pattern of growth includes an exponential growth stage, an asymptotic growth stage, and the point of inflection at which these two stages are connected (West* et al.* 1999). At the point of inflection, an organism displays maximum growth per unit time. After the substitutions, where *a* is asymptotic (limiting) growth, *b* is a parameter related to the growth at *t* = 0, and *r* is the “rate constant” that determines the spread of the curve along the time axis. Richards (1959) rewrites Equation 2 as 3

Equation 3, referred to as the Richards growth model or law, has the inflection point whose coordinates are solved as [*t* = ln(*b/*(*k* − 1))/*r*, *g* = *ak*^{1/(1−}^{k}^{)}] (Richards 1959; Nath and Moore 1992; Gregorczyk 1998). By letting *k* take different values, the Richards growth model can be reduced to three well-known growth functions,

The monomolecular curve actually describes the asymptotic phase of the logistic growth curve. The logistic curve is symmetric about the inflection point at which *t* = (ln *b*)*/r* and *g* = *a*/2. Yet, the Gompertzian curve is not symmetric, with the inflection point occurring at *t* = (ln *b*)*/r* and *g* = *a/e*. These three curves, with different flexibilities, may best fit growth data collected from a different species or organ or in a different environment.

## THE STATISTICAL METHOD

### Genetic design:

There is a standard backcross design with *n* progeny, initiated with two contrasting homozygous inbred lines. Suppose a genetic linkage map covering the genome has been constructed with polymorphic markers. The marker density and positions have been known prior to QTL mapping. We consider a special case in which each of the backcross progeny is clonally replicated. It is possible to make clones for many plants such as rice, poplar, and pine.

Field trials with the backcross are conducted at *L* locations. At each location all *n* progeny are planted in a randomized complete block design with *B* clonal replicates. Each plant is measured for longitudinal growth traits, such as plant height and tiller number, at a finite set of time points (τ) during the course of plant growth in the field. Our model allows for unevenly spaced measurement time intervals and for discrepancies in measurement schedule among progenies.

Assume that there is a segregating QTL with alleles *Q* and *q* that affects growth curves or trajectories in the backcross population grown at different locations. There are two QTL genotypes in the backcross denoted by *j* (*j* = 1 for *Qq* and 0 for *qq*). A set of parameters describing the growth curve of genotype *j* at location *l* (*l* = 1, … , *L*) is denoted by 𝒢* _{lj}*. The comparisons of these parameters among the two different QTL genotypes across different locations can determine whether and how this putative QTL differently affects growth trajectories at different locations.

### The mixture model:

The growth laws described in the preceding section have been incorporated into a QTL mapping framework based on an finite mixture model composed of different components (Ma* et al.* 2002; Wu* et al.* 2002b). In such a mixture model, each observation *y* is assumed to have arisen from one of *J* groups of QTL genotypes, each group being modeled by a density from the parametric family *f* (Lander and Botstein 1989). The population density function of *y* is where ω̃ = (ω̃_{1}, … , ω̃* _{J}*) are the mixture proportions that are constrained to be nonnegative and sum to unity; φ = (φ

_{1}, … , φ

*) are the genotype-specific parameters, with φ*

_{J}*being specific to genotype group*

_{j}*j*(

*j*= 1, … ,

*J*); and ν is a parameter that is common to all genotype groups. The mixture proportions denoted as the frequencies of QTL genotypes depend on the marker genotypes of two flanking markers bracketing the QTL. The normal density functions associated with different QTL genotypes are expressed in terms of the expected value of each genotype and residual variance.

### The likelihood function:

The likelihood function of longitudinal data, **y**, measured for the backcross at different locations for the putative QTL is constructed as 4where the mixture proportions (ω̃* _{j}*'s) denote QTL genotype frequencies in the backcross,

*i.e.*, ω̃

_{1}= ω̃

_{0}=

^{1}/

_{2}, and the multivariate normal distribution of progeny

*i*measured at τ time points is expressed as 5where is a vector of longitudinal observation composed of

*B*×

*L*subvectors each for a location (

*l*) and a block (

*b*), with

**y**

*= [*

_{ibl}*y*(1), … ,

_{ibl}*y*(τ)], and component-specific parameters φ

_{ibl}*for genotype*

_{j}*j*are specified by for different locations and blocks. But since differences in growth curve are considered only at the location level, growth curves for different blocks within a location can be regarded as identical;

*i.e.*, we have with

**u**

*= [*

_{jl}*u*(1), … ,

_{jl}*u*(τ)] for genotype

_{jl}*j*at τ different points.

At a particular time *t* in block *b* and location *l*, the relationship between the observation and expected mean can be described by a linear regression model, where *x _{i}* is the indicator variable denoted as 1 if progeny

*i*carries QTL genotype

*j*and 0 otherwise;

*e*(

_{ibl}*t*) is the residual error that is iid normal with the mean of zero and the variance of σ

^{2}

_{bl}. The errors at two different time points,

*t*

_{1}and

*t*

_{2}, are correlated with the covariance of cov

*(*

_{bl}*t*

_{1},

*t*

_{2}). These (co)variances compose a (τ × τ) matrix

**∑**

*for block*

_{bl}*b*and location

*l*. A more general model is to let the errors among different blocks and locations be correlated for the same, cov

_{}, or different time points, cov

_{}. However, in our situations, it is reasonable to assume that there are no between-block and between-location error correlations because different clonal replicates of the same genotype are planted in different environments. With this assumption the (co)variance matrix

**∑**, whose elements are the common parameter ν of the mixture model (Equation 1), is simplified into a block-diagonal matrix with a

*B*×

*L*dimension and the likelihood of Equation 4 can then be rewritten as 6

The determination of the value for the indicator variable describing the genotypes of the QTL for progeny *i* is not obvious. According to the interval mapping theory, it is possible to do so if we use the segregation information of the known flanking markers that bracket the QTL. Suppose this QTL is bracketed by two flanking markers ℳ_{ℓ} (with alleles *M*_{ℓ} and *m*_{ℓ}) and ℳ_{ℓ} + 1 (with alleles *M*_{ℓ+1} and *m*_{ℓ+1}). Thus, the QTL genotype frequencies in the backcross population (denoted by ω̃* _{j}*) should be expressed as the conditional probabilities of the unknown QTL genotypes given the known marker genotypes. Table 1 tabulated these conditional probabilities, generally expressed as ω̃

_{j}_{|}

*, where*

_{i}*j*|

*i*stands for QTL genotype

*j*given a particular marker genotype for progeny

*i*. We rewrite the likelihood function of environment-specific longitudinal data (

**y**

*) and marker information genotyped for both sexes (ℳ) as 7*

_{ibl}Note that Equation 7 is different from Equation 6 as the latter does not make use of marker information whereas the former does.

### Modeling the mean vector and (co)variance matrix:

The estimation of the mean vector **u*** _{jl}* and the (co)variance matrix

**∑**

*is statistically difficult because they involve too many unknown parameters given a possible sample size. Also, such direct estimation does not take into account the biological principles of growth and development. We incorporate the universal growth law into the estimation process of the likelihood function (Equation 7). Thus, the mean value of QTL genotype*

_{bl}*j*in location

*l*at time

*t*is expressed by where the growth parameter set 𝒢

*= (*

_{jl}*a*,

_{jl}*b*,

_{jl}*r*,

_{jl}*k*) describes the asymptotic growth, initial growth, relative growth rate, and curve shape, respectively (Richards 1959). With this growth equation, we need to estimate only the growth parameters, rather than estimate genotypic means at every point, to detect genotypic differences in growth. This can significantly reduce the number of unknown parameters to be estimated, especially when the number of time points is large.

_{jl}Similarly, the (co)variance matrix can be structured with an appropriate model. Statistical analysis of longitudinal data has established a number of structural models that capture most of the information contained in the matrix (Diggle* et al.* 2002). Here, we use a first-order autoregressive [AR(1)] model to model the structure of the matrix, which is based on two assumptions, (1) the variance σ^{2} is constant over time, and (2) the correlation decays in a proportion of ρ purely with time interval. With the AR(1) model, we need to estimate only , instead of all elements in the matrix. The advantage of such a matrix-structuring model is to reduce the number of unknown parameters, without losing the information of the matrix. Many other structural models may be more advantageous over the stationary AR(1) model, but the choice of an optimal model in a particular situation should be based on statistical tests, as described in Zimmerman and Núñez-Antón (2001).

### Computational algorithms:

As classified above, the unknown parameters that build up the likelihood function (Equation 7) include the curve parameters, matrix-structuring parameters, and the QTL genotype frequencies specified by QTL position measured in terms of the recombination fractions (θ) between the QTL and its flanking markers (see Table 1). Arrayed by , these unknowns can be estimated through differentiating the log-likelihood function of Equation 7 with respect to each unknown, setting the derivative equal to zero, and solving the log-likelihood equations. This estimation process can be implemented with the EM algorithm as described below.

The log-likelihood function of growth and marker data for block *b* and location *l* based on Equation 7 is given by with the derivative with respect to any element Ω_{ξ}, where we define 8which is a posterior probability that progeny *i* carrying a particular marker genotype with measurements from block *b* and location *l* is regarded as carrying QTL genotype *j*. We then implement the EM algorithm with the expanded parameter set {**Ω**, **Π**}, where **Π** = {Π_{j}_{|}* _{ibl}*} (the E step; Equation 8). Conditional on

**Π**, we solve for 9to get the estimates of

**Ω**(the M step; Equation 9). The estimates are then used to update

**Π**, and the process is repeated between Equations 8 and 9 until convergence. The values at convergence are the maximum-likelihood estimates (MLEs) of

**Ω**. The iterative expressions of estimating

**Ω**from the previous step are given in Ma

*et al.*(2002) and Wu

*et al.*(2004b).

As usual, the QTL position parameter can be viewed as a known parameter because a putative QTL can be searched at every 1 or 2 cM on a map interval bracketed by two markers throughout the entire linkage map. The amount of support for a QTL at a particular map position is often displayed graphically through the use of likelihood maps or profiles, which plot the likelihood-ratio test statistic as a function of map position of the putative QTL (Lander and Botstein 1989).

## HYPOTHESIS TESTS

Different from traditional mapping approaches, our functional mapping for longitudinal traits allows for the tests of a number of biologically meaningful hypotheses (Wu* et al.* 2004a). These hypothesis tests can be a *global* test for the existence of significant QTL, a *local* test for the genetic effect on growth at a particular time point, a *regional* test for the overall effect of QTL on a particular period of growth process, or an *interaction* test for the change of QTL expression across times. These tests at different levels can be formulated to test the effects of QTL × location interaction on the shape of growth.

### Global test:

Testing whether specific QTL exist to affect growth trajectories is a first step toward the understanding of the genetic architecture of growth and development. The genetic control over entire growth processes can be tested by formulating the following hypotheses: 10

H_{0} states that there are no QTL affecting growth trajectories and the two genotypic curves at each location overlap (the reduced model), whereas H_{1} proposes that such QTL do exist (the full model). The test statistic for testing the hypotheses in Equation 10 is calculated as the log-likelihood ratio of the reduced to the full model, 11where **Ω̃** and **Ω̂** denote the MLEs of the unknown parameters under H_{0} and H_{1}, respectively. The log-likelihood ratio (LR) is asymptotically χ^{2}-distributed with 4(*L* − 1) d.f. An empirical approach for determining the critical threshold is based on permutation tests, as advocated by Churchill and Doerge (1994). By repeatedly shuffling the relationships between marker genotypes and phenotypes, a series of the maximum-log-likelihood ratios are calculated, from the distribution of which the critical threshold is determined.

After a significant QTL is detected, the next test is about the effect of this QTL on growth in each location. This will use the same form as shown in Equation 10, but focusing on a location. It is interesting to test whether the same QTL genotype is expressed identically across different locations to affect growth trajectories. Such a null hypothesis test can be formulated as 12which states that any two curves between the same QTL genotypes from different sexes overlap. By combining the test results from both QTL genotypes, this test can also be used to test for QTL × location interactions. However, when location-specific curves with the same QTL genotypes are approximately parallel to each other, the area under curve (*A _{jl}*) is an appropriate criterion for this QTL × location interaction test, expressed as

In this case, the null hypothesis for testing QTL × location interaction can be formulated as 13

*i.e*., the difference between the areas under curves for two alternative QTL genotypes is set equal for all the *L* locations.

In addition to testing overall genetic effects on growth trajectories, our model allows for the tests of the additive and dominant effect as well as their interaction effects with sexes. Wu* et al.* (2004a) proposed detailed procedures for making these specific tests, all of which can be directly used or modified for this study.

### Local test:

The local test can test for the significance of the genetic effect of QTL and QTL × location interaction effect on growth traits measured at a time point (*t**) of interest. For example, the hypothesis for testing the effect of QTL on growth at a given time *t** can be formulated as 14which is equivalent to testing the difference of the full model with no restriction and the reduced model with a restriction as set in the null hypothesis.

### Regional test:

Sometimes we are interested in testing the difference of growth trajectories in a time interval rather than simply at a time point. The question of how a QTL exerts its effects on a period of growth trajectories [*t*_{1}, *t*_{2}] can be tested using a regional test approach based on the areas, covered by load curves. The hypothesis test for the genetic effect on a period of growth process is equivalent to testing the difference between the full model with no restriction and the reduced model with a restriction.

### Tests for differentiation in environment-dependent genetic expression:

How are the same QTL genotypes expressed across different environments? The solution to this question helps select superior genotypes for a complex trait. We construct the following alternative hypotheses to test this: 15

The rejection of H_{0} implies that QTL genotypes are expressed throughout growth processes differently across different environments.

### Interaction test:

The effects of QTL may change with time, which suggests the occurrence of QTL × time interaction effects on growth trajectories. The differentiation of growth with respect to time *t* represents growth rate. If the growth rates at a particular time point *t** are different between the curves of different QTL genotypes, this means that significant QTL × time interaction occurs between this time point and the next.

### Test for biologically important parameters:

A number of biological parameters can be used to evaluate the developmental characteristics of growth. The logistic growth curve can be used to determine the coordinates of a biologically important point in the entire growth trajectory—*the inflection point*—where the exponential phase ends and the asymptotic phase begins (Niklas 1994). The time at the inflection point corresponds to the time point at which a maximum growth rate occurs. The difference in the coordinates between different QTL genotypes provides important information about the genetics and evolution of growth trajectories (Niklas 1994). The genotypic differences in time and growth at the inflection point of maximum growth rate can be tested. The test for their genotypic difference is based on the restriction for *t*_{Ijl}, and for *u*.

## RESULTS

Our newly developed functional mapping model is used to map dynamic QTL in a model plant—rice. A genetic linkage map was constructed with 135 RFLP and 40 isozyme and RAPD markers to cover 12 chromosomes (Figure 1) for a doubled-haploid (DH) population of 123 lines derived from semi-dwarf IR64 and tall Azucena (Huang* et al.* 1997). This map is 2005 cM long with an average distance of 11.5 cM between a pair of adjacent markers. The DH population with two genotypes at a QTL, *QQ* and *qq*, is analytically identical to a backcross population, so that the model described above can be directly used. The DH population was cloned and different clonal replicates of the same genotype were grown in a randomized complete design with two replicates at a spacing of 15 × 20 cm at two climatically contrasting locations, Hangzhou (subtropical, 30° 16′ N) and Hainan (tropical, 19° 57′ N), China. The experiments in both locations were carried out from late May to early November 1996. After 10 days of transplanting into the field trial, plant heights (from the surface of the soil to the tip of the plant) were measured every 10 days until all lines had headed (Yan* et al.* 1998).

This study used 90 DH plants whose marker and growth data are complete for both locations. Figure 2 illustrates growth curves of the means of plant heights between two blocks for these 90 DH plants separately for Hangzhou (Figure 2A) and Hainan (Figure 2B). On average, the Hangzhou plants grew slightly better than the Hainan plants, with two displaying different growth trajectories. Substantial variation in growth curve among different plants in each location suggests that specific QTL may be involved in shaping developmental trajectories.

As a general framework, we derived our functional mapping model on the basis of a more flexible Richards growth equation. However, this equation does not necessarily work better in practice than the reduced logistic equation (*k* = 2) because the former contains one more parameter than the latter. If both the equations provide a similar likelihood of the data, the parsimonious logistic equation should be used. In this example, the logistic equation is sufficient to fit our growth data and, therefore, is employed to search for growth QTL through a genome-wide scanning approach. The assumption of constant variance for the AR(1) model may not be true in our data, as indicated by increased variance with ages in both locations (Figure 2). Wu* et al.* (2004b) proposed a transformation approach, called the transform-both-sides (TBS) model by Carroll and Ruppert (1984), to reduce variance heteroscedasticity and, therefore, increase the power of the model. This TBS-based mapping model can also preserve the biological meanings of curve parameters. In this study, we incorporate the TBS-based model through log-transformation into the functional mapping framework for analyzing QTL × location interactions. As shown in Figure 2, C and D, the log-transformation can lead to relatively constant variances in plant height growth at both locations, although a more effective transformation approach should be estimated simultaneously with the model parameters (Box and Cox 1964; Carroll and Ruppert 1984).

By comparing with the genome-wide critical threshold determined by permutation tests, we detected five significant QTL (*P* < 0.001) on chromosomes 1, 3, 7, 9, and 11 that were responsible for growth trajectories of plant heights (Figure 3). The estimated locations of these QTL correspond to the genome positions showing the maximum-log-likelihood-ratio test statistics. Table 1 gives the estimates of the positions of the detected QTL each bracketed by two flanking markers on respective chromosomes, as well as the estimates of the curve and matrix-structuring parameters.

We drew the age-dependent expression profiles of two different QTL genotypes for each QTL identified from our model for plant heights grown in Hangzhou and Hainan (Figure 4). These profiles were drawn using the MLEs of growth curve parameters given in Table 2. The profiles of gene expression allow for the characterization of the developmental timing for a QTL to turn on or turn off or the period for the QTL to trigger its effect on growth. As indicated by two divergent curves each representing a different QTL genotype, *QQ* or *qq*, the detected QTL exhibited increasing effects on plant height in both locations as rice grows. The difference between the areas under the growth curves of different genotypes reflects the influence of the QTL on the overall growth process of plant height. We use Equation 13 to test whether a QTL interacts with locations to affect entire growth trajectories. The test result indicated that none of the five detected QTL displayed significant genotype × environment interactions, suggesting that they are more general growth regulators. Using the hypothesis test in Equation 15, we found that the time-specific growth QTL detected on chromosomes 1, 3, 9, and 11 are expressed consistently stronger in Hangzhou than in Hainan, whereas the QTL on chromosome 7 displayed an opposite expression pattern (Figure 4). As expected, at all the QTL, except for one on chromosome 7, the tall Azucena parent contributes height-increasing alleles in both locations, whereas the short IR64 parent contributes height-decreasing alleles.

Adult height for rice can be genetically controlled during an early stage of development by mapping the QTL determining the timing of maximum growth rate. Such early genetic manipulation for height growth can potentially preserve more energy for later reproductive propagation. The significant QTL detected for overall growth trajectories were all observed to affect the timing of maximum height growth rate (inflection point; Figure 4), all of which, except for the one on chromosome 7, displayed significant genotype × environment interaction effects on this timing (*P* < 0.001). The QTL on chromosomes 1, 3, 9, and 11 trigger significant effects on the timing of maximum growth rate for rice grown in Hangzhou, whereas no such effects were detected for Hainan.

## DISCUSSION

The advent of powerful genetic and molecular tools has made it possible to define the machinery of development in terms of gene action and interaction of individual loci or QTL (Cheverud* et al.* 1996; Yan* et al.* 1998). A number of statistical methods have been proposed to map QTL underlying complex phenotypes primarily on the basis of goodness of fit to observational data rather than on the basis of any biological mechanism (Lander and Botstein 1989). However, the degree of the successful identification of QTL depends on the power of QTL mapping techniques to analyze growth data measured at many different time points and, more importantly, on the construction of a conceptual framework for integrating biological principles of growth laws into developmental processes through powerful statistical models.

We have framed a new statistical strategy for QTL mapping through specific incorporation of biological laws behind the phenotypic expression of complex traits. The new strategy, termed *functional mapping* (Ma* et al.* 2002; Wu* et al.* 2002a), displays greater potential for improving the precision, power, and resolution of QTL mapping in any organism, as compared with traditional mapping approaches (Ma* et al.* 2002). Functional mapping grounds theoretical genetic models in integrated developmental networks or processes and, consequently, has a direct impact on the interface of development, ecology, and evolution. With this new strategy, two major challenging questions in contemporary developmental biology can be well addressed, which regard the existence of particular regulatory genes guiding growth differentiation during an entire biological process and the alteration of the genetic architecture of a complex trait over developmental times.

In this article, we extended our theoretical model for functional mapping to study dynamic patterns of genetic effects of QTL governing growth curves and unravel the genetic machinery of organismic responses to different environments during the course of growth and development. This extended model was employed to study the genetic architecture of plant height growth trajectories in a plant model system—rice. We have successfully detected five QTL on different chromosomes that exert significant impacts on growth processes, with similar patterns across sharply contrasting environments. At four QTL, the tall parents contribute most favorable alleles to their progeny at these detected QTL. All these five QTL are expressed differently throughout entire growth processes between the two environments. These results suggest that it is possible to make efficient marker-assisted selection for rice varieties with reduced-height growth and, therefore, stronger resistance to wind and rain damage and higher grain yields through allocating more resources to the reproductive organs than to vegetative tissues (Sakamoto* et al.* 2003).

In a comparison with previous results for growth QTL based on the same material by traditional interval mapping (Yan* et al.* 1998), we have several interesting findings. First, Yan *et al.* detected 11 QTL for plant height growth in rice. Many of these QTL were not detected by our functional mapping model, but do display clear, although nonsignificant, peaks in our LR profile (Figure 3). For example, Yan *et al.* detected a QTL bracketed by markers Amy3D and E-RZ66 on chromosome 8 where our LR profile displays a peak. But according to our criterion, this QTL is not significant. This discrepancy may be due to the low criterion these authors have used. Second, many QTL detected from our model were not detected by traditional interval mapping (Yan* et al.* 1998). But all of our QTL were detected to exist at similar locations for different rice materials with a larger sample size (Li* et al.* 2003). For example, Li *et al.* found a significant plant height QTL between RG345 and RG381 on chromosome 1 that fails to be detected by Yan *et al.* but is supported by our model. This may suggest that our model displays greater power to detect growth QTL of small effects. Third, as discussed in Ma* et al.* (2002) and Wu* et al.* (2004a), our functional mapping incorporating developmental principles allows for the tests of a number of biologically meaningful hypotheses at interplay between genetics and development.

Growth for all organisms undergoes complex developmental stages. For example, rice growth includes vegetative (germination to panicle initiation), reproductive (panicle initiation to heading), and grain filling and ripening or maturation (heading to maturity; Moldenhauer and Slaton 2003). Each of these stages determines grain yield by influencing the number of panicles per unit land area, the average amount of grain produced per panicle, and the average weight of the individual grains. Our model can detect the genetic control of QTL over the length of any of these stages and any development events. Coupled with the physiological and developmental changes in various stages, our model will gain better insights into the mechanistic bases of seed formation and grain yield regulation.

Our model also has great implications for evolutionary studies. The evolution of complex organisms such as animals and plants does not occur by direct transformation of adult ancestors into adult descendants. Rather, any evolutionary change includes modifications or alterations of a series of developmental events occurring at different time periods during ontogeny (Rice 1997). The synthesization of this fundamentally important argument into the *evo-devo* framework (Raff 2000; Arthur 2002) needs knowledge of how genetic factors regulate developmental processes in an organism from embryo to adult to better adapt itself to different environments. While in the past the dynamic change of genetic control over time was estimated by traditional quantitative genetic approaches that partition total genetic variances into additive, dominant, and/or epistatic variance components (Atchley 1984; Kirkpatrick* et al.* 1994; Atchley and Zhu 1997; Pletcher and Geyer 1999; Jaffrézic and Pletcher 2000), our functional mapping model demonstrates tremendous power to precisely characterize these genetic components at any particular developmental stages and relate them to the genetic control mechanisms for other life-history traits. With this model, we are closer to addressing historically difficult questions of how small genotypic modifications are translated into phenotypic changes during evolution and how microevolutionary changes contribute to macroevolutionary events. Our model presented in this article is being implemented in our web-based software, called FunMap (Ma* et al.* 2004), to facilitate other scientists' investigations of the genotype × environment interactions for growth trajectories.

## Acknowledgments

We thank three anonymous reviewers for their constructive comments on the article that led to a better presentation of it. The preparation of this article is supported by the University of Florida Research Opportunity Fund (02050259) to R.W. The publication of this article is approved as journal series R-08963 by the Florida Agricultural Experiment Station.

## Footnotes

Communicating editor: G. Gibson

- Received May 21, 2004.
- Accepted August 10, 2004.

- Genetics Society of America