Where nutrient foramina were not obscured by matrix, digital photographs were acquired with a Canon Ixus 220 HS digital camera with a scale bar or calipers placed in the frame for scale. Where nutrient foramina could be located but were still occluded, matrix was carefully removed from the canal with an airscribe and needles until the cross-sectional area ceased to decrease below the superficial funnel. Care was taken to avoid widening the nutrient canal in this process. Then, photomicrographs were taken using either a camera adapter on a stereo microscope with a digital camera or a Hirox KH 7700 digital microscope at different planes of focus. In using photography, the picture was taken from an angle that provided visibility of the deepest unobscured point of the nutrient canal.

When nutrient foramina could not be cleared sufficiently or were not visible, femora were scanned using a General Electric phoenix V|tome|x S180/240 CT scanner (Section Palaeontology, Institut für Geowissenschaften, Universität Bonn) or a General Electric phoenix V|tome|x S240 CT scanner with custom added 180 kV tube (PaleoCT lab at the Department of Organismal Biology and Anatomy, University of Chicago), using either the 180 or 240 kV tube varying according to specimen dimensions.

Based on the μCT scans, image stacks were produced with the phoenix datos|x software and volumetric models were constructed from these image stacks using VG Studio Max software by Volume Graphics GmbH. To produce virtual cross sections of the nutrient canals, the canal paths were identified in the volumetric models. For this, we first marked the nutrient foramen in a slice window. Secondly, the canal was traced to its full extent and connected to the first mark through a vector. During this step, we also controlled for the direction of the nutrient canal toward the ossification center to avoid measuring vascular canals that did not directly feed into the medullary cavity (compare Figures 1, 2). Third, the vector was aligned coaxially with the view angle and fourth, scrolling through the cross sections along the canal, we identified the point of minimal canal cross sectional area. Last, we obtained an image file with scale bar from the virtual cross section in isometric view using the snapshot tool. In three specimens, both photographic and μCT methods were used to produce image files of the nutrient canal cross sections for cross-validation of methods and resulted in identical values. Hu et al. (2020) found no significant difference between radius and cross-sectional area of nutrient canals measured from microphotographs and micro-CT in a larger sample.

Image files from 3D models and digital photographs were loaded into ImageJ 1.49k (Fiji package) image analysis software, and cross-sectional area of the nutrient canal was averaged from at least three measurements of each file to minimize tracing subjectivity. If multiple nutrient canals had been identified, their cross-sectional areas were added. Nutrient foramen radius was calculated from a perfect circle with the same cross-sectional area as the nutrient canal as needed for Seymour’s index of blood flow.

Seymour et al. (2012) derived the index of blood flow (Q_i) from the Hagen-Poiseulle equation, which gives the flow rate of blood (Q) given laminar flow through blood vessels depending on blood viscosity, blood pressure difference (P), vessel radius (r), and vessel length (L): Q = (Pπ r⁴)/(8 L η). As blood flow and viscosity are independent of body size, both terms are eliminated to get the index of blood flow equation: Q_i = r⁴/L.

Body mass (Mb) for non-mammalian synapsids was calculated with the equation of Campione and Evans (2012) from mid-diaphysis circumference of the humerus (CH) and the Femur (CF): M_b = 10^[2.749 * log₁₀ (CH + CF) − 1.104)].

In several cases, however, only isolated femora were available. For Dimetrodon, the CH/CF ratios of Shelton et al. (2013) were employed. In the case of therapsids, ratios were measured in closely related complete skeletons at GPIT (Tübingen Palaeontological Collection, Tübingen, Germany; Sauroctonus for gorgonopsids, Chiniquodon for cynodonts, and Tetragonias for dicynodonts). For other non-therapsid synapsids, the CH/CF ratio was calculated from values of closely related taxa published in Shelton (2015). The CH value was then calculated from the CH/CF quotient and used in mass estimation.

