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.

Introduction

Biofilms are dense communities, encased in a polymer matrix, attached to a surface and/or each other[1] with a high microbial diversity compared to the bulk n class="Chemical">water system[1-3]. Microbial biofilms are important in aquatic ecosystems and are useful in many biotechnological applications, such as wastewater- or drinking water treatment. In nitrogen removal from wastewater, moving bed biofilm reactors (MBBRs) are used at many wastewater treatment plants (WWTPs). Here biofilms grow on so-called carriers, which move freely in the bioreactor (Fig. 1c), where each free-floating biofilm carrier represents a local community and the reactor represents the regional community. The local communities have defined boundaries and are separate, but are linked by dispersal with all other biofilm carriers in the reactor, in this case fed with wastewater from a full-scale WWTP to form nitrifying biofilms.

Figure 1

Biofilm structure shown by EPS staining of cryosections. The biofilm-water interface is the upper side. (a) Z400 biofilm. (b) Z50 biofilm. Scale bar: 100 µm. (c): Z400 (up) and Z50 (down) biofilm carriers; a ruler in cm is shown for size comparison.

Biofilm structure shown by EPS staining of cryosections. The biofilm-n class="Chemical">water interface is the upper side. (a) Z400 biofilm. (b) Z50 biofilm. Scale bar: 100 µm. (c): Z400 (up) and Z50 (down) biofilm carriers; a ruler in cm is shown for size comparison.

The mechanisms of community assembly are central in microbial ecology and in our understanding of formation of biodiversity in all ecosystems[4,5], including microbial communities in biofilms[6]. These mechanisms of community assembly include three major ecological processes: selection, drift and dispersal[7]. Selection, i.e. the sorting of species by prevailing local abiotic and biotic conditions, is deterministic, while drift results from stochastic birth and death events[5,7]. If local communities are further exposed to stochastic dispersal from the regional species pool (as we can assume to be the case of bioreactors in WWTPs), the expected result is that the abundance of a taxon in a local community can be predicted based on its respective abundance in the regional species pool and thereby follows neutral distribution patterns[8,9]. However, dispersal can also be deterministic if microorganisms differ in their ability to disperse within the complex spatial biofilm matrix or if their propagation is affected by interactions with species already present in the biofilm.

Selection has been suggested as the major mechanism for community assembly in stream biofilms[10-12], and other biofilms[13], while for biofilms within lakes linked by dispersal, both stochastic and deterministic factors were shown to be important[14]. The importance of both stochastic and deterministic factors was shown in an elegant study using parallel microbial electrolysis cells incubated with wastewater[15]. Other studies in wasten class="Chemical">water activated sludge systems have also shown the importance of deterministic and stochastic factors[16-19].

Deterministic assembly in biofilms could be due to specific mechanisms: Firstly, diffusion limitations form steep gradients of electron donors and acceptors in biofilms, which result in structured micro-environments. Examples are found in biofilms in WWTPs used for nitrification (i.e. the oxidation of n class="Chemical">ammonium to nitrate) during the nitrogen removal process[20,21]. Here population stratification typically occurs; ammonia oxidizing bacteria (AOB) are found closest to the oxygenated water and nitrite oxidizing bacteria (NOB) below the AOB[22-25]. If oxygen is consumed, anaerobic ammonium oxidizing (anammox) bacteria can establish in the deeper parts of the biofilm[22,26,27]. Similarly, in other multispecies biofilms anaerobic sulfate reducing bacteria are found in the biofilm interior[28]. However, functions in microbial communities are not always sorted according to such a thermodynamic “redox tower” of electron acceptors[29]. This makes detailed a priori predictions of community structure in biofilms difficult. Secondly, in addition to gradients, it was realized early on that microbial biofilms are in fact complex structures and not homogenous layers of randomly distributed organisms[3] and, ever since, architecture has been viewed as an important biofilm property. The intricate biofilm architecture consists of towers, mushroom-like structures and water filled channels[1,2,6,30]. If biofilms differ in their architecture, dispersal effects could influence community assembly by changing the available colonization surface. Furthermore, microorganism with deterministic dispersal might show preference towards different types of biofilms.

Thickness will likely influence several components of biofilm architecture, such as density, shape and porosity. It will also influence redox gradients and thus generally the local biofilm environment. However, the experimental evidence for effects of thickness on architecture, community structure and function has been difficult to obtain because biofilm thickness is the result of environmental conditions such as flow[1,6,31-33], nutrient conditions[23], development age of the biofilm[33], n class="Chemical">carbon to nitrogen (C/N) ratios[24] and temperature[34]. In most experimental systems, thickness cannot easily be isolated from these environmental factors that themselves can influence the community structure and functions.

