The Introduced Fanworm, Sabella spallanzanii, Alters Soft Sediment Macrofauna and Bacterial Communities

The Mediterranean fanworm, Sabella spallanzanii, is an introduced and established “unwanted species” in New Zealand, subject to nationwide targeted surveillance in port, marina, urban and natural environments. S. spallanzanii has the potential to change soft-sediment benthic habitats due to the physical presence of the fanworm’s tube and associated biological activities, particularly suspension feeding and bio-deposition. A 6-month field experiment was conducted to investigate the impacts of S. spallanzanii on existing communities within invaded soft-sediment habitats. Macrofaunal communities were assessed using traditional sampling and taxonomy while microbial and eukaryotic communities were characterised using metabarcoding of 16S and 18S ribosomal genes, respectively. Live and mimic S. spallanzanii were transplanted at different densities (10 - 50 individuals per m2) into experimental plots with existing assemblages, to test for potential biological and/or physical effects on benthic communities. Analyses revealed consistent, significant differences in macrofaunal, eukaryote and bacterial assemblages in the presence of live S. spallanzanii and mimics, indicating that these effects are brought about by biological and physical functions associated with the worms. The presence of S. spallanzanii did not alter total abundance and taxa richness of benthic assemblages but resulted in compositional differences. Changes in the structure of native benthic communities, as indicate by this study, can potentially impact functioning of soft-sediment habitats, through alterations to nutrient cycling, bioturbation and benthic-pelagic coupling. Quantitative measurements of impacts are crucial to understand the trajectory of marine invasions, their roles in re-structuring communities, and to guide management efforts.


