Biotic and abiotic factors and the phylogenetic structure of extinction in the evolution of Tethysuchia

Abstract

Crocodylomorpha is a large and diverse clade with a long evolutionary history now restricted to modern crocodilians. Tethysuchia is a less-inclusive clade of semi-amphibious taxa that crossed two biological crises: the second Oceanic Anoxic Event (OAE 2) and the Cretaceous/Paleogene (K/Pg) crisis. Numerous studies have sought to find the driving factors explaining crocodylomorph evolution, producing contradictory conclusions. Studies of included groups may be useful. Here, we study factors driving tethysuchian evolution using phylogenetically informed statistical analyses. First, we tested the phylogenetic structure of tethysuchian extinction at the OAE 2 and K/Pg crises. We then used phylogenetic comparative methods to test the influence of intrinsic (body size, snout proportion) and extrinsic (temperature, paleolatitude) factors on the evolution of tethysuchian diversity at the OAE 2 and the K/Pg crises. Finally, we tested whether temperature influenced the evolution of body size. We conclude that (1) extinction was not random in regard to phylogeny for Tethysuchia at the OAE 2 and K/Pg crises; (2) while an important tethysuchian turnover follows OAE 2, the K/Pg crisis was followed by an explosion in diversity of tethysuchians, probably linked to the colonization of emptied ecological niches; (3) tethysuchians lived in warmer environments after the OAE 2 crisis, possibly because of both global warming and latitudinal distribution shifts; (4) there is a significant change of snout proportion after the OAE 2 and the K/Pg crises, likely caused by niche partitioning; and (5) there is a positive correlation between body size and temperature, possibly because of a longer growth season.

Non-technical Summary

Crocodylomorpha is a large group of reptiles now restricted to modern crocodilians. Among them, Tethysuchia is a small group of semi-amphibious crocodiles that crossed two biological crises: the second Oceanic Anoxic Event (OAE 2) and the Cretaceous/Paleogene (K/Pg) crisis. Numerous studies have sought to find the driving factors explaining crocodylomorph evolution, producing contradictory conclusions. Studies of smaller groups may help find new exclusive patterns. Here, we studied factors driving tethysuchian evolution using phylogenetically informed statistical analyses. First, we tested whether or not tethysuchian extinction was random across the tips of phylogeny for both crises. Then, we tested the influence of biological (body size, snout proportion) and climatic (temperature, paleolatitude) factors on the evolution of tethysuchian diversity at the OAE 2 and K/Pg crises. Finally, we tested whether temperature influenced the evolution of body size. We conclude that (1) extinction was not random in regard to phylogeny for Tethysuchia at the OAE 2 and K/Pg crises; (2) while an important tethysuchian turnover follows OAE 2, the K/Pg crisis was followed by an explosion in diversity of tethysuchians, which may be explained by the disappearance of marine competitors such as mosasaurs; (3) tethysuchians lived in warmer environments after OAE 2, possibly because of both global warming and changes in latitudinal distribution; (4) there is an ecological diversification after both crises, observable by snout reduction, probably caused by niche partitioning; and (5) there is a positive correlation between body size and temperature, possibly because of a longer growth season.

Introduction

Crocodylomorpha is a diverse clade that emerged during the Late Triassic (Irmis et al. 2013) and occupied many ecological niches (Wilberg et al. 2019). It crossed major extinction events such as the Triassic/Jurassic (T/J) crisis, after which it radiated (Toljagić and Butler 2013; Bronzati et al. 2015), and the Cretaceous/Paleogene (K/Pg) crisis. Its diversity declined during the Cenozoic, probably due to climate cooling (Markwick 1998) or to competition with mammals in the case of terrestrial crocodylomorphs (Notosuchia) until modern days, when they are limited to 26 species sharing a similar semi-aquatic ecology (Grigg and Kirshner 2015).

Among crocodylomorphs, Tethysuchia Buffetaut, 1982 is a group of semi-aquatic freshwater and marine neosuchians (Andrade and Sayão 2014) that extended from the Kimmeridgian to the Bartonian (Jouve et al. 2021). They were probably ectothermic animals (Faure-Brac et al. 2021). While tethysuchians are ancestrally freshwater organisms (Martin et al. 2014b), independent events of colonization of the marine environment have been reported (Jouve et al. 2005a,b; Hua et al. 2007; Wilberg et al. 2019; Jouve 2021). This group is composed of two clades (Jouve 2021): Pholidosauridae Zittel and Eastman 1902, which extended from the Kimmeridgian (Mones 1980) to the Danian (Jouve and Jalil 2020) and Dyrosauroidea Jouve et al., 2021, which extended from the Barremian (Buffetaut and Hutt 1980) to the Bartonian (Buffetaut 1978). Tethysuchians faced two major extinction events. The second Oceanic Anoxic Event (OAE 2) occurred during the Cenomanian/Turonian transition, coinciding with intense volcanic activity, especially in the Caribbean, which produced metallic nutrients (Turgeon and Creaser 2008). These nutrients increased primary production, leading to greater oxidation of organic matter, stripping the ocean of oxygen, causing anoxia (Bralower 2008; Turgeon and Creaser 2008). Coupled with this, an important greenhouse effect may have been generated by the volcanic CO₂, leading to a stratified ocean that hampered oxygen delivery to deep waters (Bralower 2008; Turgeon and Creaser 2008). This event was linked to the extinction of ichthyosaurs (Fischer et al. 2016). The second event that tethysuchians crossed was the K/Pg crisis. That event coincides with a meteoric impact in the Gulf of Mexico (Hildebrand et al. 1991) and important volcanism in the Deccan Traps (Courtillot 1990). The timing and importance of each event remain heavily discussed (Schoene et al. 2019; Sprain et al. 2019). The K/Pg crisis was linked to the extinction of non-avian dinosaurs (Novacek 1999); pterosaurs (Barrett et al. 2008); and many marine reptiles, including mosasaurs and plesiosaurs (Bardet 1995).

Numerous studies have tried to identify factors driving crocodylomorph evolution. Martin et al. (2014a) suggested that Sea-Surface Temperature (SST) was positively correlated with crocodylomorph diversity, as well as with the marine colonization by tethysuchians, but they did not find a correlation between SST and tethysuchian diversity drops. Jouve et al. (2017) questioned the reliability of these results, stating that they were heavily affected by minor taxonomic updates. Mannion et al. (2015) found that diversification patterns for crocodylomorphs tracked environmental variations, but contrary to Martin et al. (2014a), no significant correlation between diversity and temperature was found for marine taxa. Jouve and Jalil (2020) found a significant positive correlation between paleotemperature and diversity during the Oxfordian–Cenomanian time interval followed by a significant negative correlation during the Turonian–Thanetian period. Bronzati et al. (2015) found that crocodylomorph diversification shifts were patchy and restricted to small intervals, whereas no such diversification shifts were found for tethysuchians. On the other hand, Jouve (2021) found an important diversification event for longirostrine (i.e., long-snouted) crocodylomorphs following the K/Pg crisis, especially regarding dyrosaurid tethysuchians. Godoy et al. (2019) did not find significant correlations between mean body size and temperature for crocodylomorphs, except for the period that extends from the Late Cretaceous to recent times. As for tethysuchians, the authors found different results depending on the body-size proxy and the paleotemperature data used. More recently, Stockdale and Benton (2021) found a significant correlation between mean body size and paleotemperature for crocodylomorphs. However, Benson et al. (2022) contested these results, pointing out the absence of log transformation before the statistical analyses. To sum up, no clear diversification driver has been found at the phylogenetic level of Crocodylomorpha. A wide ecological diversity, marked by many different lifestyles among crocodylomorphs (terrestrial, semi-aquatic, fully marine; see Wilberg et al. 2019) may explain these problems. Studies on less-inclusive groups, such as Tethysuchia, may help in finding new patterns and resolving this issue. Such studies, however, remain scarce. A new approach coding extinction/survival as a binary variable was applied recently to Notosuchia, a group of largely terrestrial crocodylomorphs (Aubier et al. 2023). These authors tested the phylogenetic structure of extinction during the K/Pg crisis and used Phylogenetic Logistic Regression (PLR) to test the factors influencing survival during the K/Pg crisis. These analyses revealed a phylogenetic structure in notosuchian extinction at the K/Pg crisis and an evolutionary trend toward larger body sizes after this crisis. This last trend was tentatively explained as being the outcome of a dietary shift (Aubier et al. 2023).

