Four different types of PBS-based polymer sheets (15 cm × 15 cm) were manufactured by SCION™ (Rotorua, New-Zealand) and used as ‘polymer type’ treatments in the experiment ( Figure 1A ): (i) OYS20 – PBS containing 20% of oyster shell content; (ii) OYS40 – PBS containing 40% of oyster shell; (iii) OYS60 - PBS containing 60% of oyster shell; (iv) Control-PBS only, with no additives added. Polybutylene succinate (PBS, PTT MCC biochemical, grade FD92PM) was purchased from Chemiplas NZ (Auckland, New Zealand). Pellets were dried at 60°C overnight before use. Oyster shells were ground and sieved. Only the fraction between 125 and 250 µm was kept. The oyster shell powder was dried at 50°C overnight before use.

The extrusion process consisted of:

For deployment, three PVC structures (120 cm × 120 cm) were fabricated, each comprising polymer sheets attached to PVC pipes ( Figure 1A ). The sheets were arranged on the structure following a Latin square design, which allowed for randomization of the bio-polymer and control sheets across the structure (Gao, 2005). Experimental frames were deployed on a floating pontoon at 0.5 m depth (polymers remained submerged at the same depth throughout the experiment) in the Nelson Marina (Nelson, New Zealand) (S -41.26°; W 173.27°), on the 29^th of January 2019, and the bio-polymer sheets were sampled 1, 3, 6, 12 and 18 weeks after deployment ( Figure 1C ) for environmental DNA (eDNA) analysis. Temperature data loggers (HOBO^® Pendant temperature loggers, Onset, USA) were attached to each of the three deployment structures, and the temperature was recorded every 15 min for the duration of the experiment. The temperature at the study site ranged between 23.3°C in January 2019 and 13.7°C in May 2019 ( Figure S1 ).

At each sampling time, three sheets of each bio-polymer (including the control) were collected (i.e., n = 12 sheets per sampling point). Biofouling isolation was undertaken in a non-destructive manner in the field, wearing sterile gloves for each sample. At each sampling time, a new polymer sheet was collected. On weeks one and three, early biofilm samples were collected using sterilized sponges (Whirl-Pak™, Speci-sponges™, Nasco, WI, USA) and thoroughly rubbing ¼ of the surface of the bio-polymer sheets on the front and back side separately. The sponge and collected material from each side were isolated into individually sealed sterile plastic bags provided by the same manufacturer. In the subsequent weeks (when mature biofouling has been established), sterile stainless-steel surgical blades (Swann-Morton, Sheffield, United Kingdom) were used to carefully detach biofouling from ¼ of the surface from the front and back side of the bio-polymer sheets. Each sample was then isolated into a sterile 50 mL centrifuge tube. Collected samples were transported on ice to the Cawthron Institute laboratory located 5 min away and processed immediately for DNA extraction Figure 1B .

In the laboratory, biofilm samples collected with sponges and blades were processed as described in Zaiko et al. (2016) and von Ammon et al. (2019), respectively. Briefly, RNA/DNA free water (40 ml; Life Technologies, Carlsbad, CA, USA) was added to each plastic bag containing a sponge sample, and the bag was crushed using a Colworth 400 laboratory stomacher press (Seward, Worthing, UK) for 2 minutes at maximum speed. Excess liquid was then squeezed from the sponges, and the resulting suspensions were pelleted by centrifugation (4,000 g, 15 min). Then, the supernatant was discarded, and pellets were used for DNA extraction. The blade harvested samples were homogenized by bead beating for 2 min (1,600 strokes/min; 1600 MiniG Spex SamplePrep), using stainless steel beads (SPEX™ SamplePrep Stainless-Steel Grinding Balls 3/8”, 9.5 mm diameter, Fisher Scientific, New Hampshire, United-States), centrifuged (15,000 g, 5 min) and the supernatant discarded. A single subsample of 0.2 g was taken from the homogenized material for subsequent DNA extraction.

DNA from all biofouling samples (n = 120) was extracted using the PowerSoil^® DNA Isolation Kit (QIAGEN, MOBIO, Carlsbad, USA), following the manufacturer’s protocol. A DNA extraction blank was included for every 40 samples (three blanks in total). The quantity and quality of the extracted DNA were assessed using a NanoPhotometer (Implen, Munich, Germany).

