Wastewater constituents impact bio ﬁ lm microbial community in receiving streams Science of the Total Environment

Swiss streams to test how bacterial and eukaryotic communities respond to wastewaterconstituents.Streambio ﬁ lmcompositionwasstrongly affectedbygeographiclocation – particularly forbacteria.However,theabundanceofcertainmicrobialcommunitymemberswasrelatedtomicropollutantsin the wastewater – among bacteria, micropollutant-associated members were found e.g. in Alphaproteobacteria , and among eukaryotes e.g. in Bacillariophyta (algal diatoms). This study corroborates several previously charac- terizedresponses(e.g. asseenindiatoms), butalsoreveals previouslyunknown community responses – suchas seenin Alphaproteobacteria .Thisstudyadvancesourunderstandingoftheecologicalimpactofthecurrentwaste- water treatment practices and provides information about potential new marker organisms to assess ecological change in stream bio ﬁ lms. © 2021 The Authors. Published by Elsevier B.V. This is an open access article under the CC BY license (http:// creativecommons.org/licenses/by/4.0/).


H I G H L I G H T S
• Stream biofilms are exposed to anthropogenic impacts from agriculture and WWTPs. • Biofilms were sequenced from 4 streams upstream and downstream of WWTPs. • Bacterial communities had a strong regional signal while eukaryotes had not. • Bacterial and eukaryotic communities shifted significantly in correlation with MPs.

• Shifted groups include diatoms and
Alphaproteobacteria, potential marker organisms.

G R A P H I C A L A B S T R A C T
a b s t r a c t a r t i c l e i n f o

Introduction
Microbial communities perform fundamentally important biological functions within organisms (i.e. host-associated microbiomes), and across ecosystems (Bardgett and van der Putten, 2014;Donaldson et al., 2015;Fuhrman et al., 2015;Grice and Segre, 2011;Guttman et al., 2014;Kwong and Moran, 2016). In stream ecosystems, microbial life is dominated by dense biofilms consisting of bacteria, archaea and microscopic eukaryotes. These biofilms form the basis of the stream food webs by performing crucial functions, such as elemental cycles, ecosystem respiration, primary production and maintenance of good water quality (Aristi et al., 2015;Battin et al., 2016). By being in contact with flowing water, stream biofilms are continuously in contact with various microorganisms and exposed to suspended nutrients and substrates. Stream biofilms are capable of resisting invasions by microorganisms to some degree and therefore have a central role in attenuating microbiological changes in stream water (Battin et al., 2007(Battin et al., , 2003Carles et al., 2021). Stream biofilms also take up substances that are present in the stream water, including nutrients and chemical pollutants (Flemming and Wingender, 2010;Singer et al., 2010), which can be toxic and/or be metabolized. Moreover, stream biofilm composition can be influenced by wastewater-born microorganisms (Carles et al., 2021). As a consequence, stream biofilms can be strongly affected by human activities such as agriculture, industry or wastewater dischargeresulting in potential alterations of the ecological functions they provide (Besemer et al., 2009).
One of the main anthropogenic influences on stream ecosystems is chemical pollution from excess nutrients and toxicants due to agricultural activities and discharged wastewater (Dudgeon et al., 2006;Schwarzenbach et al., 2006). Although wastewater treatment is aimed to assure clean water, by chemical and microbiological removal of toxic compounds and harmful microbes, wastewater discharge still contains hundreds of compounds (e.g., pharmaceuticals and pesticides), collectively known as micropollutants (Schwarzenbach et al., 2006;Stamm et al., 2016), and a large variety of microbes. These components can substantially influence the chemistry, community structure and ecosystem function of the receiving water bodies (Burdon et al., 2020Munz et al., 2017;Pascual-Benito et al., 2020;Schäfer et al., 2016;Stamm et al., 2016;Thompson et al., 2016). The mixture of micropollutants and microorganisms in wastewater can impact stream biofilm communities in complex ways, reflected for example as increased community tolerance and altered functions of algal and bacterial communities (e.g., photosynthesis and metabolic rates (Aristi et al., 2015;Corcoll et al., 2014;Tlili et al., 2017)).
Here, we apply next-generation DNA sequencing (NGS) of 16S (bacterial) and 18S (eukaryotic) rRNA genes to test how wastewater effluent affects biofilm community diversity in four Swiss streams. We also extend a recent study (Tlili et al., 2017) that found consistent changes in biofilm diversity, as well as increased tolerance to micropollutants of algal and bacterial constituents of biofilm communities, downstream of wastewater treatment plants (WWTPs). We here study the same sites and explicitly i) profile the taxonomic composition of the stream biofilm communities upstream (US) and downstream (DS) of WWTPs and ii) test whether any specific community-level or specific OTUlevel changes can be attributed to MPs or nutrients. This approach allows us to identify specific taxa responding to wastewater exposurethereby providing detailed insight into community level variation that underlie ecotoxicological and ecological responses to pollution.