INTRODUCTION
Non-indigenous species (NIS) can have profound impacts on marine coastal ecosystem functioning, by altering community structure, native species richness and ecological processes (Ruiz et al., 1999;Molnar et al., 2008). The magnitude and extent of these impacts vary across temporal and spatial scales (Carlsson et al., 2009), but can be extensive and irreversible (Edelist et al., 2013). Some NIS can regulate the availability of resources to other species by changing the physical state of the ecosystems they invade. Such species are termed "ecosystem engineers" (Jones et al., 1994) or "transformier species" (Richardson et al., 2000), which have the potential to disproportionately affect the functioning of coastal habitats (Cuddington and Hastings, 2004). Ecosystem engineers modulate invaded environments through physical and biological mechanisms. Physical mechanisms are associated with the creation of habitats and modification of the local environment by the presence of body structures of the engineers themselves. For example, introduced mussels or oysters are well-recognized ecosystem engineers that form dense beds or reefs, which can increase local diversity by providing novel hard habitat in soft-sediment environments for many species to colonize (Sousa et al., 2009). On the other hand, biological mechanisms alter the flux of resources through mechanical or chemical functions, such as filtration and bio-deposition that can shift the flow of energy between pelagic and benthic habitats. For example, the proliferation of an invasive filter-feeding clam in San Francisco Bay (USA) caused persistent changes in planktonic food webs and benthic macrofaunal community structure through alteration of benthic-pelagic coupling (Peterson, 2002;Cloern and Jassby, 2012).
The Mediterranean fan worm, Sabella spallanzanii (Gmelin, 1791), a large sabellid tube worm, is considered an ecosystem engineer owing to its capacity to form dense canopies of tubes extending to 50 cm above the sediment. The physical structure provided by the leathery tubes and spiral feeding fans of S. spallanzanii canopies can alter near-seabed hydrodynamics and reduce light penetration below. In Australia, changes in hydrodynamics caused by dense aggregations of invasive S. spallanzanii have been linked to alterations in understory invertebrate larval abundance and recruitment patterns of epifaunal assemblages (Holloway and Keough, 2002a,b). Sabella spallanzanii canopies can also provide structurally complex habitat for other species to inhabit and seek refuge from predators, thus potentially changing the structure of benthic assemblages. Additionally, S. spallanzanii is a suspension feeder that can filter large quantities of water (ca. 12 m 3 per day, Stabili et al., 2006). Its efficient removal of organic matter from the water column, and conversion to fecal material and other waste products, may alter bentho-pelagic coupling in invaded habitats. Previous studies conducted in south-eastern Australia revealed changes in benthic community structure and composition associated with S. spallanzanii (O'Brien et al., 2006;Ross et al., 2007). The impact of NIS is frequently context dependent (Didham et al., 2007;Thomsen et al., 2011), making it difficult to extrapolate findings from one location to another. This merits further efforts to understand the generality of their effects or to add context across regions that may attempt to manage their incursion.
Sabella spallanzanii is native to the Mediterranean and east Atlantic coast of Europe and has established introduced populations along the southern coast of Australia and in New Zealand, and has also been recorded in Brazil (Currie et al., 2000). Sabella spallanzanii was first detected in New Zealand in 2008 (Read et al., 2011) and has since become established in several regions around the country. It occurs subtidally at a range of densities in multiple geographic locations (100s of km apart) including on man-made structures and in natural reef and soft-sediment habitats. Although S. spallanzanii is an epifaunal species that requires a solid substrate on which to attach, it can successfully invade soft sediments by using hard substrate within the sediment as anchoring points, such as a fragment of cobble, shell or debris. The density of an introduced species often drives the magnitude of impact. Worm density is generally higher in artificial than natural habitats, and in subtidal soft sediment environments it generally increases with depth, up to a maximum of 30 m. In these natural soft sediment habitats worms typically grow in clumps (50-200 worms) attached to dead shells or live mussels and ascidians (Currie et al., 2000). Its reproductive biology, including high fecundity rates, "spermcasting" fertilization strategy, and extended reproductive season (Giangrande et al., 2000), contribute to it being found in dense aggregations as a numerically dominant, habitat-forming species (Currie et al., 2000). It has been designated an Unwanted Organism under New Zealand's Biosecurity Act 1993 and is subject to targeted surveillance in harbors, ports, marinas, and high-value natural environments around the country. However, there is still relatively scarce knowledge about the associated adverse effects of established populations. Given S. spallanzanii's wide environmental tolerance and its rapid spread potential into natural habitats, there is a need for a better understanding of its impacts on ecosystem diversity and functioning. This information is critical to understanding resource allocation and benefits of management.
This study aimed to elucidate the ecological impacts of S. spallanzanii on soft-sediment ecosystems in New Zealand. We used a manipulative field experiment to quantify density dependent effects of S. spallanzanii on the structure of soft sediment macrofaunal, eukaryote, and microbial assemblages. We manipulated the density of live and mimic S. spallanzanii to separate biological and physical effects of the worms on macrofaunal and microbial assemblages.

Experimental Design, Sample Collection, and Processing
A field experiment was conducted from September 2017 until March 2018 in subtidal soft-sediment habitats along the Rangitoto Channel, Waitemata Harbor, Auckland, New Zealand (−36.813, 174.8396). The site had a water depth of 8 m and consisted of relatively flat soft-sediment habitat. The site is situated just outside Waitemata Harbor and is located approximately 1 km outside the commercial shipping lane to the Port of Auckland. An earlier survey revealed the presence of S. spallanzanii in this area, at varying densities (0-20 individuals per m 2 ).
A two-factor experimental design was used to examine the influence of fanworms on soft sediments at the study location. A first factor, "Treatment, " had three levels: "Sabella, " corresponding to the addition (via transplants, see below) of live S. spallanzanii worms to experimental plots; "Mimic" corresponding to the addition of mimic S. spallanzanii worms; and "Control" corresponding to plots without worms. Live S. spallanzanii were obtained by divers from adjacent soft sediment invaded areas and transplanted into the plots using garden pegs. During a 2017 pilot study conducted at a nearby marina, the use of U-shaped steel garden pegs had been identified as a reliable method to anchoring and ensuring the viability of transplanted live and mimic worms (Atalah, unpubl. data). Only healthy undamaged worms with intact tubes were transplanted for this treatment. Many of the live worms collected at the study site were attached to shell fragments, sometimes with several individuals per shell, and thus treatments were created by transplanting naturally occurring clusters of worms into experimental plots to achieve the target densities. The average tube length of transplanted worms was 32.1 ± 0.66 cm (mean ± standard error), with a range of 10-57 cm. Sabella spallanzanii mimics was made from elastic bungee cord (of a length and diameter similar to real FIGURE 1 | Schematic diagram of the design and set-up of the field experiment showing the random allocation of treatments to 1 m 2 experimental plots: live Sabella spallanzanii worms (S, yellow squares), mimic S. spallanzanii worms (M, green squares) and controls (Ctrl, purple squares) without worms. Live and mimic worms were transplanted at four different densities (i.e., 10, 20, 35, and 50 worms (or mimics) per plot). Additionally, eight ambient samples (A, blue circles) were collected from adjacent uninvaded areas as procedural controls.
S. spallanzanii) with frayed ends that resembles the worms' crown structure. These were threaded onto pegs in clumps of two, three or five individuals, similar to the natural aggregations observed at that site (Supplementary Video 1).
A second factor, "Density, " had five levels: zero (control), 10, 20, 35, and 50 worms (or mimics) per plot. Twenty-four 1 m 2 plots, separated by approximately 2 m, were arranged along two perpendicular 30 m transect with a common vertex (Figure 1, n = 12 plots per transect). Plots were positioned on each side of the transects, separated by a corridor to allow divers access. Four replicate plots were assigned to Controls, two replicate plots were assigned to 10, 20, and 35 density level plots (both Sabella and Mimic), and four replicates to the 50 individual density level (Figure 1). Extra replicates in the end member treatments (i.e., four replicates each in the control and 50 density treatments, rather than two each) limited the influence of potential outliers and provided greater statistical power. To ensure interspersion, experimental plots were randomly assigned to treatment combinations. Control plots were treated in the same manner as Sabella and Mimic plots (i.e., garden pegs inserted into the sediment), but no worms were transplanted.
Experimental plots were established on 19th September 2018 and left undisturbed until 21st March 2019, when destructive sampling of macrofaunal assemblages and sediment characteristics was carried out to conclude the experiment. Two small sediment cores (3 cm internal diameter, 2 cm deep) were collected from each experimental plot: one for sediment grain size and organic matter content analysis (Gatehouse, 1971;Mook and Hoskin, 1982), the other for sediment chlorophylla and phaeopigment content (Sartory and Grobbelaar, 1984). One core (13 cm internal diameter, 15 cm deep) was collected from a randomly selected area within each of the 24 plots and the material retained on a 500 µm mesh sieve preserved in 70% isopropyl alcohol for macrofaunal community analysis. Eight additional set of samples were collected from noninvaded adjacent "Ambient" areas (i.e., without S. spallanzanii) to serve as procedural controls. Macrofauna were separated from sediment and shell hash after staining with Rose Bengal and identified to the lowest feasible taxonomic level and their abundance quantified using a binocular dissecting microscope. A random subset of samples was quality checked by an independent taxonomy expert. All samples and vouchers specimens of the identifications were kept for later verification. The macrofauna dataset was dominated by infaunal taxa, but by nature of the sampling method also encompassed surfacedwelling organisms (e.g., epifauna growing on shell hash or coarser sediment fractions). All organisms were included in the analyses, as we were interested in changes in all components of soft-sediment communities. Molecular samples were collected in parallel with macrofauna samples to enable direct comparison of the metabarcoding results (i.e., microbial and eukaryotic communities) with the traditional macrofaunal methods. Specifically, triplicate 50 g sediment samples were obtained by scraping the top 2 cm layer from haphazardly selected intact areas within each of the 24 plots using sterilized plastic vials. Molecular samples were frozen and transported to the laboratory and immediately stored at −80 • C until further processing.

High-Throughput Sequencing and Bioinformatics Analyses
The microbial and eukaryotic community compositions of the triplicate sediment samples collected at the 24 plots, represented by 10 Mimic plots (n = 30 samples), 10 Sabella plots (n = 30 samples), 4 control plots (n = 12 samples), and 7 ambient plots (n = 21 samples), were assessed using DNA metabarcoding and following a similar high-throughput sequencing library preparation method as outlined in Keeley et al. (2018). Briefly, environmental DNA was extracted from pre-homogenized 0.2 g of surface sediment using the DNeasy PowerSoil TM DNA extraction kit (Qiagen, Hilden, Germany), including DNA extraction blanks (i.e., extraction controls), followed by DNA quality and purity verification using a spectrophotometer (Eppendorf, Leipzig, Germany). A two-step tailed Polymerase Chain Reaction (PCR) amplicon procedure (Kozich et al., 2013) was applied to all DNA samples (n = 95; including pooled extraction and pooled PCR blanks), targeting ca. 400 base-pair (bp) fragments of a) the nuclear 16S rRNA bacterial gene (V3-4 region) and b) the nuclear 18S rRNA eukaryotic gene (V4 region), and using the primers, reagents and PCR thermocycling conditions described in von Ammon et al. (2018). Detailed information on sequence library preparation, demultiplexing and follow-up bioinformatics analyses are provided in Supplemental Material 2.

