Rates of morphological evolution in Captorhinidae: an adaptive radiation of Permian herbivores

Phylogeny and time calibration

The phylogeny used was that which was presented in Leibrecht et al. (in press), currently the most comprehensive cladistic analysis of captorhinids. The phylogeny was time calibrated in the R 3.1.2 (R Core Team, 2014) using the method proposed by Lloyd et al. (2016), itself an expansion of a method put forward by Hedman (2010). The method of Hedman (2010) was intended to infer confidence intervals on the age of a specific node in the tree. It is a Bayesian approach using the ages of successive straigraphically consistent outgroup taxa relative to the age of the node of interest to make inferences about the quality of sampling; large gaps between the age of the node of interest and that of the outgroups implies a poorly sampled fossil record, and therefore the age of the node of interest may be inferred to be older. Lloyd et al. (2016) designed a procedure whereby this approach could date an entire tree rather than just a specific node. In applying this method, successive outgroups are required to the total clade. The outgroups to Captorhinidae employed were: Paleothyris and Hylonomus (found to be the outgroups to Captorhinidae in the Bayesian analyses of Müller & Reisz (2005a) and Müller & Reisz (2005b)), Archaeothyris (the earliest known synapsid (Reisz, 1972)), and Westlothiana (a reptiliamorph outside the amniote crown according to Ruta & Coates, 2007). A maximum age constraint was set as 334.7 million years ago, the oldest reliable estimate using molecular dating for the origin of Amniota published within the last five years at the moment of data collection (Parfrey et al., 2011).

Uncertainty surrounding the ages of taxa was accounted for using the method of Pol & Norell (2006). For each taxon (including the outgroups), 100 first appearances and last appearances were drawn at random from a uniform probability distribution covering the full possible range of ages for that taxon. One hundred time-calibrated trees were produced from the 100 sets of ages.

Since the analysis of Leibrecht et al. (in press) produced two most parsimonious trees (MPTs), half the time calibrated trees were based on the first, and half on the second. All analyses described below were carried out on all 100 of these trees. All 100 of these trees are available in Data S1, and the age ranges allowed for each taxon in Data S2.

Reconstruction of dietary evolution

A dietary character with three states, carnivore, herbivore and omnivore, was scored for all taxa present in the phylogeny. Ancestral character states were deduced using conditional (joint) likelihood, employing the ace function in the ape package (Paradis, Claude & Strimmer, 2004) in R. This function allow three models of discrete character change: an equal rates model (transitions between all states in all directions are equally probable), a symmetrical model (transitions between two character states occur with equal probability in either direction, but different pairs of character states have different probabilities of transition) and an all-rates-different model (each transition has a different probability). In order to deduce which model was best for the data available, these models were fit to the captorhinid phylogeny using likelihood methods, employing the fitDiscrete function in the R package geiger (Harmon et al., 2008). The Akaike weights were used to deduce the best fitting model.

Since this likelihood ancestral state reconstruction produced uncertain results surrounding the character state at three nodes (see below), an alternative Bayesian approach for optimizing the dietary character at these nodes was used as an independent test in the program BayesTraits V2 (Barker, Meade & Pagel, 2007). A reversible jump Markov chain Monte Carlo (RJ MCMC) approach was used to reconstruct the ancestral states of the discrete character, allowing different combinations of character states to have different transition rates. The analysis was run for 10,000,000 iterations with 100,000 discarded as burn-in, sampling every 1,000 generations. These results are presented in the Data S1, and the analyses described below are based on the likelihood results.

Analysis of rate variation

Analysis of rate variation was carried out using the method of Lloyd, Wang & Brusatte (2012), later refined by Brusatte et al. (2014) and Close et al. (2015). Discrete morphological character scores may be taken from the matrices used in cladistic analyses, and ancestral states are deduced using likelihood. This allows the number of character changes along each branch to be counted, and rates of character change are calculated by dividing the number of changes along a branch by the branch length. The absolute value calculated for the rate of each branch, however, can be misleading due to the presence of missing data (Lloyd, Wang & Brusatte, 2012). As such it is more useful to identify branches and clades where the rates of character change are significantly higher or lower than others, rather than comparing the raw numbers. This is assessed by comparing two models using a likelihood ratio test, one where the rates of change are uniform across the whole tree and one where the branch of interest has a different rate to the rest of the tree. A similar method is used to compare rates of evolution through time and identify bins where rates of evolution are significantly high or low.

