Severe paleoclimatic change during the Toarcian (Early Jurassic) oceanic anoxic event (OAE) was characterized by a negative δ13C excursion, increased weathering, higher seawater temperatures, oceanic deoxygenation, and mass extinction. We present abundance and size data (n 36,000) for the two dominant epifaunal bivalve species from the Toarcian OAE, Yorkshire, UK. We statistically correlate the biotic data with geochemical proxies for environmental change and show that our results are comparable with changes in present-day ecosystems affected by hypoxia. Bositra radiata dominated during declining oxygen levels immediately before the Toarcian OAE sensu stricto, and shell size doubled when δ13Corg was decreasing, indicating a connection with primary productivity. Small Pseudomytiloides dubius dominated during the Toarcian OAE and varied sharply in abundance, indicating that it was highly opportunistic. P. dubius shell size is strongly related to Mo concentration, [Mo]; this indicates a relationship between size and N2-fixing primary productivity via [Mo] limitation. A secondary factor contributing to small shell size was lower oxygen levels. After the Toarcian OAE diversity increased, P. dubius was less abundant and shell size doubled, indicating that bivalve populations were less limited by resources and conditions were more favorable. Size frequency distributions show that both Toarcian species had short life spans, rapid generation times, high recruitment, and high juvenile mortality. The opportunist Mulinia lateralis is a present-day analog for P. dubius. This research provides a case study for the long-term impacts of deoxygenation upon marine ecosystems, including that being observed today.

A key impact of increasing atmospheric pCO2 and eutrophication is the development of oceanic hypoxia (with oxygen content of 1%–30% of saturation, compared to anoxia, the absence of oxygen) (Diaz and Rosenberg, 2008). Hypoxic zones have increased exponentially over the past 50 yr and now cover >245,000 km2 (∼7% of ocean area), affecting more than 400 marine systems (Diaz and Rosenberg, 2008). Expanding hypoxia is predicted to reduce marine biodiversity and biomass, resulting in decreased production, altered animal behavior and distribution, modified trophic structure, and declines in fish stocks (e.g., Falkowski et al., 2011).

Quantitative high-resolution, long-term studies of ecological responses to oceanic deoxygenation from the geological record and present day are sparse (e.g., Falkowski et al., 2011). In present-day oceans, the impacts of hypoxia are confounded by multiple synergistic anthropogenic pressures, such as marine pollution and overfishing. Investigations of the biotic effects of oceanic deoxygenation from the geological record have the advantage that they are driven entirely by natural causes and not influenced by anthropogenic factors, and the results can show how marine biota might respond to increasing deoxygenation on millennial time scales. This study quantifies, at 1000 yr resolution, the relationships between geochemical proxies for environmental change during the Toarcian (Early Jurassic) oceanic anoxic event (OAE) (ca. 183 Ma) and the size and abundance of bivalves.

The Toarcian OAE is characterized by enhanced global deposition of organic matter (e.g., Jenkyns et al., 2001), a carbon isotope excursion (δ13C) of as much as ∼−7‰ within the ocean-atmosphere system (e.g., Kemp et al., 2005; Hesselbo et al., 2007; Al-Suwaidi et al., 2010), Mg/Ca and δ18O evidence for an increase in seawater temperature of as much as 13 °C (Bailey et al., 2003), and changes in the Sr and Os isotope ratios indicating increased rates of continental weathering (Cohen et al., 2004). Cyclostratigraphic analysis indicates (1) 75 cm periodic cycles, (2) one Milankovitch forcing parameter, and (3) that the δ13C excursion is composed of 4 Milankovitch paced shifts (shifts A–D in Fig. 1A; Kemp et al., 2005, 2011). The succession was divided into Intervals 1–4 (Fig. 1A) representing different seawater oxygenation levels based on changes in the Mo isotope composition of seawater (δ98/95Mo) and Mo concentration, [Mo] (Pearce et al., 2008). Interval 1 represents the onset of marine anoxia. Interval 2 represents a global decline in the extent of oxic conditions (OAE sensu stricto), and contains four pulses when reducing conditions were at a maximum that correlate with δ13C shifts A–D. Interval 3 represents a reduction in the global extent of reducing conditions, and Interval 4 represents a return to more oxic conditions regionally. The Toarcian OAE coincided with a marine mass extinction (Little, 1996; Caswell et al., 2009).