Data Analyses
All three datasets (i.e., macrofauna, eukaryotes and bacteria) were analyzed in a similar manner using R software (R Core Team, 2018). Amplicon sequence variants (ASVs) of bacteria and eukaryotes were rarefied at a sequencing depth of 5,000 reads per sample using the "rarefy_even_depth" function in the phyloseq package (McMurdie and Holmes, 2013). This resulted in two and one triplicate samples removed for bacteria and eukaryotes, respectively, because they contained fewer reads than the sample size. Additionally, ASVs represented with a mean < 1 −4 and present <3 times in at least 20% of the samples in each dataset were removed. Rarefaction curves were calculated for each sample using the vegan package (Oksanen et al., 2018). Differences in assemblage structure between treatments were tested using permutational analysis of variance (PERMANOVA, Anderson et al., 2002) based on Bray-Curtis similarities of the log (x + 1) transformed abundance data using the "Adonis" function in the vegan package. "Treatment" was considered a fixed factor and "Density" included as a continuous covariate representing the final density of fanworms per plot. Significant terms were then investigated using pair-wise comparisons with the PERMANOVA t statistic and 999 permutations. Assemblage differences among treatment levels were visualized using canonical analysis of principal coordinates (CAP, Anderson and Willis, 2003). Indicator taxa for each treatment were determined using the R package "indicspecies" (Cáceres and Legendre, 2009). Indicator species were selected if they were significantly associated with each treatment and their relative abundance visually compared among treatments. Univariate analyses of covariances (ANCOVA) were used to test changes of taxa richness of all three datasets and macrofaunal total abundance using the same experimental design described above for the multivariate data. Tukey tests with Bonferroni adjustment for multiple comparisons were used to compare means when factor "Treatment" was found to be statistically significant in the ANCOVA analyses. In all cases the "Treatment × Density" interaction was not significant (P > 0.3), thus it was removed from the ANCOVA models. Model assumptions were checked by inspecting the model residuals. Type III sums of squares (and its direct multivariate analog) were used to analyse the unbalanced designs.
The relationship between species data and sediment environmental variables was analyzed using multivariate multiple regression (McArdle and Anderson, 2001), more specifically using the adonis routine in the vegan package. A conditional test was used where individual variables were fitted separately to test their relationship with the assemblage data (ignoring other variables) based on Bray-Curtis dissimilarities of the log (x + 1) transformed abundance data. Models were visualized using redundancy analysis ordination biplots (Borcard et al., 2011).

