The explanation of patterns in global biodiversity and taxonomic rates in the fossil record is one of the major challenges in evolutionary biology. Previous research has focused mainly on identifying (i) main trends and their consistency with alternative models of cladogenesis (Levinton 1979; Sepkoski 1984; Benton 1995, 1997; Courtillot & Gaudemer 1996) and (ii) potential causes of periods of major biodiversity turnover known as mass extinction events (Raup & Sepkoski 1982; Crowley & North 1988; Erwin 1990; Wignall & Twitchett 1996; Pope et al. 1998; Sheehan 2001; Joachimski & Buggisch 2002; White 2002; Benton & Twitchett 2003; Wignall 2005; Twitchett 2006). Only rarely have studies attempted to find statistical associations between the major global patterns and environmental variables throughout the Phanerozoic (Rothman 2001; Cornette et al. 2002).

Data on the distribution of families across geological strata (Benton 1993) represent the most recent compilation of global diversity covering both land and sea. The data were compiled under two different assumptions about the timing of fossil appearances for those families where dating was not to stage level (e.g. epoch); the ‘maximum’ assumption assumes that the lifespan of each family was from the first stage of the interval of first appearance until last stage of the interval of last appearance. The ‘minimum’ assumption is that the lifespan of a family was from the last stage of the interval of first appearance until the first stage of the interval of last appearance (Benton 1993, 1995). Data on marine genera of animals and protists (Sepkoski 2002) were compiled at stage resolution to provide maximum taxonomic inclusivity and greatest stratum comparability with Benton (1993). Three of the four major datasets we have used (the temperature, CO₂ and global family datasets) apply the same time scale (Harland et al. 1990). Therefore to maximize comparability among the different datasets and to reduce total rescaling error, because strata used in different scales are not always exactly equivalent, we rescaled the more recent genus level marine dataset to the same time scale. Taxonomic rates were recalculated using the new assumed time scale.

We included estimated atmospheric CO₂ concentrations (Berner & Kothavala 2001) as an explanatory variable in some analyses because previous analyses have shown that this may correlate with both taxonomic richness and rates (Rothman 2001; Cornette et al. 2002), and because CO₂ levels and temperature are often good predictors of each other (Royer et al. 2004). Estimates of CO₂ levels over the Phanerozoic in recent years have often been discordant (Berner & Kothavala 2001; Rothman 2002). We used those from the GEOCARB III model primarily because those data were assumed in deriving the temperature estimates we used (Royer et al. 2004). The data have also been used in previous studies on macroevolution (Cornette et al. 2002), show good convergence with other models that track carbon exchange, match well with other proxies (Royer et al. 2004) and correlate well with temperature (Royer et al. 2004), so are likely to provide a rigorous test of the effects of temperature itself.

We first used traditional measures of taxonomic diversity and rates as our response variables; total richness and the per-taxon rate (Foote 2000a). However, such measures may be susceptible to a number of biases that may obscure the true biological signal (Foote 2000a). We therefore used taxonomic rate measures that are independent of singletons: the standing richness of boundary crossing taxa (Foote 2000b), as well as the estimated per capita origination rate ‘p’ and extinction rate ‘q’ (Foote 2000a). Such measures are potentially more robust to variation in interval length and preservation rate, at least as long as this variation is small. Second, we included several control variables in our analysis to represent variation in preservation rate: numbers of sedimentary formations at epoch level were taken from columns 3 and 4 of appendix 2 of Peters & Foote (2001). Stage level data for marine formations (same original source) were obtained from Peters & Foote (2002). Data from the Correlation of Stratigraphic Units of North America used in Peters (2005) were obtained courtesy of the author. Finally, data on estimated marine genus preservation frequency in Foote (2005; from the pulsed/pulsed model) were obtained courtesy of the author. All these datasets were rescaled, as above for the Sepkoski data. The temperature and CO₂ data consist of estimates at 10Myr intervals and to provide richness and rate data at the same intervals we used the diversity and taxonomic rates of the geological interval into which each 10Myr point fell.

