TRACING ORIGIN AND COLLAPSE OF HOLOCENE BENTHIC BASELINE COMMUNITIES IN THE NORTHERN ADRIATIC SEA

The northern Adriatic Sea is a young ecosystem that originated at the onset of the Holocene, approximately 11,000 years ago, when rising sea level inundated the vast alluvial plains that had been exposed since the Würmian glaciation (Brambati 1992; Amorosi et al. 2008). In contrast to most other shallow coastal seas worldwide, it was affected by a dense human colonization of shores and hinterlands from early on in mankind's history (Turney and Brown 2007). The anthropogenic alteration of Adriatic ecosystems probably reached an initial peak as early as the time of the Roman Empire, with local exhaustion of fish stocks and increased river sediment loads due to deforestation-induced erosion (Kaligaric et al. 2006; Lotze et al. 2011). Human pressure rose again in the Middle Ages due to population growth and the expansion of coastal settlements (Lotze et al. 2006). During the last 200 years, the human impact on continental shelves has increased exponentially with the many repercussions of industrialization and a globalized economy: eutrophication, deoxygenation, pollution, coastal development, and an enormous fishing pressure at global scales including the Northern Adriatic Sea (Justić 1991; Jackson et al. 2001; Airoldi and Beck 2007; Occhipinti-Ambrogi 2007; Roberts 2007; Zillén et al. 2008; Rabalais et al. 2010; Coll et al. 2012; Tomašových and Kidwell 2017).

Several paleoecological studies have demonstrated effects of overfishing, eutrophication, hypoxia, and pollution on ecosystem composition in the northern Adriatic over the past few centuries based on historical and sediment core data (Barmawidjaja et al. 1995; Sangiorgi and Donders 2004; Ferretti et al. 2013; Kowalewski et al. 2015; Vidović et al. 2016; Gallmetzer et al. 2017; Tomašových et al. 2017, 2019). However, quantitative analyses of benthic assemblages in sediment cores capturing time intervals preceding the past 500 years are rare, and the major ecological turnovers in the northern Adriatic during the entire Holocene remain poorly explored, especially in the north-eastern part (Fiorini and Vaiani 2001; Scarponi and Kowalewski 2007; Rossi and Vaiani 2008; Lotze et al. 2011; Coll and Lotze 2016).

In view of the long-lasting human influence and in the absence of long-term stratigraphic records, it is difficult to reconstruct ecological baselines defining the ecosystem prior to major anthropogenic changes, and to disentangle the effects of long-term ecological drivers such as the post-glacial sea-level rise or climatic oscillations from human interference (Kidwell 2015; Martinelli et al. 2017; Casebolt and Kowalewski 2018).

Here, based on multiple marine sediment cores that capture the entire Holocene in the northern Adriatic Sea, we investigate down-core shifts within the molluscan community and use them as indicators for changes in the marine environment (Anderson et al. 1998; Poirier et al. 2009; Cramer et al. 2015; Gallmetzer et al. 2017). Sediment analyses for heavy metals, organic pollutants, and nutrient concentrations (proxies for human impacts) supplement the findings from the molluscan death assemblages and help in discriminating between natural and anthropogenic drivers of observed community changes. Geochronological data are based on two methods: (1) ²¹⁰Pb of sediment, which provides constraints on apparent sedimentation rates over the past ∼ 120 years, with rate estimates being sensitive to bioturbation (Peng et al. 1979), and (2) Radiocarbon-calibrated amino acid racemization (AAR) dating of mollusk shells from several layers evenly distributed along the sediment cores. These methods allow for a temporal classification of community shifts and also for an assessment of the degree of mixing between layers, one of the major challenges in inferences based on subfossil assemblages (Scarponi et al. 2013; Tomašových et al. 2014; Kosnik et al. 2015; Ritter et al. 2017; Kowalewski et al. 2018).

The four sediment cores used for this study (each 1.5 m long) are from widely separated sampling stations and are thus representative of a broad range of low-sedimentation settings in the northern Adriatic. They not only help draw a comprehensive picture of the ecological history of the whole basin, but also illustrate the heterogeneity of its soft-bottom habitats.

The long time-span embraced in the relatively short sediment columns of 1.5 m has repercussions for the temporal resolution of the ecological history preserved in the cores. In addition to low sediment accumulation rate, bioturbation complicates the detectability of shorter-term shifts by disrupting original stratification and mixing shells of considerably different age (Kosnik et al. 2007; Tomašových et al. 2017). This time-averaging, however, also implies that the environmental and ecosystem trends detected in the cores represent a strong and conservative signal, one that survived in spite of mixing (Tomašových and Kidwell 2010).

The geochronology and molluscan abundance data of the cores from the Brijuni Islands, Croatia, and from the two stations off the coast of Slovenia (stations Piran 1 and 2) were described by Schnedl et al. (2018), Mautner et al. (2018), and Tomašových et al. (2019). Here, we describe age distributions from a fourth station off the Venice lagoon, Italy, and comprehensively analyze stratigraphic (within stations) and geographic variation (among stations) in molluscan composition.

We use molluscan communities as surrogates of macrobenthic soft-bottom assemblages (Olsgard and Somerfield 2000; Kidwell 2002; Albano et al. 2016; Tyler and Kowalewski 2017), and trace long-term shifts in their composition from the beginning of the marine transgression (∼ 11,500 years BP at 40 m depth and ∼ 9,000 years BP at 20 m depth; Colantoni et al. 1979; Storms et al. 2008; Schnedl et al. 2018) through the maximum ingression (∼ 7,000 –5,500 years ago; Trincardi et al. 2011; Amorosi et al. 2017) to the highstand communities established ∼ 5,500 years ago. Controlling for centennial-scale climatic oscillations, such shelly assemblages can be considered as baseline communities, a term used here to denote the range of compositional states that characterized benthic communities at the onset of the first major human impacts in recent centuries. At the same time, our approach reveals the turning points in Holocene community developments marked by anthropogenic pressures of various kinds, and helps quantify the amount of species turnover or habitat degradation in relation to former ecosystem conditions. As such, our work also provides clues for possible restoration or conservation efforts by identifying ecological reference conditions.

The shallow northern Adriatic Sea (average depth 35 m) is delimited from the mid-Adriatic basin by a 50-m isobath extending approximately from the southern tip of the Istrian Peninsula to Ancona on the Italian coast. It is one of the few modern epicontinental seas comparable to typical Paleozoic and Mesozoic shelf environments (McKinney 2007). The inflow of the Po River and other smaller rivers along the western margin leads to high sediment and nutrient input, making this area one of the most productive in the generally oligotrophic Mediterranean (Zavatarelli et al. 1998). Because of the prevailing counterclockwise currents, nutrients and sediments are asymmetrically distributed along a west-east gradient, with meso- to eutrophic waters in the western and northern part and oligotrophic conditions in the eastern regions (Smodlaka 1986; Zuschin and Stachowitsch 2009). The trophic state fluctuates seasonally due to stable water column stratification in summer and pronounced mixing processes caused by winter storms (Giani et al. 2012). In general, the northern Adriatic seafloor is formed by fine-grained sands on the inner shelf, by silts and clays on the middle shelf, and by relict Pleistocene sands on the outer shelf, with the thickness of the Holocene muds decreasing strongly with distance from the coastline (Frignani et al. 2005). Most terrigenous muds are deposited as a belt along the western coast, where they accumulate at annual rates locally exceeding 1 cm/year (Palinkas and Nittrouer 2007; Tesi et al. 2011; Albano et al. 2017). In contrast, the eastern karst-dominated coast is sediment starved (Colantoni and Mencucci 2010) and the Holocene marine sedimentary cover in the southern Gulf of Trieste also tends to be very thin (Trobec et al. 2018).

The four sampling stations were selected to cover different sediment types and degrees of protection from bottom trawling (Fig. 1). Station “Venice” is located about 18 km south-east of the Venice Lagoon in Italian waters at a depth of 21 m (45.375233°N, 12.775783°E). The sediment is almost pure sand and is covered by large amounts of shell debris interspersed with rhodolites and occasional sessile and vagile epifauna (ascidians, sponges, hermit crabs, scallops) (Fig. 2A). Station “Brijuni” is located within the boundaries of the Brijuni Islands National Park (Croatia), at the southern tip of Veliki Brijuni Island, ∼ 400 m from shore, at a water depth of 44 m (44.885833°N, 13.747017°E). The seafloor consists of muddy sand and, at the time of sampling, showed signs of mucilage cover, Lebensspuren, and half-buried shell fragments (Fig. 2B). This location has been protected from fishing activity since 1983 when the Brijuni archipelago became a Marine Protected Area (Fatovic-Ferencic 2006).