Macrofaunal, Eukaryote, and Bacterial Diversity
Overall, the macrofaunal data comprised 75 different taxa, classified into 10 different phyla, 13 classes, 22 orders, 45 families, and 68 genera. Species level identifications were made for 27 of the 75 unique taxa. Macrofaunal assemblages were largely dominated by annelids, molluscs and crustaceans (52, 19, and 16% of all specimens, respectively, Figure 2A), with the non-indigenous bivalve Theora lubrica numerically dominant, followed by polychaetes from the Paraonidae and Cirratulidae families.
The illumina TM MiSeq sequencing run generated a total of 12,717,386 sequences, consisting of 5,592,988 bacterial 16S rRNA gene sequences and 7,124,398 eukaryotic 18S rRNA gene sequences. Following filtering (denoising, merging, and chimera removal), the number of high-quality 16S and 18S sequences retained for taxonomic assignments of Amplicon Sequence Variants (ASVs) were 1,848,310 and 3,111,615, respectively. Following reads rarefaction and the removal of ASVs represented with a mean <1 −4 and present <3 times in at least 20% of the samples, the number of ASVs retained for downstream statistical analyses for the eukaryote and bacterial datasets were 1,210 and 6,947, respectively; with rarefaction of samples indicating that all samples reached saturation (Supplemental Material 3). Eukaryotic ASVs were assigned to 29 different phyla, 64 classes, 137 orders, 232 families, 304 genera, and 339 distinct species. The output of the metabarcoding identification resulted in a somewhat similar abundance-based partitioning of higher taxa ( Figure 2B) with phyla Annelida, Mollusca and Nematoda dominating the eukaryote dataset (26, 21, 20% of all detected specimens, respectively, Figure 2B). Bacterial ASVs were assigned to 44 different phyla, 92 classes, 206 subclasses, 281 families. Bacterial assemblages were largely dominated by Protobacteria, in particular from the family Woeseiaceae and Desulfobulbaceae, followed by bacteria from the phylum Bacteroidetes ( Figure 2C). All sequence data are deposited in GenBank's Short Reads Archive, BioProject PRJNA552692, Accession number SRR9644676 to SRR9644588.