Time series (see electronic supplementary material for dataset and explanatory text) were transformed (log or square root), as appropriate, and detrended to remove long-term patterns (see table 1 in the electronic supplementary material). The detrending methodology was to fit linear and cubic polynomials, and splines increasing in flexibility to follow the patterns in the data (i.e. of increasing degrees of freedom, d.f.; Fewster et al. 2000; see figure 1 and table 1 in the electronic supplementary material). The specific detrender used was chosen to: make the series stationary; reduce the influence of outliers; make the distribution symmetrical; and maximize the signal-to-noise ratio. Low d.f. smoothers fit the trend and as the d.f. increase the smoothers pick up more of the shorter-term patterns in the data. Comparison between the fits, using autocorrelation or spectral analysis (see figure 2 in the electronic supplementary material), typically shows a pattern dominated by the trend (i.e. autocorrelation decays away very slowly), whereas very high d.f. smoothers take out all but noise around the pattern (so there are no patterns in the autocorrelation function). The appropriate smoother is one that maximizes the underlying signal. In the majority of cases, the detrender was a 5 d.f. spline (see figure 2 and table 2 in the electronic supplementary material). Detrended time series were then mean standardized.

In contrast, the per-taxon origination rates were all positively associated with temperature (table 1; figure 2), attaining significance for all families, marine families and marine genera. Again, some correlations were significantly stronger following a 10Myr lag. Per-taxon extinction rates were also significantly positively correlated with temperature for all groups (table 1; figure 3), and in no case were lagged correlations stronger than unlagged. Thus, high rates of extinction occurred at the same time as, or preceded, periods of reduced overall taxonomic diversity and increased origination of taxonomic lineages (Sepkoski 1998; Foote 2000b, 2005).

When we used measures of standing diversity, origination and extinction rate that are independent of singletons (Foote 2000a), the strength and significance of correlations were generally unaffected (table 1: boundary crossers, p, q). However, marine genus diversity became significant (figure 1), suggesting that the significant results based on traditional measures are robust to some potential biases, but also that some signal is obscured by the incorporation of singletons.

In general, associations remain in the same direction throughout the time series, though during some time intervals they differ or are weak. For example, in the Mid–Late Jurassic (ca 180–150Ma), both temperature and diversity residuals rise, giving rise to a positive association (figure 1). Extinction and temperature residuals are rather poorly associated during the Cambrian and Ordovician (ca 520–440Ma), and extinction residuals are sometimes positive during the icehouse Carboniferous phase (ca 350–280Ma; figure 3).

When atmospheric CO₂ concentrations were included as an explanatory variable in our analyses, temperature always remained significant, and CO₂ was normally not significant (table 1 in the electronic supplementary material). CO₂ was significant for both marine genus origination and extinction rate, and in the latter case was a stronger predictor than temperature. Overall, temperature was the better predictor of diversity and taxonomic rates.

In the current rapid transition from an icehouse to a greenhouse world (Huber et al. 2000; Lovejoy & Hannah 2005; IPCC 2007), Earth history may help us to estimate future effects on biodiversity (Culver & Rawson 2000). Prima facie, our results from the fossil record endorse those of ecological models (Thomas et al. 2004; Botkin et al. 2007), which demonstrate that expected future warming will adversely affect biodiversity. However, several qualifications are necessary.

A first qualification is that our results relate to the effects of residuals from the long-term trend. An increase in global temperature may therefore cause an increase in extinction rate but not necessarily an absolute decrease in biodiversity because the underlying trend is for biodiversity to increase over time.

