A General Framework for Analyzing the Genetic Architecture of Developmental Characteristics

Corresponding author: Rongling Wu, 533 McCarty Hall C, University of Florida, Gainesville, FL 32611., rwu{at}stat.ufl.edu (E-mail)

Communicating editor: M. W. FELDMAN

	ABSTRACT

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

The genetic architecture of growth traits plays a central role in shaping the growth, development, and evolution of organisms. While a limited number of models have been devised to estimate genetic effects on complex phenotypes, no model has been available to examine how gene actions and interactions alter the ontogenetic development of an organism and transform the altered ontogeny into descendants. In this article, we present a novel statistical model for mapping quantitative trait loci (QTL) determining the developmental process of complex traits. Our model is constructed within the traditional maximum-likelihood framework implemented with the EM algorithm. We employ biologically meaningful growth curve equations to model time-specific expected genetic values and the AR(1) model to structure the residual variance-covariance matrix among different time points. Because of a reduced number of parameters being estimated and the incorporation of biological principles, the new model displays increased statistical power to detect QTL exerting an effect on the shape of ontogenetic growth and development. The model allows for the tests of a number of biological hypotheses regarding the role of epistasis in determining biological growth, form, and shape and for the resolution of developmental problems at the interface with evolution. Using our newly developed model, we have successfully detected significant additive x additive epistatic effects on stem height growth trajectories in a forest tree.

THE evolution of complex organisms, such as animals and plants, does not result simply from the direct transformation of adult ancestors into adult descendants, but rather involves a cascade of developmental processes that produce the new features of each generation. An increasing number of evolutionary studies have been launched to determine the genetic or developmental changes in the rate or timing of developmental processes that must take place to derive a particular phenotype from its ancestor (RICE 1997 ; RAFF 2000 ; ROUGVIE 2001 ). A general view is that the evolution of developmental processes is affected by both the environment and many genes that act singly and in interaction with each other (LYNCH and WALSH 1998 ). However, to accurately predict the direction and rate of trait evolution, a detailed genetic architecture of how genes act and interact to control various stages of development must be quantified.

The genes predisposing for a phenotypic character that displays continuous variation among individuals are referred to as quantitative trait loci (QTL). The genetic effect or variance of QTL includes two components, additive, due to the cumulation of breeding values, and nonadditive, due to allelic (dominant) or nonallelic (epistatic) interactions. Epistatic interactions between different loci can be further partitioned into different types: additive x additive, additive x dominant (or dominant x additive), and dominant x dominant. The presence of epistasis implies that the influence of a gene on the phenotype depends critically upon the context provided by other genes. In the past, the estimation of the additive and nonadditive genetic architecture of a quantitative trait was based on the phenotypes of related individuals (LYNCH and WALSH 1998 ), although this has minimal power to detect the nonadditive genetic variances, especially epistatic variance because epistasis contributes little to the resemblance among relatives (CHEVERUD and ROUTMAN 1995 ).

The advent of DNA-based linkage maps opens a novel avenue for precisely estimating the genetic architecture of developmental traits (VAUGHN et al. 1999 ). Current statistical methods proposed to detect the main and interaction effects of QTL are based on the phenotypes of a quantitative trait measured at a limited set of landmark ages. More recently, WU et al. 2002 , WU et al. 2003A , WU et al. 2003B and MA et al. 2002 have derived a powerful functional mapping method for estimating the dynamic changes of QTL effects during a course of ontogenetic growth through the implementation of universal growth laws (WEST et al. 2001 ) and the structured residual (co)variance matrix among different time points (see KIRKPATRICK and HECKMAN 1989 ; KIRKPATRICK et al. 1990 , KIRKPATRICK et al. 1994 ; PLETCHER and GEYER 1999 ). This method has proven to be statistically more powerful and more precise because of a reduced number of parameters being estimated and to be biologically more meaningful due to the consideration of biological principles underlying trait development (WU et al. 2002 ). However, this model has not incorporated the estimation process of epistatic interactions and, thus, cannot examine the role of the entire genetic architecture in developmental trajectories.

In this article, we extend the functional mapping method to map any QTL (including additive, dominant, and epistatic) that transforms allelic and/or nonallelic effects into final phenotypes during a continuous process of development represented as ontogenetic trajectories or a path through phenotype-time space (ALBERCH et al. 1979 ; WOLF et al. 2000 ). We derive special procedures to estimate and test the impact of epistasis on trait growth because a growing body of evidence now shows that epistasis plays a more important role in determining developmental changes than originally thought (RICE 1997 , RICE 2000 ; WOLF et al. 2000 ). Epistasis can trigger an effect on the evolution of development across different levels of biological organization and these include the molecular mechanisms of gene expression and genetic architecture, the evolution of sex and recombination, the genetic coadaptation and its associated outbreeding depression, adaptive evolution, and the very process of speciation (WOLF et al. 2000 ). We use a maximum-likelihood-based method, implemented with the expectation-maximization (EM) algorithm, to estimate QTL locations and genetic effects on growth differentiation. Compared with current mapping methods, our method of incorporating growth trajectories tends to be more powerful and more precise in QTL detection and effect estimation, as demonstrated in an example using forest tree data. In practice, our method is economically more feasible than previous methods because it needs a smaller size of genotyped samples to obtain adequate power for QTL detection through the use of repeated measurements for each individual. It can be anticipated that the method proposed in this article will have potential implications for understanding the origin and evolution of development and the contributions of epistatic effects to evolutionary changes in the process of development.

	GROWTH EQUATIONS AND MIXTURE MODEL

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

