The results demonstrate a repeated association between significantly elevated prevalence and centennial scale flooding events, support the link between sea-level rise and increasing trematode activity, and can serve as historical analogues for ongoing and future anthropogenic climate change. Despite evidence for the consistent relationship between transgression and trematode prevalence, it is doubtful that a relative rise in sea level alone drove this pattern. Many factors that can influence the biota, including temperature, nutrient availability, salinity, host availability, diversity, and community structure, co-vary with sea level changes and should be tested as driving factors32, 33. Increasing temperature has been shown to increase reproductive output and infectivity of a diverse array of pathogens and parasites3, 4 (but see refs 7, 8 and 34). As parasites derive nutrition from their hosts, it is not clear that changes in nutrient availability/productivity would directly control their distribution, however biological diversity is often related to productivity and its mode of delivery across a variety of scales and systems35, 36. Diversity and productivity often increase in concert until a tipping point above which diversity begins to decline, varying with the influence of consumers and disturbance level37. In this way productivity could control the distribution and abundance of many taxa that might serve as intermediate or definitive hosts, though likely in a non-linear manner. Salinity is a primary environmental driver of mollusk turnover in the studied system (Fig. 2B), and free-swimming larval (cercaria) production and survival time tend to decrease significantly in lowered salinity regimes in paralic environments; thereby reducing infestation of intermediate and/or definitive hosts38, 39. Here, however, as in ref. 23, the lack of correlation between salinity proxy and prevalence estimates (Table 1) suggests that salinity is not a strong driving factor of trematode prevalence at this spatial and temporal scale of observation.

Table 1 Spearman rank correlation coefficients and p-values (when p < α = 0.05; otherwise indicated as ns: non-significant) between arcsine-transformed trematode prevalence values of Abra segmentum from the 204-S7 core; nMDS1 sample scores; and environmental, ecological, and taphonomic variables. nMDS1: Non metric Multidimensional Scaling axis 1. Full size table

The absence of correlation between preferred host (A. segmentum) abundance and prevalence (Table 1) rules out fluctuating host availability as a limiting factor of trematode distribution23. The median shell length of infested valves of A. segmentum was significantly larger than that of their non-infested counterparts (Mann-Whitney U, p = 2.21E-34), likely due to the accumulation of parasites through ontogeny. Prevalence values were positively and significantly correlated with host shell length (r = +0.68, p = 0.004), however there were no significant associations between shell length and either flooding pulses or nMDS1 scores (Fig. 2C and Table 1). This suggests that other environmental or ecological factors, acting as drivers of host shell length, were unlikely to have indirectly driven the temporal trend of trematode prevalence. Similarly, the lack of correlative relationships between prevalence and standardized richness, dominance, and Shannon diversity (Table 1) suggests that fluctuating biodiversity did not exhibit direct/linear control over trematode-bivalve interactions.

The role of more complex, community-level factors that may have influenced the distribution of trematode parasites can be examined by evaluating the distribution of samples and their constituent taxa in the nMDS space and assessing faunal similarity using Bray-Curtis pairwise comparisons to measure faunal turnover throughout the length of the core (Fig. 3A). Samples retrieved from brackish muds only (8.50–12.25 m core depth) display a comparable amount of turnover to that identified when comparing samples from both freshwater and brackish environments. However, the dendrogram derived from the Q-mode cluster analysis of the samples included in the nMDS ordination demonstrates that samples recording elevated or subdued trematode prevalence were distributed haphazardly across the dendrogram topology. Consequently, community structure/turnover (Fig. 3B) is unlikely to have been a driving factor of trematode prevalence within Holocene lagoonal facies.

Figure 3 Turnover and ecological similarity of assemblages across core 204-S7. (A) Turnover estimated by pairwise comparison of Bray-Curtis similarity indices and environmental distance (nMDS1 salinity). Solid circles indicate pairwise comparisons between lagoonal muddy samples from core depths of 8.50–12.25 meters and the black line indicates the ten point running average. Hollow circles indicate all other comparisons. The red line indicates the ten point running average for all comparisons. (B) Q-mode cluster analysis (UPGMA algorithm, Bray-Curtis similarity). Samples with trematode prevalence values of Abra segmentum greater than and less than the 95% CI are indicated in red and blue, respectively. Note how the samples of either exceptionally low or high prevalence values are randomly distributed across the dendrogram. Full size image

Sample nMDS1 scores were negatively correlated with the proportion of Abra valves that were too fragmented to evaluate in terms of infestation status. This pattern raises the question of how fragmentation might have influenced the parasite record (i.e., were valves with trematode pits more prone to fragmentation than non-infested valves?). All Abra valves were classified as either whole or broken, and the broken valves were further categorized into either “sufficiently complete” or “too fragmented” to determine infestation status. There was no significant difference in trematode prevalence values of whole and “sufficiently complete” broken valves (Χ 2, p = 0.16). These results suggest that the proportion of “too fragmented” valves was unlikely to represent an important confounding factor in reconstructing the stratigraphic record of trematode dynamics.

Another potential factor affecting parasite prevalence is the fluctuating availability of habitat-area for trematodes during sea-level cycles. The geologically rapid creation of new habitat during flooding pulses and their subsequent destruction during progradation could exert a first order control on trematode prevalence during high frequency cycles. As sea level continues to rise, some settings will be more strongly influenced than others. For instance, densely populated lowlands, estuarine, and riverine settings would likely display the greatest increase in trematode habitat-area during relative sea level rise as a direct effect of flooding and, indirectly, by the landward rise of the groundwater table40. Therefore, we hypothesize that gymnophallid trematode prevalence will be more strongly influenced by the creation of new habitat in brackish and freshwater settings than in shallow marine settings. Though not the direct topic of research here, an increase in wetlands created by sea level rise would generate new habitat for the gastropod intermediate hosts of Schistosoma 41, the trematodes responsible for schistosomiasis in humans.

The fossil record of the northern Adriatic points to a significant association between the prevalence of heterocious parasites and flooding events recording repeated climate-driven sea level shifts. From this historical perspective we posit that the ongoing anthropogenic warming and sea-level rise should trigger a significant upsurge in gymnophallid trematode prevalence and the expansion of wetland habitats ideal for schistosomatid intermediate hosts. The forecasted changes are expected to suppress the fecundity of common benthic organisms, exert negative impacts on ecosystems, impede ecosystem services, and, eventually, negatively affect human well-being.