Structural and functional effects of global invasion pressure on benthic marine communities—patterns, challenges and priorities

Retrospective (pre‐ vs. post‐invasion) and cross‐sectional comparisons of ecosystems exposed to high and low bioinvasion pressure, provide an alternative approach to evaluate shifts in biological communities associated with non‐indigenous species (NIS) introductions. In this study, we aimed to examine general patterns of change in community composition, structure and function in six well‐studied and globally distributed marine ecosystems that had documented histories of biological invasions.


| INTRODUC TI ON
Biological invasions are a defining characteristic of global change and the Anthropocene epoch in general (Leroy et al., 2023;Lewis & Maslin, 2015;Ricciardi, 2007).Since the middle of the twentieth century, the number of species that have been detected outside their native historical range has increased dramatically across most of the world's terrestrial and aquatic ecosystems, in association with increased global trade and travel (Bailey et al., 2020;Early et al., 2016;Seebens et al., 2018).Over that period, research on the ecological and economic consequences of non-indigenous species (NIS) has increased correspondingly (IPBES, 2019;Pejchar & Mooney, 2009;Simberloff et al., 2013).Despite decades of research on bioinvasions, documenting new incursions and the consequences of many individual NIS (Carlton, 1999;Galil, 2018;Guy-Haim et al., 2018;Katsanevakis et al., 2014;Ruiz et al., 1997;Seebens et al., 2017), we still lack understanding of the chronic, cumulative impacts of multiple NIS on communities and ecosystems (Ojaveer & Kotta, 2015;Ruiz et al., 1999).Most of our understanding of invasion biology is at the population level, focused on a single NIS, a restricted geographical area and/or timeframe (Strayer et al., 2006;Watkins et al., 2021).
There is a lack of standardised longer-term community data from the invaded environments, that can help us interpret how community assembly and ecosystem functions change over time in response to multiple invasions (Carlton, 2009).Impacts of marine NIS are known to vary from species-to-species and place-to-place, but most marine invasions have likely not reached equilibrium and the size of the invaded range for NIS is strongly associated with 'time since arrival' (Byers et al., 2015;Galil, 2021).Measuring and evaluating this dynamism at the community or ecosystem level over reasonable timeframes would help reveal the extent to which invasions contribute to marine community change and to variation in ambient community effects, including ecosystem goods and services, of this global scale phenomenon.
A great deal of our current understanding of the impacts of individual marine NIS has been derived from: (i) autecological studies, producing one-off observational datasets of the ecological profile of a NIS (e.g.Firth et al., 2021;Lutz-Collins et al., 2009;Olenina et al., 2010;Orlova et al., 2004); (ii) autecological/synecological experimental studies, where interactions with native species and assemblages are deduced by manipulative experiments by adding, removing, altering abundance or using physical mimic models of natives or NIS (e.g.Atalah et al., 2019;Floerl et al., 2004;Giddens et al., 2014;Hollebone & Hay, 2008); and (iii) presence-absence synecological studies, usually observational and/or experimental comparisons within a site, when diversity, structure and function of the communities inside and outside the NIS-affected habitat are examined (e.g.Guilhem et al., 2020;Ross et al., 2006;Zaiko et al., 2009).
However, to evaluate the chronic and cumulative effects of multiple invasions in the absence of genuine and consistent long-term datasets, retrospective studies, contrasting baseline (pre-invasion) to subsequent (post-invasion) data (e.g.Floerl et al., 2009;Forrest & Taylor, 2002;Steger et al., 2022); and cross-sectional studies, contrasting locations historically exposed to high propagule pressure and invader colonisation (e.g.due to proximity to a propagule source, Johnston et al., 2009) may offer greater insights (Strayer et al., 2006).
Trait-based analyses linking species abundances to a mixture of life history, morphological and behavioural characteristics of species confounding effects of study designs and other drivers of change.The most prominent shifts in community composition were observed in the retrospective studies, characterised by the greatest relative contribution of NIS.No uniform pattern of change in functional metrics was observed across study regions.However, functional evenness and dispersion showed a tendency to increase in systems under higher invasion pressure, refuting the hypothesis of selective accumulation of specific traits and functional homogenisation within ecosystems exposed to high invasion pressure.
Main Conclusions: Accumulation of NIS within broader communities can be a subtle process, with inherent spatial and temporal variability.Nonetheless, not only do species' proportional contributions to communities change over time in areas subjected to high bioinvasion pressure, but trait profiles can incrementally shift, which alters the original ecology of an area.Planned, long-term studies that incorporate a range of measures of environmental drivers and ecosystem response are crucial for better understanding of cumulative, community-level and ecosystem-scale change associated with biological invasions.