For visualization purposes, Q_i was plotted against body mass on a double logarithmic scale using Microsoft Excel. We performed simple linear least squares regressions for non-mammalian synapsids, non-avian dinosaurs, birds, mammals, varanids, and non-varanid ectotherms, adding calculated 95% confidence intervals to maintain comparability with the original results of Seymour et al. (2012) (see Figure 3).

First, a phylogenetic tree was constructed from the recent taxa of Seymour et al. (2012) and Allan et al. (2014) using the Timetree website¹ (Hedges et al., 2006, 2015; Kumar and Hedges, 2011) and saved in Nexus format (Newick notation). The tree file was loaded into Mesquite software (version 3.06; Maddison and Maddison, 2019) equipped with ‘Stratigraphic tools for Mesquite’ (version 1.0c, Josse et al., 2006) and manually modified to add extinct taxa. As trees had to be fully resolved for analysis with ‘nlme’ (Pinheiro et al., 2014) in R (R Core Team, 2021), polytomies were resolved to multiple dichotomies with a branch length of 0.0000001 (equaling 0.1 years). Phylogenetic relationships and divergence dates of extinct taxa were taken from Benson et al. (2014) for dinosaurs, Benson (2012) for basal synapsids, from Kemp (2012) for the interrelationships of major therapsid groups, from Fröbisch (2007) for dicynodonts, and from Botha-Brink et al. (2012) for cynodonts. Following Field et al. (2014), turtles were treated as the sister group to Archosauromorpha with the divergence dated according to the appearance of the oldest stem turtle Eunotosaurus (Bever et al., 2015, but for a contrary opinion see Lichtig and Lucas, 2021).

Statistical comparative methods generally assume independence, homoscedasticity, and normality of the data points analyzed to produce reliable results that allow reliable deductions. This is required because statistical methods test a priori for null-hypotheses that assume no relationship between variables and our regressions require normality. In data sets derived from organisms, however, said requirements may be violated due to the phylogenetic relationships between taxa (Felsenstein, 1985). Therefore, we performed phylogenetic generalized least squares regression (PGLS) (Grafen, 1989), following the methodology of Mitchell et al. (2017), to detect potential phylogenetic signal (Blomberg and Garland, 2002) in the relationships of Q_i and body mass. Additionally, PGLS confirmed significant relationships between Q_i and body mass for previous data under phylogenetic compensation (Seymour et al., 2012; Allan et al., 2014).

Q_i and body mass data from Seymour et al. (2012) as well as phylogenetic trees (Newick format) were fed into R and tested for the best fit of evolutionary model using the ‘fitcontinuous’ function of the R package ‘geiger’ (Harmon et al., 2008). A choice of Brownian Motion (Felsenstein, 1973), Ornstein-Uhlenbeck (Butler and King, 2004), Early Burst (Harmon et al., 2010), and White Noise models was available. ΔAICc (Burnham et al., 2011) gave the best-fitting model the lowest value, which was then selected.

When ‘fitcontinuous’ recommended the White Noise model, compensation for phylogenetic signal was not necessary, but in all other cases, the function ‘pgls’ of the package ‘caper’ (Orme et al., 2013) or the function ‘gls’ of the package ‘nlme’ was used with phylogenetic weights and the phylogenetic correlation structure “corMartins” of the package ‘ape’ (Paradis et al., 2004). To ensure normality in the residuals after phylogenetic adjustment, we used the Lillefors test (‘lille.test’ function) of the R package ‘nortest’ (Gross and Ligges, 2015). When the test showed p-values below 0.05, quantile-quantile plots were consulted to check for the normality of distribution and fitted values plotted against residuals to control for heteroscedasticity (see Supplementary Figures 1, 2). Outliers were then removed to produce a normal distribution of residuals and homoscedasticity in the dataset.

To confirm the distinctness of varanids from other ectothermic amniotes for Q_i vs. body mass regression under phylogenetic compensation, we performed phylogenetic ANCOVA (analysis of covariance) on the subset of recent species from Seymour et al. (2012), applying the methodology of Smaers and Rohlf (2016). We used the function ‘gls.ancova’ of the ‘evomap’ package (Smaers and Mongle, 2014) of R (R Core Team, 2021), setting Q_i as the dependent variable, higher taxon as indicator variable, and body mass as covariate (Table 2).

