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

In the northern Adriatic Sea and in most semienclosed coastal regions worldwide, hypoxia induced by eutrophication in the late 20th century caused major die-offs of coastal marine organisms. However, ecosystem responses to hypoxia over longer centennial scales are unclear because the duration of direct observations is limited to a few decades and/or the temporal resolution of sedimentary archives is compromised by slow sedimentation and bioturbation. To assess whether perturbations of ecosystems by hypoxia recurred over centuries in the northern Adriatic Sea, we evaluate the timing and forcing of past hypoxia events based on the production history of the opportunistic, hypoxia-tolerant bivalve Corbula gibba, using 210Pb data, radiocarbon dating, amino acid racemization, and distribution of foraminifers in sediment cores that capture the past 500 yr in the Gulf of Trieste. Unmixing the stratigraphic record on the basis of 311 shells of C. gibba, we show that the reconstructed fluctuations in abundance do not correlate with abundances in the raw stratigraphic record. We find that production of C. gibba has undergone major decadal-scale fluctuations since the 18th century, with outbreaks corresponding to density of more than 1000 individuals per square meter. These outbreaks represent long-term phenomena in the northern Adriatic ecosystem rather than novel states characteristic of the 20th century eutrophication. They positively correlate with centennial-scale fluctuations in sea-surface temperature, indicating that the hypoxia events were coupled with water-column stratification rather than with nutrient enrichment. INTRODUCTION 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 20th 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). Shortterm 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 210Pb of the sediment and radiocarbon-calibrated amino acid racemization (AAR) of C. gibba. METHODS Species abundances, AAR, and activities of 210Pb and 226Ra 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 Repository1). 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 14C 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). RESULTS The 210Pb 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 1 GSA Data Repository item 2017104, details on methods, and raw data on shell ages, is available online at http://www.geosociety.org /datarepository /2017/ or on request from editing@geosociety.org. *E-mails: geoltoma@savba.sk; ivo.gallmetzer@ univie.ac.at; alex.haselmair@univie.ac.at; Darrell .Kaufman@nau.edu; martin.zuschin@univie.ac.at. GEOLOGY, April 2017; v. 45; no. 4; p. 363–366 | Data Repository item 2017104 | doi:10.1130/G38676.1 | Published online 7 February 2017 © 2017 The Authors. Gold Open Access: This paper is published under the terms of the CC-BY license. Downloaded from https://pubs.geoscienceworld.org/gsa/geology/article-pdf/3550906/363.pdf by guest on 06 April 2019 364 www.gsapubs.org | Volume 45 | Number 4 | GEOLOGY equal to 0.2 (based on median shell age) to 0.24 (based on 210Pb) 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 210Pb 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-cmthick 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 hypoxiatolerant 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 210Pb 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. Figure 1. Sediment core M28, collected in A.D. 2013 at 11 m water depth in the Bay of Panzano (Gulf of Trieste). A: Radiocarbon-calibrated amino acid racemization (AAR) ages of Corbula gibba shells. Boxplots show median ages and the 25th and 75th percentiles (whiskers denote minima and maxima). One shell at 112.5 cm is 816 yr old, and 2 shells at 147.5 cm are 689 and 1637 yr old. B: Downcore changes in 210Pb and ages based on the constant flux–constant sedimentation model. C: The correlation between the 210Pb chronology and AAR chronology based on the median age of C. gibba observed at a given depth below the sediment surface. The stippled 1:1 line represents perfect agreement between the two age estimates. Figure 2. A–C: Stratigraphic trends in absolute and proportional abundance of Corbula gibba and in absolute abundance of foraminifers differing in sensitivities to hypoxia. D: Downcore changes in time averaging (interquartile age range) observed in C. gibba death assemblages in 4and 5-cm-thick increments. Black points show raw values and gray points show values corrected for calibration errors (with 95% confidence intervals). Downloaded from https://pubs.geoscienceworld.org/gsa/geology/article-pdf/3550906/363.pdf by guest on 06 April 2019 GEOLOGY | Volume 45 | Number 4 | www.gsapubs.org 365 In contrast, taking into account mixing and burial, the record reconstructed on the basis of estimated ages of C. gibba shells shows decadalscale fluctuations peaking at 1980, 1890, and a double peak between 1810 and 1780 (Fig. 3C), rather than in the late


INTRODUCTION
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). Shortterm 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 210 Pb of the sediment and radiocarbon-calibrated amino acid racemization (AAR) of C. gibba.

METHODS
Species abundances, AAR, and activities of 210 Pb and 226 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 1 ). 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 14 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).

RESULTS
The 210 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 210 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 210 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-cmthick 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 hypoxiatolerant 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 210 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 decadalscale 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 2 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 r = 0.009, p = 0.96) or proportional abundances (Spearman r = 0.14, p = 0.46), indicating that raw stratigraphy does not reliably estimate the timing of C. gibba outbreaks.

Detecting Hypoxia in the Stratigraphic Record
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.

Forcing of Hypoxia
Eutrophication following from a major increase in agricultural and urban pollution (Justić et al., 1987), coupled with seasonal watercolumn 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 seasurface temperature reconstruction from the Gulf of Taranto (Versteegh et al., 2007) shows a longterm increase punctuated by three temperature maxima (A.D. 1780, 1880, and 1950Fig. 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.

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