The V3-V4 regions of the bacterial 16S ribosomal RNA (16S rRNA) gene and the V4 region of the eukaryote nuclear 18S ribosomal RNA (18S rRNA) gene were amplified by Polymerase Chain Reaction (PCR), using respectively the bacteria-specific primers 341F: 5’-CCT ACG GGN GGC WGC AG-3’ and 805R: 5’-GAC TAC HVG GGT ATC TAA TCC-3’ (Herlemann et al., 2011; Klindworth et al., 2013) and eukaryotic-specific primers Uni18SF: 5’-AGG GCA AKY CTG GTG CCA GC-3’ and Uni 18SR: 5’ -GRC GGT ATC TRA TCG YCT T-3’ (Zhan et al., 2013). Both primer sets contained an Illumina overhang adapter for NEXTERA indexing, as in Pochon et al. (2019). Polymerase Chain Reaction analyses were done in an Eppendorf Mastercycler (Eppendorf, Hamburg, Germany) in a total volume of 50 µL using MyFi™ PCR Master Mix (sample concentration ranging from 18.9 to 45.82 ng/µL) (Bioline Meridian Bioscience, Memphis, Tennessee, USA), 2 µL of each primer (10 mM stock), and 1 µL of template DNA. The PCR cycles for the 16S rRNA gene amplification were as follows: 95°C for 3 min followed by 33 cycles of 94°C (20 s), 52°C (30 s) and 72°C (30 s), with a final extension at 72°C for 5 min. Amplifications for the 18S rRNA gene were 95°C for 3 min followed by 37 cycles of 94°C (20 s), 54°C (20 s), and 72°C (30 s) with a final extension at 72°C for 7 min. Negative (no-template) PCR controls were included in each PCR run. Amplicon PCR products were purified using AMPure XP PCR Purification beads (Agencourt, Beverly, MA, USA), quantified using a Qubit Fluorometer (Life Technologies, Carlsbad, CA, USA), and diluted to 3 ng/µL. Water control was added to test for potential contamination during the sequencing workflow. Normalized PCR products and controls (n = 133; 120 samples, four extraction blanks, seven PCR blanks, and two water blanks) were sent for library preparation and sequencing on an Illumina MiSeq™ platform at Auckland Genomics, University of Auckland, New Zealand, following the Illumina 16S rRNA metagenomics library preparation manual (D’Amore et al., 2016). Sequencing adapters and sample-specific indices were added to each amplicon via a second round of PCR using a Nextera Index kit. After that, five μl of each indexed sample was pooled, and a single clean-up of pooled PCR products was undertaken. A bioanalyzer was used to check the quality of the library, which was then diluted to 4 nM and denatured. The library was diluted to a final loading concentration of 7 ρM with a 15% spike of PhiX. Paired-end sequences (2 x 250 bp) were generated on an Illumina MiSeq instrument. Raw sequences were deposited in the NCBI short read archive under accession: PRJNA811748.

Raw sequence reads (with Illumina adapter sequences removed by sequencing facility) were trimmed using cutadapt v2.10 to remove primer DNA sequences (Martin, 2011), with no primer mismatch allowed. Then quality filtering, denoising, merging pair-end sequences, and calling amplicon sequence variants [ASVs] were performed using the DADA2 version 1.20.0 package, implemented in R version 4.0.5 (Callahan et al., 2016). Following exploration of the DNA sequence quality plots, the following parameters were applied for quality filtering and denoising for 16S rRNA and 18S rRNA datasets, respectively: “truncLen=c(220,220), maxEE=c(2,6), truncQ=2, ndmaxN=0” and “truncLen=c(150,150), maxEE=c(2,6), truncQ=2, ndmaxN=0”. Singleton ASV data were discarded, and the remaining paired-end reads were merged with a minimum overlap of 25 bp and one mismatch allowed in the overlap region. Chimera removal was performed using the default (consensus) method and the de-noised ASVs were taxonomically classified against the SILVA 132 database for 16S rRNA (Quast et al., 2012) and PR2 database for 18S rRNA (Guillou et al., 2012) using DADA2 “assignTaxonomy” command, based on the RDP classifier (Wang et al., 2007).