The character data used was from the matrix of Leibrecht et al. (in press). The time bins used to examine rate variation through time were substages, dividing the international stages into two bins: early and late. The analysis was carried out in R using functions from the package Claddis (Lloyd, 2016) on all 100 of the time calibrated trees. The data matrix is presented in Data S3.

Due to the uncertainty surrounding the optimisation of the dietary character, a stochastic mapping approach was used to examine rate heterogeneity in the different dietary classes. For each of the 100 time calibrated trees, the dietary character containing three states (carnivore, omnivore and herbivore) was mapped onto the tree using likelihood. Using the character state probabilities identified at each node, 100 possible evolutionary histories of diet in that tree were generated for each of the 100 phylogenies following the procedure outlined by Bollback (2006), giving a total of 10,000 stochastic maps. The mean rate of herbivorous branches, carnivorous branches and omnivorous branches were calculated in each stochastic map, along with the mean rate of a randomly selected set of branches with a sample size equal to the number of herbivorous branches in that map.

Taxonomic jack-knifing was used as a sensitivity analysis to examine the robusticity of the results. In each iteration, three randomly selected taxa (the maximum number that could be deleted leaving at least one member of each dietary regime) were dropped prior to running the analyses of rate variation described above.

Disparity

The character matrix of Leibrecht et al. (in press) was also used to examine morphological diversity (disparity). Morphological distances between taxa were calculated using the Maximum Observable Rescaled Distance (MORD) distance measure of Lloyd (2016), which was shown to perform better in datasets with large amounts of missing data. Following the suggestion of Brusatte et al. (2011), the internal nodes of the phylogeny were treated as data points, with their character scores inferred using ancestral state reconstruction, in order to account for the incomplete sampling of the fossil record; these data points represent ancestral taxa that may have possessed character combinations not observed in sampled taxa.

Having generated a distance matrix, once again the stochastic mapping approach was used to compare disparity in different dietary classes. For each of the 10,000 evolutionary histories generated, each taxon (both tip and node) was assigned a dietary class, and the mean MORD distance for each of the three dietary classes was calculated.

Disparity through time was investigated by subjecting the MORD distance matrices to a principal coordinate analysis. Disparity in each time bin was calculated as the sum of variances of the PC scores of each taxon in that bin. An attempt was also made to incorporate both ghost lineages and internal branches into the analysis using a novel method illustrated in Fig. 1. Taxon A is present in time bin 3, and its ancestral node is inferred to be in time bin 1. Therefore there must be a ghost lineage present in time bin 2 (Fig. 1A), which would be ignored in the disparity analysis under the method of Brusatte et al. (2011), wherein only node and tip morphologies were included. The morphology inferred in time bin 2 will depend on which model of evolution is preferred; under a gradualistic model of evolution, assuming no change in rate along the branch (Fig. 1B), the principal coordinate score in time bin 2 may be inferred by calculating the rate of change in the principal coordinate along that branch, and the amount of time between the ancestral node and the midpoint of time bin 2. Alternatively one may assume a punctuated model of evolution, where the morphological change occurs rapidly at the time of speciation in time bin 1, and the lineage experiences morphological stasis for the remaining time; thus the PC scores inferred in time bin 2 will be identical to that of the tip in time bin 3 (Fig. 1C). Both methods are used here to compare the results. Again, the stochastic mapping approach was used to assign a diet to each branch, allowing the comparison of patterns of disparity through time in each of the dietary regimes. The function used to infer trait disparity through time while including ghost lineages is presented in Data S5.

Lineage density

