Two 17 × 26 mm flat pieces of each of the three materials were used, polyethylene (PE) from plastic food bags (Ziploc), polyethylene terephthalate (PET) from transparent drinking water bottles (Coca-Cola, 1.5 L) and glass microscope slides (Marienfeld). The surfaces were tied with a fishing line to a straw and secured with plastic beads to create the “mobile” structures. Five mobiles were then secured along the dock in Herzliya Marina, Israel (32° 09′ 38.8″ N 34° 47′ 35.0″ E) such that the surfaces were submerged ∼0.5 m below the water surface (Figure 1). After 1 month each material from the mobile was sampled for DNA extraction and metabarcoding or microscopy.

Each of the mobile materials were sampled in its entirety, gently washed three times for 5 min with filtered seawater to remove unbound material and separately processed in each of the subsequent assays and procedures described below. Seawater was sampled in proximity to each of the mobiles (five times) using sterile sampling bottles. 0.5 L of the sampled water was filtered on 0.22 μm polyethersulfone membrane (Millipore) using a 20 L/min pump (MRC). DNA was extracted using the phenol–chloroform extraction method (Debeljak et al., 2017; Supplementary Material, Protocol 1).

One sample of each material (water filters, PE, glass and PET) from each mobile (five repeats) were subjected to PCR amplification with three sets of primers to amplify three barcode regions. The complete 16S rRNA gene was amplified using 27F and 1492R primers (Weisburg et al., 1991), with an expected product size of ∼1.5 kb. Primers 566F and 1289R (Hadziavdic et al., 2014) were used to amplify the V4 and V5 regions of the 18S rRNA gene, with an expected product size of ∼0.7 kb, and ITS86F and ITS4R (Op De Beeck et al., 2014) were used to amplify the fungal internal transcribed spacer 2 (ITS2) loci, with an expected product size of ∼0.4 kb. The amplification parameters and primer details are listed in Supplementary Material 1, Table 1. The reaction volume was 50 μL with 25–75 ng of template sample DNA. The PCR products were cleaned with QIAquick PCR Purification kit (QIAGEN) to meet the criteria of the MinION nanopore library preparation protocol (Karamitros and Magiorkinis, 2018).

The sequencing libraries were prepared using the Native barcoding amplicons protocol with EXP-NBD104 and SQK-LSK109 kits (Oxford Nanopore Technologies) according to Davidov et al. (2020). The 16S rRNA and 18S rRNA gene sequencing libraries were loaded to the same MinION flow cell in two batches: the first included the amplified 16S rRNA and 18S rRNA gene DNA barcodes of replicates 1, 3, and 5 (12 multiplexed libraries in total), and the second included the amplified 16S rRNA and 18S rRNA gene DNA barcodes of replicates 2 and 4 (8 multiplexed libraries in total). The two sequencing runs were separated by a washing step using EXP-WSH002 Kit (Oxford Nanopore Technologies). Each library was loaded onto to the MinION Nanopore Spot-on flow cell (FLO-MIN106D, version R9) and sequenced until reaching ∼7 Giga nucleotides (∼4 M reads). The ITS2 sequencing library was loaded to a new MinION flow cell and sequenced until reaching ∼1.3 Giga nucleotides (∼1 M reads). Base-calling for all libraries were done by the Guppy base calling software 3.3.3, using MinKNOW program with the “high accuracy” option. Raw reads were obtained in FAST5 and FASTq formats from which “pass” quality reads were subjected to further analysis.

