Effects of Phylogenetic Signal on Ancestral State Reconstruction

One of the standard tools used to understand the processes shaping trait evolution along the branches of a phylogenetic tree is the reconstruction of ancestral states (Pagel 1999). The purpose is to estimate the values of the trait of interest for every internal node of a phylogenetic tree based on the trait values of the extant species, a topology, and, depending on the method used, branch lengths and a model of trait evolution (Ronquist 2004). This approach has been used in a variety of contexts such as biogeography (e.g., Nepokroeff et al. 2003; Blackburn 2008), ecological niche evolution (e.g., Evans et al. 2009; Smith and Beaulieu 2009), and metabolic pathway evolution (e.g., Gabaldón and Huynen 2003; Christin et al. 2008). Investigations of the factors affecting the accuracy with which ancestral character states can be reconstructed have focused in particular on the choice of statistical framework (Ekman et al. 2008) and the selection of the best model of evolution (Cunningham et al. 1998; Mooers et al. 1999). However, other potential biases affecting these methods, such as the effect of tree shape (Mooers 2004), taxon sampling (Salisbury and Kim 2001), as well as reconstructing traits involved in species diversification (Goldberg and Igić 2008), have also received specific attention. Most of these studies conclude that ancestral character state reconstruction is still not perfect, and that further developments are necessary to improve its accuracy (e.g., Christin et al. 2010). Here, we examine how different estimations of branch lengths affect the accuracy of ancestral character state reconstruction. In particular, we tested the effect of using time-calibrated versus molecular branch lengths and provide guidelines to select the most appropriate branch lengths to reconstruct the ancestral state of a trait. Current maximum likelihood (ML) and Bayesian methods used to reconstruct ancestral character states incorporate information about branch lengths during the inference. Branch lengths can be specified either in units of divergence (i.e., expected number of substitutions per site per unit of time) or in units of absolute time (e.g., millions of years). In practice, time-calibrated trees (or chronograms) are nearly always chosen over trees depicting molecular changes (or phylograms; e.g., Crespi and Teo 2002; Jones et al. 2009a; Friedman et al. 2009; Skinner and Lee 2010). The implicit argument for this practical choice is that branch lengths estimated using the DNA markers employed for phylogenetic inference should not be assumed to be correlated with the rate of phenotypic evolution (Bromham et al. 2002). Evidence of substitution rate heterogeneity has been found in different taxonomic units (e.g., rodents: Spradling et al. 2001; angiosperms: Smith and Donoghue 2008), and changes in phenotypic evolutionary rate are not uncommon between lineages (O'Meara et al. 2006). However, it has recently been shown that molecular and phenotypic rates of change can be correlated, such as in flowering plants, where part of the variation in phenotypic rates could be explained by the DNA substitution rate (Davies and Savolainen 2006). Furthermore, herbaceous plant species exhibit faster substitution rates than tree and shrub species (Smith and Donoghue 2008), which suggests an important effect of life-history characteristics, and generation times in particular, on rates of molecular evolution. Similar patterns have also been found in mammals (Nikolaev et al. 2007) and invertebrates (Thomas et al. 2010). Therefore, if the rate of molecular evolution explains part of the phenotypic variation between species, we argue that inferring ancestral character states on time-calibrated trees could have deceiving results as it may not appropriately represent the evolution of species traits. However, no simple guideline to choose between chronograms and phylograms has been proposed so far. In order to investigate the effects of molecular and phenotypic rate heterogeneity on the accuracy of ancestral state reconstructions, we performed computer simulations in R (R Development Core Team 2011) using the ape (Paradis et al. 2004), geiger (Harmon et al. 2008), and picante (Kembel et al. 2010) packages. We first simulated under a uniform birth–death process 5000 trees with 64 terminal taxa whose branch lengths corresponded to relative time after rescaling the root depth to 1 (Fig. 1a). We then multiplied each branch length by a value drawn from a normal distribution with mean equal to the parent value and standard deviation equal to 0.2 (Fig. 1b), the root value being set to one. By introducing substitution rate heterogeneity, we transformed the initial ultrametric trees into phylograms where rate variation is phylogenetically conserved among lineages (Sanderson 1997). Schematic view of the tree used during the simulations with Branch B shown in dotted line subtending Clade C shown in gray. a) The original tree is simulated using a birth–death process resulting in an ultrametric tree, which is scaled by setting the root node to 1. b) Mutation rates are then applied to each branch following an autocorrelation process. c) A change in generation time is applied to Clade C after the appearance of a key innovation along Branch B indicated by a dotted line (see text for more details). To mimic a change in substitution rates, for example, due to a change in life history, a Branch B (dotted line in Fig. 1) was randomly selected and the substitution rates of every branch descending from Branch B were modified as follows. We assumed that the change in substitution rate was not instantaneous but moved toward an optimum value (defined as 10 in all simulations) following an Ornstein–Uhlenbeck process (Butler and King 2004) with a selective force of 3 and a variance of 0.2. In practice, this meant that the mean of the normal distribution used to draw the substitution rate assigned to each branch was not equal to the ancestor rate but moved toward the optimum stochastically. This whole process converted the original tree (Fig. 1a) into a phylogram with the Clade C originating from Branch B, showing an increase in substitution rate (Fig. 1c). Penalized likelihood (Sanderson 2002) was then used to transform the phylogram into an ultrametric tree using different penalties corresponding to (i) complete parametric estimation (penalty = 0), (ii) some likelihood penalty on rate change (penalty = 10), and (iii) a global molecular clock (penalty = 10,000). To account more directly for molecular rate heterogeneity, we also used BEAST (Drummond and Rambaut 2007) to infer chronograms. As BEAST, in contrast to penalized likelihood, requires DNA data to perform the dating analysis, we simulated sequences of 1000 nucleotides using the phylogram (Fig. 1c) of each simulation replicate and the F81 model (κ= 2.0). We then ran BEAST on this DNA matrix with the HKY model of nucleotide evolution, a root age of 1, and a fixed topology for 106 generations sampling every 1000. The first 2000 trees were removed as burn-in, and we constructed a maximum clade credibility tree with TreeAnnotator (Drummond and Rambaut 2007). The burn-in level and parameters of the analysis were set based on a preliminary analysis to assure an optimal convergence of the MCMC chains. Because of the intensive computational effort required by the Bayesian analyses, only 300 of the 5000 trees simulated were analyzed with BEAST. We generated continuous characters using a Brownian motion model of character evolution with a variance parameter of 1 for every node of the phylograms (Fig. 1c) and chronograms (reconstructed with penalized likelihood and BEAST). This created a data set of characters that evolved following either the molecular or ultrametric branch lengths. We also simulated characters that evolved in a manner uncorrelated to either the molecular or ultrametric branches to accommodate situations where the evolution of a character is linked to neither the substitution rates nor the time. This was done by taking the original tree (Fig. 1a) and reassigning for each branch length a value drawn from a normal distribution centered on the original branch length. This was done 4 times for each simulated tree with 4 different standard deviations (0.05, 0.1, 0.3, 0.5) to create 4 trees having increasingly uncorrelated branch lengths. We then reconstructed the characters on the phylograms and chronograms using a widely used ML algorithm (Pagel 1994) implemented in geiger (Harmon et al. 2008). For each simulated data set, we summed over the differences between the inferred and the true ancestral states at each node to obtain one value per tree. This gave a relative measure of the error in the character state reconstructed, with higher values meaning larger differences between the inferred and the true values. We also measured the tree imbalance using the Colless index (Colless 1995), the size of the clade originating on Branch B, and the phylogenetic signal using both the K (Blomberg et al. 2003) and λ (Pagel 1999) indices. In total, we investigated 30 different scenarios. First, data were simulated on 6 different phylogenetic trees (1 chronogram, 1 phylogram, and 4 trees with uncorrelated branch lengths; see earlier description). Second, the simulated data sets were analyzed on 5 different phylogenetic trees (4 different chronograms and 1 phylogram; see earlier description). Each scenario contained 5000 simulated data sets (i.e., one per tree created; 300 selected at random for the BEAST reconstructions). We summarized the general trends in our simulations by calculating the mean error overall trees for each scenario. Our results show that the error was significantly lower when characters were reconstructed on the same tree as they were simulated rather than on another set of branch lengths (t-test, t = − 4.2, df = 8.3, P = 0.003). This is perfectly logical and demonstrates the potential discrepancy that can arise if a character is reconstructed on a tree with branch lengths that do not represent its tempo of evolution. For characters that evolved on a tree with uncorrelated branch lengths, there was no difference in accuracy between reconstructing ancestral character states on phylograms or chronograms (t-test, t = − 0.18, df = 4.6, P = 0.87). Furthermore, neither the size of Clade C nor the tree imbalance explained the reconstruction accuracy (slope = − 1.15, adjusted R2= 0.002, for clade size; slope =−0.0007, adjusted R2= 0.0003, for tree imbalance). Regardless of the tree, the phylogenetic signal, measured by K or λ, always strongly explained the error in ancestral character state reconstruction (Fig. 2; K: slope =−2.379, adjusted R2 = 0.53; λ: slope = − 20.56, adjusted R2= 0.91). This result is highly intuitive as it states that, by using branch lengths having a better fit to the character evolution (stronger phylogenetic signal), the ancestral state reconstruction will be more accurate. Relationship between mean error in ancestral character state reconstruction and phylogenetic signal. Closed circles show the reconstructions done on phylograms; triangles, the reconstructions done on chronograms (3 penalized likelihood estimations); and diamonds, the reconstructions done on BEAST chronograms. Shades of gray represent the mode of character evolution (black: simulated on the phylogram; dark gray: simulated on the chronogram; light gray: simulated on uncorrelated branch lengths). To test the effect of phylogenetic signal on ancestral character state reconstruction, we used a published phylogeny containing 230 primate species (Arnold et al. 2010) and reconstructed the ancestral states of 25 continuous characters (Jones et al. 2009b) using an ML method on both published chronogram and phylogram. The phylogeny was originally built in a Bayesian framework, allowing us to replicate the analysis by using 100 pairs of phylogram/chronogram drawn from the posterior distribution of trees. For each character, species with missing data were pruned from the tree before the analysis and we standardized the value of the characters. We then summed the absolute difference between the simulated and reconstructed character states for each node inferred on the two sets of branch lengths. This difference was divided by the total number of nodes in the tree to obtain the mean difference per node. As we do not know the correct assignment of ancestral character states, we used this mean difference between chronograms and phylograms instead of the true error as in the simulations. If, as shown in the simulations, phylogenetic signal is linked to ancestral state reconstruction accuracy, reconstructions done on pairs of trees having similar phylogenetic signals should have a smaller mean difference per node than those done on a pair of trees showing very different phylogenetic signals. Results showed that both K and λ explain the mean difference in ancestral character state reconstruction, with λ showing a stronger relationship (Fig. 3; K: slope = 0.014, adjusted R2= 0.002; λ: slope = 0.21, adjusted R2= 0.19). Another important factor was taxon sampling. The mean difference between reconstructions was smaller when more taxa were present in the analysis (slope =−7.3 × 10−5, adjusted R2= 0.02). It has been shown that for maximum parsimony methods, larger taxon sampling does not necessarily lead to a more accurate reconstruction (Li et al. 2008). On the contrary, our results confirm previous findings that taxon sampling matters for ML and Bayesian methods for accurate ancestral character state reconstruction (Heath et al. 2008). Relationship between mean differences in ancestral state reconstruction and phylogenetic signal in the primate data set. Gray dots show for each pair of tree the difference in λ between the phylogram and chronogram, whereas black dots are for K. The black and gray lines show slopes of the relationship. The results of the simulations and empirical analyses show that using a phylogram or a chronogram can lead to a very different conclusion about the ancestral states of a character. This was true regardless of the character evolutionary history. Usually, genes used to build phylogenetic trees are selected to reflect mostly neutral evolution, which is supposed to have no effect on the phenotype. For comparative analysis such as ancestral character state reconstructions, this would mean that using a phylogram or a chronogram should not matter. But this reasoning misses one very important point. Substitution rate is a function of many variables, such as population size, generation time, DNA repair efficiency, metabolism, and mode of reproduction (Ballard 2000; Foltz et al. 2004; Fontanillas et al. 2007; Thomas et al. 2010). With many of these variables under selection, it is expected that variation will appear between organisms. As mutations are needed for any novelty to evolve, a higher substitution rate will increase the chance of appearance and subsequent evolution of a trait. Under this paradigm, it is thus most probable that changes in many phenotypic traits followed closely changes in the rate of DNA substitution and not a clock-like evolutionary rate (Davies and Savolainen 2006; Smith and Beaulieu 2009; Seligmann 2010). In such case, ancestral character state reconstruction will without doubt be less accurate if done on chronograms than on phylograms. The phylogenetic signal of a character measures the level of dependency in the trait value of species that is due to their phylogenetic relatedness (Revell et al. 2008). For both indices used in this study, a value of 1 is expected for a character evolving following a Brownian motion. The Brownian motion was introduced in phylogenetics to calculate independent contrasts (Felsenstein 1985). It is now used extensively in comparative methods, including estimation of ancestral character states, due to its simplicity and its connection with the genetic drift theory in population genetics (Martins 1994). Because of its inherent neutrality, it is increasingly seen as the standard null model for continuous character evolution (Butler and King 2004, Salamin et al. 2010). By choosing to reconstruct ancestral states with the tree that has the strongest phylogenetic signal one is simply fitting more accurately to the assumptions of the ancestral state reconstruction method. We showed that this was also true for characters not evolving under Brownian motion (i.e., the characters simulated on uncorrelated branch lengths). Although our simulations were done solely using continuous characters, a direct link between discrete and continuous characters exists in comparative methods (Felsenstein 2005). We thus expect the trends obtained in this study to hold for discrete characters as well. Two common indices to measure phylogenetic signal were used in this study. In both simulations and empirical studies, λ had a stronger relation with the inference accuracy. We believe that, as λ is optimized solely on the internal branches, it may relate more to ancestral state inference in contrast to K, which is measured on the whole tree. This suggests that using λ would be more appropriate when choosing between a chronogram and a phylogram to infer ancestral states. Recently Skinner (2010) showed that taking phenotypic rate heterogeneity into account improved reconstructions of ancestral limb morphology in the scincid lizard clade Lerista. But as the model Skinner proposed uses branch length to infer the probability of changes, choosing the right sets of branches is still of primary importance. Nowadays, it is becoming more common to build large phylogenies to infer patterns of evolution (Smith and Donoghue 2008). These data sets will, without doubt, encompass species with significant genotypic and phenotypic rate heterogeneity, complicating the choice between chronogram and phylogram. In this case, we recommend that inference of ancestral characters should be performed on the tree that shows the strongest phylogenetic signal as it is expected to be more informative and thus gives more accurate reconstructions. This work was supported by the Swiss National Science Foundation (grant numbers 3100AO-116442 and CRS 113-125240). This work received support from the Vital-IT facilities from the Swiss Institute of Bioinformatics. We thank Pascal-Antoine Christin, Socrates Litsios, and two anonymous reviewers for valuable comments on a previous version of the manuscript.