Treatment Effects on Assemblage
We found no differences between control plots and ambient (procedural) sediment samples (i.e., the effects of garden pegs and manipulation) on macrofaunal total abundance (t-test, P > 35), richness (t-test, P > 0.25) or assemblage structure (PERMANOVA, P > 0.05) for all three communities (i.e., macrofauna, eukaryotes, and bacteria, Supplemental Material 4). Therefore, results measured from ambient plots were excluded from all further analyses. There was a significant effect of treatment on assemblage structure consistent across macrofaunal, eukaryote and bacterial datasets (P < 0.05, Table 1), although there was no main or interactive effect of the factor Density for all three datasets (P > 0.1, Table 1). The CAP plots illustrate this pattern, showing a clear separation between assemblages from the three treatments in all three datasets (Figure 3). For the macrofaunal data, pairwise comparisons distinguished between Sabella and Control (P < 0.05); and Sabella and Mimic (P < 0.05), but not between Control and Mimic (P > 0.1). The most significant indicators for Control plots included the brushworm Phylo novazealandiae, amphipods, the burrowing crab Pilumnus novaezealandiae, the burrowing holothurian Taeniogyrus dendyi, and the polychaete Aglaophamus sp. Significant indicator taxa for Mimic plots included the bivalves Monia zelandica and Limaria africana, the polychaetes Armandia maculata, unidentified worms from the family Spionidae and Heteromastus filiformis ( Figure 4A). Sabella plots were associated with polychaete family Paraonidae, nemerteans and the polychaete Prionospio sp. (Figure 4A).
For the eukaryote communities, pair-wise comparisons found significant differences among all three treatments (P < 0.05). Controls were significantly associated with higher abundances of polychaetes from the family Lumbrineridae and nematodes from the family Microlaimidae ( Figure 4B). Most significant indicator eukaryote taxa for the Mimic treatment included polychaetes from the families Paraonidae and the non-indigenous Serpulidae tubeworm Hydroides elegans and ascidian Molgula manhattensis, whereas tellinid bivalves, harpacticoid copepods and dorvilleid polychaetes were more abundant on Sabella plots ( Figure 4B).
For bacterial assemblages, pair-wise comparisons also revealed significant differences among all three treatments (P < 0.01 for Control vs. Mimic and Mimic vs. Sabella, and P < 0.05 for Control vs. Sabella) and 55 indicator taxa were identified as being significantly associated with these differences. Most significant indicator taxa for the Control treatment included families in the phyla Actinobacteria, Bacteroidetes and Proteobacteria (Figure 4C), whereas for the Mimic treatment they included bacteria in the phyla Acidobacteria, Gemmatimonadetes, Nitrospirae, and Proteobacteria ( Figure 4C). Significant indicator bacterial taxa for the Sabella treatment most notably included families in the phyla Kiritimatiellaeota (e.g., Kiritimatiellaceae), Spirochaetes (e.g., Spirochaetaceae) and Proteobacteria (Figure 4C), including significant higher abundance of Desulfobacteraceae.
Total abundance and richness of macrofaunal and bacterial assemblages did not vary significantly among Treatment or Density (Figure 5 and Table 2), although there was a nonsignificant trend for higher numbers of macrofauna taxa richness 1 | Results of PERMANOVA analyses to test the effect of Treatment ("Sabella," "Mimic," and "Control") and Sabella spallanzanii Density on macrofaunal, eukaryote and bacterial assemblages.

