Does a species thermal bias predict its removal under climate change?
Regional species loss is often correlated with climatic changes 13 but considering climate change relative to species’ thermal niches leverages additional information to assess population vulnerability 28,31(this study). We present empirical evidence from the fossil record that immigration, persistence, extirpation, and likely the extinction of species form an escalating response gradient linked to species suitability to regional conditions, as estimated by their thermal bias. A species’ thermal bias is thus a useful attribute to predict its likely response to warming, though it is likely insufficient alone (e.g. here R2 = 18%), with magnitude of regional warming and taxonomic membership explaining additional variation. Thermal biases of individual species were highly variable across responses despite the spatiotemporal extent and focal species being well-sampled, which supports the validity of responses such as extirpation. In some times and regions, seawater anoxia likely overruled the importance of temperature in determining assemblage membership. The otherwise consistent temperature–response relationship supports cautious use of habitat suitability models based on temperature, ideally alongside additional niche variables 31, to provide evidence on extinction risk as a statistical tendency over multiple species 38.
For simplicity, especially given the large time scales, we assumed a linear relationship between a species’ thermal bias and its regional occupancy response, from immigration, through persistence, to extirpation detailed in 10, and potentially extinction. This assumed an open thermal landscape over which species could disperse, which was relatively well supported, though sampling was bookended between approximately 15 and 34°C (Fig. S4), with geographical barriers to the north (discussed below). Non-temperature habitat (and sampling) heterogeneity is likely to dominate at scales beneath our effective spatial resolution of 1000s km. As climate changes, species distributions should move through a region, following shifting isotherms according to their thermal tolerance, other habitat variables permitting (see below) 10. Though we observed occupancy responses, these are likely to have additional dimensions that also escalate according to thermal bias, such as larger body sized species being regionally replaced by smaller, opportunist species 39 and being more likely to go extinct 29.
The influence of thermal bias on occupancy response was more pronounced in rhynchonellid brachiopods than bivalves, the latter having a lower mean vulnerability to extinction during rapid warming 29,39,40. Thermal performance has ecophysiological underpinnings 25,41,42, with some organisms having more specific or limited physiological adaptations. Evidence is mounting that different ecophysiological adaptations among taxa lead to different performance outcomes, including extinction risk 25, though quantitative comparisons of the thermal performance of brachiopods and bivalves are scarce 43. Our results therefore support the view that ecophysiology predisposes vulnerable taxa and traits to greater species immigration and extirpation at multiple scales, and their extinction risk is predictable via habitat loss 42,44,45. Groups vulnerable to warming may thus be more likely to show strong range shift responses, where habitat permits, such as bony fish 4,25, relying on escape rather than tolerance. Other vulnerable groups, such as reef corals, may be more restricted in their rate of habitat tracking (though see 22). Identification of vulnerable clades or traits, alongside spatial projections of their habitat loss from physiological principles and environmental factors 42, may aid understanding of the pressures regional warming will place on species 44–46.
Linking climate-driven range shifts to extinction risk
Species may avoid climate-induced extinction by colonising newly suitable habitat 13, hence marine fauna are expected to consistently trace their thermal preferences during climate change 21. Therefore, an escalating occupancy response in line with a species’ thermal bias may be a null expectation for a region under warming 10,15, leading species with particularly cool thermal biases to be vulnerable to local and global extinction. However, there are modern examples of how disequilibria between ambient temperatures and a population or assemblage can be stable rather than a dispersal failure, especially when observed at finer spatiotemporal scales than sampled here 15,28. Even at our scales, we observed an especially large variation in thermal bias for persisting species, suggesting either that finer scale thermal heterogeneity played a substantial role in permitting species to persist (i.e. refuges), or that many species were temperature generalists. Nevertheless, our finding that regional warming increased the slope of the relationships between thermal bias and response implies that greater magnitudes of warming on average increase the cost of disequilibria between species and climate. Furthermore, the overwhelming negativity of thermal biases across responses during warming phase 2, which coincided with the highest ratio of extirpations and extinctions to persisting species, may also have been amplified by the warming on-top-of warming climatic context, which increases extinction risk 47.
The largely overlapping thermal bias values for species being extirpated and going extinct, as observed here, may betray a cause of extinction. The mean thermal bias for extirpated species fell below the expectations of our linear model and instead fell within 95% confidence intervals for species going extinct. Meanwhile, observed thermal biases for immigrating and originating species aligned well with linear expectations. Either the thermal bias values for extirpations were unusually high, or the thermal bias values of the extinctions were unusually low. The first option assumes the true relationship was linear and thus the thermal bias expectations for extinctions were valid, but anomalous values for extirpated species alone are difficult to explain. The second option implies a non-linear true relationship, with the thermal bias values for extirpated species being valid but there being more extinctions observed than expected given their thermal bias. Given the anoxia of northern waters of the north-western Tethys, especially during the T-OAE (discussed below) 26, we suspect the anoxia overruled temperature in habitat suitability, leaving species going extinct with unusually low observed thermal bias values. Poor dispersal capabilities and/or dispersal barriers can lead to a species’ failure to lessen its population thermal biases by shifting distributions, thereby shrinking its geographical range 48, and making it vulnerable to global extinction 45,49. Mechanisms of and limitations to habitat tracking should be explored during other intervals with changes in climate, sea level, and geography e.g. 49.
Overriding effects of anoxic waters and terrestrial runoff
The well-sampled, oceanic-influenced regional clusters of north and east of Iberia best supported thermal determination of assemblage membership, where any deoxygenation prior to the T-OAE 50 apparently did not preclude a signal of thermal bias. Analyses of well-oxygenated environments such as outcrops from the south-west of Europe implicate Early Toarcian warming as the main regional driver of species loss, changes in bivalve-brachiopod assemblage structure, and their body size 39,40,51. During peak T-OAE (‘warming 2’ into ‘transitional’ phases), support for thermal determination of assemblage membership dwindled in the Germanic and British regions alongside the number of occurrences. Although aquatic deoxygenation can amplify the influence of warming on ectotherm performance 25,52, bottom water anoxia is likely to supersede the ecological influence of increased temperature. Several regions during the T-OAE are characterized by black shale deposition, where hypoxic and anoxic waters have long been associated with faunal turnover and extinctions 26. Accordingly, benthic macrofaunal recovery only began after seafloor ventilation resumed, and remained incomplete in the British region by the end of our study 53. During the T-OAE, the northern waters may have essentially been unavailable as habitat for species tracing their thermal niche. This may exemplify how species ranges can be compressed as they trace thermal preferences. Although fully marine (see Table S8), the more restricted northern waters likely had greater terrestrial influence, such that bottom-water anoxia was probably dependent on productivity, as nutrients were delivered from warming-enhanced weathering 54, rather than simply temperature-dependent deoxygenation. The HadCM3 model estimated slightly lower salinity in the Germanic and British basins also 55, ranging between 33.3 and 34.6ppt across scenarios, than the other regions, while salinity was always highest east of Iberia, ranging between 34 and 35.6ppt (see SM). The semi-enclosed setting, especially of the Germanic and British basins, also likely increased the influence of local processes that global models are unlikely to capture, with the reality likely being warmer and more seasonal than estimated by our models 56. Alongside changes in sea-level-dependent seafloor ventilation 53, water density differences from freshwater input may have also encouraged stratification 55,57. Enhanced capacity to model biochemical processes and extract variables, such as oxygen levels, should expand the ability of predictive models to account for additional niche requirements. While modern oxygen minimum zones continue to spread 24, our results show how regional-scale processes can complicate the predictability of assemblage responses to temperature change.
Besides temperature and salinity, other habitat requirements for a benthic species include suitable water depth and substrate conditions, which also dictate the conditions under which a species can be sampled. The northern regions were the only ones dominated by siliciclastic substrates, which could have blocked the immigration (alongside anoxia, see previous paragraph) of carbonate-affinity species. The largest and most consistent non-temperature change occurred at the Spinatum-Tenuicostatum transition, when substantial sea-level rise33 led to increases in the frequency of deep habitat occurrences from 20–50% to 90–96% regional share. However, species thermal bias remained highly significantly associated with its occupancy response through different statistical treatments to explore the importance of this spatiotemporal scenario (see Tables S4 and S5). Being 100s km across, our regions tended to cover substantial substrate and depth variation, such that finer scale analyses may be needed to detect the influence of non-temperature habitat variables. Our focus on two-timer species also emphasised longer-term changes of the more common and better-preserved species, of which our analyses support temperature change being a key driver at broad spatial scales.
Temporal and spatial scaling
Temporal and spatial resolutions in our study were ~ 1 million years and ~ 2000 km respectively, which need appreciation to compare our results with other studies. Finer scale variations were averaged out, such as the warming at the Pliensbachian–Toarcian stage-boundary 58, although permanent ecological changes such as extinctions from short-term pulses remain. Despite the myriad of factors influencing a species occurrence at fine spatial scales, climate is expected to be one of the dominating factors at broader scales 15,59. Significant effects of thermal bias have been assessed for modern assemblages at spatial scales from surveyed sites 11 to biogeographic ‘ecoregions’, more similarly sized to our regions 28,59. At intermediate spatial scales, Flanagan et al. 31 found larger thermal biases of fish assemblages over decadal scales than inter-annual scales, which might encourage expectations that marine communities rapidly maintain equilibrium with temperature, despite evidence often to the contrary 31. At much longer time scales and with spatially coarse temperature estimates, our data also supported equilibrium between assemblage mean thermal optima (CTI) and environment temperatures (Fig. S5). However, geographical context affected observations of thermal equilibrium in a study of planktonic foraminifera over thousands of years: mid latitude assemblages tracked climate change by turnover, but decreasing assemblage turnover at high latitudes under warming and low latitudes under cooling accumulated assemblage thermal bias 23. Regions of high climate velocity, such as the tropics and poles, are likely to demand faster species’ niche-tracking than lower climate velocity regions, which is more likely to push populations of multiple species nearer to their thermal niche edges 45. However, increasing thermal bias may only increase extirpations and extinctions when changes exceed species recent climatic experience 47. Temporal resolution is not a problem per se for the application of paleontological insights to modern issues 17, but limits the mechanisms for which we can observe evidence. Future work should be directed to understanding the mechanisms underlying observed palaeontological patterns and the transferability of those mechanisms to modern climate change and the current biodiversity crisis 17.
Based on our results for the northwestern Tethys, we expect regional species extirpations and especially immigrations to be already considerable (~ 5% and ~ 26%, respectively) at + 3°C warming, such as forecast under high emissions scenarios (RCP 8.5) by the end of the century for the North Sea 35. These extirpation and immigration values are similar to projections of a paleo-validated biodiversity model for the shelf seas of Europe by 2100 60. Although an application of our results to modern warming ignores the very different time scales (= observed rates of change), the loss of a species’ thermally suitable habitat can respond directly to the magnitude of warming, regardless of the rate of warming, such as supported by empirical patterns of high latitude extinctions during hyperthermal events 45. Rates of ancient climate changes may have been sufficiently slow for most species to track habitat availability but the extremely rapid anthropogenic rates of change are likely to divide response severity between species with greater and lesser dispersal abilities 45. This may be especially the case in the tropics where climate velocities are highest 61, leaving paleobiological extrapolations most likely as underestimates.
Rare species, both range-restricted or locally uncommon, are unlikely to make it into the fossil record and thereby into our analysis. If rare species are at higher extinction and extirpation risk or tend to have narrower thermal tolerances, the overall magnitude of assemblage change including rare species can be expected to be higher than we predict. Again, this implies that inferences based on paleobiology will tend to give underestimates of whole community responses.