To compare distinctness of regression slope and intercept for the Q_i vs. body mass plots of our non-mammalian synapsids to previously published recent datasets (Seymour et al., 2012; Allan et al., 2014), we performed pairwise ANCOVA (analysis of covariance) with the ‘aov’ function of the ‘stats’ package of R, designating Q_i as the dependent variable, higher taxon as factor, and body mass as the covariate. We applied a Bonferroni correction to counteract the problem of multiple comparisons, lowering the desired significance level to α = 0.05/9 = 0.0056. If F and p-values indicated no significant slope difference, we used the function ‘anova’ of ‘nlme’ package to test for significant differences in elevation (Table 3).

We calculated phylogenetic eigenvector maps from the phylogeny using the R-package “MPSEM” (Guénard et al., 2013). For this we used the reduced dataset of 19 recent taxa with MMR values (see section “Construction of Phylogenetic Tree”) in combination with our extinct synapsid taxa. AICs corrected for small sample sizes (Hurvich and Tsai, 1989) was used to decide between the predictive models, with and without the co-predictor Q_i. Cross-validation was provided by leave-one-out cross-validation. Selecting the co-predictor Q_i (AICc 0.9407704, R² 35.52039; without co-predictor AICc 0.87, R² 42.272), we used the predictive model to estimate MMR for the extinct synapsid taxa and 95% confidence interval boundaries. To test for the influence of Chelonia mydas (with exceptionally high MMR for a turtle) on the reconstructed values, we reran the analysis without this taxon.

Phylogenetic Eigenvector Maps (PEM) belongs to a family of phylogenetic comparative methods aimed at inferring missing values in samples where values are known for the great majority of species (Molina-Venegas et al., 2018). Because of the large number of predicted values (31) relative to the smaller number of known values from which the prediction model was calculated (19), we performed a series of PEM analyses using subsets of extinct taxa and the full extant data set. Values from all subset analyses were identical to the predictions for the entire dataset.

Mass-independent MMR values are given in mL O₂ h^–1 g^–0.⁶⁷. Although maximum aerobic metabolism is empirically shown to correlate with body mass to the power of 0.82 in ectothermic amniotes (Bennett and Dawson, 1976) or 0.86 in mammals (White and Seymour, 2005), different from 0.76 for BMR (Withers, 1992), many metabolic processes in animals are dependent on surface area of body membranes, which correlate with body mass to the power of 0.67 (Hemmingsen, 1960). For this reason and to provide better comparability with previous studies of metabolic rates in synapsids, we chose the lower exponent.

The sample of non-mammalian synapsids spans two orders of body size magnitude, with the lightest animal estimated at 1.86 kg (Galesaurus) and the heaviest at 503 kg (Dinodontosaurus). Q_i-values of non-mammalian synapsids reach from a minimum 2.653E-5 mm³ for the second-smallest specimen of Clepsydrops to 3.918E-3 mm³ for the longest non-therapsid synapsid femur of Dimetrodon grandis. These two specimens illustrate that Q_i correlates strongly with body mass in our synapsid dataset as previous studies have found in other groups (Seymour et al., 2012; Allan et al., 2014).

In plotting the index of blood flow against body mass (Figure 3), non-therapsid synapsids and non-mammalian therapsids are statistically indistinguishable from recent endotherms (mammals and birds), but also from varanids. Varanids have been shown (Seymour et al., 2012) to be statistically different from mammals in slope, but not in intercept. Like recent endotherms and varanids, stem mammals are significantly different from recent non-varanid ectotherms. While varanids differ in slope from other non-avian sauropsids (Seymour et al., 2012), endotherms and non-mammalian synapsids differ significantly in intercept from non-varanid ectotherms. Like certain bird and mammal species, two therapsids and one non-therapsid synapsid overlap with a non-varanid sauropsids in the range of values. These do not impact the statistical significance of the results.

