Abstract
The effects of and the interplay between natural and anthropogenic influences on the composition of benthic communities over long time spans are poorly understood. Based on a 160-cm-long sediment core collected at 44 m water depth in the NE Adriatic Sea (Brijuni Islands, Croatia), we document changes in molluscan communities since the Holocene transgression ~11,000 years ago and assess how they were shaped by environmental changes. We find that (1) a transgressive lag deposit with a mixture of terrestrial and marine species contains abundant seagrass-associated gastropods and epifaunal suspension-feeding bivalves, (2) the maximum-flooding phase captures the establishment of epifaunal bivalve-dominated biostromes in the photic zone, and (3) the highstand phase is characterized by increasing infaunal suspension feeders and declining seagrass-dwellers in bryozoan-molluscan muddy sands. Changes in the community composition between the transgressive and the highstand phase can be explained by rising sea level, reduced light penetration, and increase in turbidity, as documented by the gradual up-core shift from coarse molluscan skeletal gravel with seagrass-associated molluscs to bryozoan sandy muds. In the uppermost 20 cm (median age <200 years), however, epifaunal and grazing species decline and deposit-feeding and chemosymbiotic species increase in abundance. These changes concur with rising concentrations of nitrogen and organic pollutants due to the impact of eutrophication, pollution, and trawling in the 20th century. The late highstand benthic assemblages with abundant bryozoans, high molluscan diversity, and abundance of soft-bottom epi- and infaunal filter feeders and herbivores represent the circalittoral baseline community largely unaffected by anthropogenic impacts.
Keywords
Introduction
Investigations on modern marine ecosystem responses to natural and anthropogenic changes are typically limited to annual or decadal time scales. However, marine ecosystems have experienced major human impacts over the past few hundred years (Jackson et al., 2001; Lotze et al., 2006), and consequences of such impacts are difficult to capture using ecological data, which are limited to a few years or decades (Kidwell, 2015). Death assemblages (DAs) of shelly fauna preserved on the seafloor or in sediment cores provide a window to the past (Edgar and Samson, 2004; Gallmetzer et al., 2017). They allow us to reconstruct community composition before and after the onset of human interference, and to assess community responses both to long-term natural environmental changes and human impacts (Gilad et al., 2018; Mautner et al., 2018; Rick et al., 2016). Therefore, the investigation of DAs can identify baseline community states that can serve as a reference point for management and restoration of sensitive, endangered marine areas (Dietl et al., 2015; Jackson et al., 2001; Kidwell, 2013; Rodriguez et al., 2001).
The northern Adriatic Sea is a prime example of an ecosystem strongly affected by centuries of human pressure (Barmawidjaja et al., 1995; Gallmetzer et al., 2017; Huntley and Scarponi, 2015; Kowalewski et al., 2015; Lotze et al., 2011; Mautner et al., 2018; Scarponi et al., 2017a). This semi-enclosed shallow basin with high riverine input, extensive soft bottoms, susceptibility to seasonal water-column stratification, and propensity to eutrophication due to high primary production in its western part is classified as particularly sensitive (Stachowitsch, 1984) and is considered one of the most degraded marine ecosystems worldwide (Lotze et al., 2006). Benthic mortalities caused by oxygen depletion and mucilage events have occurred here periodically for centuries (Crema et al., 1991; Danovaro et al., 2009; Tomašových et al., 2017) but became more frequent and widespread during the last 40 years (Degobbis et al., 2000; Faganeli et al., 1985; Justić, 1991; Stachowitsch, 1991). The increase in the frequency of hypoxia was promoted by anthropogenic eutrophication through agriculture and wastewater disposal (Barmawidjaja et al., 1995; Mavrič et al., 2010; Nerlović et al., 2011; N’siala et al., 2008; Vidović et al., 2016). Overfishing and pollution negatively affected ecosystem functioning, leading to changes in food-web structure and to reduced efficiency of bentho-pelagic coupling (Barausse et al., 2009; Coll et al., 2007). Exploitation by fishermen started early in history in the Neolithic age approx. 4000 BC (European Commission, 2017) as fish, mussels and shellfish were popular trading goods and important food sources for early settlers (Forenbaher and Miracle, 2005; Horvat, 2015; Jackson et al., 2001). Harbour, ship traffic and insufficient waste water management are the main causes for pollutants in the northern Adriatic Sea (Howarth, 2008; Pearson and Rosenberg, 1978). However, in spite of extensive paleoecological research in the northern Adriatic Sea, the Holocene history of shallow-shelf benthic ecosystems off the Istrian peninsula along the NE Adriatic shelf, which are associated with very low sedimentation rate, remains poorly known (Mautner et al., 2018).
In this study, we use molluscan DAs from a sediment core collected at 44 m water depth to reconstruct the history of the north-eastern Adriatic macrobenthic fauna off southern Istria since the Holocene transgression (~11,000 years ago). First, we hypothesize that large-scale environmental processes (i.e. sea level rise and changes in sediment input and turbidity) lead to long-term changes in species composition and diversity of molluscan assemblages, as well as in the sedimentary record of pollutants and organic matter. Second, we assess whether human impacts can be detected in the stratigraphic record, and suggest that the most rapid decrease in diversity and total molluscan abundance occurs in the youngest core layers that correspond to the past few centuries.
Material and methods
Study area
The northern Adriatic was flooded ~11,000 years ago at sites presently located at ~40 m (assuming subsidence of ~0.8 mm/year, Antonioli et al., 2009) and ~10,000 years ago at sites located at 30–35 m water depth (Antonioli et al., 2007; Storms et al., 2008; Vacchi et al., 2016), and maximum ingression was attained in the NW Adriatic ~7000–6000 years ago (Amorosi et al., 2016, 2017). Its surface circulation is primarily thermohaline and cyclonic (counterclockwise), with the south-eastward Western Adriatic Current in large part driven by the Po River, the largest river in the area and the main source of freshwater, sediments, and nutrients (McKinney, 2007). The productivity in this shallow sea (average water depth 35 m) is among the highest in the generally oligotrophic Mediterranean (Zavatarelli et al., 1998). The basin is characterized by a west-eastern gradient in nutrients (from eutrophic to oligotrophic) and sediment composition (from fine-grained muds to coarse-grained bioclastic sediments) (Zuschin and Stachowitsch, 2009). Holocene sediment cores from the north-eastern Adriatic show a succession from fluvial to brackish and marine facies (Cattaneo and Trincardi, 1999; Orogelec et al., 1991). High concentrations of Hg, Cu and Zn characterize modern deposits that accumulated during the past thousand years in the Gulf of Trieste (Covelli et al., 2006; Mautner et al., 2018; Vidović et al., 2016). Temporary sediment storage and resuspension mainly occur during the winter storm period caused by north-easterly cold winds (bora), but also during strong summer winds from southeast (sirocco) (Frignani et al., 2005). Sediments are resuspended by bottom trawling and dredging (Kaiser et al., 2002; Thrush and Dayton, 2002), which occurred in the Adriatic Sea for more than a century (De Madron et al., 2005; Palanques et al., 2001). Moreover, due to water column stratification from summer to mid-autumn, short-term hypoxic and anoxic events, frequently coupled with mucilage events, causing mass mortalities of benthic communities occur frequently in the northern Adriatic Sea (Crema et al., 1991; Stachowitsch, 1991; Tomašových et al., 2017). These events tend to occur at almost annual frequency in the NW Adriatic off the Po Delta (Alvisi and Cozzi, 2016). However, they also affect the Gulf of Trieste and NE Adriatic off Istria, although with smaller frequency and less regularly (Djakovac et al., 2015; Harding et al., 1999).
Sampling and sediment analysis
Sampling took place in the waters of the Brijuni Islands National park, which is situated close to the SW coast of the Istrian peninsula in Croatia (Figure 1). The archipelago is spread across a total area of 7.42 km2 with two major islands, Mali Brijun and Veliki Brijun (Fatovic-Ferencic, 2006; Šoštarić and Küster, 2001). The first major anthropogenic impact on these islands occurred in the late 19th century when swamps were dewatered to promote tourism (Bralic, 1990; Šoštarić and Küster, 2001). In the 20th century, the area became a presidential residence and was declared a national park in 1983. Thus, its land and waters remained largely excluded from extensive commercial and public use in the late 20th century (Fatovic-Ferencic, 2006).

