Non- indigenous molluscs in the Eastern Mediterranean have distinct traits and cannot replace historic ecosystem functioning

Aim: A large body of ecological theory predicts that non- indigenous species (NIS) are successful invaders if their niches overlap little with native taxa. Native– non-indigenous trait dissimilarity, however, may also be observed if NIS have outcompeted ecologically similar native species. Discriminating these scenarios is essential for assessing invasion impacts but requires baseline assemblage data that are fre-quently unavailable. We overcome this impediment


| INTRODUC TI ON
In today's globalized world, introductions of non-indigenous species (NIS) occur at unprecedented rates (Seebens et al., 2017), with often profound impacts on native biodiversity and ecosystem functioning (e.g., Pyšek et al., 2020;Strayer, 2012). Assessing the potential for such impacts requires knowledge on how the traits of NIS relate to those of native species, as species' functional attributes are important determinants of biotic interactions and ecosystem processes (e.g., Engelhardt et al., 2009;Hooper et al., 2005;McGill et al., 2006;Pearson et al., 2012).
To successfully establish and increase in abundance, NIS do not only require traits that allow them to cope with the local abiotic conditions, but also to acquire sufficient resources to sustain growth and reproduction (e.g., Byers, 2000;Gallien & Carboni, 2017). While some NIS may achieve the latter by being superior competitors compared to functionally similar native species (e.g., Byers, 2000), a large and influential body of ecological theory suggests that niche differences to native species are key determinants of invasion success (see Catford et al., 2009 for a review). Such differences may enable NIS to effectively exploit available resource opportunities (sensu Shea & Chesson, 2002) and alleviate direct competition with resident species, as predicted by concepts such as limiting similarity (MacArthur & Levins, 1967), Darwin's naturalization hypothesis (Darwin, 1859), the 'empty niche' hypothesis (Hierro et al., 2005), and related theoretical frameworks (see Catford et al., 2009 for a review). While aiming at mechanistic explanations of invasions success, these hypotheses also yield two major predictions regarding the impact of biological invasions: first, due to complementary resource use, NIS are unlikely to competitively displace native species; second, biological invasions should cause major functional shifts at the ecosystem scale because successful (i.e., abundant) NIS will substantially alter the trait structure of recipient assemblages (Pearson et al., 2012).
Studies exploring the functional dissimilarity between NIS and native species -either directly via traits or indirectly via phylogenetic distances -have produced contradictory results (e.g., Diez et al., 2008). However, this may largely reflect differences in the conceptual frameworks and spatial scales involved (Thuiller et al., 2010).
Indeed, considering that resource partitioning and competition are most relevant at small spatial scales and advanced stages of the invasion process (e.g., Schaefer et al., 2011), several recent empirical works found that abundant or invasive NIS often tend to differ either phylogenetically, ecologically and/or morphologically from native species (Azzurro et al., 2014;Diez et al., 2008;Divíšek et al., 2018;Pearson et al., 2012).
Such findings, however, can be confounded by the lack of a historical perspective on community dynamics in systems with long invasion histories. In such systems, any differences in the trait composition of present-day native and non-indigenous components of local assemblages could represent the outcome of two contrasting scenarios. First, abundant NIS may indeed be those whose niches overlap little with native species; hence, trait dissimilarity was present since the onset of their introduction. Second, present-day trait patterns may represent 'shifted baselines' (Pauly, 1995): successful NIS have outcompeted native species with similar traits in recipient assemblages, resulting in the trait segregation observed today. In this latter scenario, non-indigenous assemblage components should share high trait similarity with historical native assemblages.
The limited temporal coverage of observational data for many ecosystems relative to their invasion history, particularly in the marine realm (Ojaveer et al., 2018), therefore poses a major limitation to the interpretation of empirically observed trait patterns (Thuiller et al., 2010) and unbiased assessments of invasion impacts. Death assemblages -the accumulations of taxonomically identifiable organism remains encountered in a landscape or seafloor -are faithful natural archives of community composition that can help overcome this impediment (see Kidwell & Tomasovych, 2013 for a review). Due to the durability of skeletal elements such as molluscan shells, foraminifera tests and fish otoliths, and generally low sedimentation rates in open shelf environments, remains of multiple generations accumulate in death assemblages over time-scales of decades to millennia (Kidwell & Tomasovych, 2013). This extensive temporal mixing ('time-averaging') renders death assemblages compositionally inert to any recent directional change in living communities, and underlies their ability to record long-term ecological baselines and pre-impact assemblage composition. Capitalizing on these properties, death assemblages have successfully been exploited to reveal unrecorded shifts in taxonomic composition in habitats with long histories of anthropogenic alteration (e.g., Tomašových & Kidwell, 2017), and can onset of the invasion, suggesting that competition was unlikely the primary driver of the regional-scale native biodiversity loss. Our findings, however, also imply that NIS cannot functionally compensate for native species disappearance. Instead, the transition towards increasingly NIS-dominated assemblages has profoundly altered ecosystem functioning, with unknown consequences.

