Edited by David Jablonski, The University of Chicago, Chicago, IL, and approved July 30, 2012 (received for review January 18, 2012)
The geographic distribution of life on Earth supports a general pattern of increase in biodiversity with increasing temperature. However, some previous analyses of the 540-million-year Phanerozoic fossil record found a contrary relationship, with paleodiversity declining when the planet warms. These contradictory findings are hard to reconcile theoretically. We analyze marine invertebrate biodiversity patterns for the Phanerozoic Eon while controlling for sampling effort. This control appears to reverse the temporal association between temperature and biodiversity, such that taxonomic richness increases, not decreases, with temperature. Increasing temperatures also predict extinction and origination rates, alongside other abiotic and biotic predictor variables. These results undermine previous reports of a negative biodiversity-temperature relationship through time, which we attribute to paleontological sampling biases. Our findings suggest a convergence of global scale macroevolutionary and macroecological patterns for the biodiversity-temperature relationship.
Beyond small geographical scales, biodiversity consistently decreases with latitude (1⇓–3), reflecting a strong, probably causal, association between warmer climates and standing richness in both the terrestrial, water availability permitting (4), and marine realms (5, 6). Our understanding of the association between biodiversity and warm climates in space contrasts strongly with our models of how climate explains global diversity through time (7). One analysis of compendia of fossil taxa suggests that biodiversity declines with increasing global temperatures (8), but the focus on temperature as the driver, without reference to other variables, has drawn criticism (7). Analyses at smaller scales (geographic, temporal, or taxonomic) are equivocal (9). How can it be that warm temperatures should apparently have negative effects on biodiversity through time while also having positive effects across space (10) (i.e., contrasting macroevolutionary and macroecological patterns)? Recently, however, our understanding of past changes in biodiversity has been transformed by the application of techniques to control for sampling bias in paleodiversity data (11). Here we apply more-robust measures of fossil diversity, origination, and extinction through time to reevaluate the role of temperature in the context of other potential environmental drivers.
Much effort to understand macroevolutionary changes through the Phanerozoic has focused on marine invertebrates, first through Sepkoski’s genus-level compendium (12) and latterly via the Paleobiology Database Project (PaleoDB) (11). Putative factors proposed to drive temporal fluctuation in biodiversity include biotic drivers, such as competition between taxa (13, 14) and predation intensity (15). Alternatively, abiotic variables such as sea level change (16, 17), nutrient inputs and shelf redox conditions (17, 18), plate tectonic events (19), volcanism (20), bolide impacts (21), and global climate (8) have been invoked. The competing paradigms are labeled the “biotic” Red Queen and the “abiotic” Court Jester (22). However, as with the analogous debate in population biology, both paradigms probably have an explanatory role in macroevolution (23).
Current anthropogenic climate change, and the need to quantify its effects, has brought to the fore the long-recognized role of global climate in driving change in taxonomic richness (9). Over restricted geographic, temporal, or taxonomic scales, a variety of associations between temperature and richness have been reported (e.g., 23⇓⇓⇓–27). One analysis in deep time, using global-scale fossil compendia (12, 28), has found that high temperatures are associated with low taxonomic richness but high origination and extinction rates (8). This mirrors associations between atmospheric CO2 concentrations and macroevolutionary rates (29, 30), and a c.140 Myr cycle in Sepkoski’s compendium (12), alongside a similar periodicity of global climate modes (31). Despite being widely reported (e.g., refs. 32⇓–34), it is not easy to reconcile why warm temperature should apparently have negative effects on biodiversity through time while also having positive effects across space (i.e., contrasting macrovolutionary and macroecological patterns), and several counterexamples are known (7, 26, 35).
This paradox, a negative correlation between temperature and richness through time, rests on the assumption that range-through compendia adequately characterize fluctuations in biodiversity through time. Doubts have been raised by extensive recent work that shows strong associations between the amount of sedimentary rock and paleodiversity in fossil compendia (8, 36⇓–38). In short, variation in paleodiversity may be influenced by variation in the amount of preserved rock, how intensively sampled those rocks are, and the rate of publication of fossil lists from different time intervals. Most recent studies of the paleodiversity now attempt to account for bias (39), and the adoption of standardized subsampling techniques has been a major innovation (11).
Here we use time series analysis to test the association between global temperature and marine invertebrate macroevolution over the Phanerozoic using sample-standardized data. We first address whether, after removing long-term patterns, estimates of shallow-sea temperature and atmospheric CO2 concentrations remain robust predictors of these standardized macroevolutionary measures as they are for unstandardized measures. Second, we address whether a temperature proxy, δ18O, remains a significant predictor in analyses using a broad suite of potential biotic and abiotic explanatory variables. If it does, we can more strongly infer a causal role for temperature rather than simply being correlated with other variables.
After detrending all variables to remove long-term patterns (see Methods and Table S1), genus-level richness of marine invertebrates (using two commonly applied standardized subsampling techniques) is positively correlated with seawater temperature [Fig. 1D; item quota subsampling (IQS): r = +0.482, 95%CI +0.232 to +0.666; shareholder quorum subsampling (SQS): r = +0.289, 95%CI +0.010 to +0.519]. Richness is also positively correlated with atmospheric CO2 concentrations for IQS (r = +0.422, 95%CI +0.142 to +0.617), but not SQS (r = +0.192, 95%CI -0.044 to +0.422), while for both subsampling methods only temperature remains significant in multivariate models (IQS: b = +0.482, 95%CI +0.215 to +0.774; SQS: (b = +0.290, 95% CI +0.043 to +0.557). These positive relationships stand in contrast to the negative relationships found using richness measures from unstandardized data, even when rock record measures are used to control for sampling bias (8) (Fig. 1B and Table S2). Short-term fluctuations in the unstandardized and standardized richness series mostly coincide in time (Fig. 1 A and C) (17) but are of different magnitude, producing broader peaks and troughs that are out of phase, explaining why correlations with temperature are of opposite sign for the different measures.
Associations between tropical sea surface temperature and estimates of marine invertebrate richness, shown both through time series plots (A and C, N = 51) and through correlations between those time series (B and D), respectively. Time series show temperature (black circles) and standing richness (white and gray circles) for boundary crossers in Sepkoski’s compendium (A) and using item quota subsampling (C). Gray circles on all plots are the five mass extinctions of ref. 52. The residuals plotted are mean-standardized, after detrending following data transformation where necessary.
Community evenness is positively correlated with temperature through time (r = +0.451, 95%CI +0.171 to +0.654). Evenness and CO2 show a similar relationship (r = +0.332, 95%CI +0.100 to +0.532). Only temperature remains significant in multivariate models (b = +0.451, 95%CI +0.138 to +0.713).
Standardized origination rates (λ in ref. 13) are positively correlated with temperature (r = +0.378, 95%CI +0.137 to +0.598; Fig. 2B) and CO2 (r = +0.276, 95%CI +0.030 to +0.517), but only temperature remains significant in multivariate models (b = +0.375, 95%CI +0.112 to +0.615). Standardized extinction rates (μ in ref. 13) are positively correlated with temperature at a 10 Myr lag (r = +0.475, 95%CI +0.210 to +0.663), Fig. 2D) and CO2 (r = +0.404, 95%CI +0.152 to +0.618), and both remain significant in multivariate models as does their interaction (temperature b = +0.326, CO2 b = +0.176, interaction b = -0.205, 95% CI: -0.010 to -0.405). Similar relationships are found using unstandardized data (8) (Table S2).
Associations between tropical sea surface temperature and estimates of marine invertebrate origination rate (A and B) and extinction rate (C and D), shown both through time series plots (A and C, N = 51) and through correlations between those time series (B and D), respectively. Time series show temperature (black circles) and origination or extinction rates (white and gray circles), using the variables λ and μ from ref. 13. Gray circles on all plots are the five mass extinctions of ref. 52 (C) or rebounds from them (A). The extinction data in D are lagged by10 Myr. The residuals plotted are mean-standardized, after detrending following data transformation where necessary.
Modeled seawater temperature was generally a better predictor of the standardized macroevolutionary measures than δ18O (Table S3). To be conservative, δ18O was used in multivariate analysis. Significant predictors of the standardized macroevolutionary measures in bivariate correlations (Table S3) were generally quite similar and in the same direction as in multivariate analyses (Table 1). In bivariate correlations, high standing richness (IQS) is significantly predicted by low δ18O (i.e., high temperatures) and low extinction in the previous time step (Table S3). In multivariate models (Table 1), high δ13C (i.e., high productivity) is also associated with high diversity, and in lagged models high δ34S and high sea level are also associated with high diversity (Table 1). Similar results hold for standing richness using SQS, although notably high 87Sr/86Sr (i.e., high inorganic nutrient input) contributes to the unlagged models (Table 1). Mean evenness is also associated with low δ18O and low 87Sr/86Sr. Global standing richness and evenness measures are positively correlated (11) (Table S3). In bivariate analyses, origination rate is significantly correlated with extinction rate in the preceding interval, δ34S (i.e., high nutrient input), and 87Sr/86Sr. However, in multivariate analyses, δ18O is significant, in the same direction as for standing richness and for evenness (Table 1). In bivariate analyses, extinction rate is correlated with preceding-interval standing richness and evenness but also with low sea level (Table S3). In multivariate analyses, high δ34S was also significantly correlated with high extinction rates, but δ18O is not. However, when δ18O is replaced with seawater temperature in lagged models, temperature is retained (b = +0.334, 95%CI +0.0795 to +0.6574). Though the above results report experiment-wise significance, the strongest associations are robust to family-wise errors: e.g., the correlations between evenness and previous δ18O, and extinction and previous temperature (Table S3), remain significant even after strict Bonferroni correction across Table S3 (P < 0.000625).
Linear models of marine invertebrate macroevolution
The evidence for a relationship between temperature and richness in deep time, similar to that already demonstrated in relation to space, is equivocal. Some previous studies reported a positive relationship (26, 35, 40⇓–42), but other studies have found other associations (17, 24, 25, 43). At the largest scales, using range-through measures, the association between temperature and biodiversity had appeared to be negative (8), in apparent contradiction to broad spatial associations. Reconciliation is achieved when sample-standardized measures of the marine fossil record through time are used: The relationship between temperature and richness through time apparently becomes positive instead of negative, and hence temperature becomes related to diversity in similar ways in space and time.
Finding an association between temperature and standardized fossil subsamples was unsurprising given a 140 Myr periodicity in the IQS richness measure (44) matching the long-term cycle of climate modes (45, 46). While an overall association apparently remains, its sign has changed from negative to positive due to differences in the amplitude of shorter term peaks and troughs. (Fig. 1 A and C). This illustrates that, while range-through and sample-standardized curves can agree in many details (17), longer-term dynamics can differ (11), altering interpretations of causal mechanisms (13, 14, 47). Our finding that temperature is a statistically significant correlate of richness contrasts with (17). While the methods developed and deployed in ref. 17 are powerful tools for inferring causality in nonlinear systems (48), our findings are not exactly comparable: There are differences in the preprocessing of data, variables included, range of data in time and space, and in ability to consider more than three nonadditive variables in multivariate analyses. We note that δ34S approached significance in a number of our analyses (Table 1 and Table S3), along with other variables that could be driven by tectonics. A final issue relates to our ability to use results to predict the effects on biodiversity: The ecological significance of δ34S is more contentious (17, 18), whilst temperature has a number of clear, mechanistic links to the generation of biodiversity that are amenable to further testing.
Rock record metrics are increasingly used to control for geological biases (17, 39, 49). One problem with this approach is that it attempts to correct faults in the data. Ideally, we should use data without these faults. Furthermore, both the rock and fossil records could be driven by a third, common cause (17). Many possible rock record measures could be used to correct these data, which may not correlate well among themselves or which can be used in various ways (49, 50). The efficacy of these different methods remains undetermined. Our results lead us to urge caution about findings obtained through rock-record sampling correction, because conclusions will not necessarily be in close agreement with those using sample-standardized data.
It is tempting to claim that our results may help to predict the effects of current climate change (e.g., refs. 9, 32, 33). Our work suggests that warming of the oceans increases their potential to support biodiversity on geological time scales, but it would be misguided to conclude the same for short ecological time scales. During mass extinction intervals, the temperature-biodiversity relationship can break down (Fig. 1 C and D) and is influenced by compensating originations, which will not occur over short time scales, and it remains possible that high temperatures contribute to high extinction rates (Table S3). Caution is also appropriate because some of the correlations are weak, no consensus exists about the drivers of Phanerozoic diversity, and because weak associations are less robust to family-wise errors.
Possible drivers of the temperature-macroevolution associations may be partly inferred from the present results, which show that origination rates are also elevated with higher temperatures. Originations may be increased by temperature through numerous widely discussed ecological or microevolutionary mechanisms generally inferred from spatial associations (51) such as high rates of molecular evolution or increased habitat provision. Temperature changes may also sort lineages with different intrinsic turnover rates, as found for marine invertebrates during late Paleozoic glaciations (26, 40).
While extinction rates are also enhanced at high temperatures (Fig. 2D), causality is more equivocal; δ18O is not retained in models when other variables are controlled for statistically (Table 1), although modeled temperature estimates are retained in lagged models. This and previous work indicate that elevated richness is followed by increased extinction rates (13). This potentially complex network of effects may make the detection of independent drivers difficult but also predicts that increases in richness due to temperature are generally accompanied by taxonomic turnover (26, 40). Despite one or more mass extinction events (52) being plausibly precipitated by high temperatures (20, 53), some short-term temperature spikes were greater but are not matched by similar-sized extinction spikes (Fig. 2 C and D), while large reductions in richness are not always tied to major extinction events (Fig. 1C) (13).
Other abiotic predictors explain variation in macroevolutionary parameters in close accord with previous work (14, 16, 17, 54). These variables may partly explain intervals when the temperature-richness relationship appears more negative than expected (Fig. 1C); e.g., during the Great Ordovician Biodiversification Event (55), Late Jurassic/Early Cretaceous (12), and Neogene (12, 56) (Fig. 1C), intervals when richness increased in the context of a cooler climate, but when other environmental variables (e.g., sea-level, nutrient inputs, continental dispersion, shifting inter-tropical conversion zones) plausibly promoted richness increases (12, 55, 56). In general, our results suggest that warm, nutrient-rich, productive shelf seas at high eustatic sea levels are more conducive to the generation and maintenance of marine biodiversity over long geological time scales, although such settings are vulnerable to anoxia (57) and other abiotic stresses (58) that prevent a species-area relationship emerging. Multivariate analyses (Table 1) support mixed (Red Queen plus Court Jester) models for standing richness, evenness, and extinction, which contrasts with the notion (22) that Court Jester predominates at larger spatial and temporal scales but agrees with a recent analysis of Cenozoic marine Foraminifera (23).
In conclusion our results suggest that temperature is one of several abiotic and biotic variables that can enhance marine biodiversity, but also taxonomic turnover, across geological stages (26, 40). These findings raise the prospect of a greater integration between our understanding of macroevolutionary and macroecological patterns, heralding the possibility of a more mechanistic understanding of biodiversity patterns in deep time (59).
We assembled two datasets (details in SI Text): One to assess the effect of temperature and atmospheric CO2 concentrations on macroevolution, to compare with previous findings; the other to investigate more broadly the effects of environmental proxies and biotic variables on macroevolution.
In the first dataset (8) we used 10 Myr interval estimates of tropical shallow-seawater temperatures (45) and atmospheric CO2 concentrations (60). Marine invertebrate richness and evenness came from (figures 1 and 2 in ref. 11) using item quota subsampling (IQS) of fossil occurrences and also shareholder quorum subsampling (SQS) (14). Data on origination and extinction rates came from the measures λ and μ of ref. 13, which correct for pseudo-origination and extinction. For comparison, richness, origination and extinction from Sepkoski’s genus-level compendium were assembled using estimators (61), which control for preservation rates and interval duration. As a control for the availability of sedimentary rock through time, we took measures of the area of regional sedimentary rock record from Europe, Australia, and both combined (49). Data were standardized to the 2004 time scale (62). To ensure temporal matching in the sampling intervals, Akima interpolation splines (63) were applied, using the aspline function in R (64).
The second dataset built on ref. 18, containing marine invertebrate origination rates “λ” (13), eustatic sea level (65⇓–67), and isotopic proxies for environmental variables: δ18O (inverse proxy for temperature) (68); δ13C (proxy for biological activity) (68); 87Sr/86Sr (proxy for inorganic nutrient inputs) (69); δ34S (proxy for organic nutrient inputs or shelf redox conditions) (70). Data were added on marine invertebrate extinction rates (μ) (13), IQS standing richness and evenness data as above (11), and SQS richness (14) (see above), and modeled estimates of sea-water temperature (45). Akima interpolation splines (63) ensured temporal matching of sampling intervals.
Associations between variables were tested by Pearson correlation and linear modeling. Variables were transformed if necessary to allow parametric analysis and then detrended using smoothing splines to remove long-term patterns (Table S1). Residuals were standardized to a mean of zero and unit standard deviation (z-scores). Linear model simplification was performed by stepwise removal from a full model containing main effects (both datasets) and interactions (first dataset only), using the step function in R (64). Step removes parameters based on a comparison of AIC scores of all possible models with one less term. As the time series are serially autocorrelated, we tested significance (experiment-wise) through bootstrapping of the test statistic using the function boot in R, calculating confidence intervals on the test statistic using the bias corrected and accelerated (bca) technique (71).
We thank John Alroy and Dana Royer for providing data. A.J.M. was supported by a Royal Society of Edinburgh/Scottish Government Research Fellowship. This is Paleobiology Database publication 165.
Author contributions: P.J.M., M.A.B., T.G.B., and A.J.M. designed research; P.J.M., T.G.B., and A.J.M. performed research; P.J.M., M.A.B., T.G.B., and A.J.M. analyzed data; and P.J.M., M.A.B., T.G.B., and A.J.M. wrote the paper.
The authors declare no conflict of interest.
This article is a PNAS Direct Submission.
This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1200844109/-/DCSupplemental.
