Ecological baselines in the Eastern Mediterranean Sea shifted long before the availability of observational time series

Native biodiversity loss and invasions by nonindigenous species (NIS) have massively altered ecosystems worldwide, but trajectories of taxonomic and functional reorganization remain poorly understood due to the scarcity of long‐term data. Where ecological time series are available, their temporal coverage is often shorter than the history of anthropogenic changes, posing the risk of drawing misleading conclusions on systems' current states and future development. Focusing on the Eastern Mediterranean Sea, a region affected by massive biological invasions and the largest climate change‐driven collapse of native marine biodiversity ever documented, we followed the taxonomic and functional evolution of an emerging “novel ecosystem”, using a unique dataset on shelled mollusks sampled in 2005–2022 on the Israeli shelf. To quantify the alteration of observed assemblages relative to historical times, we also analyzed decades‐ to centuries‐old ecological baselines reconstructed from radiometrically dated death assemblages, time‐averaged accumulations of shells on the seafloor that constitute natural archives of past community states. Against expectations, we found no major loss of native biodiversity in the past two decades, suggesting that its collapse had occurred even earlier than 2005. Instead, assemblage taxonomic and functional richness increased, reflecting the diversification of NIS whose trait structure was, and has remained, different from the native one. The comparison with the death assemblage, however, revealed that modern assemblages are taxonomically and functionally much impoverished compared to historical communities. This implies that NIS did not compensate for the functional loss of native taxa, and that even the most complete observational dataset available for the region represents a shifted baseline that does not reflect the actual magnitude of anthropogenic changes. While highlighting the great value of observational time series, our results call for the integration of multiple information sources on past ecosystem states to better understand patterns of biodiversity loss in the Anthropocene.


| INTRODUC TI ON
Human activities have profoundly reshaped global biodiversity by driving unprecedented losses of native species and promoting the introduction of nonindigenous ones (IPBES, 2019;McCauley et al., 2015;Seebens et al., 2017).Although no region is left in an entirely pristine state (Halpern et al., 2008), some systems have been modified to a degree that renders returning to, or the restoration of, historical baseline conditions unfeasible, so-called "novel ecosystems" (Hobbs et al., 2013;Morse et al., 2014).The number and areal extent of such massively altered systems are expected to increase globally (e.g., Bianchi, 2007;Morse et al., 2014;Vergés, Steinberg, et al., 2014); however, the limited availability of long time series often prevents tracking, quantifying, and fully understanding the trajectories of taxonomic and functional reorganization, as well as their consequences for ecosystem functioning (Morse et al., 2014;Peleg et al., 2020;Steger et al., 2022;Vergés, Steinberg, et al., 2014).This is particularly problematic as conservation and restoration targets are still set more according to our perception of pristineness-itself often based on already degraded systems-rather than on objective information on ecosystem states before the onset of major human pressures, the socalled "shifting baseline syndrome" (Pauly, 1995).
The southeastern Mediterranean is probably the most striking example of an emerging novel marine ecosystem.It is the hottest sector of the Mediterranean Sea and represents the trailing distributional edge of native Atlanto-Mediterranean species, which often live close to their thermal tolerance limits and thus tend to be particularly sensitive to anthropogenic stressors (Bates et al., 2014;Rilov, 2016;Yeruham et al., 2015).For many of these organisms, the naturally extreme environment poses physiological challenges, reflected in often significantly reduced adult body sizes compared to Western Mediterranean populations (Comay et al., 2012;Por, 1989;Sharir et al., 2011;Sonin et al., 2007).Human-mediated climate warming has further exacerbated the abiotic regime in recent decades (Ibrahim et al., 2021;Ozer et al., 2017;Rilov, 2016): in the most affected parts of the basin, seawater temperatures have been rising at a rate 4-10 times the global ocean average, whichtogether with increasingly frequent and severe marine heatwaves (Ibrahim et al., 2021)-has caused the largest regional-scale collapse of native biodiversity in the modern oceans (Albano, Steger, Bošnjak, et al., 2021;Ozer et al., 2017;Rilov, 2016).Additionally, the opening of the Suez Canal approximately 150 years ago enabled the massive introduction of nonindigenous species (NIS) from the Red Sea, the so-called "Lessepsian invasion" (Por, 1978;Rilov & Galil, 2009).
The establishment and spread of thermophilic NIS are currently further facilitated by warming (Albano, Steger, Bošnjak, et al., 2021;Amarasekare & Simon, 2020;Givan et al., 2018) and have led to a progressive "tropicalization" of shallow-water ecosystems in the region (Bianchi, 2007;Bianchi & Morri, 2003;Peleg et al., 2020;Vergés, Steinberg, et al., 2014).Jointly, these processes resulted in the emergence of an unprecedented mixture of native and nonindigenous biotas and of ecological phase shifts that are predicted to expand to increasingly larger parts of the Mediterranean Sea (Bianchi, 2007;Vergés, Steinberg, et al., 2014;Vergés, Tomas, et al., 2014).Still, the timing and dynamics of these changes remain poorly known.
Here, using a unique 18-year-long time series of shelled mollusks from the shallow subtidal along the entire c. 180-km-long Israeli coastline, we follow the temporal development of individual density, assemblage composition, taxonomic richness, and functional diversity of this rapidly evolving novel ecosystem.Mollusks, in this context, are a particularly informative target as they are ubiquitous, taxonomically and ecologically diverse, and good surrogates for entire benthic assemblages (Smith, 2005;Tyler & Kowalewski, 2017).
We aim at pinpointing the timing and dynamics of the native biodiversity collapse in the region, at assessing trait overlap-and thus the potential for resource competition-between native species and NIS, and at quantifying changes in ecosystem functional properties.To determine to which extent the observational assemblage data represented a shifted ecological baseline, we further reconstructed decades-to centuries-old baselines from radiometrically dated molluscan death assemblages, the time-averaged accumulations of taxonomically identifiable empty shells on the seafloor (Kidwell, 2007;Kidwell & Tomašových, 2013).Due to the great durability of these calcareous remains and generally low sedimentation rates in open shelf environments, death assemblages integrate skeletal elements of multiple non-contemporaneous generations over decadal to millennial timescales (Kidwell, 2009(Kidwell, , 2013)).This significant temporal mixing levels out the short-term demographic fluctuations and small-scale patchiness typical of living populations, and renders death assemblages compositionally inert to recent directional changes in living assemblage composition, making them valuable natural archives of past community states (see Kidwell & Tomašových, 2013 for a review).The consideration of information from the young geohistorical record enabled us to considerably extend the temporal scope of our analysis, and thereby obtain important quantitative baseline data on past assemblages that otherwise would have remained elusive.Our results thus not only provide a better understanding of the temporal development of the novel Eastern Mediterranean benthic ecosystem but also highlight the value of integrating multiple sources of ecological data when studying biodiversity trajectories in the Anthropocene.Triplicate sediment samples were taken annually in July-August at 13 stations located at 8-13 m water depth (Figure 1; Table S1 in Supporting Information S1), using a 0.11 m 2 van Veen grab (KahlSico WA265/SS214, 32 × 35 cm, 20 L).Seven of these stations were discontinued after 2016 (Figure 1), following a reallocation of resources to also include sampling in deeper water (data not considered herein).After collection, bulk sediments were sieved over a 0.25-mm mesh and the retained material preserved in 4% formaldehyde solution (years 2005-2012) or 99% ethanol (years 2013-2022, to enable potential future barcoding of specimens).

