Diversity and Disparity of Therocephalia: Macroevolutionary Patterns through Two Mass Extinctions (original) (raw)

Introduction

Therocephalians were a diverse group of non-mammalian eutheriodont therapsids from the synapsid lineage of amniotes, first documented in the fossil record in the early Wordian1. They most likely originated in southern Gondwana, since early therocephalian fossils such as Glanosuchus macrops were found in the Beaufort Group of South Africa1 and later spread throughout the world2. The clade includes a great variety of insectivorous (e.g. Bauriamorpha), herbivorous (Bauriidae), and carnivorous species with greatly enlarged canine teeth in both jaws (e.g. Scylacosauridae)3. They are one of three therapsid lineages to survive the end-Permian mass extinction (EPME), the largest mass extinction in Earth’s history, thought to have wiped out up to 96% of all marine species and 70% of terrestrial vertebrate species4,5,6,7,8. The surviving vertebrates went through a series of diversifications afterwards to fill new vacant niches6,7,8,9.

The other two synapsid lineages which survived the EPME, the anomodonts and cynodonts, have been the subject of recent analyses examining their macroevolutionary patterns across this event10,11, with particular focus on the patterns of diversity (species richness) and disparity (morphological diversity) through time. In anomodonts, the most diverse and abundant clade of herbivores at the time, diversity and disparity were found to be decoupled, with substantial loss of diversity across the EPME, but a slow and continuous decline in disparity beginning earlier in the Permian10. In cynodonts, the lineage which eventually gave rise to mammals during the Triassic, conflicting results are found depending on the disparity measure employed: variance-based measures of disparity conflict with the diversity estimates, while range-based measures show a significant correlation11.

Many studies have found that species richness and disparity are often not correlated12,13,14. As a direct result of specialization and speciation, one could expect increasing disparity early in an evolutionary history, since clades tend to diversify along ecomorphological lines in a new environment15, leading to “early bursts” of morphological diversity independent of species richness10,16,17. It has been suggested that such decoupling can become more pronounced during periods of mass extinction, as if the extinction is non-selective or targets less specialised forms, disparity may remain high while diversity falls17,18.

We used discrete morphological characters to determine patterns of morphological variation within Therocephalia, allowing an examination of how it has changed throughout their history. By quantifying three different disparity metrics, we illustrate the pattern of therocephalian evolution and extinction across the Permian-Triassic boundary.

Results

Phylogenetic Relationships of Therocephalia

The analyses described below, both relating to diversity and disparity, incorporate a phylogenetic hypothesis. To account for uncertainty in the relationships of therocephalians, three sets of phylogenies were produced. All were based on the same character and taxon dataset19 and were inferred and time calibrated using the same fossilised birth death (FBD) model under the same parameters (see material and methods section), but different topology constraints were applied: (1) an entirely unconstrained analysis; (2) relationships constrained to those found in the parsimony analysis of Kammerer and Masyutin19 (although topology was free to vary within the polytomies found by that study) (3) relationships constrained to those found in the parsimony analysis of Liu and Abdala20 (again, topologies were free to vary within the polytomies identified by that analysis).

Although it is not the purpose of this study to examine phylogenetic relationships of Therocephalia in detail (and in any case the different topology constraints appear to have had minimal effect on the diversity and disparity estimates), there are some differences between the topologies of the maximum clade credibility trees resulting from these analyses that warrant discussion.

The parsimony analysis of Liu and Abdala20 was unable to resolve the relationships among the Scylacosauridae, an assemblage of basal taxa found early in the fossil record of therocephalians. The modifications made by Kammerer and Masyutin19 found a monophyletic Scylacosauridae containing six species under parsimony analysis. The maximum clade credibility tree from reanalysis of this data matrix here using the FBD model also finds a monophyletic Scylacosauridae when no topology constraints are applied (Fig. 1A). However, when the FBD model is applied with the topology constrained to those found by Liu and Abdala20, Syclacosauridae is found to be paraphyletic (Fig. 1C), with Pristerognathus and Glanosuchus found more closely related to Eutherocephalia than other scylacosaurids. This result is somewhat paradoxical: one might expect the analysis where no topology constraints are applied to correspond more with the topology where scylacosaurid relationships are less constrained (that of Liu and Abdala).

Figure 1

Maximum clade credibility trees of Therocephalia identified by the fossilised birth death skyline model. (A) Tree produced by analysis with no topology constraint. (B) Topology constrained to that found by Kammerer and Maysutin19. (C) Topology constrained to that found by Liu and Abdala20.

Full size image

In other aspects, the early diverging nodes of the unconstrained analysis diverge from the topology constrained to that of Kammerer and Masyutin. The unconstrained analysis and that of Liu and Abdala found Scylacosuchus to lie outside the Eutherocephalia. The parsimony analysis of Kammerer and Masyutin19 was unable to resolve the early divergences of the Eutherocephalia but found Scylacosuchus within this clade (albeit in a polytomy with most of the whaitsioids). Perplexisaurus, however, was found by Kammerer and Masyutin to be outside Eutheroocephalia, but Liu and Abdala found it nested deep within this clade, a member of the family Chthonosauridae. The unconstrained FBD analysis again agrees with Liu and Abdala.