K E Y W O R D S
benthic communities, Biological invasions, cumulative effects, historical datasets, nonindigenous species, trait analysis in a community (Bremner, 2008) may improve our understanding of community functioning among greatly varying locations at broad spatial scales (Statzner et al., 2001).Much of ecological theory predicts that species that share the same environment differ in their trait characteristics to reduce competition pressure and facilitate long-term coexistence (Pianka, 1978;Schoener, 1974).Recent mass invasions (Sax et al., 2007), however, have proved the apparent 'unsaturation' of communities indicating our generic lack of understanding on basic assembly rules of ecological systems.Successful settlement and further dispersal of non-indigenous species are determined by their ability to colonise and retain a niche in the environment (MacDougall et al., 2009).There exist many alternative views on how species traits and species invasiveness are related (Duncan & Williams, 2002;Marvier et al., 2004;Pires-Teixeira et al., 2021;Steger et al., 2022;Strauss et al., 2006).However, the combination and diversity of traits that determine successful establishment will only become apparent after the establishment of a permanent NIS population in a new region (Jiménez-Valverde et al., 2011).
Thus, the analyses of impacts of invasive species should not only focus on community diversity, but the functional trait diversity of the community and its relationship to ecosystem processes (Diaz & Cabido, 2001).Expansion or contraction of the functional space as a result of bioinvasion may signal substantial change in ecosystem functioning, susceptibility to further invasions, sustained biodiversity and provisioning of ecosystem services (Funk et al., 2008;Milanović et al., 2020;Wen et al., 2019).
In this study, we examine patterns of structural and functional community-level change in a range of well-studied marine ecosystems with documented histories of bioinvasion.We hypothesised that accumulation of NIS in ecosystems exposed to high bioinvasion pressure triggers structural divergence in the affected communities which in turn leads to functional changes.Specifically, we expect that in assemblages with higher exposure to NIS colonisation pressure: (i) substantial structural changes with higher relative NIS contributions will be observed (Bradley Bethany et al., 2019;Olenin et al., 2007); (ii) novel functions are likely to be introduced, thus more niche space will be occupied (Parker et al., 1999;Thomsen et al., 2011); (iii) the average functional dissimilarity among sites will decrease due to selective accumulation of certain species and traits across biological communities (McDowell & Byers, 2019;Smart et al., 2006).We evaluate trends within and among ecosystems by comparing paired measures of community and functional structure in either space or time.We acknowledge that these are not fully controlled measures since they can respond to both accumulation of NIS and other drivers of change, such as climate change, habitat loss, pollution, sedimentation or overfishing.
By considering a range of regional datasets and different sampling approaches, we therefore test for concordance and robustness of any detected patterns to help build our understanding of ambient community change associated with globally observed levels of invasion pressure.
This improved understanding can assist resource managers with justification and implementation of conservation strategies.Perhaps more importantly, it can also provide critical guidance for future global environmental reporting such as that undertaken by many nations to fulfil their obligations under the Convention of Biological Diversity or other environmental commitments (Department of Conservation, 2019;Dobson, 2005;Lehtiniemi et al., 2015).

| Compilation of species datasets
Two workshops attended by 20 international marine invasion experts were held in 2016 (Australia) and 2018 (Argentina) to investigate the availability of temporal (long-term) and cross-sectional datasets for global coastal (marine and estuarine) locations associated with biological invasions.We identified six regions with extensive bay-scale datasets on native and non-indigenous benthic species assemblages, allowing paired comparisons in time (years to decades; retrospective datasets) or space (high vs. low proximity to hotspots of NIS introductions within a region; cross-sectional datasets) representing differences in bioinvasion pressure (Table 1).The magnitude of bioinvasion pressure was approximated by the extent of relevant drivers and mechanisms of environmental change (Oesterwind et al., 2016) that can lead to the accumulation of NIS.Assignment of the locations to a low-or high-pressure category was a case-by-case expert decision based on a qualitative assessment of indirect drivers likely contributing to the introduction and establishment of non-native species in an ecosystem (e.g.shipping intensity, mariculture activities, tourism and recreation).Where available, quantitative data sources were utilised to underpin those decisions.The selected ecosystems were as follows: (1) coastal waters of British Columbia, Canada (BC); (2) San Francisco Bay, USA (SF); (3) Ilha Grande Bay, Brazil (BR); (4) North-Eastern Baltic Sea, Estonia (BS); (5) estuaries of New South Wales, Australia (AU) and (6) Waitematā Harbour, New Zealand (NZ).These six regions are generally at mid-to higher-latitudes, except for one low-latitude BR.The ecosystems encompassed surveys of benthic communities: fouling assemblages (BC and AU), subtidal reefs (BR and BS) and soft-sediment benthos (SF, BS and NZ).
Two cross-sectional datasets (BC and AU; Table 1) used single-time designs (sensu Wiens & Parker 1995) in which multiple locations that have historically been subjected to high NIS colonisation pressure from shipping or other vectors (high-pressure areas) were contrasted with nearby reference locations in settings not adjacent to significant shipping ports and thus subjected to lower bioinvasion pressure (low-pressure areas) (Johnston et al., 2009;Ruiz et al., 2000).The highpressure areas typically had major ports of entry for international vessels, whereas low-pressure ones were characterised by limited vessel movement (predominantly domestic).The locations were selected so that both high-and low-pressure areas featured similar environmental conditions (i.e.salinity, temperature regimes, habitat types).The other four regions (BR, SF, BS and NZ) provided data for paired contrasts of temporal datasets along the NIS accumulation curve, that is, a historic survey and similar or identical contemporary survey, with variation in the period of exposure to NIS (Table 1).For consistency with the crosssectional studies, hereafter, these are referred to as high-pressure and TA B L E 1 Overview of the datasets considered in this study.low-pressure datasets respectively.We acknowledge that high-pressure and low-pressure is not a strict unambiguous dichotomy, as it is wellknown that small estuaries and harbours lacking international shipping may be highly invaded, due to significant contribution of both intraregional transport and the non-shipping vectors (Wasson et al., 2001;Zabin et al., 2014).