| Study area, field sampling, and assemblage datasets
In the laboratory, samples were stained with ethanolic Rose Bengal (until 2012) or eosin solutions (since 2013), live-collected individuals picked under a stereomicroscope, identified to the lowest possible taxonomic level based on shell morphology, and counted (Lubinevsky et al., 2019).This study is based on the resultant 35,752 specimens representing 110 species, of which 66 were native, 39 nonindigenous, and 5 had an uncertain native/ nonindigenous status (Supporting Information S2).
To quantify the taxonomic and functional richness of the novel molluscan assemblage relative to historical communities, we analyzed data on two death assemblages collected in 2016 by an independent field campaign nearby our study sites (see Albano, Steger, Bošnjak, et al., 2021;Steger et al., 2022).One death assemblage was sampled off Ashqelon (station SG10 of Steger et al., 2022, 1).While available death assemblages had not been collected directly at our study sites (Figure 1), the biotic homogeneity of the studied sandy bottom habitat along the c. 130 km-long shelf stretch between Ashqelon and Haifa Bay (Lubinevsky et al., 2019), and the reduction in small-scale compositional variability by time-averaging (Tomašových & Kidwell, 2009) ensured they were fully representative for the close-by National Monitoring stations.Death assemblages were obtained by sieving grab-collected sediment samples over a 0.5-mm mesh, and picking quantitative splits of the retained material under a stereomicroscope until approximately 1000 shell elements had been extracted.Only fragments representing more than half of a complete shell element, with either aperture or apex (in gastropods) or the hinge (in bivalves) preserved, were considered, to avoid counting any of them more than once.Raw counts of loose bivalve valves were divided by 2, and those of polyplacophoran plates by 8, to correct for the greater contribution of shell elements by these taxa compared to gastropods and scaphopods.The shell age distribution of each death assemblage was assessed by radiocarbon-dating 15 randomly selected valves of F I G U R E 1 Map of the 13 National Monitoring stations on the Mediterranean Israeli shelf.Grey dots indicate stations discontinued after year 2016, black dots those sampled until 2022.Orange crosses mark death assemblage sampling sites NG10 and SG10 of Albano, Steger, Bošnjak, et al. (2021) and Steger et al. (2022).Station codes and localities (north to south): HM2.1-Akko; HM10-Kfar Masarik; HM23.1-KiryatHayim; HM27-Kishon-Haifa port; H3-Dado (off Haifa); NG10-off Haifa; H7-Taninim; H11-Alexander River mouth; H41-Poleg; H13-Herzliya; H16-Yarkon River mouth; H19-Soreq; H24-Ashdod; and SG10 and H28-Ashqelon.For further details and geographic coordinates of sampling stations, see Table S1.Map lines delineate study areas and do not necessarily depict accepted national boundaries.a common native bivalve.Detailed information on the selected species and dating-including laboratory protocols and data calibration procedures, as well as ages of individual valves-are available from Albano et al. (2020), Albano et al. (2021, Supplemental Material) and Steger et al. (2022, Supporting Information Appendix S1).The native death assemblage component at Ashqelon had a median age of 763 calendar years before sample collection (i.e., AD 2016).Further, more than 53% of the dated shells were "pre-Lessepsian", that is, older than AD 1869-when the Suez Canal was opened-confirming it represented a solid baseline to assess changes in benthic assemblages that occurred during the 20th century (Steger et al., 2022).
The death assemblage collected off Haifa, in contrast, had a much younger median age of 24 years before sample collection and contained only 7% pre-Lessepsian shells (Steger et al., 2022).While such a death assemblage does not reflect preinvasion baseline conditions, it can-together with the much older and geographically distant sample from Ashqelon-provide important insights into the timing and spatial extent of any shifts in ecological baselines.Specifically, if both live-dead comparisons yield qualitatively similar results, this suggests not only that similar biodiversity changes affected large portions of the shelf but also that they (largely) happened during the 20th century, that is, after the onset of major anthropogenic impacts, rather than in response to mostly natural environmental changes in the more distant past.

| Trait dataset
We functionally characterized molluscan species based on five different traits: maximum adult body size, feeding habit, environmental position, substrate affinity, and host association.
These traits reflect several important aspects of molluscan ecology (Nawrot et al., 2015), and reliable information is available for most species in the Eastern Mediterranean (Steger et al., 2022;Steger, Dunne, et al., 2021).Details on the ecological relevance of the selected traits and definitions of scored modalities (i.e., trait categories) are available from Table S2.Trait information was compiled from the literature, online databases, observations on specimens, and expert knowledge-a complete reference list is available in Supporting Information S3.We used a semi-quantitative fuzzy coding scheme for trait scoring (Chevenet et al., 1994) that takes into account that species may have affinities to more than one modality per trait.We employed the widely used system with discrete numeric scores from 0 to 3, where 0 corresponds to "no affinity", 1 to "low affinity", 2 to "high affinity", and 3 to "exclusive affinity" (Bremner et al., 2003;Degen et al., 2018).In the case of similar affinity to two or more modalities of a trait, each was scored with "2" by convention.If no information on a particular trait was available for a species (concerning 2 out of 224 species, or 0.9% of taxa), all of its modalities were assigned a score of "0" to avoid introducing biases into subsequent analyses (Chevenet et al., 1994).
Prior to analysis, raw modality scores of each trait were converted to proportions adding to 1; this way, traits were weighted equally, irrespective of the number of scored modalities (Oug et al., 2012).
The standardized species × trait modality matrix is available in Supporting Information S4.

