Non-technical Summary

Organismal metabolic rate is linked to environmental temperature and oxygen consumption and thus may be a useful predictor of extinction risk, especially during major climate-driven extinctions, given the combined stressors of warming and hypoxia. However, metabolic rate can be quantified in different ways, each with differing implications. Here, we estimate resting whole-body and mass-specific metabolic rates in bivalves using body size and temperature to assess whether metabolic rate predicts survival, particularly during intervals of rapid climate warming, or “hyperthermals.” Bivalves with lower total calorific needs are generally more likely to go extinct, which reflects an overall evolutionary increase in body size and energetics over geologic time. Conversely, bivalves with higher metabolic rates per gram of tissue are more likely to go extinct during hyperthermals relative to baseline conditions, particularly in well-sampled tropical environments, emphasizing the impact of warming on active animals living close to their thermal limits. When standardized by mass, metabolic rate may represent a useful metric with which to predict the effects of anthropogenic climate change on modern marine faunas.

Thermal Tolerance and Metabolism in Animals: The Conundrum

Tolerance of climate-related stress in experiments is highly variable across taxa and ecologies (Reddin et al. Reference Reddin, Nätscher, Kocsis, Pörtner and Kiessling2020), and explaining this variability is important for predicting how taxa will respond to anthropogenic warming. A variable often invoked to account for some of this variance is metabolic rate; however, metabolism is quantified in multiple ways in living organisms and can only be approximated for those in the fossil record, making evaluation of a potential link with survivorship during warming difficult to assess. An organism’s activity level, which can be observed or inferred for living and fossil taxa, is intuitively related to metabolic rate, because being active requires expending more energy. Peck et al. (Reference Peck, Clark, Morley, Massey and Rossetti2009) semiquantitatively measured activity levels through aspects of feeding and movement, finding that more active (and smaller) Antarctic bivalves are more likely to survive short-term warming during experiments. Clapham (Reference Clapham2017) found that more active marine organisms were also more likely to survive during the Permian–Triassic and Triassic–Jurassic extinctions, both trends associated with ocean warming, acidification, and hypoxia.

Metabolism is more specifically quantified in terms of an organism’s measured metabolic rate—its rate of energy expenditure per unit time. At higher temperatures, aerobic metabolism proceeds faster and requires more oxygen, with the temperature-dependent balance between oxygen supply and demand likely setting the thermal and hypoxic tolerances of organisms (Pörtner Reference Pörtner2010; Pörtner et al. Reference Pörtner, Bock and Mark2017). To project how anthropogenic warming might change the geographic distribution of metabolically viable habitat, Deutsch et al. (Reference Deutsch, Ferrel, Seibel, Pörtner and Huey2015) defined the metabolic index ( $ \Phi $ ) as the ratio of the temperature-dependent O₂ supply in the environment to an organism’s resting O₂ demand. Conditions where demand exceeds supply (i.e., $ \Phi $ is less than 1) are considered to be metabolically nonviable (Deutsch et al. Reference Deutsch, Ferrel, Seibel, Pörtner and Huey2015). Although organisms can survive these conditions temporarily, they must employ tolerance mechanisms, including anaerobic metabolic pathways, heat shock proteins, and suppression of aerobic activity, which are not energetically sustainable in the long term (Pörtner et al. Reference Pörtner, Bock and Mark2017). The metabolic index predicts, for example, that polar taxa are especially vulnerable to warming due to experiencing greater losses of metabolically viable habitat (Deutsch et al. Reference Deutsch, Ferrel, Seibel, Pörtner and Huey2015), a finding consistent with observations of increased extinction intensity at higher latitudes during the end-Permian extinction (Penn et al. Reference Penn, Deutsch, Payne and Sperling2018; Reddin et al. Reference Reddin, Kocsis and Kiessling2019). If metabolic rate and oxygen demand increase with warming while oxygen supply decreases, one might predict more active organisms to be more susceptible to warming, as their higher oxygen needs become increasingly difficult to meet. However, this prediction runs counter to observed patterns of selectivity based on activity level, where high-activity (and presumably high metabolic rate) taxa preferentially survive warming events (Peck et al. Reference Peck, Clark, Morley, Massey and Rossetti2009; Clapham Reference Clapham2017). However, more metabolically active organisms generally have higher CO₂ tolerance (Melzner et al. Reference Melzner, Gutowska, Langenbuch, Dupont, Lucassen, Thorndyke, Bleich and Pörtner2009), potentially making them less vulnerable to hypercapnia as a result of CO₂ buildup in the oceans during rapid climate warming (Knoll et al. Reference Knoll, Bambach, Canfield and Grotzinger1996). The multiple environmental stressors associated with hyperthermal conditions are complex, not independent, and can act on different temporal scales (Reddin et al. Reference Reddin, Nätscher, Kocsis, Pörtner and Kiessling2020). As such, the links between activity level, metabolic rate, and extinction risk during periods of rapid climate warming need greater clarity to help identify the mechanisms leading to organism vulnerability.

