Article: Archosauromorph extinction selectivity during the Triassic–Jurassic mass extinction

Many traits have been linked to extinction risk among modern vertebrates, including mode of life and body size. However, previous work has indicated there is little evidence that body size, or any other trait, was selective during past mass extinctions. Here, we investigate the impact of the Triassic–Jurassic mass extinction on early Archosauromorpha (basal dinosaurs, crocodylomorphs and their relatives) by focusing on body size and other life history traits. We built several new archosauromorph maximum‐likelihood supertrees, incorporating uncertainty in phylogenetic relationships. These supertrees were then employed as a framework to test whether extinction had a phylogenetic signal during the Triassic–Jurassic mass extinction, and whether species with certain traits were more or less likely to go extinct. We find evidence for phylogenetic signal in extinction, in that taxa were more likely to become extinct if a close relative also did. However, there is no correlation between extinction and body size, or any other tested trait. These conclusions add to previous findings that body size, and other traits, were not subject to selection during mass extinctions in closely‐related clades, although the phylogenetic signal in extinction indicates that selection may have acted on traits not investigated here.

Mass extinction events have played a major role in shaping macroevolutionary trends through time (Barnosky et al. 2011). While selection is important in determining which taxa become extinct during both background extinction and mass extinction events (McKinney 1997), analysis of the fossil record indicates that the traits involved in selection have varied through time, during both extinction regimes and across different clades (Jablonski 2005; Finnegan et al. 2017). The high rates of extinction today are comparable to those during past mass extinctions, and as a result, understanding the processes involved in extinction during those events may prove useful for modern conservation (Barnosky et al. 2011).

Several approaches can be taken when trying to determine the influence of selection on extinction. One is to assume that the phylogenetic clustering of extinctions is a logical consequence of selection acting on phylogenetically conserved traits within closely related taxa (McKinney 1997). This approach is particularly valuable in clades for which trait data are lacking, but phylogenies are relatively robust (Soul & Friedman 2017). Previous studies of this phenomenon have indicated that results vary across clades and mass extinction events; for example, while bivalve extinctions during the end‐Cretaceous mass extinction were phylogenetically clustered (Roy et al. 2009), brachiopod extinctions during the end‐Ordovician mass extinction were not (Krug & Patzkowsky 2015). Given that ecological traits are often phylogenetically conserved, these results may reflect the contrasting ecological severities of these two mass extinction events (McGhee et al. 2013).

Previous studies have found little evidence of correlations between particular individual traits and extinction during mass extinctions (Jablonski & Raup 1995; Smith & Jeffery 1998). Body size is often linked to extinction risk in extant animals (Gaston & Blackburn 1995), as larger animals often possess traits linked to extinction vulnerability, such as small population sizes and low reproductive rates (Cardillo et al. 2005). In contrast, palaeontological studies focusing on mass extinction events have failed to find evidence of this correlation (Jablonski & Raup 1995), particularly when phylogenetic relationships are taken into account (Friedman 2009; Puttick et al. 2017a). Other traits relating to extinction risk during mass extinctions include diet and motility (Payne & Clapham 2012; Song et al. 2012), with evidence that ecological specialists and slow dispersers cannot tolerate rapid environmental change (Erwin 1998; Jablonski 2005). Geographical factors such as range size and occupied latitude have also been suggested, as these traits may influence survival during spatially heterogeneous disturbances (Erwin 1998; Powell 2007; Jablonski 2008). However, some authors have suggested that few, if any, traits reduce extinction risk during mass extinctions (Jablonski 2005; Dunhill & Wills 2015). In addition, most of these studies have focused on marine taxa, and extinction selectivity in the terrestrial realm in deep time has only recently come under investigation.

The Triassic–Jurassic mass extinction event (201 Ma; Cohen et al. 2015) was the third most ecologically severe mass extinction of the Phanerozoic (McGhee et al. 2013), resulting in a major biotic turnover in the latest Triassic (Preto et al. 2010). In the seas, reef ecosystems, bivalves, ammonites and scleractinian corals declined dramatically (Bond & Grasby 2017) and on land many terrestrial tetrapods became extinct (Lucas & Tanner 2015). While it is broadly accepted that the Central Atlantic Magmatic Province eruptions were responsible for these extinctions (Whiteside et al. 2010; Bond & Grasby 2017), the mechanisms remain uncertain, but probably involved global warming and ocean acidification and anoxia, as in other such events (Preto et al. 2010).

The Triassic–Jurassic mass extinction played an important role in the evolution of the Archosauromorpha, the clade comprising crocodiles and birds and their ancestors. Having originated in the middle Permian, the archosauromorphs represent a classic evolutionary radiation (Benton et al. 2014; Ezcurra et al. 2014). The radiation occurred in three phases: the first in the Early and Middle Triassic, following the devastation of ecosystems by the Permian–Triassic mass extinction; the second following the Carnian Pluvial Episode, around 232 Ma, when dinosaurs in particular radiated explosively (Bernardi et al. 2018); and the third following the Triassic–Jurassic mass extinction (Langer et al. 2010).

Archosauromorphs were highly disparate and diverse by the end of the Triassic (Brusatte et al. 2010; Stubbs et al. 2013), reaching a cosmopolitan distribution (Nesbitt 2011). While there has been considerable recent interest in the group, the quality of the early archosauromorph fossil record is relatively poor, and our knowledge of their evolution is incomplete (Nesbitt 2011; Ezcurra 2016). During the Triassic, archosauromorphs replaced basal therapsids as the dominant large‐bodied animals of the terrestrial realm (Sookias et al. 2012). Their increase in absolute body size during the Triassic and Jurassic (Turner & Nesbitt 2013) was the product of passive, rather than actively directional, evolution (Sookias et al. 2012). At the Triassic–Jurassic mass extinction, archosaurs declined, with the aetosaurs, phytosaurs and rauisuchians all becoming extinct (e.g. Nesbitt 2011). This resulted in a major transition within the group, as the crocodile‐line archosaurs (Pseudosuchia) diminished and the bird‐line archosaurs (Avemetatarsalia; most notably dinosaurs and pterosaurs) rose to dominance (Brusatte et al. 2008). Dinosaurs subsequently became one of the most successful groups of the Mesozoic, possessing a range of adaptations, occupying many different ecological niches and exhibiting a wide variety of body sizes (Langer et al. 2010; Benton et al. 2014).