Growth equations:
When growth g is plotted against time t, different forms of growth curves will appear. Among these forms, a logistic growth curve (also referred to as the sigmoid curve of growth; NIKLAS 1994 ) is one of the most ubiquitous, having been derived from fundamental physiological and physical principles (WEST et al. 2001 ). The logistic growth curve can be mathematically described by

(1)

where a is the asymptotic value of g when t , a/(1 + b) is the value of g at t = 0, and c is the relative rate of growth (BERTALANFFY 1957 ). The logistic growth curve consists of two phases, exponential and asymptotic. The overall form of the curve is determined by different combinations of parameters a, b, and c.

In evolutionary biology, a question of how a population evolves on a logistic curve is determined by how selection acts on the growth and by the local geometry of the curve itself. Some geometric properties of the growth curve have straightforward biological interpretations. For example, the slope of the logistic curve at any given time point measures the degree to which the value of growth is sensitive to a change in age:

(2)

Such a slope represents the rate of growth at a given time. Thus, if the slope at a point is low, then that value of growth is locally buffered against age changes. The rate of growth drops off linearly as the overall size approaches some limit.

From a growth curve, we can derive the timing (t_I) of the inflection point, at which the exponential phase ends and the asymptotic phase begins (NIKLAS 1994 ). For the logistic curve, t_I is derived as

(3)

The inflection point is thought to play an important role in shaping ontogenetic growth and development.

The area under the logistic curve at an interval [t₁ t₂] describes the capacity of a given organism to grow over time. Such an area is the integral of the logistic curve, expressed as

(4)

Quantitative genetic model:
We start with a simple F₂ population of size n derived from two homozygous lines. Consider two segregating QTL responsible for a quantitative trait, _k and _l, with three genotypes, Q_kQ_k, Q_kq_k, q_kq_k, and Q_lQ_l, Q_lq_l, q_lq_l, respectively. Three genotypes at a QTL are denoted by j_k that takes 2, 1, or 0 depending on the number of capitalized alleles. The genotypic value of a two-QTL genotype (j_kj_l) can be expressed by a linear model,

(5)

where µ is the overall mean; _k, ß_k and _l, ß_l are the additive and dominant effects of the two QTL, respectively; and i, i_ß, i_ß, and i_ßß are the epistatic effects between the two QTL due to additive x additive, additive x dominant, dominant x additive, and dominant x dominant interactions, respectively. The dummy variables in Equation 5 are defined as

with w = x_kx_l, w_ß = x_kz_l, w_ß = z_kx_l, and w_ßß = z_kz_l. These two QTL can be mapped using a genetic linkage map constructed from molecular markers. There are two possibilities for the locations of the two QTL: (1) the two QTL are located on two different marker intervals or (2) the two QTL are located on the same marker interval. Consider QTL _k bracketed by two flanking markers _u and _u+1. The recombination fractions are denoted by r_u, r_k1, and r_k2, respectively, between the two markers, between _u and _k, and between _k and _u+1. The conditional probability of a given QTL genotype, conditional upon the marker genotypes for F₂ progeny i, can be generally expressed as

which depends on the location of the QTL on the marker interval, characterized by r_k1 and r_k2. Considering all possible two-marker genotypes and QTL genotypes, _ijk forms a (9 x 3) matrix. If a second QTL _l is located on a different marker interval [_v, _v+1], the conditional probabilities of a two-QTL genotype given marker intervals for progeny i are the product of the corresponding probabilities of a one-QTL genotype, i.e.,

(6)

which forms a (81 x 9) conditional probability matrix. If two QTL are located on the same marker interval, the conditional probability is expressed as

(7)

In Equation 7, denote r_k1, r_k2, and r_k3 to be the recombination fractions between marker _u and QTL _k, between QTL _k and _k+1, and between _k+1 and marker _u+1, respectively. For two QTL on the same interval, _ijkjk+1 forms a (9 x 9) matrix.

