Inner ear biomechanics reveals a Late Triassic origin for mammalian endothermy

Endothermy underpins the ecological dominance of mammals and birds in diverse environmental settings1,2. However, it is unclear when this crucial feature emerged during mammalian evolutionary history, as most of the fossil evidence is ambiguous3–17. Here we show that this key evolutionary transition can be investigated using the morphology of the endolymph-filled semicircular ducts of the inner ear, which monitor head rotations and are essential for motor coordination, navigation and spatial awareness18–22. Increased body temperatures during the ectotherm–endotherm transition of mammal ancestors would decrease endolymph viscosity, negatively affecting semicircular duct biomechanics23,24, while simultaneously increasing behavioural activity25,26 probably required improved performance27. Morphological changes to the membranous ducts and enclosing bony canals would have been necessary to maintain optimal functionality during this transition. To track these morphofunctional changes in 56 extinct synapsid species, we developed the thermo-motility index, a proxy based on bony canal morphology. The results suggest that endothermy evolved abruptly during the Late Triassic period in Mammaliamorpha, correlated with a sharp increase in body temperature (5–9 °C) and an expansion of aerobic and anaerobic capacities. Contrary to previous suggestions3–14, all stem mammaliamorphs were most probably ectotherms. Endothermy, as a crucial physiological characteristic, joins other distinctive mammalian features that arose during this period of climatic instability28. The functional morphology of the fluid-filled semicircular ducts of the inner ear is adapted to body temperature and behavioural activity and can be used to investigate the evolution of endothermy.

help to solve this problem. Semicircular ducts monitor head rotation and are filled with endolymph, a fluid whose viscosity affects function and depends on T b . Perception of head motion is essential for navigation 22 , balance 21 , vision 20 , motor coordination 21 and spatial awareness 21 . Therefore, SDS function is expected to attune to the spectrum of head rotations (that is, angular velocities and ordinary frequencies) experienced by an organism 27 . In this context, endotherms not only need to compensate for the decreased response speed (that is, decreased upper corner frequency) induced by the lower endolymph viscosity 23,24 associated with elevated T b 23,24 , but may also require a more efficient SDS than ectotherms (Supplementary Methods) because they are more active on average 1,2,25,26,33,34 . Optimal SDS function can be achieved by modifying endolymph chemistry to increase its viscosity and/or altering SDS morphology.
To track these adaptations in the fossil record, we developed the thermo-motility index (TMI), a biomechanically informed proxy derived primarily from functionally relevant morphometric parameters of the SDS and corrected for allometry (Supplementary Methods). The TMI is defined so that for a given body size, relevant SDS function = TMI + a temperature term (inversely related to T b ). T b does not contribute to the calculation of the TMI. Nevertheless, measuring the TMI should enable prediction of T b , because assuming that SDS function is attuned to natural head motion, an increase in T b (reducing the temperature term) should be proportionally compensated by an increase in the TMI. Similarly, because the TMI is positively related to maximum angular head motion (Methods and Supplementary Methods), when behavioural activity increases, the TMI may have to increase to improve SDS function. As the SDS does not fossilize, we assembled a new dataset of 50 membranous labyrinths of various vertebrates to establish morphofunctional correlations between semicircular ducts and bony canals (Supplementary Data 2 and Supplementary Note 2). Subsequently, we computed the TMI of 362 specimens, including 68 extinct synapsids, from bony labyrinths visualized through microtomography (Fig. 1, Extended Data Fig. 1 and 2 and Supplementary Data 2 and 3).