Here, we construct several new supertrees, including 184 species of Triassic and Jurassic archosauromorphs, to test for the impacts of the Triassic–Jurassic mass extinction on the group using phylogenetic comparative methods. These supertrees are used to examine the degree of phylogenetic clustering in extinctions, and to test for correlation between several traits (including body size) and extinction, during this event.

Material and method

Supertree construction

We applied a two‐step process to tree selection for topologies to be used in subsequent macroevolutionary analyses: first, composite phylogenies were generated, representing eight hypotheses of archosauromorph relationships; then the relative fit of these hypotheses were compared to phylogenies estimated from cladistic matrices. For each composite topology we obtained a likelihood score and p‐value that indicated how well the informal tree was supported by the formally estimated trees based on cladistic data. The relative measure of fit was the Robinson–Foulds metric (Robinson & Foulds 1981) and all taxa from the source trees were present in each of the informal composite trees.

We sourced and re‐analysed six character matrices for published phylogenies of Triassic and Jurassic archosauromorphs: Early to Middle Triassic archosauromorphs (Ezcurra 2016), Late Triassic to Middle Jurassic archosauromorphs (Nesbitt 2011), pterosaurs (Andres et al. 2014), basal dinosaurs (Cabreira et al. 2016), phytosaurs (Kammerer et al. 2016) and aetosaurs (Parker 2016). Each of the character matrices from these phylogenetic studies was re‐analysed individually using the Mkv + Γ model in a maximum‐likelihood approach in IQ‐TREE (Nguyen et al. 2015). As the original studies used different phylogenetic methods and presented trees using different summary metrics, the trees were re‐analysed to make them comparable.

Eight informal composite supertrees, containing the 184 archosauromorph species included in the input phylogenies from the Middle Triassic to Early Jurassic, were assembled in Mesquite (Maddison & Maddison 2017). The two main input phylogenies formed the basis of the eight supertrees (A–H): Ezcurra (2016) for trees A, C, E and G, and Nesbitt (2011) for trees B, D, F and H. To these two backbone topologies, we added a combination of published dinosaurian and pseudosuchian relationships. For dinosaurs, we used a topology based on Cabreira et al. (2016) for trees A, B, E and F, and a topology based on Baron et al. (2017) for trees C, D, G and H. For trees A–D, we preserved the pseudosuchian relationships within the backbone trees (Ezcurra (2016) in trees A and C and Nesbitt (2011) in trees B and D), but used a pseudosuchian topology based on Brusatte et al. (2010) in trees E–H. The eight informal supertrees therefore incorporated all possible combinations of two contrasting archosauromorph topologies, two dinosaur topologies, and three pseudosuchian topologies (Fig. 1; Table 1).

Figure 1 Open in figure viewerPowerPoint The first of eight proposed Triassic and Jurassic archosauromorph supertrees (A). An asterisk indicates groups which are paraphyletic in this tree, and therefore should be considered basal members of their clade. Silhouettes from http://www.phylopic.org/ (clockwise from ‘Archosauromorpha’): Prolacerta T. Michael Keesey; Chanaresuchus ‘Smokeybjb’; Eudimorphodon Steven Traver; Coelophysis Emily Willoughby; Paleorhinus Scott Hartman; Riojasuchus Nobu Tamura & T. Michael Keesey; Stagonolepis Scott Hartman; Ticinosuchus Nobu Tamura; Polonosuchus Dmitry Bogdanov & T. Michael Keesey; Protosuchus Nobu Tamura & T. Michael Keesey. Table 1. The published phylogenies on which the eight tested supertree topologies are based Proposed supertree Basis of relationships Dinosaurian relationships Pseudosuchian relationships A Ezcurra (2016) Cabreira et al. (2016) Ezcurra (2016) B Nesbitt (2011) Cabreira et al. (2016) Nesbitt (2011) C Ezcurra (2016) Baron et al. (2017) Ezcurra (2016) D Nesbitt (2011) Baron et al. (2017) Nesbitt (2011) E Ezcurra (2016) Cabreira et al. (2016) Brusatte et al. (2010) F Nesbitt (2011) Cabreira et al. (2016) Brusatte et al. (2010) G Ezcurra (2016) Baron et al. (2017) Brusatte et al. (2010) H Nesbitt (2011) Baron et al. (2017) Brusatte et al. (2010)

The relative likelihood of each informal composite supertree based on the source trees was judged by comparing the maximum likelihood fit of each source tree to each supertree, and then summarizing across all tree‐wise likelihood scores to get the likelihood of each supertree. The fit of each source tree was assessed using the maximum‐likelihood metric of Steel & Rodrigo (2008), based on an exponential distance from the informal composite supertree using Robinson–Foulds distances (Robinson & Foulds 1981). Using this approach, a maximum‐likelihood score was estimated for each supertree using L.U.St (Akanni et al. 2014). The supertrees were then ranked according to their likelihood scores, with the significance of each tree based on the approximate unbiased (AU) test (Shimodaira 2002). The AU test is an extension of the Kishino–Hasegawa (KH) and Shimodaira–Hasegawa (SH) tests, which determine whether there is a significant difference between the likelihood fits of a set of trees. In all tests, the null hypothesis is found by bootstrapping and re‐centring the likelihood differences, but the AU test uses multiple lengths of these likelihood differences during bootstrapping, which reduces vulnerability to type‐1 errors and selection biases present in the KH and SH tests respectively (Shimodaira 2002).