Macrofauna
Eukaryote Bacteria  in the Mimic compared to the Control plots. Treatment had a significant effect on eukaryote taxa richness, with Mimic having higher number of taxa compared to Control plots (P < 0.05, Figure 5 and Table 2).

Relationship Between Soft-Sediment Assemblages and Environmental Variables
The sediment at the site was relatively poor in terms of organic content (mean 2.5% LOI ± 0.1 S.E.) and was dominated by very fine sand (63-500 µm, 46.6% ± 0.9), silt particles (< 63 µm, 29.3% ± 1.1) and clay (11% ± 1.1) grain size fractions. Mean chlorophyll-a and phaeopigments contents were on average 1.92 µg g −1 (± 0.1) and 4.5 µg g −1 (± 0.3), respectively. Sabella treated plots had higher fine sand, chlorophyll-a, organic and phaeopigments compared to Mimic and Control plots (Figure 6). On the other hand, Mimic and Sabella plots had higher content of medium and coarse sediment in relation to the controls (Figure 6). The multivariate multiple regression analysis showed that gravel, medium and fine sediment had a significant relationship with macrofaunal assemblage structure, explaining 8, 7, and 9% of the variation in the data, respectively.
Eukaryote assemblage structure was significantly related to sediment chlorophyll-a concentration and marginally to the fine sediment fraction and total organic content; whereas bacterial assemblage structure was not significantly related to sediment characteristics. Although no significant relationships were detected between bacterial assemblages and sediment characteristics, there was a trend of bacterial associations with the fine sediment fraction and phaeopigment concentration (P < 0.1). The dbRDA plots (Figure 7) illustrate that all three benthic assemblages were distinguishable among treatments constrained to the variation in sediment characteristics. In all three cases Sabella plots assemblages were generally well separated from the other treatments and correlated with higher phaeopigment, chlorophyll-a, organics and the fine sediment fraction, whereas Mimic plots were related to larger sediment fractions (i.e., medium and gravel).