A second qualification is that the coarse time scale of our data does not allow us to make short-term predictions, although short-term effects also cannot be excluded. Related to this, in predictive models relating extinction to climate change, it is the changes in climate rather than the future climates themselves that are generally held responsible (Thomas et al. 2004; Lovejoy & Hannah 2005; Botkin et al. 2007; Williams et al. 2007). Our results suggest that long-term average temperatures may exert a separate effect, independent of the rate of temperature change, and thus may be related to an entirely different mechanism (see below). Furthermore, the associations we have shown are only moderate in strength (correlations ≤0.5), and time intervals are apparent when associations can be reversed.

Finally, although we have shown an association between temperature and both biodiversity and taxonomic rates, this association may not be causative. Deducing causation from correlation is, of course, difficult. The lags shown in some of our analyses suggest that temperature is affecting biodiversity and evolutionary rates, but well-known links between organisms and geophysical processes suggest we should not yet rule out the opposite direction of causation (Rothman 2001). The periodic cycle of taxonomic richness and rates seen here (ca 140Myr, figures 1–3) is also potentially associated with cosmic ray flux, the age of meteorites and possibly sea-level changes (Rohde & Muller 2005). However, the age of meteorites is likely to be linked causally to cosmic ray flux (Rohde & Muller 2005), and the latter (Shaviv & Veizer 2003) and sea-level changes (Hallam & Wignall 1999) are likely to be linked to temperature. Clearly, many of these variables are likely to exert independent or correlated effects, and both our response and explanatory variables could also be responding independently to other variables not considered here.

Previous work at the scale of the Phanerozoic has suggested associations between atmospheric CO₂ concentrations and taxonomic richness and rates (Rothman 2001; Cornette et al. 2002). However, those studies were more taxonomically restricted, and solely used traditional richness and rate measures that are known to be susceptible to biases in the fossil record (Foote 2000a). Because it still remains significant in some of our analyses, our results suggest that CO₂ may still play some direct role on biodiversity and taxonomic rates independent of temperature, though for most of our analyses temperature is more important.

Despite the above provisos, our results demand that we speculate on causative links between temperature and both biodiversity and taxonomic rates. Previous work has suggested that extinction rates primarily drive changes in Phanerozoic diversity (Foote 2000b) and that origination rates rise following mass extinctions (see Sepkoski 1998). Thus it is plausible that the associations between temperature and both biodiversity and origination are primarily driven by changes in extinction rates with temperature. The lagged associations between temperature and origination rates, but not extinction rates, are consistent with this interpretation.

The five largest positive residuals for extinction in our time series correspond to previously identified mass extinctions (Raup & Sepkoski 1982; Benton 1995) in the end Ordovician, Late Devonian, end Permian, Early Triassic (a continuance of the former) and end Cretaceous (figure 3). Of these, four correspond well with peaks in global temperature in the time series (figure 3). Thus, our analyses suggest that these mass extinctions can be viewed as part of a wider trend seen across the fossil record as a whole. A number of recent studies have suggested a link between periods of rapid global warming, perhaps driven by large igneous province eruptions, and marine crises that often accompany extinction events (Wignall 2005; Twitchett 2006). It is plausible that there is an interaction between the long-term average global temperature and the onset of such marine crises, such that the effects of rapid global warming on marine systems are exacerbated in a world that is already in a greenhouse state, increasing the probability of oceanic anoxia or breakdowns in oceanic circulation. However, the associations shown in this paper apply to terrestrial as well as marine systems, and it is currently unclear how such marine crises might affect terrestrial systems (Wignall 2005).

The risk of future extinction through rapid global warming is primarily expected to occur through mismatches between the climates to which organisms are adapted in their current range and the future distributions of those climates (Thomas et al. 2004; Botkin et al. 2007; Williams et al. 2007), a mechanism that has doubtless also been relevant in the past (Twitchett 2006). Because taxa adapted to more tropical environments may be more vulnerable to short-term changes in climate in this way (Stanley 1986; Joachimski & Buggisch 2002; Williams et al. 2007), it is plausible that high long-term global temperatures may increase the general vulnerability of species to rapid climate change, and this may also explain the associations we find.