The 16S rRNA dataset was filtered to exclude any ASVs classified as Eukaryota in the rank Kingdom, Chloroplast in the rank Class, and Mitochondria in the rank Family using the “subset_taxa” command implemented in R package phyloseq (McMurdie and Holmes, 2013). The maximum number of ASVs found across negative controls was subtracted from the corresponding ASVs to offset potential contamination noise (Bell et al., 2019; Clark et al., 2020). Then, rarefaction curves were plotted for both 16S rRNA and 18S rRNA datasets using the ‘ggrare’ function in R (package ranacapa; Kandlikar et al. (2018)). The 16S rRNA rarefaction curves indicated that the sequencing depth attained per sample adequately captured biodiversity (i.e., the curves have reached a plateau). In the 18S rRNA dataset, four samples yielded an extremely low post-filtering number of reads (<2,000) and were removed from further analyses ( Table S1 ). For calculating diversity metrics (ASV richness and Shannon index), metabarcoding data were normalized by rarefying at equal depth (10,000 and 2,000 reads for 16S rRNA and 18S rRNA data, respectively). The community structure analyses were performed on the unrarefied datasets transformed into proportional read abundance. Amplicon Sequence Variant (ASV) tables generated by our bioinformatic pipeline were uploaded into the Plymouth Routines in Multivariate Ecological Research (PRIMER 7) v7.0.13 software (Anderson, 2001; Clarke and Gorley, 2015). Permutational Multivariate Analysis of Variance (PERMANOVA) was used to investigate whether bacterial and eukaryotic community structure responds differently across experimental factors ‘Type’ (i.e., different bio-polymer concentrations), ‘Sampling’ (i.e., sampling time), and ‘Side’ (i.e., front or back of the bio-polymer sheets). Bray-Curtis dissimilarity matrices were calculated for bacterial and eukaryotic datasets (at ASV level) (Bukin et al., 2019). Bray-Curtis dissimilarity based on relative abundance of sequence reads was chosen as one of the most commonly used in metabarcoding studies (Bukin et al., 2019) and not affected by the high frequency of zero counts typical for HTS datasets (Pawlowski et al., 2016). Two-way permutational ANOVAs (PERMANOVA) with ‘Sampling time’ as a random factor and ‘Type’ as a fixed factor was performed to explore differences in bacterial and eukaryotic composition across experimental factors. Differences P ≤ 0.05 were deemed statistically significant (Bland and Altman, 1995). In addition, tests of homogeneity of dispersions (PERMDISP), were performed to quantify the variability in community structure between treatments for the bacterial and eukaryotic datasets, using "Sampling" as a group factor. 999 permutations and calculating distances to centroids. Variation in community structures among the five sampling times (1 week, three weeks, six weeks, 12 weeks, and 18 weeks) were ordinated using a Principal Coordinate Analysis (PCoA) for both bacterial and eukaryotic datasets to visualize distribution patterns across ‘Types’ and ‘Sampling.’ One-way analysis of variance (ANOVA) was performed on ASV richness and Shannon-Wiener index, separately on the bacterial and eukaryotic dataset for the group factors ‘Sampling’ and ‘Type’ using the ‘Vegan package in RStudio (Oksanen et al., 2013). Before that, Levene’s and Shapiro-Wilk’s tests were used to confirm that assumptions for normality and heterogeneity are met. Differences at P ≤ 0.05 were deemed statistically significant. Following ANOVA analysis, Tukey’s honestly significant test (Tukey’s HSD) was performed in RStudio to test differences for alpha diversity indices among sampling times with the R package agricolae. Additionally, statistical related differences in the bacterial and eukaryotic communities’ composition were investigated using a pairwise PERMANOVA design. The ten most abundant phyla (bacteria) and class (eukaryote) across datasets related to the ‘Type’ and ‘Sampling’ factors were visualized using bar plots based on standardized, average abundance across replicates.

Comparison of relative alpha diversity indices for the bacterial dataset indicated significant shifts across sampling times (relative ASV richness and Shannon Index, both p ≤ 0.01, One way ANOVA) ( Table S2 ). The three-week samples yielded the highest relative bacterial diversity estimates ( Table S2 and Figure 2 ), significantly different from other sampling times for ASV richness (Tukey’s HSD test p ≤ 0.01) and from weeks 12 and 18 samples for Shannon diversity (Tukey’s HSD test p = 0.006 and p = 0.0008, respectively). The lowest overall bacterial ASV richness was observed for 18-weeks samples, significantly less than that during sampling weeks 1-6 (Tukey’s HSD p ≤ 0.0001).

Principal coordinate analysis (PCoA) displayed a clear temporal shift in the bacterial community structure ( Figure 3A ), with increasing dispersion towards the later sampling times (i.e., week-6 and later) and convergence of bacterial communities at the two latest sampling time points (weeks 12 and 18). PERMANOVA analysis confirmed the significant effect of time (i.e., ‘Sampling time’) (p ≤ 0.001). Besides, neither the effect of the ‘Side’ factor (p = 0.458) nor the marginal effect of the ‘Type’ factor (p = 0.0804) on bacterial community structure was detected ( Table S3 ). A test of homogeneity of dispersion (PERMDISP) on the bacterial dataset revealed significant differences in the average distances to centroids among samples grouped by sampling times (p= 0.0016, F = 5,4394, PERMDISP), supporting the observed time-related increasing dispersion in PCoA plots.

In the eukaryotic dataset, a comparison of alpha diversity indices revealed significant changes across sampling time points for both ASV richness and Shannon index (One way ANOVA p ≤ 0.001) ( Table S2 ). The highest diversity estimates were observed in the 6-week samples, differing significantly from 18-week samples based on ASV richness (Tukey’s HSD test p ≤ 0.005), and from week-1, week-3, and week-18 samples based on Shannon index (Tukey’s HSD test p ≤ 0.005). Minimum diversity was observed in the week-1 samples, with significant differences in ASV richness compared with week-18 samples (p ≤ 0.05) and significant differences among week-6, week-12, and week-18 samples based on Shannon index values (p ≤ 0.0001).