Map of the Brijuni Islands, Croatia, and exact location of the sampling station within the marine protected area of the national park Brijuni (NP Brijuni). The grey dotted line marks the boundaries of the national park.
In summer 2013, at a water depth of 44 m (coordinates: 44.8857667 N and 13.747E), four sediment cores (length 160 cm, diameter 9 cm) were extracted with the Slovenian research vessel
Dating of shells
Core chronology was based on the 14C-calibrated amino acid racemization (AAR) of 305 specimens of the bivalve
In total, 17 bivalves and the single plant remnant were analysed for 14C at the AMS facility at the Poznan Radiocarbon Laboratory (Table 1). 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% H2O2 for 10 min in an ultrasonic bath. The remaining carbonate was dissolved with concentrated H3PO4 in a vacuum line. 14C was measured with a ‘Compact Carbon AMS’ (Goslar et al., 2004). Conventional 14C ages were calculated using correction for isotopic fractionation (Stuiver and Polach, 1977), on the basis of ratio 13C/12C measured in the AMS spectrometer simultaneously with the ratio 14C/12C. Conventional ages of bivalves were converted to calendar years using Calib7.1 (Stuiver and Reimer, 1993), the Marine13 calibration curve (Reimer et al., 2013), and a regional marine reservoir correction (∆R) in the NE Adriatic (Rovinj) equal to = -61 years (standard deviation = 50 years) (Siani et al., 2000). Conventional ages of the plant remnant were converted to calendar years with the IntCal13 calibration curve.
Numerical ages of 13 dead-collected and three live-collected specimens of
Small portions of 305 dead shells of
We calibrated the rate of AAR in