The present paper is aimed at elucidating the phylogenetic structure of extinction and identifying the biotic and abiotic factors driving the evolution of tethysuchian paleobiodiversity. More precisely, we tested the phylogenetic structure of tethysuchian extinction at the OAE 2 and K/Pg crises. Then, we tested the effect of intrinsic (body size, snout proportion) and extrinsic (paleolatitude, paleotemperature) factors on the evolution of tethysuchians at both crises. As paleotemperature seems to play a varying role in tethysuchian diversity depending on the time period considered (Jouve and Jalil 2020; Jouve 2021), we expect temperature to be significantly associated with the probability of belonging to the post-OAE 2 fauna. As there seems to be an overall increase in mean body size in crocodylomorphs through time (Godoy et al. 2019), we expect this overall trend to remain unaffected by the crises and body size to be correlated with the probability of belonging to the fauna that existed after the OAE 2 and the K/Pg crises. Finally, we tested whether paleotemperature is linked to body-size evolution. Previous studies did not find significant correlations between these variables in crocodylomorphs (Godoy et al. 2019; Benson et al. 2022). However, mixed results were obtained when focusing on Tethysuchia (Godoy et al. 2019). As most of these results were not significant, we do not expect to find a correlation between size and temperature in tethysuchians.

Materials and Methods

Data Acquisition