Similar to the bacterial dataset, PCoA revealed clear shifts in eukaryotic community structure for the ‘Sampling time’ factor ( Figure 3B ). Eukaryotic communities from the first two sampling weeks are closely related to each other and are distinct from other sampling time points, followed by a noticeable shift in community composition for week-6 samples. Eukaryotic communities at weeks 12 and 18 are most closely related and overlap in the multivariate space. PERMANOVA analysis further confirmed significant differences in samples distribution for the ‘Sampling’ factor (p ≤ 0.001), whereas no statistical differences were found for the ‘Side’ (p = 0.225) or ‘Type’ (p = 0.1191) factors ( Table S3 ). Statistical analysis of community dispersion (PERMDISP) also revealed significant differences in homogeneity of dispersion for the factor ‘Sampling’ (p= 0.001, F = 23,796, PERMDISP).

The overview of taxonomic succession on the substrates revealed higher temporal consistency in the bacterial community data across the tested bio-polymers compared with eukaryotes ( Figure 4 ). Proteobacteria and Bacteroidetes were the most abundant phyla in the bacterial dataset regardless of sampling time or concentration of oyster shells, comprising 66.7% and 20.1% of the overall community, respectively. The high-level taxonomic composition of the bacterial assemblage was similar across polymer types for the first two sampling time points (i.e., weeks 1 and 3), with some changes observed in week-6 samples, with a 3-6% increase in the relative abundance of Epsilonbacteraeota and Firmicutes in the OYS40 and OYS60 treatments. This observation is consistent with the type-related statistical differences reported in the pairwise PERMANOVA analysis, with significant statistical differences in community structure observed between the OYS40 and OYS60 week-6 samples ( Table S4 and Figure 4 ). At the two latest sampling times (i.e., weeks-12 and -18), bacterial community composition was similar regardless of polymer types. A minor change was observed in the control and OYS20 samples at 18 weeks, where the contribution of Actinobacteria was sensitively higher compared to other polymer types ( Figure 4 ).

The eukaryote dataset was dominated by Metazoa (84.1% of the overall read abundance). Among the ten most-represented classes ( Figure 4 ), the majority were common biofoulers: Cnidaria (40.5%), Urochordata (35.1%), and algae (e.g., Florideophyceae, Dinophyceae, 12.3%). On week-1 and week-3, Cnidaria dominated eukaryotic assemblages in all treatments, comprising 83%-91.6% of the overall read abundance, and with a noticeable increase in Arthropoda (8.5%) on week-3. The greatest variation in community structure (at the ASV level) across bio-polymer types was detected on week 3, as indicated by the pairwise PERMANOVA analyses ( Table S2 ). On week-6, a clear type-dependent shift in dominant taxa and community structure was observed. Thus, the control samples were dominated by Urochordata (40%), with Cnidaria, Arthropoda, and Plathelminthes comprising a substantial portion of the community (16.4%, 9.1%, and 18.2%, respectively). The OYS20 were dominated by Dinophyceae (27.7%), Urochordata (21.7%), unknown eukaryotes (16.2%), and Cnidaria (14.4%). In the OYS40 treatment, there was a notable decrease in Cnidaria’ relative abundance on week 6 (from 81.6% at week-3 to 23.1% at week-6), as well as an increasing contribution of Dinophyceae, Urochordata, Plathelminthes, and Syndiniales, 25.8%, 9.8%, 11.7%, and 11.2%, respectively. At this time, a substantial proportion of the eukaryotic community in the shell-enriched treatments (OYS40 and OYS60) was comprised of undertermine Syndiniales, a diverse parasitic group of Dinoflagellates (11.2% and 15.1%, respectively). On week–6, tunicates became dominant on OYS60 (59.2%), and later their relative contribution increased across all polymer types, peaking on week-18 in the control and OYS20 treatments (80.6% and 83.4%, respectively). In the latest weeks, the most pronounced difference between bio-polymer treatments was a comparatively higher contribution of red algae (Floridophyceae) in the OYS60 (from 3.9% to 13.4% on week-12 and -18, respectively) and a lower contribution of Nematoda in the control samples (from 3.1% to 2.1% on week-12 and -18 respectively). At the ASV level, no significant difference in eukaryotic community structure was detected between polymers at week-18 ( Table S4 ).