Relationships elsewhere in the three maximum clade credibility trees are largely consistent with each other, although there is some disagreement regarding the earliest divergences of the Bauriodea. Liu and Abdala17 were unable to resolve a monophyletic Ictidosuchidae (Ictidosuchus, Ictidosuchoides and Ictidosuchops) but did find a monophyletic Regisauridae (Regisaurus and Urumchia). Conversely, Kammerer and Masyutin16 found a monophyletic Ictidosuchidae but a paraphyletic Regisauridae. The unconstrained analysis here provides a novel set of relationships: Regisauridae is found within a paraphyletic Ictidosuchidae. This result is perhaps driven by the incorporation of tip ages into the FBD analysis: the relationships of Ictidosuchidae and Regisauridae in the unconstrained analysis are the most stratigraphically consistent, with the earlier appearing Ictidosuchidae forming successive outgroups to the Triassic Regisauridae.

Diversity

Phylogenetic diversity estimates calculated using the maximum clade credibility (MCC) trees identified by the three FBD analyses (Fig. 1) all indicate that diversity remained low until the late Capitanian (the latest stage of the Guadalupian) (Fig. 2). Diversity then increased rapidly, before remaining relatively constant throughout the Lopingian. When the tree produced by the unconstrained tip dating analysis is used, a slight further increase in diversity during the late Wuchiapingian is indicated (Fig. 2A), although this peak is absent in the trees produced by the topologically constrained analyses (Fig. 2B,C). This inconsistency is likely due to the fact that both the MCC trees produced in the constrained analyses support an earlier radiation of the Baurioidea during the late Guadalupian, while the unconstrained MCC indicates a later diversification of this group.

Figure 2

Phylogenetic diversity estimates of Therocephalia. The thick black line represents the diversity estimated using the maximum clade credibility tree, and the thin grey translucent lines indicate diversity estimated from 1000 trees drawn at random from the posterior distribution. (A) Diversity estimate based on trees produced by analyses with no topology constraints. (B) Based on trees constrained to topology found by Kammerer and Maysutin19. (C) Based on trees constrained to topology found by Liu and Abdala20.

Full size image

1000 randomly selected trees from the three posterior distributions produced by the tip dating analyses (Supplementary Data S1) also support an increase between the Guadalupian and the Lopingian, although there is variation in the pattern identified; whether the increase was stepwise or gradual is unclear (Fig. 2). However, a substantial increase between the early and late Capitanian is universally identified. The large range of diversity estimates obtained from the different posterior distributions during the late Guadalupian and early Lopingian may be due to the uncertainty surrounding the relationships of Scylacosauridae and early Eutherocephalia already discussed; it should be noted that the ranges are lower in the results from the trees with constrained topologies (Fig. 2B,C).

At the Permian-Triassic boundary, the phylogenetic diversity estimates calculated from all 1001 chosen phylogenies from all three tip dating analyses become considerably less variable (Fig. 2). The universal pattern identified is a decrease in species richness beginning across the Permian-Triassic boundary and continuing until the early Olenekian. While diversity does recover slightly in the late Olenekian, the Permian richness is never again reached. During the Anisian, diversity again declines, and therocephalians disappear from the fossil record in the Ladinian.

Disparity

Disparity was calculated using three metrics: the sum of variances (SOV), sum of ranges (SOR) and mean distance from centroid (DFC). The disparity curves calculated incorporated information on the ancestral morphologies deduced from the phylogenies (see Materials and Methods for more detail), and so disparity curves were calculated using all three phylogenetic hypotheses. However, the differences in the three sets of disparity curves were minimal (Figs 3 and S1, S2), so for the rest of the paper the disparity curves under discussion are those inferred from the tree produced by the unconstrained FBD analysis

Figure 3

Disparity through time estimated from the maximum clade credibility tree of the unconstrained analysis. The thick black line represents the raw disparity estimate, the thin translucent grey lines represent disparity estimates from 100 taxonomic bootstrapping replicates. (A) Disparity measured using the sum of variance; (B) Disparity measured using sum of ranges; (C) Disparity measured using mean distance from the centroid.

Full size image

All three disparity metrics, the sum of variances (SOV), sum of ranges (SOR) and mean distance from centroid (DFC) indicate a low initial disparity, which rises rapidly in the late Capitanian and then plateaus for much of the Lopingian (Fig. 3). When taxon bootstrapping is applied, this pattern is consistently observed when using the SOR and DFC measures of disparity (Fig. 3B,C), but the SOV measure shows much greater variation in the patterns identified (Fig. 3A).

Across the EPME, the raw SOV remains relatively stable, with only a slight decrease persisting into the Triassic (Fig. 3A). Again, bootstrapping the data leads to considerably more uncertainty, with a variety of trends observed and a wide range of disparity estimates. The SOR and DFC measures show similar changes across the EPME, which are robust to bootstrapping (Fig. 3B,C). Immediately prior to the extinction, in the late Changhsingian, both curves show a peak in disparity. Across the Permian-Triassic boundary both exhibit a significant decrease, but in the late Induan both rise again. During the Triassic, the SOR and DFC measures show a slight rise to a plateau in the Olenekian, but both decrease in the late Anisian prior to the extinction of Therocephalia.

Statistical tests

The strength of the correlations between the diversity and disparity varies depending on the metric applied. The correlation when SOV is used to measure the disparity is weak, but significant correlations are observed when SOR and DFC are used (Table 1).

Table 1 Results of statistical comparisons.

Full size table

When the bootstrapped disparities before the EPME are compared to those after, the significance of the change also varies depending on the disparity metric. The SOV indicates no significant change (possibly due to the large variation in patterns observed from the bootstrap data), but highly significant decreases are identified when the SOR and DFC are used.