A primary set of taxa was gathered using the Paleobiology Database (PBDB; https://paleobiodb.org). To account for potential errors, we consulted the primary literature to ensure the reliability of the data on various aspects (location, age, taxonomy, etc.). As most of the fossil record consists of skulls (Buffetaut and Hutt 1980; Hastings et al. 2011; Jouve et al. 2021), we initially chose three cranial variables: skull length (SL; from the anterior tip of the premaxilla to the posterior end of skull table), skull width (SW) at mid-orbital length, and snout proportion (SP; from the tip of the premaxilla to the anterior margin of the orbits, relative to SL). If measurements were not available from the literature or not explicitly stated to be the same as defined, we measured them, using Photofiltre software (see Supplementary File 1 for details) on published figures. As complete tethysuchian remains are scarce (Sereno et al. 2001; Jouve et al. 2006), we could not directly compare body sizes. Therefore, one of our cranial measurements had to be selected as a proxy for body size. O'Brien et al. (2019) mentionned that SW at the quadrates is a good proxy for body size for extant crocodilians. In their study, this proxy seemed to provide accurate results for Sarcosuchus imperator De Lapparent De Broin and Taquet, 1966. However, lateral compression and poor preservation of the specimens only allowed measurements at mid-orbital length. On the other hand, SW at mid-orbital length remains a missing variable in most of our sample (see Supplementary File 1 for more information). As a result, SW was excluded from further analyses. SL is the most available skull metric and has previously been used as a proxy for body size (e.g., Godoy et al. 2019; Aubier et al. 2023). However, studies have criticized this variable as subject to biases linked to group differences and have tried to address this problem using various methods (Young et al. 2011; Stockdale and Benton 2021). Most recently, Stockdale and Benton (2021) have applied principal component analysis (PCA) using various body size indicators to distinguish independent components linked to body size. However, Benson et al. (2022) pointed out that the first principal component was still heavily linked to SL. Therefore, their analyses were still heavily biased by this metric. Furthermore, specimens included in this study are not sufficiently preserved to provide satisfying measurements with respect to the PCA analysis. Another approach is to use equations to estimate the total body size using long bones (Farlow et al. 2005; Vandermark et al. 2007). However, most of these equations are based on extant crocodilians, particularly Alligator mississipiensis (Daudin, 1802), and using them for fossil species would rely on the assumption that there is not much difference in allometry between extant and extinct crocodylomorphs. However, Young et al. (2011) considered this assumption unlikely and devised an entirely new equation for metriorhynchids to counter the problem. As we lack complete tethysuchian remains, we cannot test this assumption for Tethysuchia. Therefore, we chose to keep SL as a body-size proxy. In the case of Meridiosaurus vallisparadisi Mones, 1980 and Sabinosuchus coahuilensis Shiller et al., 2016, SL measurements were not available, but rather estimations based on the length from the tip of the premaxilla to the last maxillary tooth (Fortier et al. 2011) and on the mandible length (Shiller et al. 2016), respectively. We coded their SLs accordingly and then conducted another set of analyses that excluded these estimations (see details in Supplementary File 1). Before any analysis, we log-transformed the measurements, as advised by Benson et al. (2022).

We gathered the paleoenvironments of analyzed taxa using Jouve (2021). It can be hypothesized that some Tethysuchia could possibly move between fresh and salt water, like some modern crocodilians (Grigg and Kirshner 2015). However, modern crocodilians capable of this behavior can only stay in salt water for a limited period and need at least partly functional salt glands to deal with various osmolarity problems (Grigg and Kirshner 2015). Although some tethysuchian species have been described as living in a “marine-influenced” environment that has intermediate levels of salinity, the capacity to move “at will” between marine and freshwater environments seems unlikely. Indeed, most of the species included in this study are buried either in freshwater-only or marine-only localities (Jouve 2021). Therefore, we considered that the various specimens were buried in their preferred environments and were coded accordingly. Dakotasuchus kingi Mehl, 1941 is the only exception, as it was found in one marine and two freshwater localities (Jouve 2021). It seems more likely that it had been transported from freshwater to a marine environment than the opposite, so we considered D. kingi a freshwater species. On another note, MHNM-kh01 is a heavily damaged specimen in the abundant and well-preserved Ouled Abdoun Basin which is marine (Jouve and Jalil 2020). Its state of preservation is striking compared with the other Tethysuchia from this formation (Jouve et al. 2005a,b, 2006, 2008b). Such a difference in preservation suggests transport from a freshwater to a marine locality (Jouve 2021). Therefore, we considered MHNM-kh01 to be a freshwater organism.

Moreover, the local maximum absolute paleolatitude recorded for each species was gathered using the PBDB, and local paleotemperatures were inferred using latitudinal temperature gradients from the literature considering the aforementioned paleolatitudes. However, we lack latitudinal temperature gradients for freshwater environments. Freshwater temperatures have been proposed to be close to the terrestrial ones (Newton and Mudge 2003; Pouech et al. 2014). Furthermore, tethysuchians spent time out of the water, as they had a semi-aquatic lifestyle (Andrade and Sayão 2014). Therefore, we used terrestrial temperature gradients for the species inferred as coming from freshwater environments. For marine species, we used SST gradients gathered from the literature (Frakes et al. 1994; Amiot et al. 2004; Pouech et al. 2014; Alberti et al. 2017; Zhang et al. 2019; Laugié et al. 2020; see details in Supplementary File 1). Note that no extensive latitudinal temperature gradient study has been made for the Danian. As δ¹⁸O levels between the Maastrichtian and Selandian are rather similar (Prokoph et al. 2008), we considered the mean value between Campanian–Maastrichtian and Selandian–Thanetian to be a proxy for the value of the Danian.

Supertree

Because phylogenetic comparative methods (PCMs) require a phylogeny, we decided to use the topology from Jouve (2021) as a reference. It includes the largest tethysuchian sample (n = 35) and provides an extensive review of phylogenetic relationships among both pholidosaurids and dyrosaurids. We added Brachiosuchus kababishensis Salih et al., 2022, which has been recovered as the second-earliest diverging dyrosaurid (Salih et al. 2022). Dakotasuchus kingi has a debated phylogenetic position (Jouve and Jalil 2020). It is considered to be part of the clade including Terminonaris robusta (Mook, 1934), Terminonaris browni (Osborn, 1904), and MHNM-kh01, a poorly preserved Danian specimen (Jouve and Jalil 2020) or the sister species of Pholidosaurus Meyer, 1841. We constructed a supertree for each of these two hypotheses, subsequently named Jouve 1 and Jouve 2, respectively. We also tested the topologies obtained by Sachs et al. (2021), the only ones with a satisfying Tethysuchia sample (i.e., more than 20 species), although this phylogeny was initially designed for testing phylogenetic relationships among crocodyliforms and not specifically Tethysuchia. Its most striking difference with Jouve 1 and 2 was that Vectisuchus leptognathus Buffetaut and Hutt, 1980 and Elosuchus De Lapparent De Broin, 2002 are considered pholidosaurids. Sachs et al. (2021) yielded two topologies: the first one retrieves Pholidosaurus schaumburgensis Meyer, 1841 in a clade with Oceanosuchus boecensis Hua et al., 2007 and T. robusta. The second, on the other hand, retrieves P. schaumburgensis as a sister species of the clade including O. boecensis, T. robusta, Chalawan thailandicus (Buffetaut and Ingavat, 1980) and Sarcosuchus De Lapparent De Broin and Taquet, 1966. These topologies are subsequently named Sachs 1 and Sachs 2, respectively. Other species listed in the PBDB could have been added but were excluded for various reasons. Anglosuchus geoffroyi (Owen, 1884) and Anglosuchus laticeps (Owen, 1884) are considered Bathonian. However, their ages remain doubtful, and they closely resemble Pholidosaurus purbeckensis (Mansel-Pleydell, 1888), so they may be synonyms of the latter (Jouve and Jalil 2020). The pholidosaurids Pholidosaurus milwardi Roxo, 1929 and Pholidosaurus meyeri Dunker, 1843 and the dyrosaurids Tilemsisuchus lavocati Buffetaut, 1980, Congosaurus compressus (Buffetaut, 1980), and Rhabdognathus acutirostris Bergounioux, 1955 combine poor information on their anatomy, locality, age, and/or phylogenetic position.

As mentioned earlier (see previous section), stratigraphic data were gathered using both the PBDB and primary literature. For taxa restricted to a single formation, we considered their FAD (first appearance datum) and LAD (last appearance datum) to match the stratigraphic extent of the formation. For species having occurrences in multiple formations, we considered their FADs and LADs to be as restrictive as possible: we selected the shortest time interval in which the species could be present in all of its recorded localities. However, some adjustments had to be made. Phu Kradung Formation (Thailand), where C. thailandicus was recovered, has been traditionally considered as Kimmeridgian–Tithonian according to vertebrate data (Buffetaut and Suteethorn 2007). However, recent palynology studies suggested a Berriasian age (Racey and Goodall 2009). Therefore, we considered C. thailandicus to be of Berriasian age. Hyposaurus natator Troxell, 1925 is noted as being Maastrichtian in the PBDB (Cope, 1866; Marsh, 1870). However, reviews argued that there was probably a reworking caused by bioturbation that caused Danian fossils to be trapped in an apparent Maastrichtian site (Landman et al. 2007; Wiest et al. 2016). Therefore, in our analyses, we considered it to be Danian.

Topologies were dated using the timePaleoPhy function on the paleotree package (Bapst 2012) in R v. 4.2.2 (R Core Team 2013). We used the firstLast dating method, which considers the FAD–LAD interval as a positive presence of the taxa. The nodes were dated using the mbl (minimum branch length) method, which considers the age of a node to be the same age as the FAD of the oldest fossil of the node. Therefore, FADs and LADs remain the only range data used. We must consider that this method may generate zero-length branches (ZLBs), which are intractable for many PCMs (Soul and Wright 2021). A minimal branch length can be selected to prevent ZLBs (Laurin 2004; Wang and Lloyd 2016). Here, we set it to 1 Myr using the “vartime” argument. The complete dated supertrees include 36 Tethysuchia for the phylogenies adapted from Jouve 1 and 2 and 25 Tethysuchia for Sachs 1 and 2 (see Fig. 1 and Supplementary File 2). The complete dataset, R script, and generated nexus trees are in Supplementary Files 3–5.

Figure 1. Supertree of Tethysuchia, the topology shown here is Jouve 1. The green spot indicates the Pholidosauridae; the red spot, Dyrosauroidea; the orange spot, Dyrosauridae; the yellow spot, Phosphatosaurinae; and the black spot, Hyposaurinae. The alternative topologies can be observed in Supplementary File 2.

Faunal Attribution

Each species was assigned to a fauna depending on whether its stratigraphic interval extended before or after the OAE 2 and the K/Pg crises. For the OAE 2 crisis, 15 taxa from Jouve 1 and 2 in the Kimmeridgian–Turonian time bin are referred to as “pre-OAE 2 fauna” (12 taxa for Sachs 1 and 2). The other 21 taxa (13 in Sachs 1 and 2) extend from the Campanian to the Ypresian and are referred to as “post-OAE 2 fauna.” Regarding the K/Pg crisis, 18 taxa extend from the Kimmeridgian to the Maastrichtian and are defined as “pre-K/Pg fauna” (13 in Sachs 1 and 2). The other 18 taxa (12 in Sachs 1 and 2) extend from the Danian to the Ypresian and are defined as “post-K/Pg fauna.” Thus, each crisis separates two large time bins. These time bins will be used to test differences between pre- and postcrisis faunae (see following sections) rather than to analyze the evolution of a trait through time as previous studies have done (this last procedure requires a larger sample to infer evolutionary rates; see Stockdale and Benton 2021). Here, these faunae are assumed to be homogeneous, a strong assumption considering the long time bins involved.

D-statistic

To check whether the extinction across the OAE 2 and K/Pg has a phylogenetic structure or not, we used the D-statistic (Fritz and Purvis 2010). This method measures the randomness of the extinction distribution across the tips of a given tree. More precisely, it compares the observed distribution of a binary variable (in this case, extinction vs. survival, coding each species in the “precrisis” fauna as 0 and each species in the “postcrisis” fauna as 1) with two other distributions: one that simulates the evolution of the binary trait under a Brownian model of evolution and one that simulates the evolution of the same trait under a random model of evolution. The analysis generates a D-value. If this value is equal to 1, extinction is not considered to be phylogenetically structured (i.e., the observed distribution is the same as the one produced under the simulated random evolutionary model). If extinction is clustered in the phylogeny as if it followed a Brownian evolutionary model, the D-value would equal zero. D-values can fall outside this range. This method has been used before to check extinction risk for extant organisms (Fritz and Purvis 2010; Yessoufou et al. 2012; Fontana et al. 2021) or extinction selectivity in the fossil record (Allen et al. 2019; Wilke et al. 2020; Aubier et al. 2023). We used the phylo.d function of the caper package (Orme et al. 2013) in R v. 4.2.2 (R Core Team 2013), selecting 1000 permutations (i.e., repetition of the simulations to scale D-values), as suggested by Fritz and Purvis (2010). This function provides the D-value, as well as the probability of obtaining this D-value if extinction was phylogenetically random and if it was phylogenetically structured. We performed four analyses depending on the phylogenies (Jouve 1 and 2, Sachs 1 and 2; see previous section). We excluded V. leptognathus, as it caused heteroscedasticity in the phylogenetic generalized least squares (PGLS) analysis (see “PGLS and Variation Partitioning”). The complete dataset and script can be found in the Supplementary Files 3 and 6.

PLR

We used PLR (see Ives and Garland 2010) to test whether body size, SP, paleolatitude, and paleotemperature affected the probability of belonging to the post-OAE 2 or post-K/Pg faunae. We used the phyloglm function from the phylolm R package (Tung Ho and Ané 2014) in R v. 4.2.2 (R Core Team 2013). The PLR allows the production of predictive models for a binary dependent variable using a set of explanatory variables and the phylogeny. As observations between organisms are not independent (Felsenstein 1985), we included the dated trees (see “Supertree”). PLR has been used to infer the probability of endothermy in tetrapods (Cubo et al. 2023; Faure-Brac et al. 2024) and the probability of survival after the K/Pg crisis in Notosuchia (Aubier et al. 2023), similar to our study. We used the same coding as in the D-statistic (see previous section). We performed four sets of analyses depending on the phylogenies (Jouve 1 and 2, Sachs 1 and 2; see “Supertree”). In each set, we tested four models: log-transformed SL (model A), SP (model B), paleotemperature (model C), and paleolatitude (model D). For the latter two, we considered that closely related species have a tendency to live in proximity and/or share similar habitats, following Phylogenetic Niche Conservatism (PNC; Ackerly 2003; Cooper et al. 2010). However, the alternative may be possible. Therefore, we also tested the influence of paleolatitude and paleotemperature using non-phylogenetic logistic regressions with Generalized Linear Models (GLMs). Furthermore, for each set, we considered an alternative hypothesis that excluded SL estimations for M. vallisparadisi and S. coahuilensis (see “Data Acquisition”). In total, each set comprised 24 analyses with 12 per biological crisis (8 comprising PNC, and 4 discarding it). Early analyses suggested that V. leptognathus caused heteroscedasticity in the PGLS analysis (see next section). Therefore, it was subsequently removed from both PLR and PGLS analyses. The complete dataset and script can be found in Supplementary Files 3 and 7.

PGLS and Variation Partitioning

Many studies have previously tried to find a correlation between body size and paleotemperature, producing mixed results (see “Introduction”). Here, we used the PGLS method (see Grafen and Hamilton 1989) to test whether temperature affected log-transformed SL. We used the pgls function from the caper R package (Orme 2013) in R v. 4.2.2 (R Core Team 2013). We performed four sets of analyses depending on the phylogenies (Jouve 1 and 2, Sachs 1 and 2; see “Supertree”). In each set, we tested the relationship for tethysuchians as a whole, as well as for pholidosaurids and dyrosauroids separately. As mentioned earlier, temperature may be independent from phylogeny (see previous section); therefore, we also used Generalized Least Squares (GLS) and classic Linear Models (LMs) that do not take phylogenies into account. Furthermore, for each set, we considered an alternative hypothesis that excluded SL estimation for M. vallisparadisi and S. coahuilensis (see “Data Acquisition”). Each set had a total of 18 analyses (6 per group). Shapiro-Wilk tests (Shapiro and Wilk 1965) were used to test the normality of the residual distribution. To check for homoscedasticity, we used the Breusch-Pagan test, which measures the regression error variance (Breusch and Pagan 1979). Homoscedasticity was not respected if V. leptognathus, which had a very short skull (Salisbury and Naish 2011) and lived in very cold temperatures (Frakes et al. 1994), was included. Thus, it was removed from the sample. Then, to test the quality of the paleotemperature-influenced model, we calculated the corrected Akaike Information Criterion (AICc) using the AICc function from the AICcmodavg package (Mazerolle 2013) in R v. 4.2.2 (R Core Team 2013) and compared it with a null model (i.e., no influence). We used the same sets of analyses as in the test of correlation between body size and temperature (see above). Finally, to estimate the impact of PNC on log-transformed SL variation, we used the variation partitioning method, which allows quantification of the relative impact of various components on an explanatory variable (Borcard et al. 1992). This method has been further developed to consider phylogeny as a component using a matrix of principal coordinates representing phylogeny (Desdevises et al. 2003; Peres-Neto et al. 2006; Montes et al. 2007; Piras et al. 2009; Sakamoto et al. 2010). We used the varpart function from the vegan R package (Dixon 2003) in R v. 4.2.2 (R Core Team 2013). We analyzed the variation of log-transformed SL using two components: ecology, which can be represented either by paleotemperature or paleolatitude; and phylogeny. For the latter, we retained a set of axes that contributed for more than 80% of the total variation of the phylogenetic distance matrix. As a result, we obtained four different partitions (Fig. 2): a fraction corresponding to a strictly ecological impact on log-transformed SL variation (partition A), a fraction corresponding to a strictly phylogenetic impact on log-transformed SL variation (partition B), a fraction corresponding to a combined effect of ecology phylogeny on log-transformed SL variation (partition C), and a partition corresponding to the unexplained variation (partition D). We can test the significance of partitions using redundancy analysis, except for partitions C and D. The complete dataset and script can be found in Supplementary Files 3, 8, and 9.

Figure 2. Representation of variation partitioning for a dependent variable, the gray rectangle represents all of the variation of the dependent variable. Four different partitions are proposed: partition A corresponds to the strictly ecological impact on variation, partition B corresponds to the strictly phylogenetic impact on variation, partition C corresponds to the common impact of phylogeny and ecology (Phylogenetic Niche Conservatism), and partition D corresponds to the unexplained part of variation.

Results

Testing the Phylogenetic Structure of Extinction at the OAE 2 and the K/Pg Crises

Similar results were provided by all four analyses. Indeed, in the topology Jouve 1, which considers that Dakotasuchus kingi belongs to the clade including Terminonaris and MHNM-kh01, we observe D-values of −1.004 for the OAE 2 crisis and −0.751 for the K/Pg one (Table 1, Jouve 1). These values mean that the distribution of the extinction is more phylogenetically structured than that obtained from the simulations performed under a Brownian evolutionary model. Likewise, negative D-values were yielded for the topology Jouve 2, which considers that D. kingi belongs to the clade including Pholidosaurus (D _OAE2 = −1.037 and D _K/Pg = −0.723; Table 1, Jouve 2). This was also the case for Sachs 1 and 2, which respectively consider that Pholidosaurus schaumburgensis is in a clade comprising Oceanosuchus boecensis and Terminonaris robusta (D _OAE2 = −1.388 and D _K/Pg = −0.704, Table 1, Sachs 1), and that P. schaumburgensis is the sister species of a clade including O. boecensis, T. robusta, Chalawan thailandicus, and both Sarcosuchus species (D _OAE2 = −1.331 and D _K/Pg = −0.74; Table 1, Sachs 2). These results show that extinction was not phylogenetically random at both of the studied crises. Rather, they show a phylogenetic structure of extinction (i.e., closely related species went extinct during both of the studied crises) that is robust enough to be independent from the phylogenetic placement of one or two species depending on topologies.

Table 1. Results from the D-statistic analysis for second Oceanic Anoxic Event (OAE 2) and Cretaceous/Paleogene (K/Pg) crisis. The first topology is the same as in Fig. 1. The second topology shows Dakotasuchus kingi in a clade including Pholidosaurus cherves, Pholidosaurus purbeckensis, and Pholidosaurus schaumburgensis. The third topology shows P. schaumburgensis in a clade with Oceanosuchus boecensis. The fourth topology retrieves P. schaumburgensis as a sister clade of the clade including O. boecensis, Terminonaris robusta, Sarcosuchus, and Chalawan thailandicus. These alternative topologies are provided in Supplementary File 2.

Testing the Effect of Biotic and Abiotic Factors on the Evolution of Tethysuchia after the OAE 2 and K/Pg Crises

For the first Jouve topology (Table 2, Jouve 1), the probability of belonging to the post-OAE 2 fauna is significantly explained by SP (model B) and paleotemperature (model C) but not by the log-transformed SL (model A) or paleolatitude (model D). The probability of belonging to the post-K/Pg crisis fauna is significantly explained by SP (model B) but not by paleotemperature (model C), log-transformed SL (model A), and paleolatitude (model D). The second topology produced similar results (Table 2, Jouve 2). The alternative analysis excluding Meridiosaurus vallisparadisi and Sabinosuchus coahuilensis (see “Materials and Methods”) produced similar results for both topologies (Supplementary File 10, PLR, Jouve 1 and 2). The paleotemperature estimate is positive and significant for each analysis testing its effect on the probability of belonging to the post-OAE 2 fauna, which means the variable is positively correlated with the probability of belonging to the post-OAE 2 fauna. Tethysuchians in the post-OAE 2 fauna are more likely to live in warmer climates. According to analyses using Jouve's (2021) topologies. The coefficient for SP is significant and negative for each analysis testing its effect on the probability of belonging to the post-OAE 2 and post K/Pg faunae, which means the variable is negatively correlated with the probability of belonging to the post-OAE 2 and the post-K/Pg faunae. Tethysuchians belonging to the postcrisis faunae are more likely to be short-snouted according to analyses using Jouve's (2021) topologies. However, analyses using Sachs 1 and 2 topologies yielded different results: the probability of belonging to the post-OAE 2 or to the post-K/Pg fauna is not affected by any of our models regardless of topology (Table 2, Sachs 1 and 2). The alternative hypothesis that excludes M. vallisparadisi and S. coahuilensis yields similar results (Supplementary File 10, PLR, Sachs 1 and 2). Finally, GLMs that discard PNC for paleotemperature and paleolatitude (see “Materials and Methods”) retrieved a positive effect of paleotemperature on the probability of belonging to the post-OAE 2 and K/Pg faunae regardless of topology (Table 2, Supplementary File 10, PLR). Paleolatitude also has a negative effect on the probability of belonging to the post-OAE 2 fauna. Discarding PNC shows that Tethysuchia are more likely to live in lower latitudes and warmer environments after OAE 2 and in warmer environments after the K/Pg crisis (Table 2). The alternative hypothesis that excludes M. vallisparadisi and S. coahuilensis yielded similar results (Supplementary File 10, PLR). To sum up, analyses using Jouve's (2021) hypotheses indicate a trend to warmer climates after OAE 2 and shorter snouts after K/Pg, whereas analyses using Sachs et al.'s (2021) hypotheses indicates no trend, if PNC is taken into consideration.

Table 2. Results from the phylogenetic logistic regression (PLR) and generalized linear model (GLM) analyses; significant p-values are lower than 0.05. The first topology is the same as in Fig. 1. The remaining topologies are in the same order as in Table 1. *p < 0.05; **p < 0.01; ***p < 0.001.

Testing the Correlation between Body Size and Temperature

Both Jouve topologies yielded similar results (Table 3, Jouve 1 and 2, and Fig. 3). We found a significant positive correlation between paleotemperature and log-transformed SL for tethysuchians and pholidosaurids. On the other hand, we found no significant correlation for dyrosauroids. In both Sachs topologies, we found a significant positive correlation between paleotemperature and log-transformed SL for tethysuchians (Table 3, Sachs 1 and 2). However, it should be noted that in the latter two, residuals did not follow a normal distribution. Therefore, these results are not statistically definitive. The alternative analysis considering M. vallisparadisi and S. coahuilensis SL as missing provides different results: we find no significant correlation between paleotemperature and log-transformed SL for any groups and topologies (see Supplementary File 10, PGLS). The paleotemperature model has a lower AICc than the null model for Tethysuchia in both Jouve topologies (Table 4, Jouve 1 and 2). However, the null model has a lower AICc than the paleotemperature model for dyrosauroids and pholidosaurids. For both Sachs topologies, we see close AICc values between the null model and the paleotemperature model for Tethysuchia (<0.5; see Table 4, Supplementary File 10, AICc), which indicates that the models are not different. However, for dyrosaurids and pholidosaurids, AICc is generally lower in the null model. The alternative hypothesis excluding M. vallisparadisi and S. coahuilensis yields similar results for the Jouve 1 and 2 topologies (Supplementary File 10, AICc). However, for the Sachs 1 and 2 topologies, AICc is lower in the null model for Tethysuchia.

Table 3. Results from the phylogenetic generalized least squares (PGLS), generalized least squares (GLS), and linear models (LM) analyses, significant p-values are lower than 0.05. The first topology is the same as in Fig. 1. The remaining topologies are in the same order as in Table 1. *p < 0.05; **p < 0.01; ***p < 0.001; ¹Nonnormal, p = 0.038; ²Nonhomogenous, p = 0.032.

Figure 3. Phylogenetic generalized least squares (PGLS) curve for tethysuchians (blue), pholidosaurids (green), and dyrosauroids (red). The circles correspond to Pholidosauridae species, and the triangles correspond to Dyrosauroidea species.

Table 4. Comparison of corrected Akaike information criterion (AICc) between a paleotemperature-influenced model and a null model for the phylogenetic generalized least squares (PGLS), generalized least squares (GLS), and linear models (LM) analyses. The topologies are in the same order as in Table 1.

If we discard PNC (i.e., if we rely on GLS and LM), no correlation is found between log-transformed SL and paleotemperature regardless of topologies and coding for M. vallisparadisi and S. coahuilensis (Table 3, Supplementary File 10, PGLS). Using this assumption, AICc is always lower in the null model (Table 4, Supplementary File 10, AICc). To sum up, paleotemperature has a positive correlation with log-transformed SL in Tethysuchia only if PNC is considered and if M. vallisparadisi and S. coahuilensis are not excluded. The paleotemperature model generally has lower AICc values than the null model (and therefore is the better model) in Jouve's (2021) topologies for Tethysuchia if PNC is considered and if M. vallisparadisi and S. coahuilensis are not excluded.

Finally, when testing for variation partitioning, both Jouve topologies yield similar results (Table 5, Jouve 1 and 2). Regardless of the explanatory variable composing the ecological component (i.e., paleolatitude or paleotemperature), we observe that partition B accounts for around 5% of the variation. However, we note that most of the variation remains unexplained. Neither partition A nor partition B is significant when tested with redundancy analysis. The alternative hypothesis that excludes M. vallisparadisi and S. coahuilensis provides slightly different results. Partitions A and B remain nonsignificant in redundancy analyses, but partition B provides a negative R ², while partition D accounts for around 100% of the variation for each of the analyses (Supplementary File 10, Variation Partitioning, Jouve 1 and 2). These results suggest an important effect of the two removed species on the results. In both Sachs topologies, partition D accounts for around 100% of the variation and the R ² values for the other partitions are either negative or up to 2% of the variation (Table 5, Sachs 1 and 2). Similar results can be observed when M. vallisparadisi and S. coahuilensis are excluded from the analysis: partition C contributes to around 2% of log-transformed SL when paleotemperature is the ecological component and less than 1% if paleolatitude is the ecological component. The rest of the variation is unexplained (Supplementary File 10, Variation Partitioning, Sachs 1 and 2).

Table 5. Results from the variation partitioning analyses, adjusted R ² is noted along with p-values, if possible, within parentheses.

Discussion

A Differential and Phylogenetically Structured Response to Biotic Crises

The first major peak of tethysuchian diversity occurs during Cenomanian (Jouve and Jalil 2020; Jouve 2021; Fig. 1, Supplementary File 2). This period corresponds to the highest temperature and sea level of the Mesozoic (Vérard et al. 2015; Scotese et al. 2021), which may explain the important tethysuchian fossil record, because high sea level has long been considered a factor of enhanced diversity (Martin et al. 2014a; Mannion et al. 2015; Tennant et al. 2016). During the Cenomanian/Turonian transition, Tethysuchia experienced a major diversity drop corresponding to OAE 2 (Jouve and Jalil 2020; Fig. 1, Supplementary File 2). Because half of the tethysuchians at the time were marine (Jouve 2021), they were probably heavily affected by this event, which was also linked to the extinction of ichthyosaurs (Fischer et al. 2016) and the diversification of mosasaurs (Bardet 1995). These patterns suggest an important marine faunal turnover previously mentioned in the literature (Kauffman 1995; Wan et al. 2003; Caron et al. 2006; Monnet 2009). This turnover is supported by the D-statistic analysis, which shows a phylogenetic structure of extinction at OAE 2 (Table 1). Indeed, most pholidosaurids do not survive the crisis (Fig. 1, Supplementary File 2) and all known Dyrosauridae De Stefano, 1903 appear after the crisis. OAE 2 marks a transition from pholidosaur- to dyrosaurid-dominated faunae.

Following the OAE 2, a gap in the tethysuchian fossil record occurs from the Coniacian to the Santonian (Jouve and Jalil 2020; Fig. 1, Supplementary File 2). The only known putative tethysuchian remains during this period are a partial maxilla fragment from the In Beceten Formation of Niger that is described as being similar to Tethysuchia, although no phylogenetic analysis is possible because of its fragmentary nature (Buffetaut 1974; Meunier and Larsson 2018). The Coniacian–Santonian interval coincides with a marine regression (Jouve and Jalil 2020), which can explain this drop in diversity. However, if we look at crocodylomorphs as a whole, most of the fossil record during the Coniacian–Santonian consists of fragmentary remains (Puértolas-Pascual et al. 2016; Meunier and Larsson 2018). Therefore, some of these crocodylomorph elements may have belonged to tethysuchians but have not been identified as such because the material is too fragmentary to provide a more precise taxonomic attribution. Tethysuchian biodiversity may also have been further underestimated due to sampling biases: Coniacian–Santonian formations may suffer from a lack of interest compared with other Late Cretaceous periods that are closer to major events such as OAE 2 and the K/Pg crisis. The next tethysuchian occurrences are recorded during the Campanian and the Maastrichtian (Halstead 1975; Shiller et al. 2016; Jouve and Jalil 2020; Salih et al. 2022). Most Late Cretaceous tethysuchians lived in freshwater environments (Jouve 2021). These environments were relatively spared during the K/Pg crisis, as increased potential for dormancy (i.e., a metabolically slowed or inactive state in response to harsh conditions that limits starvation), faster production recovery, more abundant detrital food sources, and the presence of eventual thermal refuges in those environments may have helped stabilize the trophic networks (Robertson et al. 2013).

Following the K/Pg crisis, an explosion in diversity occurs (Jouve 2021; Fig. 1, Supplementary File 2). Most Cenozoic tethysuchians lived in marine environments (Jouve 2021). This colonization from freshwater to marine environments may have been made possible because tethysuchians took over the niches vacated by mosasaurs and plesiosaurs that became extinct during K/Pg (Barbosa et al. 2008; Jouve et al. 2008a,b; Bardet et al. 2017; Jouve 2021). As shown by the D-statistic analysis, this diversification was phylogenetically structured (Table 1), because most of early-diverging Dyrosauridae do not cross the K/Pg boundary. Because extinction is phylogenetically structured, so is the subsequent diversification. Indeed, the postcrisis diversification affects mostly dyrosaurids, especially Hyposaurinae Nopcsa, 1928 (Fig. 1, Supplementary File 2) that heavily colonized the marine realm during the Paleogene (Jouve 2021).

Both crises had an impact on tethysuchian biodiversity: OAE 2 caused a turnover in tethysuchian diversity, likely by destabilizing the marine food chain, causing top predators such as marine pholidosaurids to become extinct (Jouve and Jalil 2020), whereas the K/Pg crisis made tethysuchian diversity explode, likely as a result of the colonization of niches vacated by mosasaurs and plesiosaurs. After a thriving period during Paleocene, Tethysuchia's evolutionary history ends with their extinction during the Bartonian (Jouve 2021). The factors explaining their extinction are uncertain (Amoudji et al. 2021). The Bartonian coincides with the beginning of the late Eocene–Oligocene cooling (Scotese et al. 2021), which may have impacted tethysuchians (Jouve 2021). Another hypothesis involving competition with new predators, including cetaceans, was mentioned by Hastings (2012). It was considered unlikely by Martin et al. (2014a), but was still mentioned by Stubbs et al. (2021). These hypotheses remain to be tested. Most recently, Scott and Anderson (2023) have tested, under the postulate that competition increases as morphological similarity increases, the competitive interactions between gnathostomes and agnathans during the late Silurian–Devonian using distance-based morphometrics. However, we lack fossil sites bearing both Tethysuchia and cetaceans to support the competition. Therefore, such an assumption cannot be tested.

An Adaptation to Warmer Temperatures and Morphological Changes after the Biotic Crises

PLR analyses showed that post-OAE 2 tethysuchians lived under higher temperatures than those pre-OAE 2 if we follow Jouve's (2021) hypotheses. We observe similar results if we discard PNC for Sachs et al.'s (2021) hypotheses (Table 2). Climatic data suggest that, except for the Cenomanian, post-OAE 2 mean temperature values were generally higher than those of the pre-OAE 2 periods (Scotese et al. 2021). Therefore, these results could be explained by an overall global warming after OAE 2. A change in latitudinal distribution between the two faunae could also explain these results. However, we find no significant difference between them if we consider PNC (Table 2). On the contrary, if we discard PNC, we observe that Tethysuchia were more likely to live at lower latitudes after the OAE 2. Both pre- and post-OAE 2 faunae have a wide latitudinal range (11.6°–53.9° and 3.5°–40.6°, respectively; see Fig. 4). However, two post-OAE 2 tethysuchians (Sabinosuchus coahuilensis and Hyposaurus natator, located in Mexico and the eastern United States, respectively) have positions relatively isolated from the others. These are not clear outliers; however, if they are excluded, the latitudinal range of the post-OAE 2 fauna is highly reduced (3.5°–27°; see Fig. 4B). Thus, the width of the latitudinal range of this fauna is largely due to only two species. Therefore, temperature differences between both faunae may be caused by the combined effect of an overall temperature increase and a generally more restrictive latitudinal distribution (although not necessarily statistically different for the latter). GLM analyses that discard PNC show that post K/Pg Tethysuchia also lived in warmer environments. Although paleolatitude is similar between both faunae, the literature shows an overall warming after K/Pg, most notably during the end of the Danian and the Paleocene–Eocene thermal maximum (see Scotese et al. 2021). We do not find any correlation between paleotemperature and the probability of belonging to the post K/Pg fauna if we consider PNC. However, the recorded fauna during the Campanian and Maastrichtian, which represent the period right before the K/Pg crisis, is still quite limited (n = 3). Therefore, a larger Campanian–Maastrichtian sample could heavily impact the statistical analyses. Further analyses may be needed to test whether the formation of paleocurrents may influence local temperature; especially for Paleogene, where marine forms are numerous (Jouve 2021). Indeed, a proto–Gulf Stream has been suggested in literature (Watkins and Self-Trail 2005). It could explain the presence of H. natator and S. coahuilensis in high latitudes during the Late Cretaceous–Paleogene, as there were warm currents on North America's eastern coast (Jouve 2021). On the other hand, colder currents have been predicted near the European islands (Pucéat et al. 2005; Herman and Spicer 2010; Herman 2013). These cold currents may have excluded tethysuchians from Europe, as there are no consensual occurrences of this clade in this region during the end of the Cretaceous–Paleogene.

Figure 4. Distribution map of tethysuchians from the (A) pre- and (B) post-OAE 2 (second Oceanic Anoxic Event) faunae. The red polygon shows the repartition without Sabinosuchus coahuilensis and Hyposaurus natator. Map generated from the Paleobiology Database.

PLR analyses using Jouve's (2021) topologies showed that the post-OAE 2 fauna was more prone to brevirostry than the pre-OAE 2 one (Table 2). These cases of snout reduction have been described in dyrosauroids, especially during the Paleogene (Jouve et al. 2005a, 2021; Hastings et al. 2010). SP and shape classification in crocodyliforms have been widely discussed in the literature, with proposals of differing categories for various clades. However, a consensus emerges, suggesting that longer and thinner snouts are generally associated with a mainly piscivorous diet and shorter snouts are generally associated with a more durophageous diet (Brochu 2001; Pierce et al. 2009; Drumheller and Wilberg 2020). Early dyrosauroids such as the opportunistic predator Elosuchus have a much longer snout than Cenozoic forms such as Chenanisuchus lateroculi Jouve et al., 2005a, Anthracosuchus balrogus Hastings et al., 2015 and Rodeosuchus machukiru Jouve et al., 2021. This pattern is absent when using Sachs et al.'s (2021) topologies. This might, at least partly, be explained by the fact that 11 species included in Jouve (2021) are not present in Sachs et al. (2021). Indeed, among these missing species are Terminonaris browni, a longirostrine pre-OAE 2 pholidosaurid; Dorbignysuchus niatu Jouve et al., 2020, a brevirostrine (i.e., short-snouted) post-K/Pg dyrosaurid; and many post-K/Pg dyrosaurids that have a mesorostrine (i.e., a medium-sized) snout. The inclusion of these 11 species in Sachs et al.'s (2021) sample would be of interest for testing whether phylogeny, sample, or both are affecting the results. We note that Sachs et al.'s (2021) matrix is designed to test crocodyliform relationships and not intraclade relationships. In contrast, Jouve (2021) provides a matrix designed for Tethysuchia. Different statistical results between topologies may also be caused by these differing approaches. We consider that a significant variation in snout length after a crisis may indicate a selective extinction of a particular diet and/or diversification caused by character displacement, both being characteristic of niche partitioning (Brown and Wilson 1956). Longirostrine Tethysuchia are still very abundant after K/Pg, with species such as Atlantosuchus coupatezi Buffetaut and Wouters, 1979 and Luciasuchus lurusinqa Jouve et al., 2021, among others. Furthermore, brevirostrine Tethysuchia are only known after the K/Pg crisis, which marks the extinction of mosasaurs and plesiosaurs (Bardet 1995; Jouve et al. 2008b, 2021; Hastings et al. 2015). As mentioned earlier (see “A Differential and Phylogenetically Structured Response to Biotic Crises”), dyrosaurids may have taken the mosasaur's ecological position after the extinction of the latter. Colonization of now-empty environments may have allowed cases of niche partitioning. Niche partitioning has been described in thalattosuchians (De Andrade et al. 2010), eusuchians (Hastings and Hellmund 2017), and marine Mesozoic squamates (Bardet 2012; Bardet et al. 2015). Such a pattern is also present in dyrosaurids. Indeed, Paleogene dyrosaurid-bearing formations often include longirostrine, brevirostrine, and mesorostrine forms, each associated with a specific diet (piscivorous, durophagous, and generalist, respectively). We observe such a pattern for the formations of the Ouled Abdoun Basin (Paleocene–Ypresian) in Morocco (Bardet et al. 2010), the Cerrejon Formation (Paleocene) in Colombia (Hastings et al. 2015), and the Santa Lucia Formation (Paleocene) in Bolivia (Jouve et al. 2021).

A Trend toward Larger Body Sizes and Warm Climates?

A relationship between body size and temperature has already been tested for tethysuchians, yielding contrasting results depending on how body size is measured (Godoy et al. 2019). However, these analyses excluded PNC and used ordinary least squares and GLS instead of phylogenetic comparative methods. A significant relationship has been found for the crocodilian crown-group using the same methodology (Godoy et al. 2019; Godoy and Turner 2020), suggesting that larger body sizes are associated with cooler climates. Similar results were found when only extant crocodylomorphs were analyzed (Lakin et al. 2020). However, these relationships were tested using GLMs and not PGLS (Lakin et al. 2020). By contrast, the PGLS performed here shows that tethysuchians were larger in warmer climates (Table 3). The subsequently created paleotemperature-influenced model has a better linear fit than the null model that postulates no correlation (Table 4).

Tethysuchia were probably ectotherms (Faure-Brac et al. 2021). Ectothermic organisms are known to have a cyclic growth linked to seasonality that can be recorded in bone histology. Indeed, we can observe periods of rapid growth (zones) and periods of slow (annuli) or arrested (lines) growth (De Buffrénil and Quilhac 2021). Such a pattern has been identified in extant crocodilians (Hutton, 1987) and the crocodylomorph fossil record (Castanet et al. 1977; De Buffrénil et al. 2021). Various histological sections of tethysuchians show the presence of lines of arrested growth (Andrade and Sayão 2014; De Buffrénil et al. 2021; Faure-Brac et al. 2021). Furthermore, some extant vertebrate ectotherms have been known to have a preferential season of growth during warm periods and hence have larger sizes in warmer climates (Hjernquist et al. 2012). Therefore it is possible that Tethysuchia living in warmer temperatures had a longer preferred growth season. With a longer growth season, they may have grown larger, thus explaining the results shown by the PGLS. However, extant crocodilians, which are also ectothermic, have been noted to follow the opposite pattern (Godoy and Turner 2020; Lakin et al. 2020). Both groups share an overall similar thermophysiology; hence, this difference in body-size distribution cannot be explained by the thermometabolism. However, temperature tolerance differences have been noted to exist among extant crocodilians and have been suggested for extinct crocodilians (Jouve et al. 2019) or between dyrosaurids and gavialoids (Jouve 2021), the latter being present in warmer climatic zones than dyrosaurids. Therefore, there may be a different response in growth to paleotemperature for tethysuchians compared with modern crocodilians. Further exploration is required to identify the origins of these differences. We also tested the correlation for dyrosauroids and pholidosaurids separately. We observe a similar correlation for pholidosaurids if we follow Jouve's (2021) topologies. However, because the pholidosaurid sample is very small (n < 8) and the null model generally has a better linear fit than the paleotemperature-influenced one, we consider that there is not enough statistical evidence to prove any correlation for pholidosaurids separately. Finally, the null model remains the better linear fit, and no correlation is found for dyosauroids after Vectisuchus leptognathus is excluded from the analysis. Vectisuchus leptognathus is a clear outlier in our sample, as it caused heteroscedasticity in PGLS analyses if it was not excluded. Furthermore, it is both the smallest known tethysuchian and the one that lives in the coldest environment (Frakes et al. 1994; Salisbury and Naish 2011). Its unique specimen has been found in the Upper Wessex Formation (Barremian) of England, and it lived among many other crocodylomorphs, mostly goniopholidids (Salisbury and Naish 2011). Its small size may result from niche partitioning with these other crocodylomorphs: smaller-sized species do not feed from the same resources as larger species. Such a pattern has been observed in Metriorhynchidae: species with similar ecologies have a wide size range and different prey (Young et al. 2011). However, V. leptognathus lived at a high paleolatitude, which may explain its low inferred paleotemperature (see Supplementary File 1). Such a temperature cannot be explained by paleolatitude alone, as it is not the highest paleolatitude of our sample (see Supplementary File 1 for details). Furthermore, during the Barremian, many crocodylomorphs were living at similarly high paleolatitudes (Salisbury and Naish 2011). The notably low paleotemperature we inferred for V. leptognathus may result from cold environmental conditions during the Barremian. Indeed, its stratigraphic extent occurs during the Tithonian–Early Barremian cool interval, which is the coldest period in the Mesozoic (Scotese et al. 2021). Therefore, a combination of a high latitude during a notably cold period explains its low paleotemperature, and niche partitioning may explain its small size. Finally, no significant correlation is found if the SL estimations of Meridiosaurus vallisparadisi and S. coahuilensis are excluded. Hence, all of these elements suggest that our results must be treated with caution, because changing the interpretation for one or two species heavily affects the results. This caution is strengthened by the results of variation partitioning that suggest different results, indicating that paleotemperature and phylogeny had a nonsignificant influence on log-transformed SL variation. According to the analyses, most of the variation remains unexplained. These differing results may be explained by the relative scarcity of SL data. Indeed, some species were excluded from both variation partitioning and PGLS analyses because they had no available SL. These missing values may have impacted the results differently depending on the methodology used. We note that variation partitioning may show that, apart from temperature, various other factors may explain log-transformed SL variation. One of these factors may be dietary differences. Indeed, dietary shifts have been shown to explain body-size variation in Canidae and Notosuchia (Van Valkenburgh et al. 2004; Aubier et al. 2023). Another possible component may be species competitiveness: species with a relatively similar ecology may limit competition for the same resources because of their larger range of body sizes. Therefore, they feed on different prey while having a similar ecology. This pattern has been suggested in metriorhynchid crocodylomorphs (Young et al. 2011). However, we lack tethysuchian fossil data to further test both of these assumptions. Finally, other poorly understood or yet undiscovered biological factors may explain more log-transformed SL variation.

Conclusion

Tethysuchians crossed two biological crises, the OAE 2 and K/Pg, during their evolutionary history. Extinction was phylogenetically structured in both of them. These crises had differential effects on paleobiodiversity: first, the OAE 2 crisis was followed by a turnover of tethysuchian diversity with a pholidosaurid-dominated fauna replaced by a dyrosaurid-dominated one. Second, the K/Pg crisis was followed by increased biodiversity, especially regarding dyrosaurids, which remained high until the Eocene. Post-OAE 2 tethysuchians lived in warmer environments than the pre-OAE 2 fauna thanks to an overall global warming, possibly combined with a more restricted lower-latitude extension. The possible colonization of new ecological niches, likely left vacant by the extinction of mosasaurs and plesiosaurs, may also have allowed morphological diversification regarding the SP and shape in the same formations. This niche partitioning is shown by the co-occurrences of multiple tethysuchians associated with diverse diets. Finally, unlike other studies (Godoy et al. 2019; Godoy and Turner 2020; Lakin et al. 2020), we found a positive correlation between body length (using the log-transformed SL as proxy) and temperature. These results may be explained by the difference in a preferential season of growth duration. Nevertheless, these results must be treated with caution, as the fossil record of tethysuchians is scarce, most notably during the Late Cretaceous. These results also depend heavily on the size estimations from two taxa in our sample, suggesting that the SL sample may be an issue. Finally, variation partitioning suggested that other factors may explain body-size variation in Tethysuchia. Therefore, further exploration is required to uncover body-size evolution in Tethysuchia.