From the eukaryote dataset, 344 ASVs were assigned with high confidence (>99% identity) to 33 putative NIS (species or genus level), contributing up to 87.7% of the total sequence abundance at the first sampling time point. Overall, putative NIS dominated eukaryotic assemblages across all sampling time points ( Figure 5 ). Among these, tunicates (e.g., Ascidiella sp. Botryllus planus, Ciona intestinalis, Ciona savignyi, Clavelina meridionalis, Styela plicata) were the most prevalent, followed by red algae (e.g., Heterothamnion sp., Neosiphonia harveyi, Polysiphonia schneideri), green algae (Chaetomorpha crassa, Ulva sp. Cladophora glomerata) sponges (Hymneniacidon sinapium, Leucosolenia sp. Mycale sp.) and polychaete worms (Spirobranchus caraniferus, Pomatoceros lamarckii) (see Table S5 ). One-way ANOVA revealed a significant effect of the ‘Sampling’ factor (p = 0.0001) on NIS relative read abundances but no effect of the ‘Type’ factor (p = 0.242). However, there was an apparent difference in relative sequence read contribution of individual NIS between bio-polymer treatments. Thus, following week-1 and week-3 timepoints during which eukaryotic communities were heavily dominated by the invasive hydrozoan Bougainvillia sp., the total NIS contribution decreased substantially across all bio-polymer types. This change coincided with the overall community shift mentioned above and was characterized by a higher divergence in NIS contribution across polymer treatments ( Figure 5 ). Thus, on week-6, C. intestinalis and Bougainvillia sp. comprised a substantial part of the control and OYS20 communities (25.8% and 15.1% for C. intestinalis and 8.2% and 10.5% for Bougainvillia sp., respectively). The OYS40 samples were characterized by a high contribution of Bougainvillia sp. (20.3%) and C. intestinalis (7.1%). In the OYS60 treatment, invasives tunicates S. plicata and C. intestinalis contributed up to 52% of the overall diversity. On week-12, the control assemblage again became dominated by NIS, with B. planus, S. plicata, and C. intestinalis comprising 49.3%, 13.7%, and 9.8%, respectively ( Figure 5 ). Interestingly, the relative contribution of invasive tunicates increased over time across all bio-polymer treatments, particularly at the two latest sampling time points (i.e., 12 and 18 weeks), with B. planus being the most abundant, followed by S. plicata and C. intestinalis. The red algae N. harveyi was exclusively detected in week-18 OYS60 samples, contributing 12.1% of the entire diversity ( Figure 5 ).

In total, 41 bacterial families belonging to the phyla Proteobacteria, Firmicutes, Actinobacteria, and Bacteroidetes were identified with potential plastic-degrading capabilities in the bacterial dataset. Our analysis shows that within this dataset, Proteobacteria was the most frequently detected bacterial phylum of plastic degraders, comprising Gammaproteobacteria (54.1%), Alphaproteobacteria (16.7%), Betaproteobacteria (16.7%), Deltaproteobacteria (8.3%). Output from the PlasticDB tool revealed six families reported to degrade PBS: Streptomycetaceae (Actinobacteria); Thermoactinomycetaceae, Paenibacillaceae, and Bacillaceae (Firmicutes); Comamonadaceae and Pseudomonadaceae (Proteobacteria). In addition, the remaining 35 bacterial families detected in the dataset were reported to degrade other types of plastics ( Figure 6 and Figure S2 ). Among the taxa reported to degrade PBS, we observed that Streptomycetaceae was exclusively associated with shell-enriched bio-polymers. All five other families reported to degrade PBS were detected in all bio-polymer samples, both shell-enriched and controls. Among the families that have been identified by the PlasticDB tool as having species reported to degrade other plastic types: 32 families were detected in all shell-enriched and control samples; Six families were detected only in shell-enriched samples (OYS20, OYS40, and OYS60), and three families were detected in the control samples and some but not all shell-enriched samples (i.e. Sphingobacteriaceae and Halomonadaceae detected in control, OYS20, and OYS40 samples; Caulobacteraceae detected in control, OYS40 and OYS60 samples) ( Figure 6 ). PERMANOVA analysis revealed a significant effect of the factor ‘Sampling’ (p = 0.001) on the relative read abundance of putative bacterial plastic degraders, but no effect of the ‘Type’ factor (p = 0.085).

The community change of putative plastic-degrading bacterial families across polymer types revealed a similar pattern for the control, OYS20, and OYS40 samples ( Figure 7 ). On week-1, we observed a substantial contribution of putative degraders, followed by a progressive decline on weeks 3 and 6. Interestingly, the temporal evolution was different in the OYS60 treatment, with a notable increase in the contribution of putative bacterial plastic degraders from week-1 to week-3. However, a general decline was observed on week six across all polymer types, coinciding with the major community shift described above. At the family level, Illumatobacteraceae was the most persistent bacterial family overall, being detected in all polymer types and across all sampling time points. The second most abundant family was Flavobacteriaceae, particularly during the first two sampling time points, with an increased contribution on week-3 and a substantial drop later on (particularly in OYS60 treatment). The family Spongiibacteraceae was exclusively detected during the first two sampling time points across polymer types, contributing up to 30% of the overall diversity in the control, OYS20, and OYS40 treatments on week 1, decreasing afterward. In the OYS60 samples, the relative contribution of Spongiibacteraceae increased from week-1 to week-3. The relative contribution of the Kordiimonadaceae increased towards the later sampling time in the control while being underrepresented in the shell-enriched treatments.

