Stratigraphic unmixing reveals repeated hypoxia events over the past 500 yr in the northern Adriatic Sea

Over recent decades marine coastal habitats globally have been affected by eutrophication and hypoxia, leading to the development of coastal dead zones (Gooday et al., 2009; Rabalais et al., 2010). The ecosystem of the northern Adriatic Sea was affected by eutrophication, algal blooms, mucilage blooms, and hypoxia during the second half of the 20^th century (Degobbis et al., 2000; Djakovac et al., 2015). However, archives of past hypoxia events in the northern Adriatic Sea are either limited to the past 150 yr (Puškarić et al., 1990; Barmawidjaja et al., 1995; Sangiorgi and Donders, 2004; Lotze et al., 2011) or are characterized by low stratigraphic resolution due to slow sedimentation and bioturbation (Van der Zwaan, 2000). Short-term archives of hypoxia can underestimate the frequency of older natural hypoxia events that intermittently occurred for centuries and millennia in marine ecosystems now affected by eutrophication (Osterman et al., 2008). To disentangle the contributions of eutrophication and climate-driven water-column stratification to past hypoxia and to determine whether past communities in the northern Adriatic coastal system responded to hypoxia the same way that the modern community does, it is necessary to trace the duration and frequency of such events over longer time scales.

Here we reconstruct a 500-yr-long record of hypoxia events in the northern Adriatic Sea by tracing changes in production of an opportunistic, hypoxia-tolerant bivalve species, Corbula gibba, and analyze molluscan and foraminiferal assemblages in two 1.5-m-long sediment cores (M28 and M29, 16 cm diameter) collected in A.D. 2013 at 11 m water depth in the Bay of Panzano (Gulf of Trieste; Figs. 1 and 2). C. gibba survives seasonal hypoxic events and, following mass mortality affecting most of the benthic macrofauna, it achieves density of thousands of individuals per square meter and can represent >80% of individuals of the bivalve assemblage (Hrs-Brenko, 2006; Nerlović et al., 2011). To account for mixing of different age cohorts, we analyze ²¹⁰Pb of the sediment and radiocarbon-calibrated amino acid racemization (AAR) of C. gibba.

Species abundances, AAR, and activities of ²¹⁰Pb and ²²⁶Ra by gamma spectroscopy were analyzed in 2-cm-thick intervals in the upper 20 cm and in 5-cm-thick intervals below (see the methods in the GSA Data Repository¹). The 2 cm intervals were pooled to 4 cm intervals in analyses of AAR and abundances (Fig. 2). AAR was used to estimate the postmortem age of 311 specimens of C. gibba in core M28 (Fig. 1A). The rate of AAR was calibrated with 11 shells dated by ¹⁴C and 3 live-collected individuals using Bayesian fitting (Allen et al., 2013; Dominguez et al., 2013). Age-frequency distributions (AFDs) are right skewed in the upper 4 cm and shift to normal-shaped AFDs in increments below (Fig. 1A). To account for mixing of the stratigraphic record (Fig. 3A) and to estimate the population density of C. gibba, we reconstruct the AFD of the entire core on the basis of shell ages in directly dated increments (white bars in Fig. 3B). First, we interpolate the shape of AFD of each undated increment by pooling AFDs of the directly dated increments stratigraphically located below and above. The mean and standard deviations of such pooled distribution are assumed to capture the AFD of the undated increment. Second, we resample shells to the total number of C. gibba in each increment in both cores, either from directly estimated or from interpolated AFDs, and tally the number of resampled shells in 10 yr cohorts. The resulting distribution of shell ages represents the AFD of all shells in the core (gray bars in Fig. 3B). To account for burial of shells (i.e., older shells can be buried below the bottom of the cored interval and thus under-represented in the entire sample of shells), we model burial with the exponential distribution and divide the whole-core AFD by its survival function (Tomašových et al., 2016), where mean time to burial (354 yr) is approximated from whole-core shell burial rate (0.42 cm/yr; Fig. 3C).

