Responses of ecosystems to environmental changes vary greatly across habitats, organisms and observational scales. The Quaternary fossil record of the Po Basin demonstrates that marine communities of the northern Adriatic re-emerged unchanged following the most recent glaciation, which lasted approximately 100 000 years. The Late Pleistocene and Holocene interglacial ecosystems were both dominated by the same species, species turnover rates approximated predictions of resampling models of a homogeneous system, and comparable bathymetric gradients in species composition, sample-level diversity, dominance and specimen abundance were observed in both time intervals. The interglacial Adriatic ecosystems appear to have been impervious to natural climate change either owing to their persistence during those long-term perturbations or their resilient recovery during interglacial phases of climate oscillations. By contrast, present-day communities of the northern Adriatic differ notably from their Holocene counterparts. The recent ecosystem shift stands in contrast to the long-term endurance of interglacial communities in face of climate-driven environmental changes.
Anthropogenic climate and environmental changes present increasing threats to marine ecosystems and their biodiversity. Diverse human-induced stressors impact marine life at multiple scales by impairing physiology of individual organisms, disturbing ecological interactions, modifying or nullifying biogeographic barriers, and inducing extirpation and extinction events [1–3]. Resulting changes can lead to ecosystem degradation  and the emergence of new ecosystem states [2,5,6]. A historical perspective is a prerequisite for calibrating the significance of anthropogenic changes [7–11].
Climate oscillations that took place during the Quaternary offer an opportunity to quantify long-term responses of marine ecosystems to naturally occurring environmental changes and provide a reference system for assessing the severity of anthropogenic processes. However, the Quaternary fossil record has provided mixed answers regarding the magnitude of marine ecosystem responses to natural environmental change. For example, Pleistocene reef communities appear to have remained stable in many areas of the world [12,13] and Quaternary molluscs of coastal California suffered unexpectedly low extinction rates . By contrast, differential biogeographic effects could be observed across mollusc species with different body size  and the western Atlantic mollusc faunas experienced high species turnover during the Quaternary . Here, using the rich fossil record of marine benthic communities from coastal and deltaic habitats of Po Basin (Italy), the ecosystem changes across two consecutive interglacial time intervals have been evaluated quantitatively and contrasted against comparable data available for present-day communities of the northern Adriatic.
2. Material and methods
(a) Study area
The distal part of the Po Basin is filled with alternating packages of Quaternary continental and shallow-marine deposits (figure 1) that record interplay between recurrent climate-driven oscillations in sea level and changes in local sediment supply . The penultimate (Late Pleistocene) package of marine deposits represents the previous interglacial highstand (marine isotope stage 5e), which lasted from approximately 130 000 to 116 000 years ago . During the subsequent glacial phase, the regional sea level dropped by about 120 m and the Adriatic Sea retreated approximately 270 km southeastward . This prolonged, stepwise phase of sea-level fall lasted from approximately, from 116 000 to 19 000 years ago and is recorded proximally by a thick package of predominantly alluvial deposits, with intervening lagoonal and paludal clays (figure 1). Renewed sea-level rise took place between 15 000 and 14 000 years ago in response to the incipient melting of the West Antarctic Ice Sheet . In the Po delta region, the resulting youngest marine sequence records the last approximately 10 000 years : the Holocene. The Late Pleistocene and the Holocene marine sequences, recorded in the subsurface of the Po coastal plain roughly at 0–30 and 95–120 m core depth intervals (figure 1), are nearly identical in terms of facies architecture: each represents a wedge-shaped succession of transgressive (TST) and highstand (HST) systems tracts dominated by coastal, shallow-marine and deltaic deposits (electronic supplementary material, data S1).
(b) Sample processing and data filtering
The fossil dataset of mollusc specimens was obtained from 16 cores drilled in the Po Plain (figure 1). A total of 611 bulk samples (approx. 375 cm3) were processed to create a sample-level, species-level data matrix (333 species and 131 780 specimens). The data were filtered by removing specimens that could not be identified to a unique species and excluding samples from non-marine habitats. For those analyses in which sample sizes and species rarity can generate substantial volatility, data were filtered further by removing small samples (n < 20) and species occurring in one sample only (see figure/table captions, and the electronic supplementary material, data S1). Because only some of the cores penetrated into the older marine sequence, the Late Pleistocene dataset is notably smaller than the Holocene dataset. After excluding non-marine samples and unidentified species, the final dataset (electronic supplementary material, table S1) included 221 species, 125 558 specimens and 453 samples (414 for Holocene and 39 for Pleistocene).
The dataset for Modern communities was derived from a compilation of marine benthic surveys of the northern Adriatic . These data were restricted to samples that match the environmental and bathymetric range estimated for the fossil dataset derived from cores. The Modern samples include counts of dead specimens to make them comparable to time-averaged samples derived from cores (electronic supplementary material, data S1). The Modern dataset includes 78 samples, 91 species and 91 552 specimens (electronic supplementary material, table S1). To make the Holocene dataset maximally analogous to the Modern dataset, which represents HST settings only, the subset of the Holocene dataset restricted to HST samples only was also derived (the Holocene-HST dataset includes 263 samples and 129 species; electronic supplementary material, table S1).
(c) Analytical methods
Multiple indirect ordination techniques, applied previously to explore relationships among the core samples in terms of their faunal composition, yielded consistent results . Only non-metric multidimensional scaling (NMDS) is reported here (electronic supplementary material, figure S1). The NMDS ordination demonstrated that the preferred bathymetry of species (independently estimated for the region using present-day ecological data) correlated strongly with species score on the first ordination axis [23,24]. Similarly, the first axis scores appear to be a robust predictor of water depth for the core samples . These water depth proxies  were used here as estimators of water depth for each core sample (figure 3; electronic supplementary material, figure S1). The bathymetry of samples included in the Modern dataset was estimated directly.
Diversity analyses included the sample-level (α) standardized diversity, evenness, the turnover β-diversity (where pairwise Bray–Curtis distances between samples were evaluated as a function of their environmental distance estimated by difference in sample water depth), and gamma diversity estimated separately for each dataset. The Late Pleistocene, Holocene and Modern rank abundance distributions (RADs) were evaluated in terms of commonly invoked theoretical RAD models. The RADs were also compared using non-parametric methods (Kolmogorov–Smirnov D and Spearman rank ρ).
The effects of unbalanced sampling were evaluated using four resampling models. These models simulate the expected outcomes under the null hypothesis that the evaluated datasets represent the same ecosystem. The randomization model (RDM) was based on random reassignment of sample membership, where Late Pleistocene and Holocene samples were pooled and then randomly reassigned without replacement to either ‘Holocene’ or ‘Pleistocene’. The same protocol was implemented to compare Holocene-HST and Modern datasets. The subsampling model (SSM) was based on random subsetting of Holocene samples, where only Holocene (or Holocene-HST) samples were used to create random subsets that mimic the Late Pleistocene (or Modern) dataset. The adjusted randomization model (ARDM) was derived by modifying the output of the RDM simulation to correct for both differences in the number of samples and the mean sample size. To further explore the sensitivity of estimates to the choice of the resampling methodology, a double standardization model (DSM) was implemented to standardize for both the number of samples and the mean sample size. In this two-step approach, all datasets were first subsampled with replacement down to 39 samples: the size of the smallest of the compared (Late Pleistocene) datasets. In the second step, all large samples included in the random set of 39 samples, were subsampled without replacement to make mean sample sizes comparable across the datasets. DSM was used to compare the DSM subsampled datasets (Holocene, Holocene-HST and Modern) to the Late Pleistocene dataset (the comparison of the Late Pleistocene dataset to its resampled replicates served as a null model). All models attempt to mimic a homogeneous system with no extirpations, originations (new species arrivals) or changes in species abundance. See the electronic supplementary material, data S1, for further information regarding the study area, sampling protocols, analytical details and methodological justifications.
The faunal composition of samples appears to be a strong correlative of water depth: multivariate ordination scores calibrated using present-day ecological data provide bathymetric estimates for each sample to less than 3 m [23,24]. The Holocene and Late Pleistocene samples cover comparable bathymetric ranges, from lagoon/shoreface to offshore transition, although the Late Pleistocene samples, owing to less extensive core coverage, represent a slightly narrower depth range than the Holocene samples (electronic supplementary material, figure S1). Depth range represented by samples included in the Modern dataset was restricted to match the range estimated for the Holocene dataset.
The total species diversity is much higher in the Holocene than in the Late Pleistocene (electronic supplementary material, table S1). However, the resampling models (electronic supplementary material, table S2) suggest consistently that the observed Pleistocene species diversity (97 species) does not differ significantly from resampling estimates predicted for a single system: (RDM: mean = 85.2; median = 85; lower quartile Q1 = 77; upper quartile Q3 = 93; p = 0.167). The SSM and ARDM models approximate closely all RDM estimates (electronic supplementary material, table S2 and figure S2). Only 11 out of 97 Pleistocene species were not found in the Holocene (3.4% of all species). These 11 apparently extirpated species were all exceedingly rare in the Late Pleistocene (electronic supplementary material, table S3). One of those species is an extant freshwater mollusc represented by one specimen that was probably transported post-mortem from more proximal habitats. The remaining species are still present in the Adriatic Sea today (electronic supplementary material, table S3). The applicable models (RDM and ARDM) predicted that, on average, six species (Q1 = 4 and Q3 = 8 for both models) were expected to be missing under the null model of 0% extirpations and the observed number of apparent extirpations was not significantly greater (p = 0.103 for both models). Out of 210 species reported from the Holocene samples, 59% were missing from the Late Pleistocene samples (124 apparent species originations). The models consistently suggested that at least as many species originations should be observed due alone to the undersampling of the Pleistocene time interval (RDM: mean = 135.8; median = 128, Q1 = 90, Q3 = 144, p = 0.167; electronic supplementary material, table S2).
By contrast, comparisons of Holocene-HST and Modern datasets indicate that observed differences depart notably from those predicted by the models. The number of observed originations (o = 38) exceeds significantly the RDM predictions (mean = 11.0, median = 11, Q1 = 8, Q3 = 13, p < 0.0001). Also, RDM predicts 64.4 species extirpations (Q1 = 59, Q3 = 70), whereas 76 Holocene-HST species are absent in Modern dataset, a nearly significant difference (p = 0.0760). Moreover, the apparently extirpated species are not limited to extremely rare species (as was the case for Pleistocene–Holocene comparison). The most abundant of the apparently extirpated Holocene-HST species (Bittium reticulatum) is represented by 451 specimens. RDM indicates that a species with such high abundance is unlikely to be missing from the Modern dataset just due to sampling (mean = 47.42, median = 24, Q1 = 18, Q3 = 95, p = 0.0002). Whereas this species is extant in the northern Adriatic, it is rarely observed in the study area. Similarly, the mean abundance of apparently extirpated Holocene-HST species (15.1 specimens) exceeds significantly RDM predictions (mean = 2.53, median = 2.25, Q1 = 1.82, Q3 = 3.27, p = 0.0001).
The RDM and SSM models may be biased owing to the fact that mean sample size varies substantially across the datasets (electronic supplementary material, table S2): it is highest for the Modern dataset, intermediate for the Pleistocene dataset, and smallest for the Holocene and Holocene-HST datasets. However, the ARDM model which attempts to correct for this bias by excluding iterations where this bias is strongly manifested (note the correct approximation of mean sample sizes of the actual datasets for ARDM in electronic supplementary material, table S2), yielded consistent (often numerically identical) outcomes (electronic supplementary material, table S2). The DSM model, which represents a different standardizing strategy, yielded congruent results. Namely, when all datasets are compared to the Late Pleistocene dataset, notably higher origination and extirpation rates are estimated for the Modern dataset relative to the Holocene and Holocene-HST datasets (electronic supplementary material, figure S3).
The same species (Lentidium mediterraneum) dominates the Holocene and Late Pleistocene samples, and its overall dominance is nearly identical for the two time intervals: Late Pleistocene (66.3%) and Holocene (65.9%). Eight out of the top 10 species recovered from the Late Pleistocene samples also belong to the top 10 species recovered from the Holocene samples (figure 2) and the species rank abundances of Holocene and Late Pleistocene datasets are positively and significantly correlated (r = 0.58; p < 0.0001; figure 2, inset). Maximum-likelihood estimates of the species RADs indicate that the Zipf model offers the best fit to both the Holocene and Late Pleistocene distributions (figure 2), using either Akaike information criteria or Bayesian information criteria. The Zipf coefficients are nearly identical for the two datasets (electronic supplementary material, table S4), and the value observed for the Late Pleistocene molluscs approximates the expected values predicted by the resampling models (electronic supplementary material, figure S2 and table S2). Because the Zipf model is a descriptive model , its fit to the data should not be interpreted as evidence for structuring importance of biotic interactions or a neutral mechanism (i.e. multiple processes can conform to the same model). However, RADs that fit the Zipf model tend to be interpreted as representing relatively more complex ecosystem structures . The two RADs are indistinguishable statistically (Kolmogorov–Smirnov D = 0.083, p = 0.750). In fact, the observed D-value is notably lower than predicted by resampling (electronic supplementary material, table S2 and figure S2), suggesting that observed RADs are more similar to each other than would be expected by chance. This apparent hyper-homogeneity is not significant and may simply reflect resampling process (in many iterations, one of the two compared randomized datasets will omit one of the few samples dominated by the most abundant species resulting in pairs of RAD curves that yield higher D-values than observed in the actual data). Conversely, only a tiny fraction of replicate datasets produced D-values that are comparable to or lower than the observed value. In contrast to Holocene–Pleistocene comparisons, Modern samples differ in multiple ways from the most comparable subset of core samples (Holocene-HST). The RAD distribution is best fit by lognormal model (figure 2; electronic supplementary material, table S4) and differs significantly from the Holocene-HST RAD (D = 0.449, p < 0.0001). The two datasets appear uncorrelated in terms of species abundance (r = 0.085, p = 0.274). The dominant species in Modern communities is different (Varicorbula gibba) and represents only 30.6% of the total specimen abundance.
The turnover β-diversity  can be quantified for each time interval in terms of bathymetric distance between pairs of samples. Pairwise Bray–Curtis similarities revealed comparable bathymetric gradients for Holocene and Pleistocene sample pairs (figure 3a). The decay-distance exponential models  indicate equivalent turnover rates (figure 3c) with similar strength of relationship: Holocene (λ = −0.23, R2 = 0.33) and Pleistocene (λ = −0.22, R2 = 0.29). The relative frequency distributions of Bray–Curtis similarity values are also congruent for Holocene and Pleistocene sample pairs (figure 3b). By contrast, Modern sample pairs show a much less pronounced gradient, pointing to relatively higher within-habitat heterogeneity (figure 3) and notably lower turnover rates (λ = −0.12, R2 = 0.08). Consistent with these results, NMDS ordination indicates a complete overlap of the Pleistocene and Holocene samples (electronic supplementary material, figure S1), and only partial overlap for Holocene-HST and Modern samples (electronic supplementary material, figure S4). Finally, permutation-based methods (ANOSIM, PERMANOVA and Mantel), used here comparatively (and not as statistical tests), all consistently suggest (electronic supplementary material, table S5) that multivariate differences between Holocene-HST and Modern datasets are notably greater than those between Late Pleistocene and Holocene datasets.
When core samples are plotted along the depth gradient, the standardized species diversity, per-sample specimen abundance and dominance reveal comparable patterns for Holocene and Late Pleistocene datasets: diversity reaches its maximum between 13 and 8 m water depth, whereas dominance and abundance peak at approximately 3 m (electronic supplementary material, figure S5). In addition, bathymetric distributions of the three most abundant Late Pleistocene species are replicated in the Holocene, including the observed preferred depth and relative specimen abundance (electronic supplementary material, figure S6).
When interpreted in the context of the resampling models, the results indicate that the Holocene coastal ecosystem is indistinguishable from its Late Pleistocene counterpart, including γ diversity, species rank abundance, β turnover diversity, bathymetric distributions of the most common species, and spatial trends in sample-level diversity, dominance and specimen abundance. By contrast, comparisons with the Modern dataset suggest that modern communities differ notably, and often significantly, from their Holocene and Late Pleistocene counterparts.
4. Potential confounding factors
Several spurious factors may have potentially contributed to the observed patterns. First, the analysed core samples are time-averaged over centennial-to-millennial time scales  potentially resulting in homogenization of core samples, either across or within time intervals. However, the similarities between the two time intervals cannot be attributed to reworking of Pleistocene fossils into Holocene sediments. A thick package of alluvial sediments (figure 1) separates the two marine sequences and all individually dated mollusc specimens from the Holocene samples lived during the last 10 000 years . Second, individual samples vary predictably in faunal composition indicating that spatial heterogeneities within each time interval, primarily a correlative of bathymetry, can be recovered despite time-averaging. This is consistent with live–dead comparative studies in modern ecosystems, which also suggest that time-averaged shell accumulations record ecological and bio-inventorying information . Third, because the Modern samples are dominated by shells of dead specimens dredged from the uppermost sediment layer, they represent time-averaged assemblages comparable to those recovered from the core samples. Fourth, time-averaging varies across the sampled transgressive–regressive cycles, as demonstrated by direct dating of mollusc shells . However, the Late Pleistocene–Holocene comparisons are based on comparable data (61.5% and 63.5% of HST samples, respectively). Similarly, the Holocene-HST–Modern datasets are restricted to HST samples only. Fifth, the observed difference between the Holocene-HST and Modern datasets could reflect selective taphonomic loss of small, thin-shelled, aragonitic specimens owing to compaction and early diagenetic dissolution. This is unlikely given that core samples are dominated by small, thin-shelled, aragonitic taxa. Finally, the Po Plain dataset is limited to shelly molluscs, and thus, provides an incomplete picture of original communities. However, molluscs represent a substantial fraction of macrobenthic communities of the present-day northern Adriatic  and multi-taxon tests in other regions suggest that molluscs may be a reliable proxy for entire macrobenthic communities .
Neither the recurrence of the interglacial communities nor the subsequent ecosystem shift can be explained by differences in sampling coverage, taphonomic filters, differential time-averaging or exclusion of non-mollusc faunas. Both the stability of the communities and their subsequent shift are consistently demonstrated using various metrics that are expected to be differentially affected by confounding factors. The high concordance of those metrics is more parsimonious with ecosystem changes than spurious drivers.
5. Discussion and conclusion
High congruence of the consecutive interglacial mollusc associations of the Po Plain may reflect long-term resilience (i.e. reassembly of analogous coastal communities during successive highstands) or long-term persistence (i.e. back-and-forth migrations of stable communities in concert with sea-level changes) of the regional coastal ecosystem. Those two rival hypotheses cannot be evaluated until adequate data for the lowstand coastal ecosystems of the northern Adriatic become available. However, as is the case for many regions of the world , the fossil record of coastal ecosystems from the last glacial phase is difficult to access. It has been drowned offshore in the Mid-Adriatic Deep (water depths >100 m) . Regardless of which of those two hypotheses is correct, the results indicate that interglacial benthic communities of Po Plain either endured or recurred despite long-term, climate-driven environmental changes. Large-scale sea-level changes were not enough to perturb the system permanently.
These results are consistent with those reported previously for the Caribbean reef communities  and reinforce some of the previous mollusc-centred studies, which were either based on occurrence data  or focused on short-term, small-magnitude climate oscillations . The results also parallel, to a various extent, long-term analyses of marine communities in the older fossil records , including empirical examples of comparably strong correlations in shared species rank abundance (r > 0.5) observed across successive stratigraphic sequences . Similarly, recurrence of faunal gradients has been demonstrated in the fossil record [32,34], although the rank order of taxa did not necessarily repeat through time as faithfully as documented here. However, the community recurrence should not be necessarily equated with coordinated community responses [32,35]. The results contrast with those Quaternary studies that document substantial changes in marine ecosystems during long-term climate oscillations [16,36]. Similarly, since the Late Pleistocene, many terrestrial communities have experienced differential responses of species or species guilds, substantial ecosystem shifts or homogenization [37–39]. The results thus reinforce the notion that ecosystem responses to naturally occurring climate changes may vary across ecosystem types .
The ecosystem shift observed between the Holocene-HST and modern communities suggest a relatively recent timing of those changes. Whereas the base of the highstand (Holocene-HST) succession dates back to approximately 5 thousand years ago , a majority of the samples are substantially younger. Owing to accelerating net accumulation rates during the highstand, younger samples are over-represented in the sampled cores . In fact, based on numerical dating of mollusc shells from the sampled cores , as much as 90% of specimens in the Holocene-HST dataset represent the last millennium and 74% represent the last 500 years (electronic supplementary material, table S6). The observed ecosystem shift is thus unlikely to have started more than a few centuries ago.
Whereas specific processes behind the recent changes cannot be readily evaluated given the temporal resolution of our data, various ecosystem changes have been documented in the region previously and linked explicitly to diverse anthropogenic processes . The observed timing is also consistent with relatively recent shifts observed for mollusc-dominated communities in other regions, as suggested by discordances between live communities and surficial death assemblages [41,42].
The timing is relevant for assessing the regional timeframe of the Anthropocene (i.e. the onset of notable anthropogenic changes). The beginning of the Anthropocene is typically placed around the end of the eighteenth century, coincident with the Industrial Revolution [43–45]. However, other authors suggested that anthropogenic changes may have begun in the Early Holocene  or even the latest Pleistocene , and a compelling case has been made for the pre-industrial roots of anthropogenic alterations in the global marine ecosystem . Although major anthropogenic alterations of the region date back to the Roman times , these changes have accelerated most dramatically in the last two centuries . The comparative results presented here reaffirm those assessments. However, the earlier onset of finer-scale changes for specific marine species of human interest  is not precluded by the ecosystem-scale analyses presented here.
Marine benthic communities of the northern Adriatic changed recently, after withstanding environmental fluctuations for at least 100 000 years. This shift is consistent with previous claims that anthropogenic restructuring of marine communities is not a trivial extension of pre-human changes  and the ongoing resetting to a new ecosystem state may last beyond the current interglacial warming. Because the post-industrial (late global market) intensification of anthropogenic impacts in the region  partly postdates the northern Adriatic surveys of modern communities used here, the ecosystem shift documented here may have escalated further in the most recent decades.
Data available from the Dryad Digital Repository: http://dx.doi.org/10.5061/dryad.t7g11.
The work was funded by the National Science Foundation (grant no. EAR-0920075) and the Jon A. and Beverly L. Thompson Endowment Fund (University of Florida).
The authors thank Raffaele Pignone (Geological, Seismic and Soil Survey of Regione Emilia Romagna) and Giovanni Gabbianelli (Bologna University) for providing access to some of the cores used in this study. Three anonymous reviewers provided constructive and illuminating criticisms that improved greatly the quality of this report.
- Received December 8, 2014.
- Accepted January 19, 2015.
- © 2015 The Author(s) Published by the Royal Society. All rights reserved.