Overall biofilm diversity
Sequencing of the 16S and 18S genes provided, after quality control, 2.8 million paired end 16S and 0.7 million paired end 18S reads. Phylogenetic classification of the OTUs (Greengenes for 16S, SILVA for 18S) revealed that these stream biofilm communities are most commonly dominated by the bacterial phyla Bacteroidetes and Proteobacteria, and the eukaryotic phyla Chlorophyta, Ochrophyta and Metazoa (Fig. 1).
For both bacteria and eukaryotes, alpha diversity measures Chao1 and Shannon were significantly different among streams (stream main effect: all p < 0.001; factorial ANOVA, Fig. S1). There was no significant overall difference between US/DS locations (US/DS location main effect: all p > 0.5; factorial ANOVA). However, significant stream × US/DS location interactions (for both bacterial and eukaryotic communities) indicated that the effect of wastewater input on alpha diversity was stream specific (all p < 0.001; factorial ANOVA, Fig. S1).

Biofilm community responses
To test the community differences across the streams and sampling locations, we performed principal coordinate analysis (PCoA) of the UniFrac distances of the proportional OTU counts across the samples. For bacteria, the UniFrac-PCoA-Axis 1 16S explained 77% and the UniFrac-PCoA-Axis 2 16S 7.2% of variation in multivariate community diversity ( Fig. 2A). For UniFrac-PCoA-Axis 1 16S there was a highly significant stream main effect, but no significant US/DS location or stream × US/DS location interaction (Fig. S2). Pairwise Tukey tests on UniFrac-PCoA-Axis 1 16S further showed that Steinach and Herisau differed strongly from the other two streams (Steinach-Buttisholz and Steinach-Hochdorf: both p < 0.001, Herisau-Hochdorf: p < 0.005 and Herisau-Buttisholz: p < 0.05; factorial ANOVA). The significance of the stream effect was further confirmed by permutational ANOVA (p < 0.001; 999 permutations). There were no significant effects of stream, US/DS location or their interaction on UniFrac-PCoA-Axis 2 16S . These results indicate that the multivariate diversity of bacterial communities in these stream biofilms is determined primarily by the source stream and are not detectably influenced by wastewater input.
For eukaryotes, UniFrac-PCoA-Axis 1 18S explained 50.7% and UniFrac-PCoA-Axis 2 18S 22.7% of variation. For both axes, there was a significant stream main effect (both axes p < 0.001; factorial ANOVA, Fig. 2B). The US/DS location main effect was not significant for either axis, but a significant stream × US/DS location interaction (both axes p < 0.001; factorial ANOVA) suggests that community shifts of eukaryotes (in response to wastewater input) are specific to each stream (i.e. indicating context dependency; Fig. S2). Visual inspection of the UniFrac-PCoA-Axes indicates no consistent patterns in Axis 1 18S , whereas inspection of Axis 2 18S indicates that DS locations become relatively similar but US locations remain distinct across streams (Fig. 2B).

Biofilm composition responses
To test whether micropollutants (MPs) in the wastewater effluent affected the biofilm composition, we performed redundancy analyses (RDA) using the proportional bacterial or eukaryotic OTU counts as response variables and nutrients and selected MPs (measured in all samples, including upstream and downstream of the WWTPS) as explanatory variables, providing 21 independent observations. Community changes in the stream biofilms were significantly associated with MPs for both bacteria and eukaryotes (RDA and permutational ANOVAs; p < 0.001), with the first two RDA axes explaining 56.5% of the total bacterial and 53.3% of the total eukaryotic variation (Fig. 3). In contrast, biofilm community changes were not significantly related to nutrient levels in the effluent -either for bacteria or eukaryotes (RDA and permutational ANOVAs p > 0.2; Fig. S03). For bacteria, the 10 most strongly affected OTUs, identified by RDA, included various unidentified species belonging to Rhodobacters as well as a Gemmata sp. and a Rhodococcus fascians. For eukaryotes, the 10 most strongly affected OTUs, identified by RDA, included various Diatom species, as well as a golden alga Chrysochaete britannica, an ulvophyte Pseudendoclonium sp. and a spirotrich Aspidisca sp.