| Traits database
A set of species traits relevant to important ecosystem functions was identified and modalities for each trait category (representative of different types of benthic invertebrates) were defined by experts attending the two workshops and were guided by the approach  m AquaNIS Editorial Board (2015).
o Hayward et al. (1997). p The Ministry for the Environment and Statistics New Zealand (2014).
TA B L E 1 (Continued) adopted in Marchini et al. (2008) (Table 2).Binary (0/1) trait information was compiled for all species in the datasets (only adult or benthic stages were considered for species with a complex life cycle).Where available, peer-reviewed scientific publications were used for the collation of trait data for each species.In all other cases, available information was retrieved from relevant on-line databases (e.g.www.corpi.ku.lt/ datab ases/ aquan is/ , http:// www.mollu sca.co.
eu/ ), grey literature and expert judgement (where no published information was available).

TA B L E 2
Categories of species traits and their modalities considered in this study.

Life form
Zoobenthos-animals living on or in the seabed

Phytobenthos-algae and higher plants living on or in the seabed
Demersal-animals living on or near seafloor, able to move about in water Parasite-an organism intimately associated/dependent on another living organism Symbiont-an organism living mutually with another species without harming it

Trophic position
Autotroph-an organism obtains metabolic energy from light by a photochemical process such as photosynthesis

Mixotroph-an organism both autotrophic and heterotrophic
Suspension Feeder-an organism feeds on particulate organic matter from the water column Deposit Feeder-an organism feeds on fragmented particulate organic matter from the substratum Omnivore-an organism feeds on a mixed diet including plant and animal material Herbivore specialist-a herbivore that feeds on specific type of plant material Herbivore generalist-a herbivore that feeds on a variety of different types of plant material Predator specialist-a predator that feeds on a specific type of animal prey Predator generalist-a predator that feeds on a specific type of animal prey (includes scavengers)

Mobility
Sessile encrusting-attached to substrate, cover with a crust or thin coating

Sessile turfing-low growing erect or filiform organisms
Sessile tubiculous-forms a structure/tube in which it lives Sessile reef-builder-forms consolidated biogenic habitat on the seabed or shore

Sessile erect-upright in position or posture
Swimmer-an organism capable of moving through the water by means of fins, limbs or appendages Crawler-an organism that moves along the substrate Burrower-an organism capable of digging in sediment/soft substrate Borer-an organism capable of penetrating a solid substrate by mechanical scraping or chemical dissolution

Body surface
Robust-heavily calcified or leathery, unlikely to be damaged by physical impacts Fragile-lightly calcified, easily damaged as a result of physical impact or pressure

Rigid-chitinous endo-or exo-skeleton
Soft-yields to the touch or pressure b Species with an assigned size modality <1 mm were excluded for SF dataset, as different sieve sizes (0.5 and 1 mm) were used in the low-pressure and high-pressure datasets respectively. c The habitat modification category was not considered for BC and AU datasets as it was considered irrelevant for settlement plate communities.
The trait categories included in the analyses were: body size, life form, trophic position, mobility, body surface, habitat modifi- excluded from the regional datasets and downstream analyses, as no reliable generalisations on biological traits could be made at that level.
In total, nine taxa were removed across all datasets.

