1 School of Biological Sciences, Flinders University, PO Box 2100, Adelaide, South Australia 5001; 2 School of Biological Sciences, University of Adelaide, Adelaide, South Australia 5005; and 3 Earth Sciences Section, South Australian Museum, North Terrace, Adelaide 5000, Australia
Search for other works by this author on: Michael S. Y. Lee Michael S. Y. Lee *1 School of Biological Sciences, Flinders University, PO Box 2100, Adelaide, South Australia 5001; 2 School of Biological Sciences, University of Adelaide, Adelaide, South Australia 5005; and 3 Earth Sciences Section, South Australian Museum, North Terrace, Adelaide 5000, Australia
*Correspondence to be sent to: School of Biological Sciences, University of Adelaide, Adelaide, South Australia 5005; E-mail: Mike.Lee@samuseum.sa.gov.au.
Search for other works by this author on:Associate Editor: Richard Glor
Systematic Biology, Volume 64, Issue 3, May 2015, Pages 532–544, https://doi.org/10.1093/sysbio/syv005
22 January 2015 29 October 2014 Revision requested: 22 December 2014 14 January 2015 22 January 2015Benedict King, Michael S. Y. Lee, Ancestral State Reconstruction, Rate Heterogeneity, and the Evolution of Reptile Viviparity, Systematic Biology, Volume 64, Issue 3, May 2015, Pages 532–544, https://doi.org/10.1093/sysbio/syv005
Navbar Search Filter Mobile Enter search term Search Navbar Search Filter Enter search term SearchVirtually all models for reconstructing ancestral states for discrete characters make the crucial assumption that the trait of interest evolves at a uniform rate across the entire tree. However, this assumption is unlikely to hold in many situations, particularly as ancestral state reconstructions are being performed on increasingly large phylogenies. Here, we show how failure to account for such variable evolutionary rates can cause highly anomalous (and likely incorrect) results, while three methods that accommodate rate variability yield the opposite, more plausible, and more robust reconstructions. The random local clock method, implemented in BEAST, estimates the position and magnitude of rate changes on the tree; split BiSSE estimates separate rate parameters for pre-specified clades; and the hidden rates model partitions each character state into a number of rate categories. Simulations show the inadequacy of traditional models when characters evolve with both asymmetry (different rates of change between states within a character) and heterotachy (different rates of character evolution across different clades). The importance of accounting for rate heterogeneity in ancestral state reconstruction is highlighted empirically with a new analysis of the evolution of viviparity in squamate reptiles, which reveal a predominance of forward (oviparous-viviparous) transitions and very few reversals.
The evolution of viviparity (live-birth, including “ovoviviparity”) from oviparity (egg-laying) has long been a topic of great scientific interest, particularly in squamate reptiles (lizards, snakes, and amphisbaenians), which show much greater lability in reproductive mode than do other amniotes. Studies have focused on a range of aspects including physiology, ecology, and evolutionary history ( Blackburn 2000; Sites et al. 2011; Shine 2014). Viviparity is present at a variety of phylogenetic levels within squamates, including some reproductively bimodal species such as Lerista bougainvillii ( Qualls and Shine 1998) and Saiphos equalis ( Smith et al. 2001). In most taxa, viviparity involves retention of shell-less eggs, with embryos sustained mainly by a large yolk. However, some skinks, including the South American genus Mabuya, and the African genera Trachylepis and Eumecia, have chorioallantoic placentae convergent with those of eutherian mammals, allowing the embryos to be sustained by direct transfer of nutrients across the placenta rather than by the yolk ( Blackburn et al. 1984; Blackburn and Flemming 2009).
One important aspect of the evolution of viviparity is the question of whether oviparity can evolve from viviparity. Oviparity has been assumed to be usually or invariably ancestral to viviparity, given the wide distribution of oviparity in amniotes and the idea that reversals to oviparity are unlikely due to the need to re-evolve complex structures and enzymes involved in the deposition of egg-shell. There have been numerous attempts to reconstruct parity mode evolution in a phylogenetic context. Blackburn (1995, 1985) identified a minimum of 92 oviparity ⇒ viviparity (O ⇒ V) transitions, and noted that approximately 20% of squamate species are viviparous. Also noted was the heterogenous distribution of viviparous groups across the phylogeny, with Scincidae, Iguanidae, Viperidae, and Anguidae containing a large proportion of the inferred origins of viviparity. These studies implicitly used parsimony to identify the minimum number of origins required by then-current taxonomy and phylogeny; they also did not seriously consider the possibility of viviparity ⇒ oviparity (V ⇒ O) “reversals.” Fraipont et al. (1996) suggested that reversals were possible and frequent, and indeed that egg guarding often evolves from viviparity. However, this study was criticized due to methodological problems ( Blackburn 1999; Shine and Lee 1999). A more explicitly phylogenetic study using parsimony to optimize reproductive mode onto a composite phylogeny ( Lee and Shine 1998) identified five possible reversals, but noted that none were well supported, requiring only a small number of extra steps to overturn. Until recently, the best evidence of a V ⇒ O reversal was provided by the oviparous erycine snake Eryx jayakari, which is deeply nested within viviparous species and further lacks the egg tooth typical of oviparous squamates ( Lynch and Wagner 2010).
The recent publication of the most comprehensive phylogenetic analysis of squamates to date (>4000 species: Pyron et al. 2013) allowed a large-scale analysis of squamate parity mode evolution ( Pyron and Burbrink 2014). That analysis employed Binary State Speciation and Extinction (BiSSE) ( Maddison et al. 2007). This model has up to six parameters: independent speciation and extinction rates for taxa possessing each of the two character states, and independent forward and reverse rates between the two character states (i.e., asymmetry). This allows ancestral state reconstruction while accounting for state-dependant diversification and extinction (e.g., viviparity might increase speciation rate). Surprisingly, the optimal model showed strong support for viviparity across many basal squamate nodes, and strong asymmetry, with huge numbers of V ⇒ O “reversals” which outnumbered forward O ⇒ V transitions. This was a particularly surprising and novel conclusion, as it strongly contradicts entrenched ideas ( Pyron and Burbrink 2014). It was acknowledged, however, that parsimony reconstructions yielded a very different and more conventional pattern: later origins and few reversals.
Importantly, the BiSSE analyses—like most model-based comparative studies of discrete characters—assumed that rates of evolution of viviparity were homogenous across squamates. However, analyses of discrete characters are being performed on increasingly large phylogenies, where it is likely that evolutionary rates will vary greatly between different clades ( Beaulieu and O'Meara 2014). Although methods for identifying evolutionary rate shifts in continuous characters have attracted recent attention (e.g., Eastman et al. 2011; Stack et al. 2011; Rabosky 2014), analogous methods for discrete characters have only very recently been developed ( Beaulieu et al. 2013, Beaulieu and O'Meara 2014) and have yet to be widely used. Similarly, there have been few attempts to incorporate among-lineage rate heterogeneity (heterotachy) into ancestral state reconstructions of discrete characters ( Skinner 2010; Beaulieu et al. 2013). The former applied a model that assumed rates were completely uncorrelated even across adjacent branches; our tests with the uncorrelated lognormal clock in BEAST showed that such models essentially mimic parsimony, and cannot identify (for instance) clades with accelerated rates of parity mode evolution. Yet, it is widely recognized that reproductive mode is more labile in some squamate clades than others ( Blackburn 1982, 1985; Shine 1985). Here, we show that accounting for this strong heterotachy in parity mode evolution across squamates gives a very different result to standard BiSSE. Random local clock (RLC) models implemented in BEAST, split BiSSE, and the hidden rates model strongly reinstate oviparity as the ancestral state in basal squamate nodes, and oviparity ⇒ viviparity as the predominant transformation (∼115 origins of viviparity). However, the analyses also robustly support the existence of a small number of reversals (in agreement with Pyron and Burbrink 2014), the best supported of which are Lachesis, E. jayakari, and Liolaemus calchaqui.
Ancestral state reconstruction was first performed in BEAST ( Drummond et al. 2012). These analyses used the time-calibrated tree of 3951 species with reproductive mode data as in the previous study ( Pyron and Burbrink 2014). In the BEAST xml script (Supplementary File S1, available on Dryad http://dx.doi.org/10.5061/dryad.cj2fg), topology and branch lengths were fixed by setting this tree as the starting tree and removing relevant operators, and parity mode was analyzed as an “alignment” of one character. The analyses used a stochastic Mk model ( Lewis 2001), with forward and reverse rates either assumed to be equal/symmetrical (Mk1) or allowed to be unequal/asymmetrical and estimated (Mk2). Parity mode could then be analyzed with clock models normally used in molecular phylogenetics (see below). Four independent runs of each analysis were performed, each of 100 or 200 million generations, with 1000 samples saved. The prior on (mean) rate of evolution was set as an exponential distribution with mean of 0.00126 (the value obtained by dividing the number of changes inferred by parsimony—116—by total branch length). This is deliberately a rather uninformative (flat) prior, which permits very fast rates at relatively high probability; the 95% upper bound entails three times as many changes (rate = 0.00378, ∼348 changes). To further ensure adequate exploration of “faster” rates, the starting value of this parameter in each analysis was 0.1 (i.e., ∼79× the mean prior). Further, increasing the rate prior 10-fold barely changed the retrieved rates. Results were viewed in Tracer ( Rambaut et al. 2014) and convergence of the four runs confirmed by superimposition of likelihood and posterior probability traces, and effective sample size (ESS) of >100 for all parameters. We recorded the total number of forward and reverse transitions, the branch on which each transition occurred in each MCMC sample, as well as the ratio of forward:reverse evolutionary rates, ancestral state for each node, and local clock rate for each branch (see below). The post burn-in samples from the four runs were combined to produce the consensus tree, and associated distributions of the above parameters. Bayes factors (sensu Kass and Raftery 1995) for alternative models were calculated by comparing marginal likelihood estimated using the stepping stone method ( Xie et al. 2011) as implemented in BEAST 1.8.
(1) First, a strict clock was used, thereby enforcing a homogenous rate of trait evolution across all branches; an asymmetric model—permitting different relative rates for forward (O ⇒ V) and reverse (V ⇒ O) transitions—was favored.
(2) The RLC was then used ( Drummond and Suchard 2010) to accommodate heterotachy, thus relaxing the assumption of a homogenous rate across all branches; the number, phylogenetic position, and magnitude of rate changes are estimated using MCMC. All RLC analyses were supported over the strict clock by Bayes factors. The focal analysis had an exponential prior of mean 1 on the number of rate changes (thus favoring few rate changes), and a gamma prior on the distribution of relative rates (shape 0.5 and scale 2; thus favoring smaller rate changes, which were equally likely to be either accelerations or decelerations). This prior (which places most weight on fewer, smaller rate changes) will favor results more similar to a strict clock, and is thus conservative in this context (i.e., testing for the effects of relaxing the strict clock). Three additional analyses with different priors on the number and magnitude of rate changes were also performed, details of which are given as Supplementary File S2 on Dryad. The priors in analysis 2 were deliberately permissive of large numbers of rate changes. One further analysis, with another set of priors, was run for 500 million generations as it appeared to be cycling slowly between two very different likelihood peaks (see Discussion).
The xml file for the focal analysis, including the stepping-stone calculations, is included as Supplementary File S1 (Dryad); the other analyses require simple changes to relevant parameters in this file.
To make our results directly comparable with Pyron and Burbrink (2014), reproductive mode evolution was inferred using the Binary State Speciation and Extinction model ( Maddison et al. 2007), a six-parameter model that estimates state-dependent speciation, extinction, and character transition rates. The data have previously been analyzed using standard BiSSE ( Pyron and Burbrink 2014); we repeated this analysis, confirming the surprising conclusions: an early origin of viviparity and multiple reversals to oviparity (see Results).
We then used split BiSSE ( FitzJohn 2012), which permits particular clades, selected a priori, to have different rates of trait evolution compared with their “background.” Appropriate clades were identified by the position of large rate shifts in the consensus tree from the focal BEAST analysis (see Results), and tested in various combinations. These analyses were performed using the R package Diversitree ( FitzJohn 2012). The split BiSSE models were constrained such that only trait transition rates (O ⇒ V and V ⇒ O) varied across the clades, while state-dependent speciation and extinction rates were homogenous across the whole phylogeny. Alternative root states were weighted according to the relative probabilities of generating the observed data from that root state, given the model parameters and tree ( FitzJohn et al. 2009). The best-fit model was also tested with equal weighting of both root states, but this had no effect on parameter estimation (Supplementary Table S3).
In addition to the standard maximum-likelihood analyses, the best-fit model was also run in a Bayesian framework using the MCMC function in Diversitree, with an exponential prior of rate = (1/(2 * state-independent diversification rate)) on all parameters. Four independent runs of 1500 generations were performed, with convergence confirmed by superimposition of likelihood and posterior probability traces, and ESS of >1000 for all parameters. The post burn-in samples for the four runs were combined for summary statistics.
The hidden rates model ( Beaulieu et al. 2013) is another method that accommodates heterotachy in ancestral state reconstruction. This model divides each character into a number of rate classes (known as “hidden states”). In contrast to the previous two methods, it does not specifically identify clades with different rates, although this is implicit in rate category assignment. We analyzed the same tree and data using this model in the R package corHMM ( Beaulieu et al. 2013). Models with between two and five rate categories for each state were analyzed and compared using AIC.
The BEAST RLC has been tested on molecular data and found to yield highly consistent results ( Drummond and Suchard 2010). However, as this method represents a novel approach to inferring ancestral states for discrete binary characters, simulations were carried out to assess its performance. Trait evolution incorporating asymmetry and/or heterotachy was carried out in R using the phytools, geiger, and ape packages ( Paradis et al. 2004; Harmon et al. 2008; Revell 2012; see Dryad Supplementary File S4). Four sets of five replicate simulations were analyzed, each with simulated trees of 500–600 taxa (tips) and a height of five time units, and an increase in trait evolution rate occurring at a random node that had between 200 and 300 descendent taxa (“fast clade”). The four sets were as follows: (1) As a control, simulations with only asymmetry in transition rates: forward rate 0.2, backward rate 0.01; (2) heterotachy only: 0.005 forward and backward rate, with an increase in the “fast clade” to 0.2; (3) heterotachy and asymmetry: 0.005 forward and 0.00025 backward rate, with an increase in the “fast clade” to 0.2 forward rate and 0.01 backward rate; (4) heterotachy and variable asymmetry, that is, the degree of asymmetry differing in the fast clade compared with the background: 0.005 forward and 0.001 backward rate, with a rate increase to 0.2 forward rate and 0.01 backward rate. The last model is more complex than any available BEAST model, and was done to assess performance when the real data are more complex than the RLC model allows. The simulated character states were exported and analyzed in BEAST using: (1) a strict clock with no asymmetry (Mk1: simplest model); (2) strict clock with asymmetry estimated (Mk2: the “correct” model for simulation set 1); (3) RLC with no asymmetry (RLC+Mk1, “correct” model for simulation set 2); and (4) RLC with asymmetry constant but estimated (RLC+Mk2, “correct” model for simulation set 3, and nearest model for set 4). Mean evolutionary rate had a deliberately flat, uninformative prior (exponential with mean 0.01). Convergence of runs was assessed as before, with two runs of each analysis being performed due to computational demands. To compare BEAST with traditional likelihood methods, ancestral state reconstructions were also performed on the simulated character states using the ace function in APE ( Paradis et al. 2004), with both equal rates (Mk1) and unequal rates (Mk2) models. Note that some simulation replicates gave results that were trivial to reconstruct (e.g., a shift to viviparity near the base of the tree typically resulted in a single large clade of viviparous taxa and a very straightforward optimization under all methods); we thus analyzed five non-trivial replicates for each of the four simulation scenarios. Models were compared using Bayes factors for the BEAST models (with marginal likelihoods calculated using the stepping stone method, Xie et al. 2011), and likelihood ratio tests for the ace models.
The performance of each ancestral state reconstruction was assessed qualitatively. After viewing the results, it was found that each analysis could be clearly placed in one of three categories: (1) those in which there was strong support for the incorrect root state; (2) correct root state but reconstruction inaccurate, in general this involved strong support for the derived state at basal nodes within the fast clade, and incorrect estimates of direction of asymmetry (Supplementary Table S4); (3) reconstruction accurate, with the correct state at basal nodes in the fast clade, and the correct estimates of direction of asymmetry. In Table 3, those reconstructions in the last category are highlighted by shading.
Assuming a strict clock on parity mode evolution—that is, a homogenous rate across the tree resulted in a reconstruction consistent with Pyron and Burbrink (2014). An asymmetrical model with separate forward (O ⇒ V) and reverse (V ⇒ O) rates was favored over a symmetrical model with equal forward and reverse rates (Bayes factor 32.8 sensu Kass and Raftery 1995). In the favored model, V ⇒ O “reversals” evolve at 6.9× the rate of O ⇒ V forward transitions (Supplementary Fig. S1 on Dryad). Accordingly, across the sampled trees, approximately 89.7 V ⇒ O reversals were reconstructed, compared with only approximately 56.7 O ⇒ V transitions. There is an early evolution of viviparity in squamates (above dibamids and gekkotans); however, support for either ancestral state was equivocal in all basal squamate nodes (Supplementary Fig. S2). A possible explanation is that high rates of trait evolution in particular clades elevated the estimated “tree-wide” rate, resulting in uncertainty in reconstructions in deeper nodes (see Discussion).
Implementing a RLC—and thus accommodating heterotachy—further improved model fit (Bayes factor 80.24 over the asymmetrical rate model) and also changed the inferred pattern drastically ( Fig. 1). The focal RLC analysis robustly restores oviparity as the ancestral state in the common ancestor of squamates (pp = 0.9997) and in all basal squamate nodes (pp >0.99). The evolutionary rate of forward (O ⇒ V) transitions greatly exceeds the reverse rate (ratio = 5.56:1), in accordance with conventional views. Accordingly, the absolute numbers of these transitions exhibit a similar ratio. There are an average of 118.73 forward transitions in the sampled trees, with origins of viviparity reconstructed on 115 branches in the consensus tree for the focal analysis (2 in Gekkota, 39 in Scincoidea, 3 in Lacertoidea, 8 in Anguimorpha, 26 in Iguania, and 37 in Serpentes). Note that this does not include origins within reproductively bimodal species, which had to be treated as unknown due to constraints in BEAST. The majority of these 115 branches have an origin with >0.5 posterior probability; however, origins with pp 0.5; see Table 1) reconstructed with reversals: all are in liolaemids and snakes (Booidea and viperids): Liolaemus calchaqui, L. chiliensis, Anomochilus leonardi, Eryx Jayakari, Calabaria reinhardtii and Lachesis. Other reversals with lower support (pp <0.5) are also mostly in these clades (Liolaemus platei, L. reichei, and Trimeresurus borneensis) except for one in Trachylepis (Scincidae). All six well-supported reversals potentially occur within a relatively narrow time “window” after the preceding origin of viviparity (see Discussion).
Reproductive mode in squamate reptiles, as reconstructed by Bayesian methods accommodating heterogeneity in evolutionary rates across different clades (heterotachy). Consensus tree from BEAST random local clock analysis. Branch color represents the most probable (pp >0.5) reproductive mode in the descendant node; the tree reveals a preponderance of forward (oviparity ⇒ viviparity) and few reverse (viparity ⇒ oviparity) changes. A higher-resolution electronic version of this figure, with names and reproductive states for terminal taxa, is available on Dryad ( Fig. 1 expanded).
List of candidate viviparity ⇒ oviparity reversals, based on the consensus tree of the BEAST analyses with the following statistics: number of sampled species, mean posterior probability of the reversal over the four main BEAST analyses, and minimum and maximum inferred duration that each lineage was viviparous before the reversal in the focal analysis (see main text)
| . | Number of species . | mean PP over four BEAST analyses . | Min time viviparous (Ma): focal analysis . | Max time viviparous (Ma): focal analysis . |
|---|---|---|---|---|
| Liolaemus chiliensis | 1 | 0.5961 | 2.5 | 15 |
| Liolaemus calchaqui | 1 | 0.8665 | 1.2 | 7.9 |
| Calabaria reinhardtii | 1 | 0.4789 | 1.9 | 78.1 |
| Eryx jayakari | 1 | 0.9414 | 22.1 | 78.1 |
| Anomochilus leonardi | 1 | 0.5487 | 10.8 | 72.94 |
| Lachesis | 2 | 0.9890 | 8.7 | 26.1 |
| . | Number of species . | mean PP over four BEAST analyses . | Min time viviparous (Ma): focal analysis . | Max time viviparous (Ma): focal analysis . |
|---|---|---|---|---|
| Liolaemus chiliensis | 1 | 0.5961 | 2.5 | 15 |
| Liolaemus calchaqui | 1 | 0.8665 | 1.2 | 7.9 |
| Calabaria reinhardtii | 1 | 0.4789 | 1.9 | 78.1 |
| Eryx jayakari | 1 | 0.9414 | 22.1 | 78.1 |
| Anomochilus leonardi | 1 | 0.5487 | 10.8 | 72.94 |
| Lachesis | 2 | 0.9890 | 8.7 | 26.1 |
Notes: Strongly supported reversals (PP >0.5 in all four analyses) are in bold. The exact parameters of the four BEAST analyses are in Supplementary File S2.
List of candidate viviparity ⇒ oviparity reversals, based on the consensus tree of the BEAST analyses with the following statistics: number of sampled species, mean posterior probability of the reversal over the four main BEAST analyses, and minimum and maximum inferred duration that each lineage was viviparous before the reversal in the focal analysis (see main text)
| . | Number of species . | mean PP over four BEAST analyses . | Min time viviparous (Ma): focal analysis . | Max time viviparous (Ma): focal analysis . |
|---|---|---|---|---|
| Liolaemus chiliensis | 1 | 0.5961 | 2.5 | 15 |
| Liolaemus calchaqui | 1 | 0.8665 | 1.2 | 7.9 |
| Calabaria reinhardtii | 1 | 0.4789 | 1.9 | 78.1 |
| Eryx jayakari | 1 | 0.9414 | 22.1 | 78.1 |
| Anomochilus leonardi | 1 | 0.5487 | 10.8 | 72.94 |
| Lachesis | 2 | 0.9890 | 8.7 | 26.1 |
| . | Number of species . | mean PP over four BEAST analyses . | Min time viviparous (Ma): focal analysis . | Max time viviparous (Ma): focal analysis . |
|---|---|---|---|---|
| Liolaemus chiliensis | 1 | 0.5961 | 2.5 | 15 |
| Liolaemus calchaqui | 1 | 0.8665 | 1.2 | 7.9 |
| Calabaria reinhardtii | 1 | 0.4789 | 1.9 | 78.1 |
| Eryx jayakari | 1 | 0.9414 | 22.1 | 78.1 |
| Anomochilus leonardi | 1 | 0.5487 | 10.8 | 72.94 |
| Lachesis | 2 | 0.9890 | 8.7 | 26.1 |
Notes: Strongly supported reversals (PP >0.5 in all four analyses) are in bold. The exact parameters of the four BEAST analyses are in Supplementary File S2.
The focal RLC analysis converged on a mean of 12.87 shifts in evolutionary rate (95% HPD = 7–17). The consensus tree reveals certain clades have greatly elevated rates of trait evolution ( Fig. 2). The fastest rates are within Sceloporus (∼0.0145, i.e., 1.45% probability of change per lineage per million years) and Liolaemus (∼0.0109) and are an order of magnitude faster than the mean rate across the tree (∼0.00125). Both genera have previously been suggested as being highly labile for reproductive mode ( Benabib et al. 1997; Schulte et al. 2000; Lambert and Wiens 2013; Pincheira-Donoso et al. 2013). Other groups with fast rates are: Scincoidea, Anguimorpha, Alethinophidia (excluding Colubridae sensu lato), and Viperidae.
Among-lineage variability in rates of evolution (heterotachy) of reproductive mode in squamates. Consensus tree from BEAST random local clock analysis colored according to rate; seven major clades/grades with rates different to their “background” are identified. Fastest rates are in a Sceloporus subclade and Liolaemus, and are an order of magnitude faster than the overall rate in squamates. A higher-resolution electronic version of this figure, with names and reproductive states for terminal taxa, is available on Dryad ( Fig. 2 expanded). Exact rates for all branches are presented in Dryad file S3, viewable in Figtree ( Rambaut 2014).
The other three BEAST analyses were implemented with different priors on the number and magnitude of rate shifts. Priors that were more permissive regarding the number of rate changes essentially allowed extra smaller-scale changes, which did not change the overall pattern. The fastest rates were always identified in the same six clades, there was a preponderance of forward transitions over reversals, and all basal nodes in squamates were robustly reconstructed as oviparous. In addition, the individual reversals identified above appeared in each analysis, although only Lachesis, E. jayakari, and L. calchaqui have posterior probability of a reversal >0.5 across all four analyses.
One further BEAST RLC analysis (Dryad Supplementary File S2, analysis 5) suggested the existence of two optimality peaks (Supplementary Fig. S3): each of the four runs spent some time sampling one optimum with fewer (∼4) rate changes and a bias toward V ⇒ O reversals, but three of the four runs eventually favored a slightly higher optimum with more (∼9) rate changes and a bias toward forward O ⇒ V transitions. Because of slow cycling between these optima, this analysis cannot be considered to have converged. However, it is consistent with the results of the other BEAST analyses using strict clock and the random local clock: assuming less (or no) rate heterogeneity in parity mode evolution favors reconstructions with a preponderance of V ⇒ O reversals.
Standard BiSSE reconstructed an early origin of viviparity in squamates and preponderance of V ⇒ O “reversals”, exactly as previously proposed ( Pyron and Burbrink 2014). However, split BiSSE models, which accommodate heterotachy by permitting different (prespecified) clades to have different evolutionary rates for parity mode, gave very different reconstructions.
The tested model with the best AIC support allowed all the following clades to have separate rates of parity mode evolution (compared with their background rate): Scincidae, Anguimorpha, Sceloporus (subclade bracketed by Sceloporus grammicus and S. adleri), Liolaemidae (Phymaturus–Liolaemus subclade), Alethinophidia, “Henophidia” (Booidea of Pyron et al. [2013] plus Bolyeriidae), and Viperidae. Clade-specific evolutionary rates retrieved in this analysis were highly congruent with results from the BEAST RLC analysis. In all partitions the rate of O ⇒ V transitions was higher than the rate of V ⇒ O reversals ( Table 2). The rate of V ⇒ O reversals is extremely low (
Repeating the split BiSSE analyses with alternative combinations of clades revealed two alternative patterns (Fig. 3): models with a greater number of split clades, a bias toward O ⇒ V transitions, and strong support, versus models with fewer split clades, a bias toward V ⇒ O reversals and weaker support. The results in Figure 3 therefore fall into two clusters, with AIC values that do not overlap. These seem to reflect two optimality “peaks” similar to those suggested by one of the BEAST analyses: a higher peak (with a greater number of rate shifts, and a preponderance of O ⇒ V transitions), and a slightly lower peak (with fewer or no rate shifts and a preponderance of V ⇒ O reversals). Table 2 summarizes the results from four key models. Details of the other models used are shown in Supplementary Table S2.
The results for 14 split BiSSE models (different numbers and combinations of clades chosen) as well as standard BiSSE. Here, model support (ΔAIC) is plotted against the bias in transition rate parameters. Transition rate bias is calculated as log(rate of O ⇒ V transition/rate of V ⇒ O reversals) for the “background” partition. The best-supported model has ΔAIC of 0, a greater ΔAIC means lower support. The results fall into two groups: (1) lower right-shaded quadrant—models with a bias toward O ⇒ V forward transitions and higher support and (2) upper left-shaded quadrant—those models which have a bias toward V ⇒ O reversals and lower support. The arrows point to the best-supported model in our analysis (with seven clade-specific rates: Table 2), and to the standard BiSSE model, which assumes a uniform rate across all squamates ( Pyron and Burbrink 2014). Since there is no overlap in ΔAIC of the two quadrants, this provides strong support for oviparity as the ancestral state and O ⇒ V as the predominant transition.
Comparison of model parameters and likelihood for the best-fit split BiSSE model (seven split clades, 20 parameters), another split BiSSE model (5 split clades, 16 parameters), and the standard BiSSE (0 split clades, 6 parameters) used by Pyron and Burbrink (2014)
 |
|
Notes: The predominant direction of change is in bold. *V ⇒ O reversals predominate: it can be seen that a bias toward reversals in Alethinophidia in the model with five split clades can be overturned by further isolating two clades of snakes that exhibit high rates of parity mode evolution (model with seven split clades). The fourth column shows that setting extinction rates and all V ⇒ O rates except Liolaemidae, “Henophidia” and Viperidae to 0 reduces the number of parameters with no effect on model fit.
Comparison of model parameters and likelihood for the best-fit split BiSSE model (seven split clades, 20 parameters), another split BiSSE model (5 split clades, 16 parameters), and the standard BiSSE (0 split clades, 6 parameters) used by Pyron and Burbrink (2014)
|
|
Notes: The predominant direction of change is in bold. *V ⇒ O reversals predominate: it can be seen that a bias toward reversals in Alethinophidia in the model with five split clades can be overturned by further isolating two clades of snakes that exhibit high rates of parity mode evolution (model with seven split clades). The fourth column shows that setting extinction rates and all V ⇒ O rates except Liolaemidae, “Henophidia” and Viperidae to 0 reduces the number of parameters with no effect on model fit.
As in Pyron and Burbrink (2014) there is support for higher speciation rate in viviparous lineages, although there is no evidence for higher extinction rates. Trait-dependent speciation rate parameters are therefore not very sensitive to different ancestral reconstructions and different estimated rates and symmetries of trait evolution. However, we hesitate to draw conclusions from this result due to the problems associated with assuming homogenous diversification dynamics across trees with incomplete taxon sampling ( Rabosky 2010; Quental and Marshall 2010; Stadler and Bokma 2013) and diversification asymmetries driven by factors other than the trait under analysis ( Rabosky and Goldberg 2015). BiSSE furthermore has problems estimating extinction rates and can retrieve unrealistically low rates ( Maddison et al. 2007; Davis et al. 2013) with speciation estimates thus representing a greatly time-averaged approximation of net diversification.
Use of four or five rate categories in the hidden rates model produced a reconstruction broadly similar to the results from the BEAST RLC analysis, with oviparity as the ancestral state for squamates and a preponderance of forward transitions. In contrast, use of only two or three categories led to a reconstruction with a viviparous ancestor and numerous reversals to oviparity. AIC supported the use of extra rate classes, and therefore an oviparous ancestor (AIC = 1149.03, 1009.845, 991.28, and 1000.697 for 2, 3, 4, and 5 rate categories, respectively). Supplementary Figure S5 shows the result for four rate categories. This further supports the conclusion that incorporating heterotachy improves the accuracy of ancestral state reconstructions of squamate parity mode.
Summaries of the performance of each method in the simulations are shown in Table 3 for BEAST and Supplementary Table S5 for ace. Figure 4 shows the simulated character states and BEAST ancestral state reconstructions for one example where the simulation had both asymmetry and heterotachy. Full results for all the simulations are in the Supplementary Information (File S5).
Simulated data and reconstructions under different models, for simulation 11. a) Simulated data with both asymmetry and heterotachy: a bias toward forward (blue/oviparous to red/viviparous) transitions, and an increase in transition rates occurring at the yellow star. b–e) Ancestral state reconstructions under four different BEAST models, which permit (b) neither asymmetry or heterotachy, (c) only asymmetry, (d) only heterotachy, and (e) both asymmetry and heterotachy. Only the last model (e) accurately reconstructs ancestral states. Notably, model (c) reconstructs incorrect states for all basal nodes in the “fast clade” and infers the wrong direction of asymmetry (preponderance of viviparous–oviparous “reversals”). The branch labels refer to the posterior probability of the reconstructed parity mode at a few key nodes. Each number refers to the node at the younger (right-hand) end of the branch.
Performance of four different ancestral state reconstruction models on simulated data with different conditions of asymmetry (forward–reverse bias) and among-lineage rate heterogeneity (heterotachy), showing model support (numbers) and quality of ancestral state reconstruction (shading)
 |