Comparing regression functions with those for non-varanid ectotherms, average Q_i-values are 13.6-fold higher for therapsids (Figure 4) and 12.9 times higher for non-therapsid synapsids (Figure 5) along the common body mass range (5.5 to 456.4 kg), with an overall stem-mammal value of 13.1 times the Q_i of sauropsids in stem mammals. These blood flow values compare closely with those of mammals (13.4 times average non-varanid ectotherm Q_i over the body mass range).

The search for the best-fitting evolutionary model found phylogenetic signal in all groups except varanids and non-mammalian synapsids, which were best described with the White Noise model. Among the groups that showed a phylogenetic signal, all were best described by the Ornstein-Uhlenbeck model of evolution and showed homoscedasticity and normality in the distribution of the residuals. The PGLS analysis of the entire dataset resulted in an Ornstein–Uhlenbeck model of evolution being selected, showing phylogenetic signal in the data. Yet, after the analysis was performed and gave a highly significant result (p < 0.001), a Lillefors test for normality returned a statistically significant result (indicating deviation from the normality in the residuals), and our Q–Q plot revealed numerous outliers among the dinosaur dataset. After the dinosaur data set (Seymour et al., 2012) and the further outlier Emeus crassus were removed, the Lillefors test passed in the phylogenetically corrected residuals.

The PGLS analysis returned a model with a much steeper slope than that of the individual groups of 0.9675. The Q–Q plot revealed a degree of heteroscedasticity caused by a low number of data points on the extreme ends of the distribution (see Supplementary Figure 1). The distribution is likely caused by the low number of sauropsids of high body mass (above 10 kg) and low body mass (below 100 g) in the data set.

Phylogenetic ANCOVA analysis returned a significant result for the distinction between varanids and non-varanid ectothermic amniotes in Q_i vs. body mass regressions (Table 1). Pairwise ANCOVA returned significant results for the distinction of varanids from other ectothermic amniotes, but we found no significant difference between varanids and endothermic groups.

Mass-specific MMR ranges over two orders of magnitude (4 to 637.7 mL O₂ h^–1 g^–0.67 in the common blue-tongued skink (Tiliqua scincoides) and the domestic goat (Capra hircus), respectively) (see Figure 6). Two Tiliqua species have the lowest values. Varanid values are higher, but the green sea turtle Chelonia mydas exceeds them in MMR. Tachyglossus, the echidna, has a lower MMR than Varanus gouldii and C. mydas. Therian mammals have the highest values, at least an order of magnitude above Tiliqua. The MMR values predicted by our analysis for non-mammalian synapsids have a much lower variation (17.5 to 37.3 mL O₂ h^–1 g^–0.67 in Lystrosaurus sp. and an undetermined dicynodont, respectively) than the complete dataset. Predicted values overlap with varanids, Tachyglossus, and C. mydas. Differences between the range of values of non-mammalian therapsids and non-therapsid synapsids are not visible.

Excluding C. mydas from the dataset resulted in lower reconstructed values in extinct synapsids. The lowest value among non-mammalian synapsids, 18, decreases to 14.2 mL O₂ h^–1 g^–0.67 in Lystrosaurus, lower than Varanus panoptes and Varanus varius, but higher than Varanus mertensi.

Evolutionary model tests detected no phylogenetic signal in either the non-therapsid synapsid or non-mammalian therapsid subsets. This pattern coincides with the wide range of Q_i values in both paraphyletic groups on the Q_i-M_b plot. The combination could represent a diverse pattern of mass-specific Q_i in the samples due to either widely differing blood flow in taxa, differing loading regimes resulting from fundamentally differing locomotion styles such as in birds (Allan et al., 2014), or the pattern could be an artifact of sampling. As blood flow in extinct taxa cannot be directly measured, the first hypothesis cannot be rejected but would require the ad hoc assumption of stem-mammal exceptionality. Therefore, we will pursue the other alternatives further. Differing femur loading regimes in major taxa influencing Q_i have been postulated for birds, but usually result in overall homogeneous values for closely related groups with similar morphology [e.g., gruiform birds in the dataset of Allan et al. (2014)]. In our non-mammalian synapsid data, however, closely related taxa can show widely differing Q_i values such as different specimens of Dimetrodon that can vary by an order of magnitude. Such wide differences in individuals from closely related populations could result from biological variability but also from sampling or measurement errors.