Discussion
The extent of microbial community diversity has long been poorly understood because of the difficulties in characterizing large numbers of taxa. With the advent of next generation DNA sequencing, high microbial diversity has been revealed in various natural systems (Bardgett and van der Putten, 2014;Donaldson et al., 2015;Fuhrman et al., 2015), including stream biofilms (Battin et al., 2016;Nega et al., 2019;Peng et al., 2018). The challenge remains to establish the determinants of microbial diversity patterns in nature and how microbial communities may respond to environmental perturbations (Lee et al., 2021;Mansfeldt et al., 2020). Here we studied microbial diversity patterns of stream biofilms in four Swiss lowland streams and tested how bacterial and eukaryotic communities respond to wastewater, which contains microorganisms, micropollutants (MPs) and nutrients. We did find clear community shifts in bacteria and eukaryotes in response to MPs in the wastewater but observed no such responses to increased nutrient levels.

Bacterial diversity
For bacteria, we found a strong geographic community shift as the two westernmost streams (Hochdorf and Buttisholz) differed strongly from the two easternmost streams (Steinach and Herisau) (Fig. 2). This geographic pattern is broadly similar to the denaturating gradient gel electrophoresis (DGGE) -based community profiling on the same biofilms (Tlili et al., 2017). Our current study allowed in addition to identify OTUs within these community shifts. The bacterial OTUs found at our study sites belong to groups previously reported in association with stream water and stream biofilms, including Actinobacteria, Proteobacteria, Bacteroidetes and Planctomycetes (Besemer et al., 2012;Levi et al., 2017;Mansfeldt et al., 2020;Zeglin, 2015). For instance, we found certain Burkholderia (belonging to Betaproteobacteria) which have the potential to degrade xenobiotics (O'Sullivan and Mahenthiralingam, 2005), and Rhizobiales (belonging to Alphaproteobacteria) which are known for The relative presence of the most common phylogenetic classes is presented as pie charts for each study location. The gray section indicates OTUs that fall below the relative abundance of 0.005% (Bacteria) or 0.01% (Eukaryota).  their role in nitrogen fixation (Newton et al., 2011) (Fig. S3). None of the observed bacterial groups, however, differed consistently between US and DS sitesindicating no obvious responses to wastewater.

Eukaryote diversity
We found a high abundance of Ciliophora (protozoans) and some fungi -particularly Cryptomycota and Bacillariophyceae, corresponding to previous reports of stream biofilms (Dopheide et al., 2008;Heino et al., 2010;Levi et al., 2017). Ciliophora (ciliates) are grazers in stream biofilms and highly important in transferring nutrients to higher trophic levels (Dopheide et al., 2009). In general, ciliates occupy a wide range of ecological niches, respond rapidly to environmental changes and are therefore used as indicators of ecosystem health (Lear et al., 2011).
The most abundant Fungi in our biofilms were Cryptomycota (Baschien et al., 2008;Jones et al., 2011), which are common in aquatic systems (Livermore and Mattes, 2013). Whereas fungi are generally known as degraders of organic matter (Miura and Urabe, 2015), the scarcity of information about the ecology and lifestyles of the microscopic Cryptomycetes complicates the interpretation of their ecological role in stream biofilms and calls for further studies to identify their role in stream biofilms. We also found certain members of Chytridiomycota, a group which is known for their parasitic lifestyles, pathogenity for aquatic taxa as well as ability to degrade organic matter (Letcher et al., 2006).
Finally, we found a high abundance and diversity of Bacillariophyta (diatoms), which are common inhabitants of river and river biofilm communities (Battin et al., 2016). Diatoms are known for their photosynthetic lifestyle, and contribute to ecosystem functioning by exuding carbohydrates and amino acids (Battin et al., 2016). Like protists (above), diatoms are sensitive indicators of pollutants and classically used in ecosystem biomonitoring as well as ecotoxicological studies, often using their morphological characteristics for identification (Keck et al., 2016). We note that NGS-based characterization of diatoms has recently been shown to be comparable with the accuracy of morphological characterization (Apothéloz-Perret-Gentil et al., 2017).
In contrast to Besemer et al. (2012), who reported abundant taxonomic groups in biofilms to be omnipresent, we found clear geographic differences in the distribution of both bacterial and eukaryotic OTUs in stream biofilms. This is seen both as substantial differences among streams (particularly for bacteria) and US/DS locations (particularly for eukaryotes). Our results support, overall, the view that stream biofilm communities vary strongly locally (Battin et al., 2016). It is important to note that these differences only became obvious at the OTU level (Fig. 3) broader phylogenetic level distribution would be in agreement with (Besemer et al., 2012) (Fig. 1), which provides clear evidence for the benefits of NGS based estimations of biofilm diversity. The specific determinants of biofilm diversity across streams are unclear, but can reflect differences in catchment characteristics (e.g. catchment land-use; ), local climatic conditions (Fierer et al., 2007;Lear et al., 2014) or historical colonization patterns (Martiny et al., 2006).

Micropollutants are affecting bacterial and eukaryotic communities
We found no significant relationship between community shifts and dissolved phosphorus, nitrogen or organic carbon within streams or in response to effluent exposure. However, we found a significant correlation between micropollutants and a shift in the stream biofilm bacterial and eukaryotic communities. Among the 44 MPs (22) which were measured in the river water of our study sites, the strongest associations of bacterial communities, revealed by constrained ordination analysis, were with atrazin, chlortoluron, penconazole, pirimicarb, tebuconazole and terbuthylazine. The strongest associations of eukaryotic communities, revealed by constrained ordination analysis, were with azoxystrobin, chlortoluron, diuron, fexofenadine, simazine and terbuthylazine.
Interestingly, the herbicides terbuthylazine and chlortoluron correlated with both bacterial and eukaryotic community shifts in our analyses, which could indicate potential direct effects of herbizides on algae and perhaps indirect effects on bacteria (i.e. via altered trophic interactions).

Affected bacterial and eukaryotic taxa
We observed a significant bacterial community shift in response to wastewater exposure. In previous studies, several MPs have been shown to cause changes in biofilm communities. For instance, (Chonova et al., 2019;Kim Tiam et al., 2014) reported changes in the diatom community within stream biofilms in response to wastewater and to a mixture of fungicides, pesticides and herbicides, respectively, and molecular fingerprinting experiments (e.g. (Proia et al., 2013;Vercraene-Eairmal et al., 2010)) have demonstrated biofilm community changes in response to specific pollutants.
Using NGS analyses, we were able to identify specific members of the biofilm communities that were negatively related to the presence of chemical pollution. For example, there was a negative response of three distinct Rhodobacter species to herbicides chlortoluron, atrazine and terbutylazine as well as insecticide pirimicarb and fungicide tebuconazole. Since chlortoluron, atrazine and terbutylazine are known inhibitors of photosynthesis (Zhu et al., 2009), this effect may be mediated through the photosynthetic lifestyle of certain Rhodobacterial groups (Ravi et al., 2019;Strnad et al., 2010). While literature regarding the responses of Rhodobacters to MPs in stream biofilms is scarce, certain studies also indicate their usefulness in bioremediation of pesticidecontaminated wastewater (e.g. (Wu et al., 2019)). The communitylevel correlations with exposure to insecticide pirimicarb, and fungicides tebuconazole and penconazole, may reflect indirect effects on grazers or fungal communities  or unknown direct action on bacterial communities, but this finding currently remains unexplained due to lack of information.
We further observed that several diatoms had complex relationships to the MPs in the wastewater. For instance, the presence of diatom Achnanthidium minutissimum coincided with a high concentration of the herbicide diuron, whereas another diatom Cocconeis sp. exhibited an opposite response. This variability in the sensitivity of diatoms is consistent with the morphology-based assays of MP sensitivity (e.g. ). Since most diatoms are photosynthetic, the observed correlations may be mediated by photosynthesis inhibition by diuron, simazine, terbutylazine and chlortoluron. Other taxa apparently correlated with MP exposure in our data were an unidentified Pseudendoclonium species (a green algae in the class Ulvophyceae), an unidentified Aspidisca species (a bacterivorous ciliate in the class Spirotrichea) and Chrysochaete britannica (an algae in the class Chrysophyceae-Synurophyceae). The community-level responses observed in correlation with exposure to antihistamine fexofenadine, antidiabetic sitagliptin and fungicide azoxystrobin remain unexplained due to our lack of understanding of the mechanisms by which they might impact microscopic eukaryotes.

In conclusion
Our NGS data of four Swiss lowland streams showed that eukaryotic community members respond to wastewater input (particularly algal community members), whereas bacterial community composition reflects more strongly the geographic location. While these patterns reflect a combination of local and regional ecological processes typical for biofilms (e.g. (Besemer, 2015)), and likely a high turnover of bacterial communities, we also show that the biofilm community composition is related to observed patterns of chemical pollution. Most importantly, the NGS approach allowed us to identify specific taxa that seem to respond to wastewater in both bacteria and eukaryotes, suggesting the potential for new bioindicator taxa. Future work should explicitly test how functional endpoints of community level assays (e.g. enzymatic activity-based assays) correlate with taxa identity.

Study locations
We profiled the stream biofilm communities upstream (US) and downstream (DS) of WWTPs in four streams in central Switzerland (Buttisholz, Herisau, Hochdorf and Steinach, Fig. 1). The streams were initially chosen so that no wastewater was present upstream of the WWTPs and that the wastewater effluent contributed at least 20% to the total stream flow at the DS locations during low flow conditions . All streams were small to moderately sized. In each stream, one US and one DS location was selected as the reference and the impacted location, respectively. The US/DS locations were chosen to be as similar as possible with regard to stream morphology, riparian land use and vegetation Munz et al., 2017). Furthermore, the DS locations were selected so that water from the effluent was completely mixed with the stream water during low flow conditions.
Our samples originate from the same sampling campaign that recently documented Pollution induced tolerance (PICT) in these streams (Tlili et al., 2017). The details of the biofilm sampling can be found in Tlili et al. (2017) and we here provide only a summary. At each sampling location, glass slides (35.5 × 13.0 cm) were fixed vertically in perforated plastic boxes (three boxes per location and three glass slides per box) and immersed at the center of each stream to allow colonization by local biofilm. After a colonization period of 6 weeks (from 15th of March to the 30th of April 2014), the glass slides were retrieved, placed individually in plastic bags containing stream water from the corresponding sampling location and transported to the laboratory at Eawag, Dübendorf, Switzerland. The samples were transported in cooling boxes within 5 h of sampling. Immediately upon arrival at the laboratory, the stream biofilm from each glass slide was carefully scraped using a polypropylene spatula and suspended in 250 mL of Evian mineral water prior to functional and structural analyses of the stream biofilmsee (Tlili et al., 2017). The samples from each of the three glass slides within a box were pooled, resulting in three independent replicates/location. For genomic DNA extraction, 2 mL of the biofilm suspension within each biological replicate were centrifuged at 14000 ×g for 30 min at 4°C. Afterward, the supernatants were removed and the pellets stored at −80°C.

