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. Reference Irmis, Nesbitt, Sues, Nesbitt, Desojo and Irmis2013) and occupied many ecological niches (Wilberg et al. Reference Wilberg, Turner and Brochu2019). It crossed major extinction events such as the Triassic/Jurassic (T/J) crisis, after which it radiated (Toljagić and Butler Reference Toljagić and Butler2013; Bronzati et al. Reference Bronzati, Montefeltro and Langer2015), and the Cretaceous/Paleogene (K/Pg) crisis. Its diversity declined during the Cenozoic, probably due to climate cooling (Markwick Reference Markwick1998) 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 Reference Grigg and Kirshner2015).
Among crocodylomorphs, Tethysuchia Buffetaut, Reference Buffetaut1982 is a group of semi-aquatic freshwater and marine neosuchians (Andrade and Sayão Reference Andrade and Sayão2014) that extended from the Kimmeridgian to the Bartonian (Jouve et al. Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021). They were probably ectothermic animals (Faure-Brac et al. Reference Faure-Brac, Amiot, De Muizon, Cubo and Lécuyer2021). While tethysuchians are ancestrally freshwater organisms (Martin et al. Reference Martin, Lauprasert, Buffetaut, Liard and Suteethorn2014b), independent events of colonization of the marine environment have been reported (Jouve et al. Reference Jouve, Bouya and Amaghzaz2005a,Reference Jouve, Iarochene, Bouya and Amaghzazb; Hua et al. Reference Hua, Buffetaut, Legall and Rogron2007; Wilberg et al. Reference Wilberg, Turner and Brochu2019; Jouve Reference Jouve2021). This group is composed of two clades (Jouve Reference Jouve2021): Pholidosauridae Zittel and Eastman Reference Zittel and Eastman1902, which extended from the Kimmeridgian (Mones Reference Mones1980) to the Danian (Jouve and Jalil Reference Jouve and Jalil2020) and Dyrosauroidea Jouve et al., Reference Jouve2021, which extended from the Barremian (Buffetaut and Hutt Reference Buffetaut and Hutt1980) to the Bartonian (Buffetaut Reference Buffetaut1978). 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 Reference Turgeon and Creaser2008). These nutrients increased primary production, leading to greater oxidation of organic matter, stripping the ocean of oxygen, causing anoxia (Bralower Reference Bralower2008; Turgeon and Creaser Reference Turgeon and Creaser2008). Coupled with this, an important greenhouse effect may have been generated by the volcanic CO2, leading to a stratified ocean that hampered oxygen delivery to deep waters (Bralower Reference Bralower2008; Turgeon and Creaser Reference Turgeon and Creaser2008). This event was linked to the extinction of ichthyosaurs (Fischer et al. Reference Fischer, Bardet, Benson, Arkhangelsky and Friedman2016). 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. Reference Hildebrand, Penfield, Kring, Pilkington, Camargo, Jacobsen and Boynton1991) and important volcanism in the Deccan Traps (Courtillot Reference Courtillot1990). The timing and importance of each event remain heavily discussed (Schoene et al. Reference Schoene, Eddy, Samperton, Keller, Keller, Adatte and Khadr2019; Sprain et al. Reference Sprain, Renne, Vanderkluysen, Pande, Self and Mittal2019). The K/Pg crisis was linked to the extinction of non-avian dinosaurs (Novacek Reference Novacek1999); pterosaurs (Barrett et al. Reference Barrett, Butler, Edwards, Milner, Hone and Buffetaut2008); and many marine reptiles, including mosasaurs and plesiosaurs (Bardet Reference Bardet1995).
Numerous studies have tried to identify factors driving crocodylomorph evolution. Martin et al. (Reference Martin, Amiot, Lécuyer and Benton2014a) 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. (Reference Jouve, Mennecart, Douteau and Jalil2017) questioned the reliability of these results, stating that they were heavily affected by minor taxonomic updates. Mannion et al. (Reference Mannion, Benson, Carrano, Tennant, Judd and Butler2015) found that diversification patterns for crocodylomorphs tracked environmental variations, but contrary to Martin et al. (Reference Martin, Amiot, Lécuyer and Benton2014a), no significant correlation between diversity and temperature was found for marine taxa. Jouve and Jalil (Reference Jouve and Jalil2020) 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. (Reference Bronzati, Montefeltro and Langer2015) 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 (Reference Jouve2021) found an important diversification event for longirostrine (i.e., long-snouted) crocodylomorphs following the K/Pg crisis, especially regarding dyrosaurid tethysuchians. Godoy et al. (Reference Godoy, Benson, Bronzati and Butler2019) 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 (Reference Stockdale and Benton2021) found a significant correlation between mean body size and paleotemperature for crocodylomorphs. However, Benson et al. (Reference Benson, Godoy, Bronzati, Butler and Gearty2022) 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. Reference Wilberg, Turner and Brochu2019) 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. Reference Aubier, Jouve, Schnyder and Cubo2023). 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. Reference Aubier, Jouve, Schnyder and Cubo2023).
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 Reference Jouve and Jalil2020; Jouve Reference Jouve2021), 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. Reference Godoy, Benson, Bronzati and Butler2019), 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. Reference Godoy, Benson, Bronzati and Butler2019; Benson et al. Reference Benson, Godoy, Bronzati, Butler and Gearty2022). However, mixed results were obtained when focusing on Tethysuchia (Godoy et al. Reference Godoy, Benson, Bronzati and Butler2019). 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 Reference Buffetaut and Hutt1980; Hastings et al. Reference Hastings, Bloch and Jaramillo2011; Jouve et al. Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021), 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. Reference Sereno, Larsson, Sidor and Gado2001; Jouve et al. Reference Jouve, Iarochène, Bouya and Amaghzaz2006), 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. (Reference O'Brien, Lynch, Vliet, Brueggen, Erickson and Gignac2019) 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, Reference De Lapparent De Broin and Taquet1966. 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. Reference Godoy, Benson, Bronzati and Butler2019; Aubier et al. Reference Aubier, Jouve, Schnyder and Cubo2023). 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. Reference Young, Bell, De Andrade and Brusatte2011; Stockdale and Benton Reference Stockdale and Benton2021). Most recently, Stockdale and Benton (Reference Stockdale and Benton2021) have applied principal component analysis (PCA) using various body size indicators to distinguish independent components linked to body size. However, Benson et al. (Reference Benson, Godoy, Bronzati, Butler and Gearty2022) 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. Reference Farlow, Hurlburt, Elsey, Britton and Langston2005; Vandermark et al. Reference Vandermark, Tarduno and Brinkman2007). However, most of these equations are based on extant crocodilians, particularly Alligator mississipiensis (Daudin, Reference Daudin1802), 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. (Reference Young, Bell, De Andrade and Brusatte2011) 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, Reference Mones1980 and Sabinosuchus coahuilensis Shiller et al., Reference Shiller, Porras-Muzquiz and Lehman2016, 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. Reference Fortier, Perea, Schultz, Pol and Larsson2011) and on the mandible length (Shiller et al. Reference Shiller, Porras-Muzquiz and Lehman2016), 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. (Reference Benson, Godoy, Bronzati, Butler and Gearty2022).
We gathered the paleoenvironments of analyzed taxa using Jouve (Reference Jouve2021). It can be hypothesized that some Tethysuchia could possibly move between fresh and salt water, like some modern crocodilians (Grigg and Kirshner Reference Grigg and Kirshner2015). 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 Reference Grigg and Kirshner2015). 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 Reference Jouve2021). Therefore, we considered that the various specimens were buried in their preferred environments and were coded accordingly. Dakotasuchus kingi Mehl, Reference Mehl1941 is the only exception, as it was found in one marine and two freshwater localities (Jouve Reference Jouve2021). 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 Reference Jouve and Jalil2020). Its state of preservation is striking compared with the other Tethysuchia from this formation (Jouve et al. Reference Jouve, Bouya and Amaghzaz2005a,Reference Jouve, Iarochene, Bouya and Amaghzazb, Reference Jouve, Iarochène, Bouya and Amaghzaz2006, Reference Jouve, Bouya and Amaghzaz2008b). Such a difference in preservation suggests transport from a freshwater to a marine locality (Jouve Reference Jouve2021). 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 Reference Newton and Mudge2003; Pouech et al. Reference Pouech, Amiot, Lécuyer, Mazin, Martineau and Fourel2014). Furthermore, tethysuchians spent time out of the water, as they had a semi-aquatic lifestyle (Andrade and Sayão Reference Andrade and Sayão2014). 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. Reference Frakes, Probst and Ludwig1994; Amiot et al. Reference Amiot, Lécuyer, Buffetaut, Fluteau, Legendre and Martineau2004; Pouech et al. Reference Pouech, Amiot, Lécuyer, Mazin, Martineau and Fourel2014; Alberti et al. Reference Alberti, Fürsich, Abdelhady and Andersen2017; Zhang et al. Reference Zhang, Hay, Wang and Gu2019; Laugié et al. Reference Laugié, Donnadieu, Ladant, Bopp and Raisson2020; see details in Supplementary File 1). Note that no extensive latitudinal temperature gradient study has been made for the Danian. As δ18O levels between the Maastrichtian and Selandian are rather similar (Prokoph et al. Reference Prokoph, Shields and Veizer2008), 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 (Reference Jouve2021) 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., Reference Salih, Evans, Bussert, Klein and Müller2022, which has been recovered as the second-earliest diverging dyrosaurid (Salih et al. Reference Salih, Evans, Bussert, Klein and Müller2022). Dakotasuchus kingi has a debated phylogenetic position (Jouve and Jalil Reference Jouve and Jalil2020). It is considered to be part of the clade including Terminonaris robusta (Mook, Reference Mook1934), Terminonaris browni (Osborn, Reference Osborn1904), and MHNM-kh01, a poorly preserved Danian specimen (Jouve and Jalil Reference Jouve and Jalil2020) or the sister species of Pholidosaurus Meyer, Reference Meyer1841. 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. (Reference Sachs, Young, Abel and Mallison2021), 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, Reference Buffetaut and Hutt1980 and Elosuchus De Lapparent De Broin, Reference De Lapparent De Broin2002 are considered pholidosaurids. Sachs et al. (Reference Sachs, Young, Abel and Mallison2021) yielded two topologies: the first one retrieves Pholidosaurus schaumburgensis Meyer, Reference Meyer1841 in a clade with Oceanosuchus boecensis Hua et al., Reference Hua, Buffetaut, Legall and Rogron2007 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, Reference Buffetaut and Ingavat1980) and Sarcosuchus De Lapparent De Broin and Taquet, Reference De Lapparent De Broin and Taquet1966. 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, Reference Owen1884) and Anglosuchus laticeps (Owen, Reference Owen1884) are considered Bathonian. However, their ages remain doubtful, and they closely resemble Pholidosaurus purbeckensis (Mansel-Pleydell, Reference Mansel-Pleydell1888), so they may be synonyms of the latter (Jouve and Jalil Reference Jouve and Jalil2020). The pholidosaurids Pholidosaurus milwardi Roxo, Reference Roxo1929 and Pholidosaurus meyeri Dunker, Reference Dunker1843 and the dyrosaurids Tilemsisuchus lavocati Buffetaut, Reference Buffetaut and Michaud1980, Congosaurus compressus (Buffetaut, Reference Buffetaut and Michaud1980), and Rhabdognathus acutirostris Bergounioux, Reference Bergounioux1955 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 Reference Buffetaut and Suteethorn2007). However, recent palynology studies suggested a Berriasian age (Racey and Goodall Reference Racey and Goodall2009). Therefore, we considered C. thailandicus to be of Berriasian age. Hyposaurus natator Troxell, Reference Troxell1925 is noted as being Maastrichtian in the PBDB (Cope, Reference Cope1866; Marsh, Reference Marsh1870). 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. Reference Landman, Johnson, Garb, Edwards and Kyte2007; Wiest et al. Reference Wiest, Buynevich, Grandstaff, Maza and Lacovara2016). Therefore, in our analyses, we considered it to be Danian.
Topologies were dated using the timePaleoPhy function on the paleotree package (Bapst Reference Bapst2012) 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 Reference Soul and Wright2021). A minimal branch length can be selected to prevent ZLBs (Laurin Reference Laurin2004; Wang and Lloyd Reference Wang and Lloyd2016). 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 Reference Stockdale and Benton2021). 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 Reference Fritz and Purvis2010). 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 Reference Fritz and Purvis2010; Yessoufou et al. Reference Yessoufou, Daru and Davies2012; Fontana et al. Reference Fontana, Furtado, Zanella, Debastiani and Hartz2021) or extinction selectivity in the fossil record (Allen et al. Reference Allen, Stubbs, Benton and Puttick2019; Wilke et al. Reference Wilke, Hauffe, Jovanovska, Cvetkoska, Donders, Ekschmitt and Francke2020; Aubier et al. Reference Aubier, Jouve, Schnyder and Cubo2023). We used the phylo.d function of the caper package (Orme et al. Reference Orme, Freckleton, Thomas, Petzoldt, Fritz, Isaac and Pearse2013) 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 (Reference Fritz and Purvis2010). 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 Reference Ives and Garland2010) 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é Reference 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 Reference Felsenstein1985), we included the dated trees (see “Supertree”). PLR has been used to infer the probability of endothermy in tetrapods (Cubo et al. Reference Cubo, Aubier, Faure-Brac, Martet, Pellarin, Pelletan and Sena2023; Faure-Brac et al. Reference Faure-Brac, Woodward, Aubier and Cubo2024) and the probability of survival after the K/Pg crisis in Notosuchia (Aubier et al. Reference Aubier, Jouve, Schnyder and Cubo2023), 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 Reference Ackerly2003; Cooper et al. Reference Cooper, Jetz and Freckleton2010). 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 Reference Grafen and Hamilton1989) to test whether temperature affected log-transformed SL. We used the pgls function from the caper R package (Orme Reference Orme, Freckleton, Thomas, Petzoldt, Fritz, Isaac and Pearse2013) 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 Reference Shapiro and Wilk1965) 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 Reference Breusch and Pagan1979). Homoscedasticity was not respected if V. leptognathus, which had a very short skull (Salisbury and Naish Reference Salisbury, Naish and Batten2011) and lived in very cold temperatures (Frakes et al. Reference Frakes, Probst and Ludwig1994), 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 Reference Mazerolle2013) 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. Reference Borcard, Legendre and Drapeau1992). This method has been further developed to consider phylogeny as a component using a matrix of principal coordinates representing phylogeny (Desdevises et al. Reference Desdevises, Legendre, Azouzi and Morand2003; Peres-Neto et al. Reference Peres-Neto, Legendre, Dray and Borcard2006; Montes et al. Reference Montes, Le Roy, Perret, De Buffrénil, Castanet and Cubo2007; Piras et al. Reference Piras, Teresi, Buscalioni and Cubo2009; Sakamoto et al. Reference Sakamoto, Lloyd and Benton2010). We used the varpart function from the vegan R package (Dixon Reference Dixon2003) 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 (Reference Jouve2021) 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 (Reference Jouve2021) 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 (Reference Jouve2021) 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.
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; 1Nonnormal, p = 0.038; 2Nonhomogenous, 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 (Reference Jouve2021) 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 2, 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 2 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 2 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 Reference Jouve and Jalil2020; Jouve Reference Jouve2021; Fig. 1, Supplementary File 2). This period corresponds to the highest temperature and sea level of the Mesozoic (Vérard et al. Reference Vérard, Hochard, Baumgartner, Stampfli and Liu2015; Scotese et al. Reference Scotese, Song, Mills and Van Der Meer2021), which may explain the important tethysuchian fossil record, because high sea level has long been considered a factor of enhanced diversity (Martin et al. Reference Martin, Amiot, Lécuyer and Benton2014a; Mannion et al. Reference Mannion, Benson, Carrano, Tennant, Judd and Butler2015; Tennant et al. Reference Tennant, Mannion and Upchurch2016). During the Cenomanian/Turonian transition, Tethysuchia experienced a major diversity drop corresponding to OAE 2 (Jouve and Jalil Reference Jouve and Jalil2020; Fig. 1, Supplementary File 2). Because half of the tethysuchians at the time were marine (Jouve Reference Jouve2021), they were probably heavily affected by this event, which was also linked to the extinction of ichthyosaurs (Fischer et al. Reference Fischer, Bardet, Benson, Arkhangelsky and Friedman2016) and the diversification of mosasaurs (Bardet Reference Bardet1995). These patterns suggest an important marine faunal turnover previously mentioned in the literature (Kauffman Reference Kauffman and Council1995; Wan et al. Reference Wan, Wignall and Zhao2003; Caron et al. Reference Caron, Dall'Agnolo, Accarie, Barrera, Kauffman, Amédro and Robaszynski2006; Monnet Reference Monnet2009). 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, Reference De Stefano1903 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 Reference Jouve and Jalil2020; 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 Reference Buffetaut1974; Meunier and Larsson Reference Meunier and Larsson2018). The Coniacian–Santonian interval coincides with a marine regression (Jouve and Jalil Reference Jouve and Jalil2020), 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. Reference Puértolas-Pascual, Blanco, Brochu and Canudo2016; Meunier and Larsson Reference Meunier and Larsson2018). 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 Reference Halstead1975; Shiller et al. Reference Shiller, Porras-Muzquiz and Lehman2016; Jouve and Jalil Reference Jouve and Jalil2020; Salih et al. Reference Salih, Evans, Bussert, Klein and Müller2022). Most Late Cretaceous tethysuchians lived in freshwater environments (Jouve Reference Jouve2021). 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. Reference Robertson, Lewis, Sheehan and Toon2013).
Following the K/Pg crisis, an explosion in diversity occurs (Jouve Reference Jouve2021; Fig. 1, Supplementary File 2). Most Cenozoic tethysuchians lived in marine environments (Jouve Reference Jouve2021). 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. Reference Barbosa, Kellner and Viana2008; Jouve et al. Reference Jouve, Bardet, Jalil, Suberbiola, Bouya and Amaghzaz2008a,Reference Jouve, Bouya and Amaghzazb; Bardet et al. Reference Bardet, Gheerbrant, Noubhani, Cappetta, Jouve, Bourdon, Suberbiola, Jalil, Vincent and Houssaye2017; Jouve Reference Jouve2021). 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, Reference Nopcsa1928 (Fig. 1, Supplementary File 2) that heavily colonized the marine realm during the Paleogene (Jouve Reference Jouve2021).
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 Reference Jouve and Jalil2020), 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 Reference Jouve2021). The factors explaining their extinction are uncertain (Amoudji et al. Reference Amoudji, Guinot, Hautier, Kassegne, Chabrol, Charruault, Johnson, Sarr, Da Costa and Martin2021). The Bartonian coincides with the beginning of the late Eocene–Oligocene cooling (Scotese et al. Reference Scotese, Song, Mills and Van Der Meer2021), which may have impacted tethysuchians (Jouve Reference Jouve2021). Another hypothesis involving competition with new predators, including cetaceans, was mentioned by Hastings (Reference Hastings2012). It was considered unlikely by Martin et al. (Reference Martin, Amiot, Lécuyer and Benton2014a), but was still mentioned by Stubbs et al. (Reference Stubbs, Pierce, Elsler, Anderson, Rayfield and Benton2021). These hypotheses remain to be tested. Most recently, Scott and Anderson (Reference Scott and Anderson2023) 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 (Reference Jouve2021) 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. Reference Scotese, Song, Mills and Van Der Meer2021). 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. Reference Scotese, Song, Mills and Van Der Meer2021). 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 Reference Jouve2021). Indeed, a proto–Gulf Stream has been suggested in literature (Watkins and Self-Trail Reference Watkins and Self-Trail2005). 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 Reference Jouve2021). On the other hand, colder currents have been predicted near the European islands (Pucéat et al. Reference Pucéat, Lécuyer and Reisberg2005; Herman and Spicer Reference Herman and Spicer2010; Herman Reference Herman2013). 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 (Reference Jouve2021) 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. Reference Jouve, Bouya and Amaghzaz2005a, Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021; Hastings et al. Reference Hastings, Bloch, Cadena and Jaramillo2010). 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 Reference Brochu2001; Pierce et al. Reference Pierce, Angielczyk and Rayfield2009; Drumheller and Wilberg Reference Drumheller and Wilberg2020). Early dyrosauroids such as the opportunistic predator Elosuchus have a much longer snout than Cenozoic forms such as Chenanisuchus lateroculi Jouve et al., Reference Jouve, Bouya and Amaghzaz2005a, Anthracosuchus balrogus Hastings et al., Reference Hastings, Bloch and Jaramillo2015 and Rodeosuchus machukiru Jouve et al., Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021. 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 (Reference Jouve2021) are not present in Sachs et al. (Reference Sachs, Young, Abel and Mallison2021). Indeed, among these missing species are Terminonaris browni, a longirostrine pre-OAE 2 pholidosaurid; Dorbignysuchus niatu Jouve et al., Reference Jouve and Jalil2020, 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 (Reference Jouve2021) 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 Reference Brown and Wilson1956). Longirostrine Tethysuchia are still very abundant after K/Pg, with species such as Atlantosuchus coupatezi Buffetaut and Wouters, Reference Buffetaut and Wouters1979 and Luciasuchus lurusinqa Jouve et al., Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021, among others. Furthermore, brevirostrine Tethysuchia are only known after the K/Pg crisis, which marks the extinction of mosasaurs and plesiosaurs (Bardet Reference Bardet1995; Jouve et al. Reference Jouve, Bouya and Amaghzaz2008b, Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021; Hastings et al. Reference Hastings, Bloch and Jaramillo2015). 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. Reference De Andrade, Young, Desojo and Brusatte2010), eusuchians (Hastings and Hellmund Reference Hastings and Hellmund2017), and marine Mesozoic squamates (Bardet Reference Bardet2012; Bardet et al. Reference Bardet, Houssaye, Vincent, Suberbiola, Amaghzaz, Jourani and Meslouh2015). 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. Reference Bardet, Suberbiola, Jouve, Bourdon, Vincent, Houssaye, Rage, Jalil, Bouya and Amaghzaz2010), the Cerrejon Formation (Paleocene) in Colombia (Hastings et al. Reference Hastings, Bloch and Jaramillo2015), and the Santa Lucia Formation (Paleocene) in Bolivia (Jouve et al. Reference Jouve, De Muizon, Cespedes-Paz, Sossa-Soruco and Knoll2021).
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. Reference Godoy, Benson, Bronzati and Butler2019). 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. Reference Godoy, Benson, Bronzati and Butler2019; Godoy and Turner Reference Godoy and Turner2020), suggesting that larger body sizes are associated with cooler climates. Similar results were found when only extant crocodylomorphs were analyzed (Lakin et al. Reference Lakin, Barrett, Stevenson, Thomas and Wills2020). However, these relationships were tested using GLMs and not PGLS (Lakin et al. Reference Lakin, Barrett, Stevenson, Thomas and Wills2020). 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. Reference Faure-Brac, Amiot, De Muizon, Cubo and Lécuyer2021). 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 Reference De Buffrénil, Quilhac, Buffrénil, De Ricqlès, Zylberberg and Padian2021). Such a pattern has been identified in extant crocodilians (Hutton, Reference Hutton1987) and the crocodylomorph fossil record (Castanet et al. Reference Castanet, Meunier and De Ricqlès1977; De Buffrénil et al. Reference De Buffrénil, Laurin, Jouve, Buffrénil, De Ricqlès, Zylberberg and Padian2021). Various histological sections of tethysuchians show the presence of lines of arrested growth (Andrade and Sayão Reference Andrade and Sayão2014; De Buffrénil et al. Reference De Buffrénil, Laurin, Jouve, Buffrénil, De Ricqlès, Zylberberg and Padian2021; Faure-Brac et al. Reference Faure-Brac, Amiot, De Muizon, Cubo and Lécuyer2021). 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. Reference Hjernquist, Söderman, Jönsson, Herczeg, Laurila and Merilä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 Reference Godoy and Turner2020; Lakin et al. Reference Lakin, Barrett, Stevenson, Thomas and Wills2020). 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. Reference Jouve, Khalloufi and Zouhri2019) or between dyrosaurids and gavialoids (Jouve Reference Jouve2021), 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 (Reference Jouve2021) 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. Reference Frakes, Probst and Ludwig1994; Salisbury and Naish Reference Salisbury, Naish and Batten2011). 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 Reference Salisbury, Naish and Batten2011). 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. Reference Young, Bell, De Andrade and Brusatte2011). 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 Reference Salisbury, Naish and Batten2011). 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. Reference Scotese, Song, Mills and Van Der Meer2021). 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. Reference Van Valkenburgh, Wang and Damuth2004; Aubier et al. Reference Aubier, Jouve, Schnyder and Cubo2023). 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. Reference Young, Bell, De Andrade and Brusatte2011). 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. Reference Godoy, Benson, Bronzati and Butler2019; Godoy and Turner Reference Godoy and Turner2020; Lakin et al. Reference Lakin, Barrett, Stevenson, Thomas and Wills2020), 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.
Acknowledgments
We thank N. Campione (University of New England), C. Brochu (University of Iowa), and two anonymous reviewers for their suggestions. We thank J. Schnyder (Sorbonne Université) for helping us with climatic data treatment. We thank J. Bardin (Sorbonne Université) for his statistical advice. We thank M. G. Faure-Brac (University of Oslo) and M. V. A. Sena (Sorbonne Université) for their help.
Competing Interests
The authors declare no competing interests.
Data Availability Statement
Data available from the Dryad and Zenodo Digital Repositories: https://doi.org/10.5061/dryad.bnzs7h4j3, https://doi.org/10.5281/zenodo.10548562.
Open access