Juvenile specimens could present outliers. The fossil dinosaurs analyzed by Seymour et al. (2012), for example, include two outliers in their distribution whose values we suspect resulted from ontogenetic effects. The Ornithomimidae indet. and Thescelosauridae indet. femora are particularly short when compared to closely related species and could thus belong to juveniles. This might explain their high Q_i values and thus large nutrient foramina due to the increased need for blood flow in a growing bone. By selecting only the largest fossil femora for a given taxon and excluding many incompletely ossified bones, we minimized this possible source of sampling error. Yet, ontogenetic influence cannot be completely eliminated except in specimens with ossification states indicating late ontogeny such as the edaphosaurs, which have predicted MMR values higher than Tachyglossus.

The preservation of most sampled bones from the Carboniferous to the Triassic presents an alternative possible source of measurement error. Fractures of bones can obliterate nutrient canals. Additionally, long bones tend to break at the point of least resistance. In femora this occurs in particular at the mid-diaphysis, the most common location of nutrient canals (compare Figure 1). For this reason, nutrient canals could not be detected in numerous fossil synapsid femora that we screened in preparation for our study. Several femora preserved multiple nutrient canals of differing size (compare Figure 2). If the largest had been destroyed in a specimen, but at least one smaller canal were present and measured, the resulting value would be much smaller than would be measured on the intact femur. Thus, we cannot exclude the possibility that the Q_i values recorded for fossil taxa are based on a fraction of the original nutrient canal cross-sectional area. Therefore, we suspect that the non-mammalian synapsid average Q_i values observed represent conservative estimates due to the caveats of data collection on fossils.

Q_i elevated over non-varanid ectothermic amniotes has been postulated for Morganucodon (Newham et al., 2020), agreeing with our findings in elevation, but since our lightest specimens have body masses at least an order of magnitude higher than the mammaliaform, comparisons with our data are difficult.

Traditional least squares regressions as well as phylogenetic and pairwise ANCOVA alone indicate that non-mammalian synapsids potentially exhibited maximum aerobic metabolic rates elevated over the level of recent ectotherms other than varanids. This hypothesis is substantiated by the reconstructed MMR values for non-mammalian synapsids using phylogenetic eigenvector maps. However, the subset of taxa with MMR measurements is limited to low sampling size, which might skew the results. The relatively high aerobic capacity of the only turtle in the dataset for which MMR values are available (Chelonia mydas) and the high abundance of varanids in the dataset used for PEM analysis skews the ancestral state reconstruction toward higher values.

The MMR value of C. mydas is higher than for other sauropsids, likely an adaptation for long-endurance swimming (Penick et al., 1996). MMR values with and without Chelonia differ, but excluding the turtle from the analysis only results in more species of varanids overlapping with non-mammalian synapsids (see Supplementary Material). The inclusion of other turtles with a lower MMR might decrease the ancestral MMR value for sauropsids but would likely not impact our main result. Additionally, endothermic archosaurs are equally excluded from the dataset, the inclusion of which might raise the ancestral value for amniotes. Overall, more sampling is needed in order to obtain more robust ancestral value reconstructions for early amniotes.

