Thickness determines microbial community structure and function in nitrifying biofilms via deterministic assembly

Microbial biofilms are ubiquitous in aquatic environments where they provide important ecosystem functions. A key property believed to influence the community structure and function of biofilms is thickness. However, since biofilm thickness is inextricably linked to external factors such as water flow, temperature, development age and nutrient conditions, its importance is difficult to quantify. Here, we designed an experimental system in a wastewater treatment plant whereby nitrifying biofilms with different thicknesses (50 or 400 µm) were grown in a single reactor, and thus subjected to identical external conditions. The 50 and 400 µm biofilm communities were significantly different. This beta-diversity between biofilms of different thickness was primarily caused by deterministic factors. Turnover (species replacement) contributed more than nestedness (species loss) to the beta-diversity, i.e. the 50 µm communities were not simply a subset of the 400 µm communities. Moreover, the two communities differed in the composition of nitrogen-transforming bacteria and in nitrogen transformation rates. The study illustrates that biofilm thickness alone is a key driver for community composition and ecosystem function, which has implications for biotechnological applications and for our general understanding of biofilm ecology.

www.nature.com/scientificreports www.nature.com/scientificreports/ surface-area 36 , they are expected to have higher species richness than thin biofilms, thus they are by chance expected to be colonized by more species; thicker biofilms would also have larger gradients of substrates and electron acceptors (e.g. due to anaerobic zones), that could allow the establishment of anaerobic taxa, increasing species richness compared to the thin biofilm.
We hypothesize that the 50 µm biofilm community would be a subset of the richer 400 µm biofilm community, due to anaerobic taxa being restricted to the thicker 400 µm biofilms, while the same aerobic taxa would occur in both the 50 and 400 µm biofilms. Alternatively, turnover, i.e. differences in species identity, could arise if biofilms of different thickness have different environmental conditions apart from gradients of substrates and electron acceptors, or due to ecological drift. We used the Baselga framework 41 to estimate how much turnover and nestedness contributed to the observed beta-diversity, and a null model approach was used to determine the importance of deterministic versus stochastic assembly processes 42 . Finally, in order to link possible differences in taxa richness (i.e. nestedness) or identity (i.e. turnover) to functional differences between biofilms of distinct thicknesses, we also measured rates of nitrogen transformations, and discuss the implications of the results for wastewater treatment.