The two “Piran” stations are located in Slovenian waters a few km NW of the city of Piran. Station “Piran 1” is about 2.3 km from shore, at 22 m depth, within the protected perimeter of an oceanographic buoy operated by Marine Biology Station Piran (45.548867°N, 13.550900°E). Station “Piran 2” lies 2 km further away from shore, at a similar depth (22.7 m), in unprotected waters (45.563200°N, 13.537033°E). Both stations are characterized by muddy sand with a high proportion of shell debris and by the presence of epifaunal multi-species clumps mainly consisting of sponges and ascidians with a rich associated vagile fauna (ophiuroids, crustaceans, mollusks) (Fig. 2C, 2D).

Changes in molluscan composition at Piran 2 and Brijuni were analyzed by Mautner et al. (2018) and Schnedl et al. (2018).

In summer 2013, at each station, three sediment cores were taken close to each other using an UWITEC piston corer with hammer action (Gallmetzer et al. 2016). Core length was ∼ 150 cm, and core diameter was 9 cm for the two cores designated for sediment and geochemical analyses, and 16 cm (or 9 cm, in the case of stiff sediments at the Venice and Brijuni stations) for the core intended for sieving and analysis of molluscan shell communities. This third core was sliced directly on board the coring vessel immediately after retrieval. The uppermost 20 cm were cut into 2-cm thick increments for higher resolution, and the rest of the core was cut into 5-cm thick increments. For better comparability of the surface section with the rest of the core, the 2-cm increments were pooled into 4-cm increments in geochronological and in some paleoecological analyses.

In the lab, core subsamples were sieved through a 1-mm mesh, and all the shells counted and identified to species level. For bivalves, only complete valves or fragments with preserved umbo were considered. The final count for each increment was achieved by adding the higher number of single valves (either right or left) to the number of articulated specimens. For gastropod shell fragments, the presence of the shell apex was the inclusion criterion. In analyses of absolute abundance, shell numbers from the 9-cm-diameter cores (stations Brijuni and Venice) were multiplied by 3.16 (the factor by which a 16-cm core increases in volume) to allow comparability of patterns in molluscan absolute abundances with the thicker cores from the Piran stations. We note that this extrapolation was not used in analyses of species richness. Similarly, in the molluscan core from Piran-2 station, where the increments from 10 to 35 cm and from 35 to 45 cm were split into quarters and half, respectively, due to high shell abundance, shell numbers were multiplied by four and by two, respectively, to account for comparable volumes. Molluscan species were assigned to functional categories with regard to feeding guild and weed-association (“weed” is used here as a general term for macroalgae and seagrasses) in order to capture down-core shifts in mollusk communities not only in terms of species turnovers and abundance trends but also on the functional level. Likewise, if possible, species were assigned to typical biocoenoses to allow for detection of biocoenotic changes in the cores (Table 1). A dataset listing down-core abundances and guild/biocoenosis assignments for each species, together with the literature used for ecological species characterization, is provided as an electronic supplement, Dataset S1, of the Online Supplemental Files.

One 9-cm-diameter core from each sampling station was used for granulometric analysis and for ²¹⁰Pb sediment dating. At the Venice station, this core was considerably shorter (only 85 cm) than the core used for mollusk shell analysis, most likely due to the presence of paleorelief between the two adjacent coring spots (fluviatile clays occur on the bottom of both cores, marking the onset of the Holocene marine ingression). Core slicing followed the same procedure as for the molluscan core. From each core increment, a 10–20 g subsample was taken and screened through a set of five sieves (1 mm, 0.5 mm, 0.25 mm, 0.125 mm, and 0.063 mm). The fractions were dried and weighed, and the fraction < 0.063 mm was further processed using a SediGraph III 5120 Particle Size Analyzer following standard procedures. We defined four grain size classes: clay (< 0.002 mm), silt (0.002–0.063 mm), sand (0.063–2 mm), and gravel (2–63 mm). The sand was largely of terrigenous origin at Venice, whereas sand and gravel fractions consisted mainly of shells, shell fragments, and other biogenic material at the other stations.

A separate 9-cm core from each station was sealed after retrieval and further processed at the Institute of Marine Sciences in Venice (ISMAR) for geochemical analysis. First, radiographic images were taken of the intact cores to detect down-core density trends and to identify the core increments most suitable for subsample extraction and geochemical analysis. Then, in each core, 10 increments (only eight in the Venice core, which was shorter) were chosen, three within the uppermost 10 cm (0–2 cm, 4–6 cm, and 8–10 cm) and seven evenly spaced along the rest of the core, at depths slightly deviating between cores from individual stations. For each of these increments, three groups of sediment parameters were measured using the protocol described in Vidović et al. (2016): concentration of metals (Hg, Cu, Cr, Ni, Pb, As, Cd, Zn in mg/kg dry weight), concentration of nutrients (total carbon, total organic carbon, and total nitrogen, given in percentage of dry weight), and concentration of organic pollutants (polychlorinated biphenyls “PCBs” and polycyclic aromatic hydrocarbons “PAHs” in ng/g dry weight) (see Online Supplemental File, Dataset S2).

Core chronology is assessed with two different methods: ²¹⁰Pb dating of sediment, and radiocarbon-calibrated amino acid racemization (AAR) dating of mollusk shells.

In the uppermost 40 cm of each core, the activities of ²¹⁰Pb and ²²⁶Ra were analyzed in 2-cm increments down to 20 cm core depth and in 5-cm increments down to 40 cm depth using a High Purity Germanium detector system. Apparent sedimentation rates were computed from the slope of excess ²¹⁰Pb according to the Constant Flux-Constant Sedimentation model (Sanchez-Cabeza and Ruiz-Fernández 2012). A surface mixed layer of 4–6 cm thickness at Piran and Brijuni with uniform or fluctuating ²¹⁰Pb-excess values was excluded from the computation of sedimentation rates. ²¹⁰Pb excess fell to supported values between 10 and 20 cm core depth where the decay profiles used in the slope estimation were terminated (Fig. 3). In the core from Venice station, the upper 10 cm (assigned here to the surface mixed layer) are characterized by ²¹⁰Pb excess values smaller than the down-core-declining values between 10–25 cm, and supported values appear for the first time below 60 cm.

From each of the four molluscan cores, a bivalve species occurring at high abundance throughout the core was used for shell dating by radiocarbon-calibrated amino acid racemization. This was the lucinid Lucinella divaricata for Venice, the venerid Gouldia minima for the two Piran stations (in addition, for Piran 2, shells of the corbulid Corbula gibba were dated), and the venerid Timoclea ovata for Brijuni. For C. gibba, the calibration is described in Tomašových et al. (2017), for G. minima in Tomašových et al. (2019), and for T. ovata in Schnedl et al. (2018). The calibration for L. divaricata is presented below.

The chronology of the Venice core is based on the shell ages of 341 specimens of Lucinella divaricata from 12 core increments (Online Supplemental File, Dataset S3), and of two specimens of the scallop Aequipecten opercularis collected in the 100–105-cm increment (Fig. 4). From 12 selected specimens of L. divaricata, shell parts were allocated for AMS ¹⁴C dating. To avoid contamination, 30% of the outer shell mass was removed prior to AMS analysis in an ultrasonic bath and in 0.5M HCl, and treated in 15% H₂O₂ for 10 minutes, again in an ultrasonic bath (Goslar et al. 2004). Conventional ¹⁴C ages were calculated using a correction for isotopic fractionation (Stuiver and Polach 2016). Conventional ages were converted to calendar years using Calib7.1 (Stuiver and Reimer 2006), the Marine13 calibration curve (Reimer et al. 2013), and a regional marine reservoir correction (ΔR) for the northeast Adriatic (Rovinj) equal to = -61 years (standard deviation = 50 years) (Siani et al. 2000) (Table 2).

Small portions of the 347 dead and of three live-collected shells of L. divaricata were analyzed for the extent of amino acid racemization (AAR) at Northern Arizona University using reverse-phase high-pressure liquid chromatography (RP-HPLC, Kaufman and Manley 1998). Specimens were leached 20% by weight with a dilute solution of HCl and the resulting solutions were hydrolyzed at 110°C for six hours to release amino acids from their peptide chains. Concentrations and D/L values of aspartic acid (Asp), glutamic acid (Glu), serine (Ser), and alanine (Ala) were measured. Asp and Glu were used in age calibrations. Six specimens were flagged as outliers as described by Kosnik and Kaufman (2008) and removed from analyses.