The index of blood flow serves as an instrument to estimate blood flow in the absence of data on blood pressure. In fact, differences of an order of magnitude between the Q_i of endotherms and non-varanid, non-avian sauropsids are translated into nearly three times higher blood pressure and a higher oxygen carrying capacity of mammal blood to produce an about fiftyfold increase in oxygenation rate of the mammalian femora over those of sauropsids other than birds and varanids (Seymour et al., 2012). No recent ectothermic amniote maintains a mean arterial blood pressure (MAP) nearing 15 kPa. Ectothermic tetrapods have incompletely separated pulmonary and systemic circulations with low blood pressures in both. In Varanus niloticus, the derived morphology of the varanid heart facilitates the highest MAP of ectothermic tetrapods at 11.6 kPa, in the lower range typical of mammals (Millard and Johansen, 1974). This indicates femoral oxygenation rates for varanids intermediate between mammals and other non-avian sauropsids, because the oxygen carrying capacity of varanid blood is no higher than that of other non-avian sauropsids (Bennett, 1973). Hence, blood pressure values are required to determine whether the blood flow rates of non-mammalian synapsids can be compared with that of recent varanids.

Seymour (2016) applied physiological principles in recent animals to determine blood pressure in extinct dinosaurs and Dimetrodon. As a continuous blood column must transfer pressure from the location of the heart to the highest point of the animal and since the density of blood is near constant (about 1000 kg/m³) across taxa, we can calculate the pressure at height of the heart necessary to maintain the blood column to the highest point of the body (P_g). In sail-backed non-therapsid synapsids such as Dimetrodon grandis and Edaphosaurus pogonias, distances between the location of the heart at the ventral side of the body cavity to the highest neural spines measure 1.5 m and 1 m, respectively (Romer and Price, 1940). In addition to this pressure, perfusion pressure is required at the highest point of the body to calculate the actual mean arterial pressure (MAP), which is unknown. Thus, mean arterial blood pressures were considerably higher than the MAP of non-avian, non-varanid sauropsids (about 5 kPa) in the two synapsids because P_g alone for Edaphosaurus (9.7 kPa) lies in the varanid range of MAP or even well within the mammalian range of MAP for Dimetrodon (15.3 kPa). As non-therapsid synapsid femora articulate in a horizontal to more parasagittal orientation at the height of the heart or below, blood pressure at the femoral nutrient foramen was as least as high as mean arterial blood pressure.

Considering the Q_i values comparable to varanids and endotherms, we conclude that the high calculated blood pressures facilitated blood flow rates equal or higher than varanids in some non-therapsid synapsids.

Huttenlocker and Farmer (2017) inferred blood cell sizes for non-eutheriodont synapsids to be larger than those of endotherms, in the range of modern non-avian sauropsids. Inversely correlating this value with relative heart size as an indicator for stamina, they suggested that elevated stamina evolved no earlier than in ancestral therocephalians and cynodonts. Yet, varanids reach much higher AAS and MMR than other ectothermic sauropsids while maintaining large blood cell sizes and low relative heart sizes (Vinogradov and Anatskaya, 2006).

However, our predictive model of MMR does not recover elevated activity metabolism in cynodonts compared to non-cynodont synapsids. Although red blood cell (RBC) sizes likely play a role in the high stamina of mammals and birds (Snyder and Sheafor, 1999; Vinogradov and Anatskaya, 2006), our prediction model for MMR does not account for RBC size and is only sensitive to the distinction between very low (i.e., turtle, lepidosaur, and crocodile) levels and intermediate to high (i.e., varanid and mammalian) levels. Therefore, our model likely underestimates MMR of cynodonts with lower RBC sizes. We hypothesize that ancestral synapsids showed elevated stamina such as in varanids and that MMR increased further in therocephalians and cynodonts. At least two phases of increase imply gradual evolution of activity metabolism in synapsids.

Higher MAP in the systemic circulation needs to be separated from the pulmonary circulation, because high blood pressures in the latter can cause fatal pulmonary edema (Smits, 1989). Varanids and pythons are the only non-archosaurian sauropsids capable of considerable pressure separation between the two circulations through specially adapted hearts (Wood et al., 1977; Burggren and Johansen, 1982; Wang et al., 2002). Recent archosaurs and mammals feature four-chambered hearts for separation of the pulmonary from the systemic blood circulation and thus have much higher MAP. Mammalian-level blood pressures in Dimetrodon imply that the four-chambered heart or a functional analog allowing systemic pressure separation was present in the earliest Sphenacodontoidea.