K E Y W O R D S
biological invasions, competition, death assemblages, ecosystem functioning, historical baselines, Lessepsian invasion, marine molluscs, non-indigenous species, traits also provide invaluable information on the trait composition of past assemblages (Miller et al., 2014).
In this study, we combine present-day and geohistorical data from living and death assemblages, respectively, to assess the impacts of the so-called Lessepsian invasion (Por, 1978) -among the largest and longest-lasting biological invasions (Rilov & Galil, 2009) -on the native biota and ecosystem functioning in the Eastern Mediterranean Sea: following the opening of the Suez Canal in 1869, hundreds of tropical Red Sea species established populations in the basin in the course of this largely unidirectional process, and have become dominant components of local communities (e.g., Edelist et al., 2012;Galil et al., 2021;Rilov et al., 2018). Concurrent with this development, massive regional-scale population collapses of formerly common native Mediterranean species have been reported and suspected to primarily represent the outcome of competitive interactions with NIS (e.g., Edelist et al., 2012;Galil, 2007).
We here test the hypothesis that native and non-indigenous components of benthic assemblages from the shallow Israeli Mediterranean shelf have differed in their trait composition since the onset of the Lessepsian invasion. We focus on molluscs, the largest and taxonomically and functionally most diverse group of Lessepsian invaders, whose shells provide rich death assemblages that enable reconstructing baselines. This extended temporal perspective, combined with a unique multi-assemblage dataset covering a wide range of shallow subtidal soft and hard substrates, allowed us to gain important insights into general patterns of native versus nonindigenous trait structure, functional diversity, and the potential of the Lessepsian invasion to impact Eastern Mediterranean ecosystem functioning.

| Field sampling and sample treatment
We analysed molluscan assemblages from 72 samples collected at 10 stations along the c. 200-km-long Israeli Mediterranean shelf between 2016 and 2018 (see Supporting Information Appendix S1, Figure S1.1 for a map). We collected 48 van Veen grab samples from subtidal soft substrates (10-41 m depth) and twenty-four 1-m 2 airlift suction samples from subtidal hard substrates (11-28 m depth).
Our analysis is based on the resulting 45,414 living individuals and 12,041 shells representing 360 species (274 native taxa and 86 NIS). Sampling was conducted in both spring and autumn to capture seasonal variation. Detailed information on sampling stations and samples is available from Supporting Information Table S1.1.
Bulk samples were sieved on a 0.5-mm mesh, the retained material fixed in ethanol and subsequently picked for molluscan individuals in the lab. For the analysis of death assemblages, quantitative splits of air-dried sample residues were picked until approximately 1,000 shell elements per station had been obtained. Only shell fragments constituting at least half of the shell, with either the apex or aperture (in gastropods) or the hinge (in bivalves) preserved, and retaining enough morphological characters to enable reliable taxonomic identification were considered. To obtain abundances comparable to living individuals, the abundance of skeletal elements was corrected for the number of elements present in living individuals, that is, dividing the raw counts of loose bivalve valves by 2, and those of polyplacophoran plates by 8 (Kowalewski et al., 2003). Shells of pelagic gastropods and non-marine taxa transported into the sampling sites were excluded from the dataset (Supporting Information Appendix S2).
To enable a rigorous temporal interpretation of the ecological information extracted from time-averaged death assemblages, we quantified their shell age distributions by radiocarbon dating (see Albano et al., 2021). At each station, we dated 10 to 15 valves of common native bivalve species that were also recorded alive in our samples (see Supporting Information Table S1.2 for further details and a list of dated species), using accelerator mass spectrometry (AMS) of powdered carbonate targets (Bright et al., 2021;Bush et al., 2013). Ages are reported in calendar years before the year of sample collection. Detailed descriptions of sample treatment, radiocarbon analysis, calibration of ages, and individual shell age data are available from Albano et al. (2021, supplementary material).