Logistic mixture model for mapping epistatic QTL:
Unlike traditional statistical models, in which marker information is associated with phenotypic values measured at one time point, our model intends to map QTL for an infinite-dimensional trait expressed as a function of time (KIRKPATRICK and HECKMAN 1989 ; KIRKPATRICK et al. 1990 , KIRKPATRICK et al. 1994 ; PLETCHER and GEYER 1999 ). Yet, the modeling of the functional relationship of a trait with time needs to be based on the measurements made at a finite number of landmark ages. It is reasonable to assume that the phenotypes of an infinite-dimensional trait measured at all time points 1, ... , m for each putative QTL genotype group follow a multivariate normal density,

where g_jkjl is the vector of the expected genotypic values of the trait measured for m times for a two-QTL genotype j_kj_l at _k and _l and is the residual variance-covariance matrix of the phenotypes measured at different times. Assuming that the two putative QTL jointly affect the growth process, g_jkjl can be modeled by a growth equation. For the logistic curve of Equation 1, we have

where each group of growth parameters (a, b, c) corresponds to a different QTL genotype. To increase the model's flexibility, the residual (co)variance matrix need to be structured using the first-order autoregressive [AR(1)] model (DAVIDIAN and GILTINAN 1995 ), expressed as

(8)

in which we assume variance stationarity, i.e., there is the same residual variance (²) for growth at different ages, and covariance stationarity, i.e., the covariance of growth between different ages, decreases proportionally (in correlation ) with increased time interval (cf. PLETCHER and GEYER 1999 ). There are two advantages when the structured matrix (8) is used. First, an explicit expression of the determinant and inverse of can be derived, which facilitates parameter estimation. Second, with such an expression, the growth-model-based mapping approach can be applied for an arbitrary number of time points.

The assumption of variance stationarity can be satisfied by transforming both sides (TBS) of the growth Equation 1, as proposed by CARROLL and RUPPERT 1984 . The transformation at the left side of Equation 1 can lead to a homogeneous variance over times, whereas the transformation at the right side of Equation 1 can preserve the biological properties of growth parameters (a, b, c). Thus, with the TBS model, the favorable advantages of structuring according to (8) can be preserved.

We formulate the likelihood function of the phenotypic data with m-dimensional measurements as

(9)

where the vector contains unknown parameters for the QTL effects, QTL position, and residual (co)variances. The position parameters r_k1 and r_l1 depend on whether the two QTL are tested on different intervals or the same interval.

The EM algorithm:
The maximum-likelihood estimates (MLEs) of the unknown parameters under a two-QTL model can be computed by implementing the EM algorithm (DEMPSTER et al. 1977 ; LANDER and BOTSTEIN 1989 ). We have incorporated the growth law (1) into the mixture-based likelihood function (9) and derived the log-likelihood equations to estimate . In the E step, calculate the expected conditional (posterior) probability of a two-QTL genotype j_kj_l given marker genotypes for progeny i,

In the M step, these posterior probabilities are used to solve the unknown parameters on the basis of the log-likelihood equations. The E and M steps are iterated until the estimates converge.

In practical computations, the QTL position parameters can be viewed as nuisance parameters because two putative QTL can be searched at given positions throughout the entire linkage map. The amount of support for the QTL at particular map positions is often displayed graphically through the use of likelihood maps or profiles, which plot the likelihood-ratio test statistics as a function of map positions of the two putative QTL.

	HYPOTHESIS TESTS

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

Different from traditional mapping approaches, our functional mapping for epistatic QTL allows for the tests of a number of biologically meaningful hypotheses. 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 ages.

Global test:
Testing whether specific QTL exist to affect the shape of growth trajectories is a first step toward the understanding of the genetic architecture of complex phenotypes. The genetic control over entire growth trajectories can be tested by formulating the following hypotheses:

(10)

H₀ states that no QTL affect growth trajectories (the reduced model), whereas H₁ proposes that such QTL do exist (the full model). The test statistic for testing the hypotheses (10) is calculated as the log-likelihood ratio of the reduced to the full model,

(11)

where and denote the MLEs of the unknown parameters under H₀ and H₁, respectively. The LR is asymptotically ² distributed with 9 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.

We can also test the global effects of different genetic components, additive, dominant, and epistatic, on the shapes of entire growth curves. The hypothesis for testing the additive effect of QTL _k on overall growth curves can be formulated as

(12)

which is equivalent to testing the difference of the full model with no restriction and the reduced model with a restriction:

(13)

Thus, the data can be fit by one less unknown parameter under the reduced model (H₀) of (12) than under the full model (H₁). An empirical approach for determining the critical threshold for the hypothesis test of (12) is based on simulation studies. Phenotypic data following a multivariate normal density are simulated for different groups of QTL genotypes whose time-dependent expected values are restricted using Equation 13. These simulated data that include no additive effect due to QTL _k are analyzed. The threshold value is determined on the basis of the distribution of the likelihood ratios (LRs) obtained from simulation replicates.

The test for the dominant effect, ß_k(t), of QTL _k is equivalent to testing the difference of the full model with no restriction and the reduced model with a restriction:

(14)

Similarly, the additive and dominant effects of QTL _l at a time point t can be tested using the following restrictions, respectively:

(15)

(16)

The test for the epistatic effects between the two QTL is equivalent to testing the differences of the full model with no restriction and the reduced model with a restriction,

(17)

for the additive x additive effect;

(18)

for the additive x dominant effect;

(19)

for the dominant x additive effect; and

(20)

for the dominant x dominant effect. In fact, these formulations of the restrictions for testing epistatic interactions of different kinds on growth trajectories are a simple extension of CHEVERUD and ROUTMAN's (1995) epistatic model for specifying the epistasis for a stationary trait. Simulation studies are performed to determine the critical thresholds for hypotheses 14–20.

Local test:
The local test can test the significance of the main (additive or dominant) effect of each QTL and the interaction (epistatic) effect between the two QTL on growth measured at a time point (t*) of interest. The tests of additive and dominant effects of individual QTL and their epistatic effects can be made on the basis of the corresponding restrictions given in Equation 13 HREF="#FD14">Equation 14Equation 15Equation 16Equation 17Equation 18Equation 19 HREF="#FD20">Equation 20. For example, the hypothesis for testing the additive effect of QTL _k on growth at a given time t* can be formulated as

(21)

which is equivalent to testing the difference of the full model with no restriction and the reduced model with a restriction:

(22)

Regional test:
It is likely that an important developmental event often occurs in a time interval rather than simply at a time point. The question of how QTL exert their effects on a period of growth process [t₁, t₂] can be tested using a regional test approach based on the areas (G_jkjl [t₁, t₂], Equation 4) covered by growth 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. The types of restriction used are similar to Equation 13Equation 14Equation 15Equation 16Equation 17Equation 18 HREF="#FD19">Equation 19Equation 20, depending on the additive effects, dominant effects, or epistatic effects of different kinds.

Interaction test:
The effects of QTL may change with age, which suggests the occurrence of QTL x age interaction effects on the growth process. The differentiation of g(t) with respect to time t represents a slope of the growth curve (growth rate, Equation 2). If the slopes at a particular time point t* are different between the curves of different QTL genotypes, this means that significant QTL x age interaction occurs between this time point and the next. The test for QTL x age interaction can be formulated with the restriction

(23)

for the additive effect of QTL _k. The tests for QTL x age interactions due to the other genetic effects can be formulated using Equation 14Equation 15Equation 16 HREF="#FD17">Equation 17Equation 18Equation 19Equation 20.

The effect of QTL x age interaction on growth can be examined during the entire growth trajectories. The global test of QTL x age interaction due to the additive effect of QTL _k can be formulated with the restriction

(24)

The restriction (24) means that there is the same slope at every time point between the two logistic curves of QTL genotypes Q_kQ_k and q_kq_k, thus suggesting that the additive effect of _k does not lead to significant QTL x age interaction on entire growth trajectories. Similar global tests of QTL x age interactions due to the other genetic effects can also be made, depending on the types of restrictions as shown in Equation 14Equation 15 HREF="#FD16">Equation 16Equation 17Equation 18Equation 19Equation 20.

Test for the timing of development:
During its ontogenetic growth, an organism would experience various developmental events. The genetic determination of the timing of development sheds light on the theoretical integration of evolution and development (RAFF 2000 ; ROUGVIE 2001 ). Using our functional mapping model, the genotypic differences in the timing (t_I) of the inflection point of maximum growth rate can be tested. According to Equation 2, the test for such a genotypic difference due to the additive effect of QTL _k is based on the restriction,

(25)

The tests of the control of other genetic components over the timing of the inflection point can be similarly made.

	A CASE STUDY

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

Materials:
The power of our statistical model for mapping QTL affecting growth trajectories was demonstrated by a case study in poplar trees. A Populus deltoides clone (designated I-69) was used as a female parent to mate with an interspecific P. deltoides x P. nigra clone (designated I-45) as a male parent (WU et al. 1992 ). A total of 450 1-year-old rooted hybrid seedlings from this cross were planted at a spacing of 4 x 5 m at a forest farm near Xuchou City, Jiangsu Province, China. The total stem heights and diameters were measured at the end of each of the growing seasons for each tree. Two parent-specific genetic linkage maps each composed of 19 linkage groups (roughly representing 19 haploid chromosomes in poplar) were constructed from restriction fragment length polymorphism, amplified fragment length polymorphism, and microsatellite markers for this hybrid progeny (YIN et al. 2002 ) and used for the genetic mapping of QTL affecting complex traits of economical importance in forest trees.

Methods:
YIN et al. 2002 used GRATTAPAGLIA and SEDEROFF's (1994) pseudo-test backcross strategy to construct two linkage maps each corresponding to a parent. Each testcross marker for these two parent-specific maps is heterozygous in one parent and null in the other. Because the two parents, I-69 and I-45, are heterozygous, there is no consistent linkage phase among dominant alleles of different markers on the same linkage group; some are in a coupling phase whereas others are in a repulsion linkage phase (YIN et al. 2002 ). Thus, unlike QTL mapping in inbred-line crosses, we need to determine the correct linkage phase between the QTL and the markers flanking it for the pseudo-test backcross mapping population (LIN et al. 2003 ). Our functional mapping proposed above was modified to incorporate the uncertainty of the QTL-marker linkage phase into the likelihood function (Appendix).