Discussion

Comparison with Previous Analyses

Huttenlocker and Smith21 carried out a preliminary assessment of therocephalian diversity through time, as well as estimating origination and extinction rates through time. They used the results of this study to support a gradual rise of therocephalian diversity throughout the Permian. Origination and extinction rates were also inferred to rise continuously throughout the Permian, indicating a pattern of constant turnover with new species constantly and rapidly replacing the old.

The diversity curve of Huttenlocker and Smith21 contradicts the curve presented here particularly regarding the early evolution of therocephalians. We find, rather than a gradual increase in species richness throughout the Permian, instead a rapid increase in the latest Guadalupian is supported by the maximum clade credibility trees. The more rapid radiation in the late Guadalupian depicted in our study is due to the extension of multiple origination events back into the Capitanian stage. Huttenlocker and Smith21 argued for a “long fuse” model of therocephalian evolution, where the origin of clades such as the whaitsioids and the baurioids occurred in the Guadalupian, but the clades did not diversify until later. The tip dating analysis, on the other hand, suggests a “shorter fuse”, at least for the whaitsioids. The bulk of their diversification appears to have occurred by the end of the Guadalupian (Fig. 1), contributing to the more rapid rise in diversity in the late Capitanian. In fact, the bauriods may also have a shorter fuse as well; while the unconstrained FBD tree suggests a Lopingian diversification of this clade (Fig. 1A), when the topologies are constrained to those found by previous studies much of the diversification of bauriods occurred in the late Guadalupian (Fig, 1B,C), over a similar timescale to that of the whaitsoids.

The difference may also reflect the different methods employed to assess diversity. While Huttenlocker and Smith21 do not describe the method used to infer their diversity estimate, comparing the richness values exhibited in their Figure 11 to the lists of taxa and their ranges presented in their Tables 2 and 3, it appears that the diversity estimate represents a taxic diversity estimate; that is, a raw count of species. Such an estimate would be heavily influenced by the substantial sampling biases prevalent in the Permian-Triassic terrestrial record22,23,24,25. The phylogenetic diversity estimate employed here endeavours to include unsampled portions of the fossil record (inferred from the phylogeny) in the diversity estimate, and has been shown by multiple simulation studies to more closely represent the true diversity than the taxic diversity estimate26,27. Moreover, a more robust tip dating approach to generating the phylogenies, incorporating estimates of rates of origination, extinction, sampling and character change (see materials and methods) provides a more reliable indication on the origination age of clades than the observed fossil record.

Radiation Following the End-Guadalupian Extinction

Adaptive zone invasion is hypothesised to coincide with increased rates of both speciation and morphological evolution. Radiation into novel or under-exploited adaptive landscapes should maximise selection for divergence28,29. It is hypothesised that, as the lineage diversified rapidly along ecological lines, morphological diversity should accumulate rapidly early in the clade’s history: an “early burst” model. An examination of palaeontological data15, suggested that early increases in disparity are the most frequently observed pattern.

Previous examinations of synapsid diversity and disparity have supported an early burst model. Anomodont evolution was characterised by an initial peak in disparity followed by a slow decline10. Cynodonts also exhibit their most rapid rates of morphological diversification early in their history11. On the other hand, eureptile clades appear to have experienced delayed radiations. Archosauromorphs were found to have low disparity and rates of evolution early in their history, but radiated in the aftermath of the EPME, presumably filling the newly vacated areas of ecospace8. Captorhinidae also show this pattern: the peak in disparity and rates of evolution occur during the Kungurian and Roadian, despite the clade originating in the Carboniferous30,31. This was attributed to the late exploration of novel regions of ecospace – the herbivorous and omnivorous diet. Therocephalians also show relatively late exploration of novel diets (herbivorous members do not appear until the Triassic), and so their macroevolutionary patterns might be expected to be more similar to that of Captorhinids than of the more closely related synapsids.

The early evolution of Therocephalia exhibits a slower increase in disparity than those observed in anomodonts and cynodonts. Although they first appear in the fossil record during the Wordian, the greatest increase in all three disparity metrics does not occur until the late Capitanian, the last substage of the Guadalupian. A plateau is reached during the early Wuchiapingian, although the SOV measure of disparity indicates a temporary trough during the late Wuchiapingian (Fig. 3A). The increase in disparity does not precede the increase in diversity as in anomodonts, but instead both increase over a similar time.

The delayed diversification of Therocephalia, both of morphology and species richness, may be linked to the middle Permian mass extinction. A late Guadalupian mass extinction has been recognised both in the marine and terrestrial realms32,33,34,35,36,37. Using data on the tetrapods of the Karoo, it was demonstrated that the extinction on land occurred at the end of the Tapinocephalus Assemblage Zone (AZ), which radiometric ages suggested was before the end of the Capitanian38. In fact, their radiometric dates suggest that the Pristerognathus AZ (the post-extinction assemblage) falls entirely within the Capitanian. The principle victims of the mass extinction were the the Dinocephalia38, including the carnivorous Anteosauridae, and the predominantly (although not exclusively) herbivorous Tapinocephalia.

The late Capitanian diversification of Therocephalia may represent an adaptive radiation following the extinction event. According to the maximum clade credibility tree produced by tip dating (Fig. 1), it is at this time that the major therocephalian clade, the Eutherocephalia, diversified. The Whaitsioidea appear and radiate during the latest Guadalupian. When one examines the changes in morphospace occupation between the early and late Capitanian (Fig. 4B,C), the increase in disparity is seen to be driven predominantly by the whaitsoids. The Akidnognathidae and Baurioidea also appear to have their roots in the late Capitanian according to the tip dating analysis, although their radiation does not occur until later. The rapid diversification of these clades, driving increases in both disparity and diversity, may represent an ecological release of large carnivorous therocephalians following the demise of the previous dominant predators, the Anteosauridae.