Processing and analysis of reads was performed using MetONTIIME pipeline and QIIME2 plugins (Bolyen et al., 2019). The MetONTIIME pipeline was executed using the script “Launch_MinION_mobile_lab.sh” and “MetONTIIME.sh”. Configuration of MetONTIIME pipeline was performed using “config_MinION_mobile_lab” R script. Launch_MinION_mobile_lab.sh script was used for performing basecalling, demultiplexing, quality filtering, adapters and PCR primers trimming. This resulted in reads, filtered by length, quality and separated by barcodes into files in fastq format. Sequences were filtered based on read quality (min_quality; 16S and 18S, 9–11; ITS2, 8), and read length was restricted (amplicon_length X, lenfil_tol Y) based on read length histograms, to give the following range: 16S, 1399–1599 nucleotides; 18S, 605–755 nucleotides; ITS2, 100–1100 nucleotides. The following steps have been performed separately for each barcode. The MetONTIIME.sh script was run for performing dereplication, clustering and taxonomic classification. Sequences were clustered into consensus sequences with default MetONTIIME pipeline parameters (de novo strategy, clustering threshold parameter perc-identity 1). Consensus sequences of the 16S rRNA and 18S rRNA genes were assigned based on SILVA 132 database (Quast et al., 2013). The following parameters were used: taxonomic classifier Blast, max_accepts 1, query coverage 0.8 and identity threshold 0.85. Consensus sequences of the ITS barcode were classified against UNITE database V8 (Nilsson et al., 2019) using taxonomic classifier Blast, max_accepts 1, query coverage 0.8 and identity threshold 0.7. We chose max_accepts 1 for faster analysis, after testing 3 and 10 values for this command with raw data subsets which came up with the same results. Then, the feature tables of the used barcodes generated using the MetONTIIME pipeline were combined using the q2-feature-table plugin of the QIIME2, using the methods: merge, merge-seqs and merge-taxa. The code used is included in Supplementary Material 2. The featured tables were imported into R using phyloseq (R Core Team, 2018).

Sequences mapped to chloroplasts, mitochondria, eukaryotes, and “unknown organisms” were removed from the 16S rRNA gene analysis, while sequences mapped to bacteria (due to mitochondrial DNA) and “unknown organisms” were removed from the 18S rRNA gene and ITS analysis. Alpha diversity was estimated by richness [observed Operational Taxonomic Units (OTUs)], together with Shannon and Pielou diversity indexes using phyloseq (McMurdie and Holmes, 2013). After alpha diversity analysis, OTUs with a read count of 1 were excluded and the remaining OTUs were normalized to relative read abundance (dividing the number of reads for each sample by the total reads count). All further analyses were based on relative read abundance. For Beta diversity we used phyloseq (McMurdie and Holmes, 2013) to perform principle coordinates analysis (PCoA). This ordinated the sequences using the Bray–Curtis distance matrix to visualize multivariate structures of the communities. To test significant differences between the groups Permutational multivariate analysis of variance (PERMANOVA) tests were calculated using MicrobiomeAnalyst (Chong et al., 2020). Stacked bar plots and heatmaps were produced using MicrobiomeAnalyst (Chong et al., 2020) and Ampvis2 (Andersen et al., 2018), respectively. We used limma (Ritchie et al., 2015) to generate linear models and the voom function to transform normalized counts to log2-counts-per-million (logCPM) and estimate mean-variance of species between surface materials, similar to Calgaro et al. (2020). As the five “mobile” replicates were placed in slightly different locations within the marina, the mean location variance was removed from all measurements so the community differences due to surface type could be better seen. OTU differential abundance was visualized using Glimma (Su et al., 2017) and OTUs that had false discovery rate (FDR) adjusted p-value < 0.05 and log-fold-change of >0 were identified as differentially more abundant. Venn diagrams were created with InteractiVenn tool (Heberle et al., 2015), based on OTUs represented by at least two reads. All Nanopore MinION filtered reads analyzed in this project were deposited in the NCBI SRA database (Agarwala et al., 2016) under BioProject PRJNA680232 (accession numbers: SRX9553811–SRX9553850, SRX10393298–SRX10393305, and SRX11261845–SRX11261863).

Scanning electron microscopy (SEM) was used to visualize the broad composition and the morphology of the biofilms on both PE and PET samples (Supplementary Material 1, Protocol 2). For fungal identification, Lactophenol Cotton Blue dye (Sigma-Aldrich) was applied to PE samples and visualized under a light microscope (Supplementary Material 1, Protocol 3). A modified crystal violet method (Supplementary Material 1, Protocol 4) was used for the indirect measurement of the relative biofilm biomass on the experimental samples.