Relations between Asp D/L and 14C calibrated shell ages for
Calibration statistics for the rate of amino acid racemization (AAR) based on paired AAR and radiometric analyses of
Sorting, identification, and counting of shells
For the investigation of the molluscan shell record, the sediment samples were washed through a 1-mm sieve and dried. The retained material was sorted under a stereomicroscope, and molluscs were identified to species level using taxonomic literature (Cossignani and Ardovini, 2011; Cossignani et al., 1992; Gofas et al., 2011a, 2011b; Huber, 2010). The final number of bivalve individuals was computed by summing up the number of double-valved specimens with the higher number of either right or left valves. Gastropod shells were counted as individuals if at least the apex was preserved, and polyplacophorans were counted as the number of plates divided by eight.
Functional groups
Molluscan species were categorized according to (1) feeding guild, (2) organism-substrate relation, (3) host association, and (4) weed association (Table 3). The categories and subcategories were treated separately in the analyses. Categorization follows Cachia et al. (1996), Chemello and Milazzo (2002), Morton and Dinesen (2011), Oporto et al. (2012), Rinaldi (2012), Rueda et al. (2009), Werner (1953, 1959), Beesley et al. (1998), Borja et al. (2000), Gofas et al. (2011a, 2011b), Haydar (2010a, 2010b), Huber (2010), Koulouri et al. (2006), and Oliver et al. (2016).
Functional categorization used to characterize the ecology of mollusc species. Subcategories may not appear in figures but are listed here to illustrate the spread of the main categories.
Statistical analysis
Species diversity is expressed as raw number of species and as the effective species richness, that is, the exponential of Shannon’s H (exp (H)). All values reported are means with 95% confidence intervals (1000 resampling iterations, randomly drawn without replacement). Differences in diversity between the six stratigraphic units defined on the basis of lithologic and geochemical composition were tested by a Wilcoxon rank test. PERMANOVA was used to assess differences in community structure between the six stratigraphic units on the basis of species composition and on the basis of functional composition related to (1) feeding guilds, (2) organism-substrate relation, (3) weed-, and (4) host-association. The uppermost unit (0–8 cm, last 200 years) was excluded from analyses due to insufficient sample size (minimum N = 50). Sample size of other molluscan assemblages varies between 133 and 1515 specimens. To analyse gradients in community composition, we performed a non-metric multidimensional scaling (NMDS) based on Bray-Curtis distances and square-root-transformed relative abundances. The influence of environmental variables on the composition of molluscan assemblages was explored by redundancy analysis (RDA) using a forward-model-selection approach. Clay, gravel, water depth, total nitrogen, as well as the heavy metals Hg, Pb, As, and Cd were used for the analysis. The variable with the highest adjusted R2 is added at each step of the model selection and its significance calculated using a permutation test (10,000 permutations) (Blanchet et al., 2008). This allows identifying the subset of environmental variables with the strongest influence. All statistical analyses were performed in R 3.2.1 (R Core Team, 2014) using the vegan package (Oksanen et al., 2015).
Results
Stratigraphic section and geochronology
The lowermost part of the core M40 directly overlies an organic-rich clay rich in plant debris, indicating the transition from terrestrial to marine conditions. The age of the plant remains (Table 1) is in accord with sea-level curves estimated for this region (Antonioli et al., 2007), indicating that the flooding occurred about 11,000 years ago. We distinguish six stratigraphic units that differ in sediment grain size, sediment and geochemical composition, and preservation (in terms of dominant components such as molluscs, bryozoans, and coralline algae). They include unit 1 formed by coarse-grained, loosely packed and moderately sorted molluscan sands at 120–160 cm, unit 2 formed by densely packed and poorly sorted molluscan gravel (shell bed) at 90–120 cm, sharply delimited by an abrupt increase in abundance of large molluscs on the bottom and by an abrupt decline in their abundance on the top, unit 3 formed by poorly sorted and loosely packed bryozoan-molluscan muddy sand at 50–90 cm, unit 4 formed by poorly sorted and loosely packed bryozoan-molluscan sandy mud at 20–50 cm, unit 5 formed by sandy mud with dispersed bioclasts at 8–20 cm with high concentrations of pollutants (PCB and PAH), and unit 6 formed by sandy mud with coralline algae and dispersed larger shells with high concentrations of total nitrogen at 0–8 cm.
In terms of the 25th and 75th age percentiles, shell ages are younger than 370 years in unit 6, range between 200–1700 years in unit 5, 1500–3600 years in unit 4, 2900–4600 years in unit 3, 4500–7500 years in unit 2, and 4800–7500 years in unit 1. Median ages of

