## Abstract

Notosuchia is a group of mostly terrestrial crocodyliforms. The presence of a prominent crest overhanging the acetabulum, slender straight-shafted long bones with muscular insertions close to the joints, and a stable knee joint suggests that they had an erect posture. This stance has been proposed to be linked to endothermy, because it is present in mammals and birds and contributes to the efficiency of their respiratory systems. However, a bone paleohistological study unexpectedly suggested that Notosuchia were ectothermic organisms. The thermophysiological status of Notosuchia deserves further analysis, because the methodology of the previous study can be improved. First, it was based on a relationship between red blood cell size and bone vascular canal diameter tested using 14 extant tetrapod species. Here we present evidence for this relationship using a more comprehensive sample of extant tetrapods (31 species). Moreover, contrary to previous results, bone cross-sectional area appears to be a significant explanatory variable (in addition to vascular canal diameter). Second, red blood cell size estimations were performed using phylogenetic eigenvector maps, and this method excludes a fraction of the phylogenetic information. This is because it generates a high number of eigenvectors requiring a selection procedure to compile a subset of them to avoid model overfitting. Here we inferred the thermophysiology of Notosuchia using phylogenetic logistic regressions, a method that overcomes this problem by including all of the phylogenetic information and a sample of 46 tetrapods. These analyses suggest that *Araripesuchus wegeneri*, *Armadillosuchus arrudai*, *Baurusuchus* sp., *Iberosuchus macrodon*, and *Stratiotosuchus maxhechti* were ectothermic organisms.

## Introduction

Notosuchia is a group of extinct, mostly terrestrial crocodyliforms. The presence of several morphological features suggests that they had an erect posture: the prominent crest overhanging the acetabulum observed in *Chimaerasuchus paradoxus* (Wu and Sues 1996), *Notosuchus terrestris* (Pol 2005), *Araripesuchus tsangatsangana* (Turner 2006), *Baurusuchus albertoi* (Nascimento and Zaher 2010), and *Stratiotosuchus maxhechti* (Riff and Kellner 2011); the straight-shafted long bones described in *Anatosuchus minor* and *Araripesuchus* spp. (Sereno and Larsson 2009); the slender limb bones with muscular insertions close to the joints reported in *Malawisuchus mwakasyungutiensis* (Gomani 1997); and the tight/stable knee joint shown in *Pissarrachampsa sera* (Godoy et al. 2016). Among extant tetrapods, only endotherms (mammals and birds) show an upright stance. This last feature has been proposed to be linked to endothermy, because it contributes to the efficiency of the respiratory system (Carrier 1987). Thus, according to this morphological evidence, we can reasonably hypothesize that Notosuchia were endothermic. However, Cubo et al. (2020) concluded that they were primitively ectothermic using two proxies: resting metabolic rate (RMR) and red blood cell size (RBC_{size}).

RMR is the minimal consumption of oxygen over time per unit of body mass measured under postabsorptive conditions during the period of normal activity of the daily cycle in resting, nonreproductive specimens (Andrews and Pough 1985; Montes et al. 2007). RMRs of extant endotherms are at least one order of magnitude higher than those of extant ectotherms of similar body mass, because the mechanisms of thermogenesis operating in the former are costly in terms of energy (Clarke and Pörtner 2010; Legendre and Davesne 2020). RMRs inferred by Cubo et al. (2020) for Notosuchia were significantly lower than the threshold separating ectotherms from endotherms.

Within extant tetrapods, RBC_{size} is lower in endotherms (mammals and birds) than in ectotherms (Amphibia, Squamata, Testudines, and Crocodylia) (Hartman and Lessler 1964; Snyder and Sheafor 1999; Soslau 2020). It has been suggested that the acquisition of lungs together with the subsequent evolution of the cardiovascular system was the driving force explaining the evolution of vertebrate RBC_{size} (Snyder and Sheafor 1999). In endotherms, thermogenetic mechanisms use a huge amount of oxygen, producing high RMRs. Considering that “Smaller capillaries [and smaller RBCs] are associated with increased potential for diffusive gas exchange” (Snyder and Sheafor 1999: 189), these features may have been positively selected in endotherms. Huttenlocker and Farmer (2017) found that RBC_{size} values are related to, and can be inferred from, bone vascular canal diameter. Cubo et al. (2020) inferred notosuchian RBC_{size} values using this last relationship and concluded that they were significantly higher than the threshold separating ectotherms from endotherms.

