Primary productivity controls on opportunistic bivalves during Early Jurassic oceanic deoxygenation

A key impact of increasing atmospheric pCO₂ 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 km² (∼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 (δ¹³C) 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 δ¹⁸O 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 δ¹³C 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 δ¹³C 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 (⁰). 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 Repository¹). 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, http://rsbweb.nih.gov/ij/index.html) 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, δ¹³C_org; Kemp et al., 2005; N isotopes, δ¹⁵N; Jenkyns et al., 2001; Mg/Ca and δ¹⁸O; 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).

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 cm²; Fig. 1A). In Interval 3, P. dubius abundance is less variable and changes more gradually.

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

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 (R² = 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).

The inverse correlation between B. radiata size and δ¹³C_org 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 δ¹³C_org. 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.

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 δ¹³C_org 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 mm²) 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 N₂-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 N₂-fixing cyanobacteria (blue-green algae), which provide fixed N to fuel primary productivity when the seawater SO₄²⁻:Mo ratio is high. Two lines of evidence indicate the likely presence of N₂ fixers during the Toarcian OAE and their importance in the trophic system. First, δ¹⁵N 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 N₂ 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 N₂ (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 N₂ fixers would have caused positive feedbacks by further depleting seawater [Mo].

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 SO₄^2–: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.