The travelling particles: community dynamics of biofilms on microplastics transferred along a salinity gradient

Microplastics (MP), as novel substrata for microbial colonization within aquatic ecosystems, are a matter of growing concern due to their potential to propagate foreign or invasive species across different environments. MP are known to harbour a diversity of microorganisms, yet little is understood of the dynamics of their biofilms and their capacity to successfully displace these microorganisms across different aquatic ecosystems typically marked by steep salinity gradients. To address this, we performed an in situ sequential incubation experiment to simulate MP transport from riverine to coastal seawaters using synthetic (high-density polyethylene, HDPE and tyre wear, TW) and natural (Wood) substrata. Bacterial communities on incubated particles were compared to each other as well as to those in surrounding waters, and their dynamics along the gradient investigated. All communities differed significantly from each other in their overall structure along the salinity gradient and were shaped by different ecological processes. While HDPE communities were governed by environmental selection, those on TW and Wood were dominated by stochastic events of dispersal and drift. Upon transfer into coastal seawaters, an almost complete turnover was observed among HDPE and TW communities. While synthetic particles displaced a minor proportion of communities across the salinity gradient, some of these comprised putatively pathogenic and resistant taxa. Our findings present an extensive assessment of MP biofilms and their dynamics upon displacement across different aquatic systems, presenting new insights into the role of MP as transport vectors.


Detailed description of sample collection
Following each incubation, particles and surrounding waters were sampled in triplicate. Particles were stored in 1 L Kautex bottles filled with water from each respective site while water samples were stored in 1L Nalgene bottles. All sampling events were conducted at outgoing tide and samples collected were transported back to the laboratory in a cooling box (approx. 4 -5 °C) on the same day for immediate processing. After each sampling event, the particles were transported to the next site under in situ conditions, where the cage containing the particles was placed into a transport container filled with water from the respective incubation sites and kept at a temperature of approximately 4 -5 °C with ice packs.