| Data analysis
We followed the regional-scale temporal evolution of molluscan living assemblages on the Israeli shelf by assessing developments in (i) individual density; (ii) taxonomic and trait composition; and (iii) taxonomic and functional diversity.To obtain larger, statistically more robust sample sizes, we merged data from the triplicate grab samples.
Analyses were conducted for complete molluscan assemblages, and separately for their native and nonindigenous components, in the R statistical programing environment (v.3.5.2;R Core Team, 2018).
Patterns of taxonomic and trait composition were visualized by nonmetric multidimensional scaling (nMDS, based on Bray-Curtis dissimilarities of square root-transformed relative abundance data) and fuzzy correspondence analysis (FCA; Chevenet et al., 1994), respectively, after removing nonrepresentative samples with fewer than 10 individuals.To obtain assemblage-level trait profiles for ordination with FCA, we weighted species' standardized trait scores [modalities × species table] by their relative abundances [species × assemblage table] using matrix multiplication (e.g., Oug et al., 2012); trait modalities not represented in any sample were removed before ordination (see Table S3 for details).To check which traits' modalities were best separated along FCA axes, we calculated correlation ratios, which represent, for each axis, the proportion of total variance explained by the individual traits (Chevenet et al., 1994).NMDS and FCA were performed using the "vegan" (v.2.5-4; Oksanen et al., 2019) and "ade4" (v. 1.7-13;Dray & Dufour, 2007) R packages, respectively.Differences in assemblage taxonomic composition among years were tested by PERMANOVA (Bray-Curtis distances, 9999 permutations, variable "station" defined as strata), and subsequent pairwise comparisons performed using the "pairwiseAdonis" function (v.0.4; Martinez Arbizu, 2020) with Holm-Bonferroni corrected p-values to account for multiple testing (Holm, 1979).PERMANOVA (Euclidean distances, 9999 permutations, "station" as strata) was further used to assess differences in trait composition among sampling years, as well as between native and nonindigenous assemblage components-a proxy of niche and resource use similarity-using their coordinates on the first three FCA axes.
As our focus was on regional biodiversity trends, we calculated taxonomic and functional diversity indices at the scale of the entire Israeli coastline, that is, on yearly data pooled across stations, considering only those sites monitored until 2022 to account for changed sampling effort post-2016 (see Figure 1).Taxonomic diversity of assemblages was quantified by calculating species richness at standardized sample completeness, or coverage (SC; Chao & Jost, 2012), using the iNEXT R package (v. 2.0.20;Hsieh et al., 2020).While observed species richness-even at constant sampling effort in terms of seafloor area-strongly depends on the number of collected individuals and species' abundance distributions (Chao et al., 2014;Chao & Jost, 2012), the use of annual estimates at equal coverage enables an unbiased assessment of diversity trajectories through time.For standardization, we used the coverage value obtained when extrapolating the least complete sample to double the observed number of individuals (SC = 0.97 for complete assemblages, SC = 0.93 for native and nonindigenous components).In contrast to standardizing to the lowest observed coverage value, this approach minimizes information loss by avoiding overly strong rarefaction of high-coverage years, while capitalizing on moderate and statistically reliable extrapolation of lower coverage samples (Hsieh et al., 2016).
Functional diversity was quantified based on the distribution of species and their relative abundances in functional (or trait) space, a multidimensional coordinate system where species are arranged by their trait similarity (Villéger et al., 2008(Villéger et al., , 2013)).To construct this space, we first grouped the species-specific trait profiles of all 224 molluscan taxa occurring in our living and death assemblage samples into 114 functional entities (i.e., unique trait profiles).We then calculated generalized Gower dissimilarities (Pavoine et al., 2009) among the entities, and performed principal coordinate analysis (PCoA) to convert these dissimilarities into Euclidean trait space coordinates, which were finally reassigned to the individual species.We built and assessed the quality of functional spaces with up to 10 dimensions using a modified version of the "quality_funct_space" R function (Maire et al., 2015), developed and scripted by S. Villéger and available from Steger, Bošnjak, et al. (2021).A four-dimensional space was selected as it represented the best trade-off between (i) faithfully representing the original trait distances among species (mean squared deviation between Gower and Euclidean distances = 0.012; see Maire et al., 2015); (ii) enabling the calculation of functional diversity indices also for years and assemblage components with a limited number of species at a given (standardized) coverage level; and (iii) the computational effort required to calculate indices for resampled assemblage datasets (next paragraph).
As functional diversity encompasses species' trait profiles and abundances, we employed a resampling strategy (1000 iterations) to obtain coverage-standardized index values for assemblages and their native and nonindigenous components throughout the time series (see Supporting Information S1, Section 1.5, for details).Based on the resampled assemblage data and the species-specific coordinates in functional space, we calculated three complementary functional alpha diversity indices: functional richness (the volume of functional space occupied by an assemblage), functional divergence (an indicator of the proportion of total abundance contributed by species with extreme trait values), and functional evenness (a quantifier of the regularity of species and abundance distribution in functional space; Mouillot et al., 2013;Villéger et al., 2008).While the latter two measures are constrained between 0 and 1 by definition, functional richness values are reported as proportions of the functional space volume occupied by the entire species pool involved in the analysis, thus enabling their intuitive interpretation (Magneville et al., 2022;Villéger et al., 2008; for more information on functional diversity indices and scaling, see also Supporting Information S1, Section 1.6).Additionally, to analyze the patterns of functional distinction between native and nonindigenous assemblage components of each year, we computed their Jaccard-like functional beta diversity, including its decomposition into additive turnover and nestedness-resultant components (Villéger et al., 2013; see Supporting Information S1, Section 1.6, for further information).All functional diversity indices were calculated using the "mFD" R package (v.1.0.3;Magneville et al., 2022) and are reported as mean values with associated 95% confidence intervals, calculated using the percentile method (Kowalewski & Novack-Gottshall, 2010).
Temporal trends in taxonomic and functional diversity across the entire study period were assessed by generalized least-squares (GLS) regression ("gls" function in the "nlme" R package, v. 3.1-137; Pinheiro et al., 2018), which allows accounting for autocorrelated errors and/ or residual heteroscedasticity.Without explicit specification of error or residual variance structures, GLS models are equivalent to standard ordinary least-squares regression (OLS; Zuur et al., 2009).All models were fitted using restricted maximum likelihood (REML) estimation; a detailed description of fitting and selection procedures is provided in Supporting Information S1, Section 1.7.Correlations among different diversity indices, for example, the trajectories of taxonomic and functional richness, were assessed by calculating Spearman's rho for time series data detrended by first-order differencing (cf.Dornelas et al., 2013;Yaffee & McGee, 2000), to avoid spurious results related to the presence of temporal trends.
Variables representing ratios (e.g., native vs. nonindigenous density or species richness) were log 2 transformed to enhance the visual interpretability of plots, with positive ratios indicating greater values of the native versus nonindigenous assemblage component, and negative values vice versa.
For our comparisons between the modern and historical (i.e., death) assemblages, we merged the yearly living assemblage data (2005)(2006)(2007)(2008)(2009)(2010)(2011)(2012)(2013)(2014)(2015)(2016)(2017)(2018)(2019)(2020)(2021)(2022) from National Monitoring stations H28 (off Ashqelon) and H3 (off Haifa) into single, temporally integrated datasets, thereby simulating time-averaging, an important feature of death assemblages (Kidwell, 2013;Kidwell & Tomašových, 2013).As historically the Israeli shelf was inhabited exclusively by Mediterranean species, only the native components of the SG10 and NG10 death assemblage datasets of Steger et al. (2022) were kept for analysis (see Supporting Information S2).To avoid an artificial inflation of (functional) diversity related to the possible presence of "outof-habitat" shells, that is, those transported into the study site by hermit crabs, storms, waves, and/or currents, we followed a conservative approach and excluded all species predominantly associated with hard bottoms, judged by their substrate affinity profiles (see Supporting Information S4; details on the applied filters are provided in Supporting Information S1, Section 1.8).We further removed records of shells of the gastropod Pirenella conica (Blainville, 1829) from the death assemblage data, as in Israel this species only lives in hypersaline coastal lagoons and salt water ponds (Barash & Danin, 1992;Taraschewski & Paperna, 1981), and thus never inhabited the open shelf.To ensure a "fair" live-dead comparison of taxonomic and functional richness, hard substrate-associated species were also deleted from the H28 and H3 living assemblage datasets.
Following these filtering steps, we determined the taxonomic coverage of both datasets, and used the lower of the two values obtained (SC = 0.96 and 0.92 for the Ashqelon and Haifa death assemblages, respectively) as the benchmark for repeatedly (n = 1000) resampling the more complete living assemblages (SC = 1.00 and 0.99 for Ashqelon and Haifa, respectively) by the algorithm for coverage standardization outlined in Supporting Information S1, Section 1.5.
Finally, based on the adjusted datasets, we calculated and compared the taxonomic and functional richness of corresponding living and death assemblages.

