The fossil record has the potential to provide valuable insights into species response to past climate change if paleontological data are combined with appropriate proxies of environmental change. Here we use a novel, multivariate approach that combines a suite of geochemical proxies with high-resolution quantitative macroinvertebrate fossil data to study responses to early Toarcian warming (ca. 183 Ma). We show that benthic and nektonic ecosystems became decoupled during warming and were driven by different environmental variables. Benthic turnover was mostly driven by variations in seawater dissolved oxygen concentration, whereas turnover among the nektonic cephalopods was primarily controlled by variations in weathering, nutrient runoff, and primary productivity. Although rapid warming has been invoked as the main trigger of this event, the paleotemperature proxy was a poor predictor of marine community dynamics, and abiotic factors indirectly linked to temperature were more important. Given that similar environmental changes characterize other episodes of global warming and are impacting present-day marine communities, we predict that similar ecological responses occurred during other past events and are a probable outcome of current changes.

The fossil record is an archive of natural experimental data documenting the responses of Earth’s ecosystems to past environmental change, and may provide insights into future responses to current climate change (Jablonski, 2004). Past biotic crises may be associated with evidence for abrupt changes in sea level, temperature, atmospheric CO2, continental weathering and run-off, and ocean circulation and oxygenation, among others (Twitchett, 2006). Mass extinction events that occurred during episodes of climate warming record a similar suite of co-occurring, interlinked environmental changes (Jenkyns, 2010), and a key challenge is to determine which of these were most influential in controlling specific aspects of extinction or ecological collapse.

A substantial improvement in understanding the role of extrinsic environmental factors in controlling past changes in biodiversity resulted from quantitative correlation of multiple environmental proxies with biodiversity and turnover metrics (Cárdenas and Harries, 2010; Hannisdal and Peters, 2011). These studies, however, focused on long-term Phanerozoic-scale trends, and a quantitative, multi-proxy approach has hitherto never been applied to a single extinction event. Furthermore, the use of individual diversity or turnover metrics is not adequate for exploring the multifaceted changes in community structure that characterize past events, and a multivariate approach is needed (Belanger and Villarosa Garcia, 2014). This is particularly important because extinction events are selective and response varies between individual organisms, guilds, or clades (Foster and Twitchett, 2014).

Here we analyze quantitative paleoecological and paleoenvironmental data with ordination and multivariate multiple regression methods to provide the first high-resolution study of how changes in a suite of environmental factors controlled benthic and nektonic community dynamics through a past warming-related biotic crisis. We studied an ∼1.7 m.y. interval spanning the early Toarcian extinction event (Early Jurassic, ca. 183 Ma) that occurred during an episode of global warming and coincided with widespread oceanic anoxia (McElwain et al., 2005; Jenkyns, 2010). It impacted both benthos (Little and Benton, 1995) and nekton (Dera et al., 2010) and has been attributed to global environmental changes brought on by eruptions of the Karoo-Ferrar large igneous province, South Africa (McElwain et al., 2005).

The most comprehensive and highest-resolution paleontological and paleoenvironmental data sets through this event derive from studies of the Whitby Mudstone Formation (Cleveland Basin, North Yorkshire, UK) (Fig. 1; Fig. DR1 in the GSA Data Repository1). The rocks are well exposed and preserved, easy to correlate, have a rich fossil record, and have been subjected to many geochemical studies. In the early Toarcian, the Cleveland Basin was situated at a paleolatitude of 30°–40°N, and formed a sub-basin of the epicontinental Laurasian Seaway (Bjerrum et al., 2001).

Facies changes record a sea-level rise through the succession, which fines up from sandstone at the base of the formation through to silty mudstones and organic-rich laminated shale at the top, deposited under anoxic conditions (Hesselbo, 2008). A negative carbon isotope excursion (Kemp et al., 2005) coincides with a negative δ18O excursion indicating a local increase in seawater temperature (Bailey et al., 2003). Changes in 87Sr/86Sr reflect increased rates of continental weathering, and changes in δ34S, changes in δ98/95Mo, and increased total organic carbon (TOC) concentrations indicate higher nutrient input, basinal restriction, and seawater anoxia (Cohen et al., 2004; McArthur et al., 2008; Newton et al., 2011).

Published records were used to construct a data set of all available quantitative paleoenvironmental proxies spanning the study interval (Fig. 2; Table DR1 and Appendix DR2 in the Data Repository). Lithology was added as a categorical variable recording facies and sea-level changes. A high-resolution paleontological data set was assembled from the benthic macroinvertebrate quantitative abundance data of Danise et al. (2013), integrated with quantitative nektonic data from the same horizons and collected with the same methods (Appendix DR1; Fig. DR1). The total data set comprises 99 samples, 66 taxa, and 11,318 individuals identified to the species or genus level (Appendix DR2).