|
Notes: Numbers are negative log marginal likelihoods. Shaded cells are those models that produced an accurate ancestral state reconstruction, as shown by comparison with the simulation. Unshaded cells represent inaccurate reconstructions. Bold denotes the optimal model(s) under Bayes factors: the best-fit model and any model within two log-likelihood units of this. The black boxes surround the models in which the parameters match those of the simulation, or in the case of Simulations 16–20, are the closest match. It is clear that when both asymmetry and heterotachy are present in the data, an accurate ancestral state reconstruction requires both of these effects to be parameterized (i.e., random local clock Mk2).
Performance of four different ancestral state reconstruction models on simulated data with different conditions of asymmetry (forward–reverse bias) and among-lineage rate heterogeneity (heterotachy), showing model support (numbers) and quality of ancestral state reconstruction (shading)
|
|
Notes: Numbers are negative log marginal likelihoods. Shaded cells are those models that produced an accurate ancestral state reconstruction, as shown by comparison with the simulation. Unshaded cells represent inaccurate reconstructions. Bold denotes the optimal model(s) under Bayes factors: the best-fit model and any model within two log-likelihood units of this. The black boxes surround the models in which the parameters match those of the simulation, or in the case of Simulations 16–20, are the closest match. It is clear that when both asymmetry and heterotachy are present in the data, an accurate ancestral state reconstruction requires both of these effects to be parameterized (i.e., random local clock Mk2).
In the case where simulations were performed with only asymmetry in forward/backward rate and no heterotachy, Mk2 models which incorporate this parameter (in ace, BEAST strict clock, BEAST random local clock) outperformed Mk1 models, as expected. However, in one simulation both the BEAST strict clock Mk2 and the BEAST RLC Mk2 models produced incorrect reconstructions, although ace mk2 performed well. All three Mk1 models always assigned the incorrect root state, likely because the high rates of evolution in the simulations led to most tips having the derived state.
Simulations involving only heterotachy and no asymmetry were found to cause significant problems to those models (incorrectly) utilizing the Mk2 model (which estimates asymmetry): ace-Mk2 always reconstructed the incorrect root state, whereas the BEAST strict clock Mk2 model did so in three out of five cases. The presence of heterotachy in the data, but asymmetry in the reconstructing model, means the model tries to explain the rate heterogeneity by incorrectly applying extreme asymmetry. For example, the variability in parity mode in the fast clade—generated by rapid changes in both directions in the simulations—is instead explained by (1) a high rate of reversals relative to forward transitions and (2) reconstructing the derived state at the base of the fast clade, with multiple reversals in this clade. This idea is supported by the estimated asymmetry for the BEAST strict clock Mk2 model in the three cases where reconstructions are inaccurate (0.07, 0.04, 0.08—i.e., strong bias toward reversals) compared with the two cases where they are better (0.44, 1.38—i.e., closer to equality). The ace Mk2 models also estimate very high asymmetry in all five cases (between 0.024 and 0.089). Further evidence is provided in Supplementary Figure S6 and the accompanying discussion. In addition, the branches between the root and the fast clade are often also reconstructed with the derived state in a manner completely at odds with parsimony, but closely resembling the parity mode reconstruction in Pyron and Burbrink (2014). Despite the inaccurate reconstructions for these Mk2 models, the extra parameter was strongly supported statistically in these cases (likelihood ratio tests for ace and Bayes factors for BEAST strict clock). Only the BEAST-RLC-Mk1 model (which most closely matched the simulation model) produced accurate reconstructions for all five simulations. The BEAST-RLC-Mk2 model produced good reconstructions for three out of five simulations, but in one case was subject to the same problem found for the other Mk2 models: high estimated asymmetry and incorrect assignment of the derived state to basal nodes.
For simulations involving both asymmetry and heterotachy, the BEAST RLC Mk2 model produced accurate reconstructions in four out of five cases. For these four, the BEAST RLC Mk2 model is strongly favored over the other BEAST models (Bayes factors). In the case where it performed badly, it appears that the simulation resulted in almost all the tips in the fast clade being in the derived state, causing all the models to be equally inaccurate. In three out of the four cases in which the BEAST RLC Mk2 model is accurate, all the other models perform badly. In one case however, those using the Mk2 model perform equally well, suggesting the simulated heterotachy has had little effect in this instance. Figure 4 shows the reconstructions for one example under these simulation conditions. Figure 4e shows that correct reconstructions are only achieved when both asymmetry and heterotachy are parameterized.
The results for the most complex simulations (set 4)—in which the degree of asymmetry differs between the background and the fast clade—resemble those in which the degree of asymmetry is the same (set 3). This suggests that, at least under certain conditions, the random local clock model can perform well even in cases where the character has evolved in a manner more complex than the model allows. In all five cases, the BEAST RLC Mk2 model produces accurate reconstructions. Bayes factors support the RLC Mk2 model over the other BEAST models in all five cases. The other models are always inaccurate, apart from one case (out of five replicates) in which the BEAST strict clock Mk2 model performs well.
Finally, we note that in all simulations showing heterotachy, the random local clock was able to accurately identify the fast clade. However, the degree of heterotachy was always underestimated, with rates overestimated in the background clade and underestimated in the fast clade (by approximately a factor of 2).
Results from this series of simulations suggest two main conclusions:
The well-supported and highly congruent results from three rather different methods provide strong evidence that oviparity characterized all basal nodes within squamates, and that reversals are infrequent. These methods all allow rates of trait evolution to vary across clades (heterotachy), an assumption that makes intuitive sense and is supported statistically in this particular instance. However, when a homogenous rate is instead assumed across the tree, models that allow asymmetry reconstruct a preponderance of “reversals” ( Pyron and Burbrink 2014). Simulations generate results consistent with the view that the latter conclusions are artifacts of an inappropriate model.
These analyses support the conventional view that viviparity evolves from oviparity and reversals are rare, but also support a small number of reversals (in agreement with Pyron and Burbrink 2014). Such reversals suggest that the loss of structures during the evolution of viviparity may initially involve switching off relevant genetic pathways while still preserving them. Reversals may be possible for a finite time window after the evolution of viviparity, before these pathways have been irreversibly degraded by mutation. As a minimum estimate, the six candidate V ⇒ O reversals occur within a 25 Ma time window after the origin of viviparity, although longer intervals are possible. However even 25 Ma is longer than most estimates for the length of time required to silence a gene through neutral evolution ( Marshall et al. 1994; Lynch and Conery 2000), although these estimates do not factor in the likely existence of pleiotropy for many genes ( Collin and Miglietta 2008).
It should be noted that these reversals are only strongly supported under the assumed tree topology. Uncertainty in the topology of trees assembled from supermatrices with extensive missing data is well-documented ( Sanderson et al. 2011; Simmons and Goloboff 2014). However the highly nested positions of Lachesis and E. jayakari within oviparous groups are also supported by finer scale analyses ( Lynch and Wagner 2010; Fenwick et al. 2012). This is not the case for L. calchaqui, which was not found to represent a reversal in a recent analysis due to differences in topology ( Pincheira-Donoso et al. 2013).
An intriguing observation is the low species diversity in the oviparous clades resulting from each putative reversal, suggesting reversals are not a pathway to diversification. A possible explanation is that reversals are “incomplete” in that not all of the advantages of oviparity can be immediately regained. Such initial constraints may limit the success of lineages that result from a reversion to oviparity, even in the habitats in which oviparity is selectively beneficial. Eryx jayakari, for example, has been proposed to have evolved oviparity as an adaptation to extremely arid deserts ( Lynch and Wagner 2010). Eryx jayakari has not regained the egg tooth ( Lynch and Wagner 2010) and this may constrain other aspects of reproductive biology such as egg-shell thickness, which is very thin in E. jayakari. It is therefore possible that the small number of reversals, and the small clades resulting from them, is due to the interaction of both ecological factors as well as developmental constraints.
Comparative analyses of trait evolution derive greater statistical power from larger phylogenies, permitting more independent contrasts. However, traits are unlikely to evolve at a homogenous rate across large swathes of the tree of life, yet the majority of comparative analyses of discrete traits make this assumption. If a character is evolving rapidly in one part of the tree, but very slowly in other parts, a model assuming a constant rate of evolution is likely to reconstruct an intermediate rate which does not match either true rate, leading to uncertain and/or erroneous reconstructions of ancestral states and forward–reverse ratios. We here demonstrate the importance of incorporating heterotachy into ancestral state reconstruction. Certain squamate clades are highly labile in reproductive mode (e.g., Sceloporus and Liolaemus are an order of magnitude faster than typical squamates), and including this variable dramatically changes inferences of the evolution of reproductive mode. Notably, parsimony analyses effectively permit infinite heterotachy ( Tuffley and Steel 1997): accordingly parsimony reconstructions of this data set also retrieve a pattern of late origins of viviparity and a preponderance of forward changes ( Pyron and Burbrink 2014).
The rate homogeneity assumption may be problematic in the case of parity mode, as discrete states are being applied to what is arguably a continuum. For example, Zootoca vivipara, which is reproductively bimodal, shows thicker eggshells and longer incubation times in Italian and Slovenian oviparous populations compared with French oviparous populations ( Heulin et al. 2002). Clades nearer the boundary of oviparity and viviparity are likely to show higher rates of transitions. Few squamates lay eggs at a developmental stage intermediate between the eggs of typical oviparous species (about stage 30) and the neonates of viviparous species (stage 40) ( Blackburn 1995). However exceptions are found in Sceloporus, oviparous members of which can exhibit facultative egg retention ( Mathies and Andrews 1996; Calderón-Espinosa et al. 2006), consistent with Sceloporus showing the highest rates of parity mode evolution in our analysis. The threshold model ( Revell 2014) is therefore an appealing fit for the evolution of viviparity. However, it was found to be computationally unfeasible on this data set.
Clades that inhabit regions conducive to the evolution of viviparity (e.g., temperate or mountainous regions [ Shine 1985]) would also be expected to have elevated rates of parity mode evolution. In particular, groups that can readily disperse between adjacent high and low altitude regions could be expected to have elevated rates of parity mode evolution. Indeed, this appears to be the case in Liolaemus in the Andes ( Schulte et al. 2000) and Sceloporus in North America (e.g., Mathies and Andrews 1995).
It is here shown that heterotachy poses major challenges for ancestral state reconstruction. This is an important and underexplored area, at least for discrete characters. Researchers performing ancestral state reconstructions on large-scale phylogenies should be vigilant for the warning signs that their analyses have been affected by the artifacts discussed in this article, including radical departures from parsimony and biological orthodoxy, reconstruction of a root state that is scored for the minority of tips, and observable patterns in tip states which indicate strong heterotachy. Particular caution should be exercised when applying models that ignore heterotachy and only allow asymmetry: such models may attempt to explain patterns caused by heterotachy with the only parameter available, asymmetry.
Supplementary material, including data files, can be found in the Dryad data repository (http://dx.doi.org/10.5061/dryad.cj2fg).
This work was supported by the Australian Research Council (to M.S.Y.L.).
The authors thank e-Research SA for grid computing facilities, A. Pyron and F. Burbrink for providing us with a timely copy of their study, R. Fitzjohn and L. Revell for help with R, F. Anderson, R. Glor, S. Lambert, and two anonymous reviewers for comments and suggestions on the manuscript.