Detection of QTL with significant epistasis:
The statistical model built upon a universal logistic growth law (WEST et al. 2001 ) is used to map epistatic QTL responsible for growth trajectories in poplars. All of the 19 linkage groups were scanned on a 2-cM scale for the existence of a pair of QTL at different genomic locations. We have successfully detected a few pairs of genomic locations at which two QTL interact to affect stem growth trajectories in poplar. Table 1 and Fig 1 illustrate an example in which a pair of QTL located at the ninth interval [TC/CTG-730, TC/CTG-735] and the thirteenth interval [AG/CTT-570, AG/CTT-595] of linkage group D16 (YIN et al. 2002 ) were found to show significant epistatic effects on stem height growth. The uppercase alleles of these two QTL were observed to be in a coupling phase with dominant alleles of their respective flanking coupling markers. The maximum (93.1) of the landscape of the log-LR test statistics across the linkage group (Fig 1), greater than the genome-wide threshold at the significance level = 0.05 (LR_T = 85.6) estimated from permutation tests, justifies the adequacy of a two-QTL model incorporating growth curves.

View larger version (68K): [in this window] [in a new window] [Download PPT slide]   Figure 1. The landscape of the log-likelihood ratio for two epistatic QTL across the same linkage group D16 containing 18 markers in poplar (YIN et al. 2002 ). Two QTL, one located on the ninth marker interval [TC/CTG-730, TC/CTG-735] and the other located on the thirteenth interval [AG/CTT-570, AG/CTT-595], were detected to epistatically determine the overall shape of stem height growth. The plot denotes the critical threshold for the existence of significant QTL. The peak of the landscape corresponding to the genomic locations of two QTL is indicated.

To show the advantages of our functional mapping model, the same data set was analyzed by traditional univariate and multivariate interval mapping approaches. Univariate interval mapping applied to the most differentiated heights at the oldest age (year 11) measured detected no QTL, whereas multivariate interval mapping for three representative ages (years 1, 6, and 11) suggested a marginal QTL at the significance level = 0.10 (results not shown). These results indicate that our functional mapping approach is statistically more powerful for detecting QTL from a given mapping population.

The pseudo-test backcross used here allows for only the significance tests of the additive effects of the two QTL and their additive x additive epistasis. The thresholds at the = 0.05 level for these tests were calculated by simulation studies with the restrictions 13 and 17, respectively. By comparing the maximum value of the LRs from the functional mapping approach with these thresholds, we found that the additive x additive epistasis has a significant impact on the overall differences of growth curve shapes in stem height growth, whereas these two QTL each display marginally significant effects.

To address a possible violation of the constant variance assumption in the matrix (8), we incorporate the TBS model (CARROLL and RUPPERT 1984 ) into our functional mapping framework. Similar results about the estimation of QTL positions and effects were obtained from the TBS-based mapping approach (data not shown). Yet, the TBS-based mapping approach provides more precise estimates of growth curve parameters, with sampling errors reduced by 20–50% compared to those from untransformed data.

The dynamic pattern of QTL effect:
The growth curves of height are drawn using the estimates of logistic parameters (Table 1) for four genotypes at two interactive QTL located on linkage group D16 (Fig 2). As described in HYPOTHESIS TESTS, our functional mapping approach can be used to test various genetic hypotheses related to the developmental process on the basis of estimated growth parameters. The additive effect of the QTL located on the thirteenth marker interval is significant throughout the entire growth process measured, but its sign is altered when trees develop into age 7–8 years (Fig 2). The QTL located on the eighth marker interval has nonsignificant additive effect on growth, but it interacts significantly with the QTL on the thirteenth interval. It is not surprising that significant QTL x age interactions are detected on height growth given the change of the signs of the additive and additive x additive effects (Fig 2).

View larger version (19K): [in this window] [in a new window] [Download PPT slide]   Figure 2. The growth curves of four different QTL genotypes drawn using parameter sets (a, b, c) in Table 1 for two QTL detected on the same linkage group D16. The coordinates of the inflection point for each curve are indicated by the horizontal and vertical lines. The differentiation pattern of growth curves beyond the maximum observed age (11 years), affected by the QTL, is represented by extended broken curves.