Micropollutant and nutrient data
We used micropollutant (MP) and nutrient data for each of the streams from Tlili et al. (2017) to characterize potential chemical exposure. This data included measurements of 44 MPs as well as dissolved organic carbon, dissolved nitrogen and dissolved phosphorus across our four study streams.

Molecular analyses
DNA extraction from the pellets was conducted using the PowerBiofilm DNA Isolation Kit (MO BIO Laboratories, CA) following the manufacturer's instructions. The isolated DNA was stored at −20°C until further processing. For the molecular genetic analyses, each biological replicate was further divided into three technical replicated to account for stochastic effects originating from PCR. However, some variation in replication resulted because of DNA recovery and technical issues. Therefore, only one biological replicate was available for the 16S for Buttisholz and Hochdorf at the DS locations and for the 18S for Buttisholz at the DS location.
After having assessed the optimal cycle number on an ABI 7500 cycler, the 16S and 18S regions were amplified resulting in 463 bp and 558 bp fragments for the 16S and 18S, respectively. The primer sequences are listed in Table S01. The correct cycle number for each sample was determined using qPCR with reaction conditions for 16S and 18S as follows: 1× Phusion HF buffer (catalogue no. M0530S), 0.2 mM dNTPs (Promega catalogue no. U1515), 0.4 μM of each primer, 1× EvaGreen (catalogue no 31000-T, 31000) and 0.5 U Phusion polymerase (catalogue no. M0530S). Thermocycling conditions consisted of an initial denaturation of 98°C for 30 s; 25 cycles of 98°C for 5 s, 60°C for 30 s, 72°C for 30 s. The reactions were performed in three technical replicates for 16S amplicons as follows: 1× GoTaq G2 (catalogue no. M7832), 0.2 mM dNTPs (Promega catalogue no. U1515) and 0.4 μM of each primer. Thermocycling conditions consisted of an initial denaturation of 95°C for 120 s; 20 cycles of 95°C for 30 s, 55°C for 45 s, 72°C for 50 s. The reactions were performed in three technical replicates for 18S amplicons as follows: 1× Phusion HF buffer (catalogue no. M0530S), 0.2 mM dNTPs (Promega catalogue no. U1515), 0.4 μM of each primer and 0.5 U Phusion polymerase (catalogue no. M0530S). Thermocycling conditions consisted of an initial denaturation of 98°C for 30 s; 23 cycles of 98°C for 5 s, 60°C for 30 s, 72°C for 30 s, with the exception of Steinach US for which, because of low DNA yield, a cycle number of 28 was used.
The three technical replicates were re-combined and subsequently run on an agarose gel to confirm expected amplification or nonamplification (for negative controls in which PCR templates were omitted). PCR products were cleaned using AMPure XP beads in a ratio of 0.8× per sample. The samples were indexed using the standard Illumina TruSeq indices and sequenced on an Illumina® MiSeq machine in a single flow cell in paired-end mode with 300 bp-read-lengths at the Genetic Diversity Center (GDC), ETH Zürich. The sequences were deposited in the National Center for Biotechnology Information (NCBI) Sequence Read Archive database under BioProject ID PRJNA414710.