Quantifying Metabolism in Fossil Organisms

Although activity levels are related to metabolic rate, a roadblock to studying how metabolism influences extinction risk is our inability to directly measure metabolic rates in extinct organisms. Resting whole-body metabolic rate, in watts (joules per second), can be estimated from temperature and biomass using a model originally derived by Gillooly et al. (Reference Gillooly, Brown, West, Savage and Charnov2001) as follows: an organism’s whole-body metabolic rate ( $ B $ ) scales with its body mass ( $ M $ ) with an exponent close to 3/4, in other words, $ B\propto {M}^{3/4} $ . The relationship between body size, temperature, and resting whole-body metabolic rate can be approximated by the expression,

(1) $$ B={B}_0{M}^{3/4}{e}^{-{E}_i/ kT} $$

where $ {e}^{-{E}_i/ kT} $ is the Boltzmann’s factor, with $ {E}_i $ being the average activation energy of rate-limiting biochemical reactions, $ k $ being Boltzmann’s constant, and $ T $ being temperature in Kelvin. $ {B}_0 $ is a taxon-specific metabolic rate scaling constant, equivalent to the average resting metabolic rate of the higher-order taxon to which the organism belongs. This assumes that resting metabolic rate is generally taxonomically conserved, an assumption supported by, for example, Saupe et al. (Reference Saupe, Hendricks, Portell, Dowsett, Haywood, Hunter and Lieberman2014), who found that thermal tolerances in mollusks were stable over millions of years covering substantial environmental change. Note that the Gillooly et al. (Reference Gillooly, Brown, West, Savage and Charnov2001) model estimates resting metabolic rate and does not explicitly account for variations in activity levels or growth rate except coarsely through the taxon-specific (generally at the ordinal or class level) adjustment of the $ {B}_0 $ coefficient. Previous applications of the Gillooly et al. (Reference Gillooly, Brown, West, Savage and Charnov2001) model to the fossil record have generally focused on the total energy use of populations and ecosystems, without reference to extinction risk. The three-quarters power scaling of metabolic rate with body mass was used by Finnegan and Droser (Reference Finnegan and Droser2008) to estimate population energy use, with temperature held constant. Finnegan et al. (Reference Finnegan, McClain, Kosnik and Payne2011) assessed changes in the energy consumption of marine gastropod assemblages over the Mesozoic marine revolution and found that mean per capita metabolic rates of gastropod assemblages rose between the Triassic and Cretaceous, largely driven by increases in body size. Payne et al. (Reference Payne, Heim, Knope and McClain2014) used the Gillooly et al. (Reference Gillooly, Brown, West, Savage and Charnov2001) model to establish that, despite being less diverse and abundant, bivalves accounted for a much greater share of energy metabolism than brachiopods in Paleozoic oceans.

Whole-body metabolic rate can be divided by the mass of the organism to yield mass-specific metabolic rate ( $ B/M $ ), the metabolic rate of one unit mass of tissue. Mass-specific metabolic rate allows for comparisons of energy expenditure across orders of magnitudes of variation in body size (Makarieva et al. Reference Makarieva, Gorshkov, Li, Chown, Reich and Gavrilov2008). Mass-specific resting metabolic rate can be approximated by the expression,

(2) $$ B/M={M}^{-1}\bullet {B}_0{M}^{3/4}{e}^{-{E}_i/ kT}={B}_0{M}^{-1/4}{e}^{-{E}_i/ kT} $$