Figure 4

Morphospace occupation of Therocephalians through time. (A) Phylomorphospace of all therocephalians, showing the position of both tips and internal nodes along Principal coordinates 1 and 2. Phylogeny is the maximum clade credibility tree of the unconstrained analysis. Colours represent affinities (those in black are not assigned to one of the five clades indicated in the legend). (B) Morphospace indicating the position of tips and internal nodes of therocephalians from the early Capitanian. Filled circles represent tips, squares represent internal nodes. (C) Late Capitanian morphospace. (D) Early Changhsingian morphospace. (E) Late Changhsingian morphospace. (F) Early Induan morphospace.

Full size image

Decoupling of Diversity and Disparity across the EPME

The impact of the end-Permian mass extinction on terrestrial vertebrates has been assessed in a large number of papers over the last decade10,24,25,39,40,41,42. These studies have highlighted several consistent patterns. First, the extinction was particularly severe for synapsids24,41, whose extinction rates were more than double those of both parareptiles and eureptiles41. Second, the recovery period of the extinction was unusually long, with continued ecosystem instability and elevated extinction rates throughout the Early Triassic41,42,43.

Among synapsids, the lineages terminated at the EPME were the two predatory lineages Biarmosuchia and Gorgonopsia. The cynodonts, of low diversity during the Permian, radiated both in terms of species richness and morphological diversity in the aftermath of the extinction11,44. The anomodonts were of extremely high species richness during the Permian, but their morphological diversity had been declining since the Guadalupian10. The decline in disparity does not seem to have been accelerated by the EPME, but species richness decreased dramatically and did not recover until the Middle Triassic.

The EPME represents the largest decrease in numbers of therocephalian species (Fig. 2). All 1001 tested phylogenetic diversity estimates support a species richness decline first across the Permian-Triassic boundary. The decline continues throughout the Early Triassic, with the number of lineages decreasing by more than half between the late Changhsingian and the early Olenekian. The Whaitsioidea, the most diverse lineage in the Lopingian, was entirely wiped out, and although the herbivorous Bauriidae radiated in the aftermath, therocephalian richness never returned to the levels observed in the Permian.

The SOR and DFC disparity curves of Therocephalia correlate significantly with the species richness curve (Table 1), and similarities may be observed before the latest Lopingian (rises in the late Capitanian to Wuchiapingian plateaus) and after the Induan (late Olenekian rises and late Anisian falls). However, diversity and disparity appear to have become decoupled in the stages surrounding the extinction events, and all three disparity curves contrast strongly with the species richness patterns. The SOV curve exhibits no significant change across the extinction event and remains high until the Middle Triassic (Fig. 3A). The SOV and DFC both increase immediately prior to the extinction in the late Changhsingian (Fig. 3B,C). There is a significant decline across the Permian-Triassic boundary (Table 1), but an almost immediate recovery during the Induan. By the Middle Triassic, both had returned to the levels of disparity observed prior to the extinction.

The contrasting patterns observed in the different disparity curves allow comment on the pattern of selectivity of the EPME. The increase in SOR and DFC immediately prior to the extinction indicates an increase in the total morphospace occupation, although the lack of change in the SOV and species richness suggests that this is due to a very small number of “specialists”; the bulk of the taxa continue to occupy the same region of morphospace. Across the EPME, there is again not substantial change in the SOV, but a significant SOR and DFC reduction. Two possible patterns are consistent with this observation: (1) the new specialists which appeared in the late Changhsingian were immediately culled by the extinction, leaving the original area of morphospace unchanged; (2) the original area was culled, and the new lineages diversified to occupy a volume of morphospace comparable to the earlier Lopingian values. The pattern of extinction observed in Therocephalia is more consistent with the latter explanation. The Whaitsioidea, which comprised the bulk of therocephalian diversity throughout the Lopingian, were eliminated in the EPME, while the Akidnognathidae and Baurioidea, which radiated in the latest Permian, survived into the Triassic (Fig. 1). This may be seen in the morphospace diagrams (although one should be careful of interpreting only two principal coordinate axes). Between the early and late Changhsingian, the baurioids and akidnognathids are the two clades that diversify, increasing morphospace occupation at both ends of PCo1 (Fig. 4D,E). During the Induan, both clades remain in similar regions of morphospace, but the whaitsoids disappear, considerably reducing disparity along PCo2 (Fig. 4F).

The Baurioidea make up the bulk of therocephalian diversity during the Triassic. Although their species richness never reaches that of the Permian, they include a wide range of morphologies (although they are restricted to small sizes), including herbivorous and potentially semi-aquatic forms43,45, and account for the morphological diversity of the Triassic being nearly equal to that of the Permian. However, by the middle of the Triassic they were very much in decline, restricted to the herbivorous Bauriidae, and all three disparity metrics show a decrease in the late Anisian (Fig. 3). Therocephalians are not known from beyond the Anisian, although their disappearance does not appear to be associated with any mass extinction event, but rather with an extended episode of elevated extinction rates and declining diversity that affected synapsids throughout the Middle and Late Triassic10,41.

Conclusions