| Trends in molluscan density
Molluscan density at the six continuously sampled National Monitoring stations on the Israeli shelf increased between 2005 (median: 117 inds.m −2 ) and 2015 (median: 1147 inds.m −2 ; Figure 2a).The median log 2 ratio of native versus nonindigenous density showed a significant global decline over time (GLS model for time series of median log 2 ratio values: slope = −0.11,p = .01),indicating an increasing share of nonindigenous individuals in molluscan assemblages (Figure S1).The ratio reached its minimum in 2020 (median: −0.58), followed by a recovery in the last two years of the time series, when median values ranged between 1.53 and 1.63 (Figure S1).
However, it remains to be seen if this development merely is a short-term fluctuation or marks the onset of a sustained longer-term increase.

| Trajectories in assemblage composition and taxonomic richness
Temporal patterns of assemblage composition as revealed by nMDS were rather weak (Figure S2) but statistically significant, though only little variance in the data was explained by the variable "sampling year" (PERMANOVA, pseudo-F = 12.90, p < .001,R 2 = .06).Samples from the beginning of the time series (2005)(2006)(2007), when log ratios of native versus nonindigenous species richness and density were particularly high (Figures S3 and S1,respectively), clustered at the periphery of the cloud comprised by later assemblages (Figure S2).Pairwise comparisons among sampling years revealed no directional change in assemblage composition over time, but confirmed that samples from 2005 and 2006 were similar to each other (PERMANOVA, pseudo-F = 2.95, p adj = .07,R 2 = .11)but distinct from almost all other years (PERMANOVA, p adj < .05for all except one comparison; Table S4).
Coverage-standardized (SC = 0.97) taxonomic richness of the regional-scale molluscan assemblage significantly increased over the 18-year study period (GLS, p < .01, Figure 3a), with a slope corresponding to the addition of 0.82 (SE: 0.27) species per year (see Table S5 for model details).A major short-term drop in richness, however, occurred in 2010 due to a transient peak in density and relative abundance of two species of mactrid bivalves that strongly depressed coverage-standardized estimates (Figure 3a).Considering native and nonindigenous assemblage components separately (at SC = 0.93), qualitatively different regional-scale temporal patterns emerged: while native species richness did not show a significant long-term trend (GLS, p = .32;Table S5), an increase was observed for the regional-scale nonindigenous assemblage component (GLS, p < .001; Figure 3b; Table S5), which drove the richness trajectory of the complete molluscan assemblage (Spearman's rho for differenced time series of NIS and complete assemblage richness = 0.65, p < .01).
This development is largely reflected also by the ratio of native versus nonindigenous species richness: until 2008, the regional-scale native assemblage component had a considerably higher coveragestandardized species richness than the nonindigenous one (log ratios >0, Figure S3), followed by a decline and stabilization at values fluctuating around zero (i.e., equal richness) by 2011.Again, a major steep decline of this ratio occurred in 2010 due to the abovementioned recruitment peak of juvenile bivalves that was largely dominated by individuals of the native species Mactra stultorum (Linnaeus, 1758) (Figure S3).

| Trait composition and functional diversity
Abundance-weighted trait profiles of molluscan assemblages from the 13 National Monitoring stations showed high temporal overlap in the FCA ordination plot (Figure S4), although significant differences among sampling years were detected by PERMANOVA (pseudo-F = 10.41,p < .001,R 2 = .05).In contrast, there was a clear functional segregation-mainly along the first FCA axis-between the native and nonindigenous components of assemblages (Figure 4 The first two axes of the global FCA explained 25.2% and 16.4% of the variation in the data, respectively (Figure 4), with the modalities of traits "maximum adult body size" (correlation ratios of 0.43 and 0.37 for axes 1 and 2, respectively), "feeding habit" (correlation ratios of 0.41 and 0.32 for axes 1 and 2, respectively), and "environmental position" (correlation ratio of 0.32 for axis 1 (but only 0.003 for axis 2)) being best separated (see Table S6 for further details on correlation ratios).Biplots depicting the association of these traits' modalities with assemblage component trait profiles revealed that small adult body size was an important feature distinguishing nonindigenous from native mollusks (Figure S6a), with size class ">2 to 4 mm" being the most relevant in numerical terms (median relative abundance 0.12 and 0.00, respectively; Figure S7).In contrast, the most abundant size class in native Mediterranean mollusks (>16 to 32 mm, median relative abundance: 0.54) was hardly represented among nonindigenous assemblage components (Figure S7).Concerning feeding types, "macrophyte herbivore", "micrograzer", and "predator" had the strongest associations with nonindigenous trait profiles (Figure S6b), but only the latter-with a median relative abundance of 0.18-was numerically common (Figure S7).Suspension feeding was much less important compared to native mollusks (median relative abundances 0.25 and 0.74, respectively), and chemosymbiotic taxa were entirely lacking, despite being relatively common among Mediterranean species (median relative abundance: 0.08).Finally, while the majority of individuals in both assemblage components lived infaunally, epifaunal lifestyle was better represented in nonindigenous than native mollusks (median relative abundances 0.29 and 0.02, respectively; Figures S6c and S7).Coverage-standardized functional richness of the regional-scale molluscan assemblage from the six continuously sampled National Monitoring stations was characterized by high interannual variability, but showed a significant increase over time (Figure 5a; GLS, p < .01;Table S5).This pattern reflected a pronounced positive trend in this index in the nonindigenous assemblage component (Figure 5d, red line; GLS, p = .001;Table S5), whereas no statistically significant change was detected for native mollusks (Figure 5d, blue line; GLS, p = .11;Table S5).Indeed, the temporal trajectory of NIS functional richness-which, until year 2010, was low but had rapidly increased already by 2011-was well correlated with (i) nonindigenous taxonomic richness (Spearman's rho for differenced time series = 0.83, p < .001)and (ii) the functional richness of the complete molluscan assemblage (Spearman's rho for differenced time series = 0.86, p < .001),suggesting that the diversification of NIS played a major role in shaping functional richness patterns on the shallow Israeli shelf.Functional divergence was high for the complete assemblage (range: 0.69-0.93; Figure 5b), as well as its native (range: 0.60-0.90)and nonindigenous (range: 0.70-0.92)components (Figure 5e), indicating that a large proportion of total abundance was contributed by species with extreme trait profiles.While the value of this index for the complete assemblage rose from c. 0.7 to c. 0.9 over the course of the time series (Figure 5b), there was no statistically significant global trend (GLS with first-order autocorrelation term, p = .14,Table S5).Similarly, no systematic long-term changes were found for both native (GLS, p = .15)and nonindigenous (GLS with first-order autocorrelation term, p = .26)assemblage components (Figure 5e; Table S5).Functional evenness of the complete molluscan assemblage was 0.41 in 2005 and reached a maximum value of 0.60 in 2018, before decreasing to 0.46 at the end of the time series (Figure 5c).Overall, and considering the autocorrelation structure of the data, no significant long-term change in this index was detected (GLS with first-order autocorrelation term, p = .27),a finding that also pertained to native (GLS with first-order autocorrelation term, p = .75)and nonindigenous (GLS, p = .20)mollusks (Table S5;

Figure 5f).
Functional beta diversity of native versus nonindigenous assemblage components sampled in the same year was remarkably high throughout the study period (range: 0.80-1.00; Figure 6, black line), and mostly related to its turnover component (amounting to 65.0%-99.8% of total beta diversity; Figure 6, orange line).This suggests that native and nonindigenous assemblage components occupied largely distinct portions of trait space.However, a significant decline in both indices occurred over time (Jaccard-like dissimilarity index: F I G U R E 5 Temporal development of functional alpha diversity indices for the complete molluscan assemblage (panels a-c; sample coverage (SC) = 0.97) and its native and nonindigenous components (panels d-f, blue and red lines, respectively; SC = 0.93), based on pooled data for the six continuously sampled National Monitoring stations.Solid lines represent mean values calculated across all simulation runs, shaded areas the 95% confidence intervals, and dashed lines in panels (a)-(c) empirically observed index values.Index values for the native assemblage component of year 2010 could not be calculated as coverage-standardized species richness was less than the minimum number of datapoints needed to define a four-dimensional trait space volume (five species) and the minimum spanning tree (three species; see Figure 3b).FRic, functional richness; FDiv, functional divergence; FEve, functional evenness.GLS, p < .01;turnover component: GLS, p < .01;Table S5), indicating increasing overlap (Figure 6).The nestedness-resultant component of functional beta diversity had low values during the entire time series, and did not show a statistically significant change over time (GLS, p = .11;Figure 6, violet line).

| Comparison of taxonomic and functional richness between modern and historical assemblages
The taxonomic richness of the native death assemblage component from Ashqelon was 39 species, its functional richness 0.25, at a sample coverage of 0.96 (Table 1).In contrast, at equal coverage, mean taxonomic and functional richness of the complete living assemblage from the nearby National Monitoring station H28 were only 8 species and 0.001 (95% CI: 7.58 × 10 −9 -0.01), respectively (Table 1).Indeed, in 6% of simulation runs, functional richness values for resampled sets of eight taxa could not even be estimated, as at least five functionally distinct, noncoplanar species are mathematically required to calculate a volume in a four-dimensional trait space (Magneville et al., 2022;Villéger et al., 2008).The live-dead comparison for the northern shelf off Haifa yielded qualitatively similar results: despite the young age of the native death assemblage from station NG10, its taxonomic richness of 50 species (at SC = 0.92) was more than three times greater than that of the temporally integrated living assemblage collected in its close vicinity (National Monitoring site H3, 16 species), and functional richness (0.29) surpassed that of modern mollusks (mean: 0.06, 95% CI: 0.02-0.12)almost fivefold (Table 1).This suggests that historical, purely native assemblages across geographically distant but abiotically similar localities likely covered a much greater range of trait values than the mixed native and nonindigenous assemblages that have inhabited the Israeli shelf for at least the past two decades.

Eastern Mediterranean
Against our hypothesis of a pronounced regional-scale loss of Mediterranean biodiversity in the past two decades, native mollusks from the shallow Israeli shelf did not show statistically significant temporal trends in any of the assessed alpha diversity indices.Instead, we observed a considerable increase in taxonomic and functional F I G U R E 6 Evolution of functional beta diversity of contemporaneous native and nonindigenous molluscan assemblage components (sample coverage = 0.93), based on pooled data for the six continuously sampled (2005-2022) National Monitoring stations.The plot shows the Jaccard-like beta diversity index (black) and its additive turnover (orange) and nestedness-resultant (violet) components.Solid lines represent mean values calculated across all simulation runs, shaded areas the 95% confidence intervals.There was a significant decline in mean values of both the Jaccard-like index and functional turnover over time (GLS, both p < .01),whereas no change was detected in the nestednessresultant component (GLS,p = .11).Index values for year 2010 could not be calculated due to the low number of native species (see caption of Figure 5 for details).Steger, Bakker, et al., 2021;Galil et al., 2021;Guarnieri et al., 2017).The taxonomic diversification went along with a pronounced increase in nonindigenous functional richness (correlation of NIS species richness with NIS functional richness, differenced time series: rho = .83,p < .001),as well as that of the entire regional molluscan assemblage (correlation of NIS species richness with complete assemblage functional richness, differenced time series: rho = .85,p < .001).This indicates that NIS tended to occupy distinct positions in trait space, that is, outside of the convex hull of assemblages from previous years (cf.Azzurro et al., 2014;Toussaint et al., 2018).The presence of trait, and hence niche, differences among NIS on the one hand, and between nonindigenous and native assemblage components on the other-the latter consistently shown also by FCA for the global dataset and individual sampling years (Figure 4; Figure S5a,b)-suggests a limited potential for direct resource competition, corroborating similar findings for fishes and mollusks in the region (Belmaker et al., 2013;Givan et al., 2017Givan et al., , 2018;;Steger et al., 2022;van Rijn et al., 2019).We identified body size, feeding habit, and environmental position as traits most relevant for distinguishing native and nonindigenous assemblage components.The latter had greater shares of small-sized (≤4 mm), predatory, and epifaunal individuals, whereas the intermediate size class ">16 to 32 mm" and suspension feeding were relatively more important among native mollusks.Our results thus match well those of Steger et al. (2022), who found that nonindigenous assemblages across a wide array of habitat conditions had larger proportions of very small-and large-sized, as well as epifaunal individuals, and feeding guild compositions markedly different from their native counterparts.
The long-term temporal stability of trait-based indices incorporating also species´ relative abundances (i.e., functional divergence and evenness)-despite the ongoing taxonomic diversification of NISfurther supported the importance of niche differences for invasion success: throughout the time series, functional divergence of complete and nonindigenous assemblages remained very high, suggesting that successful (i.e., abundant) NIS tended to be those characterized by extreme trait values, occupying the periphery of the assemblage trait space (see also Azzurro et al., 2014); stability in functional evenness indicated that (Indo-Pacific) species and their abundances did not increasingly cluster within particular regions of the trait space, hinting at an effective utilization of available resources while avoiding niche overlap.Due to the still very limited knowledge on Indo-Pacific/Red Sea molluscan assemblages, however, it remains unclear whether the observed functional differentiation primarily reflects biologically mediated trait filtering (Belmaker et al., 2013;Givan et al., 2017;Nawrot et al., 2015).Alternatively, NIS may represent a random subset of an ecologically highly diversified tropical species pool whose trait structure markedly differs from that of the native fauna of Israel, as suggested for body size in Red Sea versus Mediterranean bivalves (Nawrot et al., 2017).
Although native versus nonindigenous functional beta diversity and its turnover component remained high throughout the study period (ranges: 0.80-1.00and 0.56-1.00,respectively), both measures declined over time, suggesting that in more recent years assemblage components shared greater proportions of trait space.This development may reflect three different, not mutually exclusive, processes.
First, the initially low, but progressively increasing taxonomic richness of NIS could have led to a greater probability of species with traits similar to native ones being present for purely stochastic reasons.Second, rapid environmental changes on the Israeli shelf and the rising frequency of extreme events such as heatwaves (Garrabou et al., 2022;Ibrahim et al., 2021) may have selected for particular suites of trait values (i.e., environmental filtering; Keddy, 1992).
Third, increasing trait similarity may indicate a weakening of recipient assemblage biotic resistance, as recently suggested for fish communities, potentially even marking the onset of an "invasional meltdown," a self-reinforcing process in which already established NIS facilitate further invasions (Givan et al., 2017).Although the relative importance of these mechanisms remains unresolved, it seems reasonable to assume that if decreasing functional turnover was largely driven by changes in abiotic and/or biotic selection regimes, we would have observed a convergence of assemblage components' trait profiles over time, as species' relative abundances are expected to respond more rapidly to selective forces than regional-scale occurrence patterns (i.e., native species immigration/NIS establishment, and extirpation; Bates et al., 2014).This, however, was not the case as a pronounced segregation of native versus nonindigenous abundance-weighted trait profiles was not only apparent in the global FCA analysis including samples from all stations and years (Figure 4), but also did not show any signs of decline when yearly datasets were analyzed individually (Figure S5a,b).Still, we cannot disregard the possibility that sublethal effects on native populations, which are not captured by our abundance dataset, have gradually decreased their biotic resistance.Indeed, due to the extreme environmental conditions, today c. 60% of native molluscan species from the Israeli shallow subtidal fail to attain reproductive size, turning populations in the region into demographic sinks (Albano, Steger, Bošnjak, et al., 2021).This ongoing erosion of population fitness may have led to a gradual decline in their resource consumption, thereby increasingly creating opportunities for the establishment of NIS with similar niches (Shea & Chesson, 2002).

