Comparative analysis has become widely used during the past 20 years. It consists of comparing two or more traits across species or a trait and an environmental variable, while taking phylogenetic autocorrelation into account (Harvey and Pagel 1991). Several methods have been proposed for this type of analysis (e.g., Felsenstein 1985; Grafen 1989; Diniz-Filho et al. 1998). Westoby et al. (1995a) initiated a controversy (Ackerly and Donoghue 1995; Fitter 1995; Harvey et al. 1995a,b; Rees 1995; Westoby et al. 1995b,c; see Ricklefs 1996) about what they called “phylogenetic correction,” which is the control for phylogeny in comparative analysis. They argued that comparative methods partition the explained variation of ecological data in such a way that they “allocate the maximum possible variation in a trait to phylogeny, considering only the residual as potentially attributable to ecology” (Westoby et al. 1995a, p. 531). Indeed, these methods control for the phylogenetic component in the variables when estimating the influence of present-day ecological factors. This is justified by the principle of parsimony, because related species do share phylogenetic history, and our interest is to quantify how other factors account for the distribution of characters.

As Westoby et al. (1995a) pointed out, the phylogenetic portion of the total variance of the variable of interest may contain a phylogenetic component related to ecology, which Harvey and Pagel (1991) called “phylogenetic niche conservatism.” This concept was first proposed by Grafen (1989, p. 143). Phylogenetic niche conservatism includes the shared attributes that related species have acquired because they tended to occupy similar niches during evolutionary history. For example, two sister species could develop adaptive phenotypes in their own environment, but because they are sister species, they will tend to use the same kind of ecological niches and then develop the same kind of adaptation. Therefore, this trait will be correlated to phylogeny and ecology (environment), in possibly different proportions. This is different from a purely phylogenetic trait that is under the influence of intrinsic factors. Of course, this distinction may not always be straightforward. Westoby et al. proposed to partition the variance of the data into three portions (fig. 1 from Westoby et al. 1995a, p. 531, is equivalent to Fig. 1 of the present paper showing fractions a, b, and c): a part strictly due to ecology (a), a part strictly due to phylogeny (c), and a part due to the common influence of these two factors (b), which we call “phylogenetically structured environmental variation.” This portion encompasses phylogenetic niche conservatism. No method has been proposed to date to calculate this variation partitioning in the phylogenetic context. The purpose of this paper is to propose such a method.

We will actually partition the variation in not three but four components (Fig. 1), adding the unexplained part of the variation, fraction d. It is not our intention to defend a position against the principle of parsimony. We simply want to show how variation can be partitioned to quantify the portion of variation that can be explained jointly by the two causes considered in the study.

The Partitioning Method

Variation partitioning was proposed a few years ago for multivariate ecological data showing spatial variation (Borcard et al. 1992; Borcard and Legendre 1994). The spatial variation can be expressed as a trend-surface equation or some other statistical model, based upon the geographic coordinates of the sampling sites (Legendre and Legendre 1998). Borcard and Legendre (2002) have shown how this procedure can be generalized to model spatial structures at all spatial scales that can be perceived by the sampling design. The variation of species assemblages can be decomposed into a fraction a, which is correlated with environment variables but is not spatially structured; a fraction b, which is the spatially structured component of environmental variation; a fraction c, which is not correlated with the environmental variables used in the model but is spatially structured; and an unexplained component d.

The decomposition of phylogenetic variation is a problem of the same nature. The problem here is to take into account, in the analysis, the variance related to the phylogeny. When using the method of independent contrasts (Felsenstein 1985), phylogenetic relationships are eliminated by subtraction from the variables under study instead of being expressed on their own. One possible solution would be to express the phylogeny as a distance matrix that would be included in a multiple regression equation computed on distance matrices (Legendre et al. 1994). The distance matrix would be calculated from either the original data (obtained for instance by alignment of sequences) or a derived phylogenetic tree represented in the form of a patristic distance matrix, that is, containing distances between species computed from the tree. To express the phylogenetic variance as a distance matrix also requires the expression of all the other variables as distance matrices, which leads to work on distances instead of the raw data. Dutilleul et al. (2000, table 2) showed, however, that the correlation between two data vectors x₁ and x₂ is higher than the Mantel correlation statistic between two distance matrices D₁ and D₂ derived from x₁ and x₂. Legendre (2000) has also shown that the power of the t-test of the Pearson correlation coefficient between x₁ and x₂ is higher than that of the Mantel test between D₁ and D₂. For these reasons, it seems preferable to carry out this decomposition on rectangular data tables instead of distance matrices.

The autoregressive method (Cheverud and Dow 1985; Cheverud et al. 1985; Gittleman and Kot 1990) or a maximum-likelihood-based method (Lynch 1991) can be used to estimate the level of phylogenetic inertia, partitioning the total variation of a quantitative trait into a specific and a phylogenetic component, that is, the phylogenetic inertia. Blomberg and Garland (2002) reviewed the concept of phylogenetic inertia and discussed the biològical (and other) forces that generate it. Another approach to quantify phylogenetic inertia was proposed by Diniz-Filho et al. (1998): it expresses the phylogeny in the form of principal coordinates via a principal coordinate analysis (PCoA: Gower 1966) computed from the phylogenetic distance matrix. The same approach was used by Legendre et al. (2002) in a test for host-parasite coevolution. Using numerical simulations, Diniz-Filho et al. (1998) showed that this method provides a better estimate of phylogenetic inertia than the autoregressive method when phylogenetic autocorrelation is low.