To sum up, both proxies (RMR and RBC_{size}) suggest that Notosuchia were ectothermic organisms. Considering that paleohistological evidence (suggesting low RMR, large RBC_{size}, and ectothermy) is not congruent with morphological evidence (suggesting an erect posture, cursoriality, and endothermy), the thermophysiological status (i.e., either ectothermic or endothermic) of Notosuchia deserves further analysis.

The approach used by Cubo et al. (2020) to perform these inferences can be improved in two ways. First, notosuchian thermophysiological status inferred using RBC_{size} is based on the quoted relationship between RBC_{size} and bone vascular canal diameter (Cubo et al. 2020). This relationship was tested by Huttenlocker and Farmer (2017) using a rather small sample size (14 extant tetrapod species). Here we tested this relationship using a more comprehensive sample of extant tetrapods (31 species) and phylogenetic generalized least-squares regression (PGLS). Second, RBC_{size} estimations were performed using phylogenetic eigenvector maps (PEMs), and this method excludes a fraction of phylogenetic information. This is because PEM generates a high number of eigenvectors (*n* − 1, with *n* being the number of terminal taxa analyzed), thus requiring a selection procedure to compile a subset of eigenvectors to avoid model overfitting (Guénard et al. 2013; Legendre et al. 2016). Here we inferred the thermophysiology of Notosuchia using phylogenetic logistic regression (PLR) (Ives and Garland 2010; Tung Ho and Ané 2014), a method that overcomes this problem, because it includes all (instead of a fraction) of the phylogenetic information.

## Material and Methods