The ²¹⁰Pb ages correlate with median shell ages in the upper 30 cm (Pearson r = 0.99, p < 0.0001, Figs. 1B, 1C; Table DR1 in the Data Repository), with apparent sedimentation rates equal to 0.2 (based on median shell age) to 0.24 (based on ²¹⁰Pb) cm/yr during the deposition of the upper 30 cm (ca. A.D. 1890 to present). The upper 6 cm of surface sediment are homogenized over ∼20 yr on the basis of ²¹⁰Pb values and median ages of C. gibba. The median age of C. gibba increases from 15 yr (ca. 1998) in the uppermost 4 cm increment to ∼400 yr at 150 cm (ca. 1615; Fig. 1A), where the interquartile age range encompasses shells 500 yr old. The median date based on shell ages increases at a slower rate between 30 cm and 150 cm (from 1615 to 1864, i.e., 0.48 cm/yr) than above 30 cm (from 1864 to 2014, i.e., 0.2 cm/yr).

Absolute abundances of C. gibba in 4–5-cm-thick increments show low-amplitude gradual fluctuations through the cored interval (Fig. 2A), whereas the proportion of C. gibba relative to the number of all bivalves displays strong fluctuations, although it never exceeds 50% (Fig. 2B). Proportional abundances of hypoxia-tolerant foraminifers increase in the top 20 cm, although hypoxia-sensitive foraminifers do not show any major fluctuations through the cored interval (Fig. 2C). Proportional abundances of hypoxia-tolerant foraminifers (Bulimina and Uvigerina) correlate positively with abundances of C. gibba (Spearman coefficient r = 0.45, p = 0.01). Time averaging (interquartile age range) of C. gibba in 4–5-cm-thick increments varies between 8 and 30 yr in the upper 25 cm (as also observed on the basis of the ²¹⁰Pb profile) to ∼60–300 yr in the middle and lower part of the core (Fig. 2D). Changing the scale of the stratigraphic record from sediment depth to median calendar date, proportional and absolute abundances (Fig. 3A) attain maxima in the late 1600s, early 1800s, and late 1900s.

In contrast, taking into account mixing and burial, the record reconstructed on the basis of estimated ages of C. gibba shells shows decadal-scale fluctuations peaking at 1980, 1890, and a double peak between 1810 and 1780 (Fig. 3C), rather than in the late 1600s, early 1800s, and late 1900s, as suggested by stratigraphic records. The abundances reconstructed to 100–120 dead individuals in 10 yr cohorts and maximum lifespan equal to 5 yr in a core with cross-sectional area of 0.04 m² correspond to a standing density of 1250–1500 individuals per square meter, a density similar to times of C. gibba outbreaks observed today. Their outbreaks are thus not restricted to eutrophication-induced hypoxia events of the 20^th century, but date back to earlier centuries. The unmixed record of C. gibba does not correlate with stratigraphic records of absolute (Spearman ρ = 0.009, p = 0.96) or proportional abundances (Spearman ρ = 0.14, p = 0.46), indicating that raw stratigraphy does not reliably estimate the timing of C. gibba outbreaks.

A significantly positive correlation between abundances of hypoxia-tolerant foraminifers and C. gibba, the coincidence between the peak in abundance ca. A.D. 1980 and multiple hypoxic crises observed in the Gulf of Trieste in 1974, 1980, and 1983 (Faganeli et al., 1985), and the coincidence between the decline in abundances of C. gibba and low frequency of hypoxia in the past two decades (Mozetič et al., 2010; Giani et al., 2012), suggest that abundance peaks of C. gibba indicate past hypoxia events. Although summer hypoxia in the Gulf of Trieste occurs over a few days or few weeks, recoveries of soft-bottom communities from such events are slow and can last for several years or decades (Stachowitsch, 1991). Therefore, opportunistic species persist even after oxygen concentrations return to normal levels, and their remains carry long-term yearly or decadal memories of hypoxia. The duration of C. gibba outbreaks detected here probably reflects longer intervals in the aftermath of seasonal and/or clusters of multiple hypoxia events.

