Biomass addition alters community assembly in ultrafiltration membrane biofilms

Freshwater biofilms assemble from a pool of rare water column genotypes. Random density fluctuations and temporal species turnover of functionally equivalent potential colonizers result in compositional variability of newly formed biofilm communities. We hypothesized that stronger environmental filtering as induced by enhanced substrate levels might reduce the impact of a locally variable pool of colonizers and instead select for more universal habitat specialists. Our model were heterotrophic biofilms that form on membranes during gravity-driven ultrafiltration of lake water. In four separate experiments, biomass of the cyanobacterium Microcystis was added to the feed water of one set of treatments (BM) and the resulting biofilm communities were compared to unamended controls (CTRL). Biomass addition led to a significant shift of community assembly processes: Replicate BM biofilms were more similar to each other than by chance in 3 of 4 experiments, whereas the opposite was the case for CTRL communities. Moreover, BM communities were more stochastically assembled across experiments from a common ‘regional’ pool of biofilm colonizers, whereas the composition of CTRL communities was mainly determined by experiment-specific ‘local’ genotypes. Interestingly, community assembly processes were also related to both, physiology (aerobic vs. anaerobic lifestyle) and the phylogenetic affiliation of biofilm bacteria.

Scientific RepoRtS | (2020) 10:11552 | https://doi.org/10.1038/s41598-020-68460-x www.nature.com/scientificreports/ provided substrates may also shift the respective importance of the 'local' and 'regional' pools of potential colonizers during biofilm community assembly, e.g., by affecting functional redundancy and the balance between competitors in the source assemblage 25 . We addressed this hypothesis using biofilms in experimental devices designed to investigate the potential of gravity-driven membrane filtration (GDM) for decentralized drinking water production 26,27 . We synoptically analysed microbial community composition and assembly processes in four independent experiments that were performed over the course of several years [27][28][29] with feed water from a single lake.