Whole-body metabolic rate $ \propto $ (is proportional to) $ {M}^{3/4} $ , whereas mass-specific metabolic rate $ \propto {M}^{-1/4} $ . Assuming consistent values for $ {B}_0 $ and $ T $ , a larger organism will have a higher whole-body metabolic rate, but a smaller mass-specific metabolic rate. Whole-body and mass-specific metabolic rates should therefore be inversely correlated with one another. Given that mass-specific metabolic rates are consistently optimized across the major domains of life (Makarieva et al. Reference Makarieva, Gorshkov, Li, Chown, Reich and Gavrilov2008), we also assess mass-specific metabolic rate as a potential predictor of extinction selectivity.

Body size is the primary factor governing resting whole-body metabolic rate (Gillooly et al. Reference Gillooly, Brown, West, Savage and Charnov2001). Hence, we expect an overall pattern of extinction selectivity similar to that found by Payne and Heim (Reference Payne and Heim2020), wherein organisms with larger body sizes and thus higher whole-body metabolic rates are generally less likely to go extinct. Reddin et al. (Reference Reddin, Kocsis, Aberhan and Kiessling2021) observed extinction selectivity for body size in line with Payne and Heim (Reference Payne and Heim2020) but detected a relative shift against larger body sizes in bivalves during hyperthermal-associated stages. Peck et al. (Reference Peck, Clark, Morley, Massey and Rossetti2009) likewise found in experiments on modern marine animals that smaller species were able to survive to higher temperatures than larger ones, Thus, we expect similar relative shift against higher metabolic rates in bivalves during hyperthermal-associated stages.

In contrast, mass-specific metabolic rate may provide a unique insight into how physiology can influence extinction risk. The narrow range of variation in $ B/M $ across the major domains of life suggests maintenance by strong selective pressures (Makarieva et al. Reference Makarieva, Gorshkov, Li, Chown, Reich and Gavrilov2008), potentially including extinction. Per the oxygen- and capacity-limited thermal tolerance hypothesis, the difference between maximum and resting metabolic rate, or aerobic scope, is highest within species-specific optimum temperatures. While the metabolic rate of an organism at rest is generally far below its maximum, at high temperatures, an organism’s resting metabolic rate approaches this maximum (Pörtner Reference Pörtner2010). Therefore, organisms living in warmer environments, which force higher mass-specific metabolic rates, may already live close to their upper pejus temperatures, where physiological performance rapidly decreases (Pörtner et al. Reference Pörtner, Bennett, Bozinovic, Clarke, Lardies, Lucassen, Pelster, Schiemer and Stillman2006). With such narrow thermal safety margins, further increases in temperature such as heat waves could have disastrous effects on survival. Penn and Deutsch (Reference Penn and Deutsch2022) found that projected extirpation risk in modern marine organisms is higher where species are living close to their ecophysiological thresholds ( $ \Phi $ = 1). Additionally, due to the exponential scaling of basal metabolic rate with temperature, a 1°C increase in environmental temperature would more intensely shift the metabolic rate of an organism that already lives at a warm baseline temperature (Dillon et al. Reference Dillon, Wang and Huey2010). In sum, the implication of the temperature dependency of mass-specific metabolic rate is that a higher mass-specific metabolic rate has more severe implications for population survival during rapid warming, which increases risk of extinction during hyperthermals.

Methodological Approach

We downloaded a dataset of Phanerozoic marine animal occurrences from the Paleobiology Database (PBDB, www.paleobiodb.org) on October 18, 2023. We binned each occurrence to its respective chronostratigraphic stage, as detailed in the vignette accompanying the R package divDyn (Kocsis et al. Reference Kocsis, Reddin, Alroy and Kiessling2019), and excluded any unbinned occurrences. We then calculated three-timer sampling completeness of Alroy (Reference Alroy2008) and range-based sampling probability using divDyn (Kocsis et al. Reference Kocsis, Reddin, Alroy and Kiessling2019), removing any singletons (taxa with fewer than two occurrences). As the order Unionida is nearly completely composed of freshwater bivalves, all genera from this order were excluded from the dataset. We then calculated genus stratigraphic ranges, including last-appearance date, by stage using divDyn::modeltab().