Diversity was measured using the Simpson diversity index (1-D). Diversity curves were smoothed with the three-point moving average function, and correlated with Spearman’s correlation coefficient (rs). Ordination methods were used to visualize patterns in multivariate data. Singletons were removed from the data sets; abundances were percentage and fourth-root transformed. Environmental proxy data were normalized. Principal coordinates analysis (PCoA) was performed using a Bray-Curtis similarity measure on the macrofossil data, and a Euclidean distance measure on the environmental data. Relationships between the macrofossil data sets and the environmental variables were analyzed using nonparametric multivariate regression. Individual environmental variables were analyzed separately for their relationship with the macrofossil data (marginal test), and all together were subjected to a step-wise selection procedure (sequential test) to develop a best-fit model using the sample-size corrected Aikake’s information criterion (AICc). Distance-based redundancy analysis (dbRDA) was applied to investigate relationships between the environmental proxies and fossil assemblages.

Environmental and Biotic Change Through Time

Pre-extinction benthic communities are more diverse than coeval nektonic ones, but both record similar, in-phase trends through the Grey Shale Member (Fig. 2; rs = 0.604; p = 0.0001). Following the crisis, the trends diverge: benthic diversity remains relatively low throughout the post-extinction Mulgrave Shale Member, whereas nektonic diversity is high (rs = –0.185; p = 0.1926). The change from parallel to divergent diversity trends implies a loss of coupling between benthic and pelagic communities. Peak nektonic diversity coincides with sea-level highstand (Hesselbo, 2008), restoration of marine phytoplankton populations (Bucefalo Palliani et al., 2002), and reduction in the global extent of anoxia (Pearce et al., 2008).

PCoA allows the visualization of patterns of environmental and biotic change through time. In the PCoA of environmental data, most of the pre-extinction tenuicostatum Zone samples cluster together (Fig. 3A), indicating relative environmental stability. Significant changes occur in the uppermost semicelatum Subzone (Fig. DR2) and in the lower exaratum Subzone, with a return to stable conditions (i.e., overlapping samples) in the falciferum Subzone (Fig. 3A; Fig. DR2). Ordination of the pre-extinction benthic data resembles that of the environmental data (Fig. 3B), indicating pre-extinction stability followed by turnover at the top of the semicelatum Subzone, where diverse benthic communities occupying a variety of infaunal and epifaunal trophic niches (Danise et al., 2013) are replaced by monospecific assemblages of the epifaunal suspension-feeding bivalve Bositra radiata (“extinction interval”, Fig. 2). The post-extinction benthos also has a stable taxonomic composition, but its ordination does not resemble the environmental PCoA (Fig. 3B). The bivalve Pseudomytiloides dubius dominates, with minor occurrences of B. buchii and Meleagrinella substriata (Fig. 2; Fig. DR2), all interpreted as epifaunal, suspension-feeding opportunists tolerant to low oxygen conditions (Caswell et al., 2009).

In contrast, the nekton records radical changes in faunal composition at several horizons, and its ordination resembles the environmental PCoA plot through the entire interval (Fig. 3C). This is consistent with the high evolutionary volatility of cephalopods and their sensitivity to environmental stress (Dera et al., 2010). During the crisis interval, belemnites undergo a shift in habitat, from cold bottom waters (Passaloteuthis) to warmer surface waters (Acrocoelites), possibly in response to bottom-water anoxia (Ullmann et al., 2014). As belemnites decrease in abundance, ammonites start to dominate, with the appearance of new genera of demersal Ammonitina and epipelagic Phylloceratina (Dera et al., 2010).

Relationship Between Biotic and Environmental Change

Each of the studied paleoenvironmental proxies has a significant relationship with the paleontological data, highlighting the complex, interlinked nature of the system (Table DR2). The AICc best predictive model for the benthos is, however, a combination of just δ98/95Mo, 87Sr/86Sr, and δ13Corg, which explains 54% of the variability in the fossil data (Table 1). Removing lithology from the model improves the selection criterion, indicating that variables other than facies change better explain changes in the benthos. The constrained dbRDA ordination shows that benthic turnover during the extinction interval correlates with increasing δ98/95Mo, a proxy for marine anoxia (Pearce et al., 2008), whereas changes in the post-extinction communities correlate best with variations in δ13Corg and 87Sr/86Sr, which reflect, respectively, changes in organic matter burial and productivity, and weathering rates (Fig. 3D; Table DR1). The most important variable is, however, δ98/95 Mo, which alone explains almost 40% of the total variation (Table 1), confirming that in the Jurassic, as in the present day (Diaz and Rosenberg, 2008), dissolved oxygen was a critical control on benthic marine ecosystems (e.g., Danise et al., 2013, and references therein).