We calibrated the rate of AAR in L. divaricata on the basis of the 12 ¹⁴C-dated and three live-collected shells according to Allen et al. (2013). Asp and Glu D/L values were fit using four mathematical functions to model the relation between postmortem age and D/L values, and two uncertainty models (lognormal and gamma) using the R statistical computing language (R Core Team 2013). The four mathematical functions include time-dependent rate kinetics (TDK), constrained power-law kinetics (CPK), simple power-law kinetics (SPK), and apparent parabolic kinetics. Each model was fitted with and without a non-zero initial D/L value. The time-dependent reaction kinetic model (TDK, Allen et al. 2013) for Asp D/L was the best model with the smallest BIC (Table 3), with the initial D/L value estimated from data (TDK1). This model was used in the final calibration. The calibration equation is a*arctanh ([DL - DL₀]/[1 - DL*DL₀])^b, where DL is Asp D/L, DL₀ is Asp D/L at 0 years (here, 0.04479), and a = 43589 and b = 2.342312. The uncertainty is defined by the log-normal distribution, with the mean equal to log-transformed age estimate (in years) and the variance parameter of the log-normal distribution equal to 0.09134. AAR-calibrated calendar ages are set relative to the year of collection (2013 AD = year zero).

To summarize, between nine and 13 more or less evenly spaced core increments were chosen at each station (nine for Piran 1, eleven for Piran 2, 12 for Venice and 13 for Brijuni), and from each of these increments, ∼ 30 specimens (complete valves or fragments with preserved umbo) were randomly selected for AAR analysis. Either only right or only left valves were used per increment to avoid double-dating of the same individual. In a few increments with lower shell abundance, both right and left valves had to be used. In total, 564 specimens of G. minima (Piran 1 and 2), 232 specimens of C. gibba (Piran 2), 300 of T. ovata (Brijuni), and 341 of L. divaricata (Venice) were used in shell-based core geochronology. AAR measurements were calibrated with ¹⁴C ages of 13 specimens of G.-minima (Tomašových et al. 2019), 12 of L.-divaricata (this study), and 13 specimens of T.-ovata (Schnedl et al. 2018), respectively. Median age and interquartile age range (corrected for calibration age error according to Dominguez et al. 2016) were computed for each increment. The increments with similar grain size, lithology, and proportion of major rock-forming components, including concretions, bryozoans, coralline algae, and mollusks were grouped into stratigraphic units. For each unit, 25^th and 75^th age percentiles and median age based on shell-age frequency distributions were computed (Table 4). Age distributions significantly overlap between some units, most likely due to bioturbation (Fig. 5). Therefore, geochronology of these units is better defined by age range (bounded by the 25^th and 75^th age percentiles) rather than by a single age estimate.

To compare long-term ecological trajectories among stations, core increments (2 or 5 cm thick layers) at each station were assigned to four phases using the sequence-stratigraphic terminology of model III of Catuneanu et al. (2009) where the highstand phase terminates prior to the base-level fall (Hunt and Tucker 1992): (1) A transgressive phase (TST) ranging from the onset of the marine transgression about 10,500–11,000 years ago (at locations that are now at ∼ 40 m depth; Storms et al. 2008) or 9,000 years ago (at locations now at ∼ 20 m depth; Picone et al. 2008; Trincardi et al. 2011), correlating with the period of global sea-level rise, to the point of maximum flooding ∼ 7,000 years ago (Antonioli et al. 2007; Vacchi et al. 2016; Amorosi et al. 2017); (2) A maximum flooding zone (MFZ) marking the maximum ingression and the transition from the late transgressive to the early highstand phase (∼ 7,000–5,500 years ago; Trincardi et al. 2011; Amorosi et al. 2017) and documented by the establishment of shell beds with epifaunal bivalves at Brijuni and Piran; (3) An early highstand phase (eHST) lasting from the end of the transgressive phase to the onset of the stillstand (∼ 2,000 years ago), coinciding with the beginning progradation of the Po, Isonzo, and Tagliamento Rivers (Amorosi et al. 2008; Zecchin et al. 2015; Amorosi et al. 2017), and with the appearance of distinct signs of human interference with the marine environment; and (4) A late highstand (effectively anthropogenic) phase (lHST) still underway, marking the era of pronounced marine ecosystem changes induced by human activities in the recent past, including resource exploitation, pollution, and eutrophication (with recent oligotrophication since the 1990s; Mozetič et al. 2010; Djakovac et al. 2012).

Environmental variables used in multivariate analyses of molluscan assemblage composition are represented by: (1) percentages of clay, silt, sand, and gravel; (2) sediment concentrations of Hg, Cu, Cr, Ni, Pb, As, Cd, and Zn; (3) concentration of total organic carbon (TOC); (4) concentration of total nitrogen (NT); and (5) concentrations of organic pollutants (PAHs and PCBs). These variables were log-transformed if not normally distributed, and z-standardized to account for different scales and units. Species diversity is given as (1) the observed number of species, and (2) true diversity or effective number of species (the exponential of the Shannon entropy, exp(H)), which transforms the standard diversity index (H) into a linear scale representing the number of species in a perfectly even community (Jost 2007). Abundances were rarefied to a sample size of 150 individuals per increment, and, if necessary, extrapolated to this number according to Chao et al. (2014) using the R package “iNEXT” (Hsieh et al. 2016). Down-core differences in true diversity among the four phases are tested using the Wilcoxon rank test based on median effective species richness.

Differences in molluscan community composition between consecutive phases of individual cores are quantified by the pseudo-F statistics of permutational multivariate analysis of variance (PERMANOVA) (Anderson 2001). Non-metric multi-dimensional scaling (NMDS) based on Bray-Curtis distances (Ward and Hutchings 1996; Clarke and Warwick 2001) and fit to two dimensions is used to visualize changes or differences in mollusk community structure (1) between the four stations and (2) between five (Venice and both Piran stations) or six (Brijuni) stratigraphic units within stations (down-core changes). Bray-Curtis distances are based on square-root transformed relative abundances. Core increments with an overall mollusk shell abundance < 70 (Brijuni 0–4 cm, Venice 130–140 cm) and the increment 4–8 cm at Birjuni (outlier in the ordination of weed-association) were excluded from this analysis. Ordinations are performed not only for species abundances but also for abundances of ecological groups based on feeding type, weed-association, and biocoenotic affiliation.

To explore the sources of compositional variation in the ordination space, correlations between NMDS-axis scores and environmental variables were measured. The amount of community variation explained by environmental variables is quantified by redundancy analysis (RDA) based on a forward selection model (“ordi2step” function in R). Just as for the NMDS analyses, RDA is based on square-root transformed relative abundances—i.e., Euclidean distances used by RDA correspond to Hellinger distances, which are not sensitive to the double-zero problem, and are efficient in modeling species responses to short environmental gradients (Legendre and Gallagher 2001; Legendre and De Cáceres 2013). We note that constrained analyses of principal coordinates based on Bray-Curtis distances generate very similar results. The application of RDA, however, is limited to those ten core increments (eight in the core from Venice station) for which measurement data of environmental variables are available. To compare the proportion of compositional variation explained by stratigraphic differences within the four stations (defined by sediment depth of a given assemblage, at the scale of 5-cm increments) with the amount explained by geographic differences among the stations (defined by latitude and longitude of each station), we extend redundancy analysis to variation partitioning based on Hellinger distance (Peres-Neto et al. 2006). Statistical analyses were performed in R-studio (version 3.1.3.) using the “vegan” package (Oksanen et al. 2016). The R-code and two related input files are accessible in the supplementary material.

The sampling regions covered by the four sediment cores (NE, SE, and central-western part of the northern Adriatic Sea) represent different depositional environments and stratigraphic histories. Only the two cores from the neighboring Piran stations show similarities in the stratigraphic distribution of percentages of sand and gravel, and in the degree of shell packing.