| Trait dataset
To functionally characterize the molluscan species in our samples, we collected information on five ecologically important traitsmaximum adult body size, feeding habit, environmental position, substrate affinity and host association -using published literature, online databases and/or observations on our specimens (see Supporting Information Appendix S3 for a list of references). We chose these traits because (a) they capture several important aspects of molluscan ecology, (b) reliable information is available for the great majority of species covered in this study, and (c) they have been successfully used in studies focusing on invasion success in Lessepsian bivalves (Nawrot et al., 2015) and on native versus nonindigenous niche overlap in Eastern Mediterranean rocky intertidal molluscs . Details on the ecological relevance of the selected traits and the number of scored modalities (i.e., the different categories of a trait) are provided in Supporting Information Table S1.3.
As many species exhibit affinities to more than one modality per trait, we used a fuzzy-coding procedure (Chevenet et al., 1994) to generate the species × traits matrix. This approach allows scoring of trait modality affinity in a flexible, semi-quantitative way, providing a more realistic representation of species' functional roles in ecosystems. Following recommendations by Degen et al. (2018), we adopted the widely used numeric scoring system ranging from 0 to 3, where 0 corresponds to 'no affinity' to a given modality, 1 to 'low affinity', 2 to 'high affinity' and 3 to 'exclusive affinity'. If two or more modalities are expressed to a similar degree, each was scored with '2' by convention. In the rare case that information on a particular trait was not available for a species (concerning seven species, or 2% of taxa), all its modalities were scored with '0', so the profile would not bias subsequent analyses (Chevenet et al., 1994). After collecting trait information, raw modality scores of each trait were converted to proportions adding up to unity (e.g., Oug et al., 2012); in this way, all species and traits were weighted equally in the analyses, irrespective of the number of modalities per trait. This standardized trait matrix is provided in Supporting Information Appendix S4.

| Trait composition of native versus nonindigenous assemblage components
Patterns of multivariate trait composition for native and nonindigenous assemblage components within stations were assessed by fuzzy correspondence analysis (FCA; Chevenet et al., 1994), an extension of multiple correspondence analysis designed for fuzzycoded data. To calculate assemblage component trait profiles, we weighted standardized modality scores of species (modalities × species table) by their relative abundances (species relative abundance × assemblage component table) using matrix multiplication (see Oug et al., 2012;Steger et al., 2021); this takes into account that ecosystem functioning is mainly driven by dominant species (e.g., Gaston et al., 2018;Winfree et al., 2015). Assemblage component replicates with < 5 individuals were excluded from the analysis, including all spring-season non-indigenous components of station SG10, and station NG10 (no non-indigenous components of sufficient sample size were available). We further excluded trait modalities not represented at a given station (i.e., zero-only columns; see Supporting Information Table S1.4). FCA was performed using the 'ade4' package (v. 1.7-13; Dray & Dufour, 2007) in the statistical programming environment R (v. 3.5.2; R Core Team, 2018). Correlation ratios (Chevenet et al., 1994) were used to assess the relevance of the different traits for separating native and non-indigenous assemblage components. Differences in FCA scores of present-day native and non-indigenous assemblage components on the first three axes were tested by permutational multivariate analysis of variance (PERMANOVA; Anderson, 2001), using Euclidean distances and 9,999 permutations. For death assemblage components, the significance of the native versus non-indigenous separation was assessed by running PERMANOVA on ordination scores calculated for sets of 30 analytical replicates per station, obtained by a resampling procedure (see Supporting Information Appendix S1, section 1.5 for details). Finally, we used PERMANOVA (Bray-Curtis distances, 9,999 permutations) to test for native versus non-indigenous differences in the abundance-weighted modality composition of individual traits; to control for the influence of local environmental conditions on modality composition, sampling stations were defined as strata in this across-sites analysis. PERMANOVA was conducted using the 'vegan' R package (v. 2.5-4; Oksanen et al., 2019).