Stratigraphic distribution of shell ages of
Raw inter-quartile age ranges of
General trends in sediment composition
At the bottom of the core, dark grey, organic-rich clays with frequent plant remains mark the transgression from a paralic to a marine environment (Figure 4). Unit 1 is formed by molluscan muddy sand with a mixture of stained and pristine bioclasts. Higher up, the skeletal gravel and sand fraction increase steeply and reach a maximum of ~80% in unit 2 at 100 cm. This peak corresponds to a ~25-cm-thick, densely packed and poorly sorted bivalve shell bed with high abundance of encrusting organisms. The units 3–6 contain abundant bryozoans and molluscs and show a fining-upward trend, leading to a high proportion of clay and silt (~80%) in surface sediments (Figure 5). The heavy metals Hg, Pb and Zn attain maxima in the upper 20 cm, while Cd peaks at 50 cm at the boundary between units 3 and 4. Only Hg exceeds NOAA ERM levels during unit 2. The trend of Cr, Cu and Ni is closely correlated with the fine sediment fraction (silt and clay), while no correlation is evident for the other metals. Fe and Al show opposite trends: Fe concentration reaches 19.5 g/kg in the deepest subsample (150–152 cm), while Al concentration peaks at a depth of 23–25 cm (17.6 g/kg). High concentrations of polycyclic aromatic hydrocarbons (PAHs) occur in the uppermost 12 cm, with a maximum of 475 ng/gr at 5 cm. Polychlorinated biphenyls (PCBs) also reach their maximum (6 ng/gr) at 5 cm core depth. In the rest of the core, their concentrations are very low, in the case of PAHs close to zero. TOC and total nitrogen increase up-core (Figure 5).

(a) Paralic layer from the lowermost core increment originating from the period before the Holocene transgression, (b) bryozoan-molluscan gravel dominating the upper core parts until 95 cm depth, (c) muddy sediment surface at the sampling station with two bryozoan colonies covered by mucilage, and (d) pectinid and oyster shells from the macroid layer overgrown by serpulid tubeworms and laminar bryozoans and heavily bioeroded by boring sponges; photographs by Ivo Gallmetzer and Alexandra Haselmair.

