Biofilm thickness matters: Deterministic assembly of different functions and communities in nitrifying biofilms

Microbial biofilms are important in natural ecosystems and in biotechnological applications. Biofilm architecture influences organisms’ spatial positions, who their neighbors are, and redox gradients, which in turn determine functions. We ask if and how biofilm thickness influences community composition, architecture and functions. But biofilm thickness cannot easily be isolated from external environmental factors. We designed a metacommunity system in a wastewater treatment plant, where either 50 or 400 µm thick nitrifying biofilms were grown simultaneously on biofilm carriers in the same reactor. Model simulations showed that the 50 µm biofilms could be fully oxygenated whereas the 400 µm biofilms contained anaerobic zones. The 50 and 400 µm biofilms developed significantly different communities. due to deterministic factors were stronger than homogenizing dispersal forces in the reactor, despite the fact that biofilms experienced the same history and external conditions. Relative abundance of aerobic nitrifiers was higher in the 50 µm biofilms, while anaerobic ammonium oxidizers were more abundant in the 400 µm biofilms. However, turnover was larger than the nestedness component of between-group beta-diversity, i.e. the 50 µm biofilm was not just a subset of the thicker 400 µm biofilm with reduced taxa richness. Furthermore, the communities had different nitrogen transformation rates. The study shows that biofilm thickness has a strong impact on community composition and ecosystem function, which has implications for biotechnological applications, and for our general understanding of biofilms. IMPORTANCE Microorganisms colonize all surfaces in water and form biofilms. Diffusion limitations form steep gradients of energy and nutrient sources from the water phase into the deeper biofilm parts, influencing community composition through the biofilm. Thickness of the biofilm will affect diffusion gradients, and is therefore presumably important for biofilm composition. Since environmental factors determine thickness, studies of how thickness influences biofilm functions and community assembly, have been difficult to perform. We studied biofilms for wastewater treatment with fixed thicknesses of 50 and 400 µm during otherwise similar conditions and history. Despite growing in the same wastewater reactor, 16S rRNA gene sequencing and confocal microscopy showed the formation of two different communities, performing different ecosystem functions. Using statistical methods, we show for the first time, how biofilm thickness influences community assembly. The results help our understanding of the ecology of microbial biofilms, and in designing engineered systems based on ecological principles.

INTRODUCTION taxon in a local community can be predicted based on its respective abundance in the regional 77 species pool and thereby follows neutral distribution patterns (9, 10). However, dispersal can 78 also be deterministic if microorganisms differ in their ability to disperse within the complex 79 spatial biofilm matrix or if their propagation is affected by interactions with species already 80 present in the biofilm. Finally, we also measured rates of nitrogen transformations in the two biofilms and investigated 145 whether possible differences in ammonia-and nitrite removal rates between biofilms of distinct 146 thicknesses are primarily linked to differences in taxa richness (i.e. nestedness) or identity (i.e. 147 turnover) and discuss the implications of the results for wastewater treatment. 148 149 RESULTS 150

Two different biofilms 151
We grew biofilm communities with a maximum thickness of 50 μm and 400 μm together in the 152 same bioreactor; these communities are referred to as Z50 and Z400. CLSM images of EPS 153 stained biofilm cryosections confirmed that carrier design limited biofilm thickness ( Fig. 1A and  154 B). 155 156 Alpha and beta diversity 157 The alpha diversity parameters richness ( 0 D), first-order diversity ( 1 D) and evenness ( 1 D/ 0 D) 158 ( Fig. 2A), were all significantly higher for the thick Z400 biofilms than for the thinner Z50 159 biofilms (Welch t-test, p<0.05). We also estimated beta diversity using the presence-absence 160 based Sørensen index (βsor), which showed that Z50 and Z400 communities were different 161 (PERMANOVA βsor, p=0.002, r 2 =0.50) (Fig. 2B). 162 163 We used null modelling to estimate the standardized effect size (SES) for βsor. We observed that 164 βsor values for between-group comparisons, i.e. between Z50 and Z400, were higher than 165 expected by chance βsor values for between-group comparisons (SESβsor > +2) (Fig. 3A) 166 indicating that between-group differences were likely deterministic. On the contrary, observed 167 βsor values for within-group comparisons, i.e. between the carriers of the same type, were not 168 more different than expected by chance (|SESβsor| < 2) (Fig. 3A). In addition, estimation of the 169 quantitative RCBRAY metric, also indicated that Z50 and Z400 communities were in average more 170 dissimilar than the null expectation (between-group RCBRAY > +0.95) 171