Recently, a biofilm carrier with a defined grid wall height that defines the biofilm thickness was designed (Z-carriers, Veolia Water Technologies AB – AnoxKaldnes, Lund, Sweden)[35-37]. These carriers allow stringent experiments with different biofilm thicknesses, which have shown that thickness can affect some biofilm functions[37,38], evenness[38], biofilm architecture[36], abundance of key organisms[36,38] and functional stability after a disturbance[36]. Beside the opportunities to gain basic ecological knowledge by designing new experiments, the ability to control biofilm thickness opens for new process configurations in WWTPs. In this study, a pilot nitrifying bioreactor was filled with a mixture of Z-carriers with biofilm thicknesses of 50 and 400 µm. Thus, environmental conditions and history of the biofilms were the same. We ask if thickness, in itself or via diffusion effects or other mechanisms, is important for bacterial community structure, and if so, what the possible mechanisms of community assembly would be. The thicknesses we investigated are within the range commonly found in natural- as well as in n class="Species">man-made biofilms[24,25,31,39,40].

Differences between communities, i.e. beta-diversity could arise due to (a) turnover (species replacement) and (b) nestedness (species-loss), the latter resulting from differences in richness[41]. Accordingly, beta-diversity between thin and thick biofilms could result from selection of different species due to different environmental gradients[42], or random drift or dispersal processes. Secondly, differences in species richness could also lead to a high between-group beta-diversity[41-43]. More specifically, because thicker biofilms have larger volume and 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 wasten class="Chemical">water 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 n class="Chemical">Z400. CLSM images of EPS stained biofilm cryosections confirmed that carrier design limited biofilm thickness (Fig. 1a and b).

Alpha and beta-diversity

16S rRNA gene sequence variants (SVs) of ten Z400 and ten n class="Chemical">Z50 biofilms were analyzed. The alpha-diversity parameters richness (0D), first-order diversity (1D) and evenness (1D/0D) (Fig. 2a), were all significantly higher for the thick Z400 biofilms than for the thinner Z50 biofilms (Welch t-test, p < 0.05). We also estimated beta-diversity using the presence-absence based Sørensen index (βsor), which showed that Z50 and Z400 communities were different (PERMANOVA βsor, p = 0.002, r2 = 0.50) (Fig. 2b).

Figure 2

(a) Richness (0D), diversity (1D) and evenness (1D/0D) for the Z50 and Z400 biofilms. (b) PCoA based on the Sørensen index (βsor).

(a) Richness (0D), diversity (1D) and evenness (1D/0D) for the Z50 and n class="Chemical">Z400 biofilms. (b) PCoA based on the Sørensen index (βsor).

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 n class="Chemical">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 RCBRAY metric, also indicated that Z50 and Z400 communities were in average more dissimilar than the null expectation (between-group RCBRAY > +0.95).

Figure 3

(a) Standardized effect size for the Sørensen index (βsor); dashed lines indicate SES values of +2 and −2. (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.

(a) Standardized effect size for the Sørensen index (βsor); dashed lines indicate SES values of +2 and −2. (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), n class="Chemical">Z50 replicates (n = 10) and between the two groups.

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 is dominated by turnover rather than nestedness[44]. We observed a between-group βratio lower than 0.5 (Fig. 3c). Thus, beta-diversity between the n class="Chemical">Z50 and Z400 communities was caused by both nestedness and turnover, with turnover being more important. In addition, both βsim and βsne were significantly different between the Z50 and Z400 communities (PERMANOVA βsim, p = 0.001, r2 = 0.34; PERMANOVA βsne, p = 0.001, r2 = 0.83).

To investigate which taxa were contributing to the differences in community composition between Z50 and n class="Chemical">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 n class="Chemical">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].

Figure 4

(a) Relative read abundance of nitrifiers and anammox bacteria in Z50 and Z400. (b) Relative read abundance multiplied by total solids (TS) measurements for each carrier type. (c) Biovolume fractions of nitrifiers and anammox bacteria, as measured by qFISH.

(a) Relative read abundance of nitrifiers and anammox bacteria in Z50 and n class="Chemical">Z400. (b) Relative read abundance multiplied by total solids (TS) measurements for each carrier type. (c) Biovolume fractions of nitrifiers and anammox bacteria, as measured by qFISH.

Interestingly, we observed the anammox bacterium Brocadia in the Z400 biofilms, but it was almost absent in the n class="Chemical">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; n class="Species">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).

Figure 5

(a) FISH image of a Z400 biofilm cryosection; the water-biofilm interface is on the top. Green: Nitrosomonas. Red: Nitrospira. Yellow: Nitrotoga. Blue: Brocadia. Grey: SYTO. (b) FISH image of a Z50 biofilm cryosection; the water-biofilm interface is on the top. Green: Nitrosomonas. Red: Nitrospira. Yellow: Nitrotoga. Grey: SYTO. (c) FISH-based population distribution at different biofilm depths in Z400.

(a) FISH image of a Z400 biofilm cryosection; the n class="Chemical">water-biofilm interface is on the top. Green: Nitrosomonas. Red: Nitrospira. Yellow: Nitrotoga. Blue: Brocadia. Grey: SYTO. (b) FISH image of a Z50 biofilm cryosection; the water-biofilm interface is on the top. Green: Nitrosomonas. Red: Nitrospira. Yellow: Nitrotoga. Grey: SYTO. (c) FISH-based population distribution at different biofilm depths in Z400.