| Functional diversity
We assessed the functional diversity of native and non-indigenous assemblage components at the different stations to gain further insights into their trait structure and diversity. Functional diversity was measured based on the arrangement of species and their abundances in functional space, a multidimensional coordinate system in which species are positioned according to their trait (dis)similarity (Villéger et al., 2008). To build this space, we first grouped the 360 partially redundant species trait profiles into 151 unique functional entities (i.e., distinct trait combinations, see e.g., Villéger et al., 2011). We then calculated a functional distance matrix for these entities, using the generalized Gower dissimilarity measure (Pavoine et al., 2009), which can handle mixed variable types, including fuzzy-coded variables (Pavoine et al., 2009). The distance matrix was subsequently converted into Euclidean functional space using principal coordinates analysis (PCoA; see Legendre & Legendre, 2012). Functional spaces with up to 10 dimensions were built and their quality assessed using an updated and modified version of the 'quality_funct_space' R function (Maire et al., 2015), available from Villéger (2017). The selected 4-dimensional functional space had a mean squared deviation of 0.0147 between the Gower and Euclidean distances -confirming that it faithfully represents the original functional relationships (see Maire et al., 2015) -and was the best compromise between retaining information and enabling the calculation of functional diversity indices (see below) even for the most species-poor assemblage components. To test the sensitivity of our analysis to the choice of traits, we calculated Gower distances based on all different combinations of four out of five traits (cf. Toussaint et al., 2016). The correlation between these matrices and the distance matrix based on all five traits was assessed by Mantel tests and found to be robust against trait omission (Mantel r > .87, all p = 10 -4 ; Supporting Information Table S1.5).
We assessed the functional alpha diversity of present-day assemblage components using three complementary indices: functional richness, functional divergence and functional evenness (Villéger et al., 2008). As the latter two incorporate information on species relative abundances, these indices were calculated after assigning to all species the coordinates of their respective functional entities. Functional richness measures the amount of functional space occupied by an assemblage, that is, the volume of the minimum convex hull including all constituent species, and thus reflects the range of trait values present (Villéger et al., 2008). This index does not consider species abundances and has no theoretical upper limit.
Functional divergence measures how species abundances are distributed in assemblage trait space, and can be regarded as a proxy of functional niche differentiation (Mouchet et al., 2010); it represents the deviation of species abundances from the mean distance to the assemblage centroid (Mouchet et al., 2010;Villéger et al., 2008).
Functional divergence is independent of species richness and takes high values if species with extreme trait profiles have high relative abundances in the assemblage (Villéger et al., 2008). Functional evenness measures the regularity of species' distribution and their relative abundances in assemblage functional space, based on the minimum-spanning tree connecting all constituent species. This index therefore takes low values if abundances are clustered in certain parts of the assemblage functional space (Mouchet et al., 2010).
Finally, to explore the degree of functional differentiation between corresponding native and non-indigenous assemblage components independent of species' relative abundances, we calculated their functional beta diversity, using the Jaccard-like dissimilarity index proposed by Villéger et al. (2013). This index, ranging between 0 and 1, can be decomposed into additive turnover and nestedness- Second, we conducted paired-samples Wilcoxon signed rank tests to assess any significant differences in functional diversity index values between native and non-indigenous assemblage components; due to the limited number of hard substrate stations (n = 4), signed rank tests were not conducted for this habitat. Finally, we also tested for habitat-specific differences in the values of alpha functional diversity ratios, as well as of functional beta diversity indices, by comparing soft and hard substrate stations with Wilcoxon rank sum tests.