The stratigraphic record of C. gibba is compromised by centennial-scale time averaging of assemblages in subsurface increments, (1) leading to low to moderate proportional abundance of C. gibba (<50%) relative to dominance expected to occur at times of their outbreaks, and (2) obscuring timing of these events. The signature of past hypoxia events is diminished by bioturbational mixing with communities unaffected by hypoxia because (1) time averaging of 5 cm increments expected purely due to condensation at the observed net sedimentation rates (0.2 cm/yr in the upper part and 0.48 cm/yr in the lower part) would be just ∼10–25 yr, and (2) time averaging increases downcore (Fig. 2D) in spite of the concomitant increase in sedimentation rate.

Eutrophication following from a major increase in agricultural and urban pollution (Justić et al., 1987), coupled with seasonal water-column stratification due to particular circulation patterns, contributed to seasonal hypoxia in the northern Adriatic Sea in the late 20^th century (Djakovac et al., 2015). These changes triggered shifts in ecosystems toward dominance of fast-growing opportunistic species (Di Camillo and Cerrano, 2015). However, the outbreaks of C. gibba observed centuries ago imply that hypoxic events occurred before the 20^th century eutrophication; they were probably related to prolonged duration of seasonal stratification driven by changes in circulation rather than to anthropogenic nutrient enrichment. First, an alkenone-based sea-surface temperature reconstruction from the Gulf of Taranto (Versteegh et al., 2007) shows a long-term increase punctuated by three temperature maxima (A.D. 1780, 1880, and 1950; Fig. 3D) that coincide with abundance peaks of C. gibba. Second, mucilage blooms (gelatinous sheets and clouds of exopolymeric aggregates formed by phytoplankton exudation and degradation) tend to be stronger in warmer years (Danovaro et al., 2009) and coincide with the formation of the Istrian Coastal Countercurrent, which contributes to prolonged water-column stratification (Supić et al., 2000). The degradation of mucilage aggregates sinking to the bottom can enhance oxygen consumption and can contribute to hypoxia and benthic mass mortality (Degobbis, 1989; Sellner and Fonda-Umani, 1999).

Fluctuations in nutrient abundance driven by river discharge from the Isonzo River could also play some role (Cozzi et al., 2012), but cross-correlations between precipitation and C. gibba abundances are weak and insignificant, and C. gibba abundances correlate negatively with abundance of the foraminifer Valvulineria (Spearman r = −0.55, p = 0.001). This genus prefers organic-rich silty substrates affected by high sediment flux, and tends to be associated with colder but more humid episodes, coupled with higher water discharge from rivers in deltaic environments (Frezza et al., 2005). The peaks in abundance of C. gibba were thus probably enhanced by warming, increasing the duration of stratification, rather than by river-driven eutrophication (Crema et al., 1991).

C. gibba has been abundant on hypoxic muddy bottoms in the Paratethys and Mediterranean since the Miocene (Zuschin et al., 2014), and was abundant in the northern Adriatic Sea in the Pleistocene (Dominici, 2001; Scarponi and Kowalewski, 2013; Kowalewski et al., 2015). Semienclosed coastal ecosystems in the Paratethys were prone to water-column stratification and high riverine siliciclastic input that enhanced seasonal hypoxia and thus naturally triggered outbreaks of opportunists like C. gibba. Although recent outbreaks of C. gibba are intensified by human-induced hypoxia in the northern Adriatic Sea, they have a long history, are not novel to the 20^th century, and capture conditions analogous to those characterizing past Paratethys and Mediterranean ecosystems.

This study was funded by the Austrian Science Fund (FWF project P24901) and the Slovak Grant Agency (VEGA 0136–15). We thank T.D. Olszewski, J.W. Huntley, and an anonymous reviewer for comments, and J. Sedmak and F. Perco for help with sampling.