How to date a crocodile: estimation of neosuchian clade ages and a comparison of four time‐scaling methods

Clade ages within the crocodylomorph clade Neosuchia have long been debated. Molecular and morphological studies have yielded remarkably divergent results. Despite recent advances, there has been no comprehensive relative comparison of the major time calibration methods available to estimate clade ages based on morphological data. We used four methods (cal3, extended Hedman, smoothed ghost lineage analysis (sGLA) and the fossilized birth–death model (FBD)) to date clade ages derived from a published crocodylomorph supertree and a new neosuchian phylogeny. All time‐scaling methods applied here agree on the origination of Neosuchia during the Late Triassic or Early Jurassic, and the presence of the major extant eusuchian groups (Crocodyloidea, Gavialoidea, Alligatoroidea and Caimaininae) by the end of the Late Cretaceous. The number of distinct lineages present before the K/Pg boundary is less certain, with support for two competing scenarios in which Crocodylinae, Tomistominae and Diplocynodontinae either: (1) diverged from other eusuchian lineages before the K/Pg boundary; or (2) evolved during a ‘burst’ of diversification after the K/Pg event. Cal3 and FBD proved to be the most suitable methods for time‐scaling phylogenetic trees dominated by fossil taxa. Extended Hedman estimates are substantially older than the others, with larger standard deviations and a strong sensitivity to taxon sampling and topological changes; sGLA has similar problems. We conclude that a detailed understanding of phylogenetic relationships, tree reconstruction methods, and good taxonomic coverage (in particular the inclusion of the oldest taxon in each clade) is essential when evaluating the results of such dating analyses.