172
To determine whether the between-group beta diversity was due to nestedness or turnover, we 173 estimated the components of βsor; turnover (βsim) and dissimilarity due to nestedness (βsne) using 174 the Baselga framework (42) (Fig. 3B), to estimate the βratio (44). When βratio is smaller than 0.5, 175 beta diversity is dominated by turnover rather than βsne (44). We observed a between-group βratio 176 lower than 0.5 (Fig. 3C). The differences between the Z50 and Z400 communities due to 177 turnover were significant (PERMANOVA βsim, p=0.001, r 2 =0.34). Thus, the beta diversity 178 between the Z50 and Z400 communities was caused by both nestedness and turnover, with 179 turnover being more important. 180 181 Differences in relative abundance of taxa between Z50 and Z400 were estimated using DESeq2. 182 We found differential abundance (p(adj)<0.01, DESeq2) for 45% of the sequence variants (SVs) 183 analyzed with DESeq, while for the top 40 most abundant SVs, 32 had different abundance 184 between Z50 and Z400 (Fig. S1). Among the fraction with differential abundance, 26% of SVs 185 were more abundant in Z50, and 74% were more abundant in Z400. The effect of thickness on 186 relative abundance, if any, differed among taxa (for example see Fig. S1). 187 188 Between-group sorting of nitrifiers and anammox bacteria 189 The relative read abundance of the nitrifiers, Nitrosomonas, Nitrospira and Nitrotoga, was lower 190 in the Z400 biofilms with Nitrotoga being almost restricted to Z50 (Fig. 4A). The same trends 191 were noticed using quantitative fluorescence in situ hybridization (qFISH; Fig. 4C; Welch's t-192 test, p<0.05). It was not possible to detect by qPCR if comammox were present due to non-193 specific amplification using Nitrospira amoA primers (45). 194 195 Interestingly, we observed the anammox bacterium Brocadia in the Z400 biofilms, but it was 196 almost absent in the Z50 biofilms (Fig. 4A). This was supported by qFISH (Welch t-test 197 p<0.001) (Fig. 4C). Sorting of bacteria between thick and thin biofilms was not only limited to 198 primary producers (i.e. autotrophic nitrogen converters) but also seen among the predatory 199 Bdellovibrionales. Bacteriovorax had a higher abundance in the Z50 communities, while some 200 SVs classified as Bdellovibrio were more abundant in either Z400 or Z50 (Fig. S2). 201 202 FISH analyses of biofilm cryosections showed that the Z400 biofilm was stratified, e.g. with 203 Nitrospira being more abundant in the middle of the biofilm and the anaerobic anammox 204 bacteria being present in the deeper layers; a lack of stratification was observed Nitrosomonas 205 ( Fig. 5A and 5B). In the thin Z50 biofilms, no stratification was observed as the AOB and NOB 206 populations were located side by side (Fig. 5C). The calculated dissolved oxygen (DO) 207 concentration profiles in the biofilms are shown in Fig. 6; the results give a range of possible DO 208 concentration profiles, which are shown as shaded regions. The model predicts that Z50 biofilms 209 can be fully oxygenated but may also have anoxic regions, whereas the Z400 biofilms contain a 210 completely anoxic region in its deeper parts in all tested scenarios. 211 212

Nitrogen transformation rates 213
Two types of tests were performed separately on the Z50 and Z400 carriers; (i) actual activities 214 tested in a continuous laboratory trial, with the same incoming water as in the 0.5 m 3 reactor and 215 (ii) potential activities tested in batch trials where excess nitrogen was added. For all trials, 216 removal rates are reported per surface area and day. Actual rates of net NO3production were 217 1.4-1.5 gNO3 --N/m 2 , d for Z50 and 0.68-0.72 gNO3 --/m 2 , d for Z400. To estimate NO3 -218 production per nitrifier abundance, it is necessary to consider differences in biomass between 219 carriers. We estimate that the total nitrifier biomass per carrier surface was about the same in 220 Z50 and Z400 (Fig. 4B). Therefore, per nitrifier biomass, net NO3production was higher in Z50 221 than in Z400. 222

223
In the aerobic potential tests for net NH4 + removal (Fig. 7A), net NO3and net NO2production 224 (per carrier area) was higher for Z50 than Z400 biofilms (ANCOVA, p<0.05), while the rate of 225 net NH4 + removal was not significantly different between Z50 and Z400 (ANCOVA, p>0.05). 226 The aerobic potential removal of NO2 - (Fig. 7B) was significantly higher for Z400 than for Z50 227 (ANCOVA, p<0.05). Finally, in the anoxic potential trials, in which NH4 + and NO2were added 228 simultaneously (Fig. 7C), removal of NO2was significantly higher for Z400 than for Z50 229 (ANCOVA, p<0.05), while no significant removal or production of NH4 + was seen for either 230 Z50 or Z400. focusing on micro-pollutant degradation, also found a significant higher evenness in thicker 239 biofilms. 240 241 A null model approach was used to investigate if the differences in beta-diversity between Z50 242 and Z400 were due to deterministic or stochastic factors while accounting for the large 243 differences in richness between Z50 and Z400 (47). The results showed that the between-group 244 beta-diversity was higher than expected by chance (Fig. 3A), suggesting deterministic assembly 245 due to differences in biofilm thickness. This result was also confirmed by the fact that biofilm 246 thickness significantly affected the relative abundance of the majority of the most abundant 247 individual taxa, meaning that they showed clear preference for either thin or thick biofilms (Fig.  248   S1). Some turnover among the Z50 and Z400 replicates was observed, and was also expected 249 due to ecological drift. Low SES values (Fig. 3A) suggest stochastic assembly among replicates, 250 however, the relative importance of drift and dispersal cannot be disentangled with the 251 experimental setup used here. In addition, due to the limited number of within-group replicates, 252 these results should be interpreted with caution. Because MBBRs allow a high level of 253 replication in communities linked by dispersal, a similar setup to the one use here with higher 254 replication could be used to study stochastic assembly and to confirm the possible existence of 255 alternate states (16, 48). Overall, other studies have shown that stochastic and deterministic 256 processes can co-occur in biofilms (15, 17, 18). Our results suggest that the importance of 257 deterministic vs. stochastic assembly depends on biofilm thickness: assembly would be 258 deterministic between biofilms of different thickness, while assembly would likely be stochastic 259 among biofilms with the same thickness. 260 261 Our hypothesis was that the communities in Z50 would be an aerobic subset of the ones in Z400. 262 Thus, beta-diversity between Z50 and Z400 would largely be due to nestedness, whereas 263 turnover would have a small contribution. This was expected due to different redox profiles 264 between Z50 and Z400 biofilms (Fig. 6) which could create nestedness; oxygen in the thin Z50 265 biofilm inhibit the growth of obligate anaerobes like anammox bacteria (49). Thus, richness in 266 Z400 would be higher, because the community is a mixture of aerobic and anaerobic taxa. 267 Surprisingly, although between-group βsne was observed, the βratio was below 0.5 (Fig. 3C), 268 indicating that beta-diversity was dominated by turnover. Thus, the Z50 biofilm was not just a 269 subset of the oxic upper layers of the Z400 biofilm, but differences were primarily due to 270 turnover of taxa. For example, Nitrotoga was observed in Z50, but was nearly absent in Z400 271 (Fig, 4, S1), which cannot be easily explained by redox profiles. Independently of the 272 mechanism, it appears that thin biofilms favor the NOB Nitrotoga, which could have 273 consequences for operational strategies in wastewater treatment. 274 275 Redox profiles (Fig. 6) explain the stratification of some taxa like anammox bacteria and 276 Nitrospira in the Z400 biofilm (Fig. 5B). Nitrosomonas was the dominant population at the top 277 of the Z400 biofilm (Fig. 5A, S3C) and was also abundant in Z50. However, Nitrosomonas 278 aggregates were present throughout the Z400 biofilm, even in regions predicted to be anoxic 279 ( Fig. 5A, 5B). Furthermore, in the thin Z50 biofilm, Nitrospira was seen alongside Nitrosomonas 280 ( Fig. 5C), and here its relative abundance was actually higher than in Z400. Hence redox profiles 281 alone cannot explain the distribution of taxa in the reactor. The fact that redox is not the only 282 determinant of the distribution of microorganisms, even in strongly structured environments like 283 sediments, has been noted (30). The Z50 and Z400 biofilms also differed in their spatial 284 structure, with Z50 being more dense and having a smoother architecture, compared with the 285 Z400 ( Fig. 1) (37); furthermore extracellular nucleic acids were observed in Z400 but not in Z50 286 (data not shown). Thus, these differences could contribute to the deterministic turnover observed 287 in this study, by either selection or deterministic dispersal. Another possible mechanism for 288 between-group species turnover are biotic interactions. For instance, some SVs within the 289 predatory Bdellovibrionales were differently distributed between the biofilms (Fig. S2). It is 290 plausible that the two biofilms represented different prey communities that in turn shaped the 291 predatory Bdellovibrionales communities. Such influence on the predatory Bacteriovorax has 292 been shown, even for closely related preys (50, 51). Furthermore, Torsvik et al. suggested that 293 predation can act as a major factor driving prokaryotic diversity (52). Hence, biological 294 interactions, such as predation, could have had a large effect on these biofilm communities, as 295 shown for other wastewater biofilms (53). 296 297 Differences in community composition between Z50 and Z400 were to a larger extent 298 determined by turnover than nestedness. Therefore, we predict that differences in nitrogen 299 transformation rates among them might not necessarily be linked to the differences in richness 300 and evenness between Z50 and Z400. This is despite previous examples that have shown that 301 species richness (46) and evenness (54) may by themselves lead to higher productivity. Similar 302 to earlier studies (39), we found that the thinner biofilm had higher net NO3 --production rates, 303 despite having lower richness. This supports that species composition might be more important 304 than alpha-diversity for some processes (55), such as nitrification. Moreover, increased evenness 305 in the Z400 compared to Z50 biofilms could have resulted in lower abundance of specialized 306 taxa (56, 57), such as Nitrosomonas and Nitrospira, and thereby decrease net NO3 --production 307 rates. Despite differences in relative abundance in Nitrosomonas, their absolute abundance was 308 estimated to be the same in both Z50 and Z400. When measuring Nitrosomonas abundance, it is 309 assumed that the entire population might contribute to aerobic ammonia oxidation. However we 310 observed Nitrosomonas microcolonies throughout the Z400 biofilm depth (Fig 5A); the ones 311 living in the deeper parts of the biofilm might have little or no access to oxygen (58). These 312 Nitrosomonas cells could have low nitrification activity or they could represent strains capable of 313 anaerobic respiration (59). Therefore, Nitrosomonas abundance might not be directly correlated 314 with nitrification rates in thick biofilms. 315 316 A different ecosystem function, anaerobic NO2removal could occur via denitrification, 317 anammox or DNRA. We observed higher anaerobic NO2removal rates in Z400 than Z50. This 318 could be due to deterministic assembly, where the presence of anaerobic regions in Z400, likely 319 allowed the establishment of taxa that could use NO2as electron acceptor. This agrees with a 320 previous study (39) Between-group beta-diversity was due to both nestedness and turnover, but dominated by 330 turnover. Furthermore, based on potential and actual measurements, the two communities 331 performed ecosystem functions at different rates, which support the idea that beta-diversity in the 332 same metacommunity can lead to the emergence of multiple ecosystem functions (60). Results Since the exact concentrations of active biomass in the biofilms were unknown, the model was 395 solved for different scenarios in which the active biomass was assumed to make up 20-80% of 396 the measured total dry solids. It should be noted that the model only considers biofilm 397 heterogeneity in one dimension (the depth direction). Layers parallel to the substratum are 398 assumed to be homogenous. Real biofilms are three-dimensional structures containing channels 399 and voids, which may allow oxygen transport into deeper regions locally. See Text S1 for details. 400 401 DNA extraction and 16S rRNA gene sequencing 402 DNA was separately extracted from ten Z50 and ten Z400 carriers. DNA extraction, PCR and 403 high throughput amplicon sequencing of 16S rRNA gene was done as previously described (61)  Beta-diversity was estimated as pairwise Sørensen (βsor) dissimilarities, a presence-absence 417 metric. Principal coordinate analysis (PCoA) was used for ordination. Permutational multivariate 418 analysis of variance (PERMANOVA) (69) was used test for significant difference between group 419 centroids. The components of βsor, turnover (βsim) and dissimilarity due to nestedness (βsne), were 420 estimated as described by Baselga et al. (42) and used to calculate the beta diversity ratio (βratio) 421 as the ratio between βsne and βsor (44). If the βratio is smaller than 0.5, beta diversity is dominated 422 by turnover rather than nestedness (44). 423

424
To disentangle the contribution of stochastic and deterministic community assembly mechanisms 425 while at the same time accounting for possible differences in richness between Z50 and Z400, a 426 null model approach was used. Firstly, the standardized effect size (SES) for pairwise Sørensen 427 (SESβsor) dissimilarities were estimated in vegan using the oecosimu function. 999 null 428 communities for estimation of SESβsor were generated using the quasiswap algorithm (70), which 429 preserve species richness and species incidence. For within groups null model analyses of Z50 430 and Z400 communities, only the taxa present in Z50 or Z400 respectively were used as the 431 regional species pool. |SESβsor| > 2 was used as criteria to estimate if βsor was different than 432 expected by chance; a |SES| > 2 value is approximately a 95% confidence interval (71).    Assumed biomass distribution in the Z400 biofilm based on input from qFISH and cryosection 676 FISH images 677 Table S1. FISH probes used in this study 678