Sidlauskas (2008) highlighted that rate heterogeneity is not the only means by which different clades may have different disparities. The “efficiency” with which the taxa explore morphospace will also have an influence. A clade which continually returns to the same region of morphospace will exhibit a lower disparity than a clade that explores novel regions of morphospace, even if the rates of morphological change do not vary (Sidlauskas, 2008). In fact, it has been observed that increased disparity can occur even alongside decreases in evolutionary rates; when examining the evolution of mammals across the Cretaceous/Paleogene boundary, Slater (2013) demonstrated that, although rates of body size evolution were lower after the boundary, an increase in variance was possible due to the release of evolutionary constraint.

Sidlauskas (2008) introduced a method to examine this concept: Lineage Density. This is a measure of the amount of evolutionary change within a clade relative to the area of morphospace explored. It is calculated by dividing the sum of the morphometric branch lengths within a clade by the volume of the 95% confidence hyperellipsoid or convex hull of the clade’s morphospace. A lower lineage density indicates that the taxa are exploring novel regions of morphospace rather than continuously returning to the same morphologies.

The principal coordinate analysis described above was used as the basis for calculating lineage density. For each of the 10,000 dietary histories generated by stochastic mapping, the branch lengths within each dietary regime were calculated as the Euclidean distance between each data point i.e., the morphological distance travelled by that branch. The 95% confidence hyperellipsoid volume of each dietary regime was calculated using the functions in the R package Cluster Maechler et al. (2008). The lineage density of each dietary regime was calculated using the equation provided by Sidlauskas (2008).

Character change histories

The character list of Leibrecht et al. (in press) was divided into five categories based on the region to which the characters referred to: Skull, Palate, Mandible, Dentition and Postcranium. The functions in the package Claddis automatically calculates the most likely combination of character changes for each of the 100 time calibrated phylogenies alongside the analysis of rate variation. These character change histories were used to assess which region of the skeleton underwent the greatest change within each dietary regime. Using the 10,000 stochastic maps of the dietary character, the number of characters from each region changing within each dietary regime was counted. These counts for each region were divided by the total number of character changes occurring across the entire tree in that region to account for the fact that the characters were not evenly distributed. The list of characters and the region to which they were assigned is presented in Data S4.

Dietary evolution

Fitting of discrete models of character evolution to the dietary character indicates an equal rates model best fits the captorhinid phylogeny (Fig. 2). It should be noted that this support is not overwhelming; although the ER model is found to fit best in all 100 trees, in none does it receive an Akaike weights score above 0.8. Using this model in ancestral state reconstructions (Fig. 3) indicates that a single transition to a herbivorous diet is most probable. Labidosaurus, judged to be an omnivore by Modesto et al. (2007) on the basis of the dental morphology, is found to most likely have evolved from a herbivorous ancestor, rather than Captorhinikos chozaensis and the Moradisaurinae representing convergent transitions to herbivory from an omnivorous ancestor. There is, however, considerable uncertainty; the probability of a herbivorous ancestor is not much more than 50%. There is further uncertainty surrounding the ancestral diet of the clade containing the three species of Captorhinus, Captorhinikos chozaensis, Labidosaurus and the Moradisaurinae; while an omnivorous ancestor receives the highest likelihood, the probability is not much better than that of a carnivorous ancestor. This has implications for the transition to herbivory: the transition from carnivory to herbivory may have passed through an omnivorous phase, which was retained by the genus Captorhinus (Dodick & Modesto, 1995; Kissel, Dilkes & Reisz, 2002), or the genus Captorhinus may represent a transition to omnivory from carnivory independent of the transition to herbivory.

For the three uncertain nodes, the ancestral state reconstruction using the RJ MCMC method produced similar results. In fact this approach suggested less uncertainty surrounding the ancestral states (Fig. S1). Thus the most probable evolutionary history inferred is a single transition to herbivory via omnivory, with Labidosaurus representing a reversal to an omnivorous diet from a herbivorous ancestor.

Analyses of rate heterogeneity