Detailed description of sample preparation and amplicon sequencing
All particles within each capsule were transferred into 40 ml of 1 × phosphate buffered saline (PBS) and vortexed three times to remove any loosely-attached material, following which total biofilm DNA was extracted. The concentration of extracted DNA was quantified using the Quant-iT™ PicoGreen™ dsDNA Assay Kit (Invitrogen, Waltham, MA, USA) measured with a Tecan Infinite M200 PRO NanoQuant microplate reader (Tecan Trading AG, Switzerland).
PCR was performed in 20 µL reaction volumes containing 1 × MyTaq buffer (Bioline GmbH, Luckenwalde, Germany), 1.5 u MyTaq DNA polymerase (Bioline GmbH, Luckenwalde, Germany), 15 pmol of each forward and reverse primer (with a 10-nt barcode sequence at the 5'-end of each primer), 2 µl of BioStabII PCR Enhancer (Sigma-Aldrich Co.), and 1 µL template DNA. Thermocycling conditions were as follows: 60 s initial denaturation at 96 °C, followed by 30 cycles of denaturation for 15 s at 96 °C, annealing for 30 s at 55 °C, and elongation for 90 s at 70 °C. Amplicons were checked via gel electrophoresis, pooled, and purified using Agencourt AMPure XP beads (Beckman Coulter, Inc., IN, USA). An additional purification step was performed using MiniElute columns (QIAGEN GmbH, Hilden, Germany). Illumina libraries were constructed from the purified amplicon pools using the Ovation Rapid DR Multiplex System 1-96 (NuGEN Technologies, Inc., CA, USA) and subject to 300 bp paired-end read sequencing using an Illumina MiSeq (V3 chemistry).

Detailed description of 16S amplicon sequence processing
Forward and reverse reads were truncated at 205 bp and 245 bp, respectively, using the DADA2 plugin on the QIIME2 platform to remove low quality bases (median score < 30). Reads were dereplicated at 100 % sequence identity and denoised using a parametric error model. The resulting reads were then merged and a global alignment performed for de novo chimera detection and removal.

Details on differential abundance testing
Statistical tests were performed using the calc_diff_abund_deseq2 function offered by the metacoder package, which employs a negative binomial model and Wald test for hypothesis testing to identify differentially abundant taxa between sample types. The test was performed on raw read counts and hypothesis testing conducted at a significance level of α = 0.05, with a Benjamini-Hochberg correction applied to obtained p-values.

Details on permutational analysis of variance (PERMANOVA)
Variation partitioning was performed at a significance level of α = 0.05 using type III sum of squares and p-values obtained using 9999 permutations except for samples where too few possible permutations were available, during which Monte Carlo simulations were employed. Permuted p-values were adjusted using Benjamini-Hochberg correction.

Details on quantitative process estimates using null models
The compositional turnover of the different communities between sites was quantified through pairwise comparisons based on the ß-nearest taxon index (ßNTI). Values obtained were then compared against a null distribution, representative of stochastically-assembled communities, and differences measured in standard deviation units. Positive and negative deviations from the null distribution (ßNTI > 2, ßNTI < -2) are indicative of variable and homogenous selection, respectively. The former is defined as a high compositional turnover elicited by a shift in selective pressures whereas the latter suggests a low turnover due to exposure to relatively unvarying environmental conditions. Environmental selection is ruled out as the dominant assembly process when no significant deviations from the randomly-assembled communities are observed (-2 < ßNTI < 2), for which dispersal or ecological drift is inferred to be the primary influence. To estimate the relative contributions of these stochastic processes, the Raup-Crick metric [43] was calculated based on Bray-Curtis dissimilarities (RCbray) as a non-phylogenetic analogue to the ßNTI for the comparison of observed and expected taxonomic turnover. Similarly, deviations from the null expectation were suggestive of dispersal limitation (RCbray > 0.95) or homogenizing dispersal (RCbray > -0.95). Dispersal limitation describes a scenario where communities drift apart due to a low dispersal rate, resulting in a high compositional turnover, while homogenizing dispersal increases compositional similarities as a consequence of high dispersal rates overwhelming the effects of environmental selection. The final scenario, represented by values that fall within the null distribution (-0.95 < RCbray < 0.95), indicates a turnover undominated by selection or dispersal, otherwise referred to as ecological drift.

Fig. S1
Species richness curve generated from normalized reads based on Shannon's diversity index displaying a sufficient sampling depth. Reads were normalized to the lowest read count detected, resulting in 4036 reads per sample (indicated by the black line).

Fig. S2
Heat tree showing significant differences in the median read counts of taxa within total communities detected of different sample types. The trees are collapsed at the order level and taxa labels provided in the grey tree on the lower left. Green highlighted branches are indicative of a greater abundance of a taxon in the sample type listed in the row while taxa that are more abundant in the sample type listed in the column are highlighted in brown. Node sizes are representative of the number of ASVs detected per taxon. Water 3: particle-associated waterborne communities; Water 02: free-living waterborne communities.

Fig. S3
Relative proportions of bacterial orders detected for each sample type at each site. All orders with a mean relative abundance below 1% were grouped into one category (<1% Order). Helgoland (T) and (S) represent samples from the transferred and stationary cages, respectively. For particleassociated (Water 3) and free-living (Water 02) waterborne communities, Helgoland (T) and (S) refer to surface waters sampled from the final offshore site before and after the incubation of the transferred and stationary cages, respectively, which were pooled together.

Fig. S4
UpsetR plot displaying the total number of detected ASVs shared and unique to each sample type. The total number of ASVs detected for each sample type are shown in the horizontal bars to the left while vertical bars represent the number of ASVs shared or unique to each type highlighted by black nodes. Water 3: particle-associated waterborne communities; Water 02: free-living waterborne communities. Explained variation [%] along both Axes 1 and 2 are displayed in square brackets. Helgoland (T) and (S) represent samples from the transferred and stationary cages, respectively. For particle-associated (Water 3) and free-living (Water 02) waterborne communities, Helgoland (T) and (S) refer to surface waters sampled from the final offshore site before and after the incubation of the transferred and stationary cages, respectively, which were pooled together.  Helgoland (T) and (S) represent samples from the transferred and stationary cages, respectively. For particle-associated (Water 3) and free-living (Water 02) waterborne communities, Helgoland (T) and (S) refer to surface waters sampled from the final offshore site before and after the incubation of the transferred and stationary cages, respectively, which were pooled together.

Fig. S8
Mantel correlograms displaying significant phylogenetic signals detected across short evolutionary distances for each sample type. Filled black points are indicative of a significant relationship between the most closely-related taxa and their respective ecological niches calculated based on a combination of all environmental parameters measured (Temperature, Oxygen, Salinity, Conductivity, pH). Water 3: particle-associated waterborne communities; Water 02: free-living waterborne communities.

Fig. S9 Stacked bar plots showing the dissimilarities within A) HDPE, B) Tyre Wear, and C)
Wood communities observed between sites, partitioned into components of turnover (JTU) and nestedness (JNE) based on Jaccard dissimilarity.