Sequence processing
The raw data was processed through the UPARSE pipeline (http:// drive5.com/usearch/manual/uparse_pipeline.html; (Edgar, 2013)). For the 16S sequences, the Illumina paired reads pipeline was used. For 18S sequences, the Illumina forward reads pipeline was used because only approximately 10% of the reads remained after merging. The reads were filtered by discarding reads with a total expected error > 1 for all bases. Afterwards the sequences were dereplicated using fulllength matching. Finally, the operational taxonomic units (OTU's, (Edgar, 2013)) were clustered with a separation limit of 0.97 sequence similarity, and the chimeras filtered. The 97% sequence similarity cut off was used to reflect approximate species level divergence for microorganisms (Huse et al., 2010).
The taxonomic affiliations of the OTUs were assigned using the assign_taxonomy.py script from QIIME (www.qiime.org; (Caporaso et al., 2010)) with the Greengenes database (release 13_5; (DeSantis et al., 2006)) for the 16S sequences, and the SILVA database (release 123; (Quast et al., 2013)) for the 18S sequences. The sequences were aligned using the align_seqs.py script from QIIME with the PyNAST method (Caporaso et al., 2010) using the Greengenes, Silva and alignment databases for 16S and 18S, respectively. Phylogenetic relationships between the OTUs were determined by inferring a phylogenetic tree using FastTree (Price et al., 2010). For the 16S data, the phylum Parvarchaeota (from the kingdom of the Archaea) and for the 18S data, Homo sapiens, were used to root the tree, respectively.

Statistical analyses
All computational analyses were conducted separately for 16S and 18S data in R version 3.2.4 (R language core team 2008). We used multiple approaches to, on one hand, characterize taxonomic diversity overall and, on the other, to test the extent and context dependency of differentiation between US and DS locations. First, alpha diversity at each sampling site was estimated using Chao1 and Shannon indices, which account for rare taxa and overall number of taxa, respectively (Chao and Shen, 2003;Shannon, 1948). To gain multivariate descriptors of community change across sampling sites, we calculated weighted UniFrac distances based on OTU phylogeny and Bray-Curtis dissimilarities of the OTU tables (containing OTU counts in each sample), and conducted principal coordinate analyses (PCoA). This was performed for all eukaryotic and bacterial OTUs as well as for individual eukaryotic and bacterial phyla and classes. These analyses were performed using functions embedded in the phyloseq framework (McMurdie et al., 2013) on R. The MPs causing the strongest effects on bacterial and eukaryotic communities were identified using ordistep function embedded in the vegan package (Oksanen et al., 2019) in R. The RDA analyses for MP and nutrient effects were calculated using capscale function embedded in the vegan package on R using the bacterial or eukaryotic OTU counts as response variables and nutrients and selected MPs, measured in all samples, as explanatory variables.

Declaration of competing interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.