Dissolved oxygen (n class="Chemical">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.

Figure 6

DO concentrations profiles in the Z50 and Z400 biofilms. The shaded regions show ranges of DO concentration profiles resulting from different assumption about the fraction of the total dry solids on the carriers that is active bacteria. The dashed horizontal lines show the biofilm-liquid interface.

DO concentrations profiles in the n class="Chemical">Z50 and Z400 biofilms. The shaded regions show ranges of DO concentration profiles resulting from different assumption about the fraction of the total dry solids on the carriers that is active bacteria. The dashed horizontal lines show the biofilm-liquid interface.

Nitrogen transformation rates

Two types of tests were performed separately on the Z50 and n class="Chemical">Z400 carriers; (i) actual activities tested in a continuous laboratory trial, with the same incoming water as in the 0.5 m3 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 NO3− production were 1.4–1.5 gNO3−N/m2, d for Z50 and 0.68–0.72 gNO3−/m2, d for Z400. To estimate NO3− 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 NO3− production was higher in Z50 than in Z400.

In the aerobic potential tests for net NH4+ removal (Fig. 7a), net n class="Chemical">NO3− and net NO2− production (per carrier area) was higher for Z50 than Z400 biofilms (ANCOVA, p < 0.05), while the rate of net NH4+ removal was not significantly different between Z50 and Z400 (ANCOVA, p > 0.05). The aerobic potential removal of NO2− (Fig. 7b) was significantly higher for Z400 than for Z50 (ANCOVA, p < 0.05). Finally, in the anoxic potential trials, in which NH4+ and NO2− were added simultaneously (Fig. 7c), removal of NO2− was significantly higher for Z400 than for Z50 (ANCOVA, p < 0.05), while no significant removal or production of NH4+ was seen for either Z50 or Z400.

Figure 7

Potential conversion rates by carrier type during aerobic oxidation of NH4+ (a), aerobic oxidation of NO2− (b) and anoxic oxidation of NH4 (c) during batch tests. Significant differences between Z50 and Z400 (ANCOVA, p < 0.05) are shown with (*). Red: Z400, Blue: Z50.

Potential conversion rates by carrier type during aerobic oxidation of NH4+ (a), aerobic oxidation of n class="Chemical">NO2− (b) and anoxic oxidation of NH4 (c) during batch tests. Significant differences between Z50 and Z400 (ANCOVA, p < 0.05) are shown with (*). Red: Z400, Blue: Z50.

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 n class="Chemical">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 n class="Chemical">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 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 n class="Chemical">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 n class="Chemical">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 n class="Chemical">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 NO3− 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 NO3− production rates. Furthermore, the presence of anaerobic taxa in Z400 could have lowered net NO3− rates, via other processes like denitrification, anammox or DNRA (dissimilatory nitrate reduction to ammonia). This could also explain the observed higher anaerobic NO2− removal rates in Z400 than Z50; where the presence of anaerobic regions in Z400, allowed the establishment of taxa that could use NO2− 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 NO3− 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 n class="Chemical">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 wasten class="Chemical">water 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 m3 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 wasten class="Chemical">water (a feed with low carbon to nitrogen ratio). The average reactor load during one month before the sampling was 0.48 kg NH4+-N/m3,day and the NH4+ 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 n class="Chemical">Z400 carriers each. The incoming water was the same as the water feeding the 0.5 m3 reactor. At the time of measurement, the NH4+-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 N2-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 NO2− N and NO3− N mg/m2,day.

For the potential activity trials 3 L reactors, containing 400 carriers each, were used. The substrate consisted of NaHCO3− buffer, pH adjusted n class="Species">to 7.5 using H2SO4, with phosphorous and trace minerals added in excess[35]. Aerobic removal of NH4+ (starting concentration 35.2 NH4+-N mg/l) and NO2− (starting concentration 32.5 NO2−-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 N2-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 NH4+ and NO2− (starting concentrations 35.5 NH4+-N and 36.1 NO2−N mg/l) was measured at 30 °C and were run for 2 hours with sampling every 20 minutes. Mixing was achieved by N2-gas from the reactor bottom. Before commencing the trials, the reactor with substrate was fed with N2-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 NH4-N, NO2-N and NO3-N using standard Hach-Lange kits (LCK 303, 342 and 339, respectively).

Fluorescence 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. n class="Chemical">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 (n class="Chemical">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 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 n class="Chemical">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 rann class="Chemical">dom 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], 0D (richness) and 1D (exponential of Shannon-Wiener index). Evenness was estimated as (1D/0D). 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 n class="Chemical">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 RCbray metric[78], which is based on quantitative data, was estimated for between-group comparisons, using 999 simulated communities. |RCbray| > 0.95 values were interpreted as deviations from the random expectation[43,78].