Genetic control over the inflection point:
Equation 3 describes the coordinates of the inflection point where the exponential phase ends and the asymptotic phase begins (NIKLAS 1994 ). The difference in the coordinates between different genotypes provides important information about the genetics and evolution of growth trajectories. If different growth curves predicted by a QTL have different ages at the inflection point, this indicates that the inflection point is under genetic determination. In our example of poplars, significant additive effect due to the QTL on the thirteenth marker interval and significant additive x additive effect on the timing of the inflection point are detected. The additive effect of the QTL on the thirteenth interval delays the occurrence of the inflection point by about 0.4 year, whereas the additive x additive effect causes the inflection point to occur 0.8 year earlier (Fig 2). Because the inflection point occurs at a time of maximum growth rate, the genetic control of growth trajectory implies that it can be genetically modified to increase a tree's capacity to effectively acquire spatial resources.

	DISCUSSION

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

Increasing evidence has emerged for the role of complex genetic architecture in regulating the ontogenetic development of embryological phenotypes (CHEVERUD et al. 1983 ; ATCHLEY 1984 ; ATCHLEY and ZHU 1997 ; VAUGHN et al. 1999 ; CARLBORG et al. 2003 ) and, ultimately, shaping the evolutionary process of organismic form (WOLF et al. 2000 ). However, traditional genetic approaches are limited in estimating nonadditive effects. If epistasis is assumed to be absent, as in most quantitative genetic studies, COCKERHAM's (1963) model based on a mating design provides a nice estimate of the dominant variance. WU 1996 extended Cockerham's quantitative genetic model to estimate the epistatic variances on the basis of clonal replicates. If both parents and offspring are cloned, Wu's model can estimate the entire additive x additive epistatic variance, a partial additive x dominant epistatic variance, and a partial dominant x dominant epistatic variance.

In this article, we present a statistical approach for mapping any QTL that exert various genetic effects on growth trajectories based on a genetic linkage map. Our approach is unique in that it detects and estimates genetic effects due to allelic/nonallelic actions and interactions of QTL from physiological and developmental principles of growth. This uniqueness makes our approach advantageous in two aspects and leads us to construct a conceptual framework of evolutionary developmental biology (ARTHUR 2002 ).

First, we integrate growth equations into a statistical mapping framework to map developmental QTL that guide the trajectories of organ growth and development. Separate QTL analyses of growth at different time points are not powerful to follow the dynamics of QTL effects since relationships of growth at different times are not considered. Multivariate QTL analysis combining all time points takes into account these age-dependent relationships (KOROL et al. 2001 ), but it quickly becomes intractable when the number of time points increases. By fitting the expected genetic values at different time points by growth curves (WEST et al. 2001 ) and the residual (co)variance matrix by the AR(1) model (DAVIDIAN and GILTINAN 1995 ), our approach estimates a considerably reduced number of parameters, which thus increases significantly the power to detect epistasis. Second, because biological principles are incorporated, our approach sheds better light on the integration of development and epitasis. Our approach allows for the understanding of the genetic basis for growth and development at the cutting edge of biology (RAFF 2000 ; ARTHUR 2002 ). Growth and development are two different but related biological processes. It is likely that these two processes share some common genetic mechanisms. By mapping the timing of various developmental events, e.g., the time to first flower, on the growth curves of QTL genotypes, we can test whether the same growth genes are also involved in regulating reproductive behaviors.

The estimation precision of additive, dominant, and epistatic effects on growth at particular time points, at time intervals, or during the entire growth process depends on the estimation precision of the parameters describing growth curves. Our earlier studies using a real example have demonstrated that growth parameters of different QTL genotypes can be precisely estimated (MA et al. 2002 ; WU et al. 2002 , WU et al. 2004 ). The sampling errors of the growth parameters estimated from Fisher's information index are less than one-tenth of the parameter values when a modest sample size (90) was used. Moreover, these sampling errors can be reduced if functional mapping is based on the transform-both-sides model, as advocated by CARROLL and RUPPERT 1984 , for relaxing the variance stationarity assumption in the AR(1) model. The estimation precision of parameters in the real example is consistent with that from our simulation studies. It is inferred from these results that the estimation of the genetic architecture of a quantitative trait from our functional mapping approach should be adequately precise, although direct assessments of the sampling errors for the gene effect estimation from our approach are needed. In addition, it is possible to modify our analysis model to consider nonstationary variance-covariance structures using structured antedependent models, as proposed by NUNEZ-ANTON 1997 and NUNEZ-ANTON and ZIMMERMAN 2000 . Such a modification can be viewed as an alternative to the current method based on the variance-stationary assumption.

Genetic effects are expressed differentially during ontogeny (RICE 2000 ; WOLF et al. 2000 ), as also documented in our case study of poplar trees. Epistasis is regarded as a force to determine the direction of organismic evolution. This controversial view can now be tested and assessed, using this model and implementing the estimation process of epistasis. From an evolutionary perspective, however, our approach should be modified to consider population genetic properties of QTL alleles, QTL genotypes, and the relationship between the QTL and markers. If a mapping population is constructed using a natural population, the theory of linkage disequilibrium mapping should be integrated within the functional mapping framework. Linkage disequilibrium-based mapping provides a powerful tool for fine-scale mapping of complex traits (LOU et al. 2003 ) and, thus, the combination of this mapping strategy with our functional mapping can gain better insights into the genetic basis of development and evolution.

	ACKNOWLEDGMENTS