At the bottom of the Venice core, the Holocene transgressive surface is marked by a sudden change from dark, organic-rich fluviatile or paralic clays to marine sediments containing shells of both marine and freshwater species. The sediment consists almost exclusively of pure, well-sorted siliciclastic sands throughout the core (Fig. 6). The down-core age distribution of L. divaricata shells (Fig. 4) shows that median ages of this species range between 200 and 1,100 years in the upper 30 cm, between 900 and 2,100 years at 40–70 cm, between 3,500 and 3,700 years at 90–105 cm (¹⁴C ages of two specimens of Aequipecten are equal to 2,700 years), and between 6,300–6,800 years at 110–130 cm. The interquartile ranges of L. divaricata ages (corrected for age error) in 5-cm increments span ∼ 500 to 2,000 years. Considerable sediment mixing is suggested by similarities in median ages of 5-cm increments within the five stratigraphic units. Down-core changes in median L. divaricata ages thus indicate that net sedimentation rates were small (∼ 0.01–0.05 cm/year) during the entire Holocene period. An abrupt shift in median age between 100–105 and 110–115 cm from 3,500 to 6,300 years indicates a major hiatus in sediment accumulation. Fitting the ²¹⁰Pb excess values against sediment depth at 10–50 cm indicates that apparent sedimentation rates are 0.88 cm/year, much higher than estimates based on down-core changes in median ages of L. divaricata over the same stratigraphic interval (∼ 0.05 cm/y, with median shell age equal to 1,000 years at 45–50 cm).

The Venice core can be subdivided into four units based on down-core changes in the composition of biogenic sediments (Fig. 7): (1) In the deepest part (135–105 cm) representing the TST, in addition to the dominant molluscan and rhodolithic fragments, foraminiferal tests of the species Ammonia cf. beccarii and Massilina cf. secans and carbonate concretions are abundant. Remains of fresh-water gastropods are frequent, mainly opercula of the genus Bithynia, reflecting pre-transgressional freshwater or brackish environments; (2) Between 105 and 90 cm (MFZ), A. cf. beccarii and M. cf. secans decrease markedly in abundance; in shallower core depths, they are completely absent. Carbonate concretions are larger compared to the TST; (3) The unit between 90 and 30 cm (eHST) is characterized by a high abundance of the small irregular sea urchin Echinocyamus pusillus; and (4) In the uppermost unit (lHST), from 30 cm up to the surface, the molluscan gravel co-occurs with numerous branched rhodolites, which are also a constituent part of the living assemblage (Fig. 2A).

A dated terrestrial plant fragment from the base of the Brijuni core is ∼ 11,350 years old. The median age of Timoclea ovata shells from dated increments increases monotonically down-core, ranging from 160 years at 0–4 cm core depth to 6,100 years at 150–155 cm (Fig. 5), and a maximum shell age of 10,100 years. Net sedimentation rates based on differences in median ages of T. ovata between increments range between ∼ 0.01–0.04 cm/year. Interquartile age ranges (corrected for age error) of T. ovata in 5-cm increments range mainly between 500 and 1,000 years. The core is stratigraphically subdivided into six units: the TST at the bottom of the core (160–120 cm), the MFZ represented by a shell-bed at 120–90 cm with shells heavily coated by encrusting bryozoans (large ¹⁴C-dated shells of Arca, oysters and Glycymeris from this shell bed span between 5,500 and 7,800 years; Schnedl et al. 2018), an eHST with coarser sediment in the lower and finer sediment in the upper part (90–50 cm and 50–20 cm, respectively), and a lHST (or “anthropogenic” phase, 20–0 cm) at the core top, with the uppermost 10 cm enriched in organic pollutants, total nitrogen, and total organic carbon. Grain size composition changes strongly up-core, with a large proportion of molluscan gravel in the lower part, bryozoan-molluscan sand in the upper part, and a higher contribution of coralline algae in the uppermost part (Figs. 6, 8). The overall trend is fining upward, with sediment-starved conditions at the bottom of the core during the marine transgression phase, and a marked shell bed with encrusted large bivalves between 120 and 90 cm. As in the Venice core, shells of terrestrial (nine species) and freshwater (15 species) gastropods are abundant in the TST but decrease steadily up-core. Molluscan sand and gravel peak between 130 and 95 cm core depth, with a maximum percentage of ∼ 80%, and then gradually decrease upwards. In the uppermost 20 cm, silt and clay constitute ∼ 80% of the sediment. According to ²¹⁰Pb sediment dating, the apparent sedimentation rate during the deposition of increments between 4 and 16 cm core depth has been ∼ 0.15 cm/year, exceeding the estimates of sedimentation rate based on down-core changes in median ages of T. ovata by one order of magnitude. Sand and gravel, in this core, are mostly of biogenic origin and consist mainly of molluscan shell hash, exoskeletal parts of crustaceans (present mainly in the deeper core sections), and bryozoan colony fragments. Bryozoans increase in abundance from 90 cm up to the surface, showing a gradual up-core shift from encrusting to bushy (Cellaria spp.) growth forms. The bivalve shell bed between 120 and 90 cm depth consists of large, strongly bored shells of glycymeridids, scallops, and oysters, with common encrustation by coralline algae, serpulids, and bryozoans (Fig. 8).

Five stratigraphic units can be distinguished at the Piran stations (Mautner et al. 2018; Tomašových et al. 2019). In contrast to Venice and Brijuni, the TST unit is more expanded and the HST units are very thin. Stratigraphic mixing coupled with a temporally disjunct production of two species (e.g., Broecker et al. 1999) generates significant age offsets between median ages of Gouldia minima and Corbula gibba in the same increments (Tomašových et al. 2019). At Piran 1 and 2, the median age of G. minima shells increases steadily down-core from ∼ 1,000 years at 0–4 cm to ∼ 5,000–6,500 years at 90–100 cm, and then it remains relatively constant down-core. In contrast, the median age of C. gibba increases from ∼ 1,000 years at 0–4 cm to ∼ 7,000–8,000 years at 90–100 cm. Similar to Venice and Brijuni, the cores can be split into five major stratigraphic units based on changes in grain size: (1) the early TST at 65–150 cm formed by muddy sand; (2) the late TST formed by molluscan sand and initial oyster-rich levels at 65–35 cm; (3) the MFZ formed by the lower part of a shell bed (mainly oysters) at 35–16 cm; (4) the eHST corresponding to the upper part of the shell bed (mainly Arca noae) at 16–8 cm; and (5) the lHST represented by the topmost 8-cm-thick molluscan sandy mud.

Stations Piran 1 and Piran 2 show similar grain size composition and stratigraphy despite being ∼ 2 km apart. Both stations share a coarsening-upward trend in grain size and an up-core increase in average shell size (Fig. 6; Mautner et al. 2018). The proportion of mud in the deeper part of the core (150–75 cm) is 70-80% and decreases to less than 20% in the uppermost 25 cm in favor of a sand and gravel fraction consisting mainly of mollusk shells, with the addition of echinoderm and crustacean fragments and foraminiferal tests. At both stations, the deepest core intervals (155–135 at Piran 1 and 145–135 at Piran 2) contain a high number of Modiolus cf. adriaticus. The densely packed shell bed with large shells of epifaunal bivalves (Arca noae, Ostrea sp., Modiolus sp.) between 8 and 35 cm is a distinct feature in both cores (Fig. 9). The uppermost 8 cm consist of molluscan sands with dispersed coralline algae (Lithothamnion spp.), which locally expose the shells of the underlying shell-bed-building species at the sediment surface (Fig. 2C, 2D). At Piran 1, a loosely packed shell bed is preserved between 95 and 35 cm, formed by ostreid and mytilid shells (Ostrea sp, and Modiolus sp.). As at Venice and Brijuni stations, single-species interquartile age ranges at the Piran stations range mainly between 1,000–2,000 years (Tomašových et al. 2019). The average sedimentation rate based on the excess ²¹⁰Pb profiles between 4–14 cm amounts to 0.29 cm/y at station Piran 1 and 0.16 cm/y at Piran 2, again significantly exceeding the net sedimentation rates indicated by down-core changes in median age of G. minima and C. gibba in the upper 20 cm (Tomašových et al. 2019).

Sediment contamination with heavy metals and organic pollutants increases upwards at all stations, with marked peaks in the uppermost 20 cm (Fig. 10). The strongest contamination is observed at Brijuni, the lowest at Venice. The measured concentrations exceed NOAA benchmark values (Long et al. 1995) in several cases: Hg concentrations at Brijuni are higher than Effects Range-Median (ERM) level (at Piran 2, they are close to this threshold); Zn exceeds Threshold Effects Levels (TEL) values at all stations except for Venice. Pb, PAH and PCBs concentrations are below benchmark values at all stations, but clearly increase in the uppermost increments. TOC and TN concentrations increase in the uppermost 10–20 cm at Brijuni and Piran 1 (at Piran 2, only total nitrogen increases upwards), whereas at Venice they remain low.