Consistent with several earlier studies (Oberbeckmann et al., 2014; Zaiko et al., 2016; Oberbeckmann and Labrenz, 2020), our results revealed strong temporal patterns associated with both bacterial and eukaryotic community structure on all tested bio-polymers. The period of sample immersion in seawater was more important in shaping biofouling assemblages than the concentration of oyster shell filler, as evidenced by finding similar biological communities between oyster shell-enriched substrates and controls. Most studies that have reported plastic-specific recruitment (Zettler et al., 2013; Oberbeckmann et al., 2018; Hou et al., 2021; Martínez-Campos et al., 2021) have focused on a variety of chemically-distinct polymers (e.g., polystyrene, polyethylene, and polypropylene). Our experiment used the same base polymer (PBS) to manufacture control and treatment substrates. This might explain the lack of high-level taxonomic differentiation between bio-polymer treatments (Lakshmi et al., 2012; Muthukrishnan et al., 2019). However, differences in community structure between substrates were detected at a finer temporal scale, with significant community shifts identified on week-3 for eukaryotes and week-6 for bacteria ( Table S4 ).

Consistent with other studies (Martell et al., 2018; Enrichetti et al., 2021), hydrozoans (cnidarians) acted as primary colonizers, dominating eukaryotic assemblages on all substrate types at the first two sampling time points. In week-3, differences between bio-polymers were mainly driven by the differential abundance of arthropods, ascidians, and parasitic Syndiniales ( Figure 4 ). The community shift in the eukaryotic assemblage could be explained by the gradual increase of biomass on the bio-polymers, which in turn favored the recruitment of macrofoulers (Qian et al., 2007; Dussud et al., 2018a; Dussud et al., 2018b; Abed et al., 2019). In addition, we observed a slight increase in parasitic dinoflagellates (Syndiniales) relative abundance at week-3. Parasites in the order Syndiniales are widespread in marine environments and can infect several types of marine eukaryotes, including crustaceans, dinoflagellates, and cnidarians (Lee et al., 2008; Elifantz et al., 2013; De Tender et al., 2015). Altogether, our results suggest that the significant changes in eukaryotic community composition among substrates at week-3 were dictated by the progressive stabilization of the biofouling assemblage, together with the recruitment of fast-growing macro-invertebrates.

Proteobacteria and Bacteroidetes were the most dominant phyla across sampling time points in the bacterial dataset. Consistently with our findings, previous studies have shown that Proteobacteria act as pioneers in the early stage of biofilm formation in the marine environment, while Bacteroidetes abundance increase over time, acting as a secondary colonizer (Lee et al., 2008; Elifantz et al., 2013; De Tender et al., 2015). While previous reports have highlighted the dominance of Alphaproteobacteria colonizing plastic particles (mainly polyethylene, polypropylene and polystyrene) in the Atlantic Ocean (Debroas et al., 2017), and the Mediterranean Sea (Dussud et al., 2018a; Dussud et al., 2018b), our results revealed that Gammaproteobacteria was the dominant class in the bacterial dataset. Specific polymer-type effect or different study location could potentially explain these differences (Oberbeckmann et al., 2016). Overall, the bio-polymer type-related differences in the bacterial community composition were much more subtle than eukaryotes. The biggest shift observed on week-6 was mainly associated with was associated with increased Epsilonbacteraeota and Firmicutes in shell-enriched treatments. Interestingly, Clostridium and Arcobacter genera, which respectively belong to Firmicutes and Epsilonbacteraeota, have been described in the gut microbiome of the invasive tunicate C. intestinalis, potentially contributing to its invasiveness (Utermann, 2021). Moreover, the Arcobacter genus, which contain species with the potential for plastic degradation and pathogenesis has recently been reported in high abundance on microplastic particles in effluent sewage and estuarine environments (Kelly et al., 2021). Therefore, we assume that the increased contribution of these two bacterial taxa is either linked to high abundance of C. intestinalis detected on week-6, or could be associated with unidentified environmental sources, given the fact that the experiment was conducted in area exposed to environmental pollution.

Other finer-scale differences observed in bacterial communities between bio-polymers could be explained by the likely dominance of bacterial taxa associated with specific eukaryotic biofoulers. Marine invertebrates are holobionts: a functional ecological unit associating the host and its various microbial communities, providing host-beneficial functions (Hoffmann et al., 2009; Thompson et al., 2015). For example, numerous studies identified Proteobacteria and Bacteroidetes as dominant phyla associated with the ascidian microbiome (Erwin et al., 2013; Dror et al., 2019), cnidarians (Brown et al., 2017), and crustaceans (Datta et al., 2018). Hence, the recruitment of particular eukaryotic taxa and change in their relative abundance through time might be directly linked with the change in relative abundance observed in bacteria (Utermann et al., 2020), disguising, to some extent, bio-polymer-related changes in bacterial community composition.