Concentrations of heavy metals and organic pollutants along the core. Hg, Pb, Cd, and Zn in mg/kg; PAH and PCB in ng/g; and TOC and TN in percent dry weight. Dotted lines boundaries of age units. Hg concentrations exceed NOAA ERM benchmark (grey line).
Molluscan assemblages
Throughout the core, 23,012 individuals representing 229 mollusc species from 106 families were identified, including 24 non-marine gastropod species, which are very abundant in the deeper parts until 70 cm depth. Gastropods account for 149, bivalves for 69, polyplacophorans for 8, and scaphopods for 3 species. Gastropods are more abundant than bivalves in most of the core, but are absent in the uppermost 8 cm (last 200 years, Figure 6). Total abundance of individuals peaks in unit 1 at 130 cm (>5500 years ago) and strongly decreases up-core. Diversity changes are best illustrated by the effective number of species, which peaks between 35 and 90 cm core depth (2000–4500 years ago) (Figure 6). Apart from molluscan shells, the sediment core was also rich in calcareous skeletal remains, mostly of bryozoans, echinoderms, and crustaceans (Pifeas, 2017). Layer 95–120 cm was characterized by a high abundance of heavily bioeroded shells with thick calcareous encrustations, mainly by bryozoans. The six units differ significantly in species composition and abundance (Figure 7 and Table 4; unit 6 was excluded due to small sample size).

Abundance and diversity of molluscan taxa in the sediment core along a depth (left) and time axis (right). The unit boundaries marked by dotted lines. (a) total abundance, (b) raw number of species, and (c) effective number of species (exp (H)).

Relative abundances of the 12 most abundant species in each unit.
Differences in species composition and functional groups between age units (tested with PERMANOVA, except for exp(H), where a Wilcoxon-rank test was applied). Significant results in bold. The 0–8 cm unit was not tested because of an insufficient number of individuals.
Throughout the core, epi- and infaunal life habits dominate, and filter-feeders are the most important feeding guild. Epibionts and ectoparasites are more abundant in deeper sediment layers. The relative (and absolute) abundance of detritivores increases up-core while the abundance of seagrass-associated species decreases (Figure 8). The bivalves

Up-core shifts in the dominant substrate relations, feeding types, and weed associations since the Holocene transgression, plotted for relative (a-c) and absolute abundances (d-f). Dotted lines mark unit boundaries.

