Introduction
Predation is one of the key factors shaping marine ecosystems today (Paine Reference Paine1974; Connolly and Roughgarden Reference Connolly and Roughgarden1999; Freestone et al. Reference Freestone, Osman, Ruiz and Torchin2011; Harper and Peck Reference Harper and Peck2016) and through the Phanerozoic (Vermeij Reference Vermeij1977; Harper Reference Harper2003) and is hypothesized to be responsible for secular changes in the guild composition and energetics of benthic ecosystems (Bottjer and Ausich Reference Bottjer and Ausich1986; Bambach Reference Bambach1993; McRoberts Reference McRoberts2001; Aberhan et al. Reference Aberhan, Kiessling and Fürsich2006). Drill holes provide direct evidence of predation and are thus the most often studied predator–prey interaction in the fossil record (Kowalewski Reference Kowalewski2002). Several authors argued that the long-term increase in the predation intensity led to the replacement of epifaunal communities with sessile or poorly mobile species by infaunal benthic communities in soft-bottom habitats and that present-day epifaunal communities remain typical of habitats characterized by low predation pressure (Aronson et al. Reference Aronson, Blake and Oji1997; Gili et al. Reference Gili, Arntz, Palanques, Orejas, Clarke, Dayton, Isla, Teixidó, Rossi and López-González2006; Bowden et al. Reference Bowden, Schiaparelli, Clark and Rickard2011). However, it is unclear how species-level drilling frequencies (DFs) scale to community levels (Kelley and Hansen Reference Kelley and Hansen2006) and how local estimates scale up to regional or global scales (Hoffmeister and Kowalewski Reference Hoffmeister and Kowalewski2001). Given the nonrandom preservation of depositional environments in the rock record (Holland Reference Holland2000; Smith Reference Smith2007) and restricted outcrop area of many ancient successions, understanding these scaling relationships is crucial for interpreting time series of predator–prey interactions and distinguishing evolutionary factors (e.g., coevolution, escalation) from changes in the physical environment as their primary drivers (Leonard-Pingel and Jackson Reference Leonard-Pingel and Jackson2016).
First, global compilations demonstrate a distinct Phanerozoic increase in DF at the level of individual prey populations (e.g., Huntley and Kowalewski Reference Huntley and Kowalewski2007). However, local-scale dynamics may not directly extrapolate to larger spatial scales: the effects of predator–prey interactions on the food webs can propagate to entire regional ecosystems, but geographic variability in the strength of interactions and source-sink or mass-effect metacommunity dynamics can reduce their importance (Holt and Hoopes Reference Holt, Hoopes, Holyoak, Leibold and Holt2005). Second, global or regional patterns in local-scale DF are frequently perceived at the level of the total assemblage (Kowalewski et al. Reference Kowalewski, Dulai and Fürsich1998), which may reflect trends that are driven by varying relative abundances of prey species and differences in their susceptibility to drilling predation (Leighton Reference Leighton2002). However, most species or genera have constrained environmental and stratigraphic distribution, inhibiting assessments of species-level trends in DF across different habitats and over longer timescales. Aggregating species into higher taxa can thus average out uneven or patchy occurrences of species with narrow niches across samples and can reveal community-level ecological signals that cannot be perceived at the species level (Zuschin et al. Reference Zuschin, Nawrot, Harzhauser, Mandic and Tomasovych2017).
Here, we assess the congruence between local and regional (i.e., basin-scale) stratigraphic patterns in DFs at different levels of taxonomic resolution and their relationship to environmental proxies using sediment cores from the Holocene succession of the northern Adriatic Sea. This shelf exhibits a marked eastward decline in sea-surface primary productivity and sediment accumulation rate. McKinney (Reference McKinney2003) and McKinney and Hageman (Reference McKinney and Hageman2006) suggested that communities dominated by epifaunal suspension feeders in the northern Adriatic Sea are limited to oligotrophic habitats and that the decline in primary productivity is associated with a decline in predation pressure. In contrast, Zuschin and Stachowitsch (Reference Zuschin and Stachowitsch2009) showed that epifaunal communities in the northern Adriatic Sea are not restricted to oligotrophic settings, while predation intensity is similarly high as elsewhere in the Mediterranean Sea (Sawyer and Zuschin Reference Sawyer and Zuschin2010). However, as ecological surveys and DFs observed in surface death assemblages are affected by anthropogenic impacts (Smith et al. Reference Smith, Handley and Dietl2018), assessments of natural predation pressure require analyses of sediment cores that capture Holocene communities inhabiting the northern Adriatic shelf since its flooding ~10,000 years ago. To rigorously assess DFs in the northern Adriatic Sea during the Holocene, high-resolution field data need to be placed within a robust, temporally constrained framework that integrates facies relationships and stratal architecture (Brett Reference Brett1998; Leonard-Pingel and Jackson Reference Leonard-Pingel and Jackson2016; Patzkowsky Reference Patzkowsky2017). Therefore, we evaluated DF patterns in a sequence-stratigraphic framework at different taxonomic resolutions and correlated them with environmental factors, represented by grain size and geochemical proxies. This is the first study on predator–prey interaction in the fossil record that uses such an integrative approach.
We show that stratigraphic trends in DF at regional scales can be confounded by habitat differences in DFs and can be biased by geographic differences in the thickness of systems tracts. We also demonstrate that the frequencies of drilling predation in the oligotrophic parts of the northern Adriatic Sea were not particularly low during the Holocene as previously suggested (McKinney Reference F. K.2007: p. 243) and were similar to average Cenozoic levels (e.g., Baumiller et al. Reference Baumiller, Bitner and Emig2006; Klompmaker Reference Klompmaker2009; Sawyer and Zuschin Reference Sawyer and Zuschin2011; Chattopadhyay et al. Reference Chattopadhyay, Zuschin, Dominici and Sawyer2016; Neely et al. Reference Neely, Kelley and Friedman2021) (rather than to Paleozoic or early Mesozoic levels before the onset of the Mesozoic marine revolution; Kowalewski et al. Reference Kowalewski, Hoffmeister, Baumiller and Bambach2005).
One exception is represented by a location that was characterized by seagrass development, indicating that structured biogenic habitats like seagrass meadows can significantly affect predator–prey interactions (e.g., Orth et al. Reference Orth, Heck and van Montfrans1984; Leonard-Pingel and Jackson Reference Leonard-Pingel and Jackson2016).
Materials and Methods
Sediment Cores
Piston cores 136–160 cm in length, penetrating the Holocene succession and representing the local scale, were taken in 2013 at four stations with different sediment types and at water depths between 21 and 44 m (Fig. 1). At each station, a separate core was taken to evaluate stratigraphic variability in grain size and chemical composition at the Institute of Marine Sciences in Venice. Ten increments were chosen from each core, 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. Only eight samples in total were analyzed from the Venice core due to its shorter length. For each increment, grain size, total carbon, total organic carbon, total nitrogen, and carbonate content were measured from the sediment (Vidović et al. Reference Vidović, Nawrot, Gallmetzer, Haselmair, Tomašových, Stachowitsch, Ćosović and Zuschin2016; Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019). For faunal analysis, core increments, 4–5 cm in thickness, were sieved through a 1-mm mesh, and molluscan shells were counted and identified to the species level (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019).
The location of the Venice, Brijuni, Piran 1, and Piran 2 sampling stations in the northern Adriatic Sea. Coordinates of stations are provided in Gallmetzer et al. (Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019).
In total, 124 samples containing 70,783 molluscan individuals were analyzed. Two samples were excluded because their number of individuals fell below the threshold of 20 individuals, which resulted in a final dataset of 70,362 individuals. Muricid and naticid gastropods are abundant in our samples and are the most likely producers of predatory drill holes encountered (Bromley Reference Bromley1981; Kelley and Hansen Reference Kelley, Hansen, Kelley, Kowalewski and Hansen2003), but no attempt to distinguish between Oichnus paraboloides (made by naticids) and Oichnus simplex (made by muricids) was made. The mean percentage of naticids per sample is 0.9% with a maximum of 2.7%. For muricids, the mean is 0.7% and the maximum 3%. Naticidae (827 individuals) are strongly dominated by Euspira nitida (81.9%), followed by Euspira macilenta (15.6%), Natica stercusmuscarum (1.7%), and unidentified juvenile naticids (0.8%). Muricidae (495 individuals) are strongly dominated by Hexaplex trunculus (77.2%), followed by Bolinus brandaris (10.1%), unidentified juvenile muricids (6.5%), Ocenebra erinaceus (4.4%), Ocinebrina helleri (0.6%), Trophonospis sp. 1 (0.8%), and Typhinellus labiatus (0.4%). All naticids are apparently capable of shell drilling (e.g., Carriker Reference Carriker1981; Kabat Reference Kabat1990). For the muricids encountered in our samples, shell drilling has been documented for H. trunculus (Peharda and Morton Reference Peharda and Morton2006; Morton et al. Reference Morton, Peharda and Harper2007; Sawyer et al. Reference Sawyer, Zuschin, Riedel and Stachowitsch2009) and B. brandaris (Abidli et al. Reference Abidli, Vasconcelos, Lahbib, Barroso, Trigui El Menif and Gaspar2012; P. Vasconcelos personal communication August 2021). All juvenile muricids in our samples belong either to Hexaplex or Bolinus, and although they could not be assigned with confidence to either genus, they therefore still can count as drilling predators. Ocenebra erinaceus is also a facultative driller (Hancock Reference Hancock1960). We have no information on the predatory behavior of O. helleri, Trophonospis sp. 1, and T. labiatus, but together they only account for 1.8% of the total number of muricids. All naticids and virtually all muricids in our quantitative material can therefore count as obligatory or facultative shell-drilling predators.
To estimate DF, the total number of drilled shells was divided by the total number of bivalve, gastropod, and scaphopod individuals per core increment. Because only one valve in bivalves is usually drilled, DF was corrected using the equation of Bardhan et al. (Reference Bardhan, Chattopadhyay, Mondal, Das, Mallick, Roy and Chanda2012): DF = DV/[(RV + LV)/2 + A], where DV is the total number of drilled valves, RV the total number of right valves, LV the total number of left valves, and A the number of articulated individuals. As there are only 23 samples with minimum 20 scaphopod shells, we perform class-level analysis only for bivalves and gastropods. Family-level analyses are performed for the eight most abundant families (Table 1), which together account for 64% of the individuals in the total dataset, with a range from 17.4% (Veneridae) to 2.7% (Tellinidae). These families also have a sufficiently high number of samples with minimum 20 individuals. Similarly, a species-level analysis was restricted to the eight most abundant species, which together account for 43% of the individuals in the total dataset, with a range from 10.8% (Gouldia minima) to 2.7% (Striarca lactea). Note that the families Corbulidae and Noetiidae are each strongly dominated by a single species; therefore, the family- and species-level analyses are almost identical. A genus-level analysis was not performed, because most genera are either monospecific or numerically strongly dominated by one species. Freshwater species, which occur mostly at the base of the cores, were excluded from the analyses.
The number of samples and abundance of prey classes, dominant prey families and species, and predatory families per systems tract at the four stations.
Core Geochronology and Stratigraphic Framework
The core chronologies are based on extensive dating by radiocarbon-calibrated amino acid racemization of one or two bivalve species occurring at high abundance at each station, with supplementary dates from other taxa (documented in our former studies, see below). For the two Piran cores, G. minima (316 shells) and Corbula gibba (232 shells) were dated with additional dates from Arca noae and Ostrea sp. (Mautner et al. Reference Mautner, Gallmetzer, Haselmair, Schnedl, Tomašových and Zuschin2018; Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019; Tomašových et al. Reference Tomašových, Gallmetzer, Haselmair, Kaufman, Mavrič and Zuschin2019); for Brijuni, Timoclea ovata (305 shells) was used with additional dates from Glycymeris, Aequipecten, Chlamys, Ostrea, and a plant remnant (Schnedl et al. Reference Schnedl, Gallmetzer, Haselmair, Mautner, Tomašových and Zuschin2018; Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019); and for Venice, Lucinella divaricata (341 shells) was used with additional dates from Aequipecten (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019).
The history of relative sea-level changes since the last glacial maximum in the northern Adriatic is characterized by a rapid deepening during the transgressive phase, followed by a slow increase in accommodation space over the past 5000 years (Antonioli et al. Reference Antonioli, Anzidei, Lambeck, Auriemma, Gaddi, Furlani, Orrù, Solinas, Gaspari, Karinja and Kovačić2007; Vacchi et al. Reference Vacchi, Marriner, Morhange, Spada, Fontana and Rovere2016). Prograding highstand sediment wedges formed near river deltas (Po, Tagliamento, and Isonzo Rivers; e.g., Zecchin et al. Reference Zecchin, Gordini and Ramella2015; Amorosi et al. Reference Amorosi, Barbieri, Bruno, Campo, Drexler, Hong, Rossi, Sammartino, Scarponi, Vaiani and Bohacs2019) and in bays (e.g., Koper Bay; Novak et al. Reference Novak, Šmuc, Poglajen, Celarc and Vrabec2020), whereas transgressive lags formed in current-swept and sediment-starved environments in the Gulf of Venice, in the Gulf of Trieste, and off Istria (Cattaneo et al. Reference Cattaneo, Trincardi, Asioli and Correggiari2007; Trobec et al. Reference Trobec, Šmuc, Poglajen and Vrabec2017; Novak et al. Reference Novak, Šmuc, Poglajen and Vrabec2019). In spite of this along-strike variability in sediment supply (Trincardi et al. Reference Trincardi, Correggiari and Roveri1994; Catuneanu Reference Catuneanu2019), the timing of the maximum ingression and the shift in the stacking from retrogradation to progradation are relatively isochronous in the whole region at millennial timescales. Although beyond the seismic resolution, Holocene sediments in our cores can be subdivided at the centimeter or decimeter scale of resolution into distinct lithologic units that can be correlated (on the basis of geochronological data) with the late-transgressive phase, the phase of maximum flooding, and the highstand phase (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019). The starved locations were in the aggradational regime even during the highstand phase due to a slow increase in accommodation space (and locally even under an erosional regime), but the geochronological dating indicates that sediments in the uppermost decimeters of cores collected at Venice, Brijuni, and Piran are dominated by bioclastic (molluscan) particles produced over the past two millennia rather than by skeletal remains reworked from the transgressive sediments (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019; Tomašových et al. Reference Tomašových, Gallmetzer, Haselmair, Kaufman, Mavrič and Zuschin2019, Reference Tomašových, Gallmetzer and Zuschin2022). These deposits thus represent thin condensed sediment veneers that are isochronous with thick prograding, highstand wedges in parts of the Northern Adriatic Sea characterized by high terrigenous sediment supply. Distinct shell beds at Piran and Brijuni formed by epifaunal bivalves are located in units that were deposited during the maximum ingression (~5000–7000 years ago), in accordance with the predicted occurrence of such beds at the maximum flooding surfaces (separating the TST from the HST) due to the lack of dilution and higher productivity of epifauna under slow sedimentation rates (e.g., Kidwell Reference Kidwell1988, Reference Kidwell1993; Kondo et al. Reference Kondo, Abbott, Kitamura, Kamp, Naish, Kamataki and Saul1998; Fürsich et al. Reference Fürsich, Alberti and Pandey2021). High abundance of large epifaunal mollusks, high degree of skeletal alteration, and taphonomic feedback triggering high abundance of encrusters confirmed that the maximum flooding zone (MFZ) shell bed at Brijuni formed under limited clastic supply that was compensated by high carbonate production of mollusks, leading to temporally invariant long-term sediment accumulation (Tomašových et al. Reference Tomašových, Gallmetzer and Zuschin2022). A similar, although less distinct shell bed with scallops and concretions marks the MFZ at the Venice station (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019: fig. 7).
To compare long-term ecological trajectories among stations, core increments were assigned on the basis of grain-size, taphonomic, and paleoecological data to four intervals following the sequence-stratigraphic terminology of model III of Catuneanu et al. (Reference Catuneanu, Abreu, Bhattacharya, Blum, Dalrymple, Eriksson, Fielding, Galloway, Gibling, Giles, Holbrook, Jordan, Kendall, Macurda, Martinsen, Miall, Neal, Nummedal, Pomar, Posamentier, Pratt, Sarg, Shanley, Steel, Strasser, Tucker and Winker2009), in which the highstand phase terminates before the base-level fall (Hunt and Tucker Reference Hunt and Tucker1992): (1) a transgressive systems tract (TST) with transgressive lags and mixed assemblages formed by marine, brackish, and freshwater species; (2) an MFZ with densely packed molluscan shell beds and increments rich in concretions; (3) an early highstand systems tract (eHST); and (4) a late highstand systems tract (lHST). The shallower stations at Venice and Piran (21–23 m) were flooded ~9000–10,000 yr BP, and the deeper station at Brijuni (44 m) ~11,300 yr BP (Antonioli et al. Reference Antonioli, Anzidei, Lambeck, Auriemma, Gaddi, Furlani, Orrù, Solinas, Gaspari, Karinja and Kovačić2007; Storms et al. Reference Storms, Weltje, Terra, Cattaneo and Trincardi2008; Vacchi et al. Reference Vacchi, Marriner, Morhange, Spada, Fontana and Rovere2016). Thus, all cores generally preserve similar durations of the Holocene record. However, they differ in the thickness of the TST and HST units: the MFZ is located at 90–105 cm at Venice and at 90–120 cm at Brijuni, but only at 16–35 cm at both Piran stations (Fig. 1). We note that the bases and tops of the MFZ shell beds and concretion horizons that separate the TST and the HST at all stations are also geochronologically well defined, effectively represent distinct stratigraphic surfaces. However, all intervals were subjected to intense bioturbation so that transitions between systems tracts are ultimately represented by gradual transitions at the scale of several centimeters or even decimeters. Uppermost increments deposited during the late highstand (effectively anthropogenic) phase (lHST) are separated from increments deposited during the early highstand (lHST) by geochronological data and by the strong increase in organic and inorganic pollutants at most stations (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019: fig. 10).
Analyses
DFs between higher taxa, between stations, and between systems tracts were compared with the Wilcoxon rank-sum test and with the Kruskal-Wallis test (followed by the pairwise Wilcoxon test), using a posteriori Bonferroni correction for multiple comparison. To account for the unequal number of increments per sequence-stratigraphic unit at stations, two increments per station (i.e., the minimum number available) were randomly drawn from each systems tract, resulting in eight increments per unit in total. This subsampling procedure was repeated 10,000 times to obtain mean values and corresponding 95% confidence intervals (CIs) for regional-scale median DF per increment. Only increments with at least 20 individuals for each studied taxon were used in statistical analyses.
To explore the effect of scale on the identification of potential environmental correlates of DFs, the proportion of clay (PC), an indicator of hydrodynamic conditions, total organic carbon (TOC) and total nitrogen content (TN), both proxies for sediment organic enrichment, and CaCO3 content of the sediment were used in the analyses. The PC in bottom sediments typically explains a large portion of variation in the composition of benthic communities (Sanders Reference Sanders1958; Rosenberg Reference Rosenberg1995), and in addition to water energy can be also affected by proximity of mud-sourcing river deltas and by bioturbation that determines mud erodibility (Rhoads and Boyer Reference Rhoads and Boyer1982; Brückner et al. Reference Brückner, Schwarz, Coco, Baar, Boechat Albernaz and Kleinhans2021). At the regional scale, we estimated the effect of these environmental factors and of the relative abundance of the most likely predators (i.e., muricid and naticid gastropods relative to the total number of molluscan individuals) on the DF of the total assemblage and of the molluscan classes using generalized linear mixed models (GLMM) with station included as a random effect; for the stations we used generalized linear models (GLM) (O'Hara Reference O'Hara2009; Bates et al. Reference Bates, Mächler, Bolker and Walker2015). We also used generalized least squares (GLS; Pinheiro and Bates Reference Pinheiro and Bates2006) to assess whether per-station analyses are affected by temporal autocorrelation, with first-order autoregressive correlation structure. TOC and TN were excluded from model fitting because they are strongly correlated with each other and with PC (Supplementary Figure 1). The GLMM was fit using the lme4 R-package (v. 1.-26; Bates et al. Reference Bates, Maechler, Bolker, Walker, Christensen, Singman, Dai, Scheipl, Grothendieck, Green, Fox, Bauer and Krivitsky2020). Statistical analyses were performed with the software package R (v. 4.0.3.). All data used in this study are shown in Supplementary Table 1.
Results
DFs differ between localities and show distinct trends along individual sediment cores (Fig. 2). At the regional scale, the median per-increment DF in the total assemblage is 17%. DF is equally high in bivalves and gastropods (18%) and low in scaphopods (7%) (Fig. 3, Table 2).
Stratigraphic changes in drilling frequencies (DFs) in the total molluscan assemblage and molluscan classes at the four sampling stations with 95% binomial confidence intervals (CIs); vertical lines represent the regional median DF of the total assemblage. Note that Piston cores differ in length between 136 and 160 cm. eHST, early highstand systems tract; IHST, late highstand systems tract; MFZ, maximum flooding zone; TST, transgressive systems tract.
Box plots of drilling frequencies (DFs) of the total assemblage and the three classes at the regional scale. Numbers above box plots refer to number of samples with at least 20 individuals.
Results of statistical tests for significance of differences in drilling frequencies (DFs). Bold values indicate significance level < 0.01; underlined values indicate significance level < 0.05.
A spatial comparison reveals that the median per-station increment-level DF of the total assemblage is similarly high at Brijuni (17%), Piran 1 (19%), and Piran 2 (22%), but is very low at Venice (7%). DFs for bivalves and gastropods range from 15% to 25% at Piran 1, Piran 2, and Brijuni, although in changing rank order (Fig. 4, Table 2). The very low median DF of the total assemblage at Venice is shared by both bivalves (6%) and gastropods (8%) (Fig. 5, Table 2).
Box plots of drilling frequencies (DFs) of the total assemblage and the three classes at stations. Numbers above box plots refer to number of samples with at least 20 individuals.
Box plots of drilling frequencies (DFs) in systems tracts at stations (i.e., local scale) for the total assemblage and molluscan classes. Numbers above box plots refer to number of samples with at least 20 individuals. eHST, early highstand systems tract; IHST, late highstand systems tract; MFZ, maximum flooding zone; TST, transgressive systems tract.
At the local scale, a distinct trend emerges at most stations: DFs of the total assemblage are low in the TST and MFZ and high in the eHST and lHST. DFs increase from ~12% in the TST to ~25% in the HST at Brijuni, and from 16%–21% in the TST to 28%–30% in the HST at the Piran stations. At Venice, however, DFs are low across all systems tracts (Fig. 5). The stratigraphic increase in DF holds for bivalves and gastropods (Fig. 5) and for the families Rissoidae and Cerithiidae (Fig. 6). However, it is absent in most families, genera, or species (Figs. 6, 7).
Box plots of drilling frequencies (DFs) for the eight most abundant families in systems tracts at stations. Numbers above box plots refer to number of samples with at least 20 individuals. eHST, early highstand systems tract; IHST, late highstand systems tract; MFZ, maximum flooding zone; TST, transgressive systems tract.
Box plots of drilling frequencies (DFs) for the eight most abundant species in systems tracts at stations. Numbers above box plots refer to number of samples with at least 20 individuals. eHST, early highstand systems tract; IHST, late highstand systems tract; MFZ, maximum flooding zone; TST, transgressive systems tract.
At the regional scale, when data from all stations are considered together, the median DFs hardly differ between systems tracts at the level of the total assemblages (16%–20%) (Fig. 8, Table 2) and of bivalves (17%–19%), gastropods (15%–22%), and scaphopods (5%–12%) (Fig. 8). Standardizing the number of increments per station and systems tract demonstrates, however, that DFs are significantly higher in the HST than in the TST and MFZ, both for the total assemblage and for individual classes (Fig. 8).
Box plots of drilling frequencies (DFs) in systems tracts at the regional scale for the total assemblage and molluscan classes. Results are shown for raw data (box plots) and after resampling to account for different number of samples per systems tracts at individual stations (mean value and 95% confidence intervals [Cis] across 10,000 iterations). Numbers above box plots refer to number of samples with at least 20 individuals. eHST, early highstand systems tract; IHST, late highstand systems tract; MFZ, maximum flooding zone; TST, transgressive systems tract.
At the regional scale, the abundances of predators, PC, and CaCO3 significantly and positively correlate with DFs in the total assemblage (Fig. 9). The results are similar for bivalves and gastropods, although PC has no significant influence on DFs in bivalves and abundance of muricids has no significant influence on DFs in gastropods. In contrast, at the local scale of individual stations, the relationships between DFs and these environmental variables are divergent, mainly owing to between-station differences in the raw relationship between PC and CaCO3. The Spearman rank relation between PC and CaCO3 is negative at both Piran stations (r = −0.38, p = 0.28 at Piran 1, r = −0.9, p = 0.0009 at Piran 2) and at Venice (r = −0.6, p = 0.1), but it is strongly positive at Brijuni (r = 0.77, p = 0.01). Differencing these variables leads to low and insignificant correlations between PC and CaCO3 at the scale of stations. However, when accounting for temporal autocorrelation in GLS analyses with single predictors, the contrast between the local and regional scales persists: the effects of PC on DFs are negative at both Piran stations (slope = −0.07, p = 0.0001 at Piran 1, slope = −0.03, p = 0.09 at Piran 2), positive at Brijuni (slope = 0.04, p = 0.08), and negligible at Venice (slope = −0.003, p = 0.7). The effects of CaCO3 on DFs are consistently positive at all stations. Within-station GLM analyses with four predictors show that DFs in the total assemblage at Piran 1 increase with decreasing PC and with increasing abundances of muricids, while at Piran 2 they increase with PC, CaCO3, and abundance of naticids. At Brijuni, DFs of the total assemblage increase with PC, and no significant relations are evident for DFs at Venice. These results are largely similar for bivalves and gastropods (Fig. 9, Table 3). Accounting for temporal autocorrelation in GLS with four predictors, DFs in the total assemblage increase with decreasing PC at Piran 1 and with increasing CaCO3 at Piran 2 and Brijuni. However, as mentioned earlier, the within-station collinearity between PC and CaCO3 affects these analyses and hides the observation that DFs decline with increasing clay content at the Piran stations, whereas they increase with increasing clay content at Brijuni.
Coefficient plot showing results of the fitted generalized linear mixed models (GLMM) and generalized linear models (GLM) for the drilling frequencies (DFs) in the total assemblage and classes in the region and at each station. Effect sizes for predictors are shown with 95% confidence intervals (CIs). For significances and R 2 of the regional GLMM, see Table 3.
The fitted generalized linear mixed models (GLMM) and generalized linear models (GLM) for the total assemblage and molluscan classes at regional and local scale. Bold values indicate significance level < 0.01, underlined values indicate significance level < 0.05. R2m, marginal R squared; R2c, conditional R squared; AIC = Akaike information criterion.
Discussion
Differences in DF between Regional and Local Scales
Our results demonstrate that spatial heterogeneity in DFs is very high at local scales and between taxa, making extrapolations to entire communities and larger spatial scales problematic. At the local scale, assemblage-level DFs tend to increase during the postglacial sea-level rise (from the TST to HST) at three out of four stations, but regional-scale comparisons suggest that the median DF remained stable within the basin throughout the Holocene. This pattern can be explained by strong habitat variability within systems tracts typical for higher-order sequences like the one studied here (e.g., with the development of seagrass communities at Venice station and level-bottom communities at other stations at similar depths; Scarponi and Kowalewski Reference Scarponi and Kowalewski2007; Zuschin et al. Reference Zuschin, Harzhauser, Mandic, McGowan and Smith2011) and by the differences in thickness of systems tracts between the individual stations.
Spatial distribution and abundance of both predator and prey taxa vary between habitats (e.g., Sawyer and Zuschin Reference Sawyer and Zuschin2010; Leonard-Pingel and Jackson Reference Leonard-Pingel and Jackson2016), which in turn shift laterally through time in response to changes in sea level and available accommodation space (Holland Reference Holland2000), producing systematic changes in facies—and consequently DFs—observed at stratigraphic sections. As different stations within the basin differ in their environmental histories and thicknesses of systems tracts, the local trends in DF can be averaged out regionally and thus do not necessarily scale up to the entire basin.
Such distinct stratigraphic patterns in facies and unit thickness are evident among the studied cores. For example, owing to earlier marine flooding, the HST at Brijuni (fining-upward trend) and Venice (sandy throughout) is expanded, in contrast to the condensed, coarse-grained HST at Piran 1 and 2 (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019). Consequently, the HST is represented by many increments with low DFs at Brijuni and especially Venice and few increments with consistently high DFs at Piran. Pooling data from these stations together without accounting for the differences in numbers of increments dampens the regional DF. Conversely, the large number of increments combined with relatively high DFs in the stratigraphically expanded TST at the two Piran stations results in higher regional DF. Therefore, even though each station is represented by a similar number of equally spaced samples (Figs. 1, 2), the regional patterns are still biased by the variable thickness of systems tracts and thus by the number of samples representing each stratigraphic unit. As a result, local and regional stratigraphic patterns in DFs appear decoupled. However, a consistent temporal trend of increasing DFs emerges at both scales when the same number of samples from each station and systems tracts is used (Fig. 8). Therefore, integration of information on sequence-stratigraphic architecture and facies changes at individual locations/sections across a basin provides the most useful framework for tracking DF patterns at regional scales.
The Advantages of Analyzing DF on Assemblage and Higher-Taxon Levels
Our results suggest that the regional within-sequence trends in DFs and along environmental gradients manifest primarily at higher taxonomic levels or at the level of whole assemblages, not at the level of species or genera. Lower taxa are typically confined to specific habitats, which vary in their geographic extent through time. Therefore, most species, genera, and even families do not range through the studied Holocene sequence at sufficient abundance to allow meaningful comparison of DF patterns across the basin or between systems tracts.
Abundances of predators and environmental factors do not correlate with local DFs in a consistent way. At the regional scale, however, DF largely increases with increasing PC (and correlated increase in organic enrichment and/or decrease in hydrodynamic energy driven by the increase in water depth), with increasing predation pressure (as reflected in the relative abundance of naticid and muricid gastropods), and with increasing CaCO3 content. These relationships suggest that over the timescales of thousands of years, in which the regional species pool remained largely stable, an increase in organic enrichment and an increase in water depth enhances predator–prey interactions in the metacommunity (e.g., Menge and Sutherland Reference Menge and Sutherland1987; Leibold Reference Leibold1989).
No Evidence for Unusually Low Predation Pressure in the Northern Adriatic Sea
The northern Adriatic Sea was considered to be an area of overall unusually low predation pressure, a background condition that would allow sessile epibenthos to flourish in the oligotrophic waters of the eastern part of this basin (McKinney and Hageman Reference McKinney and Hageman2006; McKinney et al. Reference McKinney, Hageman and Jaklin2007). By pooling all increments at the regional scale, we show, however, that DFs during the Holocene have remained at levels very typical for warm-temperate seas in the Cenozoic (Yochelson et al. Reference Yochelson, Dockery and Wolf1983; Kelley and Hansen Reference Kelley, Hansen, Kelley, Kowalewski and Hansen2003). Similar to other regions, DFs are high for bivalves and gastropods and significantly lower for scaphopods. Our results demonstrate that DF in the northern Adriatic varies between and within systems tracts, generally increasing toward deeper habitats. Previous studies targeting surface death assemblages documented high environmental variability in DF in that basin as well. For example, Huntley and Scarponi (Reference Huntley and Scarponi2015) found that DF is highest in sediment-starved environments north of the Po delta, while it is considerably lower in high-sedimentation environments south of it. Sawyer and Zuschin (Reference Sawyer and Zuschin2010) demonstrated that DFs in subtidal environments of the Isonzo delta and Bay of Panzano were of typical modern levels, while those from tidal flats were very low. However, subtidal habitats in the northeast Adriatic Sea, now inhabited by abundant epifaunal suspension feeders at Piran and Brijuni, have consistently high DFs (>25%), contradicting the hypothesis that low sea-surface primary productivity typical for these areas is associated with reduced predation intensities. At the same time, our study also suggests the very low DFs can be a persistent feature in subtidal areas with seagrass development (see following section). However, this pattern only further underscores high spatial variability in drilling predation pressure frequently observed at regional scales (e.g., Hansen and Kelley Reference Hansen and Kelley1995; Hoffmeister and Kowalewski Reference Hoffmeister and Kowalewski2001; Sawyer and Zuschin Reference Sawyer and Zuschin2011) and does not qualify low predation intensity as a key condition in the northern Adriatic Sea.
Reduced Drilling Frequencies in Seagrass Habitats
At Venice, DFs are significantly lower than at the other three stations. This unusually low DF is possibly due to pervasive seagrass development at that site throughout the Holocene succession, as indicated by the high abundance of the chemosymbiotic lucinid bivalve Lucinella divaricata, which occurs preferentially in sediments of seagrass meadows, and by the occurrence of strictly seagrass-associated rissoid gastropods such as Rissoa monodonta and Rissoa membranacea (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019). The millennial-scale stability of seagrass ecosystems is well proven (e.g., Mateo et al. Reference M. A., Renom and Michener2010; Leiva-Dueñas et al. Reference Leiva-Dueñas, López-Merino, Serrano, Martínez Cortizas and Mateo2018; Hyman et al. Reference Hyman, Frazer, Jacoby, Frost and Kowalewski2019), but the significance of seagrass development for predator–prey interactions is controversial. Experiments in modern environments suggest that the root mat typically formed by seagrasses protects infauna from burrowing predators (Orth et al. Reference Orth, Heck and van Montfrans1984), and this is particularly true for infaunal bivalves (Peterson Reference Peterson1982; Irlandi Reference Irlandi1997; Goshima and Peterson Reference Goshima and Peterson2012). On the other hand, DFs in the Neogene fossil record of the Caribbean were shown to be uniformly higher in biogenic habitats such as seagrass meadows compared with soft sediments (Leonard-Pingel and Jackson Reference Leonard-Pingel and Jackson2016). The very low DFs at Venice are in line with the experimental studies and suggest that seagrasses provided mollusks a refuge from drilling predation over the paleoecologically relevant timescale of millennia.
Acknowledgments
Thanks to the captain of the coring vessel, J. Sedmak, to numerous students for their help with processing of samples, and to D. Chattopadhyay and P. Vasconcelos for sharing information on the ecology of muricid predators. This work was funded by the Austrian Science Fund (FWF): P24901. A.T. was supported by the Slovak Research and Development Agency (APVV 0555-17) and by the Slovak Scientific Grant Agency (VEGA 2/0169-19). Two reviewers provided helpful comments on the article.
Data Availability Statement
Data available from the Dryad Digital Repository: https://doi.org/10.5061/dryad.qfttdz0jj.
Supplementary Figure 1. Correlations between environmental variables (TN, TOC, CaCO3, PC).
Supplementary Table 1. All data used for this study.
Open access