Altogether, our results revealed a highly dynamic assemblage of both bacteria and eukaryotes at the early sampling time points (1-6 weeks). Two months after deployment (i.e., from sampling week-12 onwards), bacterial and eukaryotic communities became more stable, likely indicating the establishment of mature biofouling, significantly different from the early colonization (biofilm) stage. In an earlier study, Pollet et al. (2018) reported an increase in the stability of microbial community composition associated with plastic deployed in a coastal environment after 28 days of immersion. The authors suggested that plastic-associated biofilms might reach a ‘steady state’ over time and hypothesized that the temporal stability of bacterial communities in later succession stages could be explained with the capture of bacterial biofilms developed at the surface of macro-organisms colonizing the substrates (Pollet et al., 2018). Our results support these conclusions and indicate that in the context of the New Zealand coastal environment, the bacterial community associated with tested bio-polymers became more stable between 6 and 12 weeks after immersion.

Putative NIS dominated the eukaryotic assemblages across all sampling time points and bio-polymer types, with tunicates, algae, sponges, and polychaete worms being the most prevalent NIS taxa (see Table S5 ). These species are usual ‘inhabitants’ of the artificial substrates around port and marinas and were previously reported from the study area (Fletcher et al., 2017; Hawes et al., 2018). Ports and marinas are well-known gates for vessel-mediated biological invasions (Ardura et al., 2015; Borrell et al., 2017; Miralles et al., 2018; Floerl et al., 2021) and hotspots of the local NIS populations (Seebens et al., 2013; Lyons et al., 2020).

Previous studies have shown that the location and timing of introducing a new substrate into the marine environment can considerably affect the recruitment and succession pattern of biofouling assemblages (Jenkins and Martins, 2010; Zaiko et al., 2016). In the present study, bio-polymers were deployed during the summer-fall months (i.e., between January and May), corresponding to the spawning season of many non-indigenous marine invertebrates in New Zealand (O’Driscoll et al., 2003). Hence, the increase of propagules in the water column during this time of the year, coupled with the high abundance of NIS in the study area, might explain the dominance of non-indigenous metazoan detected in the eukaryotic dataset. Moreover, marinas are typical hotspots for the establishment and subsequent spread of non-indigenous species due to the proximity of introduction pathways (recreational watercraft, affinity to shipping ports) and availability of artificial habitats (floating pontoons, ropes, and buoys) suitable for the incursion of biofouling taxa (Afonso et al., 2020). NIS usually exhibit more “opportunistic traits” than native species (e.g., short lifespan, fast growth rates, high fitness); hence, NIS are often more successful than native species in pioneering and establishing on artificial substrates (Ferrario et al., 2017; Quell et al., 2021).

While previous studies reported diatoms (Mincer et al., 2016), Arthropoda (Delacuvellerie et al., 2021), fungi (Kettner et al., 2017), and ascidians (Vieira et al., 2018) as pioneer colonizers of artificial surfaces in the marine environment, our results show that Cnidaria, and particularly the invasive hydroid of the Bougainvillia genus, was the main taxon detected in the eukaryotic dataset during the early stages of recruitment (week-1 and -3). Bougainvillia’s explosive growth could have been favored by the chemical cues from the bio-polymers (Wyatt et al., 2005; Grey et al., 2018) and/or the lack of competition during the early stage of colonization (Rinkevich and Sakai, 2001; Lindquist, 2002). Although members of the family Bougainvilliidae are found worldwide (Vannucci and Rees, 1961), the genus has been described as colonizing artificial surfaces in coastal waters (mostly ports and marinas) from the South Pacific Ocean, particularly in Australia (Wyatt et al., 2005; Grey et al., 2018), and New-Zealand (Schuchert, 1998; Cairns et al., 2009). A recent metabarcoding study of biofouling from submerged surfaces (buoys, ropes, tires, and pontoon structures) in Tahiti revealed the highest relative abundance of Bougainvillia among detected NIS (Pearman et al., 2021). Comparison of the seasonal recruitment of fouling hydroids on settlement panels between cold and warm months in Brazil, Migotto et al. (2001) revealed that Bougainvillia spp. tend to dominate during the warmer months. Similar to our observations, they reported a low persistence of established colonies reflected by a sharp decline during the second month of immersion.