In the overwhelming majority of the 100 time calibrated trees, a significant rate increase is identified along the branch leading to the Moradisaurinae (Fig. 4), the clade containing exclusively herbivorous taxa. This result is robust with taxonomic jack-knifing: in 94% of the jack-knife iterations a significant rate increase is identified along this branch. The position of other significant increases in rate depends on the tree topology and the uncertainty in dating the taxa, but in more than half of the trees the branch leading to the clade containing the Moradisaurinae, Labidosaurus and Captorhinikos chozensis (the clade inferred to have a herbivorous ancestor) are found to exhibit a rate increase, as is the lineage leading to the clade containing Labidosaurus and Moradisaurinae in more than two thirds of the trees. Significant rate decreases are observed in the lineages leading to Saurorictus and to Labidosaurus in the majority of the 100 trees.

While the analyses did identify rate heterogeneity when comparing branches of the phylogeny, when comparing rates of evolution in different time bins, very little was identified. In all of the 100 time calibrate trees, a constant rate through time was found to have a higher likelihood than a different rate in each time bin.

Rates and disparity in different dietary regimes

Of the 10,000 stochastic maps of dietary evolution in captorhinids, the herbivores have a higher mean rate of discrete character change than the omnivores in 9,845 maps, and a higher rate than the carnivores in 9,986 (Fig. 5A). When the mean rates of herbivores are compared to an equal number of branches drawn at random, the herbivores have a higher mean rate in 9,484 of the stochastic map (Fig. 5B). This result is robust with taxonomic jack-knifing: in 98% of jack-knife iterations the mean rates of character change in herbivores are higher than the carnivores, and are higher than the omnivores in 97% of the iterations (Fig. S2). In all 10,000 stochastic maps, the mean morphological distance between the herbivorous taxa is greater than that of the omnivores and the carnivores, indicating a greater disparity (Fig. 6).

Disparity through time

When evolutionary change is assumed to be gradual (Fig. 7A), the carnivorous captorhinids show a gradual increase in morphological disparity up to a peak in the early Artinskian. Through the late Artinskian and Kungurian their disparity decreases, culminating in a fall to zero across the Kungurian/Roadian boundary, after which only one carnivorous captorhinid is included in the phylogeny (Saurorictus). The omnivorous captorhinids show a similarly gradual increase in disparity between the Asselian and Kungurian. Again, their disparity falls to zero across the Kungurian/Roadian boundary. The initial establishment of the disparity of the herbivorous lineages is more rapid than that of the carnivores, having exceeded the disparity of the carnivorous captorhinids by the early Kungurian. A disparity peak is reached in the late Kungurian, higher than the peaks observed either in the carnivorous or omnivorous curves. Herbivore disparity falls across the Kungurian/Roadian boundary, but recovers by the Wuchiapingian.

If morphological change is assumed to be punctuated, with the morphological change occurring at the speciation events (Fig. 7B), then in all three dietary classes peak morphological disparity is reached soon after that regime’s appearance and disparity remains fairly constant in the bins following. As observed when using the gradualistic model, however, peak disparity of the herbivores is higher than either the omnivores or the carnivores. Interestingly, the disparity of herbivores already exceeds that of the other two dietary regimes by the Artinskian when using the punctuated model. Moreover, the decrease in disparity observed in the herbivorous lineages across the Kungurian/Roadian boundary is of a much lesser extent and disparity has recovered by the Wordian.

Lineage density

While there is considerable variation in the absolute lineage densities calculated within each dietary regime (Fig. 8), the overwhelming majority of the 10,000 evolutionary histories examined show the same relative pattern. In 9,926 of these, the lowest lineage density is found in the herbivores, while the highest lineage density is observed in the omnivores. In only 58 of the 10,000 stochastic maps do the carnivores have a lower lineage density than the herbivores.

Character change histories

The majority of character changes in the carnivorous lineages occurred in the skull and postcranium (Fig. 9A). In most of the 10,000 stochastic maps the feeding apparatus (teeth and mandible) remain more conservative, with a lower proportion of character changes occurring in these regions. This changes with the transition to herbivory: the majority of the characters changing in herbivorous captorhinids are dental characters and, in many of the stochastic maps (but not all), mandibular characters (Fig. 9C). The postcranium and skull, the most plastic regions in the carnivorous captorhinids, show lower proportions of character change in this new dietary regime. There is little difference in the proportions of characters changing in each region in the omnivorous captorhinids (Fig. 8B).