Within 30 days from submersion in the marina waters, all substrates were heavily covered in biofilm (Figure 1C). The biomass was similar on all surfaces, as indicated by the crystal violet assay (Supplementary Material 1, Figure 1). SEM imaging revealed that the plastic surfaces were mostly covered with a likely bacterial biofilm (Figure 2A), probable exposed bacteria (Figure 2B) and multiple diatom species (Figures 2C–E). We were also able to identify single cell eukaryotes (Figure 2E), bryozoans (Figure 2F) and fungal-like fruiting bodies (Figure 2G). Furthermore, we observed multicellular filaments (Figure 2H), of which some were identified as algae based on the presence of chlorophyll and some were confirmed to be fungi hyphae by Lactophenol Cotton Blue staining (Figure 2I).

The nanopore sequencing run produced an average of 25,235 16S rRNA gene reads, 80,125 18S rRNA gene reads, and 24,777 ITS reads per sample. The average read length was 1,428 nucleotides for 16S rRNA gene, 669 for 18S rRNA gene and 369 for the ITS barcode. However, while both 16S rRNA and 18S rRNA gene amplicons had read length standard deviation of 24 and 20 accordingly, the ITS reads varied greatly in length with standard deviation of 154.5 nucleotides. The mapping rates of the ITS sequences to the reference databases were very low (3% on average) compared to the 16S rRNA and 18S rRNA gene mapping rates (both 99%), and the ITS sequences had a lower taxonomic resolution (Supplementary Material 1, Table 2).

Community complexity parameters including richness (number of observed OTUs), evenness (Pielou’s index) and diversity (Shannon’s index) were obtained for each of the sample types (water, PET, PE, and glass) based on 16S rRNA and 18S rRNA gene reads (Supplementary Material 1, Figure 2). For 16S rRNA gene sequence analyses, the water had significantly lower taxonomic evenness, contributing to lower OTU diversity (Shannon’s index) in comparison to the surface samples (Supplementary Material 1, Figure 2A). This agrees with previous environmental studies that showed the same trend (De Tender et al., 2015; Bryant et al., 2016; Didier et al., 2017). In contrary, the 18S rRNA gene analysis resulted in no significant difference among the samples in any of the diversity indexes (Supplementary Material 1, Figure 2B). Previous 18S rRNA gene based analyses report water communities as more diverse (Didier et al., 2017; Kettner et al., 2019; Dudek et al., 2020). Because of the low ITS mapping rates, we did not analyze the alpha diversity parameters for fungi.

To assess the similarities in the community composition between the samples, beta diversity analysis was performed using a PCoA analysis (Figure 3). In both 16S rRNA and 18S rRNA gene analyses, the five water samples formed distinctive clusters from the surface samples (Figures 3A–D). PERMANOVA analyses confirmed the water samples were significantly different from the other surfaces. In the ITS analysis, while the water samples showed homogeneity among themselves, they did not form a separate cluster from the other samples (Figures 3E,F). Despite this, PERMANOVA analyses confirmed the ITS water samples were significantly different from the other surfaces. When excluding the water samples, no distinct clusters were observed in any of the analyses (Figures 3B,D,F).

The filtered, high quality 16S rRNA gene sequences were clustered into an average of 2,733 operational taxonomic units (OTUs) for PET, 21,129 for PE and 4,163 for glass, corresponding to 524, 754, and 677 classified organisms, respectively. The top 10 most abundant prokaryotic genera within each treatment, that were identified based on the relative abundance of mapped read counts, contain some different genera across the surfaces (Figure 4A). Ekhidna, Muricauda, and Portibacter were included in the top 10 most abundant genera on PE and not in the other surfaces, whereas Rhodopirellula and OM60 clade were included only in the PET top 10. Among the above genera, Muricauda, and Rhodopirellula have been previously reported as hydrocarbon-degrading bacteria (Jiménez et al., 2011; Didier et al., 2017; Laso-Pérez et al., 2019; de Araujo et al., 2021).

To identify species with preference to plastic, we searched for species with significantly higher read ratios in the PET and PE samples as opposed to glass. Linear discriminant analysis of our 16S rRNA gene sequences revealed 35 OTUs, corresponding to 30 bacteria database matches, that were significantly more abundant in the PE samples than in the glass samples (Figure 4B and Supplementary Material 1, Figure 3). Of the 30 bacteria matches, 9 were resolved at the species level and 10 at the genera level (leaving 11 of higher taxonomic classification).