(a) NMDS plot of all determined marine mollusc species and (b) redundancy analysis, with Pb and total nitrogen (TN) identified as the main environmental drivers of molluscan community change.
Unit 1 (120–160 cm, 4800–10,500 years)
A relatively rapid shift in grain size occurred in this unit, from 80% clay and silt at 165 cm to 80% sand (mostly consisting of molluscs), at 120 cm (Figure 5). Total abundance of mollusc shells is low in the deepest layers, but increases steadily towards the upper unit boundary. This unit contains the sample with the highest abundance (1557 specimens at 125–130 cm) of molluscs in the entire core. The raw number of species shows an up-core increase, whereas the effective number of species remains more or less stable (Figure 6). The most dominant species are
Shell bed (90–120 cm, 4500–7500 years)
This unit corresponds to a marked shell bed with large-sized, randomly oriented bivalves. Grain size composition shifts from 80% sand and gravel to 50% silt and clay. Total abundance and species richness peak in this unit (Figure 6).
Unit 3 (50–90 cm, 2900–4600 years)
The sediment composition in this unit formed primarily by bryozoan muddy sand shifts towards higher proportions of silt and clay (60%), accompanied by a decrease in total molluscan abundance. Diversity (exp(H)) remains high (Figure 6) but is still significantly lower than in unit 2 (Wilcoxon rank test: p = 0.02). Although
Unit 4 (20–50 cm, 1500–3600 years)
The fining-upwards trend in sediment grain size (from 21% to 27% of clay) persists within this unit formed by bryozoan-molluscan sandy mud and characterized by a marked community shift in both species composition and functional groups. Total abundance and diversity (exp(H)) strongly decrease (Wilcoxon rank test: p = 0.02). This change occurs abruptly at about 30 cm (~2000 years ago, Figure 6). Concomitantly, epifauna proportionally decreases in favour of infauna, and relative abundances of detritus feeders and grazers rise (Figure 8). The most abundant species is
Unit 5 (8–20 cm, 200–1700 years)
Clay and silt constitute 78% of the sediment in this unit formed by rhodolithic-molluscan-bryozoan-sandy mud with higher concentrations of pollutants. Total abundances of molluscs are low in comparison to former units. Species richness and effective number of species reverse and increase again (Figure 6). Infaunal bivalves and epifaunal gastropods characterize this unit.
Unit 6 (0–8 cm, < 370 years)
The sediment in this unit, whose assemblage is strongly segregated from the rest of the core in ordination space, is formed by rhodolithic–molluscan–bryozoan sandy mud with organic enrichment and consists mainly of pelitic material (Figure 5). Abundances and diversity (exp(H)) drop to minimum values (Figure 6). Soft bottom filter-feeding and detritivore bivalves are dominant.
Relation between environmental parameters and community change
The RDA shows that Pb (F = 5.424, p = 0.002) and TN (F = 2.1638, p = 0.014) have highly significant effects on the variation in the species composition of molluscan assemblages. However, clay content, Hg, Zn, PAHs and PCBs significantly correlate with concentrations of Pb and TN. Therefore, effects of these variables can be masked by the effects of Pb and TN (Figure 9b).
Discussion
Our analyses show that the sediment cores collected at 44 m water depth from the Brijuni Islands reach back to the Holocene transgression about 11,000 years ago. They capture the transgressive phase coupled with shoreface erosion and reworking, followed by the establishment of epifaunal bivalve-dominated biostromes during the late transgressive/early highstand transition (maximum-flooding zone), by the shift to low-energy bryozoan-dominated biostromes during the highstand phase, and finally by strong environmental changes triggered by human impact mainly during the last century. Although molluscan assemblages are time-averaged at multi-centennial to millennial scales (Figure 3), the changing species composition gives insight into the environmental shifts occurring over longer, Holocene time-scales. The emergence of epifaunal molluscs associated with seagrasses or algae point to a habitat with marine vegetation during the transgressional phase. During highstand, bryozoans began to dominate and shape the benthic community at the expense of epifaunal molluscs. Finally, the up-core increase of fine sediments led to a dominance of infaunal molluscs.
Transgressive phase (unit 1): Mixing of freshwater and marine assemblages during the Holocene transgression
We infer that the molluscan assemblages during this phase largely formed since the flooding ~11,000 up to 7000–8000 years ago. The Holocene boundary with the underlying plant-rich muddy layer (Figure 4a) and the high abundance of land and freshwater gastropod species in the deepest increments are indicative of the terrestrial, freshwater-dominated environments (marshlands, swamps, and estuaries) occurring on the northern Adriatic coastal plains prior to the formation of barrier-lagoon systems during the marine transgression (Bortolami et al., 1977; Correggiari et al., 1996b; Trincardi et al., 1994). The mixture of land and freshwater snails and marine molluscs in unit 1 documents the transition from a terrestrial to a shallow, nearshore to upper shoreface environments on a sandy, high-energy seafloor (Figure 7). Lithologic changes and paleoecological and taphonomic attributes of assemblages in unit 1 in the lowermost core part are closely related to the flooding, initial deposition in nearshore barrier-lagoon system (comparable to similar systems located at 42 m water depth and flooded ~10,500 years ago, Storms et al., 2008) and subsequent reworking and cannibalization of nearshore deposits due to high-energy shoreface erosion. These processes led to a mixture of freshwater, brackish, and marine species in molluscan gravels preserved in unit 1 (i.e. transgressive lags or transgressive sand sheet, Cattaneo and Steel, 2003). The transgressive lag deposits consisting of molluscan shell hash, frequently with worn and stained preservation, indicates long residence times and repeated reworking (as indicated by mixing of freshwater and marine species and by interquartile age range of
Maximum-flooding zone (unit 2): Establishment of diverse epifaunal biostromes
The sediment is predominantly formed of gravel and sand-sized shells and shell fragments of large-sized scallops, glycymeridids, arcids, and oysters. The effective number of species peaks in this unit (Figure 6) indicating the establishment of a highly diverse marine fauna. In the late phase of sea level rise, the habitat changed from macrophyte-dominated towards biogenic hard bottoms either in the form of extensive bivalve biostromes (similar to some degree to an extinct, oyster-dominated shell bed documented by Mautner et al. (2018) in the southern Gulf of Trieste) or in the form of isolated multi-species clumps as present today in some parts of the northern Adriatic (Fedra et al., 1976). High abundances of epifaunal carnivorous molluscs (Mangeliidae, Conidae, Raphitomidae) as well as parasites on sponges (Triphoridae) and sea anemones (Epitoniidae) further demonstrate the presence of habitat-forming aggregations of large suspension-feeding invertebrates. Most molluscan shells (oysters, pectinids, arcids, and glycymeridids) are heavily bored and encrusted by bryozoans (Figure 4d), indicating low sedimentation rates and substrate stability, providing ideal conditions for settlement of epibionts and encrusters. The temporal coincidence of this shell bed with the timing of fastest sea level rise ~7000–6000 years ago and high degree of encrustation and bioerosion indicates the connection between shoreline retrogradation, reduced sediment supply, and shell-bed formation, with potential for taphonomic feedback. Similar examples were described from maximum-flooding zones of cores collected in offshore locations off the Po Delta (Scarponi et al., 2017b). However, grazing species (
Early highstand (unit 3): Algal-bryozoan meadows
The community in this section is characterized by an increasing importance of bryozoans, which make up a large part of the death assemblage (Figure 4b) (Pifeas, 2017), demonstrating a shift from vegetation-dominated communities and bivalve shell beds to bryozoan meadows. Together with the increasing siltation and sediment trapping, the shift in species composition towards higher abundance of infaunal suspension- and deposit-feeders and towards a decline of seagrass-related molluscs suggests that seagrass meadows vanished with increasing water depth, reduced light penetration and higher turbidity (Figure 8) (Abal and Dennison, 1996; Duarte, 1991). However, species associated with algae (
Early highstand (unit 4): Siltation and the growing importance of deposit feeders
The drop in mollusc abundance and the relative increase in abundance of infaunal deposit-feeders and infaunal suspension-feeders such as
Early and late anthropogenic impact (units 5 and 6): Pollution and eutrophication
The concentrations of pollutants and nitrogen content in these two units trace a growing anthropogenic impact on the sea, deriving from expanding coastal settlements, land use, industrialization, and concomitant pollution and eutrophication (Lotze et al., 2006). This impact is documented by the abrupt and steep increase in the concentrations of heavy metals (Hg, Pb, and Zn) and organic pollutants (PAHs and PCBs) (Figure 5). The most probable sources of this pollution are the adjacent port of Pula with its dockyard, and the use of the area as a naval port during the late 19th and most of the 20th century (Fabijanović Matijašec, 1971, 2008). The presence of these elements and organic compounds, which are unequivocally of industrial origin (Covelli et al., 2006; Rios et al., 2007; Vane et al., 2007), in sediment layers that are several centuries old is probably a consequence of sediment mixing. However, all mentioned pollutants recede again in the uppermost 5–10 cm of sediment, illustrating a trend-reversal after the late-20th century efforts to reduce environmental pollution (Lotze et al., 2011). In contrast to heavy metals and organic pollutants, the concentration of total nitrogen increases steeply and steadily up to the surface, indicating that wastewater discharge and agriculture represent the main sources for an ongoing eutrophication process in this coastal area (Howarth, 2008).
Our analyses indicate that compositional variability of molluscan assemblages is primarily related to changes in nitrogen and Pb concentration, which both correlate closely with clay content (Figures 5 and 9b). While increasing sediment siltation on the transition from the transgressive to the highstand phase can be mainly seen as a consequence of deepening and changes in circulation patterns and current regimes (Artegiani et al., 1997; Boldrin et al., 2009), during the deposition of the uppermost two core units (0–1000 years) this process was further augmented by human activities such as coastal deforestation (Kranjc, 2009) and intensified bottom trawling. Although the Brijuni archipelago is protected since 1983 (Fatovic-Ferencic, 2006), this fishing method has been widely used in the northern Adriatic since the 20th century (Bolje et al., 2014) and can be responsible for high rates of sediment resuspension (Kaiser et al., 2002; Korpinen et al., 2013; Zuschin and Stachowitsch, 2009).
The selective loss of epifaunal species is reminiscent of the loss of epifaunal communities in the southern Gulf of Trieste in the 20th century as documented for a sampling station at 20 m water depth (Mautner et al., 2018). However, on one hand, the shell bed at Brijuni was naturally replaced by bryozoan meadows owing to deepening, reduced light penetration, and increased sediment trapping. On the other hand, at the deeper Brijuni station (44 m), bryozoans and other epifaunal molluscs were also severely affected by eutrophication and pollution. At times of sampling, the sea-floor was largely formed by dead skeletal debris and isolated epifaunal clumps. Bryozoan meadows, which are still preserved at some other places off Istria (McKinney and Jaklin, 2001), were degraded and not represented by living colonies (Figure 4c). Under increasing anthropogenic pressure, the molluscan community shifted towards high abundances of infaunal deposit feeders such as Nuculidae (
In contrast to the decline of suspension feeders, chemosymbiotic bivalves (mainly the lucinid
Conclusion
The core sampled at Brijuni exemplifies the environmental transition from a terrestrial to a marine soft- and mixed-bottom habitat in the south-eastern part of the northern Adriatic Sea. Natural processes caused by sea level rise during the Holocene transgression explain changes in the composition and diversity of molluscan assemblages during the deposition of the first four units. Seagrasses and macroalgae probably characterized benthic habitats during the early transgressional phase, together with the formation of biostromes with epifaunal bivalves. Seagrass habitats disappeared with increasing water depth due to sediment siltation and reduced light penetration. Bryozoan meadows developed during the highstand phase, and molluscan community was dominated by infaunal suspension-feeders and deposit feeders, mostly bivalves. Grazing gastropods, which were deprived of their main food source, decreased strongly in abundance up-core. First human settlements appeared and small-scale deforestation and agriculture during the latest highstand phase further contributed to seafloor siltation. In contrast to the earlier core units, which were defined by large-scale oceanographic events such as sea level rise and changing current patterns, the effects of pollution and eutrophication dominate in the units 5 and 6. Therefore, humans have re-shaped the environment during the past centuries with increasing speed. The late highstand benthic assemblages with abundant bryozoans, high molluscan diversity, and abundance of soft-bottom epi- and infaunal filter feeders and algae-associated species emerge as the last baseline community largely unaffected by anthropogenic impacts and could thus be a reference point for any management or conservation efforts aiming to restore pre-impact conditions.
Supplemental Material
Supplement_data-table_S01 – Supplemental material for Molluscan benthic communities at Brijuni Islands (northern Adriatic Sea) shaped by Holocene sea-level rise and recent human eutrophication and pollution
Supplemental material, Supplement_data-table_S01 for Molluscan benthic communities at Brijuni Islands (northern Adriatic Sea) shaped by Holocene sea-level rise and recent human eutrophication and pollution by Schnedl Sara-Maria, Alexandra Haselmair, Ivo Gallmetzer, Anna-Katharina Mautner, Adam Tomašových and Martin Zuschin in The Holocene
Supplemental Material
Suppl_data-table_S02 – Supplemental material for Molluscan benthic communities at Brijuni Islands (northern Adriatic Sea) shaped by Holocene sea-level rise and recent human eutrophication and pollution
Supplemental material, Suppl_data-table_S02 for Molluscan benthic communities at Brijuni Islands (northern Adriatic Sea) shaped by Holocene sea-level rise and recent human eutrophication and pollution by Schnedl Sara-Maria, Alexandra Haselmair, Ivo Gallmetzer, Anna-Katharina Mautner, Adam Tomašových and Martin Zuschin in The Holocene
Footnotes
Acknowledgements
We are grateful to Graham Oliver for sharing his profound knowledge on marine bivalves and to Thomas Neubauer for his help in identifying terrestrial gastropod shells. Lukas Landler provided helpful comments on an earlier version of this manuscript. We thank Jernej Sedmak, captain of the research and working vessel Manta Bianca, with which we conducted the coring campaign in the northern Adriatic Sea. Daniele Cassin and Roberto Zonta from the Institute of Marine Sciences in Venice, Italy, carried out the geochemical analyses. Vlasta Cosovic organized the permits and helped with fieldwork. We thank two reviewers for comments.
Funding
This study was financially supported by the project ‘Historical ecology of the Northern Adriatic Sea’ funded by the Austrian Science Fund FWF: P24901. A.T. was supported by the Slovak Scientific Grant Agency (VEGA 0136/15).
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
