Genetics, Vol. 168, 1751-1762, November 2004, Copyright © 2004
doi:10.1534/genetics.104.031484

A Unified Statistical Model for Functional Mapping of Environment-Dependent Genetic Expression and Genotype x Environment Interactions for Ontogenetic Development

* Department of Statistics, University of Florida, Gainesville, Florida 32611
{ddagger} Agronomy Department, University of Florida, Gainesville, Florida 32611
{dagger} Department of Agronomy, Zhejiang University, Hangzhou, Zhejiang 310029, People's Republic of China

1 Corresponding author: Department of Statistics, 533 McCarty Hall C, University of Florida, Gainesville, FL 32611.
E-mail: rwu{at}stat.ufl.edu

Manuscript received May 21, 2004. Accepted for publication August 10, 2004.

ABSTRACT

The effects of quantitative trait loci (QTL) on phenotypic development may depend on the environment (QTL x 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 x 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 x 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 x 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 x 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 x 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 x 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 x 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 x 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 kth 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

Formula 1(1)
where g(t) represents the growth at time t and {eta} and {lambda} are the constants of metabolism and catabolism, respectively. For small values of k, integration of Equation 1 leads to the growth equation,

Formula 2(2)
where g0 is the growth at t = 0. This growth function is sigmoidal (S-shaped), approaching asymptotically the value ({eta}/{lambda})1/(1–k) as t -> {infty}. 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,

Formula 2
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

Formula 3(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 = ak1/(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,

Formula 3

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 ({tau}) 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 Formula 3lj. 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

Formula 3
where Formula 3 = (Formula 31, ... , Formula 3J) are the mixture proportions that are constrained to be nonnegative and sum to unity; {phi} = ({phi}1, ... , {phi}J) are the genotype-specific parameters, with {phi}j being specific to genotype group j (j = 1, ... , J); and {nu} 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

Formula 4(4)
where the mixture proportions (Formula 4j's) denote QTL genotype frequencies in the backcross, i.e., Formula 41 = Formula 40 = 1/2, and the multivariate normal distribution of progeny i measured at {tau} time points is expressed as

Formula 5(5)
where Formula 5 is a vector of longitudinal observation composed of B x L subvectors each for a location (l) and a block (b), with yibl = [yibl(1), ... , yibl({tau})], and component-specific parameters {phi}j for genotype j are specified by Formula 5 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

Formula 5
with ujl = [ujl(1), ... , ujl({tau})] for genotype j at {tau} 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,

Formula 5
where xi is the indicator variable denoted as 1 if progeny i carries QTL genotype j and 0 otherwise; eibl(t) is the residual error that is iid normal with the mean of zero and the variance of {sigma}2blFormula 5. The errors at two different time points, t1 and t2, are correlated with the covariance of covbl(t1, t2). These (co)variances compose a ({tau} x {tau}) matrix {sum}bl for block b and location l. A more general model is to let the errors among different blocks and locations be correlated for the same, covFormula 5Formula 5, or different time points, covFormula 5Formula 5. 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 {sum}, whose elements are the common parameter {nu} of the mixture model (Equation 1), is simplified into a block-diagonal matrix with a B x L dimension and the likelihood of Equation 4 can then be rewritten as

Formula 6(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 Formula 6{ell} (with alleles M{ell} and m{ell}) and Formula 6{ell} + 1 (with alleles M{ell}+1 and m{ell}+1). Thus, the QTL genotype frequencies in the backcross population (denoted by Formula 6j) 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 Formula 6j|i, where 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 (yibl) and marker information genotyped for both sexes (Formula 6) as

Formula 7(7)


View this table:
In this window
In a new window

 
TABLE 1

The conditional probabilities of QTL genotypes given genotypes of two molecular markers bracketing the QTL in a backcross population

 
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 ujl and the (co)variance matrix {sum}bl 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 j in location l at time t is expressed by

Formula 7
where the growth parameter set Formula 7jl = (ajl, bjl, rjl, kjl) 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.

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 {sigma}2 is constant over time, and (2) the correlation decays in a proportion of {rho} purely with time interval. With the AR(1) model, we need to estimate only Formula 7, 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 ({theta}) between the QTL and its flanking markers (see Table 1). Arrayed by Formula 7, 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

Formula 7
with the derivative with respect to any element {Omega}{xi},

Formula 7
where we define

Formula 8(8)
which 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 {{Omega}, {Pi}}, where {Pi} = {{Pi}j|ibl} (the E step; Equation 8). Conditional on {Pi}, we solve for

Formula 9(9)
to get the estimates of {Omega} (the M step; Equation 9). The estimates are then used to update {Pi}, and the process is repeated between Equations 8 and 9 until convergence. The values at convergence are the maximum-likelihood estimates (MLEs) of {Omega}. The iterative expressions of estimating {Omega} 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 x 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:

Formula 10(10)

H0 states that there are no QTL affecting growth trajectories and the two genotypic curves at each location overlap (the reduced model), whereas H1 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,

Formula 11(11)
where Formula 11 and Formula 11 denote the MLEs of the unknown parameters under H0 and H1, respectively. The log-likelihood ratio (LR) is asymptotically {chi}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

Formula 12(12)
which 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 x location interactions. However, when location-specific curves with the same QTL genotypes are approximately parallel to each other, the area under curve (Ajl) is an appropriate criterion for this QTL x location interaction test, expressed as

Formula 12

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

Formula 13(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 x 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

Formula 14(14)
which 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 [t1, t2] can be tested using a regional test approach based on the areas,

Formula 14
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:

Formula 15(15)

The rejection of H0 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 x 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 x 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

Formula 15
for tIjl, and

Formula 15
for uFormula 15.


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 x 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).


Figure 1
View larger version (36K):
In this window
In a new window
Download PPT slide
 
FIGURE 1.—

Genetic linkage maps constructed from 135 RFLP and 40 isozyme and RAPD markers for 123 DH plants derived from the tall Azucena and short IR64 parents.

 
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.


Figure 2
View larger version (49K):
In this window
In a new window
Download PPT slide
 
FIGURE 2.—

Plots of plant-height growth vs. time for between-block means of the 90 DH plants containing complete marker and phenotype data derived from the tall Azucena and short IR64 parents grown in Hangzhou (A) and Hainan (B). Plant heights were also log transformed to display approximately constant variances across times (C and D).

 
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 x 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.


Figure 3
View larger version (19K):
In this window
In a new window
Download PPT slide
 
FIGURE 3.—

The profiles of the log-likelihood ratios (LR) between the full model (there is a QTL) and reduced model (there is no QTL) for rice height growth trajectories throughout the rice genome composed of 12 chromosomes (10). The genomic position corresponding to the peak of the curve is the maximum-likelihood estimate of the QTL location. Tick marks on the x-axis represent the positions of microsatellite markers on each chromosome (bar, 20 cM). The map distances between two markers are calculated using the Haldane mapping function. The critical thresholds for acclaiming the genome-wide existence of a QTL are obtained from permutation tests. The 99.9th percentile (indicated at horizonal lines) of the distribution of the maximum LR values obtained from 1000 permutation tests is used as an empirical critical value to declare genome-wide existence of a QTL at {alpha} = 0.001.

 
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 x 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.


Figure 4
View larger version (18K):
In this window
In a new window
Download PPT slide
 
FIGURE 4.—

Four growth curves of the QTL genotypes for the QTL detected on different chromosomes for rice plants grown in two contrasting locations (Hangzhou, red curves, and Hainan, blue curves). The solid curves represent QTL genotypes QQ, whereas the broken curves represent QTL genotype qq. Allele Q is the allele inherited from the tall Azucena parent and allele q is inherited from the short IR64 parent.

 

View this table:
In this window
In a new window

 
TABLE 2

The MLEs of the QTL position, QTL effects described by growth parameters Formula 15jl = (ajl, bjl, rjl), residual variance ({sigma}l), and correlation ({rho}l) in two different locations

 
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 x 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 x environment interactions for growth trajectories.


ACKNOWLEDGEMENTS
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.


LITERATURE CITED

ARTHUR, W., 2002 The emerging conceptual framework of evolutionary developmental biology. Nature 415: 757–764.[Medline]

ATCHLEY, W. R., 1984 Ontogeny, timing of development, and genetic variance-covariance structure. Am. Nat. 123: 519–540.[CrossRef]

ATCHLEY, W. R., and J. ZHU, 1997 Developmental quantitative genetics, conditional epigenetic variability and growth in mice. Genetics 147: 765–776.[Abstract]

BOX, G. E. P., and D. R. COX, 1964 An analysis of transformations. J. R. Stat. Soc. Ser. B 26: 211–252.

CARROLL, R. J., and D. RUPPERT, 1984 Power-transformations when fitting theoretical models to data. J. Am. Stat. Assoc. 79: 321–328.[CrossRef]

CHEVERUD, J. M., E. J. ROUTMAN, F. DUARTE, B. VAN-SWINDEREN, K. COTHRAN et al., 1996 Quantitative trait loci for murine growth. Genetics 142: 1305–1319.[Abstract]

CHURCHILL, G. A., and R. W. DOERGE, 1994 Empirical threshold values for quantitative trait mapping. Genetics 138: 963–971.[Abstract]

DIGGLE, P. J., P. HEAGERTY, K. Y. LIANG and S. L. ZEGER, 2002 Analysis of Longitudinal Data. Oxford University Press, Oxford.

GREGORCZYK, A. R., 1998 Plant growth model. J. Agron. Crop Sci. 181: 243–247.

HUANG, N., A. PARCO, T. MEW, G. MAGPANTAY, S. MCCOUGH et al., 1997 RFLP mapping of isozymes, RAPD and QTLs for grain shape, brown planthopper resistance in a doubled haploid rice population. Mol. Breed. 3: 105–113.

JAFFRéZIC, F., and S. D. PLETCHER, 2000 Statistical models for estimating the genetic basis of repeated measures and other function-valued traits. Genetics 156: 913–922.[Abstract/Free Full Text]

JANSEN, R. C., 2000 Quantitative trait loci in inbred lines, pp. 567–597 in Handbook of Statistical Genetics, edited by D. J. BALDING, M. BISHOP and C. CANNINGS. Wiley, New York.

KIRKPATRICK, M., W. G. HILL and R. THOMPSON, 1994 Estimating the covariance structure of traits during growth and aging, illustrated with lactation in dairy cattle. Genet. Res. 64: 57–69.[Medline]

LANDER, E. S, and D. BOTSTEIN, 1989 Mapping Mendelian factors underlying quantitative traits using RFLP linkage maps. Genetics 121: 185–199.[Abstract/Free Full Text]

LI, Z. K., S. B. YU, R. LAFITTE, N. HUANG, B. COURTOIS et al., 2003 QTL x environment interactions in rice. I. Heading date and plant height. Theor. Appl. Genet. 108: 141–153.[CrossRef][Medline]

MA, C. X., G. CASELLA and R. L. WU, 2002 Functional mapping of quantitative trait loci underlying the character process: a theoretical framework. Genetics 161: 1751–1762.[Abstract/Free Full Text]

MA, C. X., R. L. WU and G. CASELLA, 2004 FunMap: functional mapping of complex traits. Bioinformatics 20: 1808–1811.[Abstract/Free Full Text]

MACKAY, T. F. C., 2001 Quantitative trait loci in Drosophila. Nat. Rev. Genet. 2: 11–20.[Medline]

MOLDENHAUER, K., and N. SLATON, 2003 Rice growth and development, pp. 7–14 in Rice Production Handbook. Cooperative Extension Service, University of Arkansas, Little Rock, AR.

NATH, S. R., and F. D. MOORE, III, 1992 Growth analysis by the first, second, and third derivatives of the Richards function. Growth Dev. Aging 56: 237–247.[Medline]

NIKLAS, K. L., 1994 Plant Allometry: The Scaling of Form and Process. University of Chicago, Chicago.

PIEPHO, H.-P., 2000 A mixed-model approach to mapping quantitative trait loci in barley on the basis of multiple environment data. Genetics 156: 2043–2050.[Abstract/Free Full Text]

PLETCHER, S. D., and C. J. GEYER, 1999 The genetic analysis of age-dependent traits: modeling the character process. Genetics 153: 825–835.[Abstract/Free Full Text]

RAFF, R. A., 2000 Evo-devo: the evolution of a new discipline. Nat. Rev. Genet. 1: 74–79.[CrossRef][Medline]

RICE, S. H., 1997 The analysis of ontogenetic trajectories: when a change in size or shape is not heterochrony. Proc. Natl. Acad. Sci. USA 94: 907–912.[Abstract/Free Full Text]

RICHARDS, F. J., 1959 A flexible growth function for empirical use. J. Exp. Bot. 10: 290–300.[Abstract/Free Full Text]

SAKAMOTO, T., Y. MORINAKA, K. ISHIYAMA, M. KOBAYASHI, H. ITOH et al., 2003 Genetic manipulation of gibberellin metabolism in transgenic rice. Nat. Biotechnol. 21: 909–913.[CrossRef][Medline]

VIEIRA, C., E. G. PASYUKOVA, Z-B. ZENG, J. B. HACKETT, R. F. LYMAN et al., 2000 Genotype-environment interaction for quantitative trait loci affecting lifespan in Drosophila melanogaster. Genetics 154: 213–227.[Abstract/Free Full Text]

VON BERTALANFFY, L., 1957 Quantitative laws for metabolism and growth. Q. Rev. Biol. 32: 217–231.[CrossRef][Medline]

WEST, G. B., J. H. BROWN and B. J. ENQUIST, 1997 A general model for the origin of allometric scaling laws in biology. Science 276: 122–126.[Abstract/Free Full Text]

WEST, G. B., J. H. BROWN and B. J. ENQUIST, 1999 The fourth dimension of life: fractal geometry and allometric scaling of organisms. Science 284: 1677–1679.[Abstract/Free Full Text]

WEST, G. B., J. H. BROWN and B. J. ENQUIST, 2001 A general model for ontogenetic growth. Nature 413: 628–631.

WU, R. L., 1998 The detection of plasticity genes in heterogeneous environments. Evolution 52: 967–977.[CrossRef]

WU, R. L., Z-B. ZENG, S. E. MCKEND and D. M. O'MALLEY, 2000 The case for molecular mapping in forest tree breeding. Plant Breed. Rev. 19: 41–68.

WU, R. L., C.-X. MA and G. CASELLA, 2002a Joint linkage and linkage disequilibrium mapping of quantitative trait loci in natural populations. Genetics 160: 779–792.[Abstract/Free Full Text]

WU, R. L., C.-X. MA, M. CHANG, R. C. LITTELL, S. S. WU et al., 2002b A logistic mixture model for characterizing genetic determinants causing differentiation in growth trajectories. Genet. Res. 79: 235–245.[CrossRef][Medline]

WU, R. L., C. X. MA, M. C. K. YANG, M. CHANG, U. SANTRA et al., 2003 Quantitative trait loci for growth in Populus. Genet. Res. 81: 51–64.[CrossRef][Medline]

WU, R. L., C.-X. MA, M. LIN and G. CASELLA, 2004a A general framework for analyzing the genetic architecture of developmental characteristics. Genetics 166: 1541–1551.[Abstract/Free Full Text]

WU, R. L., C.-X. MA, M. LIN, Z. H. WANG and G. CASELLA, 2004b Functional mapping of quantitative trait loci underlying growth trajectories using a transform-both-sides logistic model. Biometrics 60: 729–738.[CrossRef][Medline]

YAN, J.-Q., J. ZHU, C.-X. HE, M. BENMOUSSA and P. WU, 1998 Molecular dissection of developmental behavior of plant height in rice (Oryza sativa L.). Genetics 150: 1257–1265.[Abstract/Free Full Text]

ZHAO, W., C. MA, J. M. CHEVERUD and R. L. WU, 2004 A unifying statistical model for QTL mapping of genotype x sex interaction for developmental trajectories. Physiol. Genomics 19: 218–227.[Abstract/Free Full Text]

ZIMDAHL, H., T. KREITLER, C. GOSELE, D. GANTEN and N. HUBNER, 2002 Conserved synteny in rat and mouse for a blood pressure QTL on human chromosome 17. Hypertension 39: 1050–1052.[Abstract/Free Full Text]

ZIMMERMAN, D. L., and V. NúñEZ-ANTóN, 2001 Parametric modeling of growth curve data: an overview. Test 10: 1–73.[CrossRef]

Communicating editor: G. GIBSON




This article has been cited by other articles:


Home page
GeneticsHome page
R. Yang, H. Gao, X. Wang, J. Zhang, Z.-B. Zeng, and R. Wu
A Semiparametric Approach for Composite Functional Mapping of Dynamic Quantitative Traits
Genetics, November 1, 2007; 177(3): 1859 - 1870.
[Abstract] [Full Text] [PDF]