| A "novel ecosystem" in demise?
After a steady increase in molluscan abundance during the first decade of the time series, median individual density across the six continuously sampled National Monitoring stations dropped by more than 90% between 2015 and 2016.Similar major declines, though unquantified, have been reported also for other benthic organism groups on the Israeli shelf since 2017 (Lubinevsky et al., 2022).There is mounting evidence that the rapidly rising summer sea surface temperatures in the Eastern Mediterranean (e.g., Ozer et al., 2017;Raitsos et al., 2010) increasingly exceed the physiological tolerance limits of numerous native species, whose populations in the region naturally represent(ed) the low-latitude "warm" trailing edges of their geographic ranges, where organisms are particularly sensitive to anthropogenic climate change (Albano, Steger, Bošnjak, et al., 2021;Rilov, 2016;Rilov et al., 2022;Yeruham et al., 2015Yeruham et al., , 2020)).(Sub-)tropical Indo-Pacific NIS, in contrast, are generally thought to have benefitted from seawater warming, which likely facilitated their establishment and spread to other, more western parts of the Mediterranean (e.g., Bianchi, 2007;Raitsos et al., 2010;Rilov & Galil, 2009).Surprisingly, we found that nonindigenous mollusks on soft substrates suffered a similarly strong density decline as native ones, without any signs of recovery until the year 2022.In line with this observation, shallow-water populations of the formerly abundant Lessepsian tanaidacean crustacean Cristapseudes omercooperi (Larwood, 1954) suddenly disappeared at 5-15 m water depth after 2016, but the species was still present at a site located at 33 m depth (Lubinevsky et al., 2022).Further abrupt, massive declines of abundant Indo-Pacific NIS, such as the invasive mussels Brachidontes pharaonis (P.Fischer, 1870) and Perna perna (Linnaeus, 1758), have been observed in the recent past (Belmaker et al., 2021;Galil et al., 2022;Rilov, 2017), followed by phases of recovery (B.S. Galil, personal observation).Their causes remain poorly understood, but it has been suggested that at least the latter, too, may have been related to extreme summer heat (Galil et al., 2022).Also in with maxima of up to 32°C (Galil et al., 2022;Rilov, 2016), may be prohibitive to the establishment even of several tropical species (Belmaker et al., 2013).NIS that were abundantly present on the Israeli shelf at the time of the observed density decline often established decades ago under a more benign abiotic regime, but environmental conditions have rapidly become more extreme in recent years (Ibrahim et al., 2021;Rilov, 2016;Zamir et al., 2018).Furthermore, several NIS likely have been introduced directly from the Gulf of Suez-the Red Sea sub-basin geographically closest to the Suez Canal-where summer sea surface temperatures today are lower than off the Mediterranean Israeli coast (Albano, Steger, Bošnjak, et al., 2021).We acknowledge that various other stressors, such as disease outbreaks, pollution events, and direct human exploitation (e.g., demersal fisheries), can, individually and synergistically, have detrimental effects on biological assemblages.However, their predicted effects do not fit the patterns observed in our dataset.
Diseases can cause mass mortalities in both native species and NIS (e.g., García-March et al., 2020;Katsanevakis et al., 2019;Zirler et al., 2023), but affect only one or a few taxa at a time, rather than an entire, highly diverse phylum (Fey et al., 2015).Pollution events, while potentially impacting multiple taxa, are limited in their spatial extent, whereas molluscan density declined simultaneously at six stations spanning almost the entire Israeli coastline.Finally, demersal fishing-most notably bottom trawling-in Israeli waters is forbidden at depths as shallow and at distances as close to the shore as our monitoring sites, with tighter restrictions imposed in 2016 (J.Belmaker, personal communication;Edelist et al., 2011, map in Rothschild, 2020), and thus cannot explain the observed rapid regional-scale demise.Furthermore, due to a governmentfunded voluntary "buyout" program aiming at progressively decommissioning trawlers for reasons of environmental protection (Rothschild, 2020), the size of the demersal fishing fleet has gradually decreased during our study period.
Irrespective of underlying drivers, sudden density declines, even without translating into strong changes in assemblage composition, may entail major effects on ecosystem functioning, as ecosystem processes and their rates tend to be largely determined by abundant species (Gaston et al., 2018;Winfree et al., 2015).Thus, although median molluscan density in 2022 was similar to that at the beginning of the time series in 2005, our analysis demonstrated that presentday benthic assemblages on the Israeli shelf are highly volatile, and that even the nonindigenous biota may be affected by serious, unexpected declines, questioning their previously suggested ability to maintain and secure ecosystem functioning in the wake of the massive native species loss (Katsanevakis et al., 2018) and increasingly frequent and intense extreme events (Garrabou et al., 2022;Ibrahim et al., 2021;Zamir et al., 2018).