We thank two anonymous referees for their constructive comments on this manuscript. This work is partially supported by an Outstanding Young Investigators Award of the National Science Foundation of China (30128017), a University of Florida Research Opportunity Fund (02050259), and a University of South Florida Biodefense grant (7222061-12) to R.W. The publication of this manuscript is approved as journal series no. R-09586 by the Florida Agricultural Experiment Station.

Manuscript received September 24, 2003; Accepted for publication November 24, 2003.

	APPENDIX

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

ESTIMATION OF QTL PARAMETERS IN A PHASE-UNKNOWN PSEUDO-TEST BACKCROSS
Suppose an interval is flanked by two markers, _u and _u+1, whose dominant nonalleles are in a coupling phase. A putative QTL, _k, on this interval may have different phases relative to the marker alleles. Let p_k be the probability of an uppercase QTL allele (Q_k) to be in a coupling phase with the dominant marker alleles. Thus, 1 – p_k is the probability of Q_k to be in a repulsion phase with the dominant marker alleles. A second QTL, Q_l, that can be located on different marker intervals N_v and N_v+1, or on the same marker intervals M_u and M_u+1, also has two different QTL-marker phase probabilities, p_l and 1 – p_l. The combination of the two QTL leads to four possible (unknown) phase combinations with probabilities p_kp_l, p_k(1 – p_l), (1 – p_k)p_l, and (1 – p_k)(1 – p_l) if the allelic assignment of one QTL is independent of that of the second QTL. The estimation of unknown parameters should be based on a most likely phase combination (determined by p_k and p_l) inferred from the data. The likelihood functions for the four phase combinations are denoted by L₁₁, L₁₂, L₂₁, and L₂₂, respectively.

As shown in LIN et al. 2003 , the mixture-likelihood function of these four phase combinations can be written as

where

The probabilities (^.._ijkjl) of two-QTL genotypes, conditional upon marker genotypes for progeny i, are derived for different phase combinations for the two QTL located on different marker intervals (Table A11) or for the two QTL on the same interval (Table A22). The EM algorithm has been implemented to obtain the MLEs of vector and the phase probabilities of the two QTL p_k and p_l. In Table A11 and Table A22, the conditional probabilities are given for the QTL genotypes bracketed by two coupling dominant markers or repulsion dominant markers.

View this table: [in this window] [in a new window]   Table A1. Conditional probabilities of QTL genotypes given the marker genotypes for u and u+1 under different marker-QTL linkage phases

View this table: [in this window] [in a new window]   Table A2. Conditional probabilities of the QTL genotypes of two QTL, k and k+1, bracketed by two markers, u and u+1, conditional on the marker genotypes under different marker-QTL linkage phases

	LITERATURE CITED

TOP ABSTRACT GROWTH EQUATIONS AND MIXTURE... HYPOTHESIS TESTS A CASE STUDY DISCUSSION APPENDIX LITERATURE CITED

ALBERCH, P., S. J. GOULD, G. F. OSTER, and D. B. WAKE, 1979  Size and shape in ontogeny and phylogeny. Paleobiology 5:296-317.[Abstract]

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]

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

CARLBORG, O., S. KERJE, K. SCHUTZ, L. JACOBSSON, and P. JENSEN et al., 2003  A global search reveals epistatic interaction between QTL for early growth in the chicken. Genome Res. 13:413-421.[Abstract/Free Full Text]

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. and E. J. ROUTMAN, 1995  Epistasis and its contribution to genetic variance components. Genetics 139:1455-1461.[Abstract]

CHEVERUD, J. M., J. J. RUTLEDGE, and W. R. ATCHLEY, 1983  Quantitative genetics of development—genetic correlations among age-specific trait values and the evolution of ontogeny. Evolution 37:895-905.[CrossRef]

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

COCKERHAM, C. C., 1963 Estimation of genetic variances, pp. 53–94 in Statistical Genetics and Plant Breeding, edited by W. D. HANSON and H. F. ROBINSON. Pub. 982, National Academy of Science-National Research Council, Washington, DC.

DAVIDIAN, M., and D. M. GILTINAN, 1995 Nonlinear Models for Repeated Measurement Data. Chapman & Hall, London.

DEMPSTER, A. P., N. M. LAIRD, and D. B. RUBIN, 1977  Maximum likelihood from incomplete data via EM algorithm. J. R. Stat. Soc. Ser. B 39:1-38.

GRATTAPAGLIA, D. and R. SEDEROFF, 1994  Genetic linkage maps of Eucalyptus grandis and Eucalyptus urophylla using a pseudo-testcross: mapping strategy and RAPD markers. Genetics 137:1121-1137.[Abstract]