Results
Two different biofilms. We grew biofilm communities with a maximum thickness of 50 μm and 400 μm together in the same bioreactor; these communities are referred to as Z50 and Z400. CLSM images of EPS stained biofilm cryosections confirmed that carrier design limited biofilm thickness ( Fig. 1a and b).
We used null modelling to estimate the standardized effect size (SES) for β sor . We observed that β sor values for between-group comparisons, i.e. between Z50 and Z400, were higher than expected by chance (SES βsor > +2) (Fig. 3a) indicating that between-group differences were likely deterministic. On the contrary, observed β sor values for within-group comparisons, i.e. between the carriers of the same type, were not more different than expected by chance (|SES βsor | < 2) (Fig. 3a). In addition, estimation of the quantitative RC BRAY metric, also indicated that Z50 and Z400 communities were in average more dissimilar than the null expectation (between-group RC BRAY > +0.95).
To determine the contributions of nestedness and turnover to beta-diversity, we estimated the two components of β sor : β sim (turnover) and β sne (dissimilarity due to nestedness) using the Baselga framework 41 (Fig. 3b) and calculated the ratio between β sne and β sor , referred to as the β ratio 44 . When the β ratio is smaller than 0.5, beta-diversity  (b) β sor , β sne (dissimilarity due to nestedness) and β sim (turnover) values; the sum of β sim and β sne is β sor . (c) Beta diversity ratio. Values were estimated for pairwise comparisons among Z400 replicates (n = 10), Z50 replicates (n = 10) and between the two groups.
To investigate which taxa were contributing to the differences in community composition between Z50 and Z400, differences in relative abundance of taxa between Z50 and Z400 were estimated using DESeq2. We found differential abundance (p (adj) < 0.01, DESeq2) for 45% of the sequence variants (SVs) analyzed with DESeq, while for the top 40 most abundant SVs, 32 had different abundance between Z50 and Z400 ( Supplementary  Fig. S1). Among the fraction with differential abundance, 26% of SVs were more abundant in Z50, and 74% were more abundant in Z400. The effect of thickness on relative abundance, if any, differed among taxa (for example Supplementary see Fig. S1).
Between-group sorting of nitrifiers and anammox bacteria. The relative read abundance of the nitrifiers, Nitrosomonas, Nitrospira and Nitrotoga, was lower in the Z400 biofilms with Nitrotoga being almost restricted to Z50 (Fig. 4a). The same trends were noticed using quantitative fluorescence in situ hybridization (qFISH; Fig. 4c; Welch's t-test, p < 0.05). It was not possible to detect by qPCR if comammox were present due to non-specific amplification using Nitrospira amoA primers 45 .
Interestingly, we observed the anammox bacterium Brocadia in the Z400 biofilms, but it was almost absent in the Z50 biofilms (Fig. 4a). This was supported by qFISH (Welch t-test p < 0.001) (Fig. 4c). Sorting of bacteria between thick and thin biofilms was not only limited to primary producers (i.e. autotrophic nitrogen converters) but also seen among the predatory Bdellovibrionales. Bacteriovorax had a higher abundance in the Z50 communities, while some SVs classified as Bdellovibrio were more abundant in either Z400 or Z50 ( Supplementary Fig. S2).
FISH analyses of biofilm cryosections showed that the Z400 biofilm was likely stratified, e.g. with Nitrospira being more abundant in the middle of the biofilm and the anaerobic anammox bacteria being present in the deeper layers; Nitrosomonas biovolume was the same along the depth gradient ( Fig. 5a and b), but the biovolume fraction decreased with depth ( Supplementary Fig. S3). In the thin Z50 biofilms, no stratification was observed as the AOB and NOB populations were located side by side (Fig. 5c).
Dissolved oxygen (DO) concentration profiles ( Fig. 6) were calculated for one-dimensional biofilms having the same average community composition, density gradient, and nitrogen transformation rates as the real Z50 and Z400 biofilms. The results give a range of possible DO concentration profiles, which are shown as shaded regions (Fig. 6). The model predicts that Z50 biofilms can be fully oxygenated but may also have anoxic regions, whereas the Z400 biofilms contain a completely anoxic region in its deeper parts in all tested scenarios.
Nitrogen transformation rates. Two types of tests were performed separately on the Z50 and Z400 carriers; (i) actual activities tested in a continuous laboratory trial, with the same incoming water as in the 0.5 m 3 reactor and (ii) potential activities tested in batch trials where excess nitrogen was added. For all trials, removal rates are reported per surface area and day. Actual rates of net NO 3 − production were 1.4-1.5 gNO 3 − N/m 2 , d for Z50 and 0.68-0.72 gNO 3 − /m 2 , d for Z400. To estimate NO 3 − production per nitrifier abundance, it is necessary to consider differences in biomass between carriers. We estimated that the total nitrifier biomass per carrier surface was about the same in Z50 and Z400 (Fig. 4b). Therefore, per nitrifier biomass, net NO 3 − production was higher in Z50 than in Z400.
In the aerobic potential tests for net NH 4 + removal (Fig. 7a), net NO 3 − and net NO 2 − production (per carrier area) was higher for Z50 than Z400 biofilms (ANCOVA, p < 0.05), while the rate of net NH 4 + removal was not significantly different between Z50 and Z400 (ANCOVA, p > 0.05). The aerobic potential removal of NO 2 − (Fig. 7b) was significantly higher for Z400 than for Z50 (ANCOVA, p < 0.05). Finally, in the anoxic potential trials, in which NH 4 + and NO 2 − were added simultaneously (Fig. 7c), removal of NO 2 − was significantly higher for Z400 than for Z50 (ANCOVA, p < 0.05), while no significant removal or production of NH 4 + was seen for either Z50 or Z400.

Discussion
Although incubated in the same bioreactor and experiencing the same conditions and the same history, different microbial communities developed on carriers with thin and thick biofilms (Fig. 2b). The thicker Z400 biofilm had a higher richness and evenness than the thinner Z50 biofilm (Fig. 2a) and our results are therefore in agreement with known positive species-area relationships for microbial communities 46 . Moreover, similar to our results, Torresi et al. 38 , focusing on micro-pollutant degradation, also found a significant higher evenness in thicker biofilms.
A null model approach was used to investigate if the differences in beta-diversity between Z50 and Z400 were due to deterministic or stochastic processes while accounting for the large differences in richness between Z50 and Z400 43 . The results showed that the between-group beta-diversity was higher than expected by chance (Fig. 3a), suggesting deterministic assembly due to differences in biofilm thickness. Furthermore, biofilm thickness significantly affected the relative abundance of the majority of the most abundant individual taxa, meaning that they showed clear preference for either thin or thick biofilms ( Supplementary Fig. S1). Some turnover among the Z50 and Z400 replicates was observed, and was also expected due to ecological drift (Fig. 3b). Low SES values (Fig. 3a) suggest stochastic assembly among replicates of the same thickness, however, the relative importance of drift and dispersal cannot be disentangled with the experimental setup used here. In addition, due to the limited number of within-group replicates, these results should be interpreted with caution. Because MBBRs allow a high  www.nature.com/scientificreports www.nature.com/scientificreports/ level of replication in communities linked by dispersal, a similar setup to the one use here with higher replication could be used to study stochastic assembly and to confirm the possible existence of alternate states 15,47 . Overall, other studies have shown that stochastic and deterministic processes can co-occur in biofilms 14,16,17 . Our results suggest that the importance of deterministic vs. stochastic assembly depends on biofilm thickness: assembly would be deterministic between biofilms of different thickness, while assembly would likely be stochastic among biofilms with the same thickness.
Our hypothesis was that the communities in Z50 would be an aerobic subset of the ones in Z400. Thus, beta-diversity between Z50 and Z400 would largely be due to nestedness, whereas turnover would have a small contribution. This was expected due to different redox profiles between Z50 and Z400 biofilms (Fig. 6) which could create nestedness; oxygen in the thin Z50 biofilm inhibit the growth of obligate anaerobes like anammox bacteria 48 . Thus, richness in Z400 would be higher, because the community is a mixture of aerobic and anaerobic taxa. Between-group β sne was observed, suggesting nestedness, but the β ratio was below 0.5 (Fig. 3c), indicating that beta-diversity was dominated by turnover. Thus, the Z50 biofilm was not just a subset of the oxic upper layers of the Z400 biofilm. Instead, the differences between Z50 and Z400 communities were primarily due to turnover of taxa, which could be due to both ecological drift and deterministic processes. For example Nitrotoga was observed in Z50, but was nearly absent in Z400 (Fig. 4, Supplementary Fig. S1), which cannot be easily explained by redox profiles. This shows that thin biofilms favor the NOB Nitrotoga. Together with Nitrospira, Nitrotoga has recently been shown to be the predominant NOB in several WWTPs 49,50 . Changes in NOB composition may have consequences for operational strategies in wastewater treatment, for instance in systems with nitritation 51,52 , where NOB suppression is a prerequisite, since various NOB species respond differently to the suppression strategies 53,54 .
Redox profiles (Fig. 6) could explain the stratification of some taxa like anammox bacteria and Nitrospira in the Z400 biofilm (Fig. 5b). Nitrosomonas was the dominant population at the top of the Z400 biofilm (Fig. 5a,  Supplementary Fig. 3C) and was also abundant in Z50. However, Nitrosomonas aggregates were present throughout the Z400 biofilm, even in regions predicted to be anoxic (Fig. 5a,b), and thus might represent strains capable of anaerobic respiration 55,56 . Furthermore, in the thin Z50 biofilm, Nitrospira was seen alongside Nitrosomonas (Fig. 5c), and here its relative abundance was actually higher than in Z400. Hence, redox profiles alone cannot explain the distribution of taxa in the reactor. The fact that redox is not the only determinant of the distribution of microorganisms, even in strongly structured environments like sediments, has been noted 29 . The Z50 and Z400 biofilms also differed in their spatial structure, with Z50 being denser and having a smoother architecture, compared with the Z400 (Fig. 1) 36 . Furthermore, extracellular nucleic acids were observed in Z400 but not in Z50 (data not shown). Thus, these differences could contribute to the turnover of taxa between thin and thick biofilms observed in this study, by either selection or deterministic dispersal. Another possible mechanism for compositional turnover are biotic interactions. For instance, some SVs within the predatory Bdellovibrionales were differently distributed between the biofilms (Supplementary Fig. S2). It is plausible that the two biofilms represented different prey communities that in turn shaped the predatory Bdellovibrionales communities. Such influence on the predatory Bacteriovorax has been shown, even for closely related preys 57,58 . Furthermore, Torsvik et al. 59 suggested that predation can act as a major factor driving prokaryotic diversity. Hence, biological interactions, such as predation, could have had a large effect on these biofilm communities, as shown for other wastewater biofilms 60 .
Despite previous examples that have shown that species richness 46 may by itself lead to higher ecosystem function rates, differences in nitrogen transformation rates among Z50 and Z400 might not necessarily be linked to the observed differences in richness. Similar to earlier studies 38 , we found that the thinner biofilm had higher net NO 3 − production rates, despite having lower richness. This supports that species composition might be more important than richness for some processes 61 , such as nitrification. For example, increased evenness in the Z400 compared to Z50 biofilms could have resulted in lower relative abundance of specialized taxa 62,63 , such as Nitrosomonas and Nitrospira, and thereby decrease net NO 3 − production rates. Furthermore, the presence of anaerobic taxa in Z400 could have lowered net NO 3 − rates, via other processes like denitrification, anammox or DNRA (dissimilatory nitrate reduction to ammonia). This could also explain the observed higher anaerobic www.nature.com/scientificreports www.nature.com/scientificreports/ NO 2 − removal rates in Z400 than Z50; where the presence of anaerobic regions in Z400, allowed the establishment of taxa that could use NO 2 − as electron acceptor. This agrees with a previous study 38 , showing that an increase in biofilm thickness could lead to the emergence of new functions. Overall, this suggests that for nitrifying reactors, neither richness nor abundance of AOB are predictors for net NO 3 − production rates. Amplicon sequencing data as the one used in this study, is a measure of relative abundance and represents compositional data 64 . Furthermore the Z50 and Z400 biofilms differed in their richness and volumetric densities; therefore the absolute abundances are unknown. As difficulties arise when using abundance-based metrics, we used presence-absence metrics like the β sor and β sim in this study.
In summary, we show that biofilm thickness can influence bacterial biofilm community composition despite the fact that history and all other external conditions are similar. The differences in communities between thin and thick biofilms were likely deterministic, but differences could not always be easily explained just by differences in redox conditions (cf. 29 ). Between-group beta-diversity was primarily due to turnover, with nestedness having limited importance. Furthermore, based on potential and actual measurements, the two communities performed ecosystem functions at different rates, which support the idea that beta-diversity can lead to the emergence of multiple ecosystem functions 65 . Results from these and similar experiments can be used in design of new process strategies in wastewater treatment. For example, thinner nitrifying biofilms could be combined with ticker biofilms to increase the number of ecological functions 38 . Finally, multispecies bioreactors are well suited for experiments that can help disentangle factors of community assembly, as also suggested before 5,66 .
Methods the reactor. The 0.5 m 3 MBBR was located at the Sjölunda WWTP in Malmö, Sweden. The reactor was fed with effluent from a high-rate activated sludge process treating municipal wastewater (a feed with low carbon to nitrogen ratio). The average reactor load during one month before the sampling was 0.48 kg NH 4 + -N/m 3 ,day and the NH 4 + removal was 42%, at a pH of 7.4; dissolved oxygen (DO) concentration of 5 mg/L; and temperature of 17 °C. After 261 days of operation, carriers were sampled for DNA-sequencing, FISH and activity tests to determine nitrogen transformations. The reactor contained a mixture of Z50 and Z400 carriers (Veolia Water Technologies AB -AnoxKaldnes, Lund, Sweden) at a total filling degree of approximately 30%. Thickness of the biofilm in Z-carriers is limited by a pre-defined grid wall height 35 . Samples for optical coherence tomography measurements were taken on day 272 and data showed a biofilm thickness of 45 ± 17 and 379 ± 47 (mean ± S.D.) µm for Z50 and Z400, respectively 36 .
Nitrogen transformation activity tests. Actual activity was measured in 1 L reactors in duplicate: Two reactors with 100 Z50 carriers each, and two with 100 Z400 carriers each. The incoming water was the same as the water feeding the 0.5 m 3 reactor. At the time of measurement, the NH 4 + -N concentration was 19.6 mg/L, the DO was 5.5 mg/L, and the temperature was kept at 20 °C. Mixing was achieved by supplying a gas mix consisting of N 2 -gas and air to the bottom of the reactors at an approximate total flow of 3 L/min and the DO was controlled to 5.5 mg/L by adjusting the amount of air in the gas mix. Nitrification rates were measured from mass balance as NO 2 − N and NO 3 − N mg/m 2 ,day. For the potential activity trials 3 L reactors, containing 400 carriers each, were used. The substrate consisted of NaHCO 3 − buffer, pH adjusted to 7.5 using H 2 SO 4 , with phosphorous and trace minerals added in excess 35 . Aerobic removal of NH 4 + (starting concentration 35.2 NH 4 + -N mg/l) and NO 2 − (starting concentration 32.5 NO 2 − -N mg/l) were measured separately in two different trials at 20 °C for 1 hour, with sampling every 10 minutes. Mixing was achieved by supplying a gas mix consisting of N 2 -gas and air to the bottom of the reactors at an approximate total flow of 3 L/min. DO was controlled to 5.5 mg/L by adjusting the amount of air in the gas mix. Anaerobic trials of simultaneous removal of NH 4 + and NO 2 − (starting concentrations 35.5 NH 4 + -N and 36.1 NO 2 − N mg/l) was measured at 30 °C and were run for 2 hours with sampling every 20 minutes. Mixing was achieved by N 2 -gas from the reactor bottom. Before commencing the trials, the reactor with substrate was fed with N 2 -gas until the DO concentration was negligible and thereafter the carriers were added and the trials begun. Water samples were collected and filtered through 1.6 µm Munktell MG/A glass fiber filters and analyzed for NH 4 -N, NO 2 -N and NO 3 -N using standard Hach-Lange kits (LCK 303, 342 and 339, respectively).
Fluorescence in situ hybridization (FISH). FISH on cryosections and qFISH were done as previously described 36 . The FISH probes used in this study are shown in Supplementary Table S1. EPS and total nucleic acids on biofilm cryosections were stained with the FilmTracer SYPRO Ruby biofilm matrix stain and SYTO 40 (Thermo Fischer Scientific, USA), respectively. See supplementary information, for details.

Simulation of dissolved oxygen (DO) concentration profiles. A mathematical model was developed
for simulating DO concentration profiles in the biofilms. The model considered the activities of aerobic heterotrophic bacteria, AOB and NOB. The bulk liquid concentrations of substrates (DO, nitrite, ammonium, and organic compounds), the measured biofilm densities, the microbial community compositions (as determined by FISH), the distribution of different functional groups of microorganisms in the biofilm (as measured by FISH), and kinetic coefficients from the scientific literature were used as input parameters. The thickness of the liquid boundary layer that limits diffusion of soluble substrates, including DO, from the bulk liquid to the biofilm was determined by comparing the ammonium oxidation rates calculated by the model to those measured during the nitrogen transformation activity tests. Since the exact concentrations of active biomass in the biofilms were unknown, the model was solved for different scenarios in which the active biomass was assumed to make up 20-80% of the measured total dry solids. It should be noted that the model only considers biofilm heterogeneity in one dimension (the depth direction). Layers parallel to the substratum are assumed to be homogenous. Real www.nature.com/scientificreports www.nature.com/scientificreports/ biofilms are three-dimensional structures containing channels and voids, which may allow oxygen transport into deeper regions locally. See supplementary material for details. DNA extraction and 16S rRNA gene sequencing. DNA was separately extracted from ten Z50 and ten Z400 carriers. DNA extraction, PCR and high throughput amplicon sequencing of 16S rRNA gene was done as previously described 67 with some modifications. Sequence variants (SVs) were generated for finer resolution of taxa 68,69 . See supplementary information for details. Raw sequence reads were deposited at the NCBI Sequence Read Archive, no. SRP103666.
statistics. Data was analyzed in R (R Core Team 2018), using the packages Phyloseq 70 , Vegan 71 , DESeq2 72 and betapart 41 . Differential abundance of SVs was estimated with DESeq2 72,73 , without random subsampling before the analysis. After independent filtering in DESeq2, 2578 of 3690 SVs were analyzed. A p (adj) < 0.01 value (DESeq2) was used as criterion for statistical significance. Subsampling to even depth was done prior to estimation of alpha-diversity and beta-diversity. Alpha-diversity was calculated as the first two Hill numbers 74 , 0 D (richness) and 1 D (exponential of Shannon-Wiener index). Evenness was estimated as ( 1 D/ 0 D). Beta-diversity was estimated as pairwise Sørensen (β sor ) dissimilarities, a presence-absence metric. Principal coordinate analysis (PCoA) was used for ordination. Permutational multivariate analysis of variance (PERMANOVA) 75 was used test for significant difference between group centroids. The components of β sor , turnover (β sim ) and dissimilarity due to nestedness (β sne ), were estimated as described by Baselga et al. 41 and used to calculate the beta-diversity ratio (β ratio ) as the ratio between β sne and β sor 44 . If the β ratio is smaller than 0.5, beta-diversity is dominated by turnover rather than nestedness 44 .
To disentangle the contribution of stochastic and deterministic community assembly mechanisms while at the same time accounting for possible differences in richness between Z50 and Z400, a null model approach was used. Firstly, the standardized effect size (SES) for pairwise Sørensen (SES βsor ) dissimilarities were estimated in vegan using the oecosimu function. 999 null communities for estimation of SES βsor were generated using the quasiswap algorithm 76 , which preserve species richness and species incidence. For within groups null model analyses of Z50 and Z400 communities, only the taxa present in Z50 or Z400 respectively were used as the regional species pool. |SES βsor | > 2 was used as criteria to estimate if βsor was different than expected by chance; a |SES| > 2 value is approximately a 95% confidence interval 77 . Secondly, the RC bray metric 78 , which is based on quantitative data, was estimated for between-group comparisons, using 999 simulated communities. |RC bray | > 0.95 values were interpreted as deviations from the random expectation 43,78 .

Data Availability
All data generated or analyzed during this study will be available upon request to the corresponding author.