| Time range of historical baselines
Death assemblages from the different stations had median ages ranging from 24 years (at NG10) to 1,461 years (at SG30, Table 1; see also Albano et al. (2021)). The oldest median ages were observed at soft substrate stations off southern Israel (SG10, SG20, SG30 and SG40), where death assemblages contained at least 53% of pre-Lessepsian (i.e., pre-AD 1869) shells (Table 1). At stations SG20 and SG30, this share was 67 and 93%, respectively, suggesting that these death assemblages constitute good pre-invasion baselines (see also Supporting Information Figure S1.2). In contrast, death assemblages from soft substrates off northern Israel (stations NG10 and NG30), and those collected on hard substrates, were much younger (Supporting Information Figure S1.2), and had decadal-scale median ages (Table 1). Overall, there was a significant correlation between death assemblage median age and the proportion of pre-Lessepsian shells (Spearman's rho = .70, p = .02). TA B L E 1 Median ages (in calendar years before sample collection) of native death assemblages from the Israeli Mediterranean shelf (from Albano et al., 2021) and percentages of shells pre-dating the opening of the Suez Canal in AD 1869 (pre-Lessepsian shells), determined by radiocarbon dating

| Trait composition of assemblages
Native and non-indigenous components of present-day (i.e., living) assemblages consistently differed in their multivariate trait composition, forming well-defined clusters in ordination space (FCA: Figure 1; PERMANOVA: all p < .05, R 2 range: .32-.89, Table 2). This segregation was independent of occasionally pronounced seasonal variation in trait composition within groups, and mostly occurred along FCA axis 1 (i.e., the axis explaining the greatest amount of variance in the data), except for stations SG40 and Ashqelon −25 m, which separated along axis 2 (Figure 1). At SG40 (Figure 1e) Table 2) trait composition, whereas at Ashqelon −25 m (Figure 1g), this axis reflected differences between corresponding present-day and historical (i.e., death) assemblage components.
Similar to present-day assemblages, a strong trait segregation was apparent between native and non-indigenous components of almost all historical assemblages (Figure 1; PERMANOVA: all p = 10 -4 , R 2 range: .97-1.00, Supporting Information Table S1.6), irrespective of the amount of time that they captured (correlation analysis between the proportion of variance explained by the factor 'native versus non-indigenous' (PERMANOVA R 2 ) and death assemblage median age: Spearman's rho = .35, p = .35). In contrast to the clear native versus non-indigenous separation, corresponding present-day and historical assemblage components often functionally overlapped or were separated along a direction roughly perpendicular to the segregation of present-day native and non-indigenous assemblage components (Figure 1).
Traits that repeatedly had the highest correlation ratio for axis 1 were 'maximum adult body size' (SI; 56% of the nine stations) and 'environmental position' (EP; 22% of stations), whereas 'feeding habit' (FH) had the second highest correlation ratio in 89% of stations (Supporting Information Table S1.7). Considering axis 2, 'maximum adult body size' and 'feeding habit' had the highest correlation ratio in 56 and 33% of stations, respectively (Supporting Information Table   S1.7). A comparison of the trait modality distribution between presentday native and non-indigenous assemblage components revealed significant differences for all five traits (PERMANOVA, all p <. 001, R 2 range: .03-.26, Supporting Information Table S1.8), but complex station-specific patterns (Supporting Information Figures S1.3-S1.7).
However, non-indigenous components from most stations had a greater proportion of both very small (size classes > 1 to 2 mm and > 2 to 4 mm) and rather large-sized taxa (> 32 to 64 mm) compared to their native counterparts, whereas among the latter, medium-sized taxa (> 8 to 16 mm, and to a lesser extent > 16 to 32 mm) were usually better represented than in non-indigenous components (Supporting Information Figure S1.3). Additionally, non-indigenous components TA B L E 2 Results of permutational multivariate analysis of variance (PERMANOVA) conducted on fuzzy correspondence analysis (FCA) axis 1-3 scores of present-day molluscan assemblage components from the Israeli Mediterranean shelf, using the factor 'native versus non-indigenous' F I G U R E 2 Functional diversity index ratios (native versus non-indigenous components, log 10 -scale) for present-day molluscan assemblages from soft and hard substrates on the Israeli Mediterranean shelf. The percentage of ratios > 0 (i.e., greater index values of the native component) is provided below each boxplot. Neither assemblage component consistently had a greater functional richness (FRic, panel a), whereas native components on hard substrates always had a lower functional divergence (FDiv, panel b, all ratios < 0), but greater functional evenness (FEve, panel c, all ratios > 0) than their non-indigenous counterparts. For all indices, absolute values of index ratios did not differ significantly between soft and hard substrates (Wilcoxon tests, W ≤ 15, p-values as indicated on panels). Note the discontinuous y axis in (a); the dashed line marks index value equality of assemblage components were often characterized by a greater prevalence of surface deposit feeders and, particularly on hard substrates, fewer suspension feeders (Supporting Information Figure S1.4). Non-indigenous components from hard substrates were further richer in browsing carnivores/ectoparasites and had a higher prevalence of host-associated taxa. On soft substrates, non-indigenous components had a greater share of epifaunal but fewer host-associated taxa than native components (Supporting Information Figures S1.5 and S1.7, respectively). Trait modalities exclusively present in native assemblage components were the feeding types 'subsurface deposit feeder' and 'scavenger', but they were extremely rare (Supporting Information Figure S1.4). Similarly, chemosymbiotic molluscs were mostly restricted to native assemblage components, and usually found in very low proportions.

| Native versus non-indigenous functional diversity
Log  Information Table   S1.9). Native components of hard substrate assemblages always had a lower functional divergence than corresponding non-indigenous components (ratio range: −0.19 to −0.02, median: −0.09), that is, a smaller proportion of native abundance was contributed by species with extreme trait values. In contrast, no differences were found for soft substrate assemblages (ratio range: −0.18 to 0.09, median: −0.06; signed rank test, V = 7, p = .56; Figure 2b), similar to functional richness.
Native components had a greater functional evenness in all studied hard substrate assemblages (ratio range: 0.05 to 0.14, median: 0.12;

| Interpreting present-day trait patterns with a historical perspective
Native and non-indigenous trait profiles were consistently segregated in the ordinations, indicating that abundant NIS have traits F I G U R E 3 Functional beta diversity of native versus non-indigenous components of present-day molluscan assemblages from the Israeli Mediterranean shelf (a), and its decomposition into turnover (b) and nestedness-resultant (c) components. In both substrate types, beta diversity was mainly related to functional turnover. Soft substrate assemblages had significantly greater index values for total beta diversity and functional turnover than those from hard substrates (Wilcoxon tests, both W = 24, p = .01), whereas no such difference was found for the nestedness-resultant component (W = 11, p = .91) dissimilar from dominant native taxa. This result applies not only to present-day assemblages, but also to historical (i.e., death) ones, suggesting that Lessepsian invaders have occupied 'empty' niches since the onset of the invasion. Indeed, if the present-day pattern were the outcome of functional segregation over time, then we should have observed overlap between non-indigenous assemblage components and native death assemblages. This, however, was not the case: at almost all stations, historical trait profiles either overlapped with their corresponding present-day ones ( Figure 1a,c and e), or showed an offset along roughly the same direction (Figure 1b,d,f and g). This offset therefore reflects factors similarly affecting the trait composition of native and nonindigenous assemblage components, such as taphonomic bias (e.g., the selective loss of fragile shells or shell transport away from/into the sampling site), recent changes in local environmental conditions that modified both assemblage components, or a combination of both (Kidwell, 2013). In contrast, the accumulation of rare species, typical of time-averaged death assemblages (Kidwell, 2013), is unlikely to notably decrease functional live-dead fidelity as (a) abundance-weighting emphasizes the traits of dominant species, reducing the influence of rare taxa, even if they had distinct functional properties; (b) the number of scored trait modalities is limited (23 in this study) and therefore less sensitive to the effects of timeaveraging than taxonomic richness; and (c) the offset occurs in both highly (e.g., station SG30) and minimally time-averaged (e.g., hard substrate stations) death assemblages (Supporting Information Figure S1.2). This reasoning is supported by findings of high fidelity in feeding guild structure even for death assemblages whose taxonomic richness was greatly inflated compared to corresponding living assemblages (García-Ramos et al., 2016).
Importantly, the clear separation of native and non-indigenous trait profiles not only occurred across sampling sites with very different habitat conditions, but was also independent of the age of local death assemblages: for example, the ordination pattern obtained for station SG30 (Figure 1d), where the death assemblage was almost entirely (93%) pre-Lessepsian -and thus constitutes an ideal baseline (Table 1, Supporting Information Figure S1.2) -mirrored those for most hard substrate sites (e.g., Figure 1f,g and i) whose death assemblages generally encompassed only the last few decades. Notwithstanding their young age, such death assemblages still retain a very strong signature of diverse native assemblages because they mostly pre-date the collapse of native species and the onset of the massive increase in NIS prevalence that took place only during the most recent decades Edelist et al., 2012;Galil et al., 2021;Rilov, 2016). Thus, by adopting an extended temporal perspective utilizing death assemblages, we demonstrated that the Lessepsian invasion has followed a pattern consistent with concepts like limiting similarity (MacArthur & Levins, 1967), Darwin's naturalization hypothesis (Darwin, 1859), or the 'empty niche' hypothesis (Hierro et al., 2005;see Catford et al., 2009 for a review), and that present-day trait differences indeed do not represent 'shifted baselines' (Pauly, 1995).

