Introduction
The role of antagonistic interactions in the macroevolutionary process has long been debated by students of evolution. Many studies have focused on predator-prey relationships and their evolutionary significance over ecological and geological time scales (Vermeij Reference Vermeij1977; Signor and Brett Reference Signor and Brett1984; Madin et al. Reference Madin, Alroy, Aberhan, Fürsich, Kiessling, Kosnik and Wagner2006; Huntley and Kowalewski Reference Huntley and Kowalewski2007; Sallan et al. Reference Sallan, Kammer, Ausich and Cook2011, among many others). Predation has likely played an important role in the evolutionary process since the origin and early evolution of animals (Bengtson and Yue Reference Bengtson and Yue1992; Butterfield Reference Butterfield1997; Huntley and Kowalewski Reference Huntley and Kowalewski2007; Stanley Reference Stanley2008; Sperling et al. Reference Sperling, Frieder, Ramman, Girguis, Levin and Knoll2013), but what about parasites?
Most animal species are host to at least one endemic parasite in addition to many other relatively generalist parasite species. Parasitism is ubiquitous in the biosphere and yet, by contrast with predation, we know very little of its history in the fossil record. Most reports of parasitism in the fossil record are based on specimens of exceptional preservation (Moodie Reference Moodie1923; Cameron Reference Cameron1967; Fry and Moore Reference Fry and Moore1969; Conway Morris Reference Conway Morris1981, Reference Conway Morris1990; Savazzi Reference Savazzi1995; Feldmann Reference Feldmann1998; Feldmann and Brett Reference Feldman and Brett1998; Boucot and Poinar Reference Boucot and Poinar2010) and are rarely based on systematic surveys with large sample sizes (though see Brett Reference Brett1978; Ruiz and Lindberg Reference Ruiz and Lindberg1989; Baumiller and Gahn Reference Baumiller and Gahn2002; Gahn and Baumiller Reference Gahn and Baumiller2003; Schimmel et al. Reference Schimmel, Kowalewski and Coffey2012; Huntley and Scarponi Reference Huntley and Scarponi2012).
The digenean trematode–bivalve mollusk interaction holds great potential for the systematic elucidation of parasite-host interactions in deep time. Trematodes induce the growth of characteristic oval-shaped pits with raised rims within the pallial line of their bivalve intermediate hosts (Ruiz and Lindberg Reference Ruiz and Lindberg1989) (Fig. 1). Hence, these pits have a good fossilization potential, though are most easily recognized in fossils preserved in unlithified sediment.
Examples of analyzed trace morphologies. A, Trematode pits in Loripes lucinalis from Caorle. B, Trematode pits in Chamalea gallina from Montemarciano. C, Mudblister in C. gallina from Senigallia. D, E, U-shaped borings in Scapharca inaequivalvis from Cervia (external and internal views, respectively). Note how the blistering on the interior reflects a live-live interaction. F, G, Successful predatory drill holes on C. gallina from Senigallia and Caorle, respectively. H, Incomplete drill hole on C. gallina from Lido di Jesolo. The jagged outline of the inner hole indicates that the hole was not completed by the predator and was subsequently damaged by taphonomic processes. I, Edge drilling in C. gallina from Senigallia. All scale bars, 5 mm.
A number of studies have addressed the nature of geographic and environmental variation in predatory drilling frequency and prey effectiveness, highlighting their relevance for interpreting temporal trends (Vermeij Reference Vermeij1980; Hansen and Kelley Reference Hansen and Kelley1995; Hoffmeister and Kowalewski Reference Hoffmeister and Kowalewski2001; Kelley and Hansen Reference Kelley and Hansen2007; Sawyer and Zuschin Reference Sawyer and Zuschin2010, among others). In comparison to predation traces, almost nothing is known about the spatial or environmental variation of parasitic traces. Huntley (Reference Huntley2007), with a limited sample size, observed that trematode pits were more common in Protothaca staminea living in a rocky tidal creek than those from a muddy lagoon on San Juan Island, Washington, United States. Ruiz and Lindberg (Reference Ruiz and Lindberg1989) hypothesized that trematode traces should be restricted to nearshore environments because vertebrate predators (often mammals and birds, though fish as well) serve as the definitive host upon consuming the molluscan intermediate host. However, in their study on northern Adriatic/Po Plain modern and fossil bivalves, Huntley and Scarponi (Reference Huntley and Scarponi2012) found that trematode traces were restricted to the lower shoreface-offshore transition environments (~7–10m water depth) and not in the shallowest of nearshore samples. In the studied succession, trematode traces were significantly more likely to occur in bivalves preserved in deposits of the transgressive systems tracts (TST), characterized by barrier island/lagoon complexes and relatively starved nearshore settings, than in deposits of the highstand systems tracts (HST), represented by prograding delta-strandplain (see also Amorosi et al. Reference Amorosi, Colalongo, Fiorini, Fusco, Pasini, Vaiani and Sarti2004). This led to the hypothesis that trematode prevalence was controlled by differing environmental conditions among transgressive and normal regressive nearshore deposits (TST vs. HST).
In this paper we analyze geographic trends of trematodes and other types of biotic traces (Fig. 1) preserved on molluscan shells along the northern Adriatic coast of Italy, comparing two different sectors of the Adriatic Italian coast, with the Po River as the main divide (Fig. 2). Indeed, north of the Po River, coastal settings are relatively sediment starved and characterized by barrier island–lagoon complexes where sedimentation and mixing of fresh and marine waters occurs in estuaries and lagoons (i.e., TST-like; Fig. 2). Such relatively sediment starved settings were present across most of the study area during phases of Holocene transgression (see Amorosi et al. Reference Amorosi, Colalongo, Fiorini, Fusco, Pasini, Vaiani and Sarti2004). South of the Po River, on the contrary, coastal environments are characterized by strandplain complexes, especially in proximity to the current Po River outlet. The latter type of setting is characterized by stressed, high-sedimentation environments emplaced in the area after the peak of Holocene marine ingression (see also Trincardi et al. Reference Trincardi, Corregiari and Roveri1994) (Fig. 2). Hence the northern Adriatic, with its contrasting environmental characterization, is an ideal venue in which to compare the large-scale geographic variation of traces of biotic interactions. These traces include parasitic trematode pits, klepto-parasitic spionid polychaete mudblisters, and predatory gastropod drill holes. We also analyze spionid U-shaped borings that give valuable information on the geographic distribution of spionids; however, they can be produced as both ante- and postmortem traces on the host’s shell, so we interpret these traces with caution. Specifically, we test the following hypotheses: (1) trematode prevalence will be higher in relatively sediment starved, oligotrophic, and stable nearshore environments (e.g., north of the Po river; TST-like) and lower in stressed, eutrophic environments, owing to the Po fluvial discharges (especially in the proximity of the delta); (2) the distribution of biotic traces will predictably correlate with diversity structure and/or life mode of the benthic malacofauna (e.g., drilling predation intensity and trematode prevalence should be positively correlated with diversity [Paine, Reference Paine1966; Hechinger et al. Reference Hechinger, Lafferty, Huspeni, Brooks and Kurtis2007], or traces of suspension feeding spionids should be positively correlated with suspension-feeding hosts); and (3) there will be no significant difference between modern geographic variation of trematode prevalence trends (i.e., hypothesized higher trematode infestation north of the Po river and lower south of it) and the temporal variation of fossil trematode prevalence reported by Huntley and Scarponi (Reference Huntley and Scarponi2012) (high infestation in TST vs. low infestation in HST).
Location map of the 11 bulk samples collected (numbered circles), showing also the relative thickness of the Highstand Systems Tract (HST), coastal morphology, and fluvial network along the western flank of the Northern Adriatic Sea. The thickness of the prograding wedged deposits (HST) was determined by seismic profiling and is expressed in milliseconds (i.e., 8 ms~12 m). Data from Geologic Map of Italian Seas sheets: Venezia, Ravenna and Ancona (accessed May 4 2014 at http://www.isprambiente.gov.it/Media/carg/index_marine.html).
Materials and Methods
Bulk samples of skeletal material (1–4 liters depending upon density of material) were collected from randomly selected locations on the strandline (in 1-m2 increments as necessary to collect an adequate volume of material) at 11 locations along the northeast Adriatic coast of Italy in 2010 (Fig. 2). Skeletal material and sediment were scraped from the surface to no deeper than 2 cm. Each sample was randomly subsampled to a volume of 0.375 liters to approximate the volume of fossil sample sizes (Scarponi and Kowalewski Reference Scarponi and Kowalewski2004; Huntley and Scarponi Reference Huntley and Scarponi2012). Sediment was rinsed onsite with a 1-mm sieve. The samples were collected in February and October of 2010 (Supplementary Table 1; note the samples from Senigallia and Montemarciano were collected prior to the establishment of this project, but the samples were obtained in a similar manner, only less material was obtained). The locations are visited heavily by tourists in the summer; hence a possible biasing effect of shell removal by vacationers cannot be excluded (see Kowalewski et al. Reference Kowalewski, Domènech and Martinell2014). This bias should, however, affect the overall abundance of shell material on beaches, but not necessarily the relative abundances of taxa or individuals with or without traces. We therefore consider the tourist bias to be of minimal effect in this study. Among the collected material, analyses were restricted to the three mollusk classes because they were, by far, the most abundant and because they possess the most unambiguous traces of biotic interactions. All specimens were identified to the species level when possible.
All mollusks were examined for the presence of successful and failed predatory drill holes. Each bivalve was examined for the presence of trematode-induced pits on the interior of the valve within the pallial line, spionid-induced mudblisters, and spionid-induced U-shaped borings. Gastropods and scaphopods were not examined for the latter traces because they are not readily apparent on the shell. Prevalence values and drill hole frequency values were calculated only for species occurrences with sample sizes of ten or more specimens/valves.
Trematode prevalence was calculated as the number of infested valves divided by the total number of valves within a species occurrence. Trematode parasites typically induce traces in both valves of the individual they infest. In the rare instance that traces occur in only one valve (14% among Pacific Northwest venerids) there is no statistically significant preference for one valve over the other (Huntley Reference Huntley2007). Huntley and Scarponi (Reference Huntley and Scarponi2012) used the term “infestation frequency” to describe what we herein refer to as prevalence. The change of terminology is in keeping with the desire for a standardized vocabulary among parasitologists and ecologists (Margolis et al. Reference Margolis, Esch, Holmes, Kuris and Schad1982; Bush et al. Reference Bush, Lafferty, Lotz and Shostak1997; Rózsa et al. Reference Rózsa, Reiczigel and Majoros2000).
Prevalence values for spionid-induced mudblisters and spionid-induced U-shaped borings were calculated by dividing the number of infested valves by the total number of valves within a species occurrence. Drilling frequency was calculated by dividing the number of individuals with a predatory drill hole by the total number of individuals (bivalve n/2 to correct for two valves constituting one individual and one borehole required to kill them “both”) within a species occurrence. Ninety-five percent confidence intervals for these proportional data were calculated in SAS 9.3 based on binomial and F distributions (M. Kowalewski personal communication 2003).
In some cases, sample size for a given species was rather low, and it was difficult to determine if an absence of traces truly reflected an absence of the interaction or was the spurious result of a small sample size. To address this problem, a rarefaction analysis was conducted to determine the minimum number of specimens required to record an interaction. For each taxon and interaction combination, the species occurrence with the lowest non-infested to infested ratio was selected for rarefaction analysis (e.g., C. gallina from Cervia for U-shaped borings). Counts of non-infested and infested taxa were treated as separate taxa in a rarefaction analysis and the mean number of individuals (from 1000 iterations) needed to record an interaction (i.e., for the “taxa” count to equal two) was calculated. In the case of C. gallina from Cervia with U-shaped borings, a mean number of 45 specimens was needed to identify one specimen with a trace (44 non-infested and 1 infested). Therefore, the threshold value for this taxon/interaction combination was 44. The species-level prevalence or frequency data discussed here pertain only to samples that met the minimum threshold sample size for a given taxon/interaction combination.
Summed prevalence and drilling frequency values were calculated for each trace at each location to serve as an estimate for assemblage-level interaction intensity (Hechinger et al. Reference Hechinger, Lafferty, Huspeni, Brooks and Kurtis2007). These summed prevalence/frequency values were compared with one another, diversity metrics, and life mode (Supplementary Table 2) via Spearman’s rank correlation. The aforementioned threshold sample-size criterion was not applied to the summed prevalence/frequency data.
Species richness and evenness were calculated from the molluscan abundance data matrix (bivalves, gastropods, and scaphopods). Rarefaction analysis was conducted for each sample to calculate mean species richness at a standardized sample size of 90 individuals.
Trematode prevalence values were compared between modern specimens reported in this study and the fossil specimens reported in Huntley and Scarponi (Reference Huntley and Scarponi2012). The prevalence values were arcsine-transformed and subjected to a Mann-Whitney U-test. An α-value of 0.05 is assumed for all statistical analyses, which were carried out in PAST v.2.11 (Hammer et al. Reference Hammer, Harper and Ryan2001).
Results
A total of 17,299 molluscan valves were identified and represent 111 species (55 bivalves, 53 gastropods, and 3 scaphopods; Supplementary Table 3). Traces of biotic interactions were very common in the samples observed, vary spatially in their frequency (Table 1, Figs. 3–6), and are correlated with various diversity indices and life mode (Tables 2, 3). The results for each type of trace are reported below.
Trematode prevalence at the species occurrence level, with 95% confidence intervals and number of valves for taxa determined to have been adequately sampled via rarefaction analysis. See Materials and Methods for more details on rarefaction analysis and supplementary materials for data on all species occurrences.
Spionid mudblister prevalence at the species occurrence level. Details as in Figure 3.
Spionid U-shaped boring prevalence at the species occurrence level. Details as in Figure 3.
Predatory gastropod drill hole frequency at the species occurrence level. Details as in Figure 3.
Number of specimens of molluscan species with biotic traces (specimens n≥10) from the northeast coast of Italy (number of valves for bivalves and number of shells for gastropods). (G) indicates gastropod species. Trem, trematode-induced pits; Mud, spionid polychaete-induced mudblister; U, spionid polychaete-induced U-shaped boring; DH, complete predatory gastropod drill hole; Inc, incomplete predatory gastropod drill hole; Edge, complete predatory gastropod edge drill hole (on commissure of bivalve prey).

Summed prevalence/frequency and diversity indices by location.

Spearman’s rank correlation coefficients and p-values for comparisons between summed prevalence/frequency values, diversity indices, and life mode data by location. NS = not significant (α=0.05).

Trematode-Induced Pits
Trematode pits were found from eight of the 11 locations sampled, among six bivalve species: Cerastoderma glaucum, Chamelea gallina, Donax semistriatus, Loripes lucinalis, Mytilus galloprovincialis, and Spisula subtruncata (Fig. 3). The species were present in 15 out of 30 species occurrences culled by abundance (see rarefaction analysis described above and SOM).
Trematode prevalence values ranged from 0.002 (C. gallina from Lido di Jesolo) to 0.169 (M. galloprovincialis from Montemarciano) and displayed systematic large-scale geographic variation. North of the Po River (from Grado to Chioggia; Figs. 2, 3), parasitic infestation was identified in all four examined localities, among five of the six infested species (8 out of 13 species occurrences); with summed prevalence values ranging from 0.002 (Lido di Jesolo) to 0.119 (Caorle). Trematode traces are nearly absent immediately south of the Po delta (i.e., within 100 km from the main channel mouth; Figs. 2, 3). Here, only C. gallina in M. di Ravenna (prevalence=0.004) was infested out of eight occurrences of infested species at data set level (average infestation rate of 11.1%). It is noteworthy to stress that C. gallina was the most trematode-infested species among a small group of taxa reported by Huntley and Scarponi (Reference Huntley and Scarponi2012) and here, along the C. Borsetti to Cervia transect, it shows (with D. semistriatus and M. galloprovincialis) high abundance values and recurrence (Fig. 3). Summed trematode prevalence values increased again in magnitude starting from Cattolica (i.e., ~140 km southward), reaching maximum prevalence values (both at species and locality levels) in Montemarciano, which is located more than 200 km south of the Po delta main channel (Figs. 2, 3). Here parasitic infestation has been recovered from only two species (C. gallina and M. galloprovincialis Fig. 3), out of the four infested species recovered. However, along the Cattolica-Montemarciano transect the infested species C. glaucum and L. lucinalis were not retrieved (Fig. 3 and supplementary materials).
Summed trematode prevalence was negatively correlated with U-shaped borings and with mollusks of restricted mobility; it was positively correlated with the epifaunal mollusks and with molluscan herbivores on hard substrates. Furthermore, trematode traces were not significantly correlated with biodiversity measures.
Spionid-Induced Mudblisters
Mudblisters were retrieved in seven locations (18 species occurrences) in eight species: C. gallina, D. semistriatus, Lentidium mediterraneum, Ostrea edulis, Scapharca demirii, Scapharca sp., S. subtruncata, and Varicorbula gibba (Fig. 4). Their frequencies ranged from 0.001 (C. gallina from Marina di Ravenna and L. mediterraneum from Casal Borsetti) to 0.115 (S. demirii from Casal Borsetti).
In the northern part of the study area Schapharca species were not retrieved; however, a conspicuous number of specimens of infested taxa (see above) were retrieved between Grado and Chioggia (Fig. 4, supplementary materials). Of these, only C. gallina displayed spionid-induced mudblisters just north of the Po River delta (Chioggia, frequency=0.015; Fig. 4). The remaining species occurrences with mudblisters were found in all the seven locations south of the delta, with highest values and number of infested species near the Po delta between Casal Borsetti and Marina di Ravenna. In the former locality summed frequency attained the maximum registered in the dataset (i.e., 0.162), whereas in the latter locality all five recovered species were infested (Table 2, Fig. 4). No distinctive pattern emerges between the two transects of the southern study area. Frequency of infested species is similar between Casal Borsetti-Cervia and Cattolica-Montemarciano (88.1% on three localities vs. 77.5% on four localities respectively), and the abundance of infested specimens is likely driven by the presence/absence of host taxa. For instance, between S. Biagio and Montemarciano we could recover adequate sample sizes only for C. gallina and D. semistriatus (out of eight infested species) (Fig. 4, supplementary materials).
Summed mudblister prevalence was positively correlated with spionid U-shaped borings and with suspension feeders; and negatively correlated with species richness (raw and rarefied), epifaunal and herbivore mollusks, and subsurface and chemosymbiotic deposit feeders.
Spionid-Induced U-Shaped Borings
U-shaped borings by spionid worms were identified in four species (C. gallina, Scapharca inaequivalvis, Scapharca sp. and V. gibba) from five locations (eight species occurrences). As with mudblisters (with the exception of the Chioggia sample), U-shaped borings were only found south of the Po delta (Fig. 5), with highest frequency and prevalence immediately south of it (i.e., in all six species occurrences retrieved within 100 km of the Po main channel; Fig. 5). South of Cervia, U-shape borings are more patchily distributed. The notable variability in prevalence values south of the Po river seems controlled by presence/absence of Scapharca species, which show very high prevalence values wherever encountered (Table 1, Fig. 5).
Summed U-shaped boring prevalences were negatively correlated with predatory carnivorous mollusks and trematodes, and positively correlated with mudblisters.
Gastropod Drill Holes
Predation traces were recovered from six locations (16 species occurrences), among nine species: Bittium scabrum, C. gallina, Donax trunculus, Glycymeris nummaria, L. lucinalis, Lucinella divaricata, O. edulis, Turritella communis, and V. gibba. Among molluscan prey items of predatory gastropods, drilling frequencies ranged from 0.002 (V. gibba from Casal Borsetti) to 0.435 (L. divaricata from Lido di Jesolo). The majority of prey species occurrences were found well north of the Po delta from Grado-Lido di Jesolo (Figs. 2, 6) with summed drilling frequencies per locality ranging between 0.238 (Grado) and 1.056 (Lido di Jesolo) from seven species (12 species occurrences out of the 20 cumulatively retrieved north of the Po delta; 13 species occurrences including incomplete drill holes). Drilling predation was virtually not recorded in examined samples along the Po delta and within the 100 km southward (Fig. 6). Indeed, a lone instance of predation, among more than 500 specimens examined, occurred at Casal Borsetti (V. gibba, 0.002). Farther south of the delta, between Cattolica and Montemarciano, somewhat higher drilling frequencies were found. Predation was noted at Cattolica and Senigallia for three species occurrences of three different species, two of which were sampled only once north of the delta (Fig. 4). However, predation traces were quite variable, and drilling frequency seems to have been controlled by the absence of key taxa (see above). Cumulatively, between Cattolica and Montemarciano, only three out of eight preyed-upon species were found, and these were represented by five out of 32 possible species occurrences (Fig. 6). The one predated specimen of C. gallina (out of 279 specimens) from Senigallia was drilled through the commissure (edge drilling) rather than through the thicker shell wall.
A single incomplete bore hole was found in each of two specimens of C. gallina, one from Lido di Jesolo and one from Grado. Prey effectiveness, calculated as the proportion of incomplete borings among the total number of predatory attacks (complete borings+incomplete borings), was 0.091 (total number of attacks=11) for C. gallina from Lido di Jesolo and 1.000 for C. gallina from Grado (total number of attacks=1).
Summed drilling frequency was negatively correlated with sedentary mollusks and positively correlated with subsurface deposit feeders.
Trematode Prevalence Values in Modern and Fossil Bivalves
Prevalence values in modern bivalve species occurrences from this study ranged from 0.002 to 0.169 with a mean value of 0.057 (n=15) (Fig. 5). Prevalence values from comparable Pleistocene-Holocene species occurrences retrieved from core samples in the vicinity of Casal Borsetti and Marina di Ravenna (Huntley and Scarponi Reference Huntley and Scarponi2012) ranged from 0.002 to 0.450 with a mean value of 0.163 (n=7). There is no statistically significant difference in their median values (Mann Whitney U=32; p=0.16) or the shape of their distributions (Kolmogorov-Smirnov D=0.448; p=0.22).
Histogram of arcsine-transformed trematode prevalence values in modern and fossil species occurrences.
Discussion
Spatial Distribution of Traces
Analyses of biotic interactions along the Italian coast of the Northern Adriatic Sea present a striking disparity in the spatial distributions of these parasitic, klepto-parasitic, and predatory traces, which parallel the hydro-geomorphologic spatial differentiation of the Northern Adriatic nearshore environments. The main factors driving this distinction are (1) the Po fluvial discharge, which supplies sediment and nutrient-rich riverine waters to nearshore depositional environments, and (2) the geostrophic nearshore circulation, which produces an oligotrophic marine water flux that enters the Adriatic from the south along the eastern coast and flows counter clockwise through the northern Adriatic basin and drives the cold riverine waters south along the western Adriatic coast.
The northernmost sector of the Italian Adriatic coast (Grado to Chioggia) is sediment starved due to the lack of significant riverine inputs (Fig. 2). The cumulative suspended discharge of northern rivers is estimated at 0.3 * 107 t/yr (Trincardi et al. Reference Trincardi, Corregiari and Roveri1994). This, coupled with subsidence affecting much of the study area (Ferranti et al. Reference Ferranti, Antonioli, Mauz, Amorosi, Dai Pra, Mastronuzzi, Monaco, Orrù, Pappalardo, Radtke, Renda, Romano, Sanso and Verrubbi2006), has a strong effect on the geomorphologic and environmental features of this sector of the Northern Adriatic. Indeed, low fluvial discharge and local sea-level rise cause river-borne sediments to be deposited inland (e.g., alluvial-coastal plain) and fresh waters to mix with marine waters in the brackish conditions of river outlets. Hence, the morphologic outcome in the northernmost Italian Adriatic is a series of barrier island-lagoon complexes grading into a relatively sediment-starved sandy shallow shelf. Therefore, marine nearshore environments can also be characterized by relatively low turbidity, oligotrophic conditions, and a more stable marine salinity regime. These conditions are present at the sampling locations of Grado, Caorle and Lido di Jesolo, and are interpreted to be comparable to conditions along the greater part of the northern Adriatic coastline during the last two major transgressive pulses (i.e., early Late Pleistocene and early Holocene (Amorosi et al. Reference Amorosi, Colalongo, Fiorini, Fusco, Pasini, Vaiani and Sarti2004; Scarponi and Angeletti Reference Scarponi and Angeletti2008; Storms et al. Reference Storms, Weltje, Terra, Cattaneo and Trincardi2008; Wittmer et al. Reference Wittmer, Dexter, Scarponi, Amorosi and Kowalewski2014).
As our hypothesis predicted, trematode prevalence values were higher in this TST-like setting. Not only was trematode infestation more common and intense, on average, north of the Po delta than immediately south of it, but a greater diversity of host taxa were infested as well (particularly in Caorle, ~100 km north of the delta; Tables. 1, 2, Fig. 3). In contrast, along the wave-influenced Po delta and the adjacent southern strandplain (Casal Borsetti to Cervia; Fig. 2), the environmental conditions are completely different. Indeed, environmental parameters (e.g., salinity, turbidity) vary rapidly and unpredictably and sedimentation rates are extremely high (Scarponi et al. Reference Scarponi, Kaufman, Amorosi and Kowalewski2013) (Fig. 2). In short, these environments are strongly influenced by high fluvial discharge from the Po and Adriatic longshore currents (McKinney Reference McKinney2007).
Indeed, the Po delta displays an asymmetric configuration, mainly as the result of the strong shaping influence of counterclockwise surface currents, which simultaneously erode northern lobes and cause the southern drift of river-sourced fine-grained sediments (Corregiari et al. Reference Correggiari, Cattaneo and Trincardi2005; McKinney Reference McKinney2007) (Fig. 2). The total sediment load introduced to the northern Adriatic by the modern Po is 1.5 * 107 t/yr (i.e., five times that of all the northern rivers); notwithstanding the regional subsidence rate, this considerable sediment supply had, until recently, caused the delta shoreline to prograde at places more than 30 km in a few centuries (Maselli and Trincardi Reference Maselli and Trincardi2013). Freshwater discharge (and nutrients, pollution, etc.) of the Po into the Adriatic is 1500 m3 * s−1 (Nelson Reference Nelson1970) or an annual mean of ~49 km3 (Correggiari et al. Reference Correggiari, Cattaneo and Trincardi2005). The southward drift of fine-grained sediment and freshwater plumes create markedly different regimes from those found north of the delta. Here, stressed environments are widespread. As we predicted, trematode prevalence values are negligible along the delta and adjacent southern shoreline (only a few infested individuals were recovered at Marina di Ravenna and none at Casal Borsetti or Cervia; Figs. 2, 3). The moderately high trematode prevalence values for two species occurrences found in Chioggia should not come as a surprise given the proximity to the lagoon and minimal influence from the delta. Indeed, Chioggia is a special location in that it is on the northernmost flank of the delta but less than 7 km south of an inlet of the Venice Lagoon (Fig. 2).
The influence of the Po River delta (freshwater plumes) decreases with distance to the south, especially at San Biagio, Senigallia, and Montemarciano. Among these locations trematode prevalence roughly increases with distance south. We did not predict the gradient of increasing trematode prevalence with distance from delta influence, but its existence could provide further support for our hypothesis. However, an important confounding factor needs to be assessed. In this transect the geomorphologic and tectonic setting is different from the rest of the study area. Owing to tectonic uplift and higher gradient of the plain (i.e., proximity of the Apennines), nearshore deposits from the present and last interglacial are closely spaced and confined to a relatively narrow strip above sea level (coastal terraces) or within the reworking zone of fluvial/marine processes (Ferranti et al. Reference Ferranti, Antonioli, Mauz, Amorosi, Dai Pra, Mastronuzzi, Monaco, Orrù, Pappalardo, Radtke, Renda, Romano, Sanso and Verrubbi2006; Calderoni et al. Reference Calderoni, Della Seta, Fredi, Palmieri, Nesci, Savelli and Troiani2010). Hence, along the S. Biagio–Montemarciano transect, sub-fossil inputs from transgressive deposits (comparably richer in trematode traces see Huntley and Scarponi Reference Huntley and Scarponi2012) in currently forming beaches cannot be excluded and could be mixing with present dead remains. Hence, such shell beds could be time-averaged over 103 to 105 years. This could increase sample diversity and consequently increase chances of finding trematode-infested taxa. For this reason here, and for most of the discussion on biotic interactions below, we will refer to and compare only the transects north and immediately south of the Po delta, where these confounding factors are unlikely to be important. Indeed, from Cervia northward, Late Pleistocene transgressive deposits cannot contribute because they are recovered from 70–120 m core depth (Ferranti et al. Reference Ferranti, Antonioli, Mauz, Amorosi, Dai Pra, Mastronuzzi, Monaco, Orrù, Pappalardo, Radtke, Renda, Romano, Sanso and Verrubbi2006). Similarly, Holocene shallow-marine transgressive deposits are recognized only in the subsurface at several meters (>10) depth (Zecchin et al. Reference Zecchin, Baradello, Brancolini, Donda, Rizzetto and Tosi2008; Wittmer et al. Reference Wittmer, Dexter, Scarponi, Amorosi and Kowalewski2014); hence, they are well below the active mixing zone with minimal chance of introduction to modern death assemblages by river and/or sea-wave processes.
As for the rest of the study area, the environment-based hypothesis on the influence of biotic interaction was born from observations of fossil bivalves (Huntley and Scarponi Reference Huntley and Scarponi2012) and not derived from any described spatial distribution patterns of living trematodes in the northern Adriatic Sea. Most work on trematodes in the Adriatic is conducted by biologists, working with soft tissues, who are interested primarily in the infestation of piscine definitive hosts (Paradižnik and Radujković Reference Paradižnik and Radujković.2007), and occasionally molluscan intermediate hosts (Canzonier Reference Canzonier1972). Hence, the possibility of linking trematode spatial distribution to specific environmental variables is not feasible. However, Pietrock and Marcogliese (Reference Pietrock and Marcogliese2003) reported that the longevity and infectivity of marine parasites in the free-living stage are reduced when the parasites are exposed to freshwater dilutions. Also, various anthropogenic toxicants (e.g., Pb and Cd) have been shown to dramatically reduce survivorship and infectivity of free-living trematodes transitioning to and from intermediate hosts (Pietrock and Marcogliese Reference Pietrock and Marcogliese2003). The disparate pattern in trematode spatial distribution in the Adriatic could be influenced by reduced survivorship and/or infectivity of free-living trematodes in transition between egg and first intermediate host (miracidia) or between first and second intermediate hosts (cercaria) as they are exposed to freshwater or pollution. Anthropogenic pollution may explain the modern distribution of trematode traces relative to the Po delta, but, of course, it does not explain the same pattern recovered in the lower to middle Holocene and upper Pleistocene deposits.
Spionid polychaete traces in the forms of U-shaped borings and mudblisters had geographic distributions that were positively correlated to one another, not surprisingly, as these traces were likely made by congeneric, if not conspecific, individuals. Strikingly, the distribution of spionids was nearly the opposite of that displayed by trematodes. Spionids are generalists regarding feeding mode, with many taxa capable of suspension feeding and deposit feeding. Moreover, many taxa are considered to be opportunistic, thriving in habitats subjected to frequent disturbance, reduced oxygen, and\or organically enriched sediments; some consider spionid polychaetes to be potential indicators for degraded ecosystems (Dix et al. Reference Dix, Karlen, Grabe, Goetting, Holden and Markham2005 and references therein). Indeed, the Po River introduces elevated levels of nutrients, organically enriched sediments, and pollutants to the northern Adriatic which are often linked to anoxic events, especially in the summer when a decreased thermal gradient slows circulation in the basin (Stachowitsch Reference Stachowitsch1991; McKinney Reference McKinney2007 and references therein; Frontalini and Coccioni Reference Frontalini and Coccioni2008). The same unpredictable environmental conditions and stressors that excluded trematodes likely encouraged the proliferation of mudblister-inducing and boring spionid polychaetes.
Drilling frequencies of predatory gastropods on molluscan prey displayed significant geographic variation between north and south transects, seemingly related to the Po River delta as well. The vast majority of predatory drilling occurred north of the delta on seven species recovered in 12 occurrences from Lido di Jesolo to Grado (Figs. 2, 6). This pattern strongly contrasts with the one drilled species and its occurrences at Casal Borsetti. Farther south, predation traces slightly increase; however, as for trematodes, the different geological-morphological context and possible taphonomic biases (which need to be evaluated) cause us to interpret this trend with caution. In addition, here the small sample size of the species occurrences (which are scarcely represented in the northern sector; Fig. 6) make it difficult to evaluate the intensity of drilling (and could point toward taphonomic biases). Visually, the distributions of trematode pits and predatory drill holes are similar, but the summed prevalence/drill hole values were not significantly correlated. Nor were they correlated with spionids, for which drilling frequency has an almost inverse trend.
Sawyer and Zuschin (Reference Sawyer and Zuschin2010) reported significant differences in assemblage-level drilling frequencies from intertidal (1.4%) and subtidal (4–15 m depth; 27.4%) environments in the northeast Adriatic Sea. Although our summed drilling frequency metric is not directly comparable to Sawyer and Zuschin’s assemblage drilling frequency, our results are quite similar. Our strandline deposits are likely a mix of intertidal to shoreface environments.
The large-scale spatial trends in biotic interactions observed north and south of the Po delta seem strictly dependent upon the expected relationship between the structure of macrobenthic assemblages and the level of disturbance and fluvial discharge/supply. In a temporal perspective, this is influenced by base-level dynamics. Indeed, we can suggest that sediment-starved depositional environments characterized by relatively stable environmental parameters (e.g., salinity or sedimentation/turbidity) display higher prevalence and frequency values of trematodes and drilling traces; in contrast, spionids reach their maximum frequencies (especially as number of species infested) in less predictable depositional environments subject to strongly fluctuating environmental parameters. A possible confounding factor could be the differential time-averaging between the northern and southern sectors (<100 km from Po main channel). Indeed, the currently forming shell beds of the relatively starved northern sector should show greater time-averaging than those found south of the Po delta, whose time-averaging estimates are ≤200 years (Scarponi et al. Reference Scarponi, Kaufman, Amorosi and Kowalewski2013). Because time-averaging tends to average variations in faunal compositions, alpha diversity should increase, and with it chances to find higher richness (at species level) of biotic interactions in the north relative to the southern sector. However, it is unlikely that the higher prevalence/frequency of biotic interactions in the northern area is a simple product of differential faunal homogenization due to age mixing. First, age mixing occurred solely within relatively restricted time intervals, as transgressive deposits are well below the reworking zone (see above; Zecchin 2008); hence, the time of mixing is constrained within centuries to a few millennia. Second, higher single-species trematode prevalence values were found only in inaccessible TST sediments (e.g., C. gallina); thus, even if higher time-averaging could increase chances of finding infested taxa, higher time-averaging would lower the prevalence rates. Our estimates would therefore be conservative and not alter the species pattern.
Unfortunately, it is beyond the reach of this research to investigate either the relationship between biotic interactions and patterns of specific environmental variables or the dynamics of biotic interactions at smaller geographic scales (i.e., between sites). Our sampling scheme is not adequate to resolve/explain small-scale geographic patterns of biotic relations. Indeed, biotic variations between adjacent localities are influenced also by local factors (e.g., presence or not of minor rivers, local morphology, differential anthropic pressure, pollution) that can affect and drive small-scale geographic variability of analyzed biotic interactions. For example, trematode variations between adjacent localities could well be influenced by the spatial distribution of the definitive hosts of trematodes. The presence of these typically large-bodied, vagile, and low-abundance (relative to benthic intermediate hosts) vertebrates (whether bird, mammal, or fish) is required for the completion of the trematode life cycle. If they are locally absent or uncommon then trematode abundance and diversity should also be greatly reduced.
Biotic Interactions and Biodiversity
Hechinger et al. (Reference Hechinger, Lafferty, Huspeni, Brooks and Kurtis2007) suggested that abundance and diversity of trematodes within a local environment should be correlated with those of the co-occurring free-living benthic taxa. The reasoning was that many of the molluscan and arthropod taxa in the macrobenthos served as intermediate hosts for trematodes and were required for the completion of their life cycle. Their results from 32 sites in three wetlands in California and Mexico revealed significant positive correlations between trematode abundance/diversity and large-benthos abundance/diversity. The indication is that diverse and healthy ecosystems include parasites, so monitoring them is an inexpensive and reliable method for tracking ecosystem health. We expected to see similar patterns in the sub-fossil record from the northern Adriatic beaches of Italy, but a mixed approach is necessary for studying parasite traces preserved in dead bivalves. It is impossible to reconstruct the diversity of trematodes with our data, but one can estimate their abundance. This, of course, is true assuming that prevalence values based on traces are a predictor of prevalence values based on soft-tissue data, an assumption we are in the process of testing. We cannot estimate live host abundance because of the taphonomic bias inherent in strandline deposits, but we can estimate host diversity. The results of Hechinger et al. (Reference Hechinger, Lafferty, Huspeni, Brooks and Kurtis2007) suggest that summed trematode prevalence could be positively correlated with macrobenthos diversity. However, summed trematode prevalence was not correlated with either species richness (r=0.06 p=0.872) or rarefied diversity (r=0.17 p=0.627). This lack of correlation between host diversity and parasite abundance could (1) indicate a degraded ecosystem, (2) indicate an ecosystem in which parasitism is not a strong controlling factor of community structure, or (3) be the spurious result of taphonomic bias in that skeletal residues do not adequately capture the original biodiversity of an ecosystem.
In contrast, spionid polychaete-induced mudblisters are significantly negatively correlated with species richness (p=0.015) and even more so with rarefied diversity (p=0.009). Spionid U-shaped borings do not display the same patterns, but this is likely due to the smaller number of occurrences. Our results show that spionids thrive in low-diversity and, possibly, highly stressed environments. Indeed, spionids tend to be generalists that can tolerate a wide range of salinities, operate as suspension or deposit feeders, and thrive in environments that are organically enriched and prone to disturbance (Dix et al. Reference Dix, Karlen, Grabe, Goetting, Holden and Markham2005). These conditions are typical of those in the coastal environments influenced by the Po River where spionid traces are so common. The distribution of substrate type north and south of the Po delta could also influence the distribution of spionids, but this would be opposite to the pattern predicted from previous results (Huntley Reference Huntley2007).
Predatory drilling frequency was not correlated with any of the diversity indices. Predation frequency is often positively related to biodiversity, though the nature of thisrelationship is likely reliant on the spatial and temporal scale of observation (Paine Reference Paine1966; Hansen and Kelley Reference Hansen and Kelley1995; Vermeij Reference Vermeij2004; Huntley and Kowalewski Reference Huntley and Kowalewski2007). Some have characterized biotic interactions in the northern Adriatic as sluggish (Sköld and Rosenberg Reference Sköld and Rosenberg1996; Ceschia et al. Reference Ceschia, Zanchetta, Cestaro, Zambon and Sello2004; McKinney and Hageman Reference McKinney and Hageman2006; McKinney Reference McKinney2007; Huntley and Scarponi Reference Huntley and Scarponi2012), so perhaps because of the unusual ecological conditions we should not expect to find this relationship between predation and diversity. Kelley (Reference Kelley2006) reported drilling frequencies from the Lido di Venezia to be comparable to those from the Cretaceous of North America (lower than modern but higher than the Paleozoic). However, the work of Sawyer et al. (Reference Sawyer, Zuschin, Riedel and Stachowitsch2009) and Sawyer and Zuschin (Reference Sawyer and Zuschin2010) strongly suggest that predation intensity can be rather high in the northern Adriatic Sea, however patchy in distribution. Interestingly enough, Sawyer and Zuschin (Reference Sawyer and Zuschin2010) also failed to identify any correlation between diversity and predation intensity in the Gulf of Trieste.
Biotic Interactions and Molluscan Life Mode
Summed trematode prevalence is positively correlated with epifaunal mollusks. This result is unexpected, because only one host species out of six and three host species occurrences out of 15 are epifaunal taxa (Mytilus galloprovincialis). The remaining host species and species occurrences are semi-infaunal (Cerastoderma glaucum) or infaunal taxa (C. gallina, D. semistriatus, L. lucinalis, and S. subtruncata). Trematodes are negatively correlated with mollusks of restricted mobility and positively correlated with herbivores on rocky and/or coral substrates. These relationships are difficult to explain, and perhaps they are the result of ecologically, but not statistically, spurious correlations.
Summed spionid mudblister prevalence values are negatively correlated with epifaunal mollusks and positively correlated with suspension feeders. Many spionids are suspension feeders (and, in some cases, facultative deposit feeders), so it is likely much easier for a spionid to gain adequate nutrition when it has infested an epifaunal host, particularly one that is inducing a feeding current that the spionid can kleptoparasitize. Huntley (Reference Huntley2007) found that spionid-induced mudblisters occurred significantly more often in epifaunal individuals of Protothaca staminea living in a gravel-laden intertidal creek than in infaunal individuals living in a muddy-substrate lagoon near Friday Harbor, Washington, United States. However, the situation with our data is a bit different as the actual hosts of mudblister-forming spionids are not epifaunal but mostly shallow infaunal siphonate suspension feeders. Shallow infaunal siphonate suspension feeders do tend to be the most commonly available bivalves. Perhaps, though, this is the best type of host for a worm that can both suspension-feed and deposit-feed. Summed spionid mudblister prevalence values are also negatively correlated with the subsurface deposit feeders, chemosymbiotic deposit feeders, and herbivores on either rock or fine-grained sediment substrates. Subsurface and chemosymbiotic deposit feeders typically live in low oxygen substrates and do not induce feeding currents. Neither of these conditions is likely to make for well-fed spionids. Both types of herbivores retrieved here are typically gastropods whose shells are too small to host a spionid polychaete.
Summed prevalence of U-shaped borings is negatively correlated with the proportion of predatory carnivores, but one should be cautious and not overinterpret this lone significant correlation among U-shaped borings. U-shaped borings are rarer than the other traces (eight host species occurrences in five of 11 locations). Moreover, the same relationship is not found between spionid-induced mudblisters and predatory carnivores.
Summed predatory drilling frequencies are negatively correlated with sedentary mollusks. The proportion of sedentary mollusks is largely controlled by the presence of byssally attached epifaunal taxa. Prey individuals were mollusks of restricted mobility (82%), active mobility (14%), and immobility (4%). The negative correlation with sedentary mollusks is likely the result of diluted positive correlations with the other three mobility life mode categories. Summed predatory drilling frequencies are positively correlated with subsurface deposit feeders. Again, this is a difficult relationship to interpret as all prey individuals were suspension feeders (71%), chemosymbiotic deposit feeders (7%), or herbivores on plant/algal substrates (2%).
Interpreting Temporal Trends in the Light of Modern Spatial Variation
We found no significant difference between median value of trematode prevalence or shapes of distributions when comparing our data from modern northern Adriatic mollusks with the combined Holocene and Pleistocene trematode prevalence values among Po delta mollusks reported by Huntley and Scarponi (Reference Huntley and Scarponi2012). We should note that the sample sizes are rather small (n modern=15, n fossil=7) so some caution should be taken with interpretation. The Mann-Whitney U-test is rather robust to the vagaries of differences in variances; however, the Kolmogorov-Smirnov D-test requires at least eight observations within each group (Hammer and Harper Reference Hammer and Harper2006). The lack of a statistically significant difference in cumulative distribution shape could be due to small sample sizes beyond the resolving power of the Kolmogorov-Smirnov test, but we do accept the Mann-Whitney test results as robust. This spatio-temporal comparison supports the interpretation that trematode infestation is related both to the presence of a specific set of environmental conditions (see previous section) and to its dynamics in response to base-level changes.
In this specific case, trematode infestation was more widespread and intense in nearshore retrogradational deposits (TST) than in progradational comparable deposits (HST). The former were deposited in environments characterized by barrier-lagoon-estuary complexes and sediment-starved shallow marine environments that were widespread in the subsurface transgressive deposits of the central-southern part of the study area (see Wittmer et al. Reference Wittmer, Dexter, Scarponi, Amorosi and Kowalewski2014 and references therein); in contrast, barrier-lagoon-estuary complexes characterize only the northern sector of the study area in the modern (see also Fig. 2). That is, Huntley and Scarponi’s (Reference Huntley and Scarponi2012) temporal trend was not significantly different from modern spatial variation. It follows that the recurrence of comparable parasitism prevalence values through time and space suggests a stratigraphic paleobiological explanation.
The field area of Huntley and Scarponi (Reference Huntley and Scarponi2012) is located immediately south of the Po delta and the study mainly focused on the transgressive-regressive shallow-marine succession of the present interglacial. Here the relatively high and repeated occurrence of trematode infestation values through the transgressive deposits support the idea that environments amenable to trematodes (i.e., like modern Grado-Chioggia in Fig. 3) were more widespread spatially during the early middle Holocene. When the rate of relative sea-level rise slowed and the progradation took place (HST), sediment-starved oligotrophic (TST) environments vanished (or were greatly reduced) and trematode infestation was restricted from the local stratigraphic window. Hence, this spatial vs. temporal comparison suggests that the intense parasite prevalence pressure recovered in transgressive deposits is evidently not limited to this systems tract, but persists wherever there are TST-like environmental conditions where infestation intensity can be sustained (i.e., Grado-Chioggia sector of the North Adriatic).
The strong correspondence between trematode temporal (Po Plain) and spatial (northern Adriatic) dynamics in a well-defined sequence stratigraphic framework shed light on the stratigraphic mechanisms that shape trends of biotic interactions in the recent fossil record. It also highlights how shifting environments associated with base-level changes may bias the recording of long-term trends of antagonistic organismal interactions. Furthermore, this is true also for other biotic interactions analyzed here (i.e., spionid mudblisters) whose spatial distribution displays trends opposite to that of trematodes but are subject to the same constraining environmental factors (i.e., high sedimentary supply during progradation phases).
Taphonomic factors may influence the results of this fossil vs. modern comparison. Perhaps we cannot assume that the extent of time-averaging in fossil horizons is the same of that in modern beach-derived samples. Scarponi et al. (Reference Scarponi, Kaufman, Amorosi and Kowalewski2013) reported that currently forming upper shoreface shell beds and late Holocene beach ridges south of the Po delta record time-averaged assemblages up to several centuries, whereas time-averaging of fossil assemblages in comparable nearshore transgressive deposits is on centennial to millennial time scales. We do not know the extent of time-averaging in C. gallina; however, if our fossil and modern samples have experienced different degrees of time-averaging how might this affect trematode prevalence values? Huntley and Scarponi (Reference Huntley and Scarponi2012) illustrated that trematode infestation is sporadic in occurrence, but high prevalence values can be found in single species occurrences. If this is the case, then one would expect time-averaging to dilute trematode prevalence values in a sample. We see no significant difference in prevalence values between fossil and modern specimens, which could suggest minimal influence from differential time-averaging.
Conclusions
1. A total of 17,299 valves from 11 locations along the northern Adriatic coast of Italy were surveyed for traces of parasites, kleptoparasites, and predatory drilling gastropods.
2. Parasitic trematode-induced pits are more common and infest a greater diversity of host taxa north of the Po River delta, are mostly absent immediately south of the delta, and become more abundant >130 km south of the delta. These results support the hypothesis that trematodes are more common in starved, oligotrophic, and more stable or less stressed environments (TST-like conditions) than in high-disturbance fluvial-influenced settings subjected to high sedimentary input and dynamic salinity conditions.
3. Kleptoparasitic spionid-induced mudblisters and U-shaped borings are nearly absent north of the Po River delta and highly abundant south of the delta. Spionid traces are associated with lower-diversity assemblages and are typically found on shallow infaunal suspension-feeding bivalves.
4. Predatory gastropod drill holes were mostly found in prey taxa far north of the Po River delta in association with TST-like barrier island/lagoon/estuary environments. Predation intensity is not correlated with species richness, evenness, or rarefied diversity values. Predation intensities in our shoreface deposits fall between the low intertidal drilling frequencies and elevated sublittoral frequencies reported by Sawyer and Zuschin (Reference Sawyer and Zuschin2010), were generally higher than those reported by Kelley (Reference Kelley2006), and do not fully support the hypothesis that predation intensity is sluggish in the northern Adriatic.
5. There was no significant difference in trematode prevalence values between our modern samples and the fossil samples of Huntley and Scarponi (Reference Huntley and Scarponi2012). This reinforces the interpretation of Huntley and Scarponi (Reference Huntley and Scarponi2012) that the significant temporal shifts in trematode prevalence among Pleistocene–modern bivalves in the Po Plain/northern Adriatic are the result of environmental turnover associated with sea-level changes. These results highlight the importance of disentangling stratigraphic bias from biological signal when interpreting temporal trends.
Acknowledgments
The authors gratefully acknowledge the thoughtful reviews provided by M. Casey and S. Walker. P. De Muro and L. A. Huntley assisted with fieldwork. This project was funded by a Postdoctoral Research Fellowship and European Research Visit Grant from the Alexander von Humboldt Stiftung and a University of Missouri Research Council Grant (J.W.H.) and funding from RFO 2012 (D.S.). Fieldwork was made possible by a generous Ford Ka grant from the Corey and Heather Long Foundation (J.W.H.).
Supplementary Material
Supplemental materials deposited at Dryad: doi:10.5061/dryad.rt31g