Any method for assessing phylogenetic inertia can be used in the approach that we are proposing here; the important point is to obtain an estimate of the correlation of phylogeny alone with the trait(s) under study. We will detail our method by using principal coordinates to estimate phylogenetic inertia, but it should be kept in mind that the principle would remain the same with any other method of representing the phylogeny. Our aim is not to propose a new method of comparative analysis but rather a new way of using existing approaches to obtain supplementary information.

Contrary to Diniz-Filho et al. (1998), who only used the first few principal coordinates selected by reference to a broken-stick model (Barton and David 1956; Frontier 1976), we are proposing to use all the principal coordinates, extracted from the distance matrix, that are significantly related to the dependent variable(s). There is no particular reason why the broken-stick model would preferably select principal coordinates that are of importance for the explanation of the dependent variable. Because the phylogenetic distance matrix may not always be Euclidean, eigenvalues may be negative. When this is the case, it is possible to apply correction methods, described in Gower and Legendre (1986) and Legendre and Legendre (1998, section 9.2.4), to render all eigenvalues positive. Again, the phylogenetic distance matrix can be calculated either from the raw data (e.g., sequence alignments), which avoids the reconstruction of a tree, or from a patristic distance matrix representing the phylogenetic tree; no negative eigenvalues should appear in the latter case.

Our method to partition the variation of a trait is the following. Y is the dependent variable representing the trait under study, X_E represents the matrix of ecological explanatory variable(s), and PCs stands for the matrix of principal coordinates representing the phylogeny.

Example Using Real Data