The results presented here suggest that taxonomic diversity and discrete character disparity were correlated during the initial adaptive radiation of therocephalians but became decoupled during the EPME. The Permian-Triassic mass extinction only had a dramatic impact on therocephalian species richness, but morphological diversity was largely unchanged and remained high for 10 million years after (although the number of species never fully recovered). A greater impact on their morphological evolution was caused by the end-Guadalupian extinction, which allowed an “ecological release” causing the taxa to diversify both in terms of number of species and morphological diversity.

Materials and Methods

Phylogeny and Character Data

The dataset used in this study was taken from the most recent and comprehensive phylogenetic analysis of therocephalians19. The data matrix includes 49 therocephalian taxa with 136 discrete characters representing morphological variations across the entire skeleton. The first appearance date (FAD) and last appearance date (LAD) in millions of years ago for all taxa within the phylogeny were compiled from the Paleobiology Database, downloaded via the Fossilworks platform (www.fossilworks.org). Ranges of species are treated as continuous from the first to last appearance, ignoring possible gaps in the fossil record. The species ranges are presented in Supplementary Table S1.

In order to produce a phylogenetic framework for the diversity analysis, a Bayesian tip-dating approach was employed that uses the character data and the ages of the tip taxa to simultaneously infer the phylogenetic relationships and node ages: the Fossilised Birth Death (FBD) model46,47. This method was implemented in MrBayes 3.2.648. An offset exponential root prior was applied, with a minimum age 1 million years before the first appearance of therocephalians and a mean age at the Cisuralian/Guadalupian boundary. The choice of prior is based on previous analyses that have suggested a major turnover in dominant fauna at this time, along with the principal diversification of therapsids39,41,49 and so this is not an unreasonable time to expect therocephalians to have originated. In any case this age is not a hard lower bound; the exponential root prior allows root ages to be extended beyond the mean of the distribution. The analysis was run four times for 10,000,000 generations, sampling every 1000, thus producing a posterior distribution of 40,000 trees. 25% of these were discarded as burn-in. The full posterior distribution of trees before burnin is found in Supplementary Data S1, and the maximum clade credibility tree in Fig. 1A.

Due to disagreement and uncertainty regarding the relationships within therocephalians in recent studies, two further phylogenetic analyses were carried out. The parameters were identical to those described above, with the exception that the topology was constrained to represent the topologies found in two recent studies. The first topology constraint is that found by a parsimony analysis by Kammerer and Masyutin19, the second is that found by a parsimony analysis of an older version of the same matrix by Liu and Abdala20. Kammerer and Masyutin were unable to resolve the relationships of the early diverging eutherocephalian clades relative to each other, and indeed failed to resolve a monophyletic Whaitsioidea. Liu and Abdala20 were able to resolve these relationships, although they found Syclacosauridae to form a polytomy at the base of the therapsid tree, while Kammerer and Masyutin19 found this clade to be monophyletic. The maximum clade credibility trees of these analyses are shown in Fig. 1B,C, and the posterior distributions of trees are presented in Supplementary Data S2 and S3.

It should be noted that the matrix originally used by Liu and Abdala20 was not identical to that of Kammerer and Masyutin19; these authors made certain changes, including the removal of autapomorphies. Moreover, Liu and Abadala’s matrix did not include the later-described taxon Gorynychus. Because of the removal of autapomorphies, important both for tip-dating and disparity analyses50,51, it was decided not to examine disparity using Liu and Abdala’s matrix, but to instead, as described above, carry out the tip-dating and disparity analyses using the Kammerer and Masyutin19 matrix, but to test the impact of constraining the topology to that found by Liu and Abdala20.

Diversity

Phylogenetic diversity estimates were calculated as species richness estimated using the phyloDiv function as implemented in R package paleotree52. This method accounts for the incomplete sampling of the fossil record by including ghost lineages (lineages not sampled in the fossil record but inferred to have been present from the phylogeny) in the diversity estimate26,27. To account for uncertainty in the ages of lineages and phylogenetic relationships, phylogenetic diversity estimates were calculated for the three maximum clade credibility trees, as well as a random selection of 1000 trees from each of the Bayesian posterior distributions.

The period of time under study was split into 12 time bins, from the first appearance of therocephalians in the fossil record in the Wordian until their extinction in the Anisian. The time bins represent substages: each international stage was split into two, with the boundary at the midpoint of the stage.

Disparity

The 130 discrete morphological characters from the character matrix used to infer the phylogeny were also used to represent morphology in calculating disparity through time (note that as this matrix was designed for inferring phylogeny, the characters are only those considered parsimony informative). As well as the character scores provided for the tip taxa in the character matrix, character scores were obtained for the nodes within the tree with maximum-likelihood ancestral state reconstruction using functions in the R package Claddis53. The internal nodes could therefore be included as data points in the disparity calculation, thus allowing for ancestral morphologies which may have been outside the range of morphologies observed in the descendants54. We computed the pairwise morphological distances between each taxon as Gower dissimilarity, and performed a multi-dimensional scaling (principal coordinates analysis), again in Claddis. This produces a series of linearly uncorrelated continuous measurements representing each taxon’s position in morphospace.

The disparity through time was calculated using the same bins as diversity through time, using the R package dispRity55. This package allows the incorporation of information from the phylogeny into disparity studies, using the morphologies inferred from the internal nodes (described above).