The four sediment cores yielded ∼ 98,700 shells from 52,710 mollusk individuals belonging to 341 different species, including 180 gastropods, 120 bivalves, 9 chitons, five scaphopods, and 27 non-marine mollusks (17 freshwater and 10 land snail species). The highest density (individuals per 100 cm³) occurs in the eHST of the two Piran cores (upper part of the shell bed, maximum 618 ind/100 cm³ at Piran 2 and 314 ind/100 cm³ at Piran 1) and in the TST of the Venice (231 ind/100 cm³) and Brijuni cores (430 ind/100 cm³) (Fig. 11A).

Stratigraphic changes in mollusk density at the Piran stations on one hand and at Brijuni/Venice on the other hand show opposite trends: At Piran, density is relatively constant below 80 cm core depth, rises slowly from this point upwards, and experiences a burst above 40 cm. At Venice and Brijuni, the highest density occurs in the lowermost part and then decreases continuously up to the top. All cores except for the Venice station display a sudden density drop in the uppermost 10 cm. In all cores, the raw species richness (Fig. 11B, 11C) increases steeply immediately above the transgressive surface. At Venice and Brijuni, it peaks around 120 cm core depth (90–100 species). At Venice, raw species richness decreases slightly and stabilizes at ∼ 50 species around 60 cm depth without major changes up to the surface. At Brijuni, high species richness extends up-core until ∼ 40 cm depth, from where it declines, with a sharp drop in the uppermost 12 cm. At the two Piran stations, after a first steep increase in the early TST, raw species richness grows steadily up-core with minor fluctuations until the near-surface unit (maxima: 116 species at Piran 2 and 120 species at Piran 1) where, similar to the other stations, a sudden decline sets in, which is more pronounced at Piran 1 than at Piran 2.

Effective number of species (exp(H)) is highest at Brijuni station (median: 30), followed by Venice (28), Piran 2 (25), and Piran 1 (24) (Fig. 11D). Here, the most conspicuous stratigraphic trends are an extensive peak in the MFZ with subsequent sharp drop in the eHST at Brijuni, and a marked drop in the uppermost 30 cm in both Piran cores. This sudden diversity loss at Piran significantly distinguishes the HST phases from the MFZ and the TST at both stations (Table 5). At Brijuni, the effective number of species is significantly higher in the eHST than in the MFZ/TST; in the lHST, diversity drops again to low values similar to the MFZ. At Venice, diversity fluctuates considerably throughout the core without significant pairwise differences between the four units (Table 5). A diversity comparison between the same stratigraphic units in all four cores does not show any significant differences between the Piran cores and only partial differences between Piran and the two other stations and between Brijuni and Venice. Differences concern mainly the early and late HST phase and are determined by the marked diversity loss in the uppermost increments at Piran and the diversity increase during eHST at Brijuni station (Online Supplemental File Table S4).

The lucinid bivalve Lucinella divaricata, which occurs preferentially in sediments of seagrass meadows is the most frequent species throughout the core (Fig. 7). Absolute abundances of this species are much greater in deep than in shallow units, indicating particularly suitable habitat conditions during the TST and the MFZ. The existence of seagrass environments during this period is supported by the strictly seagrass-associated rissoid gastropods Rissoa monodonta and R. membranacea, which occur exclusively in the TST below 100 cm core depth. The core is further dominated by the epifaunal bivalves Anomia ephippium and Aequipecten opercularis that peak at 90–95 cm, the infaunal bivalve Donax venustus reaching a maximum in the TST, and by the tellinid Moerella distorta peaking at 60–65 cm. From the MFZ upwards, epifaunal species slightly decline and infaunal species associated with sandy bottoms increase in abundance (e.g., Tellina spp.). Parvicardium minimum becomes the second most abundant species, and the abundance of Glycymeris glycymeris, a species characteristic of sandy bottoms exposed to currents, increases from the MFZ onwards. In addition, in the intermediate core section, a slight increase of deposit feeders can be observed. In the lHST, the ratio between infaunal (the dominant type of substrate relation) and hard-bottom epifaunal species remains largely constant while, in terms of feeding guilds, deposit feeders decline in favor of filter feeders. Overall, this core is characterized by a high abundance of seagrass-associated species in the TST and by species adapted to well-sorted, nutrient-poor sands in the HST phases. This is also reflected in the biocoenotic structure, with a dominance of species from the detritic community and important proportions of representatives from the biocoenoses associated with photophilic algae and fine, well-sorted sands (AP, SRPV, and SFBC; see Table 1). The community appears relatively homogenous over wide sections of the core, and down-core changes on the functional or biocoenotic level are not particularly marked (Table 5), possibly due to pronounced mixing processes that generated highly time-averaged assemblages.

The TST and the MFZ are dominated by epifaunal (e.g., Striarca lactea) and vegetation-associated mollusks, of which a large proportion (61%) is strongly but not exclusively related to seagrasses (e.g., Alvania geryonia, Alvania cimex, Rissoa violacea) (Online Supplemental File Fig. S5). A shell bed with large shells of oysters, scallops, and glycymeridids heavily coated by encrusting bryozoans characterizes the MFZ (Fig. 8). Above this shell bed, erect bryozoans abruptly become a constitutive part of the biogenic gravel, and seagrass-associated species decline in favor of species with a preference for algae. Infaunal (e.g., Timoclea ovata, Parvicardium spp.) and epifaunal species typical of detritic bottoms (e.g., Bittium submammillatum, Jujubinus montagui, Anomia ephippium) dominate the community in the eHST. In the lHST, the proportion of epifaunal mollusks drops rapidly while infaunal species such as Nucula nucleus and T. ovata further increase. Besides filter feeders (the dominant feeding guild in this core), 20–40 % of all mollusk species are grazers. However, this feeding type experiences a drastic decline in the uppermost 30 cm, which is compensated by a rise in detritivore and filter-feeding species (Fig. 8). In biocoenotic terms, the mollusk community at Brijuni is best described as a coastal detritic community with a large proportion of vegetation-associated species and some representatives of the coralligenous biocoenosis and the current-swept, coarse sands (SGCF). Significant biocoenotic shifts are observed between the lower part of the core (TST and MFZ), the eHST1, and the upper core section (eHST2 and lHST) (Fig. 12B, Table 5), possibly due to a larger number of species characteristic of muddy sands in sheltered areas (SVMC) in the eHST1 and the lack of species typical for coastal terrigenous muds (VTC) in the bottom and top units of the core (Online Supplemental File Fig. S6).

The mollusk communities of the two Piran cores follow highly similar trends. The most abundant species are G. minima and C. gibba, two filter-feeders with an infaunal and semi-infaunal life habit, respectively, characteristic of muddy or muddy-sandy bottoms. Relative abundances of the two species show opposite up-core trends: C. gibba is more abundant in the early TST, where it constitutes up to 20% of the whole community, and G. minima in units with coarser sediments, from the late TST up to the surface (Fig. 9). In the early TST, the mollusk community is dominated by infaunal and epifaunal filter feeders and detritivores typical of muddy bottoms (in addition to C. gibba, e.g., Abra alba and Flexopecten glaber). A third important feeding guild, grazers, is represented by several gastropod species associated with algae and/or seagrasses (e.g., Pusillina spp. and Bittium spp.). Following the coarsening-upwards trend, the dominance of infaunal species recedes up-core in favor of epifaunal species (Fig. 9). At Piran 1, this trend leads to the appearance of large oyster and horse mussel shells (Modiolus sp.) in the upper early and the late TST. In both cores, a major shell bed between 35 and 8 cm depth is primarily formed by large oysters in the lower part and by Arca noae shells (still co-occurring with oysters and Modiolus sp.) in the upper part. Three dated oysters from this core interval vary in age between 5,900 and 6,800 years, showing that the shell bed initiated during the MFZ (Tomašových et al. 2019). This stratigraphic pattern indicates that the establishment of biogenic hard bottoms with epifaunal biostromes was initially patchy, and that the shell bed became spatially more continuous ∼ 6,500–6,000 years ago. In the Piran 2 core, between 12 and 14 cm, the number of fully-grown Arca shells is particularly high. Algae-associated grazers such as Bittium spp. and several rissoid species considerably increase in abundance within the shell bed. The uppermost 8 cm are clearly separated from the underlying units by (1) a decrease in biogenic gravel; (2) a marked decrease in overall abundance; and (3) a re-establishment of an infauna-dominated soft-bottom community similar to the early TST, with opportunistic species such as C. gibba proportionally increasing. On the biocoenosis level, species of the coastal-detritic community are dominant in both cores. However, the transgressive phases are well separated from the upper core mainly due to a larger proportion of species typical for coastal terrigenous muds and for unstable soft bottoms (Fig. 12B, Online Supplemental File Fig. S6).

