Abstract
To better utilize limited resources for their survival and reproduction, all organisms undergo developmental changes in both body size and shape during ontogeny. The genetic analysis of size change with increasing age, i.e., growth, has received considerable attention in quantitative developmental genetic studies, but the genetic architecture of ontogenetic changes in body shape and its associated allometry have been poorly understood partly due to the lack of analytical tools. In this article, we attempt to construct a multivariate statistical framework for studying the genetic regulation of ontogenetic growth and shape. We have integrated biologically meaningful mathematical functions of growth curves and developmental allometry into the estimation process of genetic mapping aimed at identifying individual quantitative trait loci (QTL) for phenotypic variation. This model defined with high dimensions can characterize the ontogenetic patterns of genetic effects of QTL over the lifetime of an organism and assess the interplay between genetic actions/interactions and phenotypic integration. The closed forms for the residual covariance matrix and its determinant and inverse were derived to overcome the computational complexity typical of our high-dimensional model. We used a worked example to validate the utility of this model. The implications of this model for genetic research of evo–devo are discussed.
THE success in selecting superior genotypes in agricultural breeding or in designing individualized medications for a particular disease critically relies upon our knowledge about the genetic basis for the phenotypic variation of morphology and physiology throughout the entire lifetime of an organism or for a period of its life. This, ultimately, is attributed to an understanding of how the genotype is related with the phenotype at the individual gene level (Lynch and Walsh 1998). A simplified assumption, i.e., there is a direct relationship between the genotype and final phenotype, although rarely existing, has been widely used by population and quantitative genetic approaches that model breeding and evolutionary processes as changes in allele frequencies or mean trait values at those genetic loci thought to affect the trait (Raff 2000; Wolf 2002). Rice (2002) constructs a unifying theory for the link between the intricate patterns of development and processes of phenotypic evolution. However, he did not characterize any specific developmental patterns on the basis of fundamental physiological and biochemical principles and made no attempt to detect the genetic variants that cause developmental and evolutionary changes.
A comprehensive model for phenotypic variation and evolution should be incorporated by three different biological processes, allometry, ontogeny, and plasticity (Wu et al. 2003a), each of which involves complicated machineries that relate the genotype and phenotype. Allometry describes the scaling and structural relationship between different traits and, therefore, determines the coordination and function of various parts of an organism (Niklas 1994). Ontogeny concerns a series of developmental changes of an organism from a fertilized egg to an adult, which are essential for the organism to live and reproduce (Klingenberg 1998). By altering the phenotype of a trait across a range of ecological conditions, plasticity enables the organism to buffer against environmental fluctuations (Schlichting and Pigliucci 1998). It is anticipated that the fusion of these three aspects of biology within the genetic context will help to gain insights into the control and regulation mechanisms of growth, development, and adaptation. It is difficult, however, to link these three aspects from fundamental biological principles in one single step without an in-depth understanding of the integration of each pair of them. More importantly, each such pair contributes to our understanding of phenotypic variation and evolution in different ways.
In this article, we present a hyperspace or systems model for integrating allometry and ontogeny through their underlying common genetic mechanisms using a multivariate statistical approach. In the past, the dynamic change of genetic correlations between different traits has been estimated by traditional quantitative genetic approaches that partition total genetic variances into additive, dominant, and/or epistatic variance components (Atchley and Zhu 1997). 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 quantitative trait loci (QTL) (Cheverud et al. 1996). However, the degree of this success 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 derive a conceptual model and statistical algorithm for studying dynamic patterns of genetic effects of QTL that govern different but developmentally related biological processes.
METHODOLOGY
The statistical model:
Consider two complex traits, x and y, whose phenotypes change as a function of age. These two traits are thought to be developmentally related if the expression of one trait accompanies the expression of the other during ontogeny. Figure 1 illustrates such an example for the developmental relationship between stem heights and diameters for an interspecific poplar hybrid population during the first 11 growing seasons. Each tree obeys a different height–diameter profile across ages, suggesting possible involvement of particular genetic factors in shaping the allometric development of the stem.
Consider a backcross design of size n, in which there are two different genotypes at each gene. Two dynamic traits, x and y, are measured at a finite set of times, 1, . . . , T, for each individual of the backcross. Assume that all progeny are measured at the same set of time points that are evenly spaced. Each of the two traits is controlled by two linked or unlinked QTL, although the effect sizes triggered by these two QTL can be trait dependent. The two traits are developmentally related because of pleiotropic effects or the linkage between the two QTL, or both (Lynch and Walsh 1998). Let j1, j2 = 1, 2 denote two possible genotypes for each QTL in the backcross. A genetic linkage map constructed from molecular markers for the backcross is used to map such QTL.
The phenotypic values of the two traits (xi and yi) for individual i measured at time t can be described in terms of the two putative QTL by a statistical model,
(1)where 
is the indicator variable denoted as 1 if the two-QTL genotype j1j2 is considered for individual i and 0 otherwise, 
is the genotypic value of trait k (k = x, y) for QTL genotype j1j2 at time t, and eki(t) is the residual of individual i for trait k.
The likelihood function:
The likelihood function of bivariate longitudinal observations, z = (zi) = {xi(t), yi(t)}t=1T, and marker information, M, in the backcross population, affected by the two QTL, is formulated on the basis of the mixture model (Equation 10), expressed as
(2)where Ω is an unknown vector that determines the parametric family 
. Here, 
is presented by a multivariate normal distribution with the mean vector fit by component-specific parameters, 
, and the (co)variance matrix fit by common parameters, Ωv. Thus, we have 
.
In Equation 2, θ, a parameter that describes the location of QTL, determines the proportions of different mixture normals. For a mapping population, n progeny can be classified into different groups on the basis of known marker genotypes. Thus, in each of such marker genotype groups, the mixture proportion or frequency of a joint QTL genotype, 
, is progeny specific and can be expressed as the conditional probability of QTL genotype j1j2 for progeny i given its marker genotype, which is more precisely symbolized by 
. Interval mapping commonly used in QTL identification (Lander and Botstein 1989) is to use two flanking markers, M1 and M2, at the same time to predict the segregation of the QTL bracketed by the two markers. The forms of the conditional probabilities depend on the recombination fractions between the QTL and the markers (and therefore the QTL positions described by θ1 and θ2 for two QTL, respectively). Two different situations are considered, one for two QTL each at a different marker interval and the second for two QTL at the same interval. The conditional probabilities of two-QTL genotypes are the products of the corresponding conditional probabilities for each QTL in the first situation, whereas they need to be derived on the basis of four-point analysis in the second situation.
Modeling the mean vector:
A general statistical framework, called functional mapping, has been proposed to genomewide map specific QTL that determine the developmental pattern of a complex trait (Ma et al. 2002; Wu et al. 2002a,b, 2003b, 2004a,b,c). A central premise of functional mapping lies in the connection between gene action or environmental effects and organ development by parametric or nonparametric models.
For the multivariate normal distribution of QTL genotype j1j2, functional mapping models the time-dependent mean vector by
(3)
(4)
(5)where genotypic values at different time points (Equation 3) are modeled by functions hx or hy (Equation 4) and further by hx and ℓ(hx) (Equation 5) because the two traits x and y are developmentally related.
A number of mathematical models have been established to describe the developmental process of a biological phenotype. For example, a series of growth equations have been derived to describe growth in height, size, or weight (von Bertalanffy 1957; Richards 1959) that occurs whenever the anabolic or metabolic rate exceeds the rate of catabolism. On the basis of fundamental principles behind biological or biochemical networks, West et al. (2001) have mathematically demonstrated the universality of these growth equations. Mathematical functions that describe age-dependent growth are generally expressed for trait k as
(6)where Λk is a set of parameters that fit the growth curve of trait k = x or y. If different joint genotypes at the two QTL have different combinations of the curve parameters, 
, this implies that the QTL plays a role in governing the difference of growth trajectories.
We permit the two traits under consideration to be developmentally integrated; i.e., they are correlated across ages. For example, stem height and radial growth in a forest tree are interrelated following a sigmoid or concave curve when the tree grows (Chaplain et al. 1999). In other examples concerning allometric scaling, biological variables, such as metabolic rate, universally scale as a three-quarter power of body biomass (West et al. 1997, 1999). We use a general expression to describe the developmental correlation between two traits x and y; i.e.,
(7)where Λ is a set of parameters that describe the relationship between the two traits x and y across ages. Considering Equations 6 and
7,
we model genotypic values of the two traits x and y for QTL genotype j1j2 by parameters arrayed in 
.
Modeling the covariance matrix:
We use the antedependence model, originally proposed by Gabriel (1962), to describe the dynamic features of the residual errors of these two traits. This model states that an observation at a particular time t depends on the previous ones, with the degree of dependence decaying with time lag. If an observation at time t is independent of all observations before t − m, this antedependence model is thought to be of m order. The antedependence model is extended to fit the structure of time-dependent variance and correlation, leading to the structured antedependence (SAD) model (Núñez-Antón and Zimmerman 2000).
Assuming the first-order structured antedependence [SAD(1)] model, the relationship between the residual errors of the two traits x and y at time t for individual i can be modeled by
(8)where ϕk and ψk are the antedependence parameters caused by trait k itself or by the other trait, respectively, and εxi(t) and εyi(t) are the time-dependent innovation error terms assumed to be bivariate normally distributed with mean zero and variance matrix,
where δx2(t) and δy2(t) are termed time-dependent “innovation variances” that can be described by a parametric function such as a polynomial of time,
(9)(Pourahmadi 1999), and ρ(t) is the correlation between the error terms of the two traits, which is assumed to be a function of time t, expressed as
(10)(Pourahmadi 1999; Jaffrézic et al. 2003). It is assumed that there is no correlation between the error terms of two traits at different time points; i.e., Corr(εx(tx), εy(ty)) = 0 (tx ≠ ty).
The (co)variance matrix (Σ) of phenotypic values for traits x and y with time-dependent residuals specified by a bivariate SAD(1) model (Equation 8) can be derived (see appendix). Also, the closed forms for the determinant and inverse of Σ can also be derived. On the basis of the bivariate SAD(1) model, the (co)variance matrix Σ (and therefore its determinant and inverse) are defined by a vector of parameters arrayed in Ωv = (uk, vk, wk, λk, ϕk, ψk) for k = x, y. These parameters can be estimated by using them to model the structure of the (co)variance matrix involved in the function-mapping model as described below.
The structure of Σ, modeled by a set of parameters Ωv, displays a few favorable properties. First, if we assume that the innovative variances, δx2 and δy2, are constant over all T time points, the residual variance for trait k at time t is given by
(11)For t2 > t1, the correlation function of the SAD(1) model is written as
(12)It can be seen from Equation 11 that even for constant innovation variances, the residual variance can change with time (Jaffrézic et al. 2003). Also, according to Equation 12, the correlation function is nonstationary for the SAD model because the correlation does not depend only on the time interval t2 − t1.
Second, since we assume that the innovative correlations exist only between the two traits at the same time point but not between different time points, the residual correlation between the two traits at different times, t1 for trait x and t2 for trait y, is derived as
(13)where we assume ρ(1) = ⋯ = ρ(T) = ρ. It is seen from Equation 13 that the across-correlation functions are not symmetrical for the SAD model; i.e., the correlation between trait x at time t1 and trait y at time t2 is not equal to the correlation between trait x at time t2 and trait y at time t1.
Computational algorithm:
The EM algorithm (Dempster et al. 1977) is implemented to estimate the parameters, 
, that model the mean vector and the parameters, Ωv, that model the structure of the (co)variance matrix. These unknown parameters, plus the QTL position, can be estimated by differentiating the log-likelihood function of Equation 2 with respect to each parameter, letting the derivatives equal zero, and solving the log-likelihood functions.
In practical computation, the QTL position parameter can be viewed as a fixed 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 log-likelihood-ratio test statistic for a QTL at a particular map position is displayed graphically to generate a likelihood map or profile. The genomic position that corresponds to a peak of the profile is the maximum-likelihood estimate (MLE) of the QTL location.
Alternatively, the Nelder–Mead simplex algorithm, originally proposed by Nelder and Mead (1965), can be used to estimate the parameters that model the mean (co)variance (Zhao et al. 2004a). It is a direct search method for nonlinear unconstrained optimization. It attempts to minimize a scalar-valued nonlinear function using only function values, without any derivative information (explicit or implicit). The algorithm uses linear adjustment of the parameters until some convergence criterion is met. The term “simplex” arises because the feasible solutions for the parameters may be represented by a polytope figure called a simplex. The simplex is a line in one dimension, a triangle in two dimensions, and a tetrahedron in three dimensions, respectively. Because the (co)variance matrix Σ and its determinant and inverse have their closed forms (appendix), the computational efficiency can be improved by incorporating these forms into the computational algorithm.
After the point estimates of parameters are obtained, we derive the approximate variance–covariance matrix and evaluate the sampling errors of the estimates (
). The techniques for doing so involve calculation of the incomplete-data information matrix, which is the negative second-order derivative of the incomplete-data log-likelihood. The incomplete-data information can be calculated by extracting the information for the missing data from the information for the complete data (Louis 1982). A different so-called supplemented EM (SEM) algorithm was proposed by Meng and Rubin (1991) to estimate the asymptotic variance–covariance matrices, which can also be used for the calculations of the sampling errors for the MLEs of the parameters (
).
Model selection:
Jaffrézic et al. (2003) proposed an ad hoc approach for model selection. Their strategy is to increase the antedependence order until the additional antedependence coefficient is close to zero.
Núñez-Antón and Zimmerman (2000) proposed using Akaike information criteria (AIC) to select the best model. Hurvich and Tsai (1989) showed that AIC can drastically underestimate the expected Kullback–Leibler information when only a few repeated measurements are available. Instead, they derived a corrected AIC, expressed as
(14)where 
is the white noise variance for trait k, T is the number of repeated measurements, and r is the order of the model. The number of parameters is heavily penalized in that models selected by AICC are typically much more parsimonious than those selected by AIC (Hurvich and Tsai 1989). The model corresponding to a minimum AICC is chosen.
HYPOTHESIS TESTS
One of the most significant advantages of functional mapping is that it can ask and address biologically meaningful questions at the interplay between gene actions and development by formulating a series of hypothesis tests. Wu et al. (2004a) described several general hypothesis tests for different purposes. Although all these general tests can be directly used in this study, we here propose the three most important and unique tests for the existence of growth QTL and structure QTL and the relationships of the two QTL.
Growth QTL:
The existence of a QTL that affects growth trajectories of two traits x and y should first be tested. The two hypotheses for this test are formulated as
(15)where Ωm is the set of parameters for the mean curve under the null hypothesis (there is no QTL). The log-likelihood values L0 and L1 under H0 and H1 are calculated. The test is performed with a log-likelihood-ratio (LR) statistic:
(16)An empirical approach for determining the critical threshold is based on permutation tests. 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.
The LR statistic is plotted against test locations of two QTL through the genome to form a landscape in which a peak LR value corresponds to the chromosomal positions of QTL. To test whether the two QTL have a pleiotropic effect on both traits x and y, we need first to estimate the curve parameters of growth trajectories for each trait, 
and 
, on the basis of Equation 6, and then formulate the following null hypotheses:
(17)and
(18)Only when both of these null hypotheses are rejected can the two QTL be considered to jointly pleiotropically affect each trait. Similar tests can be made to test the pleiotropic effects of individual QTL.
Allometry QTL:
Equation 7 describes the developmental relationship between two traits x and y over ontogeny. The growth of one trait with the second trait's growth in the form of this equation may be controlled by QTL. The test for the genetic control of this allometric growth can be made on the basis of the following hypotheses:
(19)
The significance test can be performed by calculating the LR value between the two hypotheses at the QTL positions detected. If the null hypothesis is rejected, this means that the two QTL jointly exert an effect on the allometric growth of two traits. Such QTL, called “allometry QTL,” may affect the structure and morphology of an organ.
Single vs. multiple QTL:
One of the mechanisms for genetic correlation is the linkage of different QTL on a nearby genomic location. In other words, if one QTL affects only trait x and a second QTL affects only trait y but they are located together on the same genomic region, traits x and y will be genetically correlated. Our model provides a way to test the relative importance of pleiotropy and linkage on genetic correlation.
If the two QTL are detected within the same interval, we need to determine whether they truly represent two different QTL for traits x and y, respectively, or if they can be better viewed actually as a single QTL. To do so, the following hypothesis is formulated:
The rejection of the null hypothesis implies that traits x and y are genetically correlated because of two different but linked QTL.
RESULTS
We use the newly developed model to analyze a real data set in forest trees. The plant materials used were derived from the triple hybridization of Populus. A Populus deltoides clone (designated I-69) was used as a female parent to mate with an interspecific P. deltoides × P. nigra clone (designated I-45) as a male parent (Wu et al. 1992). Both P. deltoides I-69 and P. euramericana (P. deltoides × P. nigra) I-45 were selected at the Research Institute for Poplars in Italy in the 1950s and were introduced to China in 1972. In the spring of 1988, a total of 450 1-year-old rooted three-way hybrid seedlings were planted at a spacing of 4 × 5 m at a forest farm near Xuzhou City, Jiangsu Province, China. The total stem heights and diameters measured at the end of each of 11 growing seasons are used in this example.
A subset (90) of genotypes randomly selected from the 450 hybrids was used to construct parent-dependent genetic linkage maps with random amplified polymorphic DNAs, amplified fraction length polymorphisms, and intersimple sequence repeats on the basis of a two-way pseudo-test backcross design (Yin et al. 2002). These maps provide a basic framework for mapping growth-specific QTL with our bivariate model. As an example, our QTL analysis is based on the genetic map constructed from heterozygous markers segregating in the parent P. deltoides.
Wu et al. (2002a) and Ma et al. (2002) plotted stem height and diameter growth against age for all mapped poplar progeny. The growth of these two stem traits can be well fit by a logistic equation expressed as
(20)where a is the asymptotic or limiting value of g when 
, a/(1 + b) is the initial value of g when t = 0, and c is the relative rate of growth (von Bertalanffy 1957). When plotting stem diameter growth against height growth across ages, we find that stem height and radial growth are interrelated in a particular way when a tree grows. Figure 1 illustrates the plot of stem diameter against height growth across different ages for 11-year-old trees from the pseudo-test backcross poplar hybrid progeny. A solid mathematical theory has been established to describe the age-dependent relationship between stem height (H) and diameter (D) in different forest trees. We consider a concave curve function,
(21)to describe the relationship between these two growth traits across ages, where α, β, and γ are the parameters that determine the shape of the across-age height–diameter profile.
The height–diameter profile across different ages for each individual tree derived from a mapping population between two Populus species.
Given considerable variation in growth trajectories (Ma et al. 2002) and in the across-age H–D relationship among different poplar individuals (Figure 1), we hypothesize that specific QTL each with two genotypes denoted by Qq and qq are involved in determining the ontogenetic genetic correlation between stem height and diameter growth. Because of a modest sample size used, we model only the effect of one QTL on growth trajectories to avoid the overparameterization of a two-QTL model. For an underlying QTL genotype j, the parameters that model two developmentally correlated traits are arrayed by 
in terms of the logistic curve for diameter growth (Equation 20) and concave curve (Equation 21).
To use the SAD model for mapping QTL that are responsible for both height and diameter growth trajectories in this Populus data set, we need to determine the best order of the model by calculating the AICC values on the basis of different orders, using Equation 14. The AICC values were found to increase with order from models SAD(1) to SAD(6). Also, for the SAD model of higher orders, the first-order antedependence coefficient is qualitatively (three- to fourfold) higher than the higher-order coefficients. Thus, we think that the SAD(1) model should be a reasonable fit to the hybrid data used in this example.
By scanning all the linkage groups, two QTL on groups 9 and 10 have been detected to affect both height and diameter growth trajectories in a concave curve shape. Figure 2 illustrates a plot of the LR values between the full (there is a QTL) and reduced model (their is no QTL) across all the linkage groups. The two detected pleiotropic QTL are located at 114 cM on linkage group D9 and at 16 cM on linkage group D10 because the LR peaks (38.5 and 32.2) at these positions far exceed the genomewide critical threshold at the significance level 0.001 (22.5; Table 1) determined from 1000 permutation tests throughout all the linkage groups. These two QTL could also be detected by previous functional mapping for single-growth traits (Ma et al. 2002), but at a lower significance level of 0.01 or 0.05.
The profile of log-likelihood-ratio (LR) values between the full and reduced model estimated from the bivariate SAD(1) model for analyzing stem height and diameter growth trajectories in an interspecific poplar hybrid progeny across all the linkage groups. The linkage map used is one constructed with heterozygous markers from the P. deltoides parent (Yin et al. 2002).
The maximum-likelihood estimates (and their standard errors in parentheses) of the parameters that model the mean (co)variance structure for different genotypes at the QTL detected to affect two developmentally related traits, stem height (H) and diameter (D), in interspecific hybrids of two poplar species
The MLEs of the curve parameters and of parameters that model the structure of the variance matrix were tabulated in Table 1, along with the approximate standard errors of these estimates estimated from Louis' (1982) approach. All the parameters can be estimated with reasonable precision. The MLEs of the curve parameters in Table 1 were used to draw the growth curves of the two genotypes at each QTL for height growth and its concave relationship with diameter growth (Figures 3 and 4). The testing results for pleiotropic effects based on the null hypotheses 17 and 18 suggest that both QTL detected pleiotropically affect the curves of stem height and diameter, but these two QTL display different effects on the concave relationship between height and diameter. The QTL on linkage group 10 displays a significant effect on the height–diameter profile (Figure 4B), whereas there is no such effect for the QTL on linkage group 9 (Figure 3B).
Two different genotypes for the QTL detected on linkage group 9. Yellow curves in the background are raw data. (A) Two genotypic growth curves over ages; (B) two genotypic height–diameter profiles across ages.
Two different genotypes for the QTL detected on linkage group 10. Yellow curves in the background are raw data. (A) Two genotypic growth curves over ages; (B) two genotypic height–diameter profiles across ages.
It appears that the two detected QTL have minimal effects on stem growth in the early stage and become operational after age 4–5 years. After this age the QTL effects increase with age (Figures 3A and 4A). This developmental pattern of the QTL action appears to be concordant with that of ecological competition that begins in the stand when trees' canopies are close at years 4 and 5 (Wu and Stettler 1996).
DISCUSSION
Because biological traits are derived from developmental processes and physiological regulatory mechanisms, complex multivariate systems that undergo such processes have been a major focus in current genetic research (Rice 2002; Wolf 2002). Considering a complex interplay among numerous interacting genetic loci, developmental processes, and environmental aspects of trait expression, we need to develop a systems approach that integrates the analytic and synthetic method, encompassing both holism and reductionism (Hood and Galas 2003). Such a systems approach enables the study of all elements in a network in response to developmental or environmental perturbations. Synthetic analyses of biological information from different elements through mathematical modeling provide new insights into the operation of the system. For example, to better capture information about stemwood growth in a forest tree, one should dissolve this system into its two elements, height and radial growth, and study the developmental trajectories of these two elements and their developmental interactions during ontogeny.
In this article, we have developed such a systems approach within the context of functional mapping, as advocated by Wu et al. (2002a,b, 2003b, 2004a,b,c) and Ma et al. (2002), for quantitative traits that are developmentally integrated at the organizational level. As compared to the published models, the model reported here has two advantages. The first advantage is in its biological relevance. It provides a quantitative testable framework for studying the phenotypic integration at the interplay between gene action and development at the molecular level. The model integrates information from ontogeny and timing of development to model the genetic architecture of ontogenetic growth and allometry at different time points (Schlichting and Pigliucci 1998; Wu 1998). The second advantage of this model is that it capitalizes on the SAD model that has proven to be powerful for studying time series data in statistical problems (Núñez-Antón 1997; Pourahmadi 1999; Núñez-Antón and Zimmerman 2000). The implementation of nonstationary residual variances and correlations through the SAD model overcomes the limitations due to the stationary assumptions used in the autoregressive (AR) model. Moreover, the derivations of the closed forms of the covariance matrix and its inverse and determinant have facilitated the computation and estimation of parameters. It should be noted that the SAD model may also be used to model the ontogenetic pattern of genetic variance-correlation structure for growth (Jaffrézic et al. 2003). But because the advantage of the SAD model can be best displayed for random effects, a structured pedigree composed of individuals with different levels of relatedness will be needed for the specification of the genetic variance–covariance matrix based on the SAD model.
The integration between ontogeny and allometry can be powerful in characterizing developmental changes of overall morphology, structure, and function of an organism (Klingenberg 1998). This integration has the potential to push evolutionary studies from a traditional view—evolution is the transformation of adult morphologies—toward a more synthetic evo–devo theory; i.e., evolution should be regarded precisely as the modification of developmental (embryonic) ontogeny (Rice 1997; Raff 2000; Arthur 2002). While our previous models provide several analytical tools for genetic mapping of ontogenetic growth (Ma et al. 2002; Wu et al. 2004a,c), the model proposed in this article offers a conceptual framework for understanding the genetic mechanisms that regulate differential rates of growth of different body parts over ontogeny. Coupled with the mathematical aspects of plastic response to different environmental stimuli (see Zhao et al. 2004b,c), this model can be further merged with an emerging eco–devo theory (Dusheck 2002) to shed light on the optimal developmental machinery of an organism across a range of biotic or abiotic environments.
In general, insight into the molecular and cellular targets of complex traits would offer the unprecedented opportunity to identify more precise mechanisms for growth and development. Ample information has been gathered about individual cellular components at various ages, but this has not allowed a clear understanding of the integrated genomic circuits that control mechanisms of development, phenotypic integration, and stress responses. With the emergence of our systems approach, coupled with functional genomics, we will finally have the opportunity to study growth and development the way it is supposed to be studied, i.e., in a comprehensive manner, and to study the dynamic network of genes that determines the physiology of an individual organism over time.
APPENDIX
In what follows, we derive the matrix Σ and the closed forms of its determinant and inverse (see also Zhao et al. 2005). According to Newton (1988, p. 359), Σ is positive definite if and only if a unique low triangle matrix, L, and a diagonal matrix, Σε, exist, such that LΣL′ = Σε (Jaffrézic et al. 2003), where
with
Given that matrix L is nonsingular and matrix Σε is positive semidefinite, it can be shown that
(A1)and
(A2)where
Two longitudinal traits xi = (xi(1), … , xi(T)) and yi = (yi(1), … , yi(T)) are combined into one vector:
The quadratic form ziΣ−1z′i can then be divided into three parts,
(A3)where
Acknowledgments
We thank two anonymous referees for their constructive comments on this manuscript.
Footnotes
Communicating editor: M. W. Feldman
- Received May 6, 2005.
- Accepted September 9, 2005.
- Copyright © 2006 by the Genetics Society of America