Because geographic range size has a known effect on extinction risk (e.g., Payne and Finnegan Reference Payne and Finnegan2007), we incorporated range size as a genus’s subsampled maximum great-circle distance (GCD) (Foote et al. Reference Foote, Ritterbush and Miller2016). Specifically, we took the mean over 100 iterations of the largest GCD between each genus’s occurrences in a random subset of 600 occurrences from a stage. Because geographic range size is strongly skewed, we followed Reddin et al. (Reference Reddin, Kocsis, Aberhan and Kiessling2021) to bin geographic range into ordinal levels based on the GCD values of each stage. Here, the first bin represents genera with no empirical occurrences during a particular stage (i.e., those interpolated by range-through occurrences), the second representing genera with maximum GCDs of 0 (i.e., those with only a single locality during the stage), the third representing everything else up to the third GCD quartile, and the fourth representing everything above the third quartile. Removing geographic range size from the analysis did not significantly affect our results.

Body-size estimates at the genus level were provided by the supplementary data files of Payne and Heim (Reference Payne and Heim2020), as measurements of the longest linear dimension (in millimeters) of the largest specimen of each genus. Maximum linear dimensions of bivalve shells are strongly correlated with both shell volume and biomass (Powell and Stanton Reference Powell and Stanton1985). Per Payne et al. (Reference Payne, Heim, Knope and McClain2014), bivalve ash-free dry mass ( $ M $ ), in grams, scales with the maximum linear dimension of the shell ( $ L $ ) via the following function:

(3) $$ M=1.0\times {10}^{-5}\times {L}^{2.95} $$

Thus, we estimated $ M $ for each genus (in grams) using measurements of maximum shell length.

Based on the compilation of Brey (Reference Brey2001), Payne et al. (Reference Payne, Heim, Knope and McClain2014) used $ {B}_0 $ constants of 1.4 × 10¹¹ W kg^−3/4 for heterodont bivalves, and 1.3 × 10¹¹ W kg^−3/4 for all other bivalves. We follow this approach, as did Strotz et al. (Reference Strotz, Saupe, Kimmig and Lieberman2018). We define heterodont bivalves as those included in the molecular phylogeny of Taylor et al. (Reference Taylor, Williams, Glover and Dyal2007) along with the extinct orders Actinodontida, Hippuritida, and Megalodontida.

Paleotemperature estimates were obtained for each taxon based on its observed distribution over paleogeographic maps. The modern latitude and longitude coordinates associated with each bivalve occurrence listed in the PBDB were used to reconstruct paleolatitude and paleolongitude coordinates for each occurrence via the R package rgplates (Kocsis et al. Reference Kocsis, Raja, Williams and Dowding2024) with the GPlates desktop application (Müller et al. Reference Müller, Cannon, Qin, Watson, Gurnis, Williams, Pfaffelmoser, Seton, Russell and Zahirovic2018), using the PaleoMAP model (Kocsis and Scotese Reference Kocsis and Scotese2021). We used these paleo-coordinates to retrieve the corresponding annual mean sea-surface paleotemperature outputs from HadCML3 climate models run every 5 Myr (data from Valdes et al. [Reference Valdes, Armstrong, Badger, Bradshaw, Bragg, Crucifix and Davies-Barnard2017]; interpolated by Kocsis et al. [Reference Kocsis, Reddin, Scotese, Valdes and Kiessling2021]). We then calculated the median of all paleotemperature estimates associated with each bivalve genus, irrespective of time, as a relative measure of its thermal preference, which can be done for nearly the entire dataset.

In our estimates of metabolic rate, we keep $ {E}_i $ (the average activation energy of rate-limiting biochemical reactions) constant, at 0.65 eV, an approach also used by Finnegan et al. (Reference Finnegan, McClain, Kosnik and Payne2011) and Payne et al. (Reference Payne, Heim, Knope and McClain2014). This is likely an oversimplification: as discussed in Deutsch et al. (Reference Deutsch, Ferrel, Seibel, Pörtner and Huey2015), the temperature sensitivity of metabolic rate (in other words, $ {E}_i $ ) does vary among taxa and physiologies. However, accounting for this variation is not feasible within the scope of this study, so we assume that $ {E}_i $ is conserved within Bivalvia.