Variation partitioning based on Hellinger distances shows that a major part of total species and ecological compositional variation observed at the four stations is explained by geographic differentiation (57% with species data, 60% with biocoenoses, 52% with feeding guilds, and 52% with weed-association data), far exceeding the explanatory power of stratigraphic variation at the scale of 5-cm increments (i.e., the explanatory variable is represented by sediment depth; 4% with species data, 5% with biocoenoses, 3% with feeding guilds, and 6% with weed-association data). However, eliminating the geographic component, the contribution of stratigraphic variation to total compositional variation is still higher than under a null hypothesis of no temporal change (F [species] = 12.1, p < 0.001, F [biocoenoses] = 12.1, p < 0.001, F [feeding guilds] = 12.1, p < 0.001, F [weed association] = 17.6, p < 0.001).

Differences in grain size between the Venice station and the eastern sites indicate that, at Venice, stronger bottom currents continuously removed finer sediment fractions. Spatial and stratigraphic trends in TOC and TN are probably determined by higher organic content at times of transgression, generating an up-core decline in TOC at all stations except Brijuni, and by differences in mixing, erasing anthropogenic signals in the uppermost core parts at Venice and at Piran 2. Multivariate analyses of species and functional composition point to significant community differences among stations (with the exception of the two Piran stations) and to stratigraphic differences within stations (Table 5).

The largest inter-site distances emerge in the ordination of species composition (Fig. 12A), without any overlaps among the three regions Venice, Brijuni, and Piran, while the cores from the two Piran stations are not separated. The separation of the three regions follows a gradient in grain size (mud-content) and sediment organic content (TOC, TN) responsible for the detached position of the Venice station with its nutrient-poor, sandy sediments on the left side of the NMDS plot. The strongly diverging species composition in the three regions is reflected in the differences among biocoenoses (Fig. 12B), which appear equally separated in ordination space and further emphasize the heterogeneity of benthic habitats in the northern Adriatic.

The ordination of the four stations according to functional species descriptors (feeding guild and weed-association) leads to a similar, albeit less pronounced compositional separation. The clear separation of Venice with regard to feeding guilds (Fig. 12C) and weed-association is caused by a high proportion of mollusks with chemosymbionts (mainly L. divaricata) and of species strictly associated with seagrasses, respectively (Online Supplemental File Figs. S7, S8). Brijuni and Piran stations are functionally more similar and show almost no separation along the first axis (Figs. 12C, 10D). Highstand assemblages at Piran (both Piran stations congruent) and Brijuni overlap in ordination space owing to the abundance of grazers. This group plays an important role throughout the core at Brijuni, while at Piran they emerge as a constitutive feeding guild only between 10 and 30 cm depth (Online Supplemental File Fig. S6). Similar overlaps of intermediate core depths between Brijuni and Piran are also shown for weed-association (Fig. 12D) due to parallel increasing trends in algae-associated species.

The distinction of different stratigraphic units and subunits within the individual cores is more accentuated in the analysis of species composition and biocoenoses than on the functional level, and more pronounced at Brijuni and Piran stations than at Venice (Table 5). At Piran, the eHST and lHST are not significantly different in any ordination but are always distinctly separated from the other units. These units are also clearly distinguished along the first axis, pointing to significant within-core community shifts (Fig. 12). At Brijuni, only the ordination of species composition significantly resolves all major phases. On the functional and biocoenotic level, eHST 2 and lHST form a homogeneous group, as well as TST and MFZ with regard to feeding guilds and weed-association. At Venice, TST and MFZ merge in all ordinations and differ significantly from the HST phases due to their high percentage of seagrass-associated species (Table 5).

The drivers of up-core mollusk community change identified by forward-selection-model based RDA are either proportions of mud or correlates of pollution and organic enrichment (concentrations of Pb, TN and PCB) that distinguish the uppermost part of the core from the rest. At the two Piran stations, which are characterized by a grain size composition strongly changing up-core, the mud content is the most important variable explaining ∼ 40% of the compositional variation in the dataset (Fig. 13C, 13D). In addition, Pb concentration is selected as a second driver raising the percentage of explained variation to ∼ 48% at Piran 2 (Table 6).

Pb and TN concentrations at Brijuni and concentrations of PCBs at Venice represent the main sources of within-core variation in species composition (Table 6). However, redundancy analyses for each station using only one variable at a time show that several variables are significantly influencing community change (Online Supplemental File Table S9). Many of the elemental concentrations and grain size percentages are strongly correlated (Online Supplemental File Table S10), indicating that the variables picked by forward selection represent two major environmental vectors: (1) correlates of sedimentary processes affecting sedimentation rate, winnowing and physical reworking (grain size), and (2) correlates of organic enrichment and pollution.

Age dating of bivalve shells reveals that at all stations, the sediment cores capture the whole or most of the environmental history since the flooding of the northern Adriatic shelf ∼ 11,000 years ago (at ∼ 40 m water depth, as directly documented at Brijuni; Schnedl et al. 2018) or ∼ 9,000 years ago (at ∼20 m water depth) (Antonioli et al. 2007; Storms et al. 2008; Vacchi et al. 2016). In all cores, the major phases of the Holocene marine history (transgressive, maximum flooding zone, early and late highstand phase) can be distinguished on the basis of stratigraphic changes in species and functional composition of mollusk assemblages, as detected by multivariate analyses (Table 5). However, we note that geographic differences overwhelm within-core stratigraphic variation in composition, most likely due to millennial-scale averaging and between-increment mixing that occurs within stratigraphic units (Fig. 7). For example, unmixing of abundance trajectories based on shell ages of two bivalve species showed that a major nineteenth century decline in abundance of Gouldia minima and a twentieth century increase in abundance of Corbula gibba are significantly muted in the stratigraphic record of highstand units at both Piran stations (Tomašových et al. 2019). The mismatches between ²¹⁰Pb-based apparent sedimentation rates and sedimentation rates based on down-core changes in median shell ages also indicate that the ²¹⁰Pb profiles are steepened by incomplete mixing reaching to 10–20 cm at Brijuni and Piran and to 50 cm at Venice, i.e., below the base of the fully-mixed surface layer (Fig. 3).

The benthic community development in the TST and MFZ follows natural environmental changes such as a deepening of the habitats due to rising sea level and the concomitant establishment of the modern circulation system in the northern Adriatic Sea (Cattaneo and Trincardi 1999), leading to a weakening of the hydrodynamic forces at Brijuni but to intensification of current activity at Piran and at Venice. After sea level stabilization in the eHST, the established communities remained compositionally stable for a few millennia. They can be considered as baseline communities that started to change markedly with the intensification of human impacts on the sea in the lHST. Depending on the sampling region, these baseline communities display specific characteristics.

The lowermost core interval (135–105 cm) encompasses the TST from ∼ 10,000 to 6000 years ago and represents a transgressive lag deposit with a mixture of terrestrial and shallow-marine species (Kidwell 1989; Fürsich and Oschmann 1993; Cattaneo and Steel 2003; Scarponi et al. 2017) with the highest mollusk shell abundances in the whole core. Besides the many opercula of the fresh-water snail Bythinia sp. evidencing the pre-transgressional limnic environments at this site (Amorosi et al. 2008), several species are characteristic for the environmental conditions in the initial phase of marine flooding. The surf clam Donax venustus (associated with sandy intertidal and shallow subtidal habitats), and the shallow-water miliolid foraminifer Massilina cf. secans display the highest abundances in this unit, whereas they are rare in the eHST and completely missing in the lHST. Another shallow-water foraminifer, Ammonia cf. beccarii, indicating shallow-water marginal environments (Felja et al. 2015) and potentially epiphytic on Posidonia oceanica (Debenay et al. 1998), is also abundant. The lucinid bivalve Lucinella divaricata, which is commonly associated with seagrass meadows (e.g., Covazzi Harriague et al. 2006), shows highest abundance in the TST, and the two strictly seagrass-associated rissoid bivalves Rissoa monodonta and R. membranacea (Gofas et al. 2011; Marina et al. 2012) are restricted to this phase. All these taxa contribute to the compositional separation of the lower core section from the upper units (Table 5). They provide evidence for the establishment of seagrass meadows, which probably reached their greatest extent close to the end of the late-Quaternary sea-level rise (Taviani 1980; Correggiari et al. 1996). Dead rhizome mats of P. oceanica still occur in a discontinuous belt along the northern Adriatic Italian coast at water depths between 15 and 23 m (Newton and Stefanon 1982; Rismondo et al. 1997). One of the many localities where such remains were found is located very close to the Venice station (Rismondo et al. 1997, fig. 2). Because these rhizome mats have not been dated so far, their age cannot be specified.