O R I G I N A T I O N time estimation is essential when determining the tempo of clade evolution, and serves as the basis for a range of downstream studies, such as correlating clade origination with climatic or palaeogeographical changes (e.g. Springer et al. 2011;Hastings et al. 2013;Patiño & Vanderpoorten 2015). Significant problems with clade age estimation can occur when the datasets in question are based largely, or solely, on either fossil or extant taxa. When fossil taxa dominate a dataset, the incorporation of no (or limited) molecular sequences can ignite major debates over tree dating and clade origination times (e.g. for placental mammals: O'Leary et al. 2013;dos Reis et al. 2014). Such morphology-dominated datasets are often difficult or impossible to analyse using node-based molecular clock methods, leaving palaeobiologists largely dependent on a literal reading of the fossil record (Turner et al. 2017). Over the past three decades, however, a number of methods suitable for estimating clade ages in morphology-based phylogenies have been developed. In some cases, these were originally created for use with molecular phylogenies, but subsequently adapted for morphological datasets. These include methods using Bayesian inference, such as the fossilized birth-death model (FBD) (Stadler 2010;Heath et al. 2014). The FBD is one of the models that allows for the incorporation of fossils as tips in phylogenetic trees and randomises both speciation and extinction probabilities over time. It was further refined to allow more flexibility for fossil assignments and model parameters by Gavryushkina et al. (2014). Although using these serially sampled birth-death models can lead to long ghost lineages (Turner et al. 2017), incorporating birth, death and sampling rates renders more accurate results than those yielded by models without these rates, because the latter tend to produce estimates that are too old (Matzke & Wright 2016).
In addition to Bayesian approaches, other methods for origination time estimation have been proposed specifically for use with phylogenies with high proportions of fossil taxa. In contrast to node-dating approaches, several a posteriori dating methods use the ages assigned to tips on a pre-defined phylogeny. The first methods to apply this approach were developed in the 1990s, as extensions and formalizations of the traditional palaeobiological minimum age estimation approach. Crucially, however, these took phylogenetic topology into account and created the concept of 'ghost ranges' for implied but missing portions of the fossil record (Norell et al. 1992;Smith 1994). These basic algorithms used the ages of taxa bracketing a node to assign a minimum age to that node. However, this can lead to the inference of long ghost ranges (Turner et al. 2017) and branches of zero length (Hunt & Carrano 2010). Several methods for scaling these branches have now been developed (e.g. Ruta et al. 2006;Brusatte et al. 2008;Laurin et al. 2009;Brusatte 2011) and these algorithms are still in use in several stratigraphic time-scaling approaches (e.g. Bell & Lloyd 2014). Following the terminology of Turner et al. (2017), this set of methods will be referred to as 'ghost lineage analysis' (GLA) henceforth, and the 'equal' option (for distributing a cluster of zero-length branches, see below) as 'smoothed ghost lineage analysis ' (sGLA;Turner et al. (2017)). These algorithms can be problematic because they do not fully take into account the stochastic uncertainty of node ages that is contained within the stratigraphic record, instead incorporating only a portion of uncertainty related to occurrence time. Newer methods for stochastic age estimation have been introduced, including cal3 (Bapst 2013) and, most recently, 'extended Hedman' (EH) (Lloyd et al. 2016), based on the extension and modification of a Bayesian method to date single nodes proposed by Hedman (2010).
A comparison of these four methods (FBD, sGLA, Cal3 and EH) has been strongly encouraged by previous authors (Bapst et al. 2016). However, although earlier works have applied two or more of these methods (Bapst & Hopkins 2017;Puttick et al. 2017;Soul & Friedman 2017;Benson et al. 2018;Allen et al. 2019;Godoy et al. 2019), to date there has been no simultaneous or indepth comparison of all four.
Neosuchia is a key vertebrate group with a high proportion of extinct species (Sereno et al. 2001;Andrade et al. 2011). This is a crocodylomorph clade that encompasses all extant species (alligators, caimans, crocodiles, gharials: <30 spp). and numerous extinct taxa (>300 spp) (Brochu 2003;Grigg et al. 2015;Mannion et al. 2019). In addition, Neosuchia contains a large number of extinct crocodylian crown-group taxa (74% of the neosuchian species in our dataset), a data source usually excluded in molecular analyses. Since the earliest clade age estimations for Neosuchia (Brochu 2000), numerous molecular and morphological studies using different dating algorithms (including node age estimations) have been published (Table 1). For most groups, there is a stark difference in the node age estimates obtained from molecular and morphological methods, with molecular estimates tending to be younger (Brochu 2000), as exemplified by the dates for Crocodylinae (Table 1). These younger molecular estimates can be explained by all terminal taxa being extant, an insufficient number of calibration points (e.g. only one calibration point within Neosuchia: Roos et al. 2007) and rather young age restrictions for older clades (e.g. 90 and 100 million years ago (Ma) for the origination of Crocodylia: Oaks 2011). Furthermore, the different placements of Gavialis and Tomistoma within competing phylogenies, and the absence of entirely extinct clades (e.g. Diplocynodontinae) and older extinct species of extant lineages (e.g. non-extant crocodyline species) from molecular studies, tend to bias clade age estimates towards younger dates.
One previous study used a GLA method, cal3 and EH to date crocodylomorph trees in order to evaluate the fit of different models of body size evolution, but because of their similarity to FBD results, these additional date estimates were not discussed further (Godoy et al. 2019). So far, no study has attempted to date the origination times of all major neosuchian clades comprehensively using different methods.
The main goals of this study are: (1) to obtain age estimates of origination times of key clades in neosuchian evolution, which can further enable the identification and testing of macroevolutionary patterns, such as assessing the clade's response to the K/Pg mass extinction; (2) to apply the four different time-scaling methods to the same set of topologies in order to compare the results of newer, stochastic time-scaling methods (cal3, EH) with more established ones; and (3) to assess the relative differences between these four methods (e.g. non-random patterns in the distributions of node age estimations) and their sensitivity to taxon sampling and topological variance. To do this, we present a new neosuchian phylogeny that updates and expands upon previous work (Groh et al. 2020). We then apply the four major a posteriori methods (cal3, EH, sGLA (one of the GLA methods) and FBD) to the resulting evolutionary trees. For the purpose of comparing our results with the age estimates based on previously published neosuchian phylogenies, we also apply the same four timescaling methods to the crocodylomorph supertree presented by Stockdale & Benton (2021). The node age estimates are then compared statistically, and we discuss them in the context of previous work on neosuchian evolution.

Phylogenetic trees
The new phylogenetic topologies used in this study were derived from analysis of a modified version of the dataset presented by Groh et al. (2020). Complete information on the original phylogeny, including details of character revision and construction according to the schemes laid out in Sereno (2007) and Brazeau (2011) (but see Sookias (2020) for an argument against splitting composite multistate characters) can be found in the aforementioned publication and are not repeated herein. Modifications to the Groh et al. (2020) dataset included the incorporation of 12 additional taxa based on their published descriptions, which were selected to improve the sampling of previously under-sampled groups and to ensure that the oldest-known taxon from each of the major clades was represented (Table 2). Furthermore, we revised the character state scores for 11 species based on additional fossil material (Groh et al. 2022, phylogenetics) We deleted six characters and added two new ones (see Groh et al. 2022 (phylogenetics) for detailed changes to the character list). The final dataset comprised 565 characters (82 continuous and 483 discrete) for 122 taxa (full character list in Groh et al. 2022, phylogenetics). We used parsimony analyses, following the protocol laid out in Groh et al. (2020) (full TNT protocol in Groh et al. 2022. In brief, two versions of the dataset were analysed: the full version with continuous and discrete characters; and a discrete version where all continuous characters were rediscretized (i.e. scored in discrete character states of '0' or '1' instead of ratios based on measurements). Character weights were multiplied by 1000 to adjust for the different weighting of continuous characters (see Groh et al. (2020) for details). All phylogenetic analyses were carried out in TNT v.1.5 (Goloboff & Catalano 2016), using New Technology Search with extended implied weighting under k = 3 (Goloboff 2014) followed by a traditional search. The analyses were not constrained for the known molecular topologies, because doing so clustered all the long-snouted taxa together in a stratigraphically unlikely clade. Additionally, we did not make use of Bayesian phylogenetic methods, because previous Bayesian analyses of the same dataset yielded trees  with very low resolution that would be unsuitable for clade age estimation (Groh et al. 2020). After the initial tree analyses, four taxa were identified in stratigraphically unlikely positions that differed in their phylogenetic affinities from previous studies. These are: Sarcosuchus imperator (a tethysuchian: Sereno et al. 2001;Hastings et al. 2010Hastings et al. , 2014 and Arambourgia gaudryi (an alligatorid: Brochu 1999; Martin & Lauprasert 2010) that were found in positions as sister taxa to all other neosuchians because of an unusual combination of characters and poor preservation in parts of the examined material; Rhamphosuchus crassidens (a Miocene tomistomine: Piras 2007) placed as an early diverging gavialoid because of the fragmentary nature of the specimen; and Susisuchus jaguaribensis (a susisuchid or early neosuchian: Fortier & Schultz 2009; Turner & Pritchard 2015) placed as an early crocodyloid because of the low number of characters that could be scored (18% of all characters including no continuous ones) and the similarity of several skull table characters with those of early crocodyloids. To test for the influence of such taxa in 'rogue' placements, the subsequent dating analyses were carried out on two tree sets: the complete trees and those with the four rogue taxa in question removed from the topology.
To compare the results derived from our new phylogeny to those based on previously published neosuchian topologies, we also analysed the most recent crocodylomorph supertree, constructed by Stockdale & Benton (2021). This topology was pruned so that it contained neosuchian taxa only, matching the same five outgroup taxa used in Groh et al. (2020). We used two versions of this trimmed tree in order to account for the different placement of Thalattosuchia in other neosuchian phylogenies (e.g. Wilberg 2015): one with Thalattosuchia included within Neosuchia, as the sister group to Tethysuchia; and one with Thalattosuchia placed outside Neosuchia. The resulting trees contained 171 and 148 neosuchian taxa, respectively.

Stem-based and node-based definitions and their influence on clade age estimation
One issue that has received little prior attention is the potential for creating artefactual discrepancies in proposed origination dates, because of the way in which phylogenetic clade names are defined. Phylogenetic definitions are often either stem-based or node-based, and our own definitions are a mix of both (Groh et al. 2022, phylogenetics). In this paper, we usually refer to 'origination times', 'node age estimates' or 'clade age estimates'. These describe the age estimates for the node that represents the most recent common ancestor (MRCA) of the currentlysampled species within a given clade X, irrespective of whether its definition is stem-based or node-based. Wherever we describe the timing of divergence of a clade X from its sister clade Y, we use the term 'divergence time' (Fig. 1)

Taxon ages
Terminal taxon ages were taken from the Paleobiology Database (PBDB) (https://paleobiodb.org; accessed 15 May 2020) and adjusted for the more recently revised ages established in the International Chronostratigraphic (ICS) chart (Cohen et al. 2013;v.2020/03 2022, timescaling). The Holocene was used as the time bin for extant taxa that currently lack a fossil record. A singleton is defined herein as a taxon that includes one or multiple specimens from a single stratigraphic horizon. The presence of such singletons can falsely limit age estimates to smaller ranges and lead to inaccurate rate estimations in methods such as cal3 (Foote 1997;Bapst 2014). The modified Groh et al. and Stockdale & Benton (2021) datasets contain 31% and 41% singletons, respectively. Thus, the ages (e.g. expressed as presence during a given geological stage) for singleton taxa do not correspond to genuine stratigraphic ranges because they are point occurrences in time, so any inferred age range (based on the total age range of their parent rock unit) represents uncertainty in their fossil record. For the purpose of this study, non-singleton taxa were treated as though their stratigraphic ranges were fully known. The difference between singletons and non-singleton taxa was addressed in the following way in our analyses: 1. For the time bins used in cal3 and sGLA, singleton species were assigned time bins that covered their full potential stratigraphic range, because both methods choose the ages randomly from the assigned ranges. Non-singleton taxa were assigned the time bins that represented the uncertainty surrounding their first appearance (e.g. a species with occurrences in early and late Maastrichtian deposits would be assigned an early Maastrichtian time bin, because it cannot have originated later than the early Maastrichtian). For the rates calculation in cal3, non-singleton taxa were given full time bins covering their entire stratigraphic ranges. 2. For FBD, the midpoint age of a taxon's potential stratigraphic range was used, because use of uniform distributions for ages (such as in Lee & Yates 2018) led to younger species being dated tens of millions of years older than their actual occurrence ranges. 3. For EH, the median time bin of a taxon's potential stratigraphic range was used. In the case of a stratigraphic range spanning an even number of time bins, the later of the two median time bins was used. A full list of taxon ages is available in Groh et al. (2022, timescaling).

Origination time estimation methods
To enable comparisons between the different approaches to tree dating, four different methods were applied to each of the phylogenetic trees produced by the protocols outlined above. All analysis scripts are available in Groh et al. (2022, timescaling).

Cal3
Cal3 is an a posteriori time-scaling method first described by Bapst (2013), and further elaborated in Bapst (2014). Instead of carrying out a single clade age estimation, cal3 employs repeated stochastic sampling of node ages from a given distribution. To create this distribution of possible node ages, three different rates were defined a priori: branching (= speciation) rate, extinction rate and sampling rate, using likelihood functions (Foote 1997). This prior rate estimation is problematic without using an already time-scaled phylogeny as a basis. In addition, it is influenced by the presence of singletons, because these can artificially increase the estimates of branching and extinction rates (which are usually assumed to be equal) (Bapst 2013(Bapst , 2014. However, only 31% or 41% of the taxa in our datasets are singletons (see above) which compares favourably with those in which the percentage is higher (pterosaurs; 71.6%; Andres 2012) or unknown (dinosaurs; Lloyd et al. 2016); this problem should therefore be less severe when working with the neosuchian fossil record.
Cal3 and the algorithms used to estimate its rates are implemented in the R package paleotree (Bapst 2012) and were executed in R v.4.0.4 (R Core Team 2013). To estimate the branching, extinction and sampling rates, the function make_durationFreqDisc was applied; this constructs models of the observed taxon duration frequencies in given discrete time intervals using the likelihood functions for rate estimation by Foote (1997) as a basis. Following Bapst's recommendations written in the annotated cal3 code, artificial time intervals were generated, dividing each time period from the ICS chart into smaller, equal-length, time bins of 2-3 million years (myr). Rate estimation using these artificial time bins yielded the following rates (in species per million years) that were employed in the subsequent dating analyses: (1) branching rate = 0.1727298; (2) extinction rate = 0.1727298; and (3) sampling rate = 0.2618742. For clade age estimation, the original ICS time units were used.
Each topology was analysed 1000 times to obtain a distribution of different node age estimates. Node ages were extracted using the functions mrca and branching. times in the R package ape (Paradis & Schliep 2019).

Extended Hedman
Whole tree extension of the Hedman algorithm ('extended Hedman' or EH) is an a posteriori time-scaling method proposed by Lloyd et al. (2016), based on a Bayesian method for the dating of single nodes first described by Hedman (2010). This approach has been modified and adapted for dating all nodes within a tree by Lloyd et al. (2016) and makes use of the ages of successive outgroups to estimate posterior distributions for node ages. As such, additional outgroup ages have to be added to enable the dating of all nodes in the target phylogenetic trees.
We implemented a modified version of the EH script supplied by Lloyd et al. (2016) in R v.4.0.4. (R Core Team 2013 to sample all node ages (instead of a selected few) and time-scale a tree topology with the mean ages (complete script in Groh et al. 2022, timescaling). Ten new outgroup ages were added to each tree (Table 3). Replacing these outgroup ages with older or younger ages did not lead to significant differences in the dates obtained (Groh et al. 2022, timescaling). The time bins in the occurrence file corresponded to official ICS time units (Cohen et al. 2013) and the actual number of occurrences for each species in every time bin (data from the PBDB) was used. Dates for each node were estimated 1000 times. Ages for specific nodes were extracted using mrca and branching.times as outlined above.
Smoothed ghost lineage analysis. 'Ghost-Lineage-Analysis' (GLA) describes the simplest type of a posteriori clade age estimation analysis. Initially, each node is assigned a minimum age equal to that of the oldest member within the clade (Norell et al. 1992). However, the presence of several closely related fossil taxa of approximately the same age, can lead to the introduction of zero length branches and multiple extremely short branches when time scaling is applied. Therefore, temporally 'smoothed' GLA (sGLA) makes use of a modified version of an algorithm originally described by Ruta et al. (2006) wherein these zero length branches are assigned lengths by distributing them evenly across the time between the nodes with which they are associated. This algorithm was implemented by Brusatte (2011) andLloyd et al. (2012). Although sGLA is influenced by clustered occurrences over time (Turner et al. 2017) and potentially biased towards minimum age estimates (Lloyd et al. 2016), it has been applied widely because it is implemented in the R packages paleotree (Bapst 2012) and strap (Bell & Lloyd 2014). In our implementation of sGLA, the time bins used were the same as in the cal3 analysis. The 'equal' option of bin_timePaleoPhy was enabled in order to produce sGLA dated trees, rather than simple GLA ones. sGLA was run with three different vartime variables (= minimum time by which the root age is adjusted): vartime = 1, 5 and 10 for the first analysis. The results proved similar for all three, with notable differences in the obtained ages (all <5% of the total age) observed only in the five nodes closest to the root (see Groh et al. 2022, timescaling). Therefore, vartime = 1 was used for the analyses evaluated here. Each analysis was looped to carry out 1000 estimates for each most parsimonious tree (MPT). ape (Paradis & Schliep 2019) was used to extract node ages. ) and a number of other software packages that implement Bayesian methods. It allows for the direct incorporation of fossils into Bayesian age estimation models, randomizing both speciation and extinction probabilities over time. Depending on the tree, the performance of FBD is influenced little by additional node calibrations (Matzke & Wright 2016) but the estimation of the root age has been shown to increase in accuracy as more fossils are incorporated close to the root (P€ uschel et al. 2020). To keep the results of this method comparable to the three other time scaling methods used here, nodes were not additionally constrained, and no tree building was carried out in BEAST2, restricting the analysis to the previously supplied MPTs. Log normal priors were used throughout the analysis, with 200.3 as offset for the original tree (corresponding to the age of the oldest taxon) and M = 21.0, S = 1.0. Diversification rate was set to 1.0, turnover and sampling rates to 0.5, and Rho was 0.7. Each MPT obtained during the phylogenetic analysis was inserted separately into the xml file T A B L E 3 . Ages of additional outgroup taxa for use with extended Hedman dating in this study, based on the phylogenetic relationships in Nesbitt (2011)

Node age comparisons
There are two main reasons why we wished to compare node age estimates: (1) to understand the implications for neosuchian evolutionary history; and (2) to determine if there were any patterns in the results that the different methods delivered relative to each other and across different tree topologies. For the first aim, we determined the node age estimates for a series of 23 named nodes that potentially represent the origination of major clades within Neosuchia. Ages were extracted after all dating analyses, using the package ape ( (2003). Detailed definitions and the position of clades in our phylogenies can be found in Figures 2 and 3, and Groh et al. (2022, phylogenetics). For the second aim (method comparison) we evaluated all nodes, not just the 23 named nodes mentioned above. We determined whether there was a statistically significant difference between the mean node age estimates, yielded by pairs of the dating methods, for a given node in a given tree topology. Thus, for example, for node X in topology Y, cal3 and EH would each provide 1000 node age estimates, and the means of these populations can then be compared using two-tailed t-tests (carried out here in Microsoft Excel 2021). This resulted in six two-tailed t-tests (cal3 compared to FBD, cal3 to EH, cal3 to sGLA, EH to FBD, EH to sGLA, FBD to sGLA) per node. If the p-value was significant (<0.05), we noted which method provided the younger and which method the older estimate.

Phylogenetic analysis
Our new phylogenetic analysis is based on the modified version of the Groh et al.

Comparison of different clade age estimation methods
The mean node age estimates and their standard deviations for each key clade highlight the differences between the four time-scaling methods used (Fig. 4). The results across all topologies from both sets of trees are very similar with only a few minor exceptions, which are detailed below.
Across all topologies, cal3 consistently delivers the statistically significantly youngest mean node age estimates for almost all nodes, compared to other methods (p < 0.05 for at least 99% of nodes compared to sGLA, 96.6% of nodes compared to EH and at least 66.7% compared to FBD). The only exceptions are nodes closer to the youngest tips, where FBD delivers significantly younger mean node age estimates than cal3. The oldest mean node age estimates are usually provided by EH, with at least 82.1% of all nodes having significantly older mean node age estimates than those produced by any other method. This pattern is only disrupted close to the root of the trees where sGLA and FBD often deliver the oldest ages. FBD and sGLA give relatively similar mean node age estimates, with 55-66% of mean node age estimates obtained by sGLA significantly older than those obtained by FBD and 4.3% of nodes showing no statistical difference in mean ages between the two methods (Groh et al. 2022, timescaling). These patterns apply to both the SB21 and G21 topologies.   The standard deviations of estimated node ages are usually low for cal3 and sGLA (always 10% or less than the total age estimate), intermediate for FBD (around 10%) and high for EH (usually above 10%). In addition, cal3 shows the lowest number of nodes that are influenced by the removal of the four rogue taxa in the G21 topologies, whereas EH, FBD and sGLA all experience great impact (Groh et al. 2022, timescaling). In the G21 topologies, removing the four rogue taxa leads to differences in mean node age estimates in five out of the twenty named clades, whereas mean ages were different for all named clades for FBD, and all but one (EH) and two named clades (sGLA) in the other two methods.
The variation in estimated clade ages increases with distance from the root, with a difference of 100% in several cases (e.g. Crocodylidae with mean node age estimates ranging from 62 Ma with cal3 to 127 Ma with EH in the G21 trees). In contrast, the difference in mean node age estimates close to the root is only 20% (e.g. Neosuchia with estimates ranging from 200.8 Ma (cal3) to 220.9 Ma (FBD)). A similar trend can be observed in the mean age estimates of unnamed nodes, with younger nodes closer to the tip showing a difference of up to 500% between youngest and oldest ages in several cases, whereas the difference is closer to 20-50% for older nodes closer to the root (Groh et al. 2022, timescaling).
For several nodes, the sGLA estimates display 'double peaks'. These bimodal distributions of age estimates do not denote differences in dating between tree topologies but occur for most crocodylian nodes in every tree, even when only a single topology is dated. We examine this anomaly in detail in the Discussion.

Origination times for key clades
Mean node age estimates for the origin of Neosuchia vary between 195 Ma (cal3) and 223 Ma (FBD) across all methods and both the G21 and SB21 phylogenies, with the majority of estimates between 200 and 205 Ma (Fig.  4). This clade age changes in only a minor way depending on whether Thalattosuchia is included in the definition of Neosuchia for both datasets (Table 4). Thalattosuchia is estimated to have originated between 193 Ma (FBD) and 203 Ma (EH) in the SB21 phylogeny.
Tethysuchia appears as two paraphyletic assemblages in both Groh et al. (2020) and G21 (Fig. 2) Summary of the different origination ages for the major neosuchian groups, as found by the four different methods for the two different topologies considered herein, excluding both rogue taxa and Thalattosuchia.

Neosuchian phylogeny
Modifying the dataset of Groh et al. . Our G21 results also identify Tethysuchia as paraphyletic, with Elosuchus and the newly added Meridiosaurus forming a sister group to Dyrosauridae and the remaining neosuchians. Elosuchus and Meridiosaurus are grouped together on the basis of several continuous characters, as well as the presence of two discrete apomorphies (a dorsally broad postorbital bar (char. 211) and long squamosal prongs (char. 231)). Our Dyrosauridae, which is similar in species composition to the group supported by Martin et al. (2016), is not clearly resolved, and is characterized by a five taxon polytomy in the consensus tree in both analyses. Susisuchus anatoceps is no longer placed as part of Goniopholididae as it was in Groh et al. (2020), but instead is found as a noneusuchian neosuchian that occupies varying positions, because of the addition of atoposaurids, paralligatoroids and Calsoyasuchus to the dataset. Goniopholididae now includes Vectisuchus and both species of Sunosuchus, because of the absence of Susisuchus as part of the clade. Goniopholidid monophyly is based on the quadratojugaljugal suture lying at the posterior angle of the infratemporal fenestra (char. 264); the vomer being exposed on the palate between the palatines (char. 346); the flat lateroventral surface of the dentary anterior to the external mandibular fenestra (char. 401); and thick and vertical osteoderm margins (char. 551).
Bernissartia fagesii forms the closest sister group to Eusuchia, a similar result to that found by previous analyses (Sweetman et al. 2014). Aegyptosuchidae has previously been placed as the sister group to Brevirostres (Delfino et al. 2008)

Comparison of methodological results
Cal3 generally delivered the youngest age estimates in this study, with the lowest standard deviations of all methods. These low standard deviations rarely encompassed the age ranges given by the other methods. Moreover, the deletion of taxa had a much lower impact on cal3 node ages compared to the other three methods. Indeed, cal3 estimates often remained similar despite differences in the tree topologies analysed: for the G21 dataset, mean node age estimates for the discrete and complete character set based topologies were the same for many nodes (Groh et al. 2022, timescaling); for the SB21 dataset, mean node age estimates were the same or within 2 myr of each other for all nodes (except those within Tethysuchia), irrespective of whether or not Thalattosuchia was removed from Neosuchia (Groh et al. 2022, timescaling). Knowing how the four different methods work and the assumptions they make can help us to understand the above results. Cal3 initially sets the youngest possible age for a node as the first appearance date (FAD) of the oldest member of the clade that is to be dated. It then determines the probability distribution of the potential node ages by moving from root to shallowest node from this point onwards, using pre-determined extinction, speciation and sampling rates (Bapst 2013). However, determining these rates can be more difficult in the presence of a large number of singleton taxa (Bapst 2013(Bapst , 2014. In addition, a number of branches exhibited the minimum branch length of 0.1 myr (set in accordance with the recommendations of Bapst 2014), which potentially introduces an arbitrary component to the results. Because cal3 bases its initial node age estimate on the oldest age within a clade, rather than its closest relatives, the tendency of cal3 to produce the youngest age estimates is to some extent an artefact of how the method works. However, the other methods have similar artefacts that lead, for example, to much older estimates (EH).
There are two types of taxon removal that can impact age estimates: (1) general taxon exclusion/absence; and (2) absence/exclusion of the oldest taxon in a clade. Cal3 is particularly sensitive to the second type of taxon removal, because the oldest taxon within the clade plays a critical role in determining the other node age estimates. This is exemplified by the case of Crocodyloidea in the G21 topologies, with mean node age estimates of 142.5 Ma (before the removal of Susisuchus jaguaribensis) and 74.1 Ma (after removal). Estimates can be shifted strongly towards younger node ages if the oldest known taxon is missing. However, cal3's estimation method is also the reason for its robustness concerning other changes in tree topology, such as the general removal of individual taxa (as long as they are not the oldest in a clade) and minor rearrangements of more deeply nested taxon interrelationships. Similarly, the removal of entire taxon groups, such as Thalattosuchia in the SB21 topology, has little effect on downstream ages when cal3 is applied, with differences of 100 000 years or less in the mean ages of the main nodes. If the order of taxa within a clade changes slightly, the overall initial minimum clade age tends to remain the same.
The results of EH stand in stark contrast to those of cal3. For the majority of nodes, it provides the oldest mean node age estimates, sometimes more than twice as old as those of cal3 (e.g. Gavialinae in G21). This effect is even more pronounced for the SB21 topologies, with EH mean age estimates often twice as old as the ones obtained by cal3, and up to four times as old (e.g. for Crocodylinae). In addition, EH standard deviation rates increase with distance from the root in both sets of topologies and are considerably larger, often 3-10 times higher than those of cal3, with minimum and maximum age estimates spanning a much wider range (Fig. 4). These SD ranges partially encompass the ranges of the other age estimates, overlapping with the cal3 means 33% of the time, with FBD estimates 61% of the time, and with sGLA estimates 63% of the time, for the G21 topology (Groh et al. 2022, timescaling). Although they did not report markedly younger cal3 estimates, similarly high standard deviations for EH were observed by Lloyd et al.
(2016) in their study of dinosaurs. In contrast to Lloyd et al. (2016), however, our analyses found that altering the ages of the additional outgroups had little influence on the obtained estimates.
EH estimates are based on the taxon ages within successive outgroups to the clade that is being dated. This dependence on the taxa contained within more than one clade produces a higher sensitivity to changes in tree topology, especially further from the root of the tree. These changes include both general taxon exclusion/ absence and removal of the oldest taxon in a clade. It also makes the method particularly sensitive to sampling issues and biases. If a clade is part of a well-sampled group with multiple species in each time bin, EH will produce younger age estimates that are closer to those of cal3, as observed previously (Lloyd et al. 2016). By contrast, the older the ages for the 'successive outgroups' to the dated clades, the older the obtained node ages. For the G21 topologies, this effect is demonstrated by Caimaininae (ages of 81.8 Ma and 107.3 Ma yielded by cal3 and EH, respectively) and Gavialoidea (ages of 74.2 Ma and 135.8 Ma yielded by cal3 and EH, respectively). Although poorer taxon sampling is likely to increase errors when applying any origination time estimation method, EH seems to be particularly influenced by both lower taxon coverage and uneven sampling of the fossil record.
In most cases, and for both sets of topologies, FBD supplied node ages younger than those given by EH and older than those obtained from cal3. The exceptions are the nodes closest to the root, where FBD estimates were the oldest of those provided by any method. FBD age estimates display intermediate standard deviations between the low and high values obtained by cal3 and EH, respectively. The ranges of the FBD age estimates often overlap with those from EH (although there is still a statistically significant difference in means: Fig. 4). FBD calculates node ages by using priors to determine the evolutionary rates along the entire tree. These rates are changeable throughout different time intervals, as per the skyline version of FBD (Gavryushkina et al. 2014). Although FBD does not rely on the ages of single taxa within clades for minimum node age estimates, unlike cal3 and EH, it was even more strongly influenced by the removal of the four rogue taxa than the former two methods. This includes both exclusion of the oldest taxon in a clade, and the general removal of taxa. By contrast, minor differences in tree topology cause less variation in node age estimates than they do in cal3.
FBD assumes constant speciation, extinction and sampling rates among phylogenetic branches over time and, accordingly, scales node ages even of successive older nodes, as can be seen with the estimates for the nodes close to the tree root. For example, whereas cal3, EH and sGLA suggest a diversification burst in early-branching neosuchians and closely related sister groups shortly after the end of the Triassic, FBD postulates a Late Triassic origin for Neosuchia for both sets of topologies, with a much more substantial amount of time passing between each node in the tree (Groh et al. 2022, timescaling). On average, there is less than a million years between the earliest neosuchian nodes in cal3: by contrast, the nodes in the phylogeny time-scaled with FBD are, on average, at least 2 myr apart.
The estimates obtained by sGLA display variable relationships with those produced by the other three methods. Estimates for nodes close to the root tend to be older, whereas more tipward node age estimates are younger and, at times, approach those of cal3. For both the G21 and SB21 topologies, sGLA mean node age estimates are usually similar to those from FBD, both in mean age and standard deviations. This is in contrast to previous studies, which found sGLA ages to be older than those obtained by FBD (Turner et al. 2017). As in EH, sGLA uses the next closest outgroup to the node being dated to establish a minimum node age. However, instead of taking multiple 'successive' outgroups into account, it relies only on the closest sister group in establishing the minimum age and scales up the estimate in accordance with the age of earlier nodes.
Our results show that sGLA is the method most sensitive to the general exclusion of single taxa (including the oldest one within a clade), although it does not respond as strongly to minor changes in topology as EH. Because sGLA adjusts node age estimates based on older taxa, it shares EH's greater sensitivity to sampling bias and coverage. Unlike cal3 and EH, sGLA is not per se a stochastic time-scaling algorithm, as can be seen in the anomalously bimodal frequency distribution of sGLA node age estimates. We have been unable to find an exact source for these bimodal distributions, as a detailed exploration of the mathematical background of the method lies beyond the scope of this paper.
Finally, it is often argued that greater taxon sampling (facilitated by the use of supertrees) is preferable when investigating a group's macroevolutionary history (e.g. Haeseler 2012). Here, the application of four dating methods to two independent sets of topologies with different levels of taxon sampling, allows us to examine this issue with regard to the best practice for establishing node ages. The different methods deliver roughly similar results using either the G21 MPTs or the SB21 supertrees; EH delivers the oldest mean node age estimates, and cal3 the youngest. As stated above, some methods (such as EH) are more influenced by poor taxonomic coverage than others, so the use of more taxon-rich supertrees might be regarded as preferable. However, our analyses demonstrate that the stratigraphic position of the sampled taxa is more important than overall taxonomic coverage. Specifically, the inclusion of the oldest taxon for each clade is paramount, because many methods determine node ages by using the oldest taxon within, or closely related to, the clades whose age is to be estimated. This phenomenon is illustrated by the G21 estimates for Crocodylinae, Diplocynodontinae and Eusuchia. All three groups were assigned unrealistically young clade ages (as judged from the fossil record) in the SB21 topology because of the omission of the oldest known fossils. In other cases, the greater taxon sampling of the SB21 supertree proved advantageous, delivering more stratigraphically consistent age estimates (e.g. for Dyrosauridae and Tethysuchia; see below).

Recommendations
Previous studies have advocated the use of stochastic time-scaling approaches such as cal3 and EH, over simpler methods like sGLA (Lloyd et al. 2016;Bapst & Hopkins 2017), as well as emphasizing the importance of complete taxon sampling (Lloyd et al. 2016) and the role of cal3 as a viable alternative to FBD (Bapst et al. 2016). They have also recommended the use of EH for groups with low sampling rates, and the use of cal3 for those with high sampling rates (Soul & Friedman 2017), as well as the theory underlying FBD (Bapst 2013). Other studies have used cal3, sGLA and EH to date trees for analyses of trait evolution, such as body mass in Dinosauria and Crocodylomorpha (Benson et al. 2018;Godoy et al. 2019), and extinction as a binary trait in non-teleostean actinopterygians, end-Permian tetrapods and archosauromorphs (Puttick et al. 2017;Soul & Friedman 2017;Allen et al. 2019). None of these studies, however, discussed the differences in estimated node ages themselves, only their impact on model selection and trait evolution coefficients.
After evaluating all four methods, we propose some recommendations for the best way of estimating clade ages in future studies. 1. It is clear that the various time-scaling methods produce different (often radically different) results. As such, it is dangerous to pick one method, apply it and take the results at face value, because this would not account for any sensitivity to the choice of method that the results might display. Applying only a single method could lead to an incomplete sampling of possible node ages and heavily skew estimates. If only one method is to be applied, then this choice should be justified based on the properties of the dataset (e.g. proportion of singletons, the extent of coverage of the earliest clade members, overall taxon sampling). Ideally, however, multiple dating methods should be applied (see Recommendation 5, below). If the latter leads to unhelpfully wide ranges for particular nodes of interest, then narrower 'preferred' node ages may be justified based on the results of the method(s) that appears best suited to the properties of the dataset. 2. Although taxon sampling as a whole is important, studies should pay special attention to the inclusion of the oldest putative members of the clades they wish to date. Without the inclusion of these species, mean node age estimates are often too young, even if the taxonomic coverage is otherwise good; as when supertrees are used, for example. 3. Extended Hedman generally yields older dates that often imply longer ghost ranges than the other methods, such as the 105 myr-long ghost range leading to Borealosuchus formidabilis in the G21 topology, compared with the 70 myr identified by cal3 and FBD (Groh et al. 2022, timescaling). The most significant problem is EH's sensitivity to various different issues, such as taxon absence/removal and sampling bias. It should thus only be applied to datasets with high taxon sampling for all the known clades within a tree but should be considered for use as an alternative to cal3 in the presence of numerous singletons. However, when doing so, particular attention should be paid to issues of tree topology and taxon sampling. 4. sGLA is the single most sensitive method to taxon exclusion and does not include stochastic estimates and is thus unsuitable for most analyses. In addition, the artefactually bimodal distributions for its age estimates make a correct age estimation difficult, especially in the presence of extant taxa. 5. We recommend the use of both cal3 and FBD for most datasets, especially those with incomplete taxon sampling. Care should be taken, however, when using methods such as cal3 that rely on prior estimation of speciation, extinction and sampling rates without phylogenetic information on datasets with a high proportion of singletons. A substantial amount of early neosuchian diversification took place during the Early Jurassic. All four methods agree that for the G21 trees, the lineages leading to Tethysuchia, Goniopholididae, Paralligatoridae, Atoposauridae and Eusuchia had diverged from each other by 180 Ma (Fig. 2). The clade ages based on the SB21 phylogeny agree, with the exception of the lineage leading to Paralligatoridae, which is estimated to have diverged during the Middle Jurassic (Fig. 3). With regard to stemeusuchian lineages and Eusuchia, our two preferred methods (cal3 and FBD) suggest two different evolutionary histories based on both the G21 and SB21 phylogenies (Table 4): (1) rapid divergence of the remaining noneusuchian neosuchian clades from the main lineage leading to Eusuchia during the Early Jurassic, between [190][191][192][193][194][195][196][197][198][199][200] Ma (cal3, also supported by sGLA); or (2) a slower diversification during the Late Triassic and Early Jurassic, 190-220 Ma (FBD, also supported by EH). We favour the scenario of a rapid early burst of diversification during the earliest Jurassic supported by cal3 (and sGLA) for both sets of topologies, because it coincides with the increased availability of ecological niches following the end-Triassic mass extinction (see below).

Clade ages and neosuchian evolutionary history
Previous node-based origination time estimates for a monophyletic Tethysuchia have ranged across the Early-Middle (Turner et al. 2017;Godoy et al. 2019) (cal3) to 84 Ma (FBD), eliminating the ghost range found in the G21 topology and providing a more stratigraphically consistent age estimate (Fig. 3). We favour the Middle to Late Jurassic tethysuchian origin dates based on the environmental changes discussed below and the stratigraphically more consistent Campanian dates for Dyrosauridae provided by the SB21 topologies.
Both  (Table 1). Our two favoured origination time estimation methods each support one of these possibilities for the G21 topology, yielding origination dates from 131 Ma (cal3) to 161 Ma (FBD). Similar to the pattern observed for the neosuchian origination date, the two suggested evolutionary scenarios for the main eusuchian clades are either gradual (FBD) or with several diversification bursts over a shorter time period (cal3). If Susisuchus jaguaribensis, a susisuchid or early neosuchian (Fortier & Schultz 2009;Turner & Pritchard 2015), is included as an early crocodyloid, these ages are skewed substantially towards older, Late Jurassic estimates. The mean node age estimates obtained for the SB21 topology are younger, ranging from 114 Ma (cal3) to 127 Ma (FBD). This is mainly because of the omission of the oldest known eusuchian, the Barremian Hylaeochampsa vectiana (Clark & Norell 1992), from the SB21 supertree, leading to unrealistically young estimates for the eusuchian clade age (Fig. 3). Therefore, we favour the scenario proposed by the G21 cal3 estimates, placing the origin of Eusuchia shortly after the Jurassic-Cretaceous boundary (Fig. 5).
In contrast to Eusuchia, crocodylian origination ages have mostly been estimated as Late Cretaceous (Salisbury  (Table 1). Dates based on the G21 dataset have been estimated as slightly earlier because of the inclusion of Hylaeochampsidae and Aegyptosuchidae as crocodylian clades (Fig. 2). As a result of the way the clade name definitions for Eusuchia and Crocodylia map onto the G21 topologies, these two clades are identical and so have a single ancestral node and clade age estimate. The most likely origin date for this G21 Eusuchia-Crocodylia clade is Late Jurassic or Early Cretaceous. Early Cretaceous origination time estimates for Crocodylia, however, have been strengthened recently by the discovery of Portugalosuchus azenhae, a potential crocodylian from the Cenomanian (Mateus et al. 2018), which is included in the G21 phylogeny. The SB21 topology does not include Portugalosuchus or Hylaeochampsa, and thus produces much younger and less likely mean node age estimates, of 85 Ma (cal3) to 101 Ma (FBD). In short, an Early Cretaceous origin of Crocodylia is most likely.
Similar to Crocodylia, Late Cretaceous clade ages have also been suggested for Brevirostres (Salisbury et al. 2006;Martin et al. 2010;Pu ertolas et al. 2011;Turner et al. 2017). The oldest possible specimen referable to Brevirostres (as well as Crocodyloidea, Crocodylidae and Crocodylinae) is hard to determine because of the ambiguous status of many fossils labelled as 'Crocodylus', with many relevant specimens either not belonging to Crocodylinae or being very fragmentary (Brochu 2000; Oaks 2011). For example, incomplete mandible fragments from the Cenomanian, described as Crocodylus selaslophensis (Molnar & Willis 2001), have recently been re-identified as belonging to the susisuchid Isisfordia (Hart et al. 2019; Hart 2020). As with Eusuchia, FBD and cal3 mean node age estimates differ substantially in the G21 topology, suggesting either: (1) a date of 87 Ma for Brevirostres (cal3) and a late Early to Late Cretaceous diversification of the main eusuchian groups (Gavialoidea, Hylaeochampsidae, Brevirostres); or (2) a date of 130 Ma for Brevirostres and a more gradual evolution of these clades mostly during the Early Cretaceous (FBD) (Fig. 5). In contrast to the G21 trees, the clade definition of Brevirostres means that the latter includes Gavialoidea in the SB21 topology, with clade age estimates for the group placed during the Late Cretaceous, similar to the cal3 estimates for the G21 dataset. It is difficult to prefer one of these two scenarios, because both can be linked to biotic and abiotic factors that might have promoted diversification (see below).
The diversification of Eusuchia during the Late Cretaceous included the origin of Gavialoidea, which occurred between 71 and 90 Ma, according to earlier studies  Table 1). This is in line with the mean node age estimates we obtained with our preferred methods, cal3 and FBD, for both the G21 and the SB21 datasets, with ages varying from 71 Ma (cal3) to 95 Ma (FBD). EH estimates are substantially older for both sets of topologies and Early Cretaceous in age (Table 3), a discrepancy that persists throughout the gavialoid part of the tree, as far as Gavialinae. Cal3 and FBD place the origin of Gavialinae after the K/Pg boundary, at 40 Ma and 50 Ma, respectively (Fig. 2). Therefore, a Late Cretaceous origin of Gavialoidea is likely, as is an early Palaeogene origin of Gavialinae, potentially promoted by abiotic factors (see below).
Establishing a firm minimum clade age for Crocodyloidea is difficult because of the above mentioned problem with the accurate identification of early 'Crocodylus' species, but most previous morphological studies suggest Our estimates for the G21 topology (there is no direct equivalent to Crocodyloidea in SB21) are strongly influenced by the presence of Susisuchus jaguaribensis. Removing it from the phylogeny leads to middle (100 Ma, FBD) or Late Cretaceous (74 Ma, cal3) dates for the origin of the crocodyloid clade (Figs 2, 4).
The current literature is divided on whether Crocodylidae was present before the K/Pg boundary (Table 1) Results from the G21 topology mirror this divide (Fig. 5): cal3 supports a diversification burst in Crocodylidae shortly after the K/Pg boundary, with Tomistominae, as well as Crocodylinae and multiple less diverse crocodyline groups diverging in quick succession. FBD, by contrast, favours a more gradual diversification (Fig. 5), taking place mostly during the Late Cretaceous after the origin of Crocodylidae during the middle Cretaceous (Table 4)  ). The history supported by the mean node age estimates for the SB21 topology is similar, although the estimates for Crocodylidae are younger than those provided by the G21 topology. The SB21 trees support a Maastrichtian (FBD) or early Eocene (cal3) origin for Crocodylidae and Mekosuchinae (Fig. 3). Both scenarios can be linked to potential causal abiotic factors (see below), making it difficult to prefer one over the other.
Crocodylinae has previously been estimated to have originated during the middle-late Palaeogene or early Neogene, by both morphological (Salisbury et al.  (Table 1). Our mean node age estimates for the SB21 topology mirror these dates, with 24 Ma (cal3) and 44 Ma (FBD) found for Crocodylinae. Our origination time estimates for the G21 topology are substantially older, with 61 Ma (cal3) and 79 Ma (FBD). This disagreement is mainly because of topological differences, namely: (1) the placement of the Paleocene family Planocraniidae as the closest sister group to Crocodylidae in G21; and (2) the definition of Crocodylinae and whether it should include the early .1 Ma) Crocodylus megarhinus. According to Brochu (2000), C. megarhinus is an early-branching crocodyline, but it is not included in the SB21 topology, leading to the younger ages supported by the latter. The clade definitions used here (Groh et al. 2022, phylogenetics) are similar to those in Brochu (2003) and Gatesy et al. (2004): Crocodylinae encompasses all taxa that do not belong to Mekosuchinae (which is not included in the G21 dataset) and that are more closely related to extant taxa such as C. niloticus than to Tomistominae. Since C. megarhinus is included in the G21 dataset, this leads to older age estimates than for the SB21 topology. Thus, we regard the G21 node age estimates to be more reliable in this instance, and therefore support a diversification burst of early crocodyline species during the Paleocene, shortly after the K/Pg mass extinction, with most modern Crocodylus species and their closest relatives emerging during the Miocene.
The two different scenarios for the timing of crocodylid evolution suggested by cal3 and FBD apply also to Tomistominae in the G21 topology (this clade is represented only by two recent taxa in the SB21 topology). Previous morphological estimates have placed its origin in the Paleocene (Salisbury et al. 2006;Pu ertolas et al. 2011;Godoy et al. 2019). Unlike the latter trees, the G21 phylogenies were not constrained by molecular hypotheses, thus our estimates for the origin of Tomistoma are older than those in molecular trees (e.g. 16-30 Ma;Oaks 2011). All four of our time-scaling methods placed the most recent common ancestor of tomistomines in the early Palaeogene, from 49 Ma (cal3) to 62 Ma (FBD). However, the timing of the divergence of Tomistominae from the crocodyline lineage is dated differently, either during the Late Cretaceous (FBD), or the Paleocene (cal3), tied to the early Palaeogene diversification burst in Eusuchia (Fig. 5).
The third of the major eusuchian clades to originate during the Late Cretaceous was Alligatoroidea (Salisbury et al. 2006;Pu ertolas et al. 2011;Lee & Yates 2018; Table 1). After the removal of Susisuchus jaguaribensis from the G21 topology, the mean node age estimates from both the latter and the SB21 topologies suggest an origination of alligatoroids either between 95 and106 Ma (FBD) or later at around [84][85][86], close to the age of the oldest known alligatoroid, Albertochampsa langstoni (Erickson 1972).
Diplocynodontinae, usually regarded as an alligatoroid clade (Brochu 1999;Delfino & Smith 2012), is placed as the sister clade to Crocodyloidea+Alligatoroidea in the G21 phylogenies, although it is still included in Brevirostres as defined herein. The earliest known putative diplocynodontine specimen, Diplocynodon sp., is from the Late Cretaceous (83. 8-72.1 Ma) (Grandstaff et al. 1992). However, Grandstaff et al. (1992) did not provide any descriptions of the fossils, nor is the specimen described elsewhere in the literature, making its status as the only known Cretaceous diplocynodontine doubtful. The oldest confirmed diplocynodontine material is of early Eocene, and potentially late Paleocene, age (Martin et al. 2014a;Rio et al. 2020). The lineage leading to Diplocynodontinae diverged from the brevirostrine lineage during the Early to middle Cretaceous. The origination date for Diplocynodontinae itself is estimated at either Late Cretaceous (86 Ma, FBD) or Paleocene (61 Ma,cal3) in the G21 dataset. The mean node age estimates from the SB21 phylogeny are substantially younger than the G21 FBD estimate, both being Palaeogene in age (49 Ma (cal3) and 62 Ma (FBD)) (Fig. 4). These SB21 node ages are here regarded as unlikely because of the age of the earliest confirmed fossil material of this lineage. Thus, we support a Palaeogene origin for Diplocynodontinae.
Origination times for Globidonta differ strongly between molecular and morphological studies, with estimates around the K/Pg boundary derived from the former (Roos et al. 2007;Oaks 2011), and Late Cretaceous provided by the latter (Salisbury et al. 2006;Pu ertolas et al. 2011;Lee & Yates 2018;Mateus et al. 2018;Godoy et al. 2019) (Table 1). The oldest confirmed globidontan taxa are from the Late Cretaceous, Stangerochampsa mccabei (72. 1-66.0 Ma;Wu et al. 1996) and Albertochampsa langstoni (83. 6-72.1 Ma;Erickson 1972). There are potential remains of Brachychampsa sp. that are apparently the same age (83. 8-72.1 Ma) (Storer 1993;Nessov 1995), but this cannot be verified from the available descriptions. Our cal3 dates (82 Ma and 83 Ma in SB21 and G21, respectively) are very similar to those obtained by morphological studies and the age of the earliest globidontans, whereas FBD mean node age estimates are slightly older (87 Ma and 102 Ma). As seen with the ages obtained for the other eusuchian lineages, cal3 suggests a Late Cretaceous diversification burst similar to that of Notosuchia (Sol orzano et al. 2019), whereas FBD supports a more gradual scenario of alligatoroid evolution (Fig. 5) for both sets of topologies. Thus, despite the 20 myr age range provided by morphology, it seems that molecular age estimates for globidontan origin are too young.
By contrast to the G21 trees, the SB21 topologies contain both Alligatoridae and Alligatorinae. Mean node age estimates from our preferred methods corroborate their suggested presence before the K/Pg boundary (Brochu et al. 2012), with estimated clade ages of 68-83 Ma for Alligatoridae and 66-74 Ma for Alligatorinae with cal3, respectively. This overlaps with node age estimates within Globidonta in the G21 topologies, making a Late Cretaceous origin of Alligatoridae and early Palaeogene one for Alligatorinae likely.
Previous molecular and morphological studies have both placed the origin of Caimaninae after the K/Pg boundary (Salisbury et al. 2006;Roos et al. 2007;Oaks 2011;Pu ertolas et al. 2011;Lee & Yates 2018), anywhere between 20 and 58 Ma (Table 1). The oldest caimanine taxa are represented by Eocaiman palaeocenicus from the late Paleocene (Bona 2007) and a number of other remains from the early Paleocene (Brochu 2012;Pinheiro et al. 2013;Bona et al. 2018). By contrast, our estimates from both topologies are mostly in the Late Cretaceous, 82 Ma (cal3, G21 topology) to 78 Ma (FBD, SB21 topology) and 85 Ma (FBD, G21 topology), because of the altered topology and inclusion of several old caimanine species. Only the cal3 estimates for the SB21 topology are marginally younger than the K/Pg boundary, 66 Ma (Table 3), making a Late Cretaceous origin and subsequent early Palaeogene diversification for the clade probable.

Biotic and abiotic factors influencing neosuchian evolutionary history
The end-Triassic mass extinction included the loss of tetrapod groups such as Phytosauria and Rauisuchia that had at least some ecological similarities to crocodylomorphs (Toljagi c & Butler 2013;Bronzati et al. 2015;Godoy et al. 2019). The proposed Early Jurassic origin of Neosuchia, and its diversification into Goniopholididae, Paralligatoridae, Atoposauridae, and the lineage leading to Eusuchia, might therefore represent an opportunistic replacement with the clade filling niches left vacant by the mass extinction. This radiation can be detected in both increased diversification rates and disparity among Early Jurassic crocodylomorphs (Toljagi c & Butler 2013; Bronzati et al. 2015;Godoy et al. 2019). The end of the Early Jurassic marked the onset of warmer climates (Sellwood & Valdes 2008), with increases in sea surface temperature and sea level, both of which have been linked to the diversification of marine and non-marine crocodylomorphs (Martin et al. 2014b;Tennant et al. 2016a). This supports the scenario of a Middle to Late Jurassic origin of Tethysuchia and the diversification of a range of younger goniopholidid species.
During the Early Cretaceous, numerous environmental perturbations caused several tetrapod lineages with niches similar to those of early eusuchians to go extinct, such as some non-marine turtles, mammals and non-marine crocodylomorphs (Tennant et al. 2016a). This was followed by a short period of elevated sea level and temperatures (Carvalho et al. 2010;Tennant et al. 2016a, b). The extinction of the above-mentioned lineages, coupled with these abiotic factors, could have led to an increase in the number of the ecological niches available for early eusuchians during the Early Cretaceous. Increased temperatures are linked to higher evolutionary and diversification rates in Crocodylomorpha (Markwick 1998;Tennant et al. 2016b;Sol orzano et al. 2019;Stubbs et al. 2021), potentially facilitating eusuchian diversification and supporting an Early Cretaceous origin of Eusuchia. Moreover, hotter  The K/Pg mass extinction witnessed severe reductions in the diversity of many tetrapod lineages, including the loss of all non-avian dinosaurs. This increased availability of niches and resources led to diversification bursts in several vertebrate clades, such as neornithine birds and placental mammals (de Celis et al. 2019). The early Palaeogene was also characterized by the Paleocene-Eocene Thermal Maximum (PETM) (R€ ohl et al. 2007) and the Early Eocene Climatic Optimum (EECO) (Lauretano et al. 2015). This combination of opportunistic replacements immediately after the K/Pg extinction, and higher temperatures during the Paleocene-Eocene, potentially led to increased speciation rates among eusuchians. These conditions promoted the evolution of Gavialinae, Tomistominae, Diplocynodontinae, Caimaninae, Mekosuchinae and early Crocodylinae. In addition, the diversification of Caimaninae was potentially linked to range expansion facilitated by an ephemeral North-South

CONCLUSION
We have used an updated phylogeny (based on Groh et al. 2020) and a supertree (Stockdale & Benton 2021) to examine the origination dates of major neosuchian clades and evaluate the behaviour of newer stochastic timescaling methods (cal3 and extended Hedman (EH)) in the presence of different types of data compared to fossilized birth-death (FBD) and smoothed ghost lineage analysis (sGLA). All four time-scaling methods across both sets of topologies agreed that the lineages leading to the main neosuchian clades (Tethysuchia, Goniopholididae, Atoposauridae, Paralligatoridae, Eusuchia) had diverged by the Early Jurassic. Similarly, the three main crocodylian superfamilies (Gavialoidea, Crocodyloidea, Alligatoroidea) appeared before the K/Pg boundary, possibly in response to the high temperatures of the Cretaceous Thermal Maximum, which occurred in the Turonian. Two different scenarios of eusuchian evolution are supported in both sets of topologies: a rapid diversification burst of Eusuchia, and in particular, Crocodylidae, after the K/Pg boundary, or a more gradual evolution of the different crocodylid lineages (including Tomistominae) beginning during the Late Cretaceous. Our results also support the inferred presence of Caimaninae before the K/Pg boundary, contrary to recent molecular dating estimates.
In general, stochastic time-scaling methods are not more reliable than others: cal3 is the most robust method with regards to its relative insensitivity to solitary taxon exclusions and minor changes in topology; however, EH is affected by these issues to a similar degree as sGLA and FBD. For the majority of nodes, EH provides the oldest age estimates (together with larger standard deviations), whereas cal3 provides the youngest ages and narrowest standard deviations, in part because it also allows zerolength branches. In addition, EH and sGLA are particularly influenced by sampling biases and taxon coverage in the dated topologies. Inclusion of the oldest-known taxon of each major taxonomic group in the tree is paramount in order to obtain more realistic estimates for all four methods, more so than full taxon coverage. As such, the use of supertrees instead of MPTs cannot be recommended without caveats, because supertrees rely on their source trees containing the oldest-known taxon of each clade. If both MPTs and supertrees contain the oldestknown taxa, the increased number of taxon sampling in supertrees does not per se provide better age estimates, unless sampling covers groups not represented in the MPTs. In addition, supertrees can potentially suffer from a number of other problems, such as non-independence of source trees, novel relationships not represented in any of the input MPTs, and longer resulting trees than those based on supermatrices (Gatesy Based on our results, we recommend the combined use of cal3 and FBD as preferred origination time estimation methods. Given the artefactual bimodal age estimate distributions of sGLA and the relatively arbitrary, linear criteria of scaling branch lengths, we suggest that this approach should be avoided (see also Bapst 2013). Application of EH should be treated with caution, given the extremely old ages that it tends to produce, and the greater sensitivity of these ages to variation in phylogenetic topology. Nevertheless, EH is potentially a better choice when working with datasets with near-complete taxon sampling and high proportions of singleton taxa, in order to avoid the problems caused by the latter when applying cal3. Generally, the application of a single method to date just one tree cannot be recommended, because this would only sample a (not necessarily accurate) subset of possible ages. Future studies that require dated morphological trees need to do more to justify their choice of dating method(s), should build trees with the needs of dating in mind (e.g. sampling of the oldest putative clade members), and should tailor their choice of method to the properties of the available trees and the nature of the target group's fossil record.
Several potentially fruitful lines of enquiry into the accuracy of origination time estimation methods can be identified. First, our conclusions are based on the application of multiple methods to only two independent datasets (and should be regarded as a case study): this raises the question of whether the phenomena noted above are unique to neosuchians, or if they are repeated across many different clades? Second, we should investigate whether clade age estimates for given nodes are converging on a consensus through historical research time. New fossil discoveries are likely to produce a general trend towards finding older origination times, but there is likely to come a point where the 'correct' answer is approached. It would be useful to pinpoint where we are currently with respect to this scenario of initially wide uncertainty/high error followed by gradually increasing consensus, and how this trend varies between different clades. Finally, our recommendations for method selection are necessarily somewhat qualitative at this stage: it would be helpful to identify, for example, the proportion of singleton taxa and/or the extent of taxon coverage that determines whether Calc3 or EH yields the more accurate origination time estimates, through both simulations and more empirical case studies.