KIRKPATRICK, M. and N. HECKMAN, 1989  A quantitative genetic model for growth, shape, reaction norms, and other infinite-dimensional characters. J. Math. Biol. 27:429-450.[CrossRef][Medline]

KIRKPATRICK, M., D. LOFSVOLD, and M. BULMER, 1990  Analysis of the inheritance, selection and evolution of growth trajectories. Genetics 124:979-993.[Abstract]

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]

KOROL, A. B., Y. I. RONIN, A. M. ITSKOVICH, J. PENG, and E. NEVO, 2001  Enhanced efficiency of quantitative trait loci mapping analysis based on multivariate complexes of quantitative traits. Genetics 157:1789-1803.[Abstract/Free Full Text]

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]

LIN, M., X.-Y. LOU, M. CHANG, and R. L. WU, 2003  A general statistical framework for mapping quantitative trait loci in nonmodel systems: issue for characterizing linkage phases. Genetics 165:901-913.[Abstract/Free Full Text]

LOU, X.-Y., G. CASELLA, R. C. LITTELL, M. C. K. YANG, and J. A. JOHNSON et al., 2003  A haplotype-based algorithm for multilocus linkage disequilibrium mapping of quantitative trait loci with epistasis. Genetics 163:1533-1548.[Abstract/Free Full Text]

LYNCH, M., and B. WALSH, 1998 Genetics and Analysis of Quantitative Traits. Sinauer, Sunderland, MA.

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]

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

NUNEZ-ANTON, V., 1997  Longitudinal data analysis: non-stationary error structures and antedependent models. Appl. Stoch. Models Data Anal. 13:279-287.[CrossRef]

NUNEZ-ANTON, V. and D. L. ZIMMERMAN, 2000  Modeling nonstationary longitudinal data. Biometrics 56:699-705.[CrossRef][Medline]

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]

RICE, S. H., 2000 The evolution of developmental interactions: epistasis, canalization and integration, pp. 82–98 in Epistasis and the Evolutionary Process, edited by J. B. WOLF, E. D. BRODIE III and M. J. WADE. Oxford University Press, London/New York/Oxford.

ROUGVIE, A. E., 2001  Control of developmental timing in animals. Nat. Rev. Genet. 2:690-701.[CrossRef][Medline]

VAUGHN, T. T., L. S. PLETSCHER, A. PERIPATO, K. KING-ELLISON, and E. ADAMS et al., 1999  Mapping quantitative trait loci for murine growth: a closer look at genetic architecture. Genet. Res. 74:313-322.[CrossRef][Medline]

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

WOLF, J. B., E. D. BRODIE, III and M. J. WADE, 2000 Epistasis and the Evolutionary Process. Oxford University Press, London/New York/Oxford.

WU, R. L., 1996  Detecting epistatic genetic variance with a clonally replicated design: models for low- vs. high-order nonallelic interaction. Theor. Appl. Genet. 93:102-109.[CrossRef]

WU, R. L., M. X. WANG, and M. R. HUANG, 1992  Quantitative genetics of yield breeding for Populus short rotation culture. I. Dynamics of genetic control and selection models of yield traits. Can. J. For. Res. 22:175-182.

WU, R. L., C. X. MA, M. CHANG, R. C. LITTELL, and S. S. WU et al., 2002  A logistic mixture model for detecting major genes governing growth trajectories. Genet. Res. 79:235-245.[CrossRef][Medline]

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

WU, R. L., C.-X. MA, W. ZHAO, and G. CASELLA, 2003b  Functional mapping of quantitative trait loci underlying growth rates: a parametric model. Physiol. Genomics 14:241-249.[Abstract/Free Full Text]

WU, R. L., C. X. MA, M. LIN, Z. H. WANG, and G. CASELLA, 2004  Functional mapping of quantitative trait loci underlying growth trajectories using a transform-both-sides logistic model. Biometrics in press.

YIN, T. M., X. Y. ZHANG, M. R. HUANG, M. X. WANG, and Q. ZHUGE et al., 2002  Molecular linkage maps of the Populus genome. Genome 45:541-555.[Medline]

This article has been cited by other articles:

				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]

				W. Zhao, H. Li, W. Hou, and R. Wu Wavelet-Based Parametric Functional Mapping of Developmental Trajectories With High-Dimensional Data Genetics, July 1, 2007; 176(3): 1879 - 1892. [Abstract] [Full Text] [PDF]

				R. Yang and S. Xu Bayesian Shrinkage Analysis of Quantitative Trait Loci for Dynamic Traits Genetics, June 1, 2007; 176(2): 1169 - 1185. [Abstract] [Full Text] [PDF]

				M. Lin, H. Li, W. Hou, J. A. Johnson, and R. Wu Modeling sequence sequence interactions for drug response Bioinformatics, May 15, 2007; 23(10): 1251 - 1257. [Abstract] [Full Text] [PDF]

R. Yang, Q. Tian, and S. Xu