In contrast, the best explanatory model for the nekton involves a combination of all paleoenvironmental variables except lithology, which explains 62% of the fossil data (Table 1). Variations in δ13Corg, TOC, and δ18O correlate with faunal changes during the extinction interval, whereas the transition from the exaratum to the falciferum Subzone assemblages is best predicted by changes in 87Sr/86Sr, δ34S, and δ98/95Mo (Fig. 3E). Among all the proxies, 87Sr/86Sr alone explains more than 37% of the total variation (Table 1). Increased 87Sr/86Sr during the early Toarcian has been interpreted as reflecting enhanced continental weathering driven by volcanic CO2 outgassing from eruption of the Karoo-Ferrar large igneous province (Cohen et al., 2004). Secular changes in the Sr cycle, related to nutrient input from continental weathering, have been shown to have a critical evolutionary impact on marine turnover during the Phanerozoic, with periods of increased nutrient availability associated with higher origination rates (Cárdenas and Harries, 2010). Thus, we infer that continental weathering and runoff were the key factors controlling nektonic communities during the early Toarcian, through their effects on nutrient supply and surface productivity. This is supported by the second most important controlling proxy, TOC (10.9%), which may reflect changes in primary productivity (Tribovillard et al., 1994).

Variations in 87Sr/86Sr are highly correlated with δ34S (r = 0.92; Table DR3) and could potentially mask the role of δ34S in the model. If 87Sr/86Sr data are excluded from the analysis, the resulting model still explains 58% of the variability of the data, with δ34S becoming the most important factor (Table DR4). Under anoxic conditions, positive δ34S excursions have been related to high rates of microbial sulfate reduction and pyrite formation (Gill et al., 2011), but marine sulfate levels also regulate the recycling of phosphorous through remineralization of organic matter by microbial sulfate reduction. If the δ34S signature at least partly reflects phosphorous recycling, which is vital for maintaining high primary productivity (Mort et al., 2007), this would provide additional support for the link between changes in the nekton and primary productivity.

The early Toarcian biotic crisis occurred during major global warming (McElwain et al., 2005), and temperature rise has been invoked as the main extinction trigger (Gómez and Goy, 2011). However δ18O, the proxy for seawater temperature during this interval (Bailey et al., 2003), is not one of the most important variables highlighted by our analyses. This may be further evidence that belemnite oxygen isotopes do not always reflect sea-surface temperatures, due to belemnite ecology or physiology, especially where there is an associated shift in habitat (Ullmann et al., 2014), or to the confounding effects of salinity on the δ18O signature (Price et al., 2013). Alternatively, if seawater warming occurred over a period of ∼0.7 m.y. as has been inferred by Bailey et al. (2003), then species may have been able to adapt through acclimatization.

During early Toarcian global warming, the benthos and nekton of the Cleveland Basin responded to different environmental drivers. Benthic change was mainly driven by hypoxia, whereas the nekton were driven by runoff, nutrient supply, and primary productivity. Given the global nature of early Toarcian warming, it is likely that the faunal response in other basins would have been similar. Other past episodes of global warming record a similar suite of environmental changes, and so we predict that they should also record a similar decoupling in response between the nekton and benthos. Furthermore, given that the Jurassic fauna is essentially “modern” in aspect, and that the key environmental controls identified here are already impacting shelf communities (e.g., Diaz and Rosenberg, 2008), we predict a similar response in modern marine ecosystems.

This project was supported by Natural Environment Research Council (NERC) grant to Twitchett (NE/I005641/1) and by a NERC Ph.D. studentship to Little. We thank W. Foster, D. Murphy, and M.-E. Clémence for help in the field, N. Thibault for discussions, and G. Lloyd for help with statistical analyses. M. Aberhan, A. McGowan, and an anonymous reviewer are thanked for their constructive insights during review.

1GSA Data Repository item 2015096, Figure DR1 (detailed stratigraphic log), Figure DR2 (detail of PCoA plots), Table DR1 (environmental variables and their interpretation), Table DR2 (marginal test), Table DR3 (Pearson Linear Correlation coefficients between environmental variables), Table DR4 (sequential test of the Nekton without 87Sr/86Sr), Appendix DR1 (methods, detailed information), and Appendix DR2 (macro-invertebrate and environmental data), is available online at, or on request from or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.