The disparity of the principal coordinate values in each time bin was assessed using three metrics: sum of ranges (SOR), sum of variances (SOV) and mean Euclidean distance of each taxon from the centroid (DFC). SOR and DFC are measures of the total morphospace occupied but are vulnerable to sample size12 and morphologically distant singular taxa56 but this may be compensated by bootstrapping the data in dispRity. SOV is a measure of the spread of the data within the total morphospace. For each of the three disparity measures, 100 taxon bootstrapping replicates were carried out. The resilience of the methods used to sampling heterogeneity may be shown by correlation tests between the bootstrapped disparity curves and the observed sample size in each time bin using the spearman’s rank correlation coefficient. Only three of the 100 SOV curves showed significant (p < 0.05) correlation with sample size (Supplementary Data S4). The range-based metrics performed less well; 15 of the SOR curves and 13 of the DFC curves showed a significant correlation with sample size (Supplementary Data S5 and S6). Nevertheless, the overwhelming majority of all curves did not, indicating that bootstrapping and incorporation of information from the phylogeny go some way to alleviating the impact of sampling heterogeneity.

Statistical Tests

The correlation between the three disparity curves and the diversity estimate were calculated using the Spearman’s rank correlation coefficient in R. Prior to the correlation test, the time series were corrected using generalised differencing57 to account for autocorrelation (the temporal non-independence of the data). To examine the significance of disparity changes across the Permian-Triassic boundary, disparity values of the 100 bootstrap replicates of the two time bins either side of the extinction event were compared using a Wilcoxan Signed Rank Test test.

Data Availability

All data is included in the Supplementary Data Files.

References

1. Johnson, M. R., Van Vuuren, C. J., Hegenberger, W. F., Key, R. & Show, U. Stratigraphy of the Karoo Supergroup in southern Africa: an overview. J Afr Earth Sci 23, 3–15 (1996).
   Article Google Scholar
2. Abdala, F., Rubidge, B. S. & Van Den Heever, J. The oldest therocephalians (Therapsida, Eutheriodontia) and the early diversification of Therapsida. Palaeontology 51, 1011–1024 (2008).
   Article Google Scholar
3. Rubidge, B. S. & Sidor, C. A. Evolutionary patterns among Permo-Triassic therapsids. Annu Rev Ecol Syst 32, 449–480 (2001).
   Article Google Scholar
4. Benton, M. J. When Life Nearly Died: The Greatest Mass Extinction of All Time. Thames and Hudson (2015).
5. Benton, M. J., Tverdokhlebov, V. P. & Surkov, M. V. Ecosystem remodelling among vertebrates at the Permian-Triassic boundary in Russia. Nature 432, 97–100 (2004).
   Article ADS CAS Google Scholar
6. Irmis, R. B. & Whiteside, J. H. Delayed recovery of non-marine tetrapods after the end-Permian mass extinction tracks global carbon cycle. Proc R Soc B 279:1310–1318 (2012).
7. Smith, R. M. H. & Botha-Brink, J. Anatomy of a mass extinction: sedimentological and taphonomic evidence for drought-induced die offs at the Permo-Triassic boundary in the main Karoo Basin, South Africa. Palaeogeogr, Palaeoclim, Palaeoecol 396, 99–118 (2014).
   Article Google Scholar
8. Ezcurra, M. D. & Bulter, R. J. The rise of the ruling reptiles and ecosystem recovery from the Permo-Triassic mass extinction. Proc R Soc B 285, 20180361 (2018).
   Article Google Scholar
9. Huttenlocker, A. K. Body size reductions in nonmammalian eutheriodont therapsids (Synapsida) during the end-Permian mass extinction. PLoS One 9, e87553 (2014).
   Article ADS Google Scholar
10. Ruta, M., Angielczyk, K. D., Fröbisch, J. & Benton, M. J. Decoupling of morphological disparity and taxic diversity during the adaptive radiation of anomodont therapsids. Proc R Soc B 280, 20131071 (2013).
    Article Google Scholar
11. Ruta, M., Botha-Brink, J., Mitchell, S. A. & Benton, M. J. The radiation of cynodonts and the ground plan of mammalian morphological diversity. Proc R Soc B 22, 20131865 (2013).
    Article Google Scholar
12. Foote, M. The evolution of morphological diversity. Annu Rev Ecol Syst 28, 129–152 (1997).
    Article Google Scholar
13. Wills, M. A. Morphological disparity: a primer. In J. M. Adrain, G. D. Edgecombe, B. S. Lieberman eds Fossils, Phylogeny, and Form – an Analytical Approach. Kluwer Academic, New York: pp 55–144. (2001).
14. Erwin, D. H. Disparity: morphological pattern and developmental context. Palaeontology 50, 57–73 (2007).
    Article Google Scholar
15. Hughes, M., Gerber, S. & Wills, M. A. Clades reach highest morphological disparity early in their evolution. Proc Nat Acad Sci 110, 13875–13879 (2013).
    Article ADS CAS Google Scholar
16. Bapst, D. W., Bullock, P. C., Melchin, M. J., Sheets, H. D. & Mitchell, C. E. Graptoloid diversity and disparity became decoupled during the Ordovician mass extinction. Proc Nat Acad Sci 109, 3428–3433 (2012).
    Article ADS CAS Google Scholar
17. Hopkins, M. J. Decoupling of taxonomic diversity and morphological disparity during decline of the Cambrian trilobite family Pterocephaliidae. J Evol Biol 26, 1665–1676 (2013).
    Article CAS Google Scholar
18. Foote, M. Contributions of individual taxa to overall morphological disparity. Paleobiology 19, 403–419 (1993).
    Article Google Scholar