We applied these tests to elucidate whether it was possible to reject statistically any of the informal composite trees prior to further analyses. This test differs only slightly from previous approaches that have used trees inferred using molecular data (e.g. Akanni et al. 2015; Arcila et al. 2017; Puttick et al. 2018) rather than morphological data. However, the approach remains comparable, as the metric is based upon Robinson–Foulds distances. This means the analyses use information from tree topology only, and therefore the relative likelihood fit of source trees can be assessed by testing any form of ‘supertree’, be it formal or informal. Although we re‐analysed all source matrices, a potential issue is that some of the source trees were used as a basis for building the composite phylogenies, as high support would be expected between a source tree and an informal composite tree based upon it. However, it is not obvious whether alternative arrangements of the informal composite topology can be rejected given the source trees, and so this requires formal testing. Additionally, it is advantageous to base macroevolutionary analyses on phylogenies that are supported by character data.

Data collection

The data are listed in Allen et al. (2018). Temporal occurrence data were obtained for all species included in the supertrees, on the basis of stage durations, from the Paleobiology Database. Stage dates were taken from Cohen et al. (2015), with an additional division arbitrarily placed at the midpoint of the Norian, because of the long duration of that stage (18.5 myr; Cohen et al. 2015).

Femur length has been shown to correlate with total body length in modern mammals (Biewener 1983; Christiansen 1999), and femur length is a commonly used proxy for total body mass in extinct vertebrates. We collected femur length data for 117 species, mainly from Sookias et al. (2012) and Turner & Nesbitt (2013). When there were measurements from multiple specimens of a species, median values were calculated. All data were log‐transformed prior to statistical analysis.

Information was also collected for diet (herbivore, omnivore, faunivore), posture (biped, quadruped, biped‐quadruped, biped‐flying), and habitat (terrestrial, semi‐aquatic); these data were designated by author opinion (MJB) following a literature search. A continuous palaeolatitude measure from the Paleobiology Database was also collected (for taxa with multiple localities, we used the palaeolatitude of the type specimen).

Tree‐scaling methods

Different tree‐scaling methods might have an impact on results from analysis of time‐trees (Bapst 2014; Soul & Friedman 2017). To find the uncertainty in the results attributable to the time‐scaling approach, each of the supertrees was time‐scaled using the Equal, Minimum Branch Length (MBL), cal3 and Hedman methods, with 100 repetitions of each method. Uncertainty in taxon age was refined to stage level, and for each of the 100 replicates a random age was chosen from uniform distributions bounded by the maximum and minimum age of stages for every taxon, so a unique age was chosen for each taxon at each iteration. The age of one species, whose range was limited to the Hettangian (the shortest stage within the covered timespan), was fixed at the midpoint of the stage, to retain a point within the tree from which absolute ages across the rest of the tree could be calculated. These methods, and all further analyses, were performed in R (v. 3.3.3) (R Core Team 2017), using the packages ape (v. 5.0; Paradis et al. 2004), paleotree (v. 3.0.0; Bapst 2012) and caper (v. 0.5.2; Orme 2018).

The Equal method estimates node ages using the age of the oldest descendant as the minimum bound. Zero‐length branches are then extended by sharing the branch length from the most immediately preceding branch that has a positive length (Brusatte et al. 2008). A root length of 5 myr was used when implementing this method. The MBL method also uses the age of the oldest descendant to determine node ages, but zero‐length branches are extended to a pre‐defined minimum length (Bapst 2012). Here, 2 myr was used as the minimum branch length.

The Hedman method is a Bayesian approach, with prior probability distributions for node ages constrained based on the ages of consecutive outgroups (Hedman 2010). Analyses using the Hedman method utilized the R code of Lloyd et al. (2016), with the necessary outgroups being basal taxa included in the Ezcurra (2016) phylogeny. Only the five species whose ranges predated those in the supertrees could be included, with the midpoints of their occurrence ranges used for the outgroup dates.

The cal3 method constrains nodes using oldest fossil occurrences, but estimates branch lengths using speciation, extinction and sampling rates based on the occurrence data (Bapst 2013). Analyses using the cal3 method were also conducted using the R code of Lloyd et al. (2016), with the sampling and extinction rate parameters calculated individually for each of the 100 trees, using a Foote (1997) approach. The speciation rate was set equal to the extinction rate (Lloyd et al. 2016). All cal3 trees contained internal zero‐length branches (Fig. 2), and therefore were not used in further analyses, as these lengths were considered to be biologically implausible.

Figure 2 Open in figure viewerPowerPoint Comparison of the effects of using the four different tree‐scaling methods, each applied to supertree A. The early Norian, late Norian, Rhaetian and Hettangian are indicated with grey shading from left to right on each tree.

Due to uncertainty in the timing and duration of the Triassic–Jurassic mass extinction event in the terrestrial realm (e.g. Lucas & Tanner 2015), analyses were carried out using six different time bins: late Norian, Rhaetian, Hettangian, late Norian–Rhaetian, Rhaetian–Hettangian, and early Norian–Rhaetian, with the Triassic–Jurassic boundary lying between the Rhaetian and the Hettangian.

Phylogenetic signal in extinction

To test for phylogenetic clustering of extinctions, the D‐statistic of Fritz & Purvis (2010) was used. It is generated by scaling and comparing the observed data to two simulated distributions, one calculated via a Brownian motion (BM) threshold model, the other produced using a randomization method, to simulate no phylogenetic signal (Fig. 3). A D‐statistic value close to 0 indicates that the collected data is consistent with the BM model, and hence that extinctions are phylogenetically clustered. A D‐statistic value close to 1 indicates that the extinctions are distributed randomly across the phylogenetic tree.