The decline and eventual disappearance of seagrasses has been primarily attributed to anthropogenic impacts such as deterioration of water transparency due to eutrophication or increased sediment load (e.g., Taviani 1980; Rismondo et al. 1997). Our findings, however, indicate that seagrasses at this station had already receded during the early HST and were largely eliminated in the late HST. The invertebrate community in the core increasingly features indicator species for barren sandy seafloors affected by strong currents (e.g., Echinocyamus pusillus (Telford et al. 1983), Glycymeris glycymeris (Royer et al. 2013), and Paphia rhomboides (Savina and Pouvreau 2004)). Species with a strict seagrass association are missing in the upper eHST and lHST record.

The decline of seagrasses can be linked to the establishment of the modern circulation pattern in the northern Adriatic after highstand conditions were reached (Trincardi et al. 1994). In the region of the Venice sampling station and generally along a 10–20 m isobath off the Italian coast, strong bottom currents developed, which may have enhanced the erosion of the meadows and increased water turbidity (Stefanon 1984; Casellato and Stefanon 2008). In addition, extreme climatic events such as exceptional storms and diseases could have further contributed to the recession of P. oceanica independent of anthropogenic disturbances (Duarte et al. 2006). Therefore, the community at Venice station towards the end of the eHST most likely consisted of an assemblage adapted to sandy seafloor influenced by regular, strong tidal currents and occasional storms, and interspersed with areas covered by meadows of P. oceanica.

Similar to Venice station, the Brijuni core contains ample evidence for the transition from a terrestrial to a marine environment and for shallow-water habitats in the early TST. The lowest core unit contains high abundances of reworked, poorly preserved terrestrial and freshwater gastropods and a large proportion of organic material originating from the underlying paralic layer. At the time of the marine transgression, the Brijuni station, which today is located a few hundred meters from the coast, was probably close to the rocky slope of the hill now forming the island of Veliki Brijuni. In accord with this, in the deepest core layers, the dominant mollusk species are characteristic of hard bottoms and marine vegetation, partly with a strong but not strict association with seagrasses, thus reflecting a mixed seagrass-hard-bottom habitat in shallow water. As seen at the Venice station, this lowermost core interval features the highest mollusk abundance in the whole core.

With increasing water depth, the mollusk community shifted to a dominance of large-sized bivalve species with epifaunal (scallops, oysters) or shallow-burrowing infaunal life habit (e.g., Glycymeris glycymeris, Timoclea ovata) and a preference for coarse and detritic bottoms. In the core, this community is preserved in a shell bed containing large and heavily encrusted Glycymeris shells (7,800 years old) that provided a settlement basis for arcids and oysters, which are present only in this core unit (Fig. 14B). This shell bed represents the transition from the late transgressive to the earliest highstand phase (MFZ).

As water depth increased and highstand conditions were reached, more fine-grained material was deposited at higher rate leading to increasing siltation of the seafloor. The high percentage of bryozoan gravel at the end of the MFZ (up from 90 cm) indicates a re-organization of the benthic community: faunal elements, particularly the erect bryozoans Cellaria sp. (the most abundant bryozoan in the core), Smittina cervicornis, and Myriapora truncata, dominate while seagrass habitats vanish (40 m of water depth represent the lower range limit of P. oceanica; Gobert et al. 2006). This trend is also reflected in the mollusk record where seagrass-associated species decline. Cellaria meadows of this kind form a characteristic facies of circalittoral detritic bottoms in the northeastern Adriatic (Seneš 1988a, 1988b; UNEP/MAP-RAC/SPA and Cerrano 2015) and were described as early as 1928 by Vatova (1928). In the vicinity of Rovinj, Croatia, such meadows occur at ∼ 35 m depth on detritic bottoms influenced by currents (McKinney and Jaklin 2000; Novosel 2005). By means of tubular rhizoids, Cellaria sp. is able to attach not only on small hard substrates but also directly in the sediment. Thus, it is well adapted to seafloor siltation, to which it can contribute by acting as a sediment trap (McKinney and Jaklin 2001; Hayward and McKinney 2002). While the overall abundance of mollusks decreased during the eHST, their diversity reached a peak, with many epifaunal and several infaunal species. This diverse, bryozoan-dominated benthos, at Brijuni, persisted relatively unchanged into the lHST (uppermost 20 cm).

At the Piran stations, water energy increased and sedimentation rate declined through time. The shell bed that was established during the MFZ continued to remain active during the highstand phase. This persistence is in contrast to Brijuni station where oceanographic changes (primarily increasing water depth and resulting sediment siltation and water turbidity) were associated with a decline in water energy and an increase in sedimentation rate, thus driving transgressive to early highstand community shifts that contributed to the loss of shell-beds. At Piran, transgressive sediments have a high proportion of silt, which decreases towards the surface in favor of sand and gravel. The coarser fraction mainly consists of shell fragments whose size increases up-core. With the establishment of the modern circulation system in the northern Adriatic at the end of the TST ∼ 7,000–6,000 years ago (Asioli et al. 1996), stronger bottom currents led to an increased winnowing of finer particles, to reduced sedimentation rates and thus to stabilization of the seafloor favoring colonization by epifauna (Mautner et al. 2018). Once established, epifauna-dominated communities were enhanced by positive feedback effects such as (1) the positive dependence of larval survivorship on the abundance of adult conspecifics and (2) taphonomic feedback processes by which dead hard parts of epifauna facilitate the settlement of other epifauna and inhibit burrowing activities of infauna at the same time (Kidwell and Jablonski 1983; Mautner et al. 2018). These developments culminate with the establishment of extensive banks with Arca and Modiolus between 4,000 and 2,000 years ago. In the cores, this community is represented as a massive shell bed from 35 to 8 cm core depth. Only at Piran 1 station is a similar but less condensed assemblage of large shells of epifaunal mollusks found at the transition between the early and late TST (95–45 cm). These densely colonized and highly diverse Arca/oyster bottoms persisted over large areas of the northern Adriatic Sea until the beginning of the twentieth century (Mancini 1929) before human impacts lead to their extirpation.

In all cores, signatures of pronounced shifts in the quality of the marine environment and in the benthic communities occur in the uppermost 10–20 cm. These shifts can be directly or indirectly related to human activities. The signatures of shifts are represented by (1) a marked decrease in density and effective species richness at Brijuni and both Piran stations, but not at Venice, and (2) a drop in the effective number of species, except for the Piran 1 station. On the functional level, these changes lead to infauna-dominated communities and in terms of feeding guilds to an increase of suspension feeders at the expense of grazers (Brijuni, Piran) and deposit-feeders (Venice). This convergent and uniform functional trend is all the more remarkable given the strong differences among stations in terms of their sediment composition and water depth.

Human activities contributed in a number of ways to the alteration of the benthic baseline community at the Venice station: (1) Early manipulations of the sedimentation regime of northern Adriatic rivers (e.g., Piave, Adige) by means of river alterations during Roman times and the subsequent collapse of these structures in the early Middle Ages, which led to enhanced sediment input into the basin and to increased water turbidity—a stressor for the growth of deeper seagrass meadows (Amorosi et al. 2008; Carton et al. 2009; Fontana et al. 2014); (2) Increasing eutrophication of the northern Adriatic Sea during the nineteenth and twentieth centuries negatively impacted seagrass growth not only by reducing water transparency but also by enhancing epiphytic leaf overgrowth and thus hampering photosynthesis (Justić 1988, 1991; Rismondo et al. 1997; Gobert et al. 2006; Gilad et al. 2018); and (3) Intensive bottom trawling, which was practiced in this area as early as in the eighteenth century by the Italian fishing fleets from the city of Chioggia (“chioggioti”), became a decisive factor in the alteration of northern Adriatic soft bottoms after the industrialization of fisheries in the twentieth century and an exponential increase of fishing effort (Krisch 1900; Fortibuoni 2009; Fortibuoni et al. 2010). Pollution seems to play a minor role at this station. Although RDA selects PCBs as a proxy for up-core community change, maximum PCB values remain below Threshold Effects Levels, and PCB co-varies positively with concentrations of nitrogen and carbon, indicating that it correlates with a general eutrophication trend.