We conduct a series of logistic regression models of bivalve genus extinction odds as a function of log metabolic rate and geographic range size, by stage, for each post-Carboniferous stage. To account for the effects of sampling bias on the empirical selectivity of extinction, following Reddin et al. (Reference Reddin, Kocsis, Aberhan and Kiessling2021), the estimated stage of extinction is considered probabilistic (survival represented by 0, extinction represented by value 0.5 < x < 1) rather than definite (survival represented by 0, extinction represented by 1), based on the mean of the three-timer sampling completeness of the temporal stage and of the clade over all stages. The remainder extinction probability (1 − x) falls into the following stage, acknowledging that the youngest known fossils typically fail to represent the true timing of extinction (Signor and Lipps Reference Signor and Lipps1982; Marshall Reference Marshall1990). The resulting extinction selectivity coefficients indicate the mean effect of increasing metabolic rate on extinction risk, in log-odds, during a given stage. Positive values indicate that increasing metabolic rate increases the risk of extinction, while negative values indicate that bivalves with higher metabolic rates are less likely to go extinct than those with lower metabolic rates.

Because clades and traits have different background rates of extinction (Stanley Reference Stanley1985; Raup and Boyajian Reference Raup and Boyajian1988), we also calculate relative hyperthermal vulnerability (RHV), following Reddin et al. (Reference Reddin, Kocsis, Aberhan and Kiessling2021), to quantify the significance of any hyperthermal-unique change in extinction rate relative to the group’s background rate of extinction. To do this, extinction selectivity coefficients from the logistic regression were synthesized via meta-regression for stages associated with the onset of hyperthermal conditions and, separately, for all other (non-hyperthermal) stages. The use of meta-regression separates variance in the relationship between metabolic rate and extinction risk within a stage and among stages. Stages associated with the continuation rather than onset of hyperthermal conditions were excluded from the calculation of RHV, as the first pulse of hyperthermal-associated extinction is likely to be the most analogous to modern climate change. For the purpose of this study, we focused on the six major post-Paleozoic hyperthermal events identified by Foster et al. (Reference Foster, Hull, Lunt and Zachos2018), these being the Permian–Triassic (~252 Ma), Triassic–Jurassic (~201.3 Ma), the end-Pliensbachian–early Toarcian (~183 Ma), the Aptian (~120 Ma), Cenomanian–Turonian (~94 Ma), and Paleocene–Eocene thermal maximum (~55.5 Ma). We calculated RHV as the difference between the average metabolic rate–based extinction selectivity of stages associated with hyperthermal onsets and of stages not associated with hyperthermals. We follow Reddin et al. (Reference Reddin, Kocsis, Aberhan and Kiessling2021) by using 90% confidence intervals when calculating RHV, as the small sample size of six hyperthermals means that the likelihood of a type II statistical error is high relative to the type I error when using 95% confidence intervals (Cohen Reference Cohen1982). Regardless, the significance of our result is preserved even when a 95% confidence interval is used instead (Supplementary Table S2). Positive RHV values indicate that selectivity against higher metabolic rates is stronger during hyperthermal conditions than during other intervals. We follow this procedure for both whole-body and mass-specific metabolic rates. Because latitudinally selective extinction patterns have been detected over hyperthermal events, only some aspects of which may be linked to thermal stress (Penn et al. Reference Penn, Deutsch, Payne and Sperling2018; Reddin et al. Reference Reddin, Kocsis and Kiessling2019, Reference Reddin, Aberhan, Raja and Kocsis2022), we repeat the analysis restricting the data solely to tropical bivalve genera (i.e., those with mean paleolatitude of occurrences falling between −30° and 30°), and again for polar genera (mean paleolatitude of occurrences falling outside −60° and 60°). All analyses were performed in R using the packages metafor (v. 3.8-1) (Viechtbauer Reference Viechtbauer2010), divDyn (v. 0.8.2) (Kocsis et al. Reference Kocsis, Reddin, Alroy and Kiessling2019), and lme4 (v. 1.1-35.1). Figures were created using the packages ggplot2 (v. 3.4.2) (Wickham Reference Wickham2016) and deeptime (v. 1.0.1) (Gearty Reference Gearty2024).