Figure 3 Open in figure viewerPowerPoint Example Phylo‐D plots for supertrees A and B, here having been scaled using the Equal method, with the late Norian–Rhaetian used as the extinction time bin.

Selectivity of traits

Correlation between femur length data and extinction was investigated using phylogenetic generalized least squares (PGLS) to account for non‐independence of the residuals arising from the shared history of traits (Felsenstein 1985). Pagel's lambda (λ; Pagel 1999), a measure of how phylogenetically conserved a trait is, was also estimated simultaneously within PGLS models to incorporate a measure of the degree of phylogenetic signal rather than assuming Brownian motion (Revell 2010).

Both the Phylo‐D and PGLS analyses were carried out using only taxa found within the designated time bin (‘contemporary’ taxa). In other words, all taxa that became extinct prior to the time bin were excluded, and lineages that passed through the bin were trimmed. For the purpose of the PGLS tests, StableTraits (Elliot & Mooers 2014) was used to estimate ancestral femur lengths for these trimmed lineages, which can generally be accurately reconstructed with fossil information (Puttick 2016).

Diet, posture, habitat, body size and palaeolatitude were analysed as predictors for extinction as a binary response, with a categorical link function with the phylogeny included as a random effect. These data were analysed in a Bayesian framework in MCMCglmm (Hadfield 2010) with the analyses run for 130 000 iterations, and with the first 30 000 discarded as ‘burn‐in’. The residual model priors were fixed to 1. For random effects parameter, expanded priors were used with the inverse Wishart parameters V and nu set to 0.1 and 1, and the prior mean set to 0 and variance 1000. Priors were chosen to allow for convergence, but not to overly bias parameter estimates, for example by selecting mean 0 with large variance (1000) for the prior mean. In these analyses, extinction was measured either by comparing taxa in the extinction bin to all other taxa in the phylogeny (from all time bins), or by comparing taxa in the extinction bin to taxa from the earliest Jurassic only. The ‘contemporary’ test used for the continuous femur data was not repeated here as it is necessary to elucidate ancestral states for this approach, and it would be difficult to marginalize over all uncertainty in the ancestral estimates for the discrete predictors of extinction. Therefore, for these analyses, the extinction bin was tested against the whole tree (to test if extinction was significantly different in this bin against all archosauromorph history considered), and between the time bin of interest and the Early Jurassic (to test whether extinction was significantly different between lineages that became extinct and lineages that survived). The model was fit to test for the difference in mean of each predictor in the model with no interactions. Convergence of the Bayesian analyses was assessed by examining trace plots of parameter estimation and inspection of effective sample sizes (target > 200).

Results

Supertrees

Of the eight archosauromorph supertrees, trees A and B have equal best likelihoods compared to the input phylogenies, and therefore have Δ likelihoods of 0. Based on these results, these two trees were used for the main analyses (Table 2), with supertree A used for the MCMCglmm procedure. Trees C, D and E, although not the best‐fit trees, did not have sufficiently low approximate unbiased (AU) values to be significantly rejected. Analyses carried out using trees C, D and E can be found in Allen et al. (2018). Trees F, G and H were all a poor fit (AU significance < 0.05) and were not used for subsequent analyses.

Table 2. Likelihood values and probabilities based on the approximately unbiased (AU) test of each supertree Proposed supertree Δ likelihood AU A 0 0.406 B 0 0.295 C 1.7 0.126 D 3.5 0.197 E 8.7 0.076 F 7.8 0.005 G 10.4 0.005 H 11.3 0.001

• AU values < 0.05 are indicated in italics; these trees were rejected prior to further analysis.

The four tree‐scaling methods produced considerable variation in branch lengths (Fig. 2). In particular, the MBL method pushed the root of the trees much deeper in time than the other methods, while the cal3 method produced many internal zero‐length branches, giving the appearance of polytomies at the base of several of the subgroups.

Phylogenetic clustering of extinction

Generally, there is evidence for phylogenetic signal in the extinction of taxa (Fig. 3; Table 3). Phylo‐D values vary considerably between time bins and tree‐scaling methods, but the corresponding probabilities indicate support for a Brownian motion model over the non‐phylogenetic random model. The lineage sample size, consisting of the number of taxa to become extinct during the time bin plus the number of lineages to survive to the end of the bin, varies considerably across time bins, in approximate accordance with the duration of the bin.