Experiments with Alligator mississippiensis suggest that high blood pressures might be the direct consequence of a four-chambered heart: Normally, crocodiles show MAP in the amniote ectotherm range (e.g., Jones and Shelton, 1993). When Eme et al. (2009) surgically removed the pulmonary systemic shunt common to all crocodilians in their test animals, blood pressures increased markedly. We suggest that once the four-chambered heart or a functional analog was acquired, subsequent blood pressures remained high in all Sphenacodontoidea, including therapsids.

The hypothesis of elevated metabolism in non-mammalian synapsids helps to explain several patterns in the fossil record inconsistent with previous hypotheses of a “reptilian-like” physiology in non-mammalian synapsids. Analyses of non-therapsid synapsid bone histology have uncovered FLB, a primary bone tissue typical of recent endotherms, in the early-branching Ophiacodon (Shelton and Sander, 2017) and incipient FLB in Dimetrodon (Shelton et al., 2013). This bone tissue is indicative of high growth rates, and only the secondarily ectothermic crocodiles (Seymour et al., 2004) can produce incipient FLB among recent non-avian sauropsids (Woodward et al., 2014). Elevated MMR in synapsids thus likely coincided with the elevated BMR needed for higher growth rates and FLB tissue in the earliest synapsids.

Non-therapsid synapsids also present morphological characters that we interpret here as indicative of a more energetic lifestyle. The evidence for a diaphragm-homologue in caseids suggests that the early synapsids could maintain breathing during locomotion (Lambertz et al., 2016). Plesiomorphically, continuous breathing necessary for elevated metabolic rates is hindered by locomotion in tetrapods as the axial body wall musculature tasked with undulatory locomotion cannot simultaneously inflate the lungs, known as Carrier’s constraint (Carrier, 1987). Just as endotherms, two highly active squamates groups (varanids and tegus) have evolved mechanisms to overcome this constraint (Owerkowicz et al., 1999; Klein et al., 2003).

A strong component of lateral undulation in locomotion is typical of non-avian sauropsids and associated with the sprawling stance typical of classical ectotherms. Mammals and archosaurs, on the other hand, evolved an erect stance, thus increasing stamina through the capacity of breathing during exercise with the help of the body wall musculature (Carrier, 1987). Traditionally, non-therapsid synapsids have been considered representing the basal amniote condition in the morphology of their locomotory apparatus (e.g., Blob, 2001), but non-therapsid synapsids in general, and Sphenacodontia in particular, show marked adaptations for a reduction of lateral undulation. Jones et al. (2021) found distinct locomotory adaptations in the vertebrae of non-mammalian synapsids such as increased stiffness of the column that might represent an early exaptation for the parasagittal gait that later evolved in mammals. Basal synapsid zygapophyses are oriented markedly different from the horizontal plane (30° in ophiacodonts, 45° in Dimetrodon), thus restricting the lateral movement of the vertebral column (Kemp, 2005: 101). Trackways of the Lower Permian ichnotaxon Dimetropus leisnerianus produced by (possibly sphenacodontid) non-therapsid synapsid trackmakers show an average pace angulation of about 110°, higher than the average angulation of the contemporary tetrapods at about 90° indicating a more erect stance [see Supplementary material of Kubo and Benton (2009)]. The trend to the parasagittal gait necessary for higher locomotion capacity [as a precondition to mammalian endothermy, compare Lovegrove (2017)] consequently began in non-therapsid synapsids.