Phylogenies in Figure 1 and Supplementary File 1 contain the tetrapod samples used in this study. Topologies were taken from Pyron and Wiens (2011) for amphibians; Meredith et al. (2011), Zurano et al. (2019), Kumar et al. (2013), and Upham et al. (2019) for mammals; Ast (2001) and Villa et al. (2018) for *Varanus*; Man et al. (2011) for crocodiles; Prum et al. (2015) for birds; and Pol et al. (2014) for Notosuchia. Both phylogenies were dated using Time Tree of Life (http://www.timetree.org). When the ages of two successive nodes collapsed, we arbitrarily added 1 Myr in between the more-inclusive and less-inclusive nodes to facilitate the graphic visualization of the topology. For Notosuchia, nodes were dated according to the last appearance datum (LAD) of the oldest fossil included in each node taken from the Paleobiology Database (https://paleobiodb.org). The age of the node Notosuchia (113 Myr) corresponds to the LAD of *Malawisuchus mwakasyungutiensis*. The age of the node *Armadillosuchus–Baurusuchus* (100.5 Myr) corresponds to the LAD of *Chimaerasuchus paradoxus.* The age of the node *Iberosuchus–Baurusuchus* (83.6 Myr) corresponds to the LAD of *Comahuesuchus brachybuccalis*, *Pehuenchesuchus enderi*, *Cynodontosuchus rothi*, and *Wargosuchus australis.* Finally, the age of the node *Stratiotosuchus–Baurusuchus* (66 Myr) corresponds to the LAD of these taxa. The latter (*Stratiotosuchus–Baurusuchus*), as well as *Armadillosuchus arrudai*, come from the Adamantina Formation, the age of which is still debated. We follow the hypothesis of a Campanian–Maastrichtian age proposed by some authors (e.g., Gobbo-Rodrigues et al. 1999; Batezelli 2017).

#### Testing the Relationship between RBC_{size} and Bone Vascular Canal Diameter Using PGLS

Supplementary File 1 contains the sample (31 species of extant tetrapods) and the phylogeny (topology and divergence times) used to test the relationships between the response variables (RBC_{width} and RBC_{area}) and the explanatory variables (femoral vascular canal diameter and femoral cross-sectional area including the medullary cavity). Thin sections of extant taxa are curated at the Vertebrate Hard Tissue Collection of the Museum national d'Histoire naturelle, Paris, and are available on request to the curator (D. Germain). RBC_{width} (defined as RBC minimum diameter) and RBC_{area} (either published values or values computed using maximum and minimum published diameters and assuming an ellipse) were taken from the literature (Supplementary File 2). Femoral vascular canal diameters (white arrowheads in Fig. 2) were computed as Can_{harmean} and Can_{min}, as defined by Huttenlocker and Farmer (2017). Can_{harmean}, Can_{min}, and femoral cross-sectional area were either quantified in this study or taken from Huttenlocker and Farmer (2017) (data available in Supplementary File 2).

The method of ordinary least-squares regression makes the assumption of no covariance between residuals obtained from the regression equation (i.e., the off-diagonals of the variance–covariance matrix are expected to contain zeros) (Symonds and Blomberg 2014). In analyses using interspecific data, this assumption is not verified because of the hierarchical, shared phylogenetic history among terminal taxa (i.e., closely related species are more similar than expected by chance). PGLS (Grafen 1989; Martins and Hansen 1996; Rohlf 2001; Symonds and Blomberg 2014) overcomes this problem by using a variance–covariance matrix in which off-diagonals correspond to the phylogenetic history shared by the two species under comparison. Symonds and Blomberg (2014) described PGLS as a “weighted regression” in which data points corresponding to closely related species are “downweighted.” We ran PGLS using the function *pgls* of the package caper (Orme et al. 2013) in R (R Development Core Team 2008).

#### Inferring the Thermophysiology of Notosuchia Using PLR

Figure 1 shows the phylogenetic relationships among extant taxa (46 species of tetrapods) used to construct the thermophysiology inference model (to infer the probability of being endothermic) and the extinct Notosuchia for which we performed paleobiological inferences. This model was constructed using femoral vascular canal diameter (Can_{harmean} and Can_{min}) and femoral cross-sectional area as explanatory variables. As noted earlier, femoral Can_{harmean}, Can_{min}, and femoral cross-sectional area were either quantified in this study or taken from Huttenlocker and Farmer (2017) (data available in Supplementary File 3). As before, thin sections of extant taxa are curated at the Vertebrate Hard Tissue Collection of the Museum national d'Histoire naturelle, Paris. Data for Notosuchia are taken from Cubo et al. (2020): *Araripesuchus wegeneri*, *Armadillosuchus arrudai*, *Baurusuchus* sp., *Iberosuchus macrodon*, and *Stratiotosuchus maxhechti.* Considering that the model was constructed using femora of extant species, we performed inferences only for those Notosuchia for which data for femora were available. Sample size was smaller for PGLS analyses, because data for RBC_{size} were not available for many species analyzed in PLR analyses. PLR is a generalized linear model explaining the probability of occurrence of the state “presence” of a binary response (dependent) variable (here the “presence of endothermy”) using continuous explanatory (independent) variables when residual variation of the former variable is phylogenetically structured (Ives and Garland 2010). The regression coefficients computed do account for phylogenetic correlation; when data are not phylogenetically structured, these coefficients are those of standard logistic regression (Ives and Garland 2010). PLR models contain two components. The first is controlled by parameters α (the transition rate) and μ (the asymptotic probability of being in state 1 [here the asymptotic probability of being endotherm]). Parameter α equals α_{1} + α_{0}; α_{1} being the probability that the response variable switches from 0 to 1 in each small time increment when it evolves up a phylogenetic tree, whereas the α_{0} parameter is the probability that it evolves from 1 to 0 (Ives and Garland 2010). The transition rate α is a measure of phylogenetic signal, because the larger the α, the quicker the evolutionary transitions and the lower the phylogenetic structure of data (Ives and Garland 2010). In the second component, the probability of occurrence of the state “presence of endothermy” is modeled using values of the independent (explanatory) variable (here bone vascular canal diameter). Parameters α and μ are estimated using an iterative process in which μ is estimated given α, using the quasi-likelihood function, and α is estimated given μ, using least squares until convergence (Ives and Garland 2010). Analyses were performed using the package phyloglm (Tung Ho and Ané 2014) in R (R Development Core Team 2008).

## Results

#### Testing the Relationship between RBC_{size} and Bone Vascular Canal Diameter Using PGLS

We used PGLS to test the relationships between RBC_{width} and RBC_{area} and the explanatory variables femoral vascular canal diameter (computed as Can_{harmean} and Can_{min}) and femoral cross-sectional area (data available in Supplementary File 2). Shapiro-Wilk normality tests showed that residuals of PGLS regression of RBC_{area} to Can_{harmean} + bone cross-sectional area and the regression RBC_{area} to Can_{min} + bone cross-sectional area do not follow a normal distribution (*p*-values of 0.0007557 and 0.001896, respectively). Thus, we performed a log transformation of all variables. After log transformation, residuals of all four PGLS regressions (RBC_{area} and RBC_{area} to the explanatory variables bone cross-sectional area and either Can_{min} or Can_{harmean}) do follow a normal distribution. All four of these PGLS regressions were significant and, in each regression, both explanatory variables (bone cross-sectional area and either Can_{min} or Can_{harmean}) were significant (Table 1).

#### Inferring the Thermophysiology of Notosuchia Using PLR

_{harmean}was used as the explanatory variable, we obtained a model with a transition rate α of 0.00144, an intercept estimate of 6.04 (

*p*-value = 0.004) and an estimate for the coefficient of Can

_{harmean}of −0.45 (

*p*-value = 0.001). The negative sign of the Can

_{harmean}coefficient indicates that the probability of being endothermic decreases as vascular canal diameter increases. Figure 3 shows the distribution of probabilities of being endothermic as a function of Can

_{harmean}variation. The corresponding equation is:

Ives and Garland (2010: p. 17) stated that “we assume that if μ_{i} <|$\;\bar{{\rm \mu }}$|, then trait Y will evolve toward 0; […] Conversely, if μ_{i} > |$\bar{{\rm \mu }}$|, then trait Y will evolve toward 1,” where |$\bar{{\rm \mu }}$| is the mean probability of being endotherm in our sample. Thus, we considered that |$\bar{{\rm \mu }}$| is the cutoff probability, so that an inferred probability higher than |$\bar{{\rm \mu }}$| would be evidence for endothermy. Conversely, a probability lower than |$\bar{{\rm \mu }}$| would be evidence for ectothermy. When Can_{harmean} was used as the explanatory variable, |$\bar{{\rm \mu }}$| = 0.59. To evaluate the predictive power of the model, we constructed a contingency table in which we inferred the thermometabolic regime of each extant species of the sample using its Can_{harmean}. Lines contain predictions (0, inferred ectothermy; 1, inferred endothermy) and columns contain true states (0, observed ectothermy; 1, observed endothermy):

The specificity (the ratio of quantity of true 1 inferred as 1 on the quantity of true 1; Sp = 27/(27 + 2)) equals 0.931. The sensitivity (the ratio of quantity of true 0 inferred as 0 on the quantity of true 0; Se = 14/(14 + 3)) equals 0.824. The classification error (the ratio (quantity of true 0 inferred as 1 + quantity of true 1 inferred as 0)/total; error = (3 + 2)/(14 + 2 + 3 + 27)) equals 0.109. This classification error is quite low, so we used the cutoff probability of 0.59 to perform paleobiological inferences using Can_{harmean} as the explanatory (predictor) variable.

_{min}was used as the explanatory variable, we obtained a model with a transition rate α of 0.00052, an intercept estimate of 2.58 (

*p*-value = 0.032), and an estimate for the coefficient of Can

_{min}of −0.49 (

*p*-value = 0.018). Figure 4 shows the distribution of probabilities of being endothermic as a function of Can

_{min}variation. The corresponding equation is:

Again we evaluated the quality of the model by constructing a contingency table in which we inferred the thermophysiological regime of each extant species of the sample using its Can_{min}. We used a cutoff probability of |$\bar{{\rm \mu }}$| = 0.32 (the mean probability of being endothermic in our sample), so that an inferred probability lower than 0.32 is evidence for ectothermy and a probability higher than 0.32 is evidence for endothermy:

The specificity (the ratio of quantity of true 1 inferred as 1 on the quantity of true 1; Sp = 20/(20 + 9)) equals 0.690. The sensitivity (the ratio of quantity of true 0 inferred as 0 on the quantity of true 0; Se = 15/(15 + 2)) equals 0.882. The classification error (the ratio (quantity of true 0 inferred as 1 + quantity of true 1 inferred as 0)/total; error = (2 + 9)/(15 + 9 + 2 + 20)) equals 0.239. This classification error is quite high, so we recomputed a new cutoff probability of 0.22 using the receiver operating characteristic curve. Then we constructed a new contingency table in which we inferred the thermophysiological regime of each extant species of the sample using its Can_{min} and considering that an inferred probability higher than 0.22 is evidence for endothermy:

With this contingency table, the classification error [the ratio (quantity of true 0 inferred as 1 + quantity of true 1 inferred as 0)/total; error = (2 + 1)/(15+ 1 + 2 + 28)] equals 0.0652. This classification error is quite low, so we used the cutoff probability of 0.22 to perform paleobiological inferences using Can_{min} as the explanatory (predictor) variable.

## Discussion

Notosuchia is an extremely diversified group of crocodyliforms. This diversity is particularly striking regarding their diet, suggesting that they occupied various ecological niches (Carvalho and Bertini 1999; Iori and Carvalho 2018). Uruguaysuchidae (Pol et al. 2014), the most basal notosuchians, of which *Araripesuchus wegeneri* (in our sample) is a representative, range from the Aptian (*Araripesuchus gomesii*) to the Maastrichtian (*Araripesuchus tsangatsangana*) (Price 1959; Turner 2006). Several species of this group have been inferred as being omnivorous, or even insectivorous, based on their dental complexity (Sereno and Larsson 2009; Soto et al. 2011; Nieto et al. 2021) and postcranial remains suggest that they had an erect posture (see “Introduction”). Uruguaysuchids were smaller than Sphagesauridae (Carvalho et al. 2010; Godoy et al. 2019). *Armadillosuchus* (also sampled by us) belongs to the large-bodied sphagesaurids group (Melstrom and Irmis 2019). The diagnosis of this clade is based on its peculiar dentition morphology (Price 1950). They show extremely complex manducatory systems, with evidence of “chewing” mechanisms, dental wear, and propalinal movements (e.g., see Ősi 2014; Iori and Carvalho 2018). The foraging abilities of some notosuchians, such as *Armadillosuchus*, *Mariliasuchus, or Malawisuchus*, to locate food or water have led some authors to propose the presence of burrowing habits (Gomani 1997; Nobre et al. 2008; Marinho and Carvalho 2009), a behavior that might play a role in thermoregulation (e.g., to search for a cooler shelter during dry periods, as in extant crocodilian species; Campos and Magnusson 2013). This behavior has also been proposed for the sebecosuchian *Baurusuchus salgadoensis* (Vasconcellos and Carvalho 2010). Sebecosuchia (to which *Iberosuchus* and *Stratiotosuchus*, also sampled by us, belong) were large predators with an erect posture (see “Introduction”) and cursorial abilities (Nascimento and Zaher 2010; Riff and Kellner 2011), feeding on large prey, including small sphagesaurids (Godoy et al. 2014). Indeed, their ziphodont teeth (unicuspidated, laterally compressed with serrated carinae) associated with the biomechanical performances of their skull allowed sebecosuchians to effectively handle prey after wounding it (Montefeltro et al. 2020). It is noteworthy that the inferred ectothermic sebecosuchians occupy a niche usually occupied by endothermic theropod dinosaurs (Benson et al. 2013; Zanno and Makovicky 2013).

The ectothermic condition of Notosuchia suggested by Cubo et al. (2020) is supported by the results of the present study using larger sample sizes of extant species and a more robust phylogenetic comparative method (PLR). First, the finding that RBC_{size} is related to bone vascular canal diameter (Huttenlocker and Farmer 2017) is supported by our results obtained using a sample size more than twice of that used by these authors. Thus bone vascular canal diameter can be used as a proxy to infer RBC_{size}, and then endothermy (because within tetrapods, RBC_{size} is lower in endotherms than in ectotherms; Snyder and Sheafor 1999). Huttenlocker and Farmer (2017) included bone cross-sectional area (in addition to bone vascular canal diameter) as an explanatory variable in models aimed at explaining the variation of RBC_{size}. However, in their study, bone cross-sectional area did not improve the explanatory power of models and was not retained (Huttenlocker and Farmer 2017; Supplemental Information, “Analysis I, Training Data Set for Extant Taxa”). Unexpectedly, our analyses (using a larger sample size) showed that bone cross-sectional area significantly improves the explanatory power of models and is retained, together with bone vascular canal diameter, in models explaining the variation of RBC_{size}. The fact that the estimate for bone cross-sectional area is always negative (Table 1) indicates that RBC_{size} decreases as bone cross-sectional area increases. Second, notosuchian thermometabolism is inferred here using a larger sample of extant tetrapods (more than three times the sample used by Cubo et al. [2020]) and logistic phylogenetic regressions (a method more reliable than those used in previous studies; see below). We are aware of the fact that Huttenlocker and Farmer (2017) showed that histological changes reflect changes in VO_{2max} better than changes in thermometabolism. However, we have found here that microstructural variation is linked to thermometabolism too. The models constructed to infer the probability of endothermy using vascular canal size as an explanatory (predictive) variable were highly significant, and the classification errors obtained are quite low (6.5% using Can_{min}). Thus we conclude that we can use these models confidently in paleobiological inference of thermometabolism.

Thermal paleophysiology is an emergent discipline (Cubo and Huttenlocker 2020). It has great potential resulting from the synergy between physiological studies aimed at deciphering the mechanisms of thermogenesis in extant amniotes (e.g., Bal and Periasamy 2020; Jastroch and Seebacher 2020; Grigg et al. 2021) and paleobiological inferences in extinct amniotes using phylogenetic comparative methods (e.g., Cubo et al. 2012, 2020, 2022; Legendre et al. 2016; Huttenlocker and Farmer 2017; Olivier et al. 2017; Fleischle et al. 2018; Cubo and Jalil 2019; Faure-Brac et al. 2021; Knaus et al. 2021). Logistic phylogenetic regressions are the third step in efforts performed during the last decade to carry out reliable inferences of thermometabolic status in extinct amniotes. A decade ago, Cubo et al. (2012) inferred bone growth rates using bone histological features and multiple linear regressions tested for significance using permutations in order to circumvent the nonindependence of the observations due to the phylogeny. Considering that bone growth rates are significantly related to RMR in amniotes (Montes et al. 2007), the former was used by Cubo et al. (2012) as a proxy to infer the thermometabolic status of extinct archosaurs. This method was used by Legendre et al. (2013) to infer the bone growth rate and the thermometabolic status of *Euparkeria*. In a second step, Legendre et al. (2016) adapted Guénard et al.'s (2013) PEMs to perform paleobiological inferences of RMR. This contribution represented significant methodological progress, because paleobiological inference models included the phylogeny (rather than circumventing its effects, as did the preceding method), assuming an evolutionary model (Molina-Venegas et al. 2018). PEMs have been widely used to infer the thermometabolic status of extinct amniotes (Legendre et al. 2016; Olivier et al. 2017; Fleischle et al. 2018; Cubo and Jalil 2019; Cubo et al. 2020, 2022; Faure-Brac and Cubo 2020; Faure-Brac et al. 2021; Knaus et al. 2021). Using logistic phylogenetic regressions is a new step in this sequence. This method improves upon the previous approach by using all of the phylogenetic information (rather than a fraction of it, as did PEMs in order to avoid model overfitting). An encouraging sign is that results are congruent in spite of the diversity of methods used to obtain them. Inferring the maximum metabolic rate of Notosuchia using the size of femoral nutrient foramina (Seymour et al. 2012) would be the next promising step to fully understand the thermophysiology of these amazing crocodylomorphs.

## Acknowledgments

We thank H. Lamrous (Sorbonne Université) for helping us with picture acquisition using bone thin sections at the Vertebrate Hard Tissue Collection of the Museum national d'Histoire naturelle (Paris). This study was partly funded by the project Emergences Sorbonne Université 2019 no. 243374 to J.C. The authors declare no competing interests.

## Data Availability Statement

Data available from the Dryad and Zenodo Digital Repositories: https://doi.org/10.5061/dryad.80gb5mktb, https://doi.org/10.5281/zenodo.6795231.