19. Kammerer, C. F. & Masyutin, V. A new therocephalian (Gorynychus_masyutinae gen. et sp. nov.) from the Permian Kotelnich locality, Kirov Region, Russia. PeerJ 6, e4933 (2018).
    Article Google Scholar
20. Liu, J. & Abdala, F. The tetrapod fauna of the upper Permian Naobaogou Formation of China: 1. Shiguaignathus wangi gen. et sp. nov., the first akidnognathid therocephalian from China. PeerJ 5, e4150 (2017).
    Article Google Scholar
21. Huttenlocker, A. K. & Smith, R. M. New whaitsioids (Therapsida: Therocephalia) from the Teekloof Formation of South Africa and therocephalian diversity during the end-Guadalupian extinction. PeerJ 5, e3868 (2017).
    Article Google Scholar
22. Benson, R. B. & Upchurch, P. Diversity trends in the establishment of terrestrial vertebrate ecosystems: interactions between spatial and temporal sampling biases. Geology 41, 43–46 (2013).
    Article ADS Google Scholar
23. Brocklehurst, N., Kammerer, C. F. & Fröbisch, J. The early evolution of synapsids and the influence of sampling on their fossil record. Paleobiology 39, 470–490 (2013).
    Article Google Scholar
24. Fröbisch, J. Vertebrate diversity across the end-Permian mass extinction - separating biological and geological signals. Palaeogeog, Palaeoclim, Palaeoecol 372, 50–61 (2013).
    Article ADS Google Scholar
25. Fröbisch, J. Synapsid diversity and the rock record in the Permian-Triassic Beaufort Group (Karoo Supergroup), South Africa. In C. F. Kammerer, K. D. Angielczyk, J. Fröbisch eds Early Evolutionary History of the Synapsida, Springer, Dordrecht: pp 305–320 (2014).
26. Lane, A., Janis, C. M. & Sepkoski, J. J. Estimating paleodiversities: a test of the taxic and phylogenetic methods. Paleobiology 31, 21–34 (2005).
    Article Google Scholar
27. Brocklehurst, N. A simulation-based examination of residual diversity estimates as a method of correcting for sampling bias. Palaeontol Electronica 18, 1–15 (2015).
    Google Scholar
28. Givnish, T. J. Ecology of plant speciation. Taxon 59, 1326–1366 (2010).
    Article Google Scholar
29. Givnish, T. J. Adaptive radiation versus ‘radiation’ and ‘explosive diversification’: why conceptual distinctions are fundamental to understanding evolution. New Phytol 207, 297–303 (2015).
    Article Google Scholar
30. Brocklehurst, N. Rates of morphological evolution in Captorhinidae: an adaptive radiation of Permian herbivores. PeerJ 5, e3200 (2017).
    Article Google Scholar
31. Romano, M., Brocklehurst, N. & Fröbisch, J. Discrete and continuous character-based disparity analyses converge to the same macroevolutionary signal: a case study from captorhinids. Sci Rep 7, 17531 (2017).
    Article ADS Google Scholar
32. Stanley, S. M. & Yank, X. A double mass extinction at the end of the Paleozoic era. Science 266, 1340–1344 (1994).
    Article ADS CAS Google Scholar
33. Rubidge, B. S. 27th Du Toit memorial lecture: re-uniting lost continents. Fossil reptiles from the ancient Karoo and their wanderlust. S Afr J Geol 108, 135–172 (2005).
    Article Google Scholar
34. Ota, A. & Isozaki, Y. Fusuline biotic turnover across the Guadalupian-Lopingian (Middle-Upper Permian) boundary in mid-oceanic carbonate buildups: biostratigraphy of accreted limestone in Japan. J Asian Earth Sci 26, 353–368 (2006).
    Article ADS Google Scholar
35. Bond, D. P. G. et al. The Middle Permian (Capitanian) mass extinction on land and in the oceans. Earth Sci Rev 102, 100–116 (2010).
    Article ADS Google Scholar
36. Stevens, L. G., Hilton, J., Bond, D. P., Glasspool, I. J. & Jardine, P. E. Radiation and extinction patterns in Permian floras from North China as indicators for environmental and climate change. J Geol Soc 168, 607–619 (2011).
    Article Google Scholar
37. Wignall, P. B. The End‐Permian mass extinction–how bad did it get? Geobiology 5, 303–309 (2007).
    Article Google Scholar
38. Day, M. O. et al. When and how did the terrestrial mid-Permian mass extinction occur? Evidence from the tetrapod record of the Karoo Basin, South Africa. Proc R Soc B 282, 20150834 (2015).
    Article Google Scholar
39. Sahney, S. & Benton, M. J. Recovery from the most profound mass extinction of all time. Proc R Soc B 275, 759–765 (2008).
    Article Google Scholar
40. Benton, M. J., Ruta, M., Dunhill, A. M. & Sakamoto, M. The first half of tetrapod evolution, sampling proxies, and fossil record quality. Palaeogeo, Palaeoclim, Palaeoecol 372, 18–41 (2013).
    Article ADS Google Scholar
41. Brocklehurst, N., Ruta, M., Müller, J. & Fröbisch, J. Elevated extinction rates as a trigger for diversification rate shifts: early amniotes as a case study. Sci Rep 5, 17104 (2015).
    Article ADS CAS Google Scholar