To determine what elevated aerobic metabolism in the earliest synapsids means for the chronology of mammalian endothermy evolution, we must assess the evidence for BMR in the mammalian stem group. Rey et al. (2017) produced geochemical evidence from three different assemblage zones in the Beaufort Group of South Africa for homeothermy in dicynodonts and cynodonts. Olivier et al. (2017) and Faure-Brac and Cubo (2020) used proxies for BMR (or the similar resting metabolic rate, RMR) from bone histology to infer mammalian-like BMR/RMR in dicynodonts. Therefore, metabolic thermogenesis evolved at the base of Neotherapsida (the last common ancestor of dicynodonts and mammals and all its descendants) no later than 270 Mya (Liu et al., 2010). de Ricqlès (1974) based his hypothesis of endothermy in all therapsids explicitly on the presence of FLB, but incipient FLB is also reported from alligators (Woodward et al., 2014). The only indicators of elevated basal metabolism in non-therapsid synapsids remain high blood pressures (Seymour, 2016) and the presence of FLB as well as incipient FLB (Shelton et al., 2013; Shelton and Sander, 2017). Their histology suggests that elevated growth rates and high blood pressures facilitated higher metabolic rates, but these facts by themselves do not justify pushing back the origin of endothermy into non-mammalian synapsids. Consequently, we must reject the thermoregulation-first hypothesis. Elevated aerobic capacity was present first in non-therapsid synapsids, while endothermy evolved only second in neotherapsids or later.

Non-therapsid synapsids had high MMR (this study) and low RMR/BMR (Faure-Brac and Cubo, 2020), whereas non-mammalian neotherapsids had high MMR (this study) and high RMR/BMR (Faure-Brac and Cubo, 2020). More precisely, Faure-Brac and Cubo (2020) showed that all non-therapsid synapsids analyzed but one (Ophiacodon uniformis) had RMR/BMR significantly lower than the threshold value separating endotherms from ectotherms. The value inferred for Ophiacodon uniformis was lower than this threshold, but the upper limit of the 95% confidence interval was slightly higher than it. As a result, we amend the foraging hypothesis of Hopson (2012) and the triphasic model of endothermy evolution of Lovegrove (2017). Our finding of high blood flow rates and inferred elevated aerobic capacity in synapsids predating the evolution of physiological thermoregulation is generally in agreement with the two models, although they suggest that elevated MMR evolved only in therapsids. Hopson’s idea of a varanid-like animal with increased aerobic scopes due to elevated MMR and low BMR as an intermediary stage to the first therapsid endotherms is a hypothesis our data supports. Judging from additional indicators of increased locomotor capacity (an exaptation for Lovegrove’s second phase of endothermy evolution) and inferred soft part anatomy capable of coping with high blood pressures as well as increased oxygen requirements, our picture of the non-therapsid synapsids shifts: away from that of passive sit-and-wait predators that rely mostly on anaerobic metabolism toward widely foraging animals with elevated aerobic metabolic rates for sustained activity. Coincidentally, the varanid analogy in metabolism used by Hopson (2012) confirms the intuitive connection made by Broili (1904) and Williston (1911) in the descriptions of the non-therapsid synapsids Varanosaurus (and Varanops) due to their morphological resemblance to varanids.

Our findings could also be interpreted to support the assimilation capacity model of Koteja (2000) since a widely foraging lifestyle could have evolved to satisfy metabolic requirements of intense parental care in the parent animal. Botha-Brink and Modesto (2007) published tentative evidence of parental care in varanopsids. However, the only conclusive evidence for the assimilation capacity model would be measurements of juvenile growth rates in non-mammalian synapsids that are significantly higher than in recent non-avian sauropsids. This could be obtained from bone histology.

Our tentative evidence for elevated MMR in the earliest synapsids is supported by the fact that ancestral amniotes most parsimoniously had complex lungs and that the simple lungs of squamates originated in an ancestral lepidosauromorph as an adaptation to extremely small body size (Lambertz et al., 2015). Also, endothermy appears to have evolved several times independently in diapsids, at least once among archosaurs (Legendre et al., 2016), and possibly three times in marine reptiles (Bernard et al., 2010; Fleischle et al., 2018). A low metabolism, ectothermic lifestyle has diverse advantages (Pough, 1980). We speculate that the taxa at the two extreme ends of the amniote metabolic spectrum, energetically efficient ectothermic non-avian sauropsids on the one hand, high throughput endothermic mammals and birds on the other, could have diverged from an ancestral state of low BMR and elevated aerobic metabolism.