The partitioning method was applied to the study of the link between body mass and population density in 79 species of mammals (Morand and Harvey 2000; the data can be found at  http://www.pubs.roysoc.ac.uk). Body mass has been shown to be inversely related to population density (Damuth 1981, 1987; Silva and Downing 1995). The negative relationship between body mass and population density has been discussed extensively in the ecological literature (Damuth 1981, 1993; Lawton 1990; Nee et al. 1991; Silva and Downing 1995; Morand and Poulin 1998) as well as its consequences in terms of the amount of energy used (Damuth 1981). Damuth emphasized the fact that the exponent linking body size to population density should be equal to −0.75 because the energy used by an individual is positively linked to its body mass with an exponent of +0.75. The energy used by a local population of a species in a community should then be independent of body size; this is the so-called energetic equivalent rule. Criticisms against the energetic equivalent rule are based on two issues: estimation of the exponent and the possible confounding role of phylogenetic inertia.

These two life-history traits, body mass and population density, are strongly linked to the phylogeny (Morand and Poulin 1998). We will estimate here what proportions of the body mass variable are correlated with population density alone, phylogeny alone, and jointly with population density and phylogeny (which corresponds, at least in part, to phylogenetic niche conservatism). We will also quantify the unexplained portion of the variation. For our calculations, we used the data compiled by Morand and Poulin (1998) (for other references see Morand and Harvey 2000).

The data for mammal body mass and population densities were subjected to a natural logarithmic transformation (ln) to linearize their relationship. From here on, ln(body mass) and ln(density) will simply be referred to as mass and density. Departure from normality did not have to be assessed because all tests of significance were carried out through permutational procedures (999 permutations) that do not assume normality.

The regression of body mass on population density, which provided fraction a + b (Fig. 1), was highly significant (P = 0.001). Population density, containing an embedded part of phylogenetic structure, explained 65% of body mass variation (R² = 0.65).

A patristic distance matrix was derived from the phylogenetic tree (Fig. 2) by considering each branch length to be equal to one unit. A principal coordinate analysis (PCoA) was then performed on this matrix using The R Package version 4.0 (Casgrain and Legendre 2000). Each principal coordinate (which is an eigenvector, called PC hereafter) represented an amount of phylogenetic variance proportional to the associated eigenvalue. The principal coordinate analysis generated n − 1 PCs for n species; the PCs were listed in decreasing order of variance, from PC1 to PC78. A forward selection procedure (program Canoco: ter Braak and Smilauer 1998) was used to select the PCs that significantly contributed to the explanation of the body mass variable. There was no need to use a more complicated selection procedure because by construct, the PCs are all orthogonal to one another.

Eight PCs were significant and were retained in the model: PC1, PC2, PC5, PC7, PC8, PC19, PC21, and PC56. The sum of their eigenvalues represented 35% of the total variance of the PCoA-transformed patristic distance matrix, and it explained 84% of the body mass variation in the 79 species of mammals (R² = 0.84, P = 0.010).

Fraction a + b + c was found by regressing body mass on population density and the eight significant PCs. It was equal to 89% (R² = 0.89, P = 0.001). Fractions a, b, c, and d were obtained by subtraction (Fig. 3): a = (a + b + c) − (b + c) = 0.89 − 0.84 = 0.05, thus 5%; b = (a + b) − (a) = 0.65 − 0.05 = 0.60, thus 60%; c = (b + c) − (b) = 0.84 − 0.60 = 0.24, thus 24%; and d = 1 − (a + b + c) = 1 − 0.89 = 0.11, thus 11%.

If one needs to obtain the fitted values corresponding to fractions a and c, two partial regression equations have to be computed. Fraction a was found by first regressing population density on all significant PCs and computing the residuals. The partial regression was obtained by regressing body mass on these residuals; the fitted values were computed from the regression equation. An alternative method would be to compute the multiple regression of body mass on population density and all significant PCs, and pick up the partial regression coefficient corresponding to population density in the regression equation. The same value as above was found for a (R² = 0.05).

Fraction c was found in a similar way: each individual PC was regressed on density and the residuals of these eight simple linear regressions were computed. Body mass was then regressed on the eight vectors of residuals, using multiple linear regression, and the coefficient of multiple determination (R²) was computed; this coefficient was equal to c. The fitted values can be computed from the regression equation. The same value as above was found for c (R² = 0.24).

The calculation procedure may appear tedious, especially if a large number of PCs are used. We repeated the calculations using canonical redundancy analysis (program Canoco: ter Braak and Smilauer 1998). The manipulation was simpler using that program and the results were the same. When there is a single dependent variable in the analysis, redundancy analysis is simply multiple regression, and partial redundancy analysis, which is directly available in Canoco, is partial multiple regression. The dependent variable (body mass) and the explanatory variable (population density) were represented as simple vectors, whereas the phylogeny was represented as a matrix containing the eight significant PCs. To obtain the partial linear regression analyses needed to obtain fractions a and c, we simply specified one of the two sets of explanatory variables, in turn, as the matrix of covariables of the analysis. The permutation tests of significance computed by Canoco in partial regression do take the number of covariables into account. The permutational P-value for fraction a was 0.040, and 0.009 for fraction c, after 999 random permutations. It is of particular interest that fraction a was statistically significant, because it represents the relationship between ln(body mass) and ln(population density) after controlling for the effect of phylogenetic autocorrelation.

The slope of the model II functional relationship between population density (controlling for the phylogeny) and body mass (also controlling for the phylogeny), computed using major axis regression, is −0.76. Considering its 95% confidence interval (−1.06, −0.54), this slope does not significantly differ from the value −0.75 proposed by Damuth (1981) as the true functional relationship between ln(body mass) and ln(population density) in mammals. Our results thus support Damuth's equivalent energetic rule.

Discussion

The expression “related to phylogeny” encompasses many potential causes of similarity among species, including phylogenetic niche conservatism and phylogenetic inertia (sensu Westoby et al. 1995a). These two mechanisms can cause related species to be similar in many traits. From a statistical standpoint, if one wants to know if the variation of trait X explains that of trait Y, then the phylogenetic dependence (autocorrelation) among the observations must be controlled for, using an appropriate comparative method, to avoid pseudoreplication (Hurlbert 1984).

Classical comparative methods (e.g., Cheverud et al. 1985; Felsenstein 1985; Grafen 1989; Garland et al. 1992; Gittleman and Kot 1990) take into account all variation linked to phylogeny. This is a sound procedure when one wants to test the significance of the correlation between traits and establish a causative relationship between them, based on a strong biological hypothesis, most statistical testing procedures requiring independence of the observations.

In a simple cross-species comparison, there is no consideration of any variation due to phylogenetic autocorrelation, thus no control for this potentially confounding variable. Regression analysis can be used if one wants to describe the pattern of covariation among several traits and then answer the question: Is trait X associated with trait Y? It is also a good procedure when one wants to predict the value of a trait from other trait(s), answering the question: What is the predictive value of Y for a given X? Comparative methods aim at answering a different question: Does trait X explain trait Y? This supposes controlling for the effects of confounding variables. In the present case, the confounding variable is the phylogeny. The method proposed here allows one to answer yet another question: What is the nature of the variation of Y explained by X?

Our method attributes the maximum amount of variation to a single well-established cause (phylogenetic autocorrelation) before invoking other independent causes (e.g., adaptations). Our proposal allows comparative biologists to quantify Harvey and Pagel's (1991) concept of phylogenetic niche conservatism for actual data tables. The need to interpret more finely the results of comparative analyses was pointed out by Westoby et al. (1995b), who argued that current phylogenetic correction methods (i.e., comparative methods) were not able to answer questions calling for a detailed partitioning of the phylogenetic variation.

In some cases, as in the example presented in this paper, most of the variation of an ecological variable can be attributed to phylogenetic niche conservatism (Fig. 3). Compared to the phylogenetic niche conservatism portion (fraction b), the functional relationship between ln(body mass) and ln(population density) that can be attributed to genetic constraints (fraction a) is small, albeit significant. The method presented in this paper allowed us to quantify these influences (fractions a and b) and obtain a more precise explanation of the variation of the trait under study.