During the observed community shift on week-6, total NIS contribution decreased substantially across all polymer types, particularly in OYS20 and OYS40 treatments, where other eukaryotes became dominant. Interestingly, species richness in OYS20 and OYS40 treatments was slightly higher than other polymer types. Higher NIS contribution in the control and OYS60 samples probably results from the presence of key taxa, potentially depressing species richness once established on the substrate. For example, Blum et al. (2007) compared the species richness between panels from which C. intestinalis was experimentally removed and panels with naturally dense populations in San Francisco Bay, revealing that species richness was negatively correlated with the tunicates abundance. The authors suggested that many species were occurring more often in the absence of C. intestinalis (e.g., Styela spp. and Didemnum spp.). In contrast, other species (e.g., Botrylloides spp.) occurred more often on C. intestinalis-dominated assemblages. Our results support these observations and indicate that a higher contribution of C. intestinalis was found in Control and OYS60 treatments, compared with OYS20 and OYS40.

In general, the relative abundance of invasive tunicates increased over time across all substrates, particularly on weeks 12 and 18. Botryllus planus was the most abundant among tunicates, followed by Styela plicata and C. intestinalis. Interestingly, the contribution of B. planus was somewhat related to the oyster shell concentration in bio-polymers ( Figure 5 ). As previously observed by Dalby Jr and Young (1993), B. planus often occurs on oyster farms, potentially impacting their growth rate. In addition, the clubbed tunicate S. plicata is often reported forming dense clumps on any submerged surfaces, including shells of the commercially important oyster Crassostrea gigas (Rodriguez and Ibarra-Obando, 2008; Marins et al., 2010). In our experiment, this tunicate seemed to recruit earlier on OYS60 bio-polymers, suggesting a specific affinity of this species to the shell-enriched substrate. Similarly, the invasive tunicate C. intestinalis has been described as a significant fouling species in shellfish farms, being able to colonize nearly 100% of the available substrate (Howes et al., 2007; LeBlanc et al., 2007). Today, C. intestinalis represents a major problem for mussel farmers on the Atlantic coast of Nova Scotia (Clancey and Hinton, 2003). Nevertheless, there was no evidence of C. intestinalis preferring the shell-enriched substrates, possibly due to interspecific competition with other colonial ascidians. For example, Paetzold et al. (2012) reported lower recruitment of C. intestinalis on aquaculture gear when other colonial Botryllus spp. tunicates were present. The authors hypothesized that pre-settled substrates might impede the recruitment of C. intestinalis through competition for space or release of negative settlement cues produced by other colonial ascidians. In this study, B. planus dominated the eukaryotic assemblages at 12 and 18 weeks, particularly on shell-enriched bio-polymers (i.e., OYS40 and OYS60), likely repulsing the attachment of C. intestinalis through chemical deterrence of the tunicate larvae (Krug, 2006).

Bacterial bio-degraders play a key role in the natural decontamination and bioremediation of environmental contaminants (Devi et al., 2016). Organisms proposed capable of plastic biodegradation are phylogenetically spread around the bacterial tree of life, with over 286 bacteria from five different phyla reported in the literature (Gambarini et al., 2021). Polybutylene succinate is a biodegradable polymer with excellent biodegradability, thermoplastic processability, and balanced mechanical properties (Xu and Guo, 2010). Many studies on PBS have reported significant biodegradation within several months in different environmental matrices such as soil (compost), water, and seawater (Nishioka et al., 1994; Hirotsu et al., 2000; Morohoshi et al., 2018). Moreover, PBS can be naturally degraded in the environment by bacteria (Teeraphatpornchai et al., 2003) and fungi (Kim and Rhee, 2003). The Comamonadaceae family was one of the six families identified with PBS degraders. A search on the PlasticDB webserver revealed only one enzyme described in the literature with depolymerase activity towards PBS. This enzyme was isolated from the bacterium Acidovorax delafieldii, which belongs to the family Comamonadaceae, indicating that bacteria belonging to this family may be associated with biodegradation of PBS in our study.

Overall, our results show that the different concentrations of oyster-shell did not affect the specific recruitment of presumed bacterial degraders, although significant changes were observed between sampling time points ( Figure 7 ). Early time points had a greater abundance of bacterial families with plastic biodegradation reports, while later time points displayed a decrease in abundance of these families. Erni-Cassola et al. (2020) reported this phenomenon for polyethylene strips incubated in situ in coastal seawater. The authors investigated the abundance of bacteria capable of biodegrading hydrocarbons. They found that taxa that might break down hydrocarbons are abundant in very early colonizing populations (two days), following a decrease in their relative abundance at a later time point (nine days). The authors hypothesized that this phenomenon was driven by the depletion of labile, plastic-derived substrates formed by weathering. After this initial depletion, non-hydrolyzable polymers substrates may not serve as an essential carbon source for the microbial biofilm attached to the plastic. Interestingly, a similar pattern was observed on PBS even though, unlike polyethylene, PBS has been described as a biodegradable polymer (Rafiqah et al., 2021). Bacterial degraders detected in the present study are probably outnumbered by later, more generalist bacterial colonizers. Hence, the development of later biofilm stages might hinder our ability to detect putative PBS degraders.