Table 3. The results of the Phylo‐D tests for supertrees A and B Time bin Method Extinct Survived D (SD) p (D > 0) Proportion p (D > 0) < 0.05 p (D < 1) Proportion p (D < 1) < 0.05 Tree A LN Equal 35 22 0.0274 (0.2726) 0.0122 0.93 0.4868 0.01 MBL 35 20 0.3505 (0.2657) 0.0366 0.84 0.2275 0.13 Hedman 34 26 −0.4194 (0.3281) 0.0010 1 0.7632 0 R Equal 10 15 −0.1545 (0.4271) 0.0298 0.81 0.5679 0 MBL 10 14 0.0547 (0.3306) 0.0362 0.78 0.4790 0 Hedman 10 17 −0.6040 (0.4847) 0.0156 0.94 0.7338 0 H Equal 5 11 0.3338 (0.3287) 0.1658 0 0.3750 0.01 MBL 5 10 0.4804 (0.2927) 0.2003 0.01 0.2932 0.04 Hedman 5 12 0.1265 (0.3543) 0.1471 0.05 0.4703 0 LN + R Equal 45 15 −0.7918 (0.1568) <0.0001 1 0.9193 0 MBL 45 14 −0.3557 (0.0701) <0.0001 1 0.7826 0 Hedman 44 17 −1.3382 (0.1886) <0.0001 1 0.9771 0 R + H Equal 15 11 0.6220 (0.3277) 0.2366 0.13 0.2057 0.1 MBL 15 10 0.8654 (0.2767) 0.3773 0.02 0.1093 0.41 Hedman 15 12 0.3799 (0.4129) 0.1721 0.22 0.3609 0.04 EN + LN + R Equal 71 15 −0.9472 (0.1225) <0.0001 1 0.9628 0 MBL 71 14 −0.5139 (0.0593) <0.0001 1 0.8854 0 Hedman 71 17 −1.4909 (0.1620) <0.0001 1 0.9907 0 Tree B LN Equal 35 22 0.0098 (0.2653) 0.0102 0.97 0.4955 0.01 MBL 35 20 0.3361 (0.2535) 0.0315 0.85 0.2312 0.15 Hedman 34 26 −0.4041 (0.3347) 0.0021 0.99 0.7536 0 R Equal 10 15 −0.1730 (0.4450) 0.0327 0.77 0.5783 0 MBL 10 15 0.1415 (0.3137) 0.0409 0.74 0.4270 0 Hedman 10 17 −0.6645 (0.4538) 0.0128 0.94 0.7579 0 H Equal 5 11 0.3787 (0.3509) 0.1813 0.04 0.3549 0 MBL 5 10 0.5351 (0.2794) 0.2131 0 0.2581 0.02 Hedman 5 12 0.1781 (0.3814) 0.1628 0.04 0.4557 0 LN + R Equal 45 15 −0.8151 (0.1645) <0.0001 1 0.9247 0 MBL 45 15 −0.3226 (0.0652) <0.0001 1 0.7664 0 Hedman 44 17 −1.4320 (0.2012) <0.0001 1 0.9822 0 R + H Equal 15 11 0.6971 (0.2666) 0.2663 0.03 0.1755 0.09 MBL 15 10 0.8931 (0.2404) 0.3848 0.02 0.0881 0.44 Hedman 15 12 0.4120 (0.3760) 0.1826 0.12 0.3443 0.01 EN + LN + R Equal 70 15 −0.9849 (0.1167) <0.0001 1 0.9662 0 MBL 71 15 −0.4784 (0.0528) <0.0001 1 0.8733 0 Hedman 71 17 −1.5676 (0.1673) <0.0001 1 0.9928 0

• Notation follows that of Fritz & Purvis (2010); ‘p (D < 1)’ denotes the proportion of D values that indicate significant phylogenetic signal, while ‘p (D > 0)’ denotes the proportion of D values that indicate significant deviation from a Brownian Motion model. Time bins: EN, early Norian; LN, late Norian; R, Rhaetian; H, Hettangian. All values are the mean of 100 repetitions, with standard deviations given in brackets for D values. Values in italics indicate p < 0.05. All values given to four decimal places except lineage extinction/survival.

The greatest contrast in results is between time bins. The most extreme contrast can be seen between the Hettangian and Rhaetian–Hettangian bins, which have positive D‐statistic means, and relatively large standard deviations of 0.2–0.4, and the late Norian – Rhaetian and the whole Norian – Rhaetian bins, which have negative D‐statistic means and relatively small standard deviations of 0.05–0.2. In addition, every tree has a probability of deviating from a Brownian Motion model of less than 0.05 for these two bins.

Considerable variation in D‐statistic values can also be seen across the different tree‐scaling methods. Generally, trees scaled with the MBL method have higher D‐statistic means and smaller standard deviations. However, trees scaled with the Hedman method have lower D‐statistic values and larger standard deviations. Trees scaled with the Equal method tend to lie between the two. However, D‐statistic values are similar across all tested topologies (Table 3; Allen et al. 2018).

Trait selectivity