The Toarcian of northwest Europe was dominated by bivalves and ammonites (0). The aims of this study were to (1) quantify the size and abundance of the only two bivalve species [Bositra radiata (Goldfuss) and Pseudomytiloides dubius (Sowerby)] that occurred in significant numbers during the Toarcian OAE, (2) correlate these biotic data with geochemical proxies for environmental changes, and (3) compare the results with present-day benthic communities.

Abundance and size data were collected for P. dubius and B. radiata from ∼25 m of strata exposed near Whitby, Yorkshire, UK (from 54°32′48.64″N, 00°45′59.50″W to 54°24′23.37″N, 00°29′19.03″W) (see Fig. DR1 in the GSA Data Repository1). We measured ∼30,000 P. dubius and ∼6000 B. radiata shells from 221 stratigraphic levels (for each level, n = 24–954). Shell length and surface area were measured using image analysis (ImageJ 1.41i, U.S. National Institutes of Health, of field photographs. The mean stratigraphic resolution of the abundance and size data is 4.5–11 cm (Table DR1 in the Data Repository). The duration of the ∼75 cm Milankovitch cycle is uncertain, but is most likely precession or obliquity (Kemp et al., 2011), thus 1 cm represents ∼250–500 yr.

Shell size and abundance data were grouped by individual stratigraphic levels, three-point moving averages, and the ∼75 cm Milankovitch cycles. All discrete groupings of data were compared using nonparametric statistical tests (Kruskal-Wallis, Jonckheere-Terpstra, and Mann-Whitney U) and reduced major axis regression. Population size structure was explored using size frequency distributions (SFDs) and survivorship patterns were modeled using the Weibull frequency distribution. Relationships between shell size and geochemical proxies (δ98/95Mo, [Mo]; Pearce et al., 2008; TOC/P, δ13Corg; Kemp et al., 2005; N isotopes, δ15N; Jenkyns et al., 2001; Mg/Ca and δ18O; Bailey et al., 2003) were investigated using linear regression and forced-entry multiple regression analyses. Data were considered over the entire event, using geochemical trends, and over Intervals 1–3, reflecting different levels of oxygenation defined by Pearce et al. (2008).

Pseudomytiloides dubius and B. radiata occur as 1–2-shell-thick monospecific shell pavements. The shells are thin, unabraded, preserved as original calcite, and represent in situ time-averaged accumulations formed by natural mortality with minimal postmortem displacement; this concept is supported by well-preserved lamination and a lack of trace fossils. A few horizons in the lower part of bed 32 contain shell fragments and/or gutter casts suggesting minor reworking. Positively skewed SFDs with large size variation suggest that no size-selective taphonomic removal occurred (Fig. DR5).

Both B. radiata and P. dubius show weak negative relationships between shell size and abundance (Figs. 1B and 1C). SFDs for both species are positively skewed with unimodal peaks in the smallest two size classes (Fig. DR5). P. dubius SFDs are more strongly skewed than B. radiata SFDs. The dominant survivorship pattern for both species shows constantly increasing mortality with age (Table 1).

Abundance Variations

Pseudomytiloides dubius and B. radiata rarely co-occur, and when δ98/95Mo and [Mo] are at their lowest, fewer than four individuals occur (three “very poorly fossiliferous intervals”; Fig. 1A; Caswell et al., 2009). B. radiata abundance varies fourfold within Interval 1, is absent in Interval 2, and briefly reappears in Interval 3. In Interval 2, P. dubius shows abrupt and large changes in abundance (∼0–500 individuals per 400 cm2; Fig. 1A). In Interval 3, P. dubius abundance is less variable and changes more gradually.

Size Variations

Bositra radiata median shell size shows statistically significant increases (Kruskal-Wallis, Jonckheere-Terpstra, p < 0.001; Fig. 1A) as δ13Corg decreases both throughout Interval 1 and within each of the δ13Corg cycles in Interval 1 (p < 0.001; Table 1; Fig. 1A). When δ13Corg is decreasing, the B. radiata SFD is unimodal and less positively skewed (Fig. DR5) than when δ13Corg is increasing.

Pseudomytiloides dubius median shell size shows statistically significant increases both throughout the Toarcian OAE (Kruskal-Wallis, Jonckheere-Terpstra, p < 0.001) and from Interval 2 to Interval 3 (Table 1; Fig. 1A). SFDs are unimodal and less positively skewed over Interval 3 than Interval 2 (Fig. DR5). Of the SFDs in Interval 3, ∼50% are multimodal (Fig. DR6). Weibull shape parameters indicate that P. dubius survivorship was significantly more age dependent during Interval 2 than Interval 3 (Table 1).

Relationships with Geochemical Proxies

Results of linear regressions of [Mo], δ98/95Mo, and TOC/P with both mean and maximum shell size are statistically significant (Table DR4; note that the mean is more representative of the whole population; the maximum is not unduly influenced by any one size class). Linear regressions show that the strongest relationships are for data averaged per ∼75 cm Milankovitch cycle. Multiple regression of P. dubius size averaged per Milankovitch cycle with δ98/95Mo, [Mo], and TOC/P shows a strong relationship (R2 = 0.53). [Mo] was the strongest and only statistically significant predictor variable in the model (Table 2). Linear regression of the data split stratigraphically based on the geochemical proxies shows the following statistically significant relationships: a strong positive relationship between [Mo] and P. dubius maximum shell size for Intervals 2 and 3, with the rate of size increase being greater during Interval 2 than Interval 3 (Fig. 1D); and, for Intervals 2 and 3, a strong positive correlation between shell size when δ98/95Mo > 1.3 (Fig. 1E). There were no statistically significant (p > 0.05) correlations with any of the other geochemical proxies. The relationships between size and δ98/95Mo could result from an indirect relationship with a third variable such as temperature.

Natural ecological variability and time averaging of the rock record have resulted in high temporal fluctuations of the raw size data from the Toarcian OAE. Abrupt and large shifts in abundance indicate that both species were highly opportunistic and their population size could increase rapidly under favorable environmental conditions such as increased food supply. These observations are comparable to present-day benthic communities that experience decreased diversity, increased density, and dominance by opportunists during long-term or severe deoxygenation (Levin et al., 2009).

The inverse relationship between body size and abundance for the two Toarcian bivalve species (Figs. 1B and 1C) is common in present-day organisms (Griffiths, 1992), and shows that smaller individuals utilize resources (e.g., food and energy) more efficiently. The regression line gradients (b) show that resource utilization by the two bivalves is comparable to that of present-day marine benthos (b = −1.08 to −1.22; Griffiths, 1992). In addition, resource utilization in relation to body size is more efficient for P. dubius (b = −1.4) than B. radiata (b = −1.0). This increased efficiency of resource utilization by smaller individuals could explain why taxa of decreased body size often occur after mass extinctions (termed the Lilliput effect).

The strong positive skew of the SFDs indicates that recruitment and juvenile mortality were high, particularly for P. dubius. Present-day seasonal hypoxia results in communities dominated by taxa of small size, with short life spans and generation times, and high larval production (Levin et al., 2009). Both of the Toarcian bivalves probably had life spans and generation times similar to those of the present-day opportunist Mulinia lateralis (see section 3 of the Data Repository), which is adapted to anoxic conditions (Shumway et al., 1993) and lives for ∼2 yr, reaching sexual maturity after just 2 months (cf. Mya arenaria, an equilibrium species with a 20 yr lifespan and sexual maturity of 2 yr).

The lack of co-occurrence of P. dubius and B. radiata is interpreted as competitive exclusion driven by differences in their preferred food source and adaptation and/or tolerance to low oxygen levels (Caswell et al., 2009). P. dubius is the only abundant bivalve found during the most reducing conditions of the Toarcian OAE and other Mesozoic periods, indicating that it was better adapted to deoxygenation. Shallower water areas could have provided new recruits for the bivalve populations (Caswell and Coe, 2012).

Interval 1

The inverse correlation between B. radiata size and δ13Corg in Interval 1 (Fig. 1A) on both short and long time scales is attributed to greater availability of food during increased primary productivity that also lowered δ13Corg. This is supported by the less positively skewed SFDs found during periods of higher primary productivity, that show that B. radiata populations had lower recruitment and juvenile mortality, attained larger size, and had more constant mortality rates.

Interval 2

The almost complete dominance of positively skewed P. dubius SFDs in Interval 2 (Fig. 1A) shows that recruitment and juvenile mortality were high and that populations consisted predominantly of one cohort. This suggests that recruitment only occurred a few times and individuals had short life spans. The four pulses of elevated reducing conditions are interpreted to have led to the very poorly fossiliferous intervals because, prior to δ13Corg shifts C and D (Fig. 1A), P. dubius abundance decreases concurrently with δ98/95Mo.

The small overall size of P. dubius in Interval 2 (∼30% of shells have a surface area >60 mm2) would have allowed more efficient use of limited resources. Shumway et al. (1993) showed that present-day small individuals (<5 mm) of M. lateralis are twice as tolerant of anoxia as larger individuals (>10 mm) at 20 °C. A weak relationship between P. dubius size and paleoredox is suggested by the statistically significant positive correlation between bivalve size and both δ98/95Mo and TOC/P.

[Mo] is an essential cofactor for key enzymes involved in the nitrogen and sulfur cycles that can, in turn, have a strong influence on primary productivity (Cole et al., 1993). We suggest that [Mo] limited primary productivity by N2-fixing cyanobacteria and/or methanotrophs during the Toarcian OAE. This process produced a strong connection between [Mo] and bivalve size, and thus links primary production with secondary production by the suspension-feeding bivalves. Throughout Interval 2 the [Mo] in seawater was unusually low because the reducing conditions severely depleted the global seawater inventory (Pearce et al., 2008; Fig. 1A). Cole et al. (1993) showed that [Mo] can limit marine N2-fixing cyanobacteria (blue-green algae), which provide fixed N to fuel primary productivity when the seawater SO42−:Mo ratio is high. Two lines of evidence indicate the likely presence of N2 fixers during the Toarcian OAE and their importance in the trophic system. First, δ15N data from the interval of reducing conditions in the Toarcian OAE indicate that the water column was partially denitrified (Jenkyns et al., 2001), a condition that facilitates blooms of N2 fixers (e.g., Conley et al., 2009). Second, as oxygen concentration started to decrease prior to the Toarcian OAE, there was a change in the primary producers from dinoflagellates and calcareous nanoplankton to the prasinophyte (green algae) Tasmanites and sphaeromorphs (Palliani et al., 2002). Both Tasmanites and sphaeromorphs appear to have been well adapted to denitrified conditions, because Tasmanites have a strong preference for reduced nitrates (Prauss, 2007) and sphaeromorphs are thought to have close affinities to methanotrophs (Palliani et al., 2002), some of which can also fix N2 (Murrell and Dalton, 1983). The high P. dubius abundance coincides stratigraphically with that of Tasmanites and sphaeromorphs, and thus it seems likely that they represent the preferred food sources for P. dubius (see section 3 in the Data Repository). In addition to the low global inventory of [Mo], periodic blooms of N2 fixers would have caused positive feedbacks by further depleting seawater [Mo].

Interval 3

The overall increase in P. dubius size and decrease in abundance in Interval 3 is interpreted to indicate greater resource availability than in Interval 2. This is supported by greater species diversity (Caswell et al., 2009; oxygenated interval in Fig. 1A) and multimodal SFDs, suggesting that multiple recruitments occurred. The sixfold [Mo] increase between Intervals 2 and 3 would have lowered the SO42–:Mo and may explain why [Mo] had a less pronounced effect on P. dubius size (via primary production) during Interval 3 (Fig. 1D). B. radiata also briefly appears in Interval 3, suggesting improved conditions.

During the Toarcian OAE, there were large and geologically rapid changes in abundance, body size, and population dynamics of P. dubius and B. radiata. Periods of higher primary productivity and greater oxygenation supported larger bivalves. When oxygen levels were at their lowest, primary productivity was limited and bivalves were virtually excluded from this area of the shelf for tens of thousands of years. As food supply and oxygen levels increased, the bivalves that colonized the seafloor were highly opportunistic, had short lifespans, rapid generation times, and high juvenile recruitment. Their small body size enabled the most efficient use of resources.

This study shows that it is possible to establish statistically robust relationships between biological and geochemical proxies for environmental change at high temporal resolution during a period of ancient global oceanic deoxygenation. Quantifying the biotic changes that occurred during other periods of oceanic deoxygenation is an important goal for improving our understanding of how present-day marine systems will respond to deoxygenation on long time scales.

This research was funded by a National Environment Research Council Ph.D. studentship to Caswell. We thank E. Maddison, D. Weatherby, A. Hunt, P. Clay, and P. Sheldon for help with field work; A. Cohen and P. Skelton for help during the initial stages; and C. Frid, W. Gosling, C. Mahaffey, C. Pearce, D. Ruban, E. Thomas, and an anonymous reviewer for comments on the manuscript.

1GSA Data Repository item 2013324, methods, data analyses, statistical results, and additional background information, is available online at, or on request from or Documents Secretary, GSA, P.O. Box 9140, Boulder, CO 80301, USA.