Among the bacterial species that showed significantly higher read representation in the PE samples compared with the glass and water samples, Muricauda sp. NBRC 112477 showed highest significance values (adjusted p-value = 0.003) (Supplementary Material 1, Table 3). We identified 3 other species of the genus Muricauda, which matches with the top 10 genera for PE (Figure 4A), and 2 species from the genus Winogradskyella, which has been previously reported as a hydrocarbon degrading genera (Wang et al., 2014). Other species that were significantly differentially abundant on PE included Maritimimonas rapanae, Robiginitalea myxolifaciens and Psychroserpens sp. MEBiC05023, all from the Flavobacteriaceae family, and Marivirga tractuosa. Of the unresolved sequences, Alcanivorax sp. is a well-known degrader of alkanes and petroleum (Yakimov et al., 2019). It was also recently shown that certain Alcanivorax species are able to grow and form biofilms when PE is the main carbon source (Delacuvellerie et al., 2019).

The high quality 18S rRNA gene sequences were clustered into an average of 10,997 OTUs for PET, 11,608 for PE, and 7,956 for glass, corresponding to 730, 735, and 506 classified eukaryotes, respectively (Supplementary Table 2). Within the 18S rRNA gene mapped reads, the genera with the highest relative abundance across all surfaces was the bryozoan Amathia. However, there were also differences between the top 10 most abundant genera among the surfaces (Figure 5A). On PE, the most abundant genera included Scytosiphon, a genus of brown seaweed, the bryozoan Bugulina, the entoprocta Barentsia and the sabellid polychaete Parasabella. On the other hand, the top 10 genera on PET included the copepods genus Acartia, as well as the diatom genus Nitzschia and the ciliate genus Dysteria.

Linear discriminant analysis of our 18S rRNA gene sequences identified five OTUs, corresponding to four species and one genus, with significantly higher read ratios in the plastic samples compared to glass samples (Supplementary Material 1, Table 4). These species included two brown algae; Petalonia fascia and Scytosiphon lomentaria, with higher read ratios in both PET and PE samples. The benthic dinoflagellate Prorocentrum emarginatum and the red alga Hydrolithon sp. had higher relative representation in the PE samples, while reads mapped to the diatom Haslea spicula were found in relative higher ratios in the PET samples (Figure 5B).

For the identification of fungi we used the ITS2 barcode. The ITS2 locus has been shown to be the most suitable taxonomic barcode for the characterization of fungal communities and has been previously used for the identification of fungi from marine plastic debris (Op De Beeck et al., 2014; Davidov et al., 2020). Due to the limitations of the available databases for marine fungi, only 3% of the reads were mapped to reference sequences. Despite this, our analysis still identified an average of 57 fungal OTUs in the PET samples and 58 in the PE samples, of which 23 and 25, respectively, were classified to the species level. We also identified an average of 87 fungal OTUs in the glass samples and 114 in the water samples, of which 29 and 41, respectively, were classified to the species level (Supplementary Material 1, Table 2).

Similar to the 16S rRNA and 18S rRNA gene based beta-diversity analyses, the fungi (barcoded with ITS) analysis did show distinction between the taxonomic composition of the water vs. the composition of the biota associated with surfaces (Figure 3E). Furthermore, within the ITS mapped reads, there were differences between the top 10 relatively abundant genera across the surfaces (Figure 6A). On PE, the most abundant genera included Neocatenulostroma, Penicillium, and Vishniacozyma. Whereas in the top 10 genera of the PET samples we identified Candida, Cyberlindnera, and Rhodosporidiobolus. Due to the limited mapping results, differential abundance analysis was not performed. Nevertheless, of all 316 OTUs identified, five were found only on PE (Figure 6B) corresponding to species Aspergillus penicillioides, Bipolaris sorokiniana, Filobasidium magnum, Knufia mediterranea and Ramichloridium cucurbitae (Figure 6C). Another five OTUs were only found on the PET samples corresponding to Cryptococcus aspenensis, Cyberlindnera jadinii, Debaryomyces vindobonensis, Pyrenochaetopsis leptospora, and Symmetrospora coprosmae. Another two OTUs corresponding to the species Candida sake, and Peniophora lycii were identified in the PET and the PE samples but not in the glass samples and the water.