| Quality assurance and pre-processing of the datasets
The six datasets were groomed by removing inconsistencies in nomenclature (standardised against the WORMS http:// www.marin espec ies.org and ITIS https:// www.itis.gov databases).For analyses where quantitative data were considered, community datasets were transformed into relative abundances (% of total community) and rank abundances.In the SF dataset, different sieve sizes (0.5 and 1 mm) were used for sample processing in 1987 and 2012 respectively.To mitigate the potential bias introduced by this methodological inconsistency, the small-bodied meiofaunal species (e.g.cumaceans or ostracods) were removed from the datasets and not considered in the analyses.
After compilation of the traits database, the data were crosschecked for inconsistencies by comparing the trait profiles of closely related species (family level) and addressing any mismatches through additional literature searches and reviews.For quantitative traits analyses, a weighted traits matrix was created for each dataset and computed as a cross-product of binary traits data and transformed species data matrices.Non-quantitative comparisons between high-and low-pressure datasets within each region were made on occurrence of individual traits and trait profiles, which are the combination of traits expressed by individual species.

| Exploratory analyses of biological assemblages and traits space
Since the datasets did not meet the normality assumptions, nonparametric Kruskal-Wallis tests were used to identify significant differences in the relative abundance of NIS in low-and high-pressure datasets for each study region (see Table 1).For each global study region, structural (species data) and functional (traits data) shifts in species assemblages across the pressure (low vs. high) and spatial (study area) factors (see Table 1) were investigated using permutational multivariate analysis of variance (PERMANOVA) based on Bray-Curtis similarity matrices (species or weighted traits).
PERMANOVA was performed using the adonis2 function of the vegan package (Oksanen et al., 2019).A crossed-factor design (except for AU and NZ datasets, where 'area' was nested within the 'pressure' factor) with 999 permutations was applied.Hereafter, we refer to the 'area' factor as a geographical domain within the study region (i.e. an estuary, or a distinct part of a larger basin).
To visualise the multivariate structure of species and weighted traits data, principal component analysis (PCA) biplots were produced for each global study region using the fviz_pca_biplot function within the factoextra package (Kassambara & Mundt, 2017).The function multipatt of the package indicspecies (De Caceres & Legendre, 2009) was used to determine subsets of traits that were indicative of either lowpressure or high-pressure datasets.This approach allows determining indicator species (traits in our case) using an analysis of the relationship between the occurrence or abundance values from a set of sampled sites and the classification of the same sites into site groups, which may represent habitat types, community types, disturbance states, etc.The Indicator Value index measuring the association between a trait and a pressure-relevant group was calculated on the weighted traits matrix for each region.The statistical significance of this relationship was then tested using a permutation test, based on 999 permutations.
We used multidimensional indices of functional composition and diversity (FRic, FEve and FDis, Table 3) to compare changes in trait space among the temporal and spatial samples and explored the contributions of NIS to the changes.All functional diversity metrics were computed in the FD package (Laliberté & Legendre, 2010).
Non-parametric Kruskal-Wallis tests were used to identify statistically significant differences in functional diversity metrics between high-and low-pressure datasets.
In addition to exploring the general functional diversity and describing the functional space using functional diversity metrics, we also considered 'trait profiles'-the full multidimensional combination of all trait modalities exhibited by a species or taxon (i.e.functional species).We assumed that shifts in the trait profiles of taxa represented in a community might have functional implications at ecosystem level.Therefore, to better understand the different levels of potential functional shift, we assessed the magnitude of overlap in both individual traits space and trait profiles for each ecosystem's low-and high-pressure datasets.
All analyses, calculations and visualisations were performed in R v.3.5 (R Core Team, 2014).

| Overview of patterns associated with biodiversity and traits distribution
The number of reported NIS varied across locations from 0 in the lowpressure NZ dataset (that is, no reported invasions in the 1920s-1930s era) to 26 in the high-pressure SF dataset.Significant differences in the average relative abundance of NIS between low-and high-pressure datasets were only detected in the SF and NZ studies, with both regions experiencing higher NIS contributions to abundance over time (Figure 1).
In most datasets, the structure of species assemblages was significantly affected by both the 'pressure' and 'area' (spatial) factors or their interaction term (Table 4).The exception was the NZ dataset, where there were only pressure-related differences.The functional structure (weighted traits) of communities was also significantly affected by both main effects of the 'pressure' and 'area' factors or their interaction (except NZ where only 'pressure' had a significant effect, Table 5).
However, the amount of unexplained variation (residual R 2 ) was high (62%-91%) for all regions, habitats, species and traits datasets.
Overall, differences in species composition between groups were likely driven by dispersion of the data, rather than a shift in their centroids (Figure 2 and Figures S1-S6).For example, in the AU, BR and NZ datasets, high-pressure communities were structurally more variable compared to the low-pressure ones, reflected in plots that had much larger 'multivariate species space' (represented by dispersion of samples) for recent surveys compared to historic ones (Figure 2).The opposite pattern was detected for SF however, where compositional dispersion of the benthic community substantially contracted in the high-pressure dataset relative to the low-pressure one (Figure 2 and Measures the amount of functional space occupied by a species assemblage and is naturally positively correlated with the number of species present (the more species there are, the larger the functional space occupied when species traits are somewhat randomly distributed).However, two communities with the same number of species may have different FRic when functional traits of species are more closely clustered in one community than in the other a .
Calculates the volume of trait space with the convex hull volume, which represents the smallest convex hull that encloses all species.With a complex algorithm, the most extreme points (vertices) can be determined and the volume encompassed by these vertices is calculated.Not weighted by species abundance d .
With accumulation of NIS, novel functions are likely to be introduced, thus the FRic will increase (more niche space is occupied).

Figure S4
). Across all regions, NIS were among the 10 species that contributed most to differentiating the principal components (highly correlating with one of two PCA main axes, Figures S1-S6) and generally were associated with high-pressure samples.Only the SF dataset had major contributions of three NIS that correlated with low-pressure samples from the southern area (Figure S4).Pressure-related changes in traits, where detected, were associated with shifts in centroids (multivariate mean) rather than dispersion, and this was most pronounced in the SF and NZ datasets (Figure 3 and Figures S4 and S6).There was substantial overlap of traits space in the four other datasets (Figure 3, Figure S1-S6).In fact, when looking into individual trait occurrence within each region, most (84%-100%) were shared between low-pressure and high-pressure datasets (Figure 4).Trait modalities unique to the low-pressure assemblages within a particular region were as follows: fragile 'body surface' (BR), parasite 'life form' (SF), borer 'mobility' (SF) and 'trophic positions' of predator specialist (BS) and herbivore specialist (NZ).
When considering the full combination of all trait modalities exhibited by a species or taxon ('trait profiles'), 10%-46% were unique for low-pressure samples and 15%-32%-for high-pressure samples (Figure 4).This suggests that while regional trait pools remained largely unchanged in the studied ecosystems, the representation of trait modalities in individual taxa shifted markedly between low-and high-pressure datasets.For example, in BC, 100% of individual traits were shared between low-and high-pressure datasets, but 76 (46%) trait profiles (that occurred in 91 taxa) were reported exclusively in the low-pressure dataset (Figure 4).
Analysis of indicator traits revealed a higher number of trait modalities significantly associated with the high-pressure (35) than the low-pressure ( 26) datasets (Table 6).No general pattern emerged for trait occurrence or change in traits across regions and suites of indicator traits were region-specific.No significant indicator traits were detected for the BC and BR datasets, which is consistent with PCA results (see Figure S1, Figure S3).There were more traits significantly associated with low-pressure datasets in SF and BS, whereas high-pressure trait associations were more common in AU and NZ.
NIS prevalence (the binary biogeographic trait) was only a statistically significant indicator in two high-pressure datasets (SF and NZ), where a significant increase in the relative abundance of NIS was also reported over the covered timeframe (Figure 1).

| Functional diversity
Functional diversity analysis returned differential responses in functional diversity metrics across regions (Figure 5).A significant increase in functional richness was only evident in AU and NZ datasets, while in SF it had substantially decreased.These results align with indicator TA B L E 5 Results of PERMANOVA (computed by R function adonis2) for trait data from six considered datasets: British Columbia (BC), New South Wales Australia (AU), Ilha Grande Bay (Brazil), San Francisco Bay USA (SF), North-Eastern Baltic Sea (BS), Waitematā Harbour New Zealand (NZ).trait associations in Table 4. Functional evenness and dispersion tended to increase across datasets (except BC), with significant change in SF (both evenness and dispersion), AU and NZ (dispersion only).Overall, the strongest shifts between low-and high-pressure datasets were detected in AU, SF and NZ, with no consistent effect for all three of these regions, and no significant changes apparent for the other three regions.SF and NZ had significant changes in NIS contributions to community abundance (Figure 2) but AU did not.The dispersion of species composition (in multivariate space) was higher in high-pressure datasets for all three regions (Figure 2) but shifts in centroids for species traits only occurred for SF and NZ (Figure 3).Indicator traits were more heavily associated with high-pressure datasets for AU and NZ, but were associated with the low-pressure dataset for SF (Table 4).SF had a loss of functional richness from low-to high-pressure datasets, while AU and NZ gained functional richness (Figure 5).Overall, out of the seven significant differences observed in functional diversity metrics, six were greater in the high-pressure datasets.

| DISCUSS ION
Our study shows that the accumulation of NIS within broader communities can be a subtle process, with inherent spatial and temporal variability that, nonetheless, can drive underlying shifts in community and ecosystem characteristics.Not only do species' proportional contributions to communities change over time in areas subjected to high bioinvasion pressure, but trait profiles can incrementally shift, which alters the original ecology of an area.A notable and somewhat unexpected finding was that the assumed bioinvasion pressure is not always unambiguously manifested in the available NIS accumulation data.This emphasises that NIS-focused surveys might not, by themselves, represent the community-level changes adequately and more holistic approaches are imperative to measure the long-term consequences of bioinvasions and other stressors for marine ecosystems and assess the bioinvasion-related status.

| Differential patterns of invasion-related change in community structure
Our analyses revealed that invasion-related change in communities was not consistent across the range of systems examined in this study.There are good reasons to expect community-level effects of NIS to be spatially and temporally variable (Bracewell et al., 2021;Clark & Johnston, 2011).Since human activity became a dominant force in global biotic exchange, the structural F I G U R E 2 Two-dimensional PCA visualisations of species composition from low-pressure and high-pressure datasets across all considered regions: BC-British Columbia coastal waters, Canada; AU-New South Wales estuaries, Australia; BR-Ilha Grande Bay, Brazil; SF-San Francisco Bay, USA; BS-North-Eastern Baltic Sea, Estonia; NZ-Waitematā Harbour, New Zealand (see Table 1 for details).The concentration ellipses cover 95% confidence interval for each group of samples.
composition of assemblages changed through the successive introduction and establishment of NIS within local species pools as well as in response to other changes and stressors in receiving environments, including coastal hardening, maritime sprawl, pollution, commercial fishing and warming temperatures (Floerl et al., 2021;Hopkins et al., 2021;Occhipinti-Ambrogi, 2007).Significant timelags in population development following introduction (Crooks & Soulé, 1996;Guastella et al., 2021), 'boom and bust' dynamics (Simberloff & Gibbons, 2004), seasonal (Schiel & Thompson, 2012) and interannual variation in recruitment (Crooks, 1996), perenniality (Thibaut et al., 2004) and facilitative effects of prior invasions (Grosholz, 2005;Zaiko et al., 2007) can all affect the outcome of NIS incursion and magnitude of manifested impact.Following ecosystems over long time periods (decades or longer, ideally) may enable detection of shifts in communities that proceed at a slower rate or involve time lags (Ojaveer et al., 2021).Therefore, we expected the retrospective studies to reflect changes in the regional species pool more reliably than the cross-sectional studies.Indeed, the most prominent shifts in community composition were observed in temporal SF and NZ datasets, with some contrasting patterns that were nonetheless characterised by the greatest relative contribution of NIS in the high-pressure dataset.
The pronounced changes in SF community structure were largely determined by high dimensional dispersion in low-pressure samples (Figure 3).This was driven mostly by differences in the relative contributions of three NIS (the polychaetes Heteromastus filiformis and Streblospio benedicti, and the bivalve Mya arenaria) in the southern basin (Nichols & Thompson, 1985;Robert, 1881).All three represent invasions that were established in the ecosystem well before the (historic) low-pressure dataset was acquired.It is likely that their long-term impact diminished over decades with ubiquitous spread of new NIS and accumulation of other anthropogenic pressures and environmental changes in the region, resulting in contracting variability between the north and south basins of the bay (Ely & Owens Viani, 2010).In contrast to SF, substantially higher dispersion was observed in the NZ high-pressure dataset compared to the low-pressure one.A significant shift in community composition was driven by two non-indigenous molluscs (Theora lubrica and Limaria orientalis), introduced into Waitematā Harbour in the 1970s (Hayward et al., 1997) with strong tolerance of sedimentation.It is likely that anthropogenic perturbations in the ecosystem over the last few decades, particularly sediment accumulation, have altered densities of large native bivalves and other functionally important benthic species, thus facilitating a shift in the prevalence of these NIS in muddy subtidal habitats F I G U R E 3 Two-dimensional PCA visualisations of species traits from low-pressure and high-pressure datasets across all considered regions: BC-British Columbia coastal waters, Canada; AU-New South Wales estuaries, Australia; BR-Ilha Grande Bay, Brazil; SF-San Francisco Bay, USA; BS-North-Eastern Baltic Sea, Estonia; NZ-Waitematā Harbour, New Zealand (see Table 1 for details).The concentration ellipses cover 95% confidence interval for each group of samples.(Lohrer et al., 2008).The relatively long timelines between low-and high-pressure data collection, supported by published research on other environmental dynamics of these regions, provide context and an understanding of some underlying mechanisms for how invasions have progressed with other factors to shift community composition of these two systems (Atalah et al., 2019;Hayden et al., 2009;Jimenez et al., 2018;Kerr et al., 2016).In both cases, benthic community baselines appear to have shifted dramatically more than the other regions based on the current data.
NIS responding to environmental change as well as driving community shifts is an inevitable characteristic of chronic invasion effects (Bauer, 2012;Didham et al., 2005;MacDougall & Turkington, 2005;Vitousek et al., 1997).NIS can take advantage of human-induced impacts on ecosystems and homogenisation of habitats, potentially contributing further to the decline of native species and losses of ecosystem functions (Byers, 2002;Piola & Johnston, 2008).For instance, the Baltic Sea has faced a major increase in anthropogenic pressures and partly human-induced regime shifts that started before the historical (low-pressure) dataset was collected (Österblom et al., 2007).Therefore, the communities examined, already subjected to substantial environmental changes in the ecosystem, exhibited a rather smoothed response to the chronic effects of accumulation of NIS over the period spanned by the datasets (Hewitt et al., 2016).This highlights the difficulty or impossibility of separating confounding factors from retrospective community analyses of invasions.
The BR community comparison had apparent shifts in dispersion of species and traits (Figures 3 and 4), as occurred in NZ, but no significant differences in relative NIS abundance (Figure 2) or functional diversity metrics (Figure 5), as were observed for the Baltic Sea.BR was considered a retrospective dataset but represents a comparatively short timeframe and primarily captures spread, rather than accumulation, of NIS in the region.The sampling in this system focused on range expansion of two sun coral species (Creed et al., 2017;Silva et al., 2014) and their possible interactions with-and effects onnative species.Although there did not appear to be major shifts or invasion-driven changes in traits, an overall expansion of community dispersion between low-and high-pressure sampling periods was noteworthy.This may be linked to the effect of Tubastraea coccinea on community structure in the high-pressure samples.It was noted previously by Guilhem et al. (2020), that an increase in sun coral cover in the invaded areas was associated with intensified turnover (i.e. higher spatial heterogeneity) of native species.This likely explains the increased dispersion in the high-pressure species data.
The cross-sectional studies considered here (AU, BC) are typical of single-time, post-impact studies that compare ecological assemblages at sites subject to a perturbation with sites that have not been exposed to the stressor ('Impact-Reference' designs) or which compare sites that have experienced different levels of exposure ('Gradient' designs; Eberhardt and Thomas, 1991;Wiens and Parker, 1995).In both cases, the presence of long-established NIS in the reference (low pressure) locations likely confounded comparison with high-pressure sites so that any invasion-related impacts are indistinguishable from other influences and local processes affecting community assembly (Davis et al., 2005).A related issue is that contemporary proxies of colonisation pressure, such as the numbers of vessel arrivals or recent ballast discharge volumes at a location, may fail to characterise historical patterns of invasion, since high-use environments may acquire NIS from other, uncharacterised pathways (Bailey et al., 2020;Ojaveer et al., 2018).The confounding effect of the pre-established NIS may partly explain comparable and even somewhat lower relative abundances of NIS in the high-pressure cross-sectional datasets (BC and AU).On the other hand, although caution was taken to include data for paired assessments only from the ecosystems with comparable environmental conditions (e.g.temperature and salinity ranges), we were unable to control for all possible biogeographical differences.For instance, the Pacific coast of Canada (represented by the Strait of Georgia dataset here) is by far more complex than the North Coast.
This is most likely due to the environmental influence of the Fraser River but also likely reflects highly variable marine use (shipping, recreational boating, aquaculture, etc.) which might have an effect on the patterns observed in NIS and wider benthic communities.
Still, both in BC and AU datasets, the pressure factor had a significant effect on variances in community composition, with major discrepancies between low-and high-pressure sites driven by NIS (Figure 2, Figures S1 and S2).

F I G U R E 4
Partitioning of the individual trait modalities (upper graph) and trait profiles (species-specific combination of all trait modalities, bottom graph) between low-pressure and high-pressure datasets from six regions: British Columbia (BC), New South Wales Australia (AU), Ilha Grande Bay Brazil (BR), San Francisco Bay USA (SF), North-Eastern Baltic Sea (BS), Waitematā Harbour New Zealand (NZ).

| Shifts in trait and functional make-up of benthic communities under invasion pressure
Our analyses revealed a significant effect of the invasion pressure factor on trait distribution within most of the regions examined (Table 5).NIS may bring new traits and novel functions to an ecosystem (Parker et al., 1999;Ruesink et al., 1995;Thomsen et al., 2011), increasing functional richness (as observed in e.g.AU and NZ datasets).Alternatively, regional ecosystems will likely favour invaders with similar characteristics to those of the recipient community, including NIS already present (Duncan & Williams, 2002).Thus, as colonisation pressure (sensu Lockwood et al., 2009) increases, the average functional dissimilarity among sites might decrease due to selective accumulation of successful invaders (biotic homogenisation), replacement of native species, or prevalence of particular biological traits (Pires-Teixeira et al., 2021;Smart et al., 2006), without any apparent change in functional richness.Decreases in functional distinctiveness, in turn, can increase vulnerability to broadscale perturbations by synchronising local biological responses (Olden et al., 2004) and removing functional redundancy or community resilience, ultimately risking prolonged loss of diversity by restricting recolonisation capacity (Clavero & Garcia-Berthou, 2005).The SF dataset was the only one with a significant reduction in functional richness, which was combined with a significant increase in functional evenness (also unique to SF; Figure 5), suggesting a simplification of benthic communities with NIS playing a major role.Contrary to our hypothesised increase in functional richness associated with elevated bioinvasion pressure, there was no uniform pattern across study regions.However, all three measures increased in high-pressure datasets in six of the seven detected significant differences and there was also a reasonable concordance in functional evenness and dispersion response to the pressure factor.Both measures showed a tendency to increase in the high-pressure datasets (except for BC), thus not supporting our hypothesised selective accumulation of certain traits and functional homogenisation in the invaded ecosystems (Table 3).
Overall, no single observed trait shift was consistent among all regions-any changes were context-specific.However, across all datasets, a substantial change in functional profiles (combinations of traits) within the largely overlapping functional space of individual traits was detected.It has been shown previously that the same environmental drivers define the ecological constitution of both native and NIS communities (Floerl et al., 2009;Pysek et al., 2020).However, certain combinations of traits might convey adaptive advantages for native species under increasing pressures and invasive success of NIS (Boltovskoy et al., 2021;Novoa et al., 2020;Nunez-Mir et al., 2019;Reichard & Hamilton, 1997).
This could explain the phenomenon observed in our study and most profound in BC (Figure 5), where no individual trait was characteristic or distinctive to low-or high-pressure communities, F I G U R E 5 Overview of three functional diversity metrics examined (see but trait profiles (the combination of traits within a species) were highly related to one or the other pressure state.The cumulative portrait for a unique species of BC changed from long-lived to short-lived species that were smaller in size with a lower signal for deposit-feeding (illustrated in Figure S2 for both low-and high-pressure locations).Such change, although quite subtle within a whole-community context, can lead to substantial shifts in ecosystem's functional characteristics and service provision.For example, cumulative effect of a small filtrator and a large deposit-feeder will not be equal to that of a small depositfeeder and a large filtrator (although the traits remain the same for both cases).This means that species identity is important and functional changes should be considered and interpreted in conjunction with community structure assessments.
In contrast to other studied regions, there was a very large overlap in multivariate trait space between low-and high-pressure BS datasets.The Baltic Sea ecosystem is known for numerous strong spatial environmental gradients as well as multiple shifts in environmental conditions over recent and geological times (Österblom et al., 2007;Zettler et al., 2014).As a consequence, its communities primarily consist of very tolerant and opportunistic species, that is, traits that are often found among invasive species (Byers, 2002;Piola & Johnston, 2008).Although the Baltic Sea hosts a great number of non-indigenous species relative to its total species richness, invader traits were not unique to the 'trait space' of resident species, and thereby long-term changes in native communities were not as severe as in many other studied ecosystems.

| Further considerations for disentangling bioinvasion effects in the context of globally changing marine ecosystems
Comprehensive functional-taxonomic community characterisations are challenging and, in the context of our study, hindered by a lack of comparable, long term, whole-community datasets that capture invasion dynamics within a range of other drivers of change.Our results show that despite the exponential growth of bioinvasion studies over the last several decades (Ojaveer et al., 2021;Ruiz et al., 2000), we are still facing major challenges in quantifying and communicating the chronic, community-and ecosystem-scale effects of introduction for one of the best studied communities-macrozoobenthos-in our global marine environment.
The challenge ahead is to source historical datasets that are in a form that would lend themselves to contemporary comparison of functional change.This is likely to be difficult because many historical datasets were natural history-type inventories (focused on richness) that lacked measures of abundance, which is important for comparisons of structural composition.In parallel, there are emerging longer-term repeated-measures datasets of marine community dynamics that highlight the interacting roles of invasion and environmental perturbations in benthic systems (Chang et al., 2018;Nygård et al., 2020;Philippe et al., 2017).In NIS-focused studies, native context is often lacking, resulting in disconnected records of native and non-indigenous biodiversity for the same time and place.Furthermore, in those datasets that are available and seemingly suitable for comparative analyses (including those considered in our study), a lack of concordance and consistency in data collection and reports impedes their consolidation and applicability for a large-scale synthesis, complicating interpretation of data.These hurdles only heighten the need for sourcing and unpacking historical information, foundational in ecology research (Swetnam et al., 1999) to develop our understanding of community baselines and invasionrelated processes.
Other common challenges in bioinvasion ecology that can impede comparative functional traits analyses are (i) precise taxonomic identification or taxonomic bias (Ojaveer et al., 2021) and (ii) correctly assigning biogeographical status for all component species in a community (Carlton, 1996(Carlton, , 2009;;Marchini & Cardeccia, 2017).
Attention is needed to address changes in species nomenclature and designation that may have occurred between two sampling events separated by several decades.Although it is not uncommon for species within the same genus or family to share similar trait modalities (Grabowski et al., 2007), sometimes a misidentified species can bring a mistaken function in the dataset, especially within the traits related to trophic behaviour, environmental tolerance and life history, or when a NIS is incorrectly identified as a native species (Costello et al., 2021;Marchini & Cardeccia, 2017).Another impediment comes from our still limited knowledge of the natural history of many marine invertebrates and thus the difficulty in defining trait profiles and their shift across life-cycle stages for many marine organisms (Cardeccia et al., 2018), which limits both the accuracy and precision of functional traits analyses.It would be valuable to establish and estimate the uncertainty parameters of traits data (e.g. depending on the information source) if assessment of functional changes is implemented in routine bioinvasion management practices.This would also allow identifying the critical knowledge gaps across traits categories and taxonomic groups.
As well as addressing these challenges in near-future research, it is worth asking how to avoid these issues going forward such that mid-and late-century evaluations of current 'baselines' can proceed with fewer problems.Much more comprehensive surveys of coastal marine systems starting now might help interpreting changes associated with multiple anthropogenically driven pressures (including NIS introductions) as well as their interactions and identifying tipping points in the future.To better disentangle pressure-response relations, establishment of harmonised, ecosystem-based monitoring programmes is crucial (Lehtiniemi et al., 2015;Ojaveer et al., 2021).Comprehensive long-term observational and experimental datasets would allow robust quantification of NIS impacts, accounting for interannual and spatial variation, as well as other environmental and anthropogenic covariates (Cleland et al., 2004;Cusser et al., 2021).Greater transparency in dataset publication is already being implemented and will greatly improve access to raw datasets.Wider implementation of paired morphological and genetic identification methods will also facilitate these investigations, along with globally harmonised species nomenclature linked to curated biological traits databases.Similarly, the concept cryptogenic species (Carlton, 1996) is well established and can help avoid mischaracterisations of the native and non-native portions of communities when reconstructing invasion histories.Finally, a stronger awareness of shifting baselines and a determination to avoid overlooking base states has emerged recently in response to Anthropocene impacts and strengthened the field of historical ecology as a sub-discipline in its own right (Kittinger et al., 2015).
Each of these components can give rise to more widespread generation of reliable and comparable large-scale community datasets globally (Lehtiniemi et al., 2015).This will undoubtedly facilitate our understanding of variation for the long-term consequences of global bioinvasion pressure at the ambient scale of regional marine ecosystems, and promote better management frameworks to address it (Ojaveer et al., 2021;Ruiz & Hewitt, 2002).Ultimately, this will help quantify the scale of change we might expect if we fail to develop effective barriers to the continued human-assisted spread of marine species.

ACK N OWLED G EM ENTS
This work is dedicated to the memory of our late colleague and co-author Professor Susan L. Williams, who was tragically killed by a negligent driver in a traffic accident on 24 April 2018.Susan played an important role in the development of the objectives and conceptual framework of the analyses presented in this paper.
We thank Georgia Thomson-Laing and Witold Ming (Cawthron) for their help with quality assurance of the assembled datasets.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors declare no conflict of interest for this study and the submitted manuscript.

PEER R E V I E W
The peer review history for this article is available at https:// www.webof scien ce.com/ api/ gatew ay/ wos/ peer-review/ 10. 1111/ ddi.
Habitat modification cCanopy-providing floating substrate by their living and dead tissues Matrix-forming-provide seafloor substrate by their living and dead tissues Substrate-modifying-modify physical/chemical properties of the habitat Temperature tolerance Wide range-tolerates wide range of temperatures Narrow range-tolerates temperatures typical for one climatic zone Longevity Short-lived-<2 years Long-lived-≥2 years a Here we referred to the maximum body size of adult individuals.

TA B L E 3
Components of functional diversity and functional diversity indices considered in this study to test for changes in community assembly along temporal and spatial gradients in relation to non-indigenous species effects.

F
Averaged relative abundance (% of total community) of NIS across the analysed datasets.The error bars indicate standard deviation and asterisks highlight a significant difference as per Kruskal-Wallis tests.TA B L E 4 Results of PERMANOVA for species data from six considered global coastal regions: British Columbia (BC), New South Wales Australia (AU), Ilha Grande Bay (Brazil), San Francisco Bay USA (SF), North-Eastern Baltic Sea (BS), Waitematā Harbour New Zealand (NZ).

TA B L E 6
Trait modalities determined as significant (p < 0.05) indicators of either low-pressure (blue cells) or high-pressure datasets (red cells) from six regions: British Columbia (BC), New South Wales Australia (AU), Ilha Grande Bay Brazil (BR), San Francisco Bay USA (SF), North-Eastern Baltic Sea (BS), Waitematā Harbour New Zealand (NZ).
Thanks to Beatriz Fleury and Projeto Coral-Sol for providing us with partially unpublished monitoring data from Brazil for use in this study.JCC acknowledges financial support from Fundação Carlos Chagas Filho de Amparo à Pesquisa do Estado do Rio de Janeiro (JCC, FAPERJ E26/201.286/2014); and Conselho Nacional de Desenvolvimento Científico e Tecnológico (JCC, CNPq-307117/2014-6).We also thank Professors Chad L. Hewitt (Murdoch University) and Marnie L. Campbell (Deakin University) for contributions to initial discussions and the study idea development.This research was funded by New Zealand Ministry for Business, Innovation and Employment (MBIE) research grant C01X1511 and internal funds contributed by the authors' organisations.JK and HO acknowledge financial support from the EEA grant 'Climate Change Mitigation and Adaptation' call I 'Ecosystem resilience increased' project 'Impacts of invasive alien species and climate change on marine ecosystems in Estonia'.Contributions of AZ, OF, ID and GAH were also supported by New Zealand Ministry of Business, Innovation and Employment funding (CAWX1904-A toolbox to underpin and enable tomorrow's marine biosecurity system).

Table 3 )
: mean values per region [British Columbia (BC), New South Wales Australia (AU), Ilha Grande Bay Brazil (BR), San Francisco Bay USA (SF), North-Eastern Baltic Sea (BS), Waitematā Harbour New Zealand (NZ)] and dataset (low-pressure vs. high-pressure) with error bars representing standard deviations.Asterisks indicate statistically significant change.