Femur length shows high phylogenetic signal across the Archosauromorpha (Pagel's λ > 0.75), but there is little evidence of a significant correlation between femur length and extinction probability. Across all time‐scaling methods, only 1.11% of the analyses show a significant relationship between femur length and extinction (Table 4). The PGLS output values are relatively consistent between supertrees and across time bins, with the most variation seen between time‐scaling methods.

Table 4. PGLS coefficients and estimated for a regression of body size against extinction with supertrees A and B Time bin Method Extinct Survived Pagel's λ (SD) R² (SD) p Proportion p < 0.05 Tree A LN Equal 24 16 0.9645 (0.0524) −0.0126 (0.0160) 0.5862 0 MBL 24 15 0.7930 (0.0897) 0.0031 (0.0272) 0.4098 0.01 Hedman 24 20 0.9833 (0.0194) −0.0034 (0.0197) 0.4575 0.01 R Equal 7 12 0.8899 (0.1054) −0.0291 (0.0351) 0.5949 0.01 MBL 7 12 0.7549 (0.2618) 0.0016 (0.0567) 0.4285 0.01 Hedman 7 14 0.9685 (0.0503) −0.0154 (0.0378) 0.5148 0 H Equal 5 8 0.8990 (0.0328) 0.0631 (0.0948) 0.2473 0.09 MBL 4 8 0.9438 (0.0264) 0.0311 (0.0627) 0.2938 0 Hedman 5 9 0.9774 (0.0290) 0.0486 (0.0625) 0.2475 0.01 LN + R Equal 32 12 0.9517 (0.0459) −0.0184 (0.0041) 0.6703 0 MBL 32 12 0.8193 (0.0696) −0.0128 (0.0054) 0.5211 0 Hedman 31 14 0.9934 (0.0139) −0.0126 (0.0047) 0.5211 0 R + H Equal 12 8 0.9388 (0.0526) −0.0100 (0.0674) 0.5224 0.04 MBL 12 8 0.8278 (0.2297) −0.0385 (0.0173) 0.6394 0 Hedman 12 9 0.9901 (0.0159) −0.0312 (0.0258) 0.6420 0 EN + LN + R Equal 48 12 0.9462 (0.0504) −0.0140 (0.0025) 0.6950 0 MBL 47 12 0.8388 (0.0600) −0.0127 (0.0025) 0.6241 0 Hedman 47 14 0.9890 (0.0227) −0.0087 (0.0024) 0.4989 0 Tree B LN Equal 24 17 0.9695 (0.0500) −0.0121 (0.0155) 0.5720 0 MBL 25 15 0.8014 (0.0787) 0.0050 (0.0264) 0.3912 0.01 Hedman 24 21 0.9902 (0.0147) −0.0081 (0.0145) 0.5184 0 R Equal 7 12 0.8716 (0.1412) −0.0227 (0.0424) 0.5704 0 MBL 7 12 0.7826 (0.2651) −0.0067 (0.0575) 0.4878 0.04 Hedman 7 14 0.9702 (0.0603) −0.0151 (0.0388) 0.5133 0.01 H Equal 5 8 0.9003 (0.0330) 0.0737 (0.0976) 0.2312 0.12 MBL 4 8 0.9441 (0.0274) 0.0303 (0.0659) 0.2967 0.02 Hedman 5 9 0.9792 (0.0261) 0.0425 (0.0541) 0.2537 0 LN + R Equal 32 12 0.9561 (0.0453) −0.0183 (0.0034) 0.6589 0 MBL 32 12 0.8247 (0.0647) −0.0104 (0.0065) 0.4838 0 Hedman 32 14 0.9958 (0.0090) −0.0131 (0.0038) 0.5291 0 R + H Equal 12 8 0.9359 (0.0533) −0.0140 (0.0506) 0.5133 0.02 MBL 12 8 0.8431 (0.2319) −0.0301 (0.0333) 0.6011 0 Hedman 12 9 0.9913 (0.0138) −0.0354 (0.0186) 0.6673 0 EN + LN + R Equal 47 12 0.9576 (0.0443) −0.0133 (0.0026) 0.6527 0 MBL 47 12 0.8403 (0.0507) −0.0104 (0.0033) 0.5501 0 Hedman 47 14 0.9929 (0.0166) −0.0080 (0.0026) 0.4799 0

• Time bins: EN, early Norian; LN, late Norian; R, Rhaetian; H, Hettangian. All values are the mean of 100 repetitions, with standard deviations given in brackets for Pagel's λ and R² values. All values given to four decimal places except lineage extinction/survival.

The multivariate models provide no evidence of a role in extinction for femur length, diet, posture, habitat or palaeolatitude, from either the full tree or comparisons of the extinction boundary and the earliest Jurassic; these results are consistent across trees and dating methods (Table 5).

Table 5. Results from the multivariate MCMCglmm of extinction against various discrete extinction predictors for supertree A Posterior mean Lower 95% CI Upper 95% CI p Proportion p < 0.05 Equal (Intercept) 3.82 (−15.6, 19.5) −25.1 (−48.5, −11.7) 33.4 (14.2, 55.7) 0.643 (0.219, 0.989) 0 Length −1.17 (−3.59, 0.803) −5.47 (−8.36, −2.98) 3.06 (0.696, 5.55) 0.587 (0.105, 1) 0 Diet (Herbivore) −3.61 (−6.65, −0.828) −12.3 (−15.9, −9.61) 5.03 (1.49, 7.86) 0.384 (0.12, 0.838) 0 Diet (Omnivore) 10.4 (−10.6, 13.4) −10.5 (−33.1, −5.52) 30.7 (10.2, 35.3) 0.331 (0.183, 0.641) 0 Latitude 0.967 (−0.734, 3.66) −1.93 (−4.87, −0.253) 3.92 (2.96, 7.77) 0.51 (0.058, 0.951) 0 Posture (Biped flying) −3.03 (−5.36, 0.131) −17.7 (−22.1, −12.9) 11.1 (8.36, 13.6) 0.666 (0.426, 0.989) 0 Posture (Biped quadruped) 6.72 (3.25, 10.4) −3.99 (−9.56, 0.259) 17.5 (13.8, 25.1) 0.21 (0.047, 0.573) 0.02 Posture (Quadruped) −2 (−4.27, 0.872) −8.99 (−12.7, −4.65) 4.68 (2.39, 7.07) 0.558 (0.249, 0.997) 0 Habitat (Terrestrial) −0.584 (−5.23, 8.7) −13.3 (−18.7, −3.66) 12.5 (7.57, 21.1) 0.802 (0.137, 1) 0 MBL (Intercept) 4.05 (−16.4, 19.8) −27.7 (−48.8, −10.4) 36.3 (14.1, 53.4) 0.67 (0.19, 0.997) 0 Length −1.33 (−3.15, 1.19) −5.73 (−8.18, −2.76) 2.98 (1.04, 5.27) 0.528 (0.156, 0.981) 0 Diet (Herbivore) −3.7 (−7.23, −1.29) −13.3 (−17.7, −10.1) 5.6 (1.56, 8.7) 0.408 (0.105, 0.744) 0 Diet (Omnivore) 9.38 (−11.9, 13.7) −12.7 (−34.3, −3.43) 29.8 (7.75, 34.9) 0.341 (0.136, 0.569) 0 Latitude 0.963 (−0.983, 3.6) −2.7 (−4.94, 0.135) 4.83 (2.52, 7.54) 0.594 (0.048, 0.997) 0.01 Posture (Biped flying) −2.41 (−4.76, 0.559) −21.7 (−26.5, −15.4) 16.2 (13.1, 20) 0.797 (0.614, 0.982) 0 Posture (Biped quadruped) 7.07 (3.14, 10.1) −4.66 (−11, −0.128) 19.1 (15.3, 25.6) 0.222 (0.053, 0.604) 0 Posture (Quadruped) −2.93 (−4.85, −0.617) −10.1 (−12.7, −7.04) 4.03 (2.22, 7.42) 0.414 (0.188, 0.885) 0 Habitat (Terrestrial) −1.05 (−5.11, 6.35) −16.5 (−22.8, −8.82) 14.3 (8.87, 23.1) 0.75 (0.39, 0.993) 0 Hedman (Intercept) 2.69 (−14.7, 18.3) −27.2 (−52.4, −11.7) 33.1 (14.9, 54.6) 0.661 (0.279, 0.977) 0 Length −0.994 (−3.16, 0.775) −5.41 (−7.7, −3.22) 3.31 (1.28, 5.07) 0.654 (0.163, 0.985) 0 Diet (Herbivore) −3.38 (−5.54, −0.958) −11.6 (−14.3, −8.33) 4.77 (2.15, 7.98) 0.388 (0.175, 0.837) 0 Diet (Omnivore) 11.6 (−10.3, 14.6) −9.66 (−32.1, −4.56) 31.9 (9.6, 36.5) 0.277 (0.15, 0.6) 0 Latitude 1.13 (−1.39, 3.21) −2.11 (−5.25, −0.00246) 4.38 (2.58, 7.55) 0.514 (0.072, 0.992) 0 Posture (Biped flying) −2.04 (−4.19, −0.247) −14.5 (−18.3, −11.7) 10.2 (8.23, 12.3) 0.738 (0.488, 0.966) 0 Posture (Biped quadruped) 7.66 (4.11, 10.4) −3.52 (−8.13, 0.23) 18.3 (13.6, 24.5) 0.171 (0.039, 0.481) 0.01 Posture (Quadruped) −1.44 (−3.68, 1.2) −8.82 (−11.5, −5.75) 5.61 (3.1, 8.58) 0.718 (0.332, 0.999) 0 Habitat (Terrestrial) −0.801 (−6.16, 7.04) −13.4 (−19.4, −4.38) 11.7 (5.95, 19.5) 0.734 (0.244, 0.993) 0

• ‘CI’ refers to credible intervals. All values are the median of 100 repetitions, with the ranges given in brackets.

Discussion

Archosauromorph extinctions during the Triassic–Jurassic mass extinction event were not randomly distributed across the phylogeny, and therefore it is likely that clade‐level selection played a role in discriminating survivors and victims (Raup et al. 1973; Raup 1981). The results here show little evidence of a correlation between body size and extinction in Archosauromorpha during this mass extinction, supporting previous work that suggested no role of body size in archosauromorph subclade extinction across the Triassic–Jurassic boundary (Turner & Nesbitt 2013). In fact, our analyses indicate that there is no link between any studied ecological trait and extinction vulnerability for these taxa at this time. The aetosaurs, phytosaurs and rauisuchians, all of which became extinct, were ecologically diverse at the time (Stubbs et al. 2013), but so too were the clades that survived the mass extinction; particularly the dinosaurs, amongst which a great variety of body sizes and ecological traits can be seen (Brusatte et al. 2008). The only archosauromorph trait considered here that was largely lost across the mass extinction event was a semi‐aquatic lifestyle, which was most commonly associated with the Phytosauria (Stocker & Butler 2013).

Archosauromorph extinctions during the Norian and Rhaetian were phylogenetically clustered, and thus loss of lineages was not random with respect to phylogeny (Raup 1981; McKinney 1997). These patterns are evidenced in all time bins except the Hettangian and Rhaetian–Hettangian, in which there was insufficient evidence to support either phylogenetic clustering or a random distribution of extinctions across the phylogeny. This may relate to the fact that the Hettangian lasted only 2 myr (Cohen et al. 2015) and only a small number of lineages were included, just 15–17 for the Hettangian and 25–27 for the Rhaetian–Hettangian (Table 3). Phylo‐D tests are considered reliable only when applied to datasets with a minimum of 50 species (Fritz & Purvis 2010). However, these results may instead reflect a difference in evolutionary dynamics between the extinctions of the late Norian and Rhaetian, which are generally attributed to the Triassic–Jurassic mass extinction event, and the extinctions of the Hettangian, which may instead reflect the demise of disaster taxa during the subsequent period of biotic recovery (Preto et al. 2010).

Body size is not linked to extinction in the Archosauromorpha during the Triassic–Jurassic mass extinction event. This is in line with previous studies on the marine fossil record, particularly those of Friedman (2009), Sallan & Galimberti (2015) and Puttick et al. (2017a), who reached similar conclusions regarding the extinction of fishes during the end‐Cretaceous, end‐Devonian and end‐Permian mass extinction events respectively. Thus, evidence is accruing to suggest that body size was not a selective factor for vertebrates during past mass extinction events in either the marine or the terrestrial realm.

The finding that body size is not a selective factor is contrary to common assumptions concerning the extinction of the most famous extinct archosauromorph clade, the Dinosauria, at the end of the Cretaceous, 66 Ma. As Raup (1986) reported, it has often been said in popular literature that dinosaurs died out because they were large, whereas their relatives, birds and crocodiles, survived because they were small. This has not been subjected to a phylogenetically‐based analysis, but results might be difficult to assess as whole clades died out (non‐avian dinosaurs), whereas extinction and survival were more complex among bird and crocodilian taxa.

None of our analysed traits is significant in determining extinction or survival. Several factors could account for this. Poor preservation and lack of data availability have resulted in small sample sizes, reducing the statistical power of the analyses (Tables 3, 4). This is particularly true of the analyses of diet, posture, habitat and palaeolatitude, in which comparisons were made between sampled species only in coarse time bins. This contrasts with our ‘contemporary’ analyses of femur length, in which all available lineages, rather than sampled species in the bins, were used in the models.

The lack of a significant signal from the analyses of body size, diet, posture, habitat and palaeolatitude may be because the sample sizes were insufficiently large to detect the underlying evolutionary processes. Alternatively, perhaps there truly was no selection acting on these traits. It is possible that other characteristics were influential, such as geographical range size, but previous studies of terrestrial vertebrates during the Triassic–Jurassic mass extinction have suggested this is not the case (Dunhill & Wills 2015). Our finding of phylogenetic clustering of extinction in the Archosauromorpha might be a result of sampling bias rather than selective extinction (Soul & Friedman 2017). For example, it could be that closely related taxa became extinct because they were all small or lived in a particular unpreserved habitat, and so were less likely to fossilize. Alternatively, it could be a result of selective preservation of certain taxa in Lagerstätten. However, the Late Triassic archosauromorphs include taxa with disparate body sizes, diets, and other ecological characteristics. Regardless, here we found no evidence for trait‐related selection at this mass extinction event; a result that conforms with many previous analyses on other taxonomic groups, hitherto largely marine (Jablonski & Raup 1995; Smith & Jeffery 1998; Friedman 2009; Puttick et al. 2017a). This is in contrast to the plant fossil record (McElwain & Punyasena 2007).

Although the six chosen input phylogenies could have been analysed individually, constructing supertrees was the preferred approach, as larger taxon samples produce more reliable phylogenies (Zwickl & Hillis 2002) and reduce uncertainty in tests of phylogenetic signal (Fritz & Purvis 2010; Münkemüller et al. 2012). While these findings emerged from comparisons of formally inferred phylogenetic trees, we are here assuming that they apply also to supertrees. Many authors faced with a similar problem have used informal composite phylogenies, where multiple smaller trees are merged, but the outcome is rarely subjected to statistical testing (Akanni et al. 2014). The use of maximum likelihood topologies is potentially problematic, as the method can create a single fully‐resolved topology that gives false precision and little accuracy (Puttick et al. 2017b; O'Reilly et al. 2018). The accuracy of these maximum likelihood trees can be improved by accounting for uncertainty using bootstrapping and subsequent collapsing of poorly‐supported branches into polytomies (Brown et al. 2017; Puttick et al. 2017c; O'Reilly et al. 2018). The use of bootstrapping could be incorporated into similar comparisons between inferred trees and informal composite trees in the future, but at present the use of bootstrapping with this method would cause problems of pseudoreplication, as the input trees are assumed to be independent, but this is not the case with data from bootstrapping as they are based on partially overlapping data. Furthermore, the approach here potentially mitigates problems of inaccurate precision with maximum likelihood, as the inferred trees are not used directly, but are compared to assumed informal composite topologies.

The methods used in this study to test the proposed supertrees performed as expected, in that supertrees A and B were expected to show the best fit to the input matrices, being most heavily based on the relationships in the Ezcurra (2016) and Nesbitt (2011) trees from which two of the matrices were obtained (Table 1). However, the methods outperformed our expectations in being able to discriminate between the other supertrees, with the differences in likelihood and AU test outputs sufficient to reject three of the eight proposed topologies.

The relative likelihoods of the eight proposed supertrees lend more support to the traditional ornithischian and saurischian dinosaur clades than the sister‐group relationship between Ornithischia and Theropoda proposed by Baron et al. (2017). This may be because the taxon sample used in that study contrasted with that used by Cabreira et al. (2016), from which the matrix analysed here was obtained. The supertree likelihoods also support the crocodylomorphs as derived rauisuchians, rather than as the sister group to the Aetosauria, as proposed by Brusatte et al. (2010). However, the pseudosuchian fossil record is relatively fragmentary, and relationships within the group have been determined using small matrices with a high proportion of missing character data (Kammerer et al. 2016; Parker 2016), making this result uncertain.

While there was considerable variation in the statistical test output parameters, generally there was much less variation in the significance of the results, regardless of the tree topology or time‐scaling method used (Tables 3 and 4). This outcome contrasts with the conclusions of many recent studies which have highlighted that the choice of tree‐scaling method can be a major source of uncertainty in phylogenetic comparative analyses (Bapst 2014; Halliday & Goswami 2016; Soul & Friedman 2017).

Conclusions

Archosauromorph extinctions during the Triassic–Jurassic extinction event were phylogenetically clustered, and this suggests a role for selection. However, this study also adds to the mounting evidence that body size is not a key trait determining extinction or survival during past mass extinction events, and analysis of other traits also showed no evidence for selective extinction. It is possible that small sample sizes or coarse fossil sampling are limiting our ability to identify these correlations, or alternatively these extinctions were indiscriminate with respect to traits. Uncertainty in phylogenetic relationships, the use of different tree‐scaling methods and focusing on multiple time periods did not change the broad conclusions of these analyses. As such, these issues should not pose barriers to the application of the methods presented here to other taxonomic groups, providing bias is carefully considered in order to construct an appropriate statistical framework.

Author contributions

TLS, MJB and MNP conceived the project. BJA completed the data collection and calculations as part of her thesis for the MSc in Palaeobiology at the University of Bristol. Statistical approaches were led by MNP, and calculations were developed by BJA, TLS and MNP. BJA led the writing of the manuscript, and all authors contributed to the final version of the paper.

Acknowledgements

We thank members of the Bristol Palaeobiology group for feedback. We also thank David Bapst, Graeme Lloyd, Emily Rayfield and Tom Williams for their advice. Funding was partly supported by Natural Environment Research Council studentship NE/L002574/1 (to BJA). We also thank Martin Ezcurra and an anonymous reviewer for feedback that improved this manuscript.