Long-term baseline data that allow tracking how predator-prey interactions have responded to intensifying human impacts are often lacking. Here, we assess temporal changes in benthic community composition and interactions between drilling predatory gastropods and their molluscan prey using the Holocene fossil record of the shallow northern Adriatic Sea, which is characterized by a long history of human transformation. Molluscan assemblages differ between the Isonzo and Po prodelta, but both show consistent temporal trends in the abundance of dominant species. Samples of mollusc prey collected at high stratigraphic resolution indicate that drilling frequencies have drastically declined in the Po prodelta since the mid-20th century, while a weaker trend in the more condensed sediments of the Isonzo prodelta is not statistically significant. The decrease in drilling predation intensity and the community turnover is related to the loss of predatory gastropods, increased relative abundance of less-preferred prey, and a rapidly progressing infaunalization during the most recent decades. Our results align with data showing the substantial depletion of marine resources at higher trophic levels in the region and indicate that the strong simplification of the food web initiated in the late 19^th century accelerated further since the mid-20th century.

Sediment cores were collected at three stations in the NAS in 2013, with two stations located in the Po prodelta at 21 m water depth (Po 3, Po 4), and one station in the Isonzo prodelta  in the Gulf of Trieste at 12 m water depth (Panzano). These two deltaic systems are characterized by the highest sediment accumulation rates in the northern Adriatic Sea and thus provide high, decadal-scale resolution of stratigraphic archives. At all three stations the seafloor is muddy. Changes in community composition and DF were assessed in two 16-cm-diameter replicate cores at each station, which were taken a few meters from each other. These replicate cores are very similar in faunal composition and were pooled in this study to increase sample size per increment when analyzing DFs and community composition. The cores were sliced into 2-cm increments for the uppermost 20 cm, and 5-cm-increments below. The smaller increments at the top may allow a higher resolution for surface strata, which are potentially more intensively affected by anthropogenic impacts. However, the thickness of the surface mixed layer (6 cm at Panzano and 16-20 cm at Po) and the decadal to centennial scale of time-averaging show that pairs of adjacent 2-cm increments can be pooled to obtain a better comparability with the 5-cm increments used for the rest of the core.
The core age models are based on radiocarbon-calibrated amino acid racemization of shells of V. gibba from cores M13 and M21 in the Po delta and core M28 in the Isonzo prodelta. Time averaging of molluscan assemblages in individual increments varies between 10-50 years at the Po prodelta and between ~20 years and a few centuries at the Isonzo prodelta.  Four distinct temporal units were distinguished using the age models based on V.gibba, including the middle to late Holocene with 40 samples, 17th to early 19th century with 21 samples, late 19th to early 20th century with 32 samples, and late 20th to 21st century, with 42 samples. 
Sedimentation rates at the station Panzano are ~2-4 mm/year and the sedimentary record of the cores covers approximately the last 500 years of ecosystem history, providing baseline data on variability in DF and community composition before the onset of strong human impacts in the 20th century. Conversely, the sediment cores from the Po prodelta provide a higher stratigraphic resolution but extend back only to the late 19th century due to high sedimentation rates of ~10-20 mm/year. Therefore, we included older samples from three sediment cores (205-S9, EM-S7 and EM-S13) drilled on the coastal Po plain near Ferrara to evaluate the significance of the recent ecological changes in the context of longer temporal trends. These older samples represent middle to late Holocene prodelta sediments (1.1 to 7.0 ky BP), deposited under similar environmental conditions but subsequently buried under the prograding delta front and coastal plain deposits, and thus currently located at 19 to 26 m below the surface of the present-day Po plain. 
All samples were sieved using a 1 mm mesh size and only complete molluscan shells and valves (>90% preserved) were counted and identified to the species level. Abundance of bivalves was calculated as the number of articulated specimens plus the combined number of left and right valves divided by two. To estimate DF (i.e., the percent of prey specimens with complete drill holes representing a proxy for the frequency of successful predatory attacks), 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 following equation: DF = DV/[(RV + LV)/2 + A], where DV is the total number of drilled valves, RV is the total number of right valves, LV the total number of left valves, and A the number of articulated individuals. Only samples with at least 20 individuals were used in the analyses. In addition to the DF of the total assemblage we also calculated DFs for bivalves, gastropods and the four most abundant species, including the bivalves Varicorbula gibba and Kurtiella bidentata and the gastropods Turritellinella tricarinata and Tritia varicosa. The four species, however, are not abundant enough in temporal bins to make a rigorous comparison of DFs within and between cores at the species level (electronic supplementary material, table S3). Whereas high-resolution analyses of DFs is not possible for individual taxa, we supplemented assemblage level analyses by coarse-resolution assessment of the four common species. In that case, core samples were pooled into pre-anthropogenic (middle-late Holocene, 17th to early 19th century, late 19th to early 20th century) and post-anthropogenic (late 20th to early 21st century) and Fisher’s test was applied to assess for statistical differences between the two sample groups. Predatory gastropods, which make up only 1.38% of the total molluscan assemblage, are strongly dominated by two species of naticids, Euspira nitida (70%) and E. macilenta (20%), while the few muricids (9%) are mostly juveniles of Hexaplex trunculus and Bolinus brandaris. We were not able to distinguish consistently between drill holes made by these two families of shell drilling gastropods.
We evaluated stratigraphic changes in the composition of molluscan assemblages using non-metric multidimensional scaling ordination (NMDS; k=2 dimensions) based on square-root-transformed proportional species abundances and Bray-Curtis dissimilarities. The differences in molluscan community composition between the four temporal units were tested with non-parametric permutational multivariate analysis of variance (PERMANOVA). Differences in DFs between the units were assessed with the Kruskal-Wallis test and the pairwise Wilcoxon rank-sum test, using a posteriori Bonferroni correction for multiple comparisons. Distance-based redundancy analysis was used to assess whether differences in species-level community composition are related to differences in DF.  The principal coordinate analysis using Bray-Curtis dissimilarity was performed on the square-root-transformed proportional species abundances, and the resulting ordination scores were subjected to redundancy analysis, using the function “capscale” in the vegan R package.
To assess stratigraphic trends in the relative abundance of drilling predators, a rank correlation between the relative abundance of predators and depth in the core was compared to a resampling distribution of Spearman rho values under a null model in which the relative abundance of predators is constant along the core. The resampling distribution is based on a Monte Carlo simulation of vertical trends generated by random drawing of drilling predators with probability given by their overall relative frequency in cores. At each core level, the total number of simulated specimens was based on the actual sample size for this level. For each simulation run, a Spearman rho value was calculated for random data. Positive rho values indicate an upward decline and negative values indicate an upward increase in predator abundance toward the top of the core. The two-tailed significance value p is based on the number of replicate runs in which the simulated absolute values of rho exceeded the observed absolute rho value. The results are based on 10000 iterations.
To assess trends in DF, a serial Monte Carlo model was used based on the DF data in a given core. In each iteration of the model, a random walk was simulated with an initial drilling frequency value based on the value observed at the base of the core (when the starting value was selected randomly from among all values observed in the core, the results (not shown) are effectively the same). Then, values were shifted by adding one first difference, using a first difference value observed between two randomly selected adjacent levels. The process was repeated for all successive levels resulting in a random time series based on the actual first differences observed in the core. If the shift resulted in an impossible negative value of DF, the shift was reversed. Alternative runs with negative values converted to 0 (not shown) yielded comparable outcomes. The process was repeated 10000 times and the resulting set of 10000 time series was used to estimate the range of possible trends expected under a random walk. The results were also used to estimate significance of changes in DF across the core. In this approach, the core was divided into two samples (upper and lower parts of the cores) starting from the third increment in the core and then successively shifting the boundary between the lower and upper parts to find the dividing point that results in the maximum difference in DF. The difference was measured using t-statistic (difference between the mean DFs between the two parts of the core divided by the standard error). The resulting trend was then compared against 95% confidence bands based on the analysis of the 10000 Monte Carlo series simulating non-directional random trends. In addition, a Spearman rank correlation was computed for drilling frequency and core depth (positive correlation indicates an upward decrease in drilling predation and negative correlation indicates an upward increase in DF). Distribution of expected values of Spearman rho was computed based on 10000 Monte Carlo runs. Monte Carlo simulations and spearman rank correlations were performed for the total assemblage, gastropods and bivalves at all three stations (except for gastropods at station Po3 because of a low number of samples with minimum 20 individuals). Statistical analyses were performed with the software package R (v. 4.3.1).