Results
Diversity and composition of biofilm communities. The complete sequencing dataset comprised 3,495 OTUs and approximately 10 6 reads. It was rarefied to the sample with the lowest read number (1.34 × 10 4 ) by random read exclusion 27 . The normalized data set used for all subsequent analyses consisted of 2,721 OTUs with a total of 0.21 × 10 6 reads. The CTRL treatment had significantly higher OTU richness (734 ± 165, mean ± standard deviation) than the BM treatment (357 ± 63) (Student's t test, n = 8, p < 0.001). The two treatments shared 786 genotypes that together formed 87% of total reads. Sixteen OTUs, representing almost one quarter of total reads, were present in every single community (Table 1).
In general, only very few OTUs on GDM biofilms were affiliated to the genotypes that are typically found in freshwater bacterioplankton 30 . Proteobacteria represented 51% of the total amount of OTUs in biofilm communities and 61% of the total amount of reads (Fig. 1). However, many OTUs within all phylogenetic groups of biofilm bacteria were treatment-specific, i.e., were either absent or had considerably lower read numbers in one of the treatments. A total of 63 OTUs significantly contributed to the distance-based separation of the BM and CTRL treatments (SIMPER analysis, p < 0.05). The numbers in brackets after the colour codes of the individual phylogenetic groups in Fig. 1 report their relative distribution. The majority of the treatment-discriminating OTUs were affiliated with Betaproteobacteria (31%), Alphaproteobacteria (24%), and Bacteroidetes (20%).
Community structure and assembly processes. More than 40% of OTUs from either treatment only occurred in a single biofilm community (Fig. 2). By contrast, only 49 and 75 OTUs were present in all eight samples of the BM and CTRL treatments, respectively. The subset of 'single sample' OTUs had significantly more reads in the CTRL than in the BM treatments, whereas the opposite was the case for those OTUs that occurred in all communities of a treatment type (Student's t tests, n = 16, p < 0.01).
Average linkage clustering of Bray-Curtis dissimilarities indicated a general separation between treatment types (Fig. 3), i.e., the communities of the BM treatments from all 4 experiments were more similar to each other than to the corresponding samples from the CTRL treatment. Similarity profile analysis showed that all phylogenetic groups except Firmicutes significantly contributed to this separation, albeit to a variable extent (Fig. 1). The clustering of only those OTUs that occurred in all biofilms (Table 1) also resulted in the perfect separation of samples according to treatment type (Bray-Curtis dissimilarity: 78%, data not shown). Five of these 16 'universal colonizers' significantly contributed to this separation, all of them with higher read numbers in the BM treatment ( Table 1).
The separation of communities according to treatment type was further corroborated by an error-free classification by Random Forest (RF) analysis (Out-of-Bag Error (OOB) = 0%, number of trees = 1,000, number of samples = 16, categories = 2). However, a RF classification of samples according to experiment was also remarkably accurate, and only 1 out of 16 samples was misclassified (OOB = 6.25%, number of trees = 1,000, number of samples = 16, categories = 4). The 20 OTUs that were most relevant for classification according to treatment type and experiment are listed in suppl. Tables S1 and S2, respectively.
Biological replicates of the BM treatments had significantly lower beta diversity (dissimilarity values: 0.39 ± 0.07) than those of the CTRL treatment (0.51 ± 0.07; one-tailed Student's t-test, n = 8, p < 0.05) (Fig. 3). Moreover, Raup-Crick indices showed that the biological replicates of the BM treatments were significantly more similar to each other than by chance (RC > 0.95) in 3 out of 4 experiments, whereas the opposite was the case for replicates of the CTRL treatments (RC < -0.95).
Null model analysis of community assembly processes indicated that the biofilms of the CTRL treatments were predominantly shaped by the 'local' pool of OTUs from the corresponding experiments (NST: 26%). By contrast, a significantly higher (p < 0.001) influence of the 'regional' pool of OTUs (i.e., that were present in several or all experiments) was found for biofilms of the BM treatments (NST: 66%). A graphic interpretation of these results is given in Fig. 4. In addition, assembly processes within the two treatments were also assessed for each of the phylogenetic lineages depicted in Fig. 1. Striking differences in treatment-specific NST values to those of the total community were observed for some phylogenetic groups, most prominently for Bacteroidetes and Firmicutes (Fig. 4, suppl. Table S3).
Aerobic and anaerobic subcommunities. 39% of all OTUs representing 78% of total reads could be assigned to either aerobic or anaerobic taxa, based on information on their closest cultured relatives. If possible, anaerobic taxa were further split into strict (272 OTUs) and facultative (62 OTUs) types. However, since both groups behaved similarly in the below described analyses, this distinction was abandoned for the sake of increased statistical power.

Discussion
It is still challenging to accurately quantify the 'rare biosphere' of potential biofilm colonizers in the bacterioplankton 31,32 and, thus, to study biofilm assembly processes by direct comparison with the source community. However, other approaches help to identify factors that modulate the respective importance of species sorting, immigration and stochasticity in biofilms, e.g., the analysis of communities in the context of the biotic and abiotic variation of the surrounding water 33,34 , or in their response to habitat manipulation 35 . The diversity (richness) of the GDM biofilm communities was significantly lower in the substrate amended BM treatments than in the CTRL. This is the opposite of what has been described for planktonic diversity in oligotrophic ponds 36 , or groundwater microbial assemblages 35 , but agrees with observations on microbial assemblages in eutrophied lake sediments or from a marine phytoplankton bloom 37,38 . Such seemingly contradictory findings might potentially be reconciled by assuming a hump-shaped relationship between microbial diversity and productivity, as observed in soil 39 : On the one hand, severe substrate limitation will impose restrictions on the survival of many genotypes and select for a small set of oligocarbophilic specialists. On the other hand, extremely high substrate levels will specifically favour the most rapidly growing 'opportunistic' genotypes, which in turn might negatively affect the growth of others by releasing waste products or toxins.
BM treatments also featured a significantly higher NST score than the CTRL. At a first glance, it appears paradoxical that these communities should be mainly shaped by neutral processes such as immigration and ecological drift at arguably more selective conditions. However, null-model-based approaches are very sensitive to regional pool pre-selection: Biofilm communities consist of a subset of the bacterioplankton source assemblage that is already strongly filtered by deterministic processes 4 . Our analysis thus assessed if the individual communities tended to be formed by bacteria that were able to colonize the biofilms and at the same time were specific to individual experiments ('deterministic' with respect to the impact of the local species pool) or by such genotypes that were present in all source assemblages ('stochastic' , i.e., mainly formed by the regional species pool) (Fig. 4). www.nature.com/scientificreports/ Thus, the high NST value of the BM treatments indicates that colonization of these biofilms was largely unrelated to priority effects, the compositional variability or the species turnover in the source water between experiments, and mainly driven by niche-related processes. The importance of environmental filtering in BM treatments is also suggested by the significantly higher similarity of the biological replicates within single experiment, i.e., lower Bray-Curtis distances (Fig. 3) and high Raup-Crick indices. Our findings agree with observations from parallel bioreactors, where higher substrate supply led to the formation of less divergent communities 40 . Competitive exclusion might be directly related to the higher availability of organic C and nutrients, e.g., by favouring more opportunistically growing genotypes 41 . Eutrophication and high organic load increased the importance of niche-driven selection in coastal bacterioplankton assemblages 42 . However, two indirect effects of biomass addition also need to be considered: For one, more extreme habitat conditions, in particular the more pronounced oxygen limitation due to higher respiration 43,44 , might also have favoured niche-driven assembly processes in BM treatments. Environmental filtering was most prominent at conditions that most strongly deviated from those of the source assemblages (i.e., draught) in the sediments of a salt pan 45 . Secondly, due to the added biomass the transmembrane flux in the BM treatments was on average between 2-3 times lower than in the CTRL treatments [27][28][29]43 . Thus, the observed assembly patterns might at least in parts be due to 'mass effects' , i.e., the successful colonizers of BM biofilms from the regional species pool (several of which are listed in Table 1) might have simply been more abundant in the source water than their competitors from the local pool.
The respective importance of assembly processes such as immigration and local selection is usually assessed at the level of entire microbial communities, implicitly assuming that these processes indiscriminately act on different community components. This has, e.g., been challenged by a separate assessment of the selection forces acting on common vs. rare members of microbial assemblages 31,32,46 . Biofilms are spatially heterogeneous habitats, featuring a complex three-dimensional architecture with pronounced horizontal and vertical gradients of oxygen or substrates 43,47 . It is conceivable that there might be differences in the assembly processes of distinct subcommunities that are either limited to particular microniches, such as obligate anaerobic bacteria 48 , or that share fundamental life style related traits 49 .
We found that niche-driven and neutral processes may selectively act on different compartments of GDM biofilm assemblages: In contrast to their aerobic counterparts, OTUs with an anaerobic metabolism were not selected from the local (i.e., experiment-specific) species pool in the CTRL treatment. Following the argument made above, we conclude that their establishment in GDM biofilms was thus mainly governed by niche-driven processes. Moreover, the clear difference in the composition of the anaerobic subcommunities between the two treatments suggests that anaerobic genotypes originated from several sources with different habitat properties, e.g. from sinking 'lake snow' particles 50 vs. from food pellets within zooplankton carcasses 51 . The contrasting www.nature.com/scientificreports/ availability of organic C in the two treatments might have been a direct selective factor, yet it probably also shaped biofilm communities indirectly by defining the spatial patterns of oxygen distribution: On the one hand, interspersed anoxic microniches within a largely oxygenated matrix would lead to high habitat heterogeneity and variable selection in the CTRL biofilms 44 . This agrees with the more than twice as high number of anaerobic OTUs that were exclusive to this treatment (159 vs. 70 OTUs). Moreover, 98% of the 605 reads affiliated to the Nitrosomonadaceae (microaerophilic lithoautotrophic ammonia oxidizers) originated from the CTRL treatment, suggesting an interplay of oxygen gradients and C limitation. On the other hand, the high loads of organic material in the BM treatments likely led to considerable lower redox conditions in large areas of these biofilms. This conclusion is supported by the observation that > 90% of reads from OTUs affiliated with sulphate reducing bacteria (total: 1536 reads) and with the strictly anaerobic Clostridiales (total: 6,413 reads) were found in this treatment. BM biofilms also featured tenfold higher read numbers of the most common 'globally occurring' OTU ( Table 1) that is closely related to the facultative iron reducer Rhodoferax ferrireducens (NR_074760, 98.8% sequence identity). In addition, there was evidence that microbes from individual phylogenetic groups did not adhere to the overall pattern of community assembly (Suppl . Table S3). In particular, OTUs affiliated with Bacteroidetes always tended to be recruited from the local species pool irrespective of treatment type, whereas the opposite was the case for Firmicutes (mainly consisting of Clostridiales classified as anaerobes). It has been argued that higher bacterial taxonomic ranks may to some extent also be ecologically meaningful units 52 , possibly because complex life style traits encoded by many genes are also phylogenetically conserved 49 . The addition of cyanobacterial biomass led to a significant shift towards a small subset of 'globally occurring' treatment-specific colonizers of GDM biofilms (Fig. 2) both, within and between experiments. Moreover, the treatment types could be unambiguously separated by only considering 5 ubiquitous genotypes that all had 5 to 10 times higher read numbers in the BM treatment (Table 1). Our observation of a few 'universally occurring' and without (CTRL) addition of cyanobacterial biomass, as deduced by null model analysis (Normalized Stochasticity Ratio). The 'regional pool' of biofilm colonizers is defined as those species that were present in the feed water of all experiments, whereas the 'local pool' refers to genotypes that were specific to the feed water of a single experiment.
Scientific RepoRtS | (2020) 10:11552 | https://doi.org/10.1038/s41598-020-68460-x www.nature.com/scientificreports/ indicator species for the addition of cyanobacterial biomass to GDM biofilms is derived from experiments conducted in 2011, 2014 and 2017, and thus appears to be rather robust with respect to variability of the source assemblages. Moreover, all 5 of the above mentioned OTUs were also present in biomass-amended GDM biofilms fed with water from a stream 28 . Two of them, Rhodoferax sp (LN870966) and Undibacterium seohonense (KC735151), also proliferated in GDM biofilms after addition of starch 43 , suggesting that their role as indicators might go beyond a particular lake and carbon source. Metabarcoding of pro-and eukaryotic marker genes from different habitats -including biofilms 53 -is increasingly regarded as a powerful tool for future environmental monitoring 54,55 . Recently, Keeley et al. reported that a small number of abundant microbial genotypes within operationally defined 'eco-groups' could accurately predict the enrichment of benthic habitats in salmon farms 56 . It is conceivable that bacterial biofilms in various technical or natural systems might also feature a small set of omnipresent reliable meters for fundamental aspects of habitat conditions such as substrate levels. Pronounced cyanobacterial blooms are a frequent phenomenon in many eutrophied lacustrine systems, resulting in concentrations of particulate organic carbon that are in the range of our experimental addition in the BM treatments 57 .
Our study might thus help to develop a better understanding of the microbial communities that develop on ultrafiltration membranes or in sand filters fed with water from such systems.

Materials and methods
Experimental systems. Microbial biofilms grew on the ultrafiltration membranes (150 kDa nominal cutoff, polyethersulfone membrane; Microdyn Nadir, Wiesbaden, Germany) of experimental gravity-driven water filtration (GDM) systems. Water from Lake Zurich, Switzerland (location: 47°19′13.24″ N, 8°33′11.86″ E) obtained from the aerobic layer at 5 m depth was used as continuous feed in four independent experiments conducted over a period of 6 years ( Table 2) [27][28][29] . Lake water was collected in sedimentation tanks and the GDM systems were supplied from the overflow. More details and a graphic depiction of the experimental systems are given in Silva et al. 27 and 28 .
In each experiment, two treatments were run in parallel with the same feed water: One set of GDMs (2-3 biological replicates per experiment) received a daily dose of destroyed biomass from an axenic culture of Microcystis aeruginosa PCC7806 (BM treatment, Table 2). Another set of GDMs was operated without additional manipulation (CTRL treatment) until biofilm maturation (2-3 weeks). The input of particulate organic carbon from lake water to the CTRL treatments was by approximately 2 orders of magnitude smaller than that of the BM treatments, as estimated from the concentrations of particulate organic matter (POC) in Lake Zurich 58 and daily flux rates through GDM systems 27,28,43 . DNA extraction and sequencing. At the end of each experiment, the biofilms were collected and DNA was extracted using the DNeasy PowerBiofilm Kit (Qiagen, Germany) according to the manufacturer's specification (but extending the removal step for inhibitors to one hour). The recovered DNA was dissolved in 10 mM Tris buffer and stored at -20 °C until further processing.
Partial 16S rRNA genes were amplified with a primer pair that excludes chloroplasts and cyanobacteria (799F-1115R) 59,60 , but using a modified reverse primer for increased coverage 27 . Sequences were obtained by Illumina MiSeq paired end (2 × 300 bp) sequencing. Sequencing data from two of the four experiments (Exp 2, 4) have been published before, albeit separately 27,28 . The stored DNA extracts from another published experiment (Exp 1) 29 were reamplified and resequenced on the Illumina platform for the purpose of this study. Samples from all experiments were processed by the same company (LGC Genomics, Germany).
Amplicon sequence data from two biological replicates of each treatment type from all four experiments were collectively re-analysed in order to produce a single, coherent set of operational taxonomic units (OTUs). The raw sequence data processing, the definition of OTUs and their taxonomic assignment according to the SILVA taxonomy (version 132) 61 were performed by an in-house pipeline as described previously 27 .
Data treatment and statistical analysis. Data analysis was carried out in R 62 . Biofilm communities were clustered by the average linkage method based on their Bray-Curtis dissimilarity. This was done for the complete dataset and separately for each of the phylogenetic groups depicted in Fig. 1. Similarity profile analysis was used to establish significant clusters (1,000 simulations, α = 0.001) and the stability of the main branches was estimated by bootstrapping (1,000 bootstraps). These calculations were performed with the R packages clustsig, vegan, and fpc [63][64][65] . Similarity percentage (Simper) analysis based on Bray-Curtis dissimilarity (R package vegan) was used to examine which of the 'core' genotypes that were present in all samples (Table 2) significantly contributed to the separation between treatments. www.nature.com/scientificreports/ Differences in various community parameters between the BM and CTRL treatments were tested for significance by Student's t-tests, following a prior testing for normality (Shapiro-Wilk test) and homoscedasticity (Levene's test) with the R package car [66][67][68] . Proportional data were logit transformed prior to statistical testing 69 . In order to account for a possible confounding effect of differences in alpha diversity between treatments 70 we also calculated the pairwise Raup-Crick indices (RC) of the biological replicates using the R package NST 71 .
A Random Forest (RF) classification (R package randomForest) was performed to assess how accurately the samples could be assigned to treatments or experiments by means of their respective community composition, and to identify the responsible OTUs 72,73 . The RF algorithm is a supervised machine learning model based on decision trees to classify data into pre-defined categories. The number of required decision trees (range: 1,000 to 50,000) was set to a value where the out-of-bag (OOB) error was stable. The Gini Impurity Metric was used as criterion to identify OTUs that were most relevant for the classification.
The influence of the 'local' and 'regional' pools of OTUs (i.e., OTUs that tended to be more or less specific for a single experiment or treatment) on biofilm community assembly in the two treatments was estimated by a null model analysis using the Normalized Stochasticity Ratio (NST; R package NST) 71 . NST values > 50% generally point to a predominance of stochastic processes on community assembly 71 . In our specific context, high NST values speak for the importance of biofilm colonizers from the 'regional pool' of OTUs that were present in similar read numbers across several or all experiments. Differences between NST values were assessed for statistical significance by bootstrap analyses (1,000 bootstraps).