Supplementary Tables
Table S1 Information on incubation sites, and samples collected and processed in this study. (A) Descriptions of the location and conditions at each site are presented along with sampling dates and times. A total of 99 samples (45 particle-associated; P01 -P45, 54 waterborne; W01 -W54) were analysed. (B) The number of particles or volume of water from which DNA was extracted, volumes of DNA sequenced, and their respective concentrations are provided for each sample. The table is appended as an Excel document.

Table S2
Mean relative abundances (%) of bacterial phyla detected per sample type. Helgoland (T) and (S) represent samples from the transferred and stationary cages, respectively. For particle-associated (Water 3) and free living (Water 02) waterborne communities, Helgoland (T) and (S) refer to surface waters sampled from the final offshore site before and after the incubation of the transferred and stationary cages, respectively, which were pooled together. The table is appended as an Excel document.

Table S3
Mean relative abundances (%) of bacterial classes detected per sample type. Helgoland (T) and (S) represent samples from the transferred and stationary cages, respectively. For particle-associated (Water 3) and free living (Water 02) waterborne communities, Helgoland (T) and (S) refer to surface waters sampled from the final offshore site before and after the incubation of the transferred and stationary cages, respectively, which were pooled together. The table is appended as an Excel document.   Table S6 Results of the Kruskal-Wallis Rank Sums test and pairwise Dunn's Test comparing the phylogenetic diversity of the different sample types based on Faith's PD. All p-values were adjusted using Benjamini-Hochberg correction (padj.). Significant values (p < 0.05) are highlighted in bold. TW: Tyre Wear; Water 3: particle-associated waterborne communities; Water 02: free-living waterborne communities.   Table S8 Results of a post hoc two-way PERMANOVA and PERMDISP comparing the different sample types collected from the transferred cage and surrounding waters at each sampling site based on taxonomic and phylogenetic dissimilarity metrics. All p-values were adjusted using Benjamini-Hochberg correction (padj.). Significance (p < 0.05) is highlighted in bold and additionally starred when corresponding pairwise PERMDISP p-values were also significant (*). TW: Tyre Wear; Water 3: particle-associated waterborne communities; Water 02: free-living waterborne communities; t: pseudo t-statistic. The table is provided as an Excel document.

Table S9
Results of a post hoc two-way PERMANOVA and PERMDISP comparing between the different sampling sites for each sample type collected from the transferred cage and surrounding waters based on taxonomic and phylogenetic dissimilarity metrics. All p-values were adjusted using Benjamini-Hochberg correction (padj.). Significance (p < 0.05) is highlighted in bold and additionally starred when corresponding pairwise PERMDISP p-values were also significant (*). TW: Tyre Wear; Water 3: particle-associated waterborne communities; Water 02: free-living waterborne communities; t: pseudo t-statistic. The table is provided as an Excel document.   Wood, and (D, E) waterborne communities between sites. The sum (∑) and mean (µ) percentage contribution of each bacterial family to differences observed between sites and their average dissimilarities (%) are shown. The number of ASVs (NASV) belonging to each family, their total and mean abundance at each site, and the difference in the mean counts (δ) detected between each pair of sites are also presented. Helgoland (T) and (S) represent samples from the transferred and stationary cages, respectively. For particle-associated (Water 3) and free-living (Water 02) waterborne communities, Helgoland (T) refers to surface waters sampled from the final offshore site before and after the incubation of the transferred cage which were pooled together. The table is provided as an Excel document.

Table S12
Results of a one-way PERMANOVA showing the effect of Sample Type on community structure differences between samples from the transferred and stationary cages based on taxonomic and phylogenetic dissimilarity metrics. p-values were adjusted using Benjamini-Hochberg correction (padj.) and significance (p < 0.05) highlighted in bold. d.f.: degrees of freedom; SS: sum of squares; Sq. root: square root.

Resemblance
Sources

Table S13
Results of a post hoc one-way PERMANOVA and PERMDISP comparing between samples collected from the transferred and stationary cages at the final offshore site (Helgoland only) based on taxonomic and phylogenetic dissimilarity metrics. All p-values were adjusted using Benjamini-Hochberg correction (padj.). Significance (p < 0.05) is highlighted in bold and additionally starred when corresponding pairwise PERMDISP p-values were also significant (*). (T) and (S) represent samples from the transferred and stationary cages, respectively. The table is provided as an Excel document.