The neutral theory of molecular evolution predicts that rates of phenotypic change are largely independent from genotypic change. A recent study by Bromham et al. (2002) confirmed this expectation, finding no evidence for correlated phenotypic and molecular evolutionary rates in animals. We reevaluate this hypothesis, sampling at different taxonomic levels in plants and animals, using Bayesian inference to reconstruct phylogenetic trees and estimate rates of molecular evolution. We use independent contrasts in branch lengths to maximize the information extracted from each of the trees and nodal posterior probabilities to assess the influence of phylogenetic error. Our results indicate that in vascular plants between 2% and 11% of the variation in phenotypic rates of change can be explained by the rate of genotypic change. These results may be explained by the idea that processes that affect general evolutionary rates, such as body size, may also be expected to influence rates of morphological change.
In addition to species numbers, two commonly used biodiversity measures are phenotypic and genetic (genotypic) variation (e.g., Mallet 1996; Williams and Humphries 1996; Owens and Bennet 2000; Purvis and Hector 2000). However, how the two measures relate to each other has been much debated. The neutral theory (Kimura 1983) predicts that the vast majority of DNA evolution is hidden from selection and is thus decoupled from morphological change. Although it is possible to trace the genetic basis for some morphological traits, the summed genetic component of such changes involve only a fraction of the genome. For example, observations of accelerated rates of molecular evolution during the rapid morphological diversification in the adaptive radiation of the Hawaiian silversword alliance were found to be largely restricted to nonsynonymous sites within a subset of genes regulating developmental processes (Barrier et al. 2001). There was no evidence for an associated acceleration of neutral mutation rates, even within the same subset of genes (Barrier et al. 2001). Because of evidence such as this, suggesting phenotypic change is controlled by few mutations in regulatory genes, for example, Bromham et al. (2002) argued that any correlation between rates of genetic and morphological change cannot be a direct consequence of adaptive change, but instead would have to be a product of a more general genomewide process.
Bromham et al. (2002) proposed two potential mechanisms that might result in a correlation between rates of molecular and morphological evolution. First, a higher substitution rate may increase the frequency of phenotypically advantageous mutations, thereby influencing the pace of morphological change. Second, an independent variable, such as population size, might influence genotypic and phenotypic rates separately. For example, if rapidly speciating lineages were characterized by small population size, we might expect accelerated rates of fixation for alleles with weak selective effect (Bromham 2003). Rapidly diversifying lineages might have smaller population sizes as a product of repeated population subdivision; alternatively, the increased rate of fixation of nearly neutral alleles may lead to hybrid incompatibility (Orr 1995) and thereby drive the process of lineage divergence.
To date, few studies have demonstrated a significant correlation between genetic and morphological rates of change (Barrier et al. 2001; but see Barraclough and Savolainen 2001), supporting expectations. Omland (1997), however, suggested that the frequency of published studies demonstrating a decoupling of genetic and phenotypic change were biased in that the taxa selected for such studies were done so on the a priori expectation that they would exhibit different rates of change. In an analysis of eight clades, encompassing a diverse array of taxa from Anas (dabbling ducks) to Sedum (a genus of succulent flowering plants), Omland found that the number of reconstructed morphological change along the branches of the phylogenetic tree (morphological branch lengths) were significantly correlated with the number of reconstructed changes in genetic characters (molecular branch lengths). However, Bromham et al. (2002) questioned whether Omland's statistical approach corrected sufficiently for the phylogenetic nonindependence of taxa needed to account for the confounding effect of evolutionary time and shared phylogenetic history.
Bromham and coworkers readdressed the relationships among morphological and molecular rates in 13 new vertebrate datasets plus three of the datasets included within Omland's original study. Three statistical approaches were adopted, namely PARE, APPLE, and GRAPE (Bromham et al. 2002). PARE estimated rates of change by dividing the number of reconstructed changes per branch by its estimated time-span, APPLE contrasted the average branch lengths for the largest possible sister group contrasts (those containing the largest number of terminal taxa), and GRAPE compared the relative branch lengths from the taxon pair with the largest tip-to-tip branch length. No significant association was found between molecular and morphological rates of change, including the reanalysis of the data within Omland's (1997) original study.
Although we certainly agree with the need to ensure rigorous nonindependence among data within any statistical analysis, we question the utility of the three tests outlined above. First, PARE uses nonparametric rate smoothing (NPRS; Sanderson 1997) to estimate relative branching times from the molecular dataset. NPRS attempts to correct for rate variation across lineages by assuming rate autocorrelation, an assumption that may be violated among datasets (Sanderson 2002). Second, APPLE and GRAPE limit the number of datapoints included within the analysis, for the most part, to one per tree, limiting the statistical power to detect a significant relationship.
Here we used independent contrasts (Felsenstien 1985) to investigate the relationship between rates of molecular and morphological change. We sampled data at several taxonomic levels in plants and animals, although we focused on vascular plants. We used Bayesian inference (Huelsenbeck et al. 2001) to reconstruct phylogenetic trees from DNA sequence data and estimate posterior probabilities supporting the relationships among clades within each tree. Our approach attempts to maximize the amount of information that may be obtained from each tree topology, while maintaining statistical independence.
Methods
Phylogenetic Reconstructions
We collated 13 phylogenetic datasets for vascular plants and one for weevils from the recent literature in which both molecular (from different sets of genes and DNA regions) and discrete morphological characters have been scored. In addition, we included a subset of the data (i.e., those that were available) originally analyzed by Omland (1997; beetles, echinoids, and salamanders) and Bromham et al. (2002; whales, crocodilians, and lizards). A list of datasets, source publications, and respective sample sizes is given in Table 1. After excluding taxa without DNA sequence data, phylogenetic trees were reconstructed from the molecular data using likelihood-based Bayesian analyses in MrBayes (ver. 2.01; Huelsenbeck and Ronquist 2001), and assuming the HKY85 (Hasegawa et al. 1985) model of DNA evolution. More complex evolutionary models tended to provide a better fit to the data (Appendix Table 1); however, for consistency within our analysis we used the HKY85 model because it provides a reasonable compromise between computational tractability (especially given the large size of many of our matrices) and generality (Soltis et al. 2002). While parameter-rich models tend to score higher likelihood values, they do not necessarily lead to increased phylogenetic accuracy (Piontkivska 2004). To evaluate the likely influence of model choice on branch lengths, we repeated the Bayesian analysis of the Nandi et al. (1998) dataset using the more complex model of DNA evolution that was estimated to provide the best fit to the data (see Appendix Table 1). If alternate models produce widely discrepant branch lengths, it may be most apparent within this dataset, because it includes the greatest number of taxa and hence the greatest number of branches in the phylogenetic tree.
For the Bayesian searches we sampled every 10th generation across four chains (three heated, one cold) for 5 × 106 generations, saving trees with branch lengths. Log likelihoods were plotted against generation to check for stationarity, and the first 1 × 106 generations were discarded as burn-in. A summary tree with clade posterior probabilities and a phylogram with mean branch lengths was output using the “sumt” command in MrBayes. All trees are available from the authors. Finally, morphological branch lengths were estimated from the Bayesian topologies in PAUP (ver. 4.0b10; Swofford 2003) using parsimony for both ACCTRAN and DELTRAN optimizations.
Branch Length Contrasts
Statistical analyses
For every node in the phylogenetic trees, we calculated node-to-tip distances for both the Bayesian molecular and morphological phylograms. For nodes deeper in the tree, this distance was estimated by successively averaging the branch lengths from the tips of the tree (see Barraclough et al. 1996). We then calculated the difference in branch lengths between the sister-clades descending from each node. Because sister clades are by definition the same age, differences in branch lengths directly reflect differences in evolutionary rates of change. We performed regression through the origin (Harvey and Pagel 1991) on the branch length contrasts in the statistical package R (ver. 1.90; freely available at http://www.r-project.org/). To evaluate the generality of our findings, we analyzed each phylogenetic tree separately and then combined contrasts within both animals and plants.
Because nodes deeper in the tree will by definition subtend a greater number of terminals, they are vulnerable to the bias introduced by the node density effect (branches broken by numerous nodes tend to recover more changes and thereby be, on average, longer than unbroken branches; Fitch and Beintema 1990). Therefore, we repeated our analysis including only nodes subtending terminals (here referred to as “terminal sisters”). However, characters in morphological datasets might have been intentionally selected to resolve relationships between taxa and thence underscore autapomorphies; consequently, terminal branch lengths may be underestimated (Bromham et al. 2002). Therefore, we randomly deleted one taxon from each sister pair and reestimated terminal branch lengths for the new sister pairs (here referred to as “pruned sisters”). Molecular branch lengths were recalculated in PAUP using the mean model parameter values derived from the Bayesian analysis after the burn-in had been discarded, and the regression analysis was repeated on the new sets of terminal contrasts.
Sensitivity analysis
We employed three approaches to assess the sensitivity of our results to violations of the assumptions of the regression model and to phylogenetic error. First, to maximize homogeneity of variance, we standardized the contrasts in both morphological and molecular branch lengths by dividing them by the average of the branch lengths leading to the respective sister pairs. Second, we used the nodal Bayesian posterior probabilities as weights in the regression analysis, thereby upweighting contrasts with greater phylogenetic support. Last, we reevaluated significant relationships using the nonparametric Spearman's rank correlation coefficient, as it removes the assumption that the data follow a normal distribution.
Estimating Type I error
To explore potential bias, we simulated branch length contrasts using the Nandi et al. (1998) phylogenetic dataset for angiosperms, our largest and most comprehensive dataset. We used a randomization procedure to shuffle morphological character states among the taxa, while keeping both the number of morphological traits and the overall frequency of character states constant. Morphological branch lengths were then calculated in PAUP with ACCTRAN from the Bayesian tree topology. We performed 1000 randomizations, producing 1000 phylograms with differing morphological branch lengths but identical branching pattern. Our null expectation is of no relationship between morphological and molecular branch lengths. A significant correlation between morphological branch lengths derived from the randomization procedure and the Bayesian molecular branch lengths would indicate Type I error.
Results
All Bayesian searches were found to have reached stationarity in log likelihoods prior to the burn-in, indicating that run length (number of generations) was adequate for sufficient sampling from the posterior probability distribution of trees (Ronquist and Huelsenbeck 2005). Comparing raw contrasts in branch lengths for matching terminal sisters, derived from analysis of the Nandi et al. (1998) dataset using alternative models of DNA evolution, indicated that model choice had little impact on the magnitude of the contrasts (r2 = 0.99; regression between the contrasts in molecular branch lengths estimated from the HKY85 and best-fit model of DNA evolution).
Our randomization procedure revealed a strong tendency for contrasts in morphological branch lengths to correlate with contrasts in molecular branch lengths when all nodes were included within the analysis. Of 1000 replicates, 999 (979 for the standardized contrasts) demonstrated a significant correlation between the branch length contrasts. However, when only terminal sister contrasts were considered in the regression model, just 54 (51 for the standardized contrasts) replicates demonstrated a significant correlation. By shuffling character states randomly among the taxa, we removed phylogenetic signal from the morphological characters, hence the number of significant results should be directly proportional to alpha (around 50 at alpha = 0.05), and this is what was found for terminal sister contrasts. The randomization results show that the analysis including all nodes has an unacceptably high Type I error rate (>95%), and we therefore report only results for terminal sister contrasts in the following section.
When each tree was analyzed separately, we found no consistent evidence for a significant correlation between contrasts in morphological and molecular branch lengths for either the raw or standardized contrasts, regardless of taxonomic group, optimization algorithm (ACCTRAN versus DELTRAN), and weighting (Tables 2, 3). However, the combined terminal sisters were found to exhibit a highly significant but weak correlation between the raw contrasts for plants (P < 0.01, r2 = 0.04; P < 0.01, r2 = 0.11, for the ACCTRAN and DELTRAN optimizations, respectively; all regression models for the combined sister contrasts are given in Table 4). Standardized contrasts were also highly significant for plants, but the strength of the correlation tended to be marginally less. For plants, weighting the contrasts by Bayesian posterior probabilities resulted in a slight decrease in the explanatory power of the regression models for the raw contrasts, but an increase in explanatory power for the standardized contrasts. Within plants the DELTRAN optimization of morphological branch lengths tended to demonstrate the stronger correlation; however, this pattern was reversed for the pruned sister contrasts. In general, the correlation between the branch length contrasts were weaker for the pruned sisters and nonsignificant when the contrasts were standardized for the DELTRAN optimization. The regressions models for the combined sister contrasts were not significant for the animal datasets (Table 4).
The Spearman's rank test supported the overall significant correlation between morphological and molecular branch lengths for both the raw and standardized contrasts within plants (Table 5). However, the significance and strength of the relationship depended upon the branch length optimization algorithm. As for the regression models, the DELTRAN optimization demonstrated the stronger relationship for the terminal sister contrasts (Spearman rank r = 0.14, P = 0.03; r = 0.19 P < 0.01, DELTRAN optimization for the raw and standardized contrasts, respectively), and the ACCTRAN optimization was the stronger for the pruned sisters (Spearman rank r = 0.56, P < 0.01; r = 0.51, P < 0.01, ACCTRAN optimization for the raw and standardized contrasts, respectively; Table 5).
Discussion
Sister-group comparisons offer a statistically rigorous approach for comparing evolutionary rates (Barraclough et al. 1998). We find that molecular rates of change can explain around 2% to 11% of the variation in morphological rates of change in plants, depending upon the set of sister taxa and alternative treatments of the underlying contrasts, but only when datasets were combined. However, we find no evidence for a similar correlation among the animal taxa included within our analyses.
Our results depend upon the accuracy with which state changes can be mapped upon the branches of the phylogenetic trees and the correct identification of sister pairs. We assume the topologies derived from the Bayesian analysis of the molecular data provide the best estimates of phylogeny. By using the mean molecular branch lengths (for the terminal sisters contrasts) from the post-burn-in posterior distribution of trees, we reduce the likelihood of a single incorrect tree topology distorting branch length estimates. Varying the evolutionary model of DNA substitution had little appreciable effect on branch length contrasts. Furthermore, comparisons between the weighted and unweighted regression models indicated that our results were relatively insensitive to phylogenetic uncertainty.
It is possible the Bayesian consensus trees do not provide the most parsimonious reconstructions for the morphological data, or that parsimony is an inappropriate model for inferring the evolutionary history of morphological characters. However, mapping of characters on independently derived phylogenetic trees is commonplace (see Pagel 1999), and the assumption that the tree represents the true tree is implicit within the majority of comparative analyses using independent contrasts (Felsenstein 1985; Harvey and Pagel 1991). If reversals and/or parallelisms are frequent, parsimony will tend to underestimate branch lengths. However, a practical alternative is presently unavailable for mapping multiple characters across large phylogenetic trees, although Markov chain Monte Carlo methods in conjunction with Bayesian inference may offer viable alternatives in the future (Pagel et al. 2004; Ronquist 2004). We believe any error in branch lengths would tend to result in loss of signal when only terminal branches are compared; hence, our analysis is conservative.
Although we found a correlation between plant rates of genotypic and phenotypic evolution when datasets were combined, this relationship was not apparent within the individual clades, consistent with the previous findings of Bromham et al. (2002). We suggest that the apparent conflict with the significant result for the combined sister contrasts in plants may partly reflect lack of statistical power in the separate clade analyses. A post hoc power analysis revealed that a sample size between 71 and 399 would be required to detect the weak correlation we observe among the combined contrasts with 80% probability (see Appendix Fig. 1). The largest single dataset within our analysis, the Nandi et al. (1998) phylogenetic tree for angiosperms, contributes just 49 contrasts; the probability of detecting a significant relationship from these data is between 16% and 64%. By extension, it is possible that the nonsignificant findings within animals and the previous analysis of Bromham and coworkers might also be a product of low statistical power (most animal datasets were significantly smaller than the plant datasets used here; see also Omland 1997). As larger phylogenetic trees become more available through initiatives such as the National Science Foundation–funded Assembling the Tree of Life consortium, it will be possible to examine this hypothesis further.
It is widely recognized that rates of molecular evolution tend to vary significantly among genes, lineages, and possibly over geological time (Soltis et al. 2002; Bromham 2003). Numerous interacting processes are likely to influence both phenotypic and genotypic diversity, for instance, generation time (Bromham et al. 1996), population size (Gillespie 1999), DNA repair, and environment (Gaut et al. 1996; Pawlowski et al. 1997; Davies et al. 2004). There are also many reasons why we might not expect an exact correspondence between phenotypic change and genetic change, even for DNA regions known to code for the trait of interest, for example, due to pleiotropic or epistatic effects. It is therefore not surprising that it has proven difficult to establish a strong link between the two. Furthermore, lineage-specific factors are likely to distort the relationship between rates of change so that generalities may not be possible at finer taxonomic scales.
It has been suggested that elevated rates of molecular evolution associated with adaptive radiations could result in overestimates of divergence times from molecular data (e.g., Lee 1998; Bromham 2003). The more species-rich families of flowering plants were found to have had faster rates of molecular evolution (Barraclough et al. 1996; Savolainen and Goudet 1998; Barraclough and Savolainen 2001), but there was no significant association between morphological rates and species numbers. Davies et al. (2004) proposed that the correlation between molecular rates and species richness might have been an artifact of both covarying with environment independently. We therefore argue that the concerns over dating adaptive radiations are largely unfounded; in any case, our results indicate any bias would be small in comparison to other sources of error within current dating techniques (Sanderson and Doyle 2001; Graur and Martin 2004).
One of the most striking aspects of our study is the fact that the genotype-phenotype correlation does not hold for animals. As mentioned above, it is possible that the correlation was simply more difficult to detect in the animal datasets due to their relatively small sizes (four of the seven animal clades sampled 12 species or fewer). It is also possible that the nature of the morphological and/or molecular characters coded in these matrices is different in the plant versus animal datasets. For instance, polygenic traits may be more common among the plant morphological characters, which may in turn be more prone to change due to background mutations, as any one of multiple genes might be affected. In contrast, animal studies may have traditionally focussed on characters that are more strictly controlled by one or a few genes, and therefore largely independently from the background mutation rate. Recent advances in our understanding of phenotypic evolution and development (see Raff 2000; Baguñà and Garcia-Fernàndez 2003) may provide the details required to explore these questions more thoroughly, once suitable data is available.
Of the molecular characters, most genes used in the animal datasets are from the ribosomal DNA (12S, 16S, 28S) or the mitochondrial genome (cytochrone b, cytochrone oxidase I), whereas in the plant datasets the DNA regions are mostly from plastid DNA or noncoding ribosomal DNA (internal transcribed spacers). Therefore the different results we obtained from plants versus animals could be due to functional constrains imposed on coding ribosomal DNA and/or the animal mitochondrion that affect substitution patterns and prohibit the detection of any linkage between molecular and morphological rates of change. By contrast, the plastid genome, which is less constrained (Savolainen et al. 2002), would provide a better record of background mutation rates and therefore more likely to correlate with morphological rate of change.
Alternatively, it remains possible that plants do truly differ from animals with regard to their relationships between rates of molecular and morphological evolution. This could be due to various reasons, for example, if morphology-controlling genes in plants are more widely spread throughout the genome and therefore more susceptible to background mutation. In animals, mutations that affect morphologies could also be either better corrected by efficient DNA repair systems or more lethal. These hypotheses are still speculative; nonetheless, our findings of a significant association between genotypic and phenotypic rates of change in plants remains.
Our results are consistent with a broadly positive relationship between genomewide substitution rates and morphological rates of change in vascular plants (which can explain 2– 11% of the variance in rates). Whether this relationship may be casual or a result of both being correlated with a third and independent variable, such as body size (Martin and Palumbi 1993; Gillooly et al. 2005), remains to be ascertained. However, we suggest evidence for covariation with rates of cladogenesis is weak, in agreement with Barraclough and Savolainen (2001). If the processes driving morphological change are sensitive to general evolutionary rates, the relationship between rates of change should be strongest for neutral substitutions (also see Omland 1997). However, at present, in both plant and animal phylogenetic datasets, molecular sampling was not sufficient to examine thoroughly differences in the patterns of neutral versus nonneutral molecular change. We await future studies that evaluate this hypothesis from the neutral theory more directly and, in particular, examine potential variation in the mode of evolutionary change between plants and animals.
Acknowledgments
We are very thankful to the authors listed in Table 1; L. Bromham and T. Givnish for kindly sending us their matrices for these analyses; and C. Barr, R. Grenyer, T. Barraclough, P. Soltis, and two anonymous reviewers for their comments and help. This work was partially funded by a grant from the Leverhulme Trust to VS.
Literature Cited
Appendix Fig. 1. Post hoc power analysis illustrating the relationship between statistical power and sample size for the upper and lower limits of the correlation coefficients (r) observed among the combined datasets in plants (alpha = 0.05, null expectation r = 0). Filled circles r = 0.14 (r2 = 0.02), open circles r = 0.33 (r2 = 0.11)
Table 1.
Phylogenetic datasets used in the analyses, listing number of taxa, number of characters, DNA regions, and source.
Table 2.
Results of the regression analysis on the raw contrasts, using both ACCTRAN (ACC) and DELTRAN (DEL) morphological branch lengths (+, positive relationship; −, negative relationship; NA, insufficient contrast to perform regression analyses; NS, not significant; *P < 0.05, **P < 0.01)
Table 3.
Results of the regression analysis on the standardized contrasts, using both ACCTRAN (ACC) and DELTRAN (DEL) mor phological branch lengths (notations as for Table 2)
Table 4.
Regression models for the combined set of raw and standardized contrasts in animals and plants, using both ACCTRAN (ACC) and DELTRAN (DEL) morphological branch lengths (w, weighted regression; ps, pruned sisters; NA, not applicable; NS, not significant; *P < 0.05, **P < 0.01)
Table 5.
Spearman's rank correlation coefficient for the combined set of contrasts in plants (notations as for Table 4)
Appendix Table 1. Best-fit model of DNA evolution for each clade, based on Bayesian information criterion within a decision theory framework, as implemented in DT-ModSel (Minin et al. 2003). This model selection procedure is most appropriate for our purposes, as it incorporates relative branch length error as a model performance measure. See Posada and Crondall (1998) for descrip tions of substitution models