| Shifted baselines in the Eastern Mediterranean
While the analysis of our unique 18-year-long observational dataset fauna, obscuring the dramatic negative impacts that environmental changes have had on regional biodiversity (e.g., Albano, Steger, Bošnjak, et al., 2021;Rilov, 2016).Furthermore, the effects of the dominated by an invasive red macroalga that revealed considerable differences in their net carbon budgets and habitat provisioning capacity (Peleg et al., 2020).Indeed, natural species losses, which tend to be nonrandom with respect to response traits, may commonly have more severe effects on ecosystem functioning than previously suggested by experimental evidence, which is usually generated by randomly removing species (Wolf et al., 2021).
In light of the sixth mass extinction (Cowie et al., 2022) Living shelled mollusks inhabiting the shallow subtidal sandy strip along the Mediterranean Israeli coast were sampled from 2005 to 2022 in the framework of the National Monitoring program of the Israel Oceanographic and Limnological Research (IOLR).
located c. 750 m northeast of National Monitoring station H28 at 11 m water depth), the other off Haifa (station NG10, situated c. 870 m south of National Monitoring station H3 at 10 m water depth; Steger et al., 2022, orange crosses in Figure The positive trend ended with a sharp decline in median density by 90.6% between 2015 and 2016, followed by a stabilization at considerably lower values during the period 2016-2022 (range: 53-133 inds.m −2 , Figure 2a; Wilcoxon rank sum test comparing yearly medians of the two time periods: W = 73, p < .001).A decomposition of total density into native and nonindigenous assemblage components revealed that their temporal development was well synchronized (Spearman's rho for detrended time series of median densities = 0.69, p < .01);however, the density drop between 2015 and 2016 was more pronounced for NIS (96.1% decline in median, Figure 2c) than for native mollusks (77.4% decline in median, Figure 2b).After this abrupt decrease, median densities of the nonindigenous assemblage component (range: 12-36 inds.m −2 ) were similar to those observed at the beginning of the time series in 2005-2006 (range: 18-39 inds.m −2 , Figure 2c).
; PERMANOVA, stations as strata, pseudo-F = 129.20,p < .001,R 2 = .30),notwithstanding samples in the ordination covered large spatial (most of the Israeli coastline) and temporal (nearly two decades of sampling) scales.Indeed, FCAs covering the 18 sampling years individually (Figure S5a,b) consistently showed an even stronger native versus nonindigenous differentiation, confirming that the small area of overlap in Figure 4 largely reflected interannual temporal variability, rather than the simultaneous presence of functionally similar trait profiles.

F
Boxplots showing the temporal development of log 2 (x + 1) transformed molluscan density at the six continuously sampled National Monitoring stations.(a) Total molluscan density; (b) native molluscan density; (c) nonindigenous molluscan density.The major decline in molluscan density between 2015 and 2016 reflects the corresponding decreases in both native and nonindigenous assemblage components.inds., individuals.

F
I G U R E 3 Temporal development of taxonomic richness of the Israeli molluscan assemblage between 2005 and 2022, based on pooled data for the six continuously sampled National Monitoring stations.(a) Richness trajectory of the complete assemblage (dashed line: observed values; solid line: estimates at sample coverage (SC) = 0.97); (b) trajectories of the native (blue) and nonindigenous (red) assemblage components (SC = 0.93).Shaded areas mark 95% confidence intervals.The richness of the complete assemblage and its nonindigenous component significantly increased over time (GLS, both p < .01),whereas no significant change was found for the native assemblage component (p = .32;see Table S5 for model details).F I G U R E 4 Fuzzy correspondence analysis ordination plot of trait profiles of native (blue) and nonindigenous (nonindig.;red) molluscan assemblage components from all 13 National Monitoring stations, sampled between 2005 and 2022.Native and nonindigenous trait profiles are strongly segregated (PERMANOVA, pseudo-F = 129.20,p < .001,R 2 = .30)across sampling stations and years.
2015, the Eastern Mediterranean Sea was affected by a massive heatwave (SenGupta et al., 2020): the core event lasted from 17 July 2015 to 21 August 2015 and was the most severe recorded for the region throughout the authors' 36-year study period(1982-  2017)  regarding intensity, duration and spatial extent, highlighting its considerable potential to impact benthic populations in very shallow water.Sampling for the National Monitoring program in that year was conducted in early August, when the heatwave had already started but not yet ended, suggesting that any major and long-lasting biotic effects might be mostly visible only in the data of year 2016, matching well the empirically observed patterns (Figure2).Although our reasoning remains hypothetical at this point, the potential sensitivity of NIS to very high temperatures was indicated also by a macroecological study that evaluated the nativerange environmental conditions of Indo-Pacific fishes, suggesting that summer heat peaks in the present-day Eastern Mediterranean, provided valuable insights into the development of the soft-bottom biota on the Israeli shelf, several independent lines of evidence clearly demonstrate that assemblages at the start of the time series represented an already significantly shifted baseline.First and most obviously, approximately half of the taxonomic richness (at SC = 0.93), and more than a third of the individuals comprising the molluscan assemblage in year 2005 were nonindigenous.Lessepsian NIS have entered the Mediterranean Sea since shortly after the opening of the Suez Canal in 1869, and their number has steadily increased over time.By 2005, more than 300 introduced species had been recorded just from the Israeli coast-the vast majority of them of Indo-Pacific origin, and one third being mollusks(Galil et al., 2021).Second, also the | 13 of 17STEGER et al.native component of our sandy-bottom assemblage has undergone massive changes already prior to the start of the National Monitoring program.By compiling information on selected, formerly abundant shallow-water species from a range of historical sources,Rilov (2016) revealed major climate-driven multi-species population collapses on the Israeli shore that were much advanced already by the turn of the millennium.These findings were recently corroborated and extended by an analysis of radiometrically dated molluscan death assemblages that demonstrated massive species losses across different shallowsubtidal habitats in the region, with only 5%-12% of native mollusks originally present still recorded alive at the time of sampling in 2016-2018(Albano, Steger, Bošnjak, et al., 2021).Third and last, the novel mixed native and nonindigenous living assemblages from National Monitoring sites off southern and northern Israel proved taxonomically and functionally much impoverished compared to nearby native death assemblages, despite the former being "enriched" by Indo-Pacific NIS.Both stations, notwithstanding their geographic distance and differences in death assemblage age, represented homogenous habitat conditions(Lubinevsky et al., 2019) and yielded qualitatively similar results, reaffirming they most likely reflect consistent regionalscale ecological patterns, rather than particularities of the local geohistorical records or effects of unrecorded past point-source disturbances.Importantly, bothAlbano, Steger, Bošnjak, et al. (2021)    and our study recovered pronounced signals of a historically much greater (native) biodiversity even from young death assemblages with median ages of just a few decades.This confirms that the majority of taxonomic and functional losses had occurred rather recently, though -according to the National Monitoring data and Rilov (2016) -largely happened before year 2005.Indeed, we noted only a few potential regional-scale population collapses when analyzing our time series, including, for example, the widespread and common Mediterranean bivalves Macomangulus tenuis (da Costa, 1778) (present at the six continuously sampled stations until 2011) and Acanthocardia tuberculata(Linnaeus, 1758(Linnaeus,  ) (present until 2015)).Assuming a constant extirpation risk per species and unit time, assemblage-level loss rates are expected to decrease with increasing taxonomic impoverishment, potentially explaining the lack of significant negative trends in native species richness in the National Monitoring time series.The evidence presented above highlights that several results of our time series analysis-despite covering a significant period of almost two decades-may easily lead to biased conclusions about the long-term development and current state of the Eastern Mediterranean benthos, unless shifted baselines are explicitly accounted for.Specifically, the overall increase in taxonomic and functional richness related to the diversification of NIS masked that the present-day mixed native/nonindigenous biota on the Israeli shelf has remained much impoverished compared to the historical native Lessepsian invasion on assemblage functional composition are grossly underestimated if only data collected since 2005 are taken into consideration: due to the consistently different trait profiles of nonindigenous assemblage components, present-day benthic assemblages almost certainly differ functionally from their historical, purely native predecessors, though the exact nature of changes remains obscure due to the lack of historical data (see alsoSteger et al., 2022).However, during the study period, no apparent directional changes in assemblage trait profiles were noted, despite an overall increasing share of Lessepsian mollusks.This was likely related to the great variability in the proportions of native versus nonindigenous individuals at the different stations in most years, itself probably indicative of the ongoing taxonomic reorganization of the Israeli shallow-water benthic assemblages.Further evidence for invasion-related functional changes in Eastern Mediterranean assemblages was recently provided also by in situ comparisons between relic patches of native brown algal (Cystoseira) forests, algal turfs (due to overgrazing by Lessepsian rabbit fishes), and habitats and the ongoing erosion of natural biogeographic patterns by humanmediated species introductions(Bernardo-Madrid et al., 2019), the need to closely follow the trajectories of ecosystem reorganization is greater than ever before.However, even in novel ecosystems shaped by massive biological invasions and unprecedented losses of native biodiversity, shifted baselines may easily lead to wrong conclusions about a systems' state and fate if not explicitly quantified.Due to the often long history of anthropogenic alterations, these biases can affect analyses even of systems where multi-decadal high-quality observational data are available.Our study on Israeli molluscan assemblages provides an example of how the impediment of limited temporal coverage may be alleviated by integrating information from the young geohistorical record (seeKidwell & Tomašových, 2013) into the analysis of ecological time series.Considering that (i) changes in several ecosystems may have primarily involved a reshuffling of abundances of native species, rather than conspicuous invasions by tropical NIS as on the Israeli shelf and (ii) the coverage by (old) scientific literature or other observational data sources may be limited, shifted baselines will often be much more difficult to detect than in the Eastern Mediterranean, underscoring the need to access diverse sources of information to reconstruct past ecosystem states.Conceptualization; data curation; formal analysis; investigation; methodology; software; validation; visualization; writing -original draft; writing -review and editing.Cesare Bogi: Data curation; investigation; writing -review and editing.Hadas Lubinevsky: Data curation; investigation; writing -review and editing.Bella S. Galil: Data curation; investigation; writing -review and editing.Martin Zuschin: Funding acquisition; resources; writing -review and editing.Paolo G. Albano: Conceptualization; data curation; funding acquisition; investigation; validation; writing -original draft; writingreview and editing.
of the nonindigenous assemblage component, which drove corresponding positive trends in index values of the entire molluscan assemblage.Expected NIS richness rose by 0.55 species per year throughout the study period, consistent with ongoing introductions and the establishment of new Indo-Pacific taxa in the region (e.g.,Albano, Note: Modern assemblages, comprised of native and nonindigenous species, were standardized to the sample coverage (SC) of corresponding native death assemblages to enable unbiased comparisons.For modern assemblages, estimates of species richness and mean functional richness are provided together with their 95% confidence intervals (CIs); for death assemblages, empirically observed values are presented.Abbreviations: FRic, functional richness; No. of inds., number of individuals; Tax.richness, taxonomic/species richness.diversity