| The role of NIS in the collapse of native species
The consistent differences in trait composition, particularly the lack of overlap between non-indigenous trait profiles and historical native ones, suggest that native assemblages on the shallow Israeli Mediterranean shelf have not been functionally displaced by NIS.
Our findings for molluscs corroborate studies on modern fish assemblages that found a low potential for direct competition between native and non-indigenous species in the Eastern Mediterranean Sea (Arndt et al., 2018;Buba & Belmaker, 2019;Givan et al., 2018). We here extend observations to another ecologically important phylum, a diverse range of shallow subtidal soft and hard substrates and, most importantly, across a long temporal perspective. Although none of the above-mentioned fish studies had data on native assemblages before the onset of the Lessepsian invasion, they analysed temporal trends using trawl datasets from two distinct time periods (1990-1994 versus 2008-2011) that saw a major turnover in the local fish communities: a marked increase in both the number and relative abundance of non-indigenous fishes, and a strong decline of several native species (Edelist et al., 2012). Together, fishes and molluscs account for 54% of the 452 multi-cellular NIS recorded in Israel to date (Galil et al., 2021), suggesting that these findings might be of general validity, although a broader coverage of Lessepsian phyla is required to draw final conclusions.
We emphasize that our assemblage-scale results must not be interpreted as dismissing the potential for significant negative interactions. For example, NIS can affect native assemblages through various direct and indirect interactions other than competition, such as predation (e.g., Green et al., 2012) or habitat modification by ecosystem engineers (e.g., Crooks, 2002). Nevertheless, it seems improbable that these factors can explain the massive (88-95% in the shallow subtidal) loss of native molluscan diversity observed at all stations covered by this study (see Albano et al., 2021). While nonindigenous predators can severely decimate vulnerable native species (e.g., Chiba & Sato, 2013;Hadfield et al., 1993), it is unlikely that they could drive to (near) extirpation most representatives of a taxonomically, ecologically and morphologically highly diversified group such as Mediterranean molluscs, in both soft and hard substrates.
For example, the formerly abundant large-sized shallow-water predatory whelk Stramonita haemastoma (Linnaeus, 1767) has undergone a major collapse along the Israeli shore in recent decades despite the absence of known non-indigenous predators or potential competitors (Rilov, 2016). This is surprising, considering that this species was shown to prefer the abundant Lessepsian mussel Brachidontes pharaonis (P. Fischer, 1870) as prey over smaller native mussels and barnacles, and therefore should even have energetically benefited from the invasion (Rilov, 2016;Rilov et al., 2002 (Vergés et al., 2014;Yeruham et al., 2020), the shrub-like macroalga Galaxaura rugosa (J. Ellis & Solander) J. V. Lamouroux, 1816 (Peleg et al., 2020), and the bivalves Chama pacifica Broderip, 1835 and B. pharaonis (Nawrot et al., 2015) are restricted to particular habitat types and depth intervals, whereas native biodiversity loss invariably occurred across all studied habitats . Furthermore, several native components of soft-bottom molluscan assemblages apparently have not changed functionally over time (Figure 1a,c and e) as would have been expected under modified local habitat conditions (see e.g., Kidwell, 2009). Together, these findings imply that regionalscale changes in abiotic conditions that broadly affect species across stations, such as warming seawater temperatures (Ozer et al., 2017), likely play a major role as driver of the massive native biodiversity loss in shallow subtidal habitats along the Israeli coast (see Albano et al., 2021;Givan et al., 2018;Rilov, 2016). Warming acts on such a broad spatial scale that as conditions become unsuitable for many native species, those with more restricted ranges may eventually go globally extinct as a small geographical range is among the best predictors of extinction risk (Chichorro et al., 2019).

| A 'novel ecosystem' in the Eastern Mediterranean
Given the distinct trait profiles of native and non-indigenous assemblage components, the faunal transition towards increasingly NIS-dominated systems (e.g., Rilov, 2016;Rilov et al., 2018) has considerably altered key functional properties of shallow-water molluscan assemblages on the Israeli Mediterranean shelf, which today has all the characteristics of a 'novel ecosystem' (sensu Hobbs et al., 2013a). Such ecosystems emerge through direct or indirect human agency that has resulted in an irreversible departure from historical baseline conditions (e.g., due to alterations related to climate warming and biological invasions); in addition, they are characterized by novel combinations of species and/or functional properties, with ramifications for important ecosystem processes (Hobbs et al., 2013b).
Our results therefore suggest that non-indigenous assemblage components cannot maintain or restore the functions lost due to the massive regional collapse of native populations (e.g., Albano et al., 2021;Rilov, 2016), despite their often considerable functional richness, which was on par with that of their native counterparts in all hard substrate assemblages, and even exceeded it at 50% of the soft substrate stations. This interpretation was further supported by the analysis of native versus non-indigenous functional beta diversity: at all stations, beta diversity was high (mean value: .75) and mainly related to functional turnover, which -on average -contributed 90% to the Jaccard-like dissimilarity index. In addition to these general patterns, native versus non-indigenous trait space occupation also showed habitat-specific differences: on hard substrates, nonindigenous molluscs always had a greater functional divergence (logindex ratios < 0) and lower functional evenness (log-ratios > 0). This may indicate a greater degree of niche differentiation and a patchier distribution of species and abundances in trait space, respectively (Mouchet et al., 2010;Mouillot et al., 2013;Villéger et al., 2008).
On the other hand, functional turnover, and hence overall beta di-  (Nawrot et al., 2017); this likely contributes to the observed greater share also of large-sized species among nonindigenous components.
Lessepsian molluscs also directly affect trophic networks in the Eastern Mediterranean by modifying feeding guild composition.
Non-indigenous assemblage components across subtidal habitats were often characterized by a greater share of surface deposit feeders, whereas suspension feeders were better represented among native components from hard substrates. The relative abundance of these functional groups can have important implications for pelagic-benthic coupling, sediment stability and biogeochemical fluxes at the sediment-water interface (Griffiths et al., 2017;Rhoads & Young, 1970;Snelgrove, 1997). Among higher-order consumers, a greater proportion of ectoparasites/browsing carnivores (and consequently host-associated individuals) was found among nonindigenous molluscan components from hard substrates. A surprisingly large number of parasitic gastropods, particularly of the family Pyramidellidae, has been introduced to the Mediterranean Sea (Oliverio & Taviani, 2003), but hardly anything is known about their hosts. While parasites may be co-introduced with their original hosts (Boussellaa et al., 2018;Galil, 2007), it has been observed that Indo-Pacific pyramidellids also feed on native Mediterranean molluscs (so-called 'parasite spillover'; Oliverio, 1994), with unknown implications for the population dynamics of these new hosts. Last but not least, we found a greater importance of epifaunal lifestyle for non-indigenous assemblage components in most soft-bottom habitats, suggesting potential implications for the vulnerability to predation, sediment oxygenation and biogeochemical processes, for example via changing rates of bioturbation (de Moura Queirós et al., 2011;Kristensen et al., 2012;Vermeij, 1987).
With ongoing warming, thermophilic NIS will clearly be the winners in the future Eastern Mediterranean , likely further increasing the functional deviation of shallow water assemblages from their historical baselines. The interplay of massively declining native biota, booming non-indigenous populations and constant novel introductions through the Suez Canal suggests that the composition of local benthic assemblages will likely continue to evolve rapidly, with potentially profound but little understood consequences for ecosystem services. and P.G.A. J.S. and P.G.A. wrote the initial draft. All authors contributed to the interpretation of the results and manuscript writing.

DATA AVA I L A B I L I T Y S TAT E M E N T
Raw data and R scripts are available from the figshare repository: