Introduction

Predation is one of the key factors shaping marine ecosystems today (Paine 1974; Connolly and Roughgarden 1999; Freestone et al. 2011; Harper and Peck 2016) and through the Phanerozoic (Vermeij 1977; Harper 2003) and is hypothesized to be responsible for secular changes in the guild composition and energetics of benthic ecosystems (Bottjer and Ausich 1986; Bambach 1993; McRoberts 2001; Aberhan et al. 2006). Drill holes provide direct evidence of predation and are thus the most often studied predator–prey interaction in the fossil record (Kowalewski 2002). 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. 1997; Gili et al. 2006; Bowden et al. 2011). However, it is unclear how species-level drilling frequencies (DFs) scale to community levels (Kelley and Hansen 2006) and how local estimates scale up to regional or global scales (Hoffmeister and Kowalewski 2001). Given the nonrandom preservation of depositional environments in the rock record (Holland 2000; Smith 2007) 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 2016).

First, global compilations demonstrate a distinct Phanerozoic increase in DF at the level of individual prey populations (e.g., Huntley and Kowalewski 2007). 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 2005). Second, global or regional patterns in local-scale DF are frequently perceived at the level of the total assemblage (Kowalewski et al. 1998), which may reflect trends that are driven by varying relative abundances of prey species and differences in their susceptibility to drilling predation (Leighton 2002). 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. 2017).

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 (2003) and McKinney and Hageman (2006) 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 (2009) 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 2010). However, as ecological surveys and DFs observed in surface death assemblages are affected by anthropogenic impacts (Smith et al. 2018), 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 1998; Leonard-Pingel and Jackson 2016; Patzkowsky 2017). 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 2007: p. 243) and were similar to average Cenozoic levels (e.g., Baumiller et al. 2006; Klompmaker 2009; Sawyer and Zuschin 2011; Chattopadhyay et al. 2016; Neely et al. 2021) (rather than to Paleozoic or early Mesozoic levels before the onset of the Mesozoic marine revolution; Kowalewski et al. 2005).

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. 1984; Leonard-Pingel and Jackson 2016).

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. 2016; Gallmetzer et al. 2019). 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. 2019).

Figure1.

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. (2019).

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 1981; Kelley and Hansen 2003), 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 1981; Kabat 1990). For the muricids encountered in our samples, shell drilling has been documented for H. trunculus (Peharda and Morton 2006; Morton et al. 2007; Sawyer et al. 2009) and B. brandaris (Abidli et al. 2012; 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 1960). 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. (2012): 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.

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. 2018; Gallmetzer et al. 2019; Tomašových et al. 2019); for Brijuni, Timoclea ovata (305 shells) was used with additional dates from Glycymeris, Aequipecten, Chlamys, Ostrea, and a plant remnant (Schnedl et al. 2018; Gallmetzer et al. 2019); and for Venice, Lucinella divaricata (341 shells) was used with additional dates from Aequipecten (Gallmetzer et al. 2019).

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. 2007; Vacchi et al. 2016). Prograding highstand sediment wedges formed near river deltas (Po, Tagliamento, and Isonzo Rivers; e.g., Zecchin et al. 2015; Amorosi et al. 2019) and in bays (e.g., Koper Bay; Novak et al. 2020), 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. 2007; Trobec et al. 2017; Novak et al. 2019). In spite of this along-strike variability in sediment supply (Trincardi et al. 1994; Catuneanu 2019), 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. 2019). 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. 2019; Tomašových et al. 2019, 2022). 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 1988, 1993; Kondo et al. 1998; Fürsich et al. 2021). 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. 2022). A similar, although less distinct shell bed with scallops and concretions marks the MFZ at the Venice station (Gallmetzer et al. 2019: fig. 7).

Table 1.

The number of samples and abundance of prey classes, dominant prey families and species, and predatory families per systems tract at the four stations.

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. (2009), in which the highstand phase terminates before the base-level fall (Hunt and Tucker 1992): (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. 2007; Storms et al. 2008; Vacchi et al. 2016). 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. 2019: 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 CaCO₃ 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 1958; Rosenberg 1995), 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 1982; Brückner et al. 2021). 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 2009; Bates et al. 2015). We also used generalized least squares (GLS; Pinheiro and Bates 2006) 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. 2020). 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).

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).

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).

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).

At the regional scale, the abundances of predators, PC, and CaCO₃ 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 CaCO₃. The Spearman rank relation between PC and CaCO₃ 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 CaCO₃ 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 CaCO₃ 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, CaCO₃, 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 CaCO₃ at Piran 2 and Brijuni. However, as mentioned earlier, the within-station collinearity between PC and CaCO₃ 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.

Figure2.

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.

Figure3.

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.

Table 2.

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.

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 2007; Zuschin et al. 2011) and by the differences in thickness of systems tracts between the individual stations.

Figure4.

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.

Spatial distribution and abundance of both predator and prey taxa vary between habitats (e.g., Sawyer and Zuschin 2010; Leonard-Pingel and Jackson 2016), which in turn shift laterally through time in response to changes in sea level and available accommodation space (Holland 2000), 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.

Figure5.

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.

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. 2019). 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.

Figure6.

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.

Figure7.

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.

Figure8.

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.

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 CaCO₃ 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 1987; Leibold 1989).

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 2006; McKinney et al. 2007). 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. 1983; Kelley and Hansen 2003). 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 (2015) 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 (2010) 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 1995; Hoffmeister and Kowalewski 2001; Sawyer and Zuschin 2011) and does not qualify low predation intensity as a key condition in the northern Adriatic Sea.

Figure9.

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² of the regional GLMM, see Table 3.

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.

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. 2019). The millennial-scale stability of seagrass ecosystems is well proven (e.g., Mateo et al. 2010; Leiva-Dueñas et al. 2018; Hyman et al. 2019), 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. 1984), and this is particularly true for infaunal bivalves (Peterson 1982; Irlandi 1997; Goshima and Peterson 2012). 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 2016). 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.