Empirical testing of TMI accuracy
We demonstrated empirically that the TMI is explained by and positively correlated to T b (Fig. 2, solid curve; phylogenetic generalized least-squares regression (PGLS), relative likelihood against null model = 1.00; Supplementary Note 2). Among the semicircular canals, the TMI of the anterior canal, which shows the highest correlation with T b , is better explained by T b for clades than species (species: n = 230, adjusted R 2 = 0.07; clades: n = 28, adjusted R 2 = 0. 86). This suggests that other explanatory variables affecting head motion and/or SDS performance (for example, behavioural ecology 35 , bauplan 27 and visual acuity 36 ) drive morphological variation at the species level and blur the correlation (Supplementary Note 2). Consequently, reversing the relationship to predict T b from the TMI gives more accurate results for clades than for species (species: n = 230, adjusted R 2 = 0.01; clades: n = 28, adjusted R 2 = 0.86; Supplementary Note 2). Notably, the relationship between the TMI and T b shows a steeper increase than would be theoretically expected from compensation only for temperature-induced changes in viscosity (Fig. 2b, dashed curve). This steeper-than-expected increase probably reflects a positive correlation between overall head motion and T b (Supplementary Note 2). This corollary is congruent with the aerobic capacity model 34 in that elevated T b in endotherms may not have evolved in isolation, but in tandem with increased aerobic capacity, translating into increased behavioural activity. At high-ranking taxonomic levels, T b may explain the size differences observed among the SDS of vertebrates 27 when body size is accounted for.
The TMI of extant endotherms is significantly higher than that of extant ectotherms (PGLS ANOVA, endotherms: n = 143, ectotherms: n = 127, adjusted R 2 = 0.06, Glass' Δ = 1.77 with 95% confidence interval [1.48-2.06], P = 2.57 10 −5 ; Supplementary Note 2). Fishes, amphibians and turtles show the lowest T b , which are reflected by their low TMI values (Fig. 2). Conversely, birds and mammals show the highest TMIs, in agreement with their higher T b . Crocodilians and lepidosaurs are intermediate between these extremes. In tetrapods, the TMI shows a significant negative correlation with canal cross-sectional thickness

Article
relative to the radius of curvature and the radius of curvature relative to body mass (Supplementary Note 2). Among extant species, mammals show the lowest values for relative cross-sectional thickness and relative radius of curvature (Extended Data Fig. 3), leading to high TMI values. However, birds have semicircular canals with large relative radii of curvature and relative cross-sectional thicknesses that are average for amniotes (Extended Data Fig. 3). In this context, to attain their high TMI values, birds responded differently to increased T b by modifying aspects of their membranous labyrinth and the physicochemical properties of their endolymph 37 (Extended Data Fig. 4) while maintaining plesiomorphic bony labyrinth dimensions. The TMI distributions of ectotherms and endotherms show little overlap (Fig. 3), indicating that this metric could be useful for assessing the thermoregulatory regime of extinct organisms.

Phylogenetic origin of mammalian endothermy
To understand the thermophysiology of fossil synapsids, we computed body temperature distributions from the TMI of fossil specimens, using Brownian motion simulations calibrated with extant data (Supplementary Note 2), and calculated the probabilities of endothermy for each fossil taxon and node by fitting a phylogenetic logistic regression to the TMI distribution of extant species (Fig. 3, Table 1, Extended Data Figs. 5 and 6 and Extended Data Tables 1-3). We obtained a cross-validated error rate of 5% by using this regression to classify the thermoregulatory regime of extant species and clades, when regarding probabilities for species (P > 0.70) and for clades (P > 0.53) indicating endothermy, and for species (P < 0.30) and clades (P < 0.47) implying ectothermy, as well as intermediate probabilities as expressing uncertainty. In this framework, maximum-likelihood ancestral state reconstruction shows that the clade threshold for endothermy is crossed at the Mammaliamorpha node (P = 0.70), and non-mammalian mammaliamorphs are the only non-mammalian synapsids classified as endotherms (Extended Data  Tables 1-3). This clade includes all descendants of the last common ancestor of tritylodontids and mammals 38 and its origin is currently calibrated at approximately 233 million years ago (Ma), during the Carnian pluvial episode. The TMI of non-mammalian mammaliamorphs (1.29 ± 0.18 (± s.d.), non-mammalian mammaliamorphs: n = 9) is intermediate between the TMI measured in extant Ornithorhynchus (1.25) and that of Tachyglossus (1.41). These monotremes are capable of producing their own body heat while inactive but show low basal metabolic rates, relatively low T b (31-34 °C) and periodic torpor 33,39 . Mammaliamorphs were probably characterized by this thermoregulatory strategy, defined as basoendothermy 1 , with a predicted T b of approximately 34 °C (95% confidence interval [31.3-36.6]; Table 1 and Extended Data Fig. 5), and behavioural activity that was probably intermediate between extant ectotherms and endotherms (Extended Data Fig. 7). Conversely, all non-mammaliamorph synapsid groups classify as ectotherms, with probabilities below the clade threshold for ectothermy (P = [0.13-0.36]), and predicted T b between 24-29 °C (95% confidence interval [21.9-31.6 °C]; Table 1 and Extended Data  Table 3), a range that overlaps with extant turtles, lepidosaurs and crocodilians.

Caveats and literature comparisons
Predictions of body temperature and thermoregulatory regime are subject to uncontrolled parameters affecting the TMI (that is, endolymph chemistry and residual variation of bony/membranous correlations). However, these parameters remained stable on average along the mammalian lineage, between the nodes of Amniota and Mammalia (Extended Data Fig. 4 and Supplementary Note 2). This suggests that thermoregulatory regimes estimated at different nodes within the mammalian lineage are accurately predicted by the bony labyrinth of fossil synapsids. However, we cannot exclude that outbranching synapsid clades (for example, bidentalian dicynodonts 6,10 ) could have independently evolved much thinner membranous semicircular ducts than expected from their bony canals (for example, Phocidae 40 ) and/ or modified endolymph chemistry toward increased viscosity (for example, European plaice 23 ) to maintain SDS function. Similarly, we cannot exclude that these organisms could have been much less active than is typical for their size (for example, Zaglossus 41 ). Each of these scenarios would have resulted in T b at the upper tail of distributions predicted for these clades (Extended Data Fig. 5), generally below 31 °C and significantly lower than predicted for Mammaliamorpha (Table 1). Our sample includes some taxa analysed in previous studies using stable oxygen isotopes 10 and osteohistologically inferred metabolic rates 6 . Our evidence is largely consistent with these studies in interpreting sphenacodontians, biarmosuchians, dinocephalians, basal dicynodonts and therocephalians as ectotherms 6,10 , and in predicting increased T b in cynognathians 10 (28.9 °C, 95% confidence interval [26.3-31.6 °C]). As mentioned previously, interpreting results from the TMI of individual fossil taxa rather than clades can be misleading. For example, although the clade Cynognathia is clearly interpreted as ectothermic in our analysis (P = 0.36), probabilities of individual species do not enable unambiguous classification (Extended Data  Tables 1 and 3).  Numbers in square brackets represent the confidence intervals of predicted T b . Numbers in bold represent probabilities significantly excluding T b higher than or equal to 31 °C, or higher than or equal to T b of non-mammal mammaliamorphs obtained during the same simulations.

Tempo and mode in evolution of endothermy
To investigate evolutionary tempo, we fitted different evolutionary models to the TMI of synapsids (for example, Brownian motion, continuous trend, early burst; Supplementary Note 2). The best-supported model corresponds to partial Brownian motion (that is, Pagel's λ = 0.5), with a major shift in evolutionary rate in the branch leading to Mammaliamorpha (Fig. 4). This shift is interpreted as the transition to endothermy, with a change in T b of +5-9 °C (Table 1) compared with non-mammaliamorph eucynodonts. Notably, all models incorporating gradual evolution of the TMI are rejected, as is a change of evolutionary rate after the end-Permian mass extinction, despite other major effects of this extinction 42 . The acquisition of endothermy in Mammaliamorpha seems to be temporally and phylogenetically consistent with a pulse in brain enlargement 17 , the development of respiratory turbinals 7 , the evolution of infraorbital canals (related to the presence of vibrissae 16 ) and functional differentiation of the vertebral column 43 , and occurred during the miniaturization process that started in probainognathians 1,44 . However, the acquisitions of parasagittal posture 3 , the diaphragm 5 , microvascularization 4 , fibrolamellar bone 13 and parental care 9 might all have predated the onset of endothermy. These features probably mark gradual stochastic increases in T b (up to +3-4 °C in eucynodonts, similar to the difference between crocodiles and geckos) and behavioural activity during the evolution of non-probainognathian synapsids. Although our sample is small (n = 3), our simulations suggest that this gradual stochastic increase stopped in non-mammaliamorph probainognathians (Supplementary Note 2), which probably reverted to plesiomorphic synapsid body temperatures (Table 1). This may be explained by the onset of miniaturization 44 , coupled with exploration of low-temperature niches (for example, nocturnality or mesic environments). Combining our results on TMI evolution with physiological and experimental information from extant analogues, we posit that endothermy emerged from the combination of heat production through sarcolipin-mediated thermogenesis and thermal insulation via pelage, which seem essential in mammals to reach and maintain high T b when exposed to cold 31,45,46 . The sln gene, which controls sarcolipin-mediated thermogenesis, is synapomorphic for vertebrates but has been found to be expressed in mammals, the opah and is suspected to be expressed in birds 46 . Fossil evidence from early Late Jurassic docodonts (for example, Castorocauda) and haramiyids (for example, Maiopatagium) suggests that pelage could have originated much earlier. Thus, although any of these features could have had an earlier origin, endothermy may have emerged from their combination in Mammaliamorpha 1,47 , providing them a selective advantage to cope with a progressively cooler Triassic climate 28 in the context of reduced thermal inertia from miniaturization 44 . Despite being spatiotemporally biased 48 , the fossil record indeed seems to suggest that local mammaliamorph species richness increased during the Rhaetian cool interval 49 , the coldest time during the Triassic period.
thermoregulatory strategies to some extent. Furthermore, they must have resorted to long-term torpor during unfavourable periods of the cold season, congruent with histological evidence of cyclical growth 13 and burrowing 50 . Similarly, during the exceptionally high temperatures of the Early Triassic period, sampled non-eucynodont therapsids were probably forced to aestivate 50 and/or to become nocturnal.

Conclusions
The thermal budget of vertebrates results from a complex interplay between internally generated heat and heat gained or lost at the body surface. Although inferring the thermoregulatory regime of extinct species can be complicated for these reasons, we use the TMI to demonstrate that a sharp shift to endothermy most probably occurred at the base of Mammaliamorpha, accompanied by expanded aerobic and anaerobic capacities and increased body temperatures. The Late Triassic emergence of the physiology and bauplan characterizing mammals was a decisive step that facilitated their evolutionary radiation during the Mesozoic Era and subsequent major ecological expansion during the Cenozoic Era.

Online content
Any methods, additional references, Nature Research reporting summaries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contributions and competing interests; and statements of data and code availability are available at https://doi.org/10.1038/s41586-022-04963-z. Publisher's note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Data acquisition and processing
The SC endocasts in the 3D bony sample were downloaded from https://www.morphosource.org/ or segmented using Amira 5.3.3 (Visage Imaging and Konrad-Zuse-Zentrum) from existing and new computed tomography (CT) images obtained with μCT, or propagation phase-contrast synchrotron microtomography (Extended Data Fig. 1). Segmentation was based on the contrast between bone and air for extant specimens, but was conducted manually for fossils. A wax endocast model of Dimetrodon sp. (FMNH PR 4976) was digitized using photogrammetry (Supplementary Note 3). The endocasts of SDs and SCs in the 3D membranous sample were segmented from μCT scans of the ear region of extant specimens, acquired from autopsies or museum collections and prepared as described previouslymembranous labyrinths of 40 extant fish. Segmentation and additional processing 19 was done in Avizo 7.1 (Visualization Science Group) and Geomagic Studio 12 (Raindrop Geomagic). Further information about the specimens and scans are given in Supplementary Data 3.

Measurements of 3D specimens
Four linear measurements were taken per SC for all specimens in the 3D bony and membranous samples: two-dimensional length of the slender portionmembranous labyrinths of 40 extant fish, cross-sectional diameter and major and minor axes of each SC (Extended Data Fig. 2). These functionally relevant measurements mirror SD metrics and were designed for ease of measurement. The length of the slender portion is taken on the SD midline, which is closer to the outermost wall of the SC 19,62,63 , from the ampullar junction of the SD to the common crus or vestibule (Extended Data Fig. 2). The distal part of the slender portion of the lateral canal is often fused with the vestibule, especially in ectothermic taxa, and ends at the common crus wall. The major and minor axes are taken perpendicular to each other, following the ellipse that best fits the SD torus. They do not follow anatomical landmarks but their endpoints are placed on the SD midline (Extended Data Fig.  2). The radius of curvature and eccentricity of each SC torus was calculated from the semi-minor and semi-major axes, respectively taken as half the measured minor and major axes (Supplementary Methods, 'Biomechanics'). Measurements were taken five times and averaged to reduce error. Measurements of major and minor axes of the SC were taken using the 'straight line' tool of ImageJ 64 , on scaled screenshots from Amira or Avizo, where the screen plane was aligned to the plane of each SC. Measurements of lengths of slender portions 19 were taken on the same screenshots, using the 'segmented line' tool of ImageJ. SC planes were approximated visually before taking screenshots. A comparison with a more accurate method, which uses landmarks placed along each canal to find the best fitting SC planes, was made in XLSTAT 2018. Morphological and functional analyses of the 3D membranous sample were performed using Ariadne Toolbox 19 . Input data included surfaces and linesets representing individual components of the duct system (for example, slender parts and ampullae, and modelled cupulae 19 ). The average height of each cilia area was measured on the μCT scans. Endolymph density was taken as 1,000 kg m −3 , Poisson ratio as 0.48 and cupula shear modulus as 1.44 Pa 19 . Output parameters included, among others, wall shape factor, 3D length and cross-sectional area of the slender portion of each SD, enclosed area of the projection of each SD torus on its maximal response plane, and a transfer factor linking endolymph volume displacement to cilia deflection (hereafter called cilia deflection factor). The cilia deflection factor corresponds to the ratio between average deflection of the cupula at the level of the height of the underlying cilia to the average endolymph volume displacement for a given angular head rotation. For specimens for which we have membranous labyrinths, the deflection of the cupula at the level of the height of the underlying cilia has been simulated through finite element analysis 19 . The cupula, generally not visible, has been modelled as an extrusion of the crista ampullaris towards the roof of the cupula. Such modelling of cupula anatomy has been confirmed to closely match the actual morphology of the cupula in Saimiri sciureus (see Ariadne Manual Data Preparation 19 , pages 31-58). Cilia deflection factors of published specimens 51-55 were estimated from average cupula thickness and cross-sectional area of the ampulla above the crista, following allometric regressions (Supplementary Note 2, models 1-6).

Measurements of 2D images
The 2D membranous sample consists of published photographs 56-61 and measurements 27 of membranous labyrinths. We used ImageJ to process photographs and measure areas enclosed by SD tori, major and minor axes, lengths of slender portions, and inner duct radii (Supplementary Data 2). Photographs were generally taken in lateral view and show the planes of the anterior and posterior SDs at an angle, projected onto the sagittal plane. The resulting distortion can be corrected by multiplying the horizontal axis of the image by 1/cosine of the angle of a specific duct, providing a good representation of undistorted size and shape. Where available, the angles of the planes of the anterior and posterior ducts to the sagittal plane were measured in 'dorsal' view. For specimens without a 'dorsal' view, we generally used an angle of 45° to undistort the views (Supplementary Data 2). The undistorted views of the anterior and posterior ducts were used to measure corresponding SDs. When a 'dorsal' view was available, measurements for the lateral SD were done on this view. For specimens without a 'dorsal' view, the maximum distance measured between any points on the lateral semicircular duct in lateral view was taken as its major axis. Missing measurements were estimated in R, using cross-validated partial least-square regressions on a set with all measurements (Supplementary Data 2).

Body temperature, body size and phylogeny
Body temperatures of extant species were obtained from the literature (Supplementary Data 2), supplemented by Traitbank (https://eol.org/ traitbank). For ectothermic taxa, preferred T b were chosen. When both activity and preferred ranges of T b were provided, we took the average of the latter. In a few instances (11 ectotherm species) we used online pet care information to estimate preferred T b . Values for fish species were obtained from FishBase (http://www.fishbase.org).
Three body size variables were used: condylobasal length, condyloanteroorbital length and body mass. Condylobasal length (posteriormost border of the occipital condyle(s) to the anteriormost tip of the snout, projected onto the sagittal plane) was used to represent cranial size and could be measured for most specimens (Supplementary Data 2). To avoid bias introduced by snout lengths that are either exceptionally long (for example, gharials) or short (for example, anomodonts), we defined the condylo-anteroorbital length as the linear distance from the posteriormost border of the occipital condyle(s) to the anteriormost border of the orbit, projected onto the sagittal plane (see measurements in Supplementary Data 2). Body mass was measured whenever possible (Supplementary Data 2). We used typical adult values rather than extremes, and male and female values were averaged for sexually dimorphic species. For fossils and some specimens of the 3D bony sample, body mass was predicted using an allometric multiple regression of body mass and condylobasal and condylo-anteroorbital lengths (Supplementary Note 2, model 7). Predictions were refined by computing phylogenetically informed body mass residuals.
The phylogeny used in this study (Supplementary Data 6) was built in Mesquite 3.5 68 , using relationships and divergence dates between extant taxa provided by TimeTree 69 (http://www.timetree.org) as a backbone to which extinct clades and species were connected. Relationships at the level of Neoaves were modified from TimeTree to fit assumptions of other published accounts 70,71 . Relationships and divergence dates of actinopterygians and chondrichthyans that were not covered by TimeTree were obtained from the literature [72][73][74] . The reasoning behind phylogenetic relationships, divergence dates, and last occurrence data is detailed in Supplementary Note 4 and follows best practices 75 .

Defining the thermo-motility index
SDs are filled with a fluid (endolymph), whose inertia leads to the deflection of the cupula during head rotation. This deflection generates a signal that is integrated in the brain and is key to essential body functions 18,19-22 . SDs are angular velocity transducers that are most effective when the signal is unsaturated 76 and in phase with angular velocity of head motion 18 . This only happens for a specific detection range of angular velocities and a specific frequency bandwidth of rotations 18 , each bracketed by lower and upper limits. Consequently, the angular velocity range and frequency spectrum of naturally occurring angular head motion is expected to fall within the detection range and frequency bandwidth of SDs. Here we particularly focus on the maxima of angular head motion and semicircular duct function (that is, upper corner frequency and saturating angular velocity, Supplementary Methods, 'Biomechanics'), because it is likely that accurately capturing high angular head motion during the most intense locomotory behaviours is most important for selective fitness of the organism.
Body temperature directly affects the viscosity of the endolymph, such that it decreases when T b increases 23,24 . For an organism with a given endolymph composition and SDS morphology, an increase in T b will result in a reduction of the upper corner frequency and saturating angular velocity. If the SDS morphology does not change or endolymph composition remains stable, the upper limits of angular head motion frequency or velocity will likely fall outside the range of effective SDS function, affecting adversely the selective fitness of the organism. Alternatively, the organism could adapt to more sluggish behaviours, but empirical evidence suggests that, if anything, increased T b correlates to increased locomotor activity 25,26,34 (Supplementary Data 1 and Supplementary Note 2). Thus, analysing the maximum SDS function of extinct species should allow us to track changes in T b , and the acquisition of endothermy. Therefore, we derive the TMI from SDS morphology and endolymph composition, using equations governing the saturating angular velocity and the upper corner frequency, but excluding T b (Supplementary Methods, 'Biomechanics'). Consequently, when T b increases, compensatory changes in the TMI, traceable in fossils, are required to keep semicircular duct function constant.
In practice, as fossils do not preserve the membranous labyrinth, we need to predict the morphology of their SDS. To do so, we first used a dataset of 50 species (Supplementary Data 2), for which we had bony and membranous labyrinths, to regress morphometrical parameters of the ducts (including their cross-sectional area) against morphometric parameters of the bony canals (Supplementary Note 2, models 8-19). Next, we used these regressions to obtain fitted values of these duct parameters using morphometric parameters of the bony canals of fossil specimens. Finally, we estimated the residual parameters of the SDs for the fossil specimens through ancestral state estimation, using the phylogeny mentioned above. Similarly, direct information on endolymph composition is lacking in fossils. Data on endolymph viscosity at a given temperature has been published for only eight species (one bird, two mammals, five euteleosts) 23,37,77-82 . Using published data, we calculated the physicochemical component of the TMI of these species (Supplementary Methods, 'Biomechanics') and found that, except for Columba livia and Pleuronectes platessa, they have endolymph with low viscosity, close to that of water (Extended Data Fig. 4). The phylogenetic distribution of the physicochemical component of the TMI in available species indicates that low-viscosity endolymph was the basal condition for Euarchontoglires and Euteleostei, thus parsimoniously the basal condition for Osteichthyes as well. This phylogenetic context allows us to assume a priori that the physicochemical component of the TMI of synapsid species did not differ from the basal condition. Following these observations, a low-viscosity endolymph was chosen for all species analysed in this study, except birds and pleuronectids. For non-avian diapsids, amphibians, and chondrichthyans, theoretical and empirical relationships between the TMI and T b are consistent suggesting, a posteriori, the retention of low-viscosity endolymph in these taxa ( Fig. 2 and Supplementary Methods). Finally, because the TMI scales allometrically, it needs to be corrected for body size when comparing between species (Supplementary Note 2, models 26-31).
A more comprehensive introduction to the biomechanics of the TMI and its complete mathematical derivations are available in Supplementary Methods, 'Biomechanics'. Moreover, a step-by-step protocol describing how to obtain the TMI from bony or membranous morphology is also provided in the Supplementary Methods. A complete description of the statistics related to the computation of the TMI and related analyses is available in Supplementary Note 2.

Statistical analyses
Computation of the TMIs and statistical analyses were done in R using the packages phytools 0.7-70 83  To determine the morphological parameters of the bony anterior SC (for example, cross-sectional radius and length of the slender portion, radius of curvature and eccentricity) that carry the most information about the transition toward endothermy in the synapsid lineage, we regressed size-corrected morphological parameters on the TMI using PGLS regressions (Supplementary Note 2, models 49-56).
To assess which model best fit the evolution of the TMI, we fitted nine evolutionary models on the TMI (for example, Brownian motion, early burst) and compared their corrected Akaike information criteria (Supplementary Note 2, models [58][59][60][61][62][63][64][65][66]. Results suggest that an evolutionary model with Brownian motion along the proposed phylogenetic tree, with branch lengths scaled using a Pagel's λ of 0.50, and a major shift in the rate of evolution along the branch directly leading to Mammaliamorpha, is best for explaining the evolution of the TMI (Supplementary Note 2, model 66).
Phylogenetic logistic regressions were fitted to predict the thermoregulatory regime of extant amniote specimens from their TMI (Supplementary Note 2, model 57). The resulting phylogenetic logistic regression was used to predict probabilities of endothermy for each node of the tree and each fossil species, using the corresponding TMI (Fig. 3, Extended Data Fig. 6 and Extended Data Tables 1 and 2). The TMI for each node of the phylogenetic tree was obtained by maximum-likelihood reconstruction of ancestral states under Brownian motion, using TMIs of all tips and the phylogenetic tree reflecting the best evolutionary model. Accuracy of predictions was assessed 100 times by randomly selecting 151 extant species, fitting a phylogenetic logistic regression to predict their thermoregulatory regime from their TMI, using the fitted phylogenetic logistic regression to predict the thermoregulatory regime of 75 extant species that were not used to produce the phylogenetic logistic regression, and comparing outcomes with real thermoregulatory regimes. The thresholds for endothermy and ectothermy were determined by using the probabilities that gave an error rate of 5%.
Multiple additional tests were conducted to assess whether the TMI is free from body size information (Supplementary Note 2, models 38-39) and robust (Supplementary Note 2, models 43-47).
Probability distributions of predicted T b of fossil species and groups were obtained (1) by using the PGLS regression between T b (via the temperature ratio) and the TMI for extant clades (Supplementary Note 2, model 41), and (2) by simulating 15,000 times the evolution of residuals of the temperature ratio using Brownian motion on a phylogenetic tree reflecting the best evolutionary model for these residuals (Supplementary Note 2, section 24).
To assess whether our small non-mammaliamorph probainognathian sample could have affected our conclusions on the presence of a shift to endothermy at the root of Mammaliamorpha, we computed the likelihood of obtaining the mean TMI we report for this group by chance alone, when randomly sampling and averaging three specimens from distributions of higher TMI means (for example, cynognathian-like mean), using standard deviations observed in well-sampled clades (for example, mammals). This experiment has been replicated 10,000 times to obtain likelihood distributions (Supplementary Note 2, section 26).

Post hoc analyses
We investigated whether high/maximum stride and impact frequencies (as a proxy for behaviourally induced angular head motion 94 ) were positively correlated to T b , and if they explained part of the variance of the TMI left unexplained by T b . To do so, we performed PGLS regressions (Supplementary Note 2, models 84-89), using stride frequency, impact frequency, body mass, and T b , sampled for 106 taxa. Impact frequencies were obtained from stride frequencies, using multiplying factors estimated from the analysis of gait types (Supplementary Data 1).
Similarly, we assessed if maximum anaerobic speed (as a proxy for behavioural activity) is higher in endotherms than ectotherms, and if it is positively correlated to T b . To do so, we performed PGLS ANOVA and PGLS regressions (Supplementary Note 2, models 76-80), using maximum anaerobic speeds, body mass, and T b , sampled in 428 taxa. Data collection followed the same protocols as above and references are provided in Supplementary Data 1.
We explored which thermoregulatory strategies (for example, aestivation) could have been adopted by fossil synapsids predicted as ectotherms, by comparing their predicted T b to the palaeoclimate of their palaeolocation (Supplementary Data 4). To do so, we used the temperature curve of Scotese et al. 28 . to obtain minimum and maximum global average temperatures (GAT) corresponding to the age range of each fossil specimen. These GATs were corrected for the palaeolatitude of each specimen (obtained from the Paleobiology Database: https:// paleobiodb.org/), using Fig. 7 of Scotese et al. 28 , to obtain the latitudinal average temperature (LAT). In addition, we reviewed the literature to compile climate types that have been proposed for the palaeolocation of each sampled specimen [95][96][97][98][99][100][101] . Furthermore, we collected temperature data (for example, summer maximum temperature) for 150 locations and 16 different Köppen climates from https://weatherspark.com/. Climate type and LAT were used in R to predict seasonal temperatures at the palaeolocation of each fossil specimen, with a range of uncertainty. Palaeolocation temperature and predicted body temperature (Supplementary Note 2) were randomly sampled from corresponding ranges of uncertainty, and compared 1,000 times to address specific questions about thermoregulatory strategies (for example, whether predicted body temperature is higher than the maximum local night temperature-if true the specimen is probably obligatory diurnal). Results from these tests allowed us to compute probabilities for given thermoregulatory strategies to occur.

Reporting summary
Further information on research design is available in the Nature Research Reporting Summary linked to this paper.

Data availability
The raw datasets used in this study are available in the Supplementary Dataset. Links for CT scan datasets and bony labyrinth 3D meshes obtained from https://www.morphosource.org/ can be found in Supplementary Data 3. Some bird skull measurements and fish lengths were obtained from https://skullsite.com/ and fishbase.org, respectively (Supplementary Data 2). Time calibrations between most extant species were obtained from timetree.org. Body mass and T b of some extant species were obtained from https://eol.org/traitbank (Supplementary Data 2). Source data are provided with this paper.

Code availability
The R scripts used in this study are available in the Supplementary Dataset.