The highest pollution values of all four sampling stations including high concentrations of several heavy metals (Hg, Zn, Pb) and organic pollutants (PAHs and PCBs) were observed in the lHST units of the Brijuni core (Fig. 10; see also Schnedl et al. 2018). The most probable sources are the industrial port of Pula located very close to the sampling station, and the use of the Brijuni islands area as a naval port during the late nineteenth and most of the twentieth century (Mlakar 1971; Fabianovic and Matijasec 2008). In addition to elevated contaminant levels, rising concentrations of total nitrogen in the uppermost core increments reflect ongoing eutrophication in this area mainly due to wastewater and agricultural discharge (Howarth 2008). Eutrophication and heavy metal contamination are highlighted in the RDA as the main factors influencing community change in surficial sediments. Furthermore, sediment siltation increases in the uppermost increments in spite of the absence of major sea-level fluctuations and can thus be attributed to human influences. Two factors are likely to contribute to siltation: (1) An increase in soil erosion due to extensive deforestation during the Roman period and in the late Middle Ages (Beug 1977; Oldfield et al. 2003; Kaligaric et al. 2006; Kranjc 2009), and (2) Excavation and construction activities in the nearby harbor of Pula and associated dispersal of sediment plumes by tidal currents. While bottom trawling is a widespread and intensely practiced fishing method in the northern Adriatic Sea, the Brijuni sampling station was protected from this disturbance at least over the last 30 years due to its location within the boundaries of the Brijuni Islands National Park, which was established in 1983 (Fatovic-Ferencic 2006). Nevertheless, the recent community shifts at this station, with rapidly declining density and diversity, a loss of grazers and increase of infaunal species in the uppermost 20 cm, illustrate the strong effects of cumulative human pressures on soft-bottom benthic assemblages and their difficulty in recovering despite decades of protection (Guarnieri et al. 2016).

While the concentration of several heavy metals (e.g., Hg and Zn) and of PAHs are high in the lHST unit of the Piran stations and may negatively affect some benthic biota (Mautner et al. 2018), the crucial factors in disrupting the extensive Arca beds were most likely the growing pressures from bottom trawling and frequent hypoxic/anoxic events in the twentieth century (Gallmetzer et al. 2017; Tomašových et al. 2017). On fishing maps for the northern Adriatic Sea from the beginning of the twentieth century, wide areas are demarcated where the seafloor was still covered by “banchi di mussoli” (Italian for “Arca noae banks”) and which are described as danger zones for bottom trawlers due to the damage this seafloor cover could cause to fishing gear (Mancini 1929). Both Piran stations fall within this area.

However, with the modernization of fisheries between the two world wars (boats were increasingly equipped with motors, and new, larger and more destructive trawls were employed), the Arca beds were soon completely eradicated, and by the 1950s the last Arca communities had vanished from northern Adriatic soft bottoms (Fortibuoni 2009; Lotze et al. 2011; Mautner et al. 2018). A compilation of catch data for Arca noae in the northern Adriatic illustrates the population crash around 1950 and the dramatic resource exploitation in the preceding years (Fig. 14A). Populations of this species in the Adriatic Sea are characterized by slow growth rate and are thus sensitive to over-exploitation (Peharda et al. 2002; Peharda et al. 2006). Similar scenarios of industrialized mollusk harvesting leading to the destruction of molluscan biostromes (particularly oyster banks) and to a profound alteration of extensive benthic ecosystems appear to have occurred worldwide in the late nineteenth and early twentieth century (e.g., Edgar and Samson 2004; Thurstan et al. 2013).

An additional stressor for the northern Adriatic benthic communities in the second half of the twentieth century are hypoxic or anoxic events whose frequency increased drastically towards the 80s and 90s leading to die-offs of benthic organisms over wide areas (Justić 1991; Stachowitsch 1991; Riedel et al. 2012; Tomašových et al. 2017). Such repeated events keep benthic communities permanently at an early stage of recovery with moderate diversity and dominance by small-sized individuals (Gallmetzer et al. 2017). Meanwhile, infaunal and opportunistic, resistant species such as C. gibba are promoted (Occhipinti-Ambrogi et al. 2005; Mavric et al. 2010; Nerlović et al. 2011; Tomašových et al. 2018), while epifauna persists predominantly in the form of isolated, patchy aggregations around fast-growing species (e.g., ascidians and sponges) (Riedel et al. 2012; Tomašových et al. 2017).

The three locations in the northern Adriatic Sea markedly differ in molluscan community composition and in important habitat-forming environmental drivers, in particular sediment composition and grain size. Because all cores reach back to the marine transgression at the beginning of the Holocene, the sedimentological and geochemical composition of each sediment column epitomizes a site-specific, long-term interaction of multiple shifting environmental factors and ecological feedbacks. These range from shifting bottom current patterns, increasing water depth, changing sediment and nutrient input, to recent anthropogenic disturbances such as pollution and sea-floor perturbation by fishing gear. Our results suggest that habitat differences between sites might have been even more pronounced in earlier phases of the Holocene before human impact peaked in modern times.

In the three sampling areas—Venice, Brijuni, and Piran—different benthic communities established in accord with the prevailing sedimentological and oceanographic conditions. These communities emerged during the eHST (Venice, Brijuni) or the MFZ (Piran) after sea level and circulation had reached their modern character, and persisted over several hundred to thousands of years.

At some sampling stations, in the course of the Holocene, very similar communities became established, persisted for a limited time, and then vanished following diverging environmental developments. For example, during the MFZ, at both the Brijuni and Piran stations, epifaunal communities were established probably due to reduced sedimentation rate coupled with high availability of large dead bivalve shells on the seafloor (preserved in the core as a shell bed). At Piran, these assemblages were at the beginning of a succession leading to the formation of extensive beds with Arca and Modiolus, which represent the baseline community at this station. At Brijuni, however, different oceanographic conditions and increasing siltation favored the development of epifaunal assemblages dominated by bryozoans.

The anthropogenic impacts causing decline or downfall of baseline communities are either historical or very recent, depending on the sampling stations. For example, human-induced siltation processes caused by soil erosion due to deforestation and land use (as visible from the Brijuni core) can reach back several millennia, whereas the impact of bottom trawling and hypoxia became crucial only in the twentieth century. Recent and direct impacts such as bottom trawling can quickly alter benthic communities already affected by long-term stressors such as increasing siltation or eutrophication (probably over a few years or decades, as at the Piran stations).

Long-term community changes over the Holocene were primarily driven by gradual changes in sediment composition and/or sea level. Anthropogenic disturbances affected the benthos either by repeated severe impacts (fishing, hypoxic or anoxic events) or gradual changes (eutrophication, pollution). These changes also led to a marked shift of the physical and structural parameters of the substrate, as is the case for bottom trawling. By removing the living cover from the seafloor (e.g., seagrass meadows at Venice station, epifaunal shelly assemblages at Piran), dredging reduces natural filtering capacity and enhances resuspension processes. Thus, benthic communities suffer the double impact of direct damage and of gradual changes of sediment texture and sedimentation regime.

At all stations, regardless of differences in sediment composition and water depth, the human-induced alteration of baseline communities involves a shift towards infauna-dominated communities and a decline of epifaunal forms. The trend is paralleled by an increase of suspension-feeding species at the expense of grazers (Brijuni, Piran) and deposit-feeders (Venice). This profound ecological change reduced species diversity and affected life habits indicating that modern soft-bottom benthic communities in the northern Adriatic Sea no longer show the high geographic variability of the pre-anthropogenic benthos.

Many thanks to the captain of the coring vessel, Jernej Sedmak, for his tireless commitment during the coring campaign in summer 2013. We thank Christa Hermann and Assistant Professor Dr. Robert Peticzka (Laboratory of Physical Geography, University of Vienna) for their help with grain size analysis, and Graham Oliver for his assistance with the identification of bivalves. We thank Thomas Olszewski and two anonymous reviewers for their helpful comments on the manuscript. Adam Tomašových was supported by the Slovak Research and Development Agency (APVV 17-0155) and by the Slovak Scientific Grant Agency (VEGA 0169/19). This work was funded by the Austrian Science Fund (FWF): P24901.

Data are available from the PALAIOS Data Archive: http://www.sepm.org/pages.aspx?pageid=332.