DISCUSSION
Soft-sediment communities in Waitemata differed significantly in the presence and absence of Mediterranean fanworms     and their mimics. This effect was consistent across benthic macrofaunal, eukaryote and bacterial assemblages, and indicates that S. spallanzanii may affect benthic community structure through both biological and physical processes and mechanisms. The measured differences among treatments were based on relative abundances of several taxa between plots with transplanted live and mimic Sabella and control plots rather than differences in overall abundance or richness. These differences in community structure can have important consequences for ecosystem functioning of soft-sediment habitats (Thrush et al., 2017) based on shifts in feeding, burrowing and sediment irrigation intensity that can drive crucial functions, such as nutrient and organic matter processing. Small mobile taxa (mainly amphipods) and larger bioturbator taxa such as ophiuroids, a burrowing crab (Pilumnus novaezealandiae), and a burrowing holothurian (Taeniogyrus dendyi) were less abundant, or absent, in treatment plots compared to controls. This may be a direct or indirect consequence of both physical presence and biological activities of fanworms. Resident communities may be affected by small-scale changes in hydrodynamics linked to the physical structures of S. spallanzanii canopies, as changes in water velocities can alter sediment stability and concentrations of key solutes (e.g., oxygen, sulfide; O'Brien et al., 2006). Mobile bioturbator taxa can substantially change the physical and chemical properties of soft sediment habitats and thus are crucial for the functioning of these ecosystems (Lohrer et al., 2004). On the other hand, S. spallanzanii plots in our study were also associated with significantly higher abundance (compared to controls) of direct deposit-feeders, such as polychaetes from the family Paraonidae, Dorvilleidae and Spionidae, and Nematoda. The biodeposition activity of the fanworm may has the potential to change nutrient cycling and the microbial community, with knock-on effects on deposit feeders. Suspension feeders, such as S. spallanzanii can increase sedimentation of organic material and produce abundant feces and pseudofeces rich in organic matter, which may serve as food sources for deposit feeding organisms (Commito and Boncavage, 1989). This biological activity may explain the observed patterns of higher concentration of chlorophyll-a, organics and phaeopigments in Sabella treated plots. As such, there has been some concern about the effect of S. spallanzanii on nitrogen cycling including changes in nutrient availability and its consequences for community structure (Harris et al., 1996).
Changes in bacterial communities in response to experimental treatments were consistent with those observed for macrofaunal and eukaryote communities. Soft-sediment dwelling bacteria play a crucial role in the functioning of these habitats and are recognized for their sensitivity to environmental change. Recent studies have highlighted the ability of S. spallanzanii to filter, accumulate and concentrate bacteria from the surrounding environment (Stabili et al., 2006;Licciano et al., 2007). Sabella spallanzanii have also been shown to exude mucus with antibacterial properties, which may influence the composition of microbial biofilms on the worm's tube and surrounding sediments (Stabili et al., 2011). In contrast, nutrient enrichment and increased primary production observed in Sabella treated plots may decrease the redox potential in sediments and favor sulfide-oxidizing and sulfate-reducing bacteria for example, as illustrated by the observed increase in the abundance of Desulfobacteraceae. Regardless of the mechanisms that may be responsible for these effects, the observed changes in bacterial communities may have important consequences for the functioning of soft-sediment habitats, most notably on nutrient and organic matter cycling.
It is often assumed that the impacts of invasive species are density-dependent, although few studies have assessed this contention (Thomsen et al., 2011) and this may not always be the case. We found no effect of density in our simulated invasion of such an engineer throughout our study, which may indicate that no real effect occurred, that the range of experimental densities was too small to have an effect, or that the duration of the study was insufficient for density-based differences to emerge. Non-indigenous filterfeeder species, such as S. spallanzanii, are more likely to affect communities over longer time scales than provide instantaneous shifts or impacts (Byrnes and Stachowicz, 2009). Nonetheless, the simple presence of fanworms or mimics (≥10 m −2 ) appeared to prompt community differences compared to plots without them. One mechanism is through secondary colonization of S. spallanzanii tubes that can provide additional habitat for a range of epifauna typically associated with hard substrates (O'Brien et al., 2006). Although no significant changes in total macrofaunal abundance were detected, the presence of S. spallanzanii and mimics was associated with higher eukaryote taxa richness and higher abundance of encrusting and non-indigenous organisms such as tube worms (Hydroides elegans), ascidians and hydroids, evidenced by the macrofaunal and eukaryote datasets. These changes were reflected in a non-significant trend for an increase in macrofauna taxa richness in mimic plots compared to the controls. Although metabarcoding identification resulted in similar taxonomic partitioning, estimated taxa richness was an order of magnitude higher than that using morphological identification. Metabarcoding provides a more holistic view of the metazoan taxonomic diversity, regardless of the size and developmental stage (e.g., including eggs or larvae). The eukaryote data not only included macrofaunal taxa that dominated morphological samples, but also small-sized meiofaunal species (<1 mm), extending the scope of the analysis. Novel epifaunal taxa facilitated by the presence of S. spallanzanii may have important implications in the functioning of invaded soft-sediment habitats, as the combined filtering capacity of these organisms is considerably enhanced compared to uninvaded habitats (Lemmens et al., 1996). Furthermore, positive and potentially synergistic interactions among invasive species may lead to accelerated impacts on native ecosystems (Simberloff and Von Holle, 1999).
The results reported here are generally consistent with previous studies looking at the impacts of S. spallanzanii in soft sediments assemblages (O'Brien et al., 2006;Ross et al., 2007). However, the effects recorded in these Australian studies from Port Phillip Bay were dependent on spatial distribution and density of S. spallanzanii, with negligible effects recorded at low densities (<1 individual per m 2 ). However, similarly to this study, localized dense clumps (15 individuals per m 2 ), caused changes in benthic infauna structure characterized by lower abundances of cumaceans, ostracods and harpacticoid copepods (O'Brien et al., 2006). Ross et al. (2007) found no effect of S. spallanzanii on the resident macrofauna, with the exception of lumbrinerid polychaetes and gammarid amphipods; however, these taxa only represented a small proportion of those present. Ross et al. (2013) found no difference in the total abundance of macrofauna between S. spallanzanii and control plots. However, the composition of assemblages did change significantly in the presence of S. spallanzanii, with a significant increase in the abundance of echinoderms (largely attributable to brittle stars) in the presence of S. spallanzanii at all three sites investigated (Ross et al., 2013).
Our study site experiences relatively high velocity tidal currents (average 1.25 m s −1 , Greig and Proctor, 1988), which may pose physical stress on organisms and are unfavorable for the accumulation of organic matter. It is likely that these prevailing hydrographic conditions and resulting sedimentary physical conditions may exert a stronger meso-scale effect on biological communities compared to the relatively subtle observed response to the experimental treatments. In this context, the generality of observed effects to a regional scale or other locations is difficult to gauge. On the other hand, conditions that promote a mixed particle size benthic surface (larger fragments within a sandy-silt matrix) and regular water flow may be ideal for S. spallanzanii settlement and growth. Sabella spallanzanii is also among the only or very few organisms in this area's benthic community that can extend as high as 50 cm above the seafloor into the water column to take advantage of passing food-laden currents.
A key factor for future work is to monitor population density, especially clusters of high density, because the magnitude of effects would undoubtedly increase if such novel canopies of this species proliferated.
Shallow subtidal soft sediment habitats are an extensive component of bays and estuaries and provide an important environment for S. spallanzanii establishment and spread. This species represents a relatively novel introduction as an epibenthic (rather than infaunal) invader that is not (almost) exclusively associated with hard-bottom habitats, like a majority of marine introduced species in New Zealand (Cranfield et al., 1998) and elsewhere (Byrnes et al., 2007;Ruiz et al., 2011). Given that these habitats are a dominant and important feature in the study area and wider region, it would be valuable to replicate this study in other areas to better understand the generality of effects of this invasion. As it is a listed "unwanted" organism under New Zealand legislation, there are active campaigns to (a) prevent its anthropogenic spread (mainly via shipping and boating) and (b) locally extirpate any nascent populations that are detected in discrete new areas using targeted removals (Read et al., 2011). These efforts, combined with monitoring hotspot populations that are not being managed, provide a precautionary approach to S. spallanzanii in New Zealand (Champion, 2018). An eradication attempt in the Auckland region is considered too uncertain and potentially expensive to initiate at present using existing tools and strategies. This view may change as marine eradication matures and if impacts of NIS on benthic communities, functioning, and maritime industries (e.g., aquaculture) proliferate (Edwards and Leung, 2009;Soliman and Inglis, 2018).
NIS are generally considered among one of the greatest threats to marine ecosystems, however our understanding about their overall impacts is often insufficient (Thomsen et al., 2011), partly because of the lack of adequate studies and context-dependency of their effects. Elucidating and predicting context-dependent impacts of NIS are crucial for assessing the associated risks to biodiversity and ecosystem functioning, and to underpin management efforts. At a global scale it has been shown that NIS have significant, although modest in magnitude, ecological effects (Anton et al., 2019). Using a combination of morphological and molecular taxonomical techniques, our study provides insights into the nature and magnitude of local-scale effects on different communities inhabiting soft-sediment habitats, which may inform the implementation of effective knowledge-based management and mitigation measures that reduce the impacts of NIS.

DATA AVAILABILITY STATEMENT
All data and R scripts used for the analyses are available at https:// github.com/jatalah/sabella_impacts_experiment.git.

FUNDING
This study was funded by the New Zealand Ministry of Business, Innovation and Employment, under the Programme What's at stake?-Enabling decision-making through better measurement, forecasting, and evaluation of the impacts of nonnative organisms in NZ's changing ocean (C01X1511).

ACKNOWLEDGMENTS
We thank Isaac Eckert and Fiona Gower (Cawthron Institute) for helping with sample processing, Samantha Happy from Auckland Council for helping with pilot study, Westhaven Marina for access to the study sites and Dr. Ian Davidson for his valuable insights comments on this manuscript. We also acknowledge MPI-New Zealand Biosecurity (Kathy Walls) for processing permits for the transplanting of Sabella spallanzanii.

SUPPLEMENTARY MATERIAL
The Supplementary Material for this article can be found online at: https://www.frontiersin.org/articles/10.3389/fevo.

2019.00481/full#supplementary-material
Supplementary Video 1 | Video footage of the field set-up of the experiment in conducted in Waitemata Harbour, Auckland, New Zealand showing the treatments to 1 m 2 experimental plots: live Sabella spallanzanii worms, mimic S. spallanzanii worms and controls without worms. Live and mimic worms were transplanted at four different densities [i.e. 10, 20, 35 and 50 worms (or mimics) per plot].