42. Roopnarine, P. D., Angielczyk, K. D., Wang, S. C. & Hertog, R. Trophic network models explain instability of Early Triassic terrestrial communities. Proc R Soc B 274, 2077–2086 (2007).
    Article Google Scholar
43. Ivakhnenko, M. F. Permian and Triassic therocephals (Eutherapsida) of eastern Europe. Paleontol J 45, 981–1144 (2011).
    Article Google Scholar
44. Lukic-Walther, M., Brocklehurst, N., Kammere, C. F. & Fröbisch, J. Diverity patterns of nonmammalian cynodonts (Synapsida, Therapsida) and the impact of taxonomic practice and research history on diversity estimates. Paleobiology In Press (2018).
45. Abdala, F., Jashashvili, T., Rubidge, B. S. & van den Heever, J. New material of Microgomphodon oligocynus (Eutherapsida, Therocephalia) and the taxonomy of southern African Bauriidae. In C. F. Kammerer, K. D. Angielczyk, J. Fröbisch eds, Early Evolutionary History of the Synapsida, Springer, Dordrecht: pp. 209–231 (2014).
46. Stadler, T. Sampling-through-time in birth-death trees. J Theor Biol 267, 396–404 (2010).
    Article ADS MathSciNet Google Scholar
47. Heath, T. A., Huelsenbeck, J. P. & Stadler, T. The fossilized birth-death process for coherent calibration of divergence-time estimates. Proc Nat Acad Sci 111, E2957–E2966 (2014).
    Article ADS CAS Google Scholar
48. Ronquist, F. & Huelsenbeck, J. P. MrBayes 3: Bayesian phylogenetic inference under mixed models. Bioinformatics 19, 1572–1574 (2003).
    Article CAS Google Scholar
49. Brocklehurst, N., Day, M. O., Rubidge, B. S. & Fröbisch, J. Olson’s Extinction and the latitudinal biodiversity gradient of tetrapods in the Permian. Proc R Soc B 284, 20170231 (2017).
    Article Google Scholar
50. Gould, S. J. The disparity of the Burgess Shale arthropod fauna and the limits of cladistic analysis: why we must strive to quantify morphospace. Paleobiology 17, 411–423 (1991).
    Article Google Scholar
51. Matzke, N. J. & Irmis, R. B. Including autapomorphies is important for paleontological tip-dating with clocklike data, but not with non-clock data. PeerJ 6, e4553 (2018).
    Article Google Scholar
52. Bapst, D. W. paleotree: an R package for paleontological and phylogenetic analyses of evolution. Methods Ecol Evol 5, 803–807 (2012).
    Article Google Scholar
53. Lloyd, G. T. Estimating morphological diversity and tempo with discrete character-taxon matrices: implementation, challenges, progress, and future directions. Biol J Linn Soc 118, 131–151 (2016).
    Article Google Scholar
54. Brusatte, S. L., Montanari, S., Yi, H.-Y. & Norell, M. A. Phylogenetic corrections for morphological disparity analysis: new methodology and case studies. Paleobiology 37, 1–22 (2011).
    Article Google Scholar
55. Guillerme, T. dispRity: a package for measuring disparity in R. Zenodo, 10.5281/zenodo.55646 (2016).
56. Korn, D., Hopkins, M. J. & Walton, S. A. Extinction space— method for the quantification and classification of changes in morphospace across extinction boundaries. Evolution 67, 2795–2810 (2013).
    PubMed Google Scholar
57. McKinney, M. L. Classifying and analyzing evolutionary trends. In McNamara, K. J. ed. Evolutionary Trends. Belhaven, London: pp. 28–58 (1990).

Download references

Acknowledgements

This study was financially supported by a Sofja Kovalevskaja Award of the Alexander von Humboldt Foundation donated by the German Federal Ministry for Education and Research to JF. We are grateful to Graeme Lloyd, Roger Close and Robin Beck for assistance with Claddis and MrBayes.

Author information

Authors and Affiliations

1. Institut für Biologie, Humboldt-Universität zu Berlin, Invalidenstraße 42, Berlin, D-10115, Germany
   Henrik Richard Grunert & Jörg Fröbisch
2. Museum für Naturkunde, Leibniz-Institut für Evolutions- und Biodiversitätsforschung, Invalidenstraße 43, D-10115, Berlin, Germany
   Henrik Richard Grunert, Neil Brocklehurst & Jörg Fröbisch
3. Department of Earth Sciences, University of Oxford, South Parks Road, Oxford, OX1 3AN, UK
   Neil Brocklehurst

Authors

1. Henrik Richard Grunert
   You can also search for this author inPubMed Google Scholar
2. Neil Brocklehurst
   You can also search for this author inPubMed Google Scholar
3. Jörg Fröbisch
   You can also search for this author inPubMed Google Scholar

Contributions

H.R.G. collected and analysed data, wrote the paper. N.B. analysed data, wrote the paper. J.F. wrote the paper.

Corresponding author

Correspondence toNeil Brocklehurst.

Ethics declarations

Competing Interests

The authors declare no competing interests.

Additional information

Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary information

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

Reprints and permissions

About this article

Cite this article

Grunert, H.R., Brocklehurst, N. & Fröbisch, J. Diversity and Disparity of Therocephalia: Macroevolutionary Patterns through Two Mass Extinctions.Sci Rep 9, 5063 (2019). https://doi.org/10.1038/s41598-019-41628-w

Download citation

• Received: 21 August 2018
• Accepted: 07 March 2019
• Published: 25 March 2019
• DOI: https://doi.org/10.1038/s41598-019-41628-w