Non-technical Summary
To better understand how ecosystems react to long-term climate changes, we can examine the Late Quaternary fossil record, a period characterized by strong climate fluctuations. Although much research has focused on species presence, it is also crucial to recognize changes in functional traits—biologically and ecologically important characteristics of organisms—because these traits influence ecosystem functions. We examined how molluscan functional traits changed across the latest Quaternary climate shifts in fluvially influenced nearshore settings of the Adriatic Sea. Frequencies of different traits were similar during the last and current warm periods (interglacials) but distinct in the last cold period (glacial). Trait diversity was higher during the last glacial, likely because lower temperatures reduced the dominance of species common in interglacials, allowing less frequently occurring species to thrive, and higher habitat diversity leading to more varied traits. The similarity between past and current interglacials suggests that in the Adriatic nearshore systems, moderate warming comparable to conditions during the last interglacial may not necessarily lead to major changes in ecosystem functioning. However, human activities such as seafloor disturbance, nutrient input, and the introduction of nonindigenous species could affect ecosystem functioning to a higher degree. By studying past trait responses to natural climate change, we established a baseline to assess the magnitude of the ongoing and future human impacts on the functioning of the shallow-marine ecosystems of the Adriatic Sea.
Introduction
Climate change is considered one of the major factors shaping marine biodiversity. While ecological monitoring, experimental data, and modeling studies can be useful to predict short-term consequences, they are insufficient to understand the long-term responses of entire ecosystems to climatic conditions that have not been experienced during the limited time frame of real-time scientific observation. In this respect, paleontological data can provide ample information on past ecosystem responses to climate perturbations (Dietl et al. Reference Dietl, Kidwell, Brenner, Burney, Flessa, Jackson and Koch2015; Kiessling et al. Reference Kiessling, Smith and Raja2023; Finnegan et al. Reference Finnegan, Harnik, Lockwood, Lotze, McClenachan and Kahanamoku2024; Malanoski et al. Reference Malanoski, Finnegan, Huang, Blake, Niocaill and Saupe2026). The Quaternary fossil record is dominated by extant taxa and therefore provides numerous examples of how climate change has affected marine ecosystems (e.g., Willis and Birks Reference Willis and Birks2006; Dietl et al. Reference Dietl, Kidwell, Brenner, Burney, Flessa, Jackson and Koch2015; Kowalewski et al. Reference Kowalewski, Nawrot, Scarponi, Tomašových and Zuschin2023). The Late Quaternary is characterized by geologically rapid fluctuations (~100 kyr cycles) between colder glacials and warmer interglacials (Pisias and Moore Reference Pisias and Moore1981; Blackburn et al. Reference Blackburn, Kodama and Piccione2024), which makes it a particularly interesting period for evaluating the natural range of variability in marine communities and their resilience to long-term environmental change.
The Late Quaternary succession of the Po-Adriatic system, located in the central-northern Mediterranean, has a well-understood stratigraphy and paleoenvironmental history (Trincardi et al. Reference Trincardi, Campiani, Correggiari, Foglini, Maselli and Remia2014; Wittmer et al. Reference Wittmer, Dexter, Scarponi, Amorosi and Kowalewski2014; Amorosi et al. Reference Amorosi, Maselli and Trincardi2016; Zecchin et al. Reference Zecchin, Donda and Forlin2017; Campo et al. Reference Campo, Pellegrini, Sammartino, Trincardi and Amorosi2024 and references therein). During the latest glacial–interglacial oscillations, sea-surface temperatures (SSTs) in the central Mediterranean differed strongly between glacials and interglacials. Available reconstructions suggest that last interglacial (Marine Isotope Stage [MIS] 5e) SSTs were locally ~2–3°C warmer than estimates for the Holocene Climate Optimum (HCO) (Siani et al. Reference Siani, Magny, Paterne, Debret and Fontugne2013; Santagati et al. Reference Santagati, Perri, Bernasconi, Borrelli, Guerrieri and Critelli2024; De Santis et al. Reference De Santis, Cerrone, Meschis, Scicchitano, Ascione and Caldara2025), although the data of Piva et al. (Reference Piva, Asioli, Andersen, Grimalt, Schneider and Trincardi2008) indicate even larger offset. By contrast, terminal last-glacial SSTs were at least 6°C lower than during the HCO (Piva et al. Reference Piva, Asioli, Andersen, Grimalt, Schneider and Trincardi2008; Bazzicalupo et al. Reference Bazzicalupo, Sicre, Checa, Maiorano, Margaritelli, Klein and Pena2022). Apart from paleoclimatic changes, the study area was also strongly affected by changing basin morphology. During the last glacial period, large parts of the present-day Adriatic, especially the shallow northern region, were subaerially exposed. Consequently, major river systems, such as the Po River, experienced extensive reorganization (e.g., Maselli et al. Reference Maselli, Trincardi, Asioli, Ceregato, Rizzetto and Taviani2014; Amorosi et al. Reference Amorosi, Maselli and Trincardi2016). During this cold period, fluvially influenced marine environments were characterized by high sedimentation rates, eutrophic waters, and lower salinity due to frequent freshwater inflows into the basin, which was at that time confined to the central and southern parts of the present-day Adriatic (Asioli et al. Reference Asioli, Trincardi, Lowe, Ariztegui, Langone and Oldfield2001; Pellegrini et al. Reference Pellegrini, Maselli, Gamberi, Asioli, Bohacs, Drexler and Trincardi2017, Reference Pellegrini, Asioli, Bohacs, Drexler, Feldman, Sweet and Maselli2018).
Scarponi et al. (Reference Scarponi, Nawrot, Azzarone, Pellegrini, Gamberi, Trincardi and Kowalewski2022) investigated how molluscan assemblages from nearshore settings responded to these major environmental turnovers. They found strong similarity in the faunal composition of assemblages from the last interglacial and the preindustrial Holocene, both of which differed significantly from their glacial counterparts. This pattern points to a resilient response of Adriatic benthic communities to major climatic fluctuations. Community composition was most likely governed by the interplay of species’ environmental tolerances and local environmental conditions, following the environmental assembly model (Jackson and Blois Reference Jackson and Blois2015), which led to the formation of non-analogous communities during the glacial but the reappearance of communities from the last interglacial when previous climatic conditions returned during the Holocene. While assemblages from the preindustrial Holocene were very similar to those of the last interglacial, they differ significantly from those of modern nearshore settings (Kowalewski et al. Reference Kowalewski, Wittmer, Dexter, Amorosi and Scarponi2015; Scarponi et al. Reference Scarponi, Rojas, Nawrot, Cheli, Kowalewski, Nawrot, Dominici, Tomašových and Zuschin2023). The response of molluscan communities to Quaternary environmental change thus informs us about their natural range of variability, a baseline against which the effects of recent multifaceted anthropogenic changes can be assessed.
The formation of non-analogous communities in response to Quaternary climate change has been described by numerous other studies (e.g., Jackson and Overpeck Reference Jackson and Overpeck2000; Tager et al. Reference Tager, Webster, Potts, Renema, Braga and Pandolfi2010; Pardi and Graham Reference Pardi and Graham2019) and anthropogenic impacts are equally leading to the emergence of novel ecosystems (Hobbs et al. Reference Hobbs, Higgs and Hall2013; Morse et al. Reference Morse, Pellissier, Cianciola, Brereton, Sullivan, Shonka, Wheeler and McDowell2014; Barnosky et al. Reference Barnosky, Hadly, Gonzalez, Head, Polly, Lawing and Eronen2017). However, it is unclear to which degree such turnovers in species composition would also result in major changes in ecosystem functioning. On the one hand, species that disappear could be replaced by ecologically similar ones, and even in the case of high biodiversity loss, functional redundancy could buffer ecosystem functioning (Fonseca and Ganade Reference Fonseca and Ganade2001; Nunes et al. Reference Nunes, Barlow, França, Berenguer, Solar, Louzada and Leitão2021; Yeeles et al. Reference Yeeles, Lach, Hobbs and Didham2025). On the other hand, some ecological roles might be performed by only a few species and therefore be vulnerable to extinction when major environmental perturbations occur (Mouillot et al. Reference Mouillot, Villéger, Parravicini, Kulbicki, Arias-González, Bender and Chabanet2014; Teichert et al. Reference Teichert, Lepage, Sagouis, Borja, Chust, Ferreira, Pasquaud, Schinegger, Segurado and Argillier2017; Carlot et al. Reference Carlot, Galobart, Gómez-Gras, Santamaría, Golo, Sini and Cebrian2025). As the retention and augmentation of ecosystem services are crucial to biodiversity conservation, there have been increasing calls for taxon-free approaches, such as studies on functional traits (Cadotte et al. Reference Cadotte, Carscadden and Mirotchnick2011; Gagic et al. Reference Gagic, Bartomeus, Jonsson, Taylor, Winqvist, Fischer and Slade2015; Barnosky et al. Reference Barnosky, Hadly, Gonzalez, Head, Polly, Lawing and Eronen2017; Guy-Haim and Bouchet Reference Guy-Haim and Bouchet2025). Such approaches provide the additional advantage of higher predictive power across regions, where taxonomic identity differs but ecosystem functioning might be very similar (Charvet et al. Reference Charvet, Statzner, Usseglio-Polatera and Dumont2000; Statzner and Bêche Reference Statzner and Bêche2010; Villéger et al. Reference Villéger, Grenouillet and Brosse2013). The same argument can be made for paleontological studies that compare assemblages across time periods where extinction, migration, and speciation have led to taxonomically distinct communities (e.g., Villéger et al. Reference Villéger, Novack-Gottshall and Mouillot2011; Dineen et al. Reference Dineen, Fraiser and Sheehan2014; Edelman-Furstenberg and Kidwell Reference Edelman-Furstenberg and Kidwell2015; Knope et al. Reference Knope, Heim, Frishkoff and Payne2015).
This study aims to assess the effects of Late Quaternary climate change on the functional diversity of nearshore molluscan assemblages of the Adriatic Sea by applying biological traits analysis (BTA). This multivariate approach uses a multitude of life-history, morphological and behavioral traits that are relevant to a species’ role in the ecosystem (Bremner et al. Reference Bremner, Rogers and Frid2006). It is an established approach in ecological studies on functional diversity (e.g., Tillin et al. Reference Tillin, Hiddink, Jennings and Kaiser2006; Paganelli et al. Reference Paganelli, Marchini and Occhipinti-Ambrogi2012; Bolam and Eggleton Reference Bolam and Eggleton2014; Floyd et al. Reference Floyd, Mizuyama, Obuchi, Sommer, Miller, Kawamura, Kise, Reimer and Beger2020) and is increasingly applied in historical ecology, conservation paleobiology, and paleoecology as well (e.g., Steger et al. Reference Steger, Bošnjak, Belmaker, Galil, Zuschin and Albano2022, Reference Steger, Bogi, Lubinevsky, Galil, Zuschin and Albano2024; Grossmann et al. Reference Grossmann, Nielsen, Rivadeneira and Valdivia2023; Tyler and Kowalewski Reference Tyler and Kowalewski2025). We focus on the functional diversity of benthic mollusks, which are a crucial component of macrobenthic ecosystems due to their high abundance and diversity (Webb et al. Reference Webb, Aleffi, Amouroux, Bachelet, Degraer, Dounas and Fleischer2009; Khan et al. Reference Khan, Manokaran and Lyla2017; Shou et al. Reference Shou, Liao, Tang, Chen, Jiang, Gao and Chen2018; Tyler and Kowalewski Reference Tyler and Kowalewski2018). Similar to other benthic organisms, they are often used as ecosystem indicators, due to their limited mobility and high sensitivity to environmental change (Nerlović et al. Reference Nerlović, Doğan and Hrs-Brenko2011; Moraitis et al. Reference Moraitis, Tsikopoulou, Geropoulos, Dimitriou, Papageorgiou, Giannoulaki, Valavanis and Karakassis2018; Fitzgerald et al. Reference Fitzgerald, Ryan, Scarponi and Huntley2024). It has been demonstrated that mollusks can serve as surrogates for the entire benthic community and that their death assemblages have high ecological fidelity due to their good preservation potential (Smith Reference Smith2005; Kidwell and Tomašových Reference Kidwell and Tomašových2013; Tyler and Kowalewski Reference Tyler and Kowalewski2017, Reference Tyler and Kowalewski2025). Our dataset comprises extant species with well-known ecology. Therefore, information on a diverse set of functional traits is available for these species (see Supplementary Data 1). We selected seven traits: feeding guild, motility, attachment type, substrate relation, maximum adult body size, type of larval development, and vegetation association. These traits (see Supplementary Table S1) affect crucial ecosystem processes, such as bioturbation, nutrient recycling, and secondary biomass production, and are responsive to environmental perturbations. Our study helps to address whether major species turnovers caused by long-term climate change are accompanied by changes in functional composition and diversity, with crucial implications for ecosystem functioning.
Materials and Methods
Samples
The samples studied herein come from the last interglacial (hereafter referred to as LIG), the last glacial (LG), and the current interglacial (Holocene, CIG). The interglacial samples are derived from sediment cores taken in the present-day Po coastal plain where LIG and CIG marine sediments are preserved under more recent alluvial and delta plain deposits (Fig. 1). Conversely, LG samples come from sediment cores taken more than 250 km southeast of the interglacial localities in the central and southern Adriatic, at the edge of the Mid-Adriatic Deep, where the shoreline was located during the late LG (Pellegrini et al. Reference Pellegrini, Asioli, Bohacs, Drexler, Feldman, Sweet and Maselli2018; Fig. 1). As with most fossil assemblages (Kidwell and Tomašových Reference Kidwell and Tomašových2013), our samples may represent different degrees of time averaging. All samples used in this study are derived from aggradational to progradational lower shoreface to foreshore deposits with varying degrees of fluvial influence, hereafter referred to as nearshore, to ensure consistency in environmental context, sedimentation rates, and taphonomic regime across time periods. By restricting the datasets to this narrow range of stratigraphic and depositional settings, we aimed to minimize variation not only in environmental context and related physical processes but also in the temporal resolution represented by individual samples. The analyzed samples are time-averaged on the scale of a few centuries (Scarponi et al. Reference Scarponi, Kaufman, Amorosi and Kowalewski2013). Therefore, our study is not designed to estimate rates of community change or transient ecological dynamics during short-term (multi-decadal to centennial) warming or cooling events but to compare functional composition and diversity across distinct climatic states. LIG samples have been limited to the period from ~125 to 110 ka cal BP (within MIS 5), LG samples from ~18 to 12 ka cal BP (within MIS 2), and CIG samples from ~7 ka cal BP to 1750 CE (within MIS 1).
The present-day Adriatic Sea and its coastline during the last glacial maximum (LGM; based on Maselli et al. Reference Maselli, Trincardi, Asioli, Ceregato, Rizzetto and Taviani2014). The sites of cores from which samples from the last interglacial (LIG; n = 21) and the current interglacial (CIG; n = 170) are derived are shown in yellow; those from which last glacial (LG) samples (n = 32) are derived are shown in blue.
Figure 1. Long description
A map of the Adriatic Sea region with a coordinate grid from 12.000 to 20.000 East and 40.000 to 46.000 North.
In the top right, an inset map shows the Mediterranean Sea with a black box highlighting the Adriatic study area.
The main map displays the present-day Adriatic Sea in white and the landmass in dark gray. A light gray shaded area surrounding the white sea represents land exposed during the L G M, bounded by a dashed line indicating the L G M coastline.
Data points are distributed as follows:
* Interglacial sites: A cluster of yellow dots is located in the Northwest, along the Italian coast near 44.500 North.
* Glacial sites: Blue dots are located further Southeast. Three dots are clustered in the central Adriatic near 43.000 North, and one dot is positioned further South along the eastern coast near 41.500 North.
A legend in the bottom left identifies yellow circles as interglacial sites, blue circles as glacial sites, the dashed line as the coastline during L G M, and light gray shading as area exposed during the L G M. A scale bar in the bottom right indicates 0, 75, and 150 kilometers. A red North arrow and the text E P S G 4326 are in the top left.
Mollusk assemblages were obtained from sediment samples processed following the protocol described by Wittmer et al. (Reference Wittmer, Dexter, Scarponi, Amorosi and Kowalewski2014). Samples were dried, disaggregated using a 4% hydrogen peroxide solution, and wet-sieved through a 1-mm mesh. All identifiable molluscan remains were picked from the sieve residue and identified to species level whenever possible. Complete specimens and unique fragments (with umbo in bivalves, with apex in gastropods) were counted.
Trait Dataset
In addition to the information on species abundances and sample properties, we compiled data from the literature and online databases on the feeding guild, motility, attachment type, substrate relation, maximum adult body size, type of larval development, and vegetation association of each species (a reference list is provided in Supplementary Data 1). The traits were scored using a fuzzy coding scheme (Chevenet et al. Reference Chevenet, Doléadec and Chessel1994), which allows species to be assigned to multiple categories or modalities (sensu Chevenet et al. Reference Chevenet, Doléadec and Chessel1994) of a trait. For example, a suspension feeder that opportunistically shifts to deposit feeding can be assigned to both modalities this way. We chose 30 modalities for all traits combined (Supplementary Table S1). Depending on the affinity of a species to a trait modality, a score from 0 to 3 was assigned, with 0 indicating no affinity; 1, low affinity; 2, high affinity; and 3, exclusive affinity. If a species had the same level of affinity to multiple modalities, it was given a score of 2 for each by default. To weigh all traits equally, we converted species scores for trait modalities to proportions adding up to unity. In the rare case that information on a trait was unobtainable for a species, each modality was assigned a score of 0. This way, the species is not taken into account for this trait during the analyses.
Data Processing and Culling
Before the analysis, all specimens that could not be identified to species level were removed from the dataset to ensure that traits were scored and treated at the same taxonomic level. Freshwater and terrestrial species, occasionally found in our samples, were also omitted from the analysis because they are allochthonous to the studied depositional settings. Counts of isolated bivalve valves were halved to approximate the number of individuals. All samples with fewer than 25 specimens and monospecific assemblages were removed. Finally, trait modalities that were not represented by any species were excluded. The final dataset comprises 35,986 specimens from 98 samples and 107 species, representing 28 modalities across all traits (interval-specific numbers of specimens, samples, and modalities in Supplementary Table S2). Data processing and all subsequent analyses were performed in R statistical environment v. 4.5.2 (R Core Team 2025).
Analysis of Functional Composition
To test for differences in functional composition between the assemblages, we applied fuzzy correspondence analysis (FCA), an ordination method representing an extension of multiple correspondence analysis, which is suitable for the analysis of discrete variables such as our trait modalities (Chevenet et al. Reference Chevenet, Doléadec and Chessel1994). To perform the FCA, we created a matrix of relative abundances of trait modalities per sample by multiplying the relative species abundances and trait modality scores for each species (relative abundances and standardized number of species of each trait modality per time interval are shown in Supplementary Figs. S1 and S2). When ordination methods are applied to explore differences in taxonomic composition, the abundance matrix is often transformed or standardized (e.g., log-transformation, Wisconsin double-relativization) to downweigh the influence of very abundant species. For the analysis of functional composition, we used relative abundances without further transformation because differences in abundance between samples are already much reduced when different species are aggregated into trait modalities and downweighing the effect of very abundant species, which can be the main drivers of ecosystem functioning, is counterproductive (Winfree et al. Reference Winfree, Fox, Williams, Reilly and Cariveau2015; Gaston et al. Reference Gaston, Cox, Canavelli, García, Hughes, Maas, Martínez, Ogada and Inger2018). Nevertheless, to test the effects of transformations on the results, we repeated the FCA with log and square-root transformation of the species abundance matrix. The results of these analyses were comparable to the FCA without prior transformation (Supplementary Fig. S3, Supplementary Table S3). FCA can depict both samples (by using the weighted average of the trait modalities represented in them) and trait modalities (by using the weighted average of the samples in which they occur) in the same ordination space. In addition to visualizing differences in functional composition, FCA depicts the separation among modalities of a trait by providing correlation ratios, which represent proportions of the total variance explained by each axis (Chevenet et al. Reference Chevenet, Doléadec and Chessel1994). Therefore, correlation ratios can be used to discern which traits differ the most between samples.
A major difference between the LG and the interglacial assemblages is the high dominance of the bivalve Lentidium mediterraneum in the latter (>85%; Scarponi et al. Reference Scarponi, Nawrot, Azzarone, Pellegrini, Gamberi, Trincardi and Kowalewski2022). Due to its very high abundance, this species might strongly influence differences in functional composition. To evaluate the impact of this single, dominant species on BTA, we repeated the FCA with a dataset without L. mediterraneum (Supplementary Table S2), using the same culling criteria for samples (≥25 specimens and ≥2 species per sample). FCA was performed using the ade4 package v. 1.7-23 (Dray and Dufour Reference Dray and Dufour2007).
To corroborate the conclusions drawn from the FCA results, we also calculated pairwise distances between samples based on Bray-Curtis (BC) dissimilarity for all three pairwise comparisons between time intervals (LIG–LG, LG–CIG, LIG–CIG) using the same sample-trait matrix as for the FCA. This resulted in distributions of BC distances for each comparison, indicating whether mean BC distance is greater for glacial–interglacial comparisons than for the comparison between the interglacials. As sample sizes differ strongly across time intervals (see Supplementary Table S2), differences in functional diversity might be affected by unbalanced sampling. To account for this, we generated null distributions by pooling all specimens, randomly reallocating them into the original sampling structure of each interval (i.e., the same number of samples and per-sample specimen counts as observed), and then recalculating pairwise BC distances and their mean for each comparison. This process was repeated 999 times, and the resulting distributions of mean BC distances were compared with the observed values. This null model allowed us to statistically evaluate whether observed differences in BC distance are the result of different functional composition or unbalanced sampling across the time intervals. BC distances were calculated using the vegan package v. 2.7-1 (Oksanen et al. Reference Oksanen, Simpson, Blanchet, Kindt, Legendre, Minchin and O’Hara2025).
Analysis of Functional Diversity
Functional diversity represents the diversity of trait profiles, thus indicating the diversity of species’ niches or functions in an ecosystem (Cadotte et al. Reference Cadotte, Carscadden and Mirotchnick2011). To test for differences in functional diversity, we assessed the distribution of the assemblages from each time interval in functional space, a multidimensional coordinate system in which species are arranged according to their trait (dis)similarity (Villéger et al. Reference Villéger, Mason and Mouillot2008; Laliberté and Legendre Reference Laliberté and Legendre2010). To create functional space, first, trait-based distances between species were calculated using the Gower measure, which can be applied to fuzzy-coded data (Gower Reference Gower1971; De Bello et al. Reference De Bello, Botta‐Dukát, Lepš and Fibich2021). Then, based on these distances, a principal coordinate analysis was computed to create a multidimensional functional space (Villéger et al. Reference Villéger, Mason and Mouillot2008; Laliberté and Legendre Reference Laliberté and Legendre2010). The selection of the dimensionality of the functional spaces is a crucial step, as it affects how closely the original distances are preserved (Maire et al. Reference Maire, Grenouillet, Brosse and Villéger2015; Magneville et al. Reference Magneville, Loiseau, Albouy, Casajus, Claverie, Escalas, Leprieur, Maire, Mouillot and Villéger2022). We chose a four-dimensional space, which had a mean squared deviation of 0.005 between Gower and Euclidean distances (Supplementary Table S5). We consider it a good compromise between the quality of the functional space and computational power (further details in Supplementary Material).
Based on functional space, we calculated four indices of assemblage-level functional diversity for each time interval (after pooling all samples of a given age together): functional dispersion, functional evenness, functional richness, and functional divergence (Villéger et al. Reference Villéger, Mason and Mouillot2008; Laliberté and Legendre Reference Laliberté and Legendre2010). Functional dispersion is the mean distance of species to the centroid of an assemblage, with both the centroid and distances weighted by abundances. Therefore, it provides information on the extent to which trait profiles and the functions provided by organisms vary on average (Laliberté and Legendre Reference Laliberté and Legendre2010). Functional evenness measures how evenly species and their abundances are distributed in functional space, based on the minimum spanning tree connecting all species in an assemblage. Therefore, it indicates whether species performing different functions have comparable abundances or not. Functional richness is the total volume occupied by species in functional space and therefore represents the range of functional profiles in a given assemblage. It indicates what range of niches is available under given environmental conditions. Functional divergence defines how species’ abundances are distributed in functional space relative to a centroid computed from species’ positions (not abundance-weighted, unlike for functional dispersion). High values are attained when dominant species are located far from the centroid, indicating that the most commonly performed functions are very distinct from the functions provided by rarer species (Villéger et al. Reference Villéger, Mason and Mouillot2008). Functional evenness and divergence are by definition constrained between 0 and 1, whereas functional dispersion and richness have no theoretical upper limit (Villéger et al. Reference Villéger, Mason and Mouillot2008; Laliberté and Legendre Reference Laliberté and Legendre2010). Therefore, to facilitate comparison across the time intervals, we divided the latter two indices by their maxima across all assemblages to obtain values between 0 and 1 (Magneville et al. Reference Magneville, Loiseau, Albouy, Casajus, Claverie, Escalas, Leprieur, Maire, Mouillot and Villéger2022). To account for unbalanced sampling, the dataset was resampled in the same way as for the analysis of BC dissimilarity. Thus, the original sampling structure of each interval was kept, specimens were redistributed freely across time intervals, and then the indices were recalculated for each assemblage. This process was repeated 999 times. This null model provided context for the observed values. Furthermore, we calculated functional richness for all three assemblages standardized to the sample size of the smallest assemblage, that of the LG (n = 1955), as this index is strongly related to species richness (Villéger et al. Reference Villéger, Mason and Mouillot2008), and consequently sample size (further details in the Supplementary Material).
Finally, we analyzed differences between the assemblages by calculating pairwise functional beta diversity between the three time intervals through a Jaccard-like index (Villéger et al. Reference Villéger, Grenouillet and Brosse2013). This index, constrained between 0 and 1, can be decomposed into turnover and nestedness-resultant components. If two assemblages show little overlap in functional space, the turnover component will attain high values. Conversely, the nestedness-resultant component will be high if one assemblage fills a small portion of the volume occupied by the other assemblage (Villéger et al. Reference Villéger, Grenouillet and Brosse2013). Differences in sample size also constitute an obstacle to this analysis. Therefore, we applied the same resampling procedure as for BC dissimilarity and indices of functional diversity, again using 999 permutations, to obtain values of beta diversity and its components, assuming that differences between assemblages are due solely to sample size. The resampling distributions were compared with observed values. The creation of functional space and the calculation of functional diversity indices and beta diversity were performed using the mFD package v. 1.0.7 (Magneville et al. Reference Magneville, Loiseau, Albouy, Casajus, Claverie, Escalas, Leprieur, Maire, Mouillot and Villéger2022).
Results
FCA
In our interpretation of the FCA, we only considered the first three axes, which cumulatively account for 85.1% of the variance in the data, explaining 49.4 %, 25.9%, and 9.8%, respectively. Glacial and interglacial samples are separated along axis 1, albeit with some overlap (Fig. 2). Conversely, samples from the LIG and CIG overlap strongly. The LG centroid is clearly distinct from those of the LIG and CIG. On the other hand, there is a greater similarity in the position of the samples from the three intervals along axis 2 and axis 3 (Supplementary Fig. S4). However, several LG samples plot outside the range of the interglacial samples and several LIG and CIG samples plot outside the range of the glacial samples along these axes as well. Overall, the distribution of samples in FCA ordination space clearly indicates differences in functional composition between the glacial and the interglacials.
Fuzzy correspondence analysis (FCA) based on a matrix of relative abundance of trait modalities. The first two axes shown here explain 75.3% of the variance in the data. Circles represent samples (last interglacial [LIG] in green, last glacial [LG] in blue, current interglacial [CIG] in red) and their size is proportional to the decadal logarithm of the sample size. The minimum convex hull for each time interval is shown in the same color as the samples, and centroids are shown as crosses. In addition to samples, the modalities of attachment type (dark gray, highest correlation ratio with axis 1), feeding guild (light gray, second-highest correlation ratio with axis 1), and maximum adult size (white, highest correlation ratio with axis 2) are shown as rectangles. Remaining traits are shown in Supplementary Fig. S5. Abbreviations: by, byssally attached; ce, cemented; cs, (chemo)symbiotic feeder; df, deposit and detritus feeder; hg, herbivore/grazer; la, large; me, medium; pa, parasitic feeder; pr, predator; sc, scavenger; sm, small; sf, suspension feeder; un, unattached; vs, very small.
Figure 2. Long description
An attribute space diagram with Axis 1 on the horizontal ranging from negative 2.5 to 0.5 (49.4 percent variance) and Axis 2 on the vertical ranging from negative 1.5 to 2.0 (25.9 percent variance).
Data Clusters:
* A large blue convex hull representing L G (last glacial) samples dominates the left and center, extending from Axis 1 value negative 2.2 to 0.0. Its centroid cross is located near Axis 1 negative 0.6.
* A red convex hull representing C I G (current interglacial) samples is concentrated on the right, overlapping with a green hull representing L I G (last interglacial) samples near the origin.
* Sample circles vary in size based on the decadal logarithm of sample size.
Trait Modality Labels (Rectangles):
* Far Left: h g (herbivore/grazer) at the top left peak of the blue hull; v s (very small) in the center-left.
* Center: c s (chemosymbiotic), d f (deposit/detritus feeder), u n (unattached), and p a (parasitic) cluster around the L G centroid.
* Bottom Center: m e (medium), l a (large), p r (predator), and s c (scavenger) are positioned below the main clusters.
* Far Right: s m (small), b y (byssally attached), s f (suspension feeder), and c e (cemented) are clustered at the positive end of Axis 1, associated with the C I G and L I G hulls.
To assess how traits vary along the FCA axes, we plotted the modalities alongside samples (Fig. 2, Supplementary Fig. S5) and defined the correlation ratios (Table 1). Attachment type has the highest correlation ratio of 0.381 with axis 1, indicating notable differences in attachment type between the time intervals. Indeed, more than 70% of the species in most LG samples have unattached live habit, while byssally attached species dominate in LIG and CIG samples (Supplementary Fig. S1C). Feeding guild has the second-highest correlation ratio with axis 1 at 0.257. Suspension feeders predominate in most samples across all time intervals, but their relative abundance tends to be higher in interglacial samples, while herbivores and deposit feeders often attain higher abundances in LG samples (Supplementary Fig. S1A). Interestingly, the herbivore-grazer modality, as well as plant and algae association, plot closely to the outlying LG samples (Fig. 2, Supplementary Fig. S5B), suggesting that they represent vegetated habitats, whereas there is no such indication for any interglacial samples. Maximum adult size shows a comparatively low correlation ratio with axis 1 at 0.161, but the highest with axis 2 (0.289). The two dominant size classes are small (4–16 mm) and medium (16–64 mm). While medium-sized species tend to be more common in LG than in LIG and CIG samples, there is very strong variation in the proportion of these size classes within time intervals (Supplementary Fig. S1E).
Correlation ratios of each trait with the first three axes of the fuzzy correspondence analysis (FCA). These ratios represent proportions of the total variance explained by each axis to depict the separation among modalities of a trait (Chevenet et al. Reference Chevenet, Doléadec and Chessel1994)
Table 1. Long description
The table consists of seven columns representing biological traits and three rows representing F C A axes.
* Axis 1 values: Attachment type 0.381, Feeding guild 0.257, Motility 0.207, Vegetation association 0.201, Substrate relation 0.200, Maximum adult size 0.161, and Larval development 0.064.
* Axis 2 values: Attachment type 0.144, Feeding guild 0.131, Motility 0.055, Vegetation association 0.084, Substrate relation 0.066, Maximum adult size 0.289, and Larval development 0.003.
* Axis 3 values: Attachment type 0.005, Feeding guild 0.005, Motility 0.060, Vegetation association less than 0.001, Substrate relation 0.216, Maximum adult size 0.005, and Larval development 0.002.
The exclusion of Lentidium mediterraneum strongly changes the results (see Supplementary Material). There is no clear separation of time intervals along any of the three FCA axes (Supplementary Fig. S6). Furthermore, feeding type and vegetation association have the highest correlation ratios with axis 1, whereas attachment type has comparatively low correlation ratios with axis 1 and 2 (Supplementary Table S4). This confirms that the functional separation of the LG and the interglacials is strongly linked to the interglacial dominance of L. mediterraneum.
BC Dissimilarity
The analysis of BC distances corroborates the conclusion that glacial and interglacial assemblages differ in functional composition (Fig. 3). The mean BC distance is higher for the comparisons between the LG and the LIG (x̄ = 0.288) or the CIG (x̄ = 0.304) than for the comparison between the interglacials (x̄ = 0.154). While no BC distance between any LIG and CIG sample is higher than 0.4, the distance of LG samples to interglacial samples can be up to almost 0.7. These results cannot be attributed to unbalanced sampling, because observed means differ significantly from those of the null distributions generated by randomly assigning specimens to the temporal intervals (Fig. 3). The observed mean BC distance is significantly higher than the mean of the resampling distribution for all comparisons, but the difference is about twice as high for the glacial–interglacial comparisons.
Bray-Curtis (BC) distances for pairwise comparison of samples from (A) the last interglacial (LIG) and last glacial (LG), (B) the LG and current interglacial (CIG), and (C) the LIG and CIG. Distributions are shown in light gray and their means as solid gray lines. Histograms shown in transparent dark gray are resampling distributions of mean BC distances created by randomly reassigning the temporal affinity of specimens (n per = 999). The means of these resampling distributions are shown as dashed gray lines. The differences between these means and the observed means (solid gray lines) ∆x̄ are given for each comparison of time intervals. The p-values estimating the significance of the observed means were calculated as the proportion of permutations that deviated more from the mean of the resampling distribution than the observed mean (p-values were significant for all three comparisons).
Figure 3. Long description
A multi-panel histogram with three vertically stacked panels labeled A, B, and C. The x-axis represents Bray-Curtis dissimilarity ranging from 0.0 to 1.0. The y-axis represents the number of pairwise comparisons from 0 to 400.
* Panel A (L I G-L G): A wide, light gray distribution of observed values is spread across the x-axis with a solid gray mean line at approximately 0.29. A very narrow, dark gray resampling distribution is clustered near 0.03 with a dashed mean line. Text indicates p equals 0.001 and Delta x-bar equals 0.257.
* Panel B (L G-C I G): The light gray observed distribution is shifted further right with a solid mean line at approximately 0.30. The narrow dark gray resampling distribution remains near 0.03. Text indicates p equals 0.001 and Delta x-bar equals 0.271.
* Panel C (L I G-C I G): The light gray observed distribution is more concentrated toward the left than the previous panels, with a solid mean line at approximately 0.15. The narrow dark gray resampling distribution is again near 0.03. Text indicates p equals 0.001 and Delta x-bar equals 0.124.
Patterns in Functional Diversity
The results of functional diversity analysis (Fig. 4, Supplementary Table S6) indicate that, in addition to differences in composition, there are also differences in diversity between the time intervals. Functional dispersion is more than three times higher in the LG assemblage than in the interglacial assemblages. This pattern becomes even more significant when we compare the observed values to the null distributions created based on the assumption that the three time periods represent the same assemblage and differences are attributable to sample size alone. According to this null model, the functional dispersion of the LG assemblage should be comparable to that of the interglacials (Fig. 4A). Functional evenness is more similar between time intervals, but it is also highest in the LG assemblage, while the observed values for the interglacials are lower than expected by the null model (Fig. 4B). Functional richness shows an apparent increase through time from the LIG to the CIG. However, functional richness of the CIG assemblage is lower than expected by the null model. While this is also true for the LIG, the observed value for the LG lies within the 95th percentile of the resampling distribution (Fig. 4C). When all assemblages are standardized to the sample size of the LG, functional richness is much higher for the LG than for the interglacials (Supplementary Fig. S8). Functional divergence is very high (>0.90) in all assemblages (Fig. 4D). This indicates that the most abundant species have very distinct trait profiles (see also Supplementary Fig. S7).
Functional diversity measured as (A) functional dispersion, (B) functional evenness, (C) functional richness, and (D) functional divergence. Observed values are shown as diamonds and are given in Supplementary Table S6. The last interglacial (LIG) is shown in green, the last glacial (LG) in blue, and the current interglacial (CIG) in red. The box plots represent resampling distributions of the indices created by randomly reassigning specimens to each time interval, keeping the number of samples and specimens of each time interval constant (n per = 999). Shaded areas behind the box plots represent the 95th percentile intervals of these distributions. Box plots and shaded areas are all shown in gray because they represent a null model that all samples come from the same assemblage.
Figure 4. Long description
A multi-panel graph with a shared Y-axis labeled Functional diversity ranging from 0.0 to 1.0. Each panel contains three categories on the X-axis: L I G (green diamond), L G (blue diamond), and C I G (red diamond), overlaid on gray box plots with shaded 95th percentile intervals.
* Panel A. Functional Dispersion: L I G and C I G diamonds are within or near the null model range at approximately 0.2, while the L G diamond is a significant outlier at 0.7.
* Panel B. Functional Evenness: All three diamonds (L I G, L G, C I G) fall within or very close to the null model box plots, centered between 0.4 and 0.5.
* Panel C. Functional Richness: L I G (0.55) and C I G (0.85) diamonds are outliers below their respective null model distributions, while the L G diamond (0.72) is within the upper range of its null model.
* Panel D. Functional Divergence: All three diamonds are positioned at the top of the scale near 1.0. L I G and L G are slightly below the null model, while C I G is at the very top, matching the null model median.
Functional beta diversity is similar for all comparisons, but they differ in the relative contributions of the components of beta diversity (Fig. 5, Supplementary Table S7). The differences between the LG and both the LIG and CIG are mainly driven by the turnover component, demonstrating that glacial and interglacial assemblages occupy distinct portions of functional space. On the other hand, beta diversity between the LIG and CIG is mainly driven by the nestedness-resultant component. Based on our null model, dominance of the nestedness-resultant component is expected in all comparisons, which is unsurprising given the much larger sample size of the CIG assemblage.
Functional beta diversity between the (A) last interglacial (LIG) and last glacial (LG), (B) the LG and current interglacial (CIG), and (C) the LIG and CIG. Observed values of beta diversity (calculated as Jaccard-like dissimilarity), its turnover, and nestedness-resultant component are shown as black diamonds and given in Supplementary Table S7. The box plots represent resampling distributions of the beta diversity and its components created by randomly assigning specimens to each time interval, keeping the number of samples and specimens of each time interval constant (n per = 999). Shaded areas behind the box plots represent the 95th percentile intervals of these distributions.
Figure 5. Long description
A multi-panel graph with a shared vertical Y axis labeled Functional beta diversity ranging from 0.0 to 0.7. The horizontal X axis for each panel contains three categories: J D, T O, and N R. Each category features a dark gray box plot with a median line, whiskers, a light gray shaded background representing the 95th percentile, and a black diamond representing observed values.
* Panel A, L I G-L G: J D has a median near 0.3 and an observed diamond at 0.4. T O has a low median near 0.05 and an observed diamond at 0.23. N R has a median near 0.23 and an observed diamond at 0.16.
* Panel B, L G-C I G: J D has a median near 0.4 and an observed diamond at 0.32. T O has a median at 0.0 and an observed diamond at 0.22. N R has a median near 0.39 and an observed diamond at 0.1.
* Panel C, L I G-C I G: J D has a median near 0.21 and an observed diamond at 0.35. T O has a median near 0.01 and an observed diamond at 0.02. N R has a median near 0.2 and an observed diamond at 0.33.
In all panels, the observed values (diamonds) often fall outside the interquartile range of the resampling distributions, particularly for the T O component.
Discussion
The results of this study contribute to understanding the effects of long-term climate change on modern marine ecosystems, as only a few studies have focused on changes in marine functional diversity during the Late Quaternary (e.g., Ávila et al. Reference Ávila, Arruda, Silva, Hipólito, Quartau, Ramalho and Madeira2025). We found that the taxonomic turnover in nearshore molluscan assemblages of the western Adriatic during the Late Quaternary led to changes in functional composition and diversity. Thus, species that replaced each other or changed in abundance at glacial–interglacial transitions were at least partly functionally distinct. One potential caveat of our study is that the ecological fidelity of fossil assemblages can be compromised by taphonomic biases that lead to differential preservation of taxa and, potentially, functional groups (Barrena Reference Barrena1993; Albano and Sabelli Reference Albano and Sabelli2011; García-Ramos et al. Reference García-Ramos, Albano, Harzhauser, Piller and Zuschin2016). However, the results of a comprehensive study of the functional diversity of all major taxonomic groups of marine macrobenthos, including mollusks, in Onslow Bay, North Carolina (USA), suggest that fossil assemblages archive functional diversity patterns with high fidelity relative to the original communities, despite preservation biases (Tyler and Kowalewski Reference Tyler and Kowalewski2025). In a literature review on the macrobenthic assemblages of other regions worldwide, Tyler and Kowalewski (Reference Tyler and Kowalewski2025) found that most have a similar set of higher taxa with comparable preservation potential to Onslow Bay. Therefore, their results are likely applicable to many benthic assemblages. Furthermore, the assemblages studied herein are dominated by small- to medium-sized, thin-shelled aragonitic bivalves with comparable durability, and thus taphonomic filtering is unlikely to have played an important or differential role in modulating their composition.
Differences in the scale of time averaging between samples could represent another confounding factor due to the positive relationship between time averaging and species richness (Kidwell and Tomašových Reference Kidwell and Tomašových2013; Tomašových et al. Reference Tomašových, Nawrot, Kowalewski, Scarponi and Zuschin2024), which in turn could affect patterns in functional richness (Villéger et al. Reference Villéger, Mason and Mouillot2008). In our case, several aspects of the study design reduce this concern: samples were restricted to fluvially influenced nearshore depositional settings, limiting variation in time averaging, and allochthonous terrestrial and freshwater species were excluded. We therefore consider it unlikely that differences in time averaging represent the main driver of the functional diversity patterns observed at the scale of glacial and interglacial time intervals.
The high fidelity of the functional structure of fossil assemblages is highly relevant, because experiments (Gagic et al. Reference Gagic, Bartomeus, Jonsson, Taylor, Winqvist, Fischer and Slade2015; Lefcheck and Duffy Reference Lefcheck and Duffy2015) and ecological surveys (Bolam and Eggleton Reference Bolam and Eggleton2014) have demonstrated that trait-based approaches are a useful tool to link community composition and ecosystem functioning. Therefore, it is crucial to understand what might have led to changes in the functional diversity of nearshore communities in the past to predict how ecosystems could respond to ongoing and future environmental change in the Adriatic.
Drivers of Higher Functional Diversity during the LG
The LG assemblage differs in functional composition from those of the interglacials (Figs. 2, 3). The decomposition of beta diversity between time intervals (Fig. 5) clearly indicates that this is due to a turnover in functional profiles, as the corresponding component is higher. Our null model, based on the assumption that the assemblage composition was the same during the three time intervals, suggests that the nestedness-resultant component should be higher if differences between assemblages were caused by uneven sampling effort. As this is not observed despite large differences in sample size, our results are robust. In addition to its distinct composition, the LG assemblage also shows higher functional diversity than the interglacial assemblages, as indicated by multiple indices. Also this difference cannot be explained by uneven sampling effort, because observed values differ strongly from the expectations derived from our null model. Functional dispersion and evenness are higher in the LG assemblage, indicating that functional trait profiles were more variable on average and that individuals were more evenly distributed among them. Functional richness is highest for the CIG assemblage, but this is due to its much larger sample size. If all assemblages are standardized to the same sample size, functional richness is significantly higher for the LG assemblage (Supplementary Fig. S8). Therefore, glacial communities also had a wider range of functional profiles. Functional divergence is high in all assemblages, indicating that they are dominated by species with distinct trait profiles (see distribution of species in functional space in Supplementary Fig. S7).
The observed differences between the glacial and interglacials are not unexpected given the patterns in taxonomic composition and diversity described by Scarponi et al. (Reference Scarponi, Nawrot, Azzarone, Pellegrini, Gamberi, Trincardi and Kowalewski2022). Assemblages from the LIG and CIG are characterized by very high dominance of Lentidium mediterraneum, which makes up more than 85% of the specimens from either interglacial. This species is also very common in shallow nearshore environments of the present-day Po Delta area, where it can reach mean population densities of 130,000 individuals/m3 and relative abundances of 95% in macrobenthic communities (Ambrogi and Bedulli Reference Ambrogi and Bedulli1981). For comparison, the most abundant species of the LG assemblage, the bivalve Spisula subtruncata, makes up 44% of individuals, while the relative abundance of L. mediterraneum, still ranking third, is one order of magnitude lower (~8%) than in the interglacial assemblages. The large disparity in the relative abundance of L. mediterraneum is reflected in the functional composition: the trait that best characterizes compositional differences between glacial and interglacial assemblages is the presence of byssal attachment, the trait modality of L. mediterraneum. When this species was removed from the dataset, this trait lost its importance for explaining variation in functional composition between samples (Supplementary Table S4). Moreover, without L. mediterraneum there was no clear distinction between glacial and interglacial intervals anymore (Supplementary Fig. S6), demonstrating how one highly abundant species can govern differences in functional traits between assemblages. While some might argue that relative abundances should be transformed before analysis to downweigh the effect of very dominant species, we deliberately opted against it, because very abundant species can be the drivers of ecosystem processes (Winfree et al. Reference Winfree, Fox, Williams, Reilly and Cariveau2015; Gaston et al. Reference Gaston, Cox, Canavelli, García, Hughes, Maas, Martínez, Ogada and Inger2018). Moreover, analyses based on transformed data still suggest distinct composition of the LG samples, although overlap with interglacial assemblages is higher (Supplementary Fig. S3).
Not all observed differences in functional composition are related to the most abundant species. Without L. mediterraneum, vegetation association and feeding type became the traits that best explain variation in functional composition (Supplementary Table S4). Plant- or algae-associated species, as well as herbivores and grazers, are more frequent and abundant in the LG assemblage, which includes a few samples that are strongly associated with these trait modalities (Fig. 2). Vegetated habitats (e.g., seagrass and macroalgal meadows) occur in shallow-marine to coastal environments of the present-day Adriatic and were even more common in the preindustrial Holocene (Lotze et al. Reference Lotze, Coll and Dunne2011). However, they are spatially patchy, especially in the northern part of the basin, as Mediterranean phanerogams and macroalgae prefer sheltered settings and areas not affected by strong river outflow and high turbidity (Short et al. Reference Short, Carruthers, Dennison and Waycott2007; Falace et al. Reference Falace, Alongi, Cormaci, Furnari, Curiel, Cecere and Petrocelli2010; Telesca et al. Reference Telesca, Belluscio, Criscoli, Ardizzone, Apostolaki, Fraschetti and Gristina2015). The detection of seagrass- or algae-associated assemblages in the LG, despite the lower number of samples representing that time interval, could reflect greater habitat heterogeneity, including more sheltered settings, in targeted nearshore environments during that time. This interpretation is consistent with paleoenvironmental reconstructions of the central Adriatic lowstand Po system, which included a wave-dominated delta positioned close to the shelf break, as well as a more complex coastal topography with large embayments located away from the main sediment entry point. Moreover, storm-enhanced hyperpycnal discharge coupled with offshore-directed flows likely promoted efficient redistribution of sediment during high-energy events, potentially limiting the persistence of fine-grained, plume-affected conditions in parts of the nearshore and enabling localized development of vegetated habitats in sheltered areas (Gamberi et al. Reference Gamberi, Pellegrini, Dalla Valle, Scarponi, Bohacs and Trincardi2020). Finally, because the seagrass Posidonia oceanica, endemic to the Mediterranean Sea, was likely climatically constrained during glacial conditions (Chefaoui et al. Reference Chefaoui, Duarte and Serrão2017), vegetated habitats during that time were likely dominated by other seagrasses (e.g., Zostera, Cymodocea) and/or macroalgae.
Another potential driver of changes in functional composition and diversity is temperature. Estimated SST was at least 6°C lower during the late LG than during the Holocene (Piva et al. Reference Piva, Asioli, Andersen, Grimalt, Schneider and Trincardi2008; Bazzicalupo et al. Reference Bazzicalupo, Sicre, Checa, Maiorano, Margaritelli, Klein and Pena2022), and even lower compared with the LIG (Piva et al. Reference Piva, Asioli, Andersen, Grimalt, Schneider and Trincardi2008). This had a profound impact on taxonomic composition of nearshore communities of the Adriatic Sea, as the relative abundance of species that are confined to warm-temperate regions today decreased substantially, whereas species reaching higher latitudes became much more common (Scarponi et al. Reference Scarponi, Nawrot, Azzarone, Pellegrini, Gamberi, Trincardi and Kowalewski2022). In concordance with the latitudinal biodiversity gradient, modern molluscan assemblages from cold-temperate environments are characterized by lower functional richness, measured as the number of distinct functional groups (Edie et al. Reference Edie, Jablonski and Valentine2018; Schumm et al. Reference Schumm, Edie, Collins, Gómez-Bahamón, Supriya, White, Price and Jablonski2019), a pattern that has also been found in the fossil record (e.g., Grossmann et al. Reference Grossmann, Nielsen, Rivadeneira and Valdivia2023). A potential explanation for this gradient is environmental filtering at higher latitudes, as lower productivity and harsher abiotic conditions, such as lower temperatures and greater seasonality, limit the set of viable functional profiles (Edie et al. Reference Edie, Jablonski and Valentine2018; Schumm et al. Reference Schumm, Edie, Collins, Gómez-Bahamón, Supriya, White, Price and Jablonski2019).
In contrast to the present-day latitudinal gradient, higher functional diversity during the colder LG than during the warmer interglacial periods was characteristic of molluscan assemblages in the Adriatic Sea. This pattern may reflect regional environmental filtering acting on a fauna shaped by Mediterranean geological and long-term climate history. Specifically, the Mediterranean molluscan fauna is dominated by species that can tolerate lower temperatures, as most warm-water species were eliminated at the transition from the Pliocene to the Pleistocene (Monegatti and Raffi Reference Monegatti and Raffi2001; Mondanaro et al. Reference Mondanaro, Dominici and Danise2024 and references therein). Consequently, few species disappeared from the Adriatic during the glacial period, although those that were dominant during interglacials but less adapted to colder climate, like the Mediterranean endemic L. mediterraneum, might have been impacted and shifted their distributions southward. The decreasing abundance and more restricted distribution of such species in the Adriatic might have enabled the establishment and spread of other species, including those with different functional profiles. Among the species that became more abundant during the last glacial were also many cosmopolitan taxa that today occur in temperate as well as the functionally more diverse (sub)tropical regions (Scarponi et al. Reference Scarponi, Nawrot, Azzarone, Pellegrini, Gamberi, Trincardi and Kowalewski2022). To conclude, it is plausible that a combination of a climate that was less favorable to dominant species of the interglacials and higher habitat heterogeneity due to basin morphology led to communities with higher taxonomic evenness during the last glacial, which in turn translated into higher functional diversity in targeted nearshore settings.
Natural Range of Functional Variability as a Baseline for Assessing Future Change
While the LG assemblage differs in functional composition, the LIG and CIG assemblages show higher overlap (Fig. 2). The LIG assemblage represents a subset of the functional diversity of the CIG assemblage (Fig. 5C, Supplementary Table S7), although this is probably due to its much smaller sample size. This high functional similarity of the CIG assemblage—which represents baseline community states of the preindustrial Holocene—to the molluscan fauna of the LIG, has important implications for forecasting future responses of Adriatic coastal ecosystems, as the LIG is considered an analogue for Earth’s future climate amid global warming (Fischer et al. Reference Fischer, Meissner, Mix, Abram, Austermann, Brovkin and Capron2018). It was estimated that SSTs in the Adriatic were ~2–3°C higher than present-day SSTs (Siani et al. Reference Siani, Magny, Paterne, Debret and Fontugne2013; Santagati et al. Reference Santagati, Perri, Bernasconi, Borrelli, Guerrieri and Critelli2024; De Santis et al. Reference De Santis, Cerrone, Meschis, Scicchitano, Ascione and Caldara2025). Even though the Mediterranean is considered a climatic hotspot, warming of the mean SST in the central Adriatic should be constrained to less than 2°C (with respect to preindustrial time) until the year 2100 (Shaltout and Omstedt Reference Shaltout and Omstedt2014) under low to intermediate emission scenarios RCP 2.6–4.5 (IPCC 2023). Our results suggest that such levels of warming alone are unlikely to induce strong changes— comparable to the turnover at glacial–interglacial transitions—in the ecosystem functioning of molluscan communities in the fluvially influenced settings investigated here, because trait structures of CIG and LIG assemblages showed very high similarity despite higher SSTs during the latter interval.
Our results contrast with the findings from different environmental and geographic settings. Ávila et al. (Reference Ávila, Arruda, Silva, Hipólito, Quartau, Ramalho and Madeira2025) compared last interglacial (MIS 5e) and present-day shallow-marine molluscan assemblages from the Azores and found that most functional diversity indices were higher in the modern assemblage, whereas functional divergence was higher in the MIS 5e assemblages. They interpreted these results as indicating that modern molluscan communities occupy a greater number of ecological niches and show greater niche specialization, while MIS 5e communities showed more efficient use of available resources—implying more oligotrophic settings—but lower resilience to invasions and climate change. Shallow-water settings in the Azores, which are dominated by hard substrates and lack major fluvial influence, differ strongly from the soft-bottom settings of the western Adriatic that we studied. The comparison between different settings shows that functional responses to warming are not uniform across shallow-marine systems and can depend on the environmental context. However, studies on Late Quaternary changes in functional diversity are lacking for most marine settings and taxonomic groups. Existing studies (e.g., Pandolfi Reference Pandolfi1996; Tager et al. Reference Tager, Webster, Potts, Renema, Braga and Pandolfi2010; Scarponi et al. Reference Scarponi, Nawrot, Azzarone, Pellegrini, Gamberi, Trincardi and Kowalewski2022) mainly focus on taxonomic composition. To better understand how long-term climate change affects marine ecosystem functioning, more studies applying BTA to the Late Quaternary fossil record are needed.
As for targeted ecosystems, similarity between the LIG and preindustrial Holocene suggests that, under natural conditions, the functional composition and diversity of Adriatic nearshore molluscan assemblages are resistant to warming of comparable magnitude to LIG conditions. Whether this remains true under the faster pace of anthropogenic warming and under additional human pressures requires direct comparison with functional data from modern settings. Indeed, coastal areas also vary in the nature and intensity of other anthropogenic stressors, with major impacts on current and future trajectories of ecosystem functioning. The Adriatic has a long history of human influence, which intensified during the second half of the twentieth century (Lotze et al. Reference Lotze, Coll and Dunne2011). Among the most relevant anthropogenic impacts in the region are bottom trawling (Eigaard et al. Reference Eigaard, Bastardie, Hintzen, Buhl-Mortensen, Buhl-Mortensen, Catarino and Dinesen2017; Pitcher et al. Reference Pitcher, Hiddink, Jennings, Collie, Parma, Amoroso and Mazor2022) and hypoxic events (Justić Reference Justić1991; Giani et al. Reference Giani, Djakovac, Degobbis, Cozzi, Solidoro and Umani2012), which can lead to changes in functional diversity. Bottom trawling can change functional composition of benthic communities by favoring specific trait modalities such as high mobility, infaunality, and scavenging (Tillin et al. Reference Tillin, Hiddink, Jennings and Kaiser2006; de Juan et al. Reference de Juan, Thrush and Demestre2007). Eutrophication and hypoxia can lead to decreased functional diversity and an increasing dominance of species with opportunistic traits (Tomašových et al. Reference Tomašových, Albano, Fuksi, Gallmetzer, Haselmair, Kowalewski, Nawrot, Nerlović, Scarponi and Zuschin2020; Bon et al. Reference Bon, Grall, Gusmao, Fajardo, Harrod and Pacheco2021; Haselmair et al. Reference Haselmair, Gallmetzer, Tomašových, Wieser, Übelhör and Zuschin2021; Harnik et al. Reference Harnik, Chao, Collins and Rillo2024). Such shifts in response to bottom trawling, eutrophication, and hypoxic events have occurred in molluscan offshore communities of the northern Adriatic mainly during the latest centuries. Epifauna declined and infauna became dominant, while feeding types such as grazers and herbivores decreased in abundance, leading to a loss in habitat heterogeneity (Gallmetzer et al. Reference Gallmetzer, Haselmair, Tomašových, Mautner, Schnedl, Cassin, Zonta and Zuschin2019; Haselmair et al. Reference Haselmair, Gallmetzer, Tomašových, Wieser, Übelhör and Zuschin2021). In a study comparing Holocene assemblages studied herein, which are all pre–1750 CE, to surficial death assemblages of present-day nearshore environments, Scarponi et al. (Reference Scarponi, Rojas, Nawrot, Cheli, Kowalewski, Nawrot, Dominici, Tomašových and Zuschin2023) found a strong imprint of anthropogenic impacts. In upper nearshore settings (0–4 m depth), a higher abundance of hard-substrate species was observed, which could be related to the construction of coastal infrastructure providing new forms of substrate. In lower nearshore settings (4–10 m depth) L. mediterraneum has decreased in abundance, and the bivalve Varicorbula gibba, an indicator species of stressed environments (Hrs-Brenko Reference Hrs-Brenko2006; Nerlović et al. Reference Nerlović, Doğan and Hrs-Brenko2011; but see also Leshno et al. Reference Leshno, Benjamini and Edelman-Furstenberg2016), and Chamelea gallina have become more widespread. The standardized species richness decreased compared with the preindustrial Holocene (CIG). This could be the result of a multitude of anthropogenic stressors, like hydraulic dredging by C. gallina fisheries, which leads to high sediment disturbance (Grazioli et al. Reference Grazioli, Guerranti, Pastorino, Esposito, Bianco, Simonetti, Rainis, Renzi and Terlizzi2022). Future studies could follow a multi-trait approach, as we did here, to determine whether the observed ecological turnovers have led to changes in functional composition and diversity.
The Adriatic has also been impacted by the spread of nonindigenous species, including invasive species that affect ecosystem functioning (Servello et al. Reference Servello, Andaloro, Azzurro, Castriota, Catra, Chiarore and Crocetta2019; Slišković et al. Reference Slišković, Piria, Nerlović, Ivelja, Gavrilović and Mrčelić2021). Among many examples, in the Po River delta, the alien bivalve Arcuatula senhousia rapidly proliferated and spread as an opportunistic colonizer forming dense byssal mats that depress the abundance of other macrobenthic species, causing a shift from a primarily deposit- to suspension-feeding community (Mistri Reference Mistri2003). A recent example in the study area is the strong impact of the Atlantic blue crab (Callinectes sapidus) on epifaunal and shallow-infaunal bivalves, which had detrimental effects on clam fisheries (Gaglio et al. Reference Gaglio, Gavioli, Turolla, Lanzoni and Castaldelli2025; Gavioli et al. Reference Gavioli, Castaldelli, Eggleston and Christian2025). Currently, the majority of the nonindigenous species in the Adriatic are cold-water species (Servello et al. Reference Servello, Andaloro, Azzurro, Castriota, Catra, Chiarore and Crocetta2019; Giacobelli et al. Reference Giacobelli, Tanduo, Orlandi, Turolla and Crocetta2025), but this is expected to change due to the northward spread of warm-water species in the Mediterranean (Albano et al. Reference Albano, Schultz, Wessely, Taviani, Dullinger and Danise2024). At the same time, range contractions and shifts are expected for many native species, including key taxa of nearshore settings like C. gallina (Gallagher and Albano Reference Gallagher and Albano2023; Schultz et al. Reference Schultz, Wessely, Dullinger and Albano2024). However, being the northernmost subregion of the Mediterranean, the Adriatic Sea could act as a climate refugium (Gallagher and Albano Reference Gallagher and Albano2023; Baldan et al. Reference Baldan, Chauvier‐Mendes, Panzeri, Cossarini, Solidoro and Bandelj2026). Nonetheless, under very strong global warming, a high species turnover in Adriatic benthic communities is likely to occur before the end of the century (Baldan et al. Reference Baldan, Chauvier‐Mendes, Panzeri, Cossarini, Solidoro and Bandelj2026). This has already happened in the southeastern Mediterranean, where many species from the Red Sea have entered the Mediterranean through the Suez Canal to become dominant elements of the shallow-shelf fauna (Rilov Reference Rilov2016; Albano et al. Reference Albano, Steger, Bošnjak, Dunne, Guifarro, Turapova, Hua, Kaufman, Rilov and Zuschin2021). It was found that the nonindigenous component of molluscan assemblages in this area was functionally different from the native one, even when historical death assemblages are considered. Thus, the high influx of nonindigenous species has led to the formation of novel ecosystems with unique functioning (Steger et al. Reference Steger, Bošnjak, Belmaker, Galil, Zuschin and Albano2022, Reference Steger, Bogi, Lubinevsky, Galil, Zuschin and Albano2024). Under extreme emission scenarios, a similar restructuring of ecosystem functioning could occur in Adriatic nearshore communities, pushing their functional composition and diversity outside the natural range of variability of the Late Quaternary.
Conclusions
We applied biological traits analysis to test whether the functional composition and diversity of molluscan fluvially influenced nearshore communities in the Adriatic changed in response to Late Quaternary climatic fluctuations. We found that the assemblage of the last glacial differed in composition from those of the last interglacial and the Holocene and displayed higher functional diversity. This is the result of a more diverse and even species assembly in nearshore settings that probably formed in response to temperature changes and increased habitat diversity. On the other hand, any differences in environmental conditions between the interglacials have not led to changes in functional diversity, as high similarity between assemblages from the last interglacial and the Holocene were observed. However, under more extreme scenarios, high species turnover involving the range contraction of native species and the spread of nonindigenous species could alter functional composition substantially. Furthermore, other anthropogenic impacts that the Adriatic has been experiencing for decades to centuries might also have had strong effects on functional diversity, and more studies using a trait-based approach are needed to discern how much ecosystems have already changed. In conjunction with the natural baseline for functional responses to long-term climate change established by our study, such insights would prove valuable to predict and guide future ecosystem functioning in these crucial coastal ecosystems.
Acknowledgments
The authors thank P. Albano and F. Ye (both Stazione Zoologica Anton Dohrn) for fruitful discussion on the methods of BTA applied in this study. This research was supported by the European Union–NextGenerationEU through the Italian Ministry of University and Research under Piano Nazionale di Ripresa e Resilienza (PNRR): Mission 4 Component C2, Investment 1.1 “Conservation of Life on Earth: The Fossil Record as an Unparalleled Archive of Ecological and Evolutionary Responses to Past Warming Events” code 2022WEZR44. The data used here were partly derived from previous projects funded by the U.S. National Science Foundation grants EAR-1559196 and EAR-0920075.
Author Contribution
L.S., R.N., M.K., and D.S. were involved in the conceptualization, methodology and visualization of the study. L.S., V.P., and D.S. were involved in the assembly of the trait dataset. L.S., R.N., M.K., and D.S. implemented the study and wrote the initial draft. All authors reviewed and contributed to the final version of the manuscript and Supplementary Material.
Competing Interests
The authors declare no competing interests.
Data Availability Statement
Supplementary Material, including text, figures, and tables supporting the conclusions of this study, as well as the datasets and R scripts needed to reproduce all analyses, figures and tables, is available on Dryad at: https://doi.org/10.5061/dryad.j0zpc86xg and on Zenodo at: https://doi.org/10.5281/zenodo.21641460.
Open access