An Alternative Source of $ {\mathbf{B}}_{\mathbf{0}} $ Values

The $ {B}_0 $ standard metabolic rate coefficient is a taxon-specific normalization constant that attempts to account for physiological differences between different higher-order clades. Payne et al. (Reference Payne, Heim, Knope and McClain2014) and Strotz et al. (Reference Strotz, Saupe, Kimmig and Lieberman2018) used $ {B}_0 $ constants of 7.87 × 10⁸ W g^−3/4 for heterodont bivalves, and 7.31 × 10⁸ W g^−3/4 for all other bivalves, based on the compilation of Brey (Reference Brey2001). An alternative source of standard metabolic rate constants for gastropods was used by Finnegan et al. (Reference Finnegan, McClain, Kosnik and Payne2011), taking the average $ {B}_0 $ of multicellular ectotherms of 9.91 × 10⁷ W g^−3/4 given by Gillooly et al. (Reference Gillooly, Brown, West, Savage and Charnov2001) and multiplying by experimentally derived standard metabolic coefficients for higher-order gastropod taxa provided by Vladimirova (Reference Vladimirova2001). Strotz et al. (Reference Strotz, Saupe, Kimmig and Lieberman2018) applied this method instead in their analysis of Pliocene to Recent gastropods. Equivalent standard metabolic coefficients for bivalves are provided by Vladimirova et al. (Reference Vladimirova, Kleimenov and Radzinskaya2003), and applying these coefficients does not affect the overall patterns of extinction selectivity based on whole-body and mass-specific metabolic rates (i.e., the trends through time and with latitude). However, estimating resting metabolic rate using $ {B}_0 $ coefficients from Vladimirova et al. (Reference Vladimirova, Kleimenov and Radzinskaya2003) removes the difference between hyperthermal and baseline extinction selectivity based on mass-specific metabolic rate. $ {B}_0 $ values sourced from Brey (Reference Brey2001) and Vladimirova et al. (Reference Vladimirova, Kleimenov and Radzinskaya2003) show no correlation (Supplementary Fig. S7), meaning that that these two sources likely measure different things.

Because whole-body metabolic rate is so heavily dependent on body mass, the sourcing of $ {B}_0 $ values was of little consequence for previous work focusing on whole-body metabolic rate. However, because the alternative metric of mass-specific metabolic rate significantly reduces the importance of body mass in the Gillooly et al. (Reference Gillooly, Brown, West, Savage and Charnov2001) equation, the influence of the $ {B}_0 $ coefficient increases. Brey (Reference Brey2001) consists of a compilation of resting metabolic rates for marine invertebrates where $ {B}_0 $ values are determined by controlling for latitude, depth, and body mass, ultimately splitting the dataset into heterodont and non-heterodont bivalves. Vladimirova et al. (Reference Vladimirova, Kleimenov and Radzinskaya2003) instead combined experimental rates of oxygen consumption for select bivalve taxa with published data on energy metabolism to obtain clade-specific mass-dependent coefficients of resting metabolic rate. $ {B}_0 $ values from Vladimirova et al. (Reference Vladimirova, Kleimenov and Radzinskaya2003) provide greater taxonomic resolution that might better reflect variation in resting metabolic rates within Bivalvia. Using these values results in the $ {B}_0 $ coefficient having a much larger effect on mass-specific metabolic rate (Table 1). However, the relatively low number of measurements combined with the limited number of orders and families represented means that we have less confidence in assigning $ {B}_0 $ values from Vladimirova et al. (Reference Vladimirova, Kleimenov and Radzinskaya2003) in our dataset.

Given the uncertainty associated with assigning $ {B}_0 $ values, we can alternatively ignore the effects of the taxon-specific normalization coefficient altogether and simply run the analysis holding $ {B}_0 $ constant for all bivalves. This approach yields a pattern of RHV resembling that obtained from our initial approach using values from Brey (Reference Brey2001) for heterodont and non-heterodont bivalves, with mass-specific metabolic rate being a significant predictor of hyperthermal vulnerability. Future analyses must take care when choosing which method to employ in estimating resting metabolic rate. Additional experimental constraints on $ {B}_0 $ for a range of taxa would be welcome in the effort to link extinction risk to metabolic rate.