Distinct growth stages controlled by the interplay of deterministic and stochastic processes in functional anammox bioﬁlms

Mainstream anaerobic ammonium oxidation (anammox) represents one of the most promising energy- eﬃcient mechanisms of ﬁxed nitrogen elimination from wastewaters. However, little is known about the exact processes and drivers of microbial community assembly within the complex microbial bioﬁlms that support anammox in engineered ecosystems. Here, we followed anammox bioﬁlm development on fresh carriers in an established 8m 3 mainstream anammox reactor that is exposed to seasonal temperature changes (~25-12 ° C) and varying NH 4 + concentrations (5-25 mg/L). We use ﬂuorescence in situ hybridization and 16S rRNA gene sequencing to show that three distinct stages of bioﬁlm development emerge naturally from microbial community composition and bioﬁlm structure. Neutral modelling and network analysis are employed to elucidate the relative importance of stochastic versus deterministic processes and synergistic and antagonistic interactions in the bioﬁlms during their development. We ﬁnd that the different phases are characterized by a dynamic succession and an interplay of both stochastic and deter- ministic processes. The observed growth stages ( Colonization, Succession and Maturation ) appear to be the prerequisite for the anticipated growth of anammox bacteria and for reaching a bioﬁlm community struc- ture that supports the desired metabolic and functional capacities observed for bioﬁlm carriers already present in the system (~100g NH4-N m 3 d − 1 ). We discuss the relevance of this improved understanding of anammox-community ecology and bioﬁlm development in the context of its practical application in the start-up, conﬁguration, and optimization of anammox bioﬁlm reactors.


Introduction
Anaerobic ammonium oxidation (anammox) involves the simultaneous oxidation of ammonium (NH 4 + ) and reduction of nitrite (NO 2 − ) under oxygen-limiting conditions, and is orchestrated by a unique lineage of bacteria (AMX), which all belong to the phylum Planctomycetes, and display a high diversity with 10 different candidate species discovered so far ( Kartal et al., 2011 ;Strous et al., 2006 ;Zhang and Okabe, 2020a ). Used under mainstream conditions in wastewater treatment plants (WWTP), anammox is also one of the most promising mechanisms of fixed nitrogen (N) elimination, since it represents an important step towards energy autarky in the treatment of municipal wastewaters ( Lotti et al., 2014 ;Siegrist, Salzgeber, Eugster and Joss, 2008 ). Currently, anammox is widely applied and represents a robust method for autotrophic N removal under sidestream conditions (high NH 4 + concentrations and temperatures) ( Lackner et al., 2014 ;Nsenga Kumwimba et al., 2020 ). However, application of the anammox processes under mainstream conditions suffers from process instabilities due to unexpected fluctuation of environmental temperature, dissolved oxygen, low NH 4 + concentrations and high C/N ratios ( Joss et al., 2011 ;Laureni et al., 2016 ;Pijuan et al., 2020 ;Wells et al., 2017 ).
AMX currently lack representatives available in pure culture. They are characterized by very slow growth rates and low cell  ( Kuenen, 2008 ), which makes sufficient retention of biomass one of the main challenges. For their application in wastewater treatment, they are grown either in biofilm reactors on various carrier materials or as granules to retain sufficient biomass ( Abma et al., 2007 ). Lower space requirements, lower sludge production, and resilience towards changes in the reactor configuration are additional benefits of biofilm-based wastewater treatment systems ( Chen and Chen, 20 0 0 ;Wilderer and McSwain, 2004 ;Zhao et al., 2019 ).
The microbial communities that ultimately form the biofilms performing anammox in wastewater treatment systems (AMX biofilms) are very diverse. AMX, although performing the anammox reaction, typically comprise only a small fraction of the community. Heterotrophic bacteria, that are crucial for the formation of the biofilm but do not contribute directly to the anammox process itself, typically represent the majority of microbes in these biofilms ( Laureni et al., 2015;Mozumder et al., 2014;Ni et al., 2012 ). Additionally, many heterotrophic bacteria in mainstream anammox systems have the functional capacity to perform denitrification (complete or single steps) and to contribute to the N removal in the reactor or scavenging of detritus and peptides produced by anammox bacteria Niederdorfer et al., 2021 ;Speth et al., 2016 ).
While there are numerous studies focusing on microbial community assembly in systems with suspended biomass ( Chu et al., 2015 ;Gonzalez-Martinez et al., 2015 ;Luo et al., 2017 ), the ecological processes underlying AMX biofilm formation on synthetic carrier materials are still poorly understood. The high diversity and internal structure indicates a high potential for internal competition, synergies and feedback mechanisms. At the same time, the system is in constant contact with the bacteria and variable environmental conditions of the inflow. Therefore both deterministic Flemming and Wuertz, 2019 ) as well as stochastic effects ( Sloan et al., 2006 ;Woodcock and Sloan, 2017 ) are expected to govern community assembly. This complicates efforts to fully understand AMX biofilm formation for engineered applications ( Lawson et al., 2019 ). Deterministic processes involve non-random, niche-based mechanisms including environmental filtering and various biological interactions between community members (e.g., competition, facilitation, mutualism, and predation) ( Chase and Myers, 2011 ;Ofi¸t eru et al., 2010 ;Vellend et al., 2014 ). Stochastic processes, (termed neutral effects in community ecology) on the other hand are random ecological processes that generate community diversity patterns. These processes typically include probabilistic dispersal, random speciation and extinction, and ecological drift (e.g., random changes in organism abundance) ( Chase and Myers, 2011 ;Zhou et al., 2013 ;Zhou and Ning, 2017 ). A neutral model posits that random variation in extinction and speciation events, coupled with limited dispersal, can account for relative abundance distribution ( Chisholm, Burgman, Hubbell and Borda-De-Água, 2004 )). Successful microbial community engineering for wastewater treatment will require better mechanistic understanding of these processes ( Kleerebezem and van Loosdrecht, 2007;Moralejo-Gárate et al., 2011;Vlaeminck et al., 2012 ).
A spatial, i.e. microscopy-based investigation of AMX bacterial biofilm provides important information about the structural and spatial requirements of AMX within such a complex community ( Almstrand et al., 2013 ;Kindaichi et al., 2007 ;Laureni et al., 2015 ). Understanding of the spatial structure as well as temporal community dynamics, is required to understand the ecological processes that lead to an established biofilm community ( Brislawn et al., 2019;Datta et al., 2016;Niederdorfer et al., 2016;Shade et al., 2014 ). Given the high microbial diversity of mainstream AMX biofilms, co-occurrence network modelling provides a useful tool to elucidate the mechanistic background of ecological interactions between individual community members ( Faust and Raes, 2012 ;Ju and Zhang, 2015 ;Röttjers and Faust, 2018 ;Widder et al., 2016 ).
In this study, we investigated the microbial community assembly on synthetic biofilm carriers in a pilot scale mainstream anammox reactor in the presence of already established mature biofilm. We studied the spatial biofilm structure using confocal microscopy and fluorescence in situ hybridization. The microbial community composition was analyzed with 16S rRNA gene amplicon sequencing. The establishment of the biofilm was monitored over the course of a year from initial colonization to mature biofilms.
We hypothesized i) that heterotrophic bacteria are the pioneering colonizers that initiate biofilm formation and create the required base layer for second-wave colonizers including AMX, and ii) that the structure of, and development within, biofilms can be divided into three distinct growth stages -colonization, succession and maturation, as observed in other systems ( Allison et al., 20 0 0;Hall-Stoodley et al., 2004 ): We expected colonization to be driven by initial colonizers capable of attaching and growing on the carrier surface. During the succession stage we expected increasingly intense competition for space and nutrients in the juvenile biofilm. During maturation, we expected convergence towards a well-adapted, quasi-equilibrium state determined by the environmental factors in the reactor. We further hypothesized iii) that the relative importance of stochastic (e.g., growth, death and immigration) versus deterministic processes (e.g., environmental factors, biotic interactions, niche differentiation, legacy effects) for community assembly changes over the different stages of biofilm development.
Using neutral, co-occurrence and checkerboard score modeling of the AMX biofilm community assembly allowed us to disentangle the deterministic and stochastic ecological factors that ultimately shape the community structure of the mature biofilm. To our knowledge, this is one of the first studies using a toolbox of various community modelling approaches combined with spatial analysis to unravel assembly mechanisms of surface-attached AMX biofilms in engineered ecosystems. For practitioners and engineers, this work provides useful information on ecological constraints governing the establishment of functional AMX biofilms under mainstream conditions, but also highlights their inherent stability.

Experimental Set-Up and Sampling
The experiment was conducted in an 8 m 3 mainstream anammox moving bed biofilm reactor, which is part of the pilot wastewater treatment plant at the Swiss Federal Institute of Aquatic Science and Technology (Eawag) in Dübendorf. A separate, SBR operated, A-stage pilot-scale reactor (8m 3 ), fed with fresh mainstream municipal wastewater from the municipality Dübendorf, provided the inflow for the anammox reactor. The influent has the following composition: 50-60 mg ·L −1 COD, 5-25 mgNH 4 + -N ·L −1 , 0-1 mgNO 2 − -N ·L −1 and 5-10 mgNO 3 − -N ·L −1 . Neither A-Stage nor anammox stage has sludge recirculation. Given seasonal temperature variations (~25-12 °C) of the inflowing wastewater and the varying NH 4 + concentrations of 5-25 mg/L our anammox reactor was exposed to mainstream conditions throughout its operational history. The reactor was operated in sequencing batch reactor (SBR) mode and was discontinuously filled and emptied every 6-7 hours, depending on the process cycle length. The reactor was operated as follows: (1) Decantation: The reactor was decanted from 8 m3 to 4.5 m3 (ca. 60 min), while the carriers were retained with a sieve.
(2) Filling: The reactor was filled with effluent from the Astage SBR reactor from 4.5 m3 to 8 m3 (60 min), i.e., 44% volume exchange. (3) Reaction phase: During the reaction phase, the reac-tor was stirred and nitrite was continuously added (3.5 M NaNO2 solution) until an ammonium concentration of < 1 mgNH 4 + -N ·L −1 was measured (120-180 min). As the reactor system is intended to simulate the anoxic AMX stage of a 2-stage process, we did not apply aeration strategies for aerobic ammonium oxidizing bacterial growth and always maintained an anaerobic environment. To prevent washout of the microbial biomass during wastewater decanting, 150 m 2 /m 3 of FLUOPUR® carrier material (WABAG Water Technology Ltd., Switzerland) fleece tiles made of synthetic polymeric fibers (12 ×12 ×0.1 mm) were used in the anammox reactor. At the start of our experiment, the AMX reactor had been operating successfully for over 2 years (~100g NH4-N m 3 d −1 ). For the colonization experiment, we added 35,0 0 0 marked (cut corners) blank biofilm carriers to the system. Carrier samples (n = 6) were collected twice per week over the first eight weeks after carrier addition, and on a weekly basis thereafter, for a total period of one year. Non-marked carriers that had stayed in the tank since June 2016 were also sampled at each time point for comparison. The samples were immediately transferred to cryo-tubes, snap frozen in liquid nitrogen and stored at -80 °C until further analysis.

Fluorescence in situ hybridization (FISH)
Six sampling points (one, three, six, nine, twelve and 24 + months) were chosen for cryosectioning and analysis by FISH. Duplicate carriers were analyzed for each time point, with an additional blank carrier as a control. The carriers were cut into three sections, with the middle section being used for fixation, cryosectioning and FISH-Confocal laser scanning microscopy (CLSM). FISH was performed according to established standard operational procedures ( Nielsen, 2009 ) using the oligonucleotide probes listed in Supplementary Table 1. The probe EUB 338, specific for the domain Bacteria, was used together with a number of species-specific 16S rRNA-directed oligonucleotide probes to identify, nitrite-oxidizing bacteria (NOB), and AMX in the cryosectioned samples ( Laureni et al., 2015 ). We use the term EUB for bacteria stained by EUB 338 and not by the specific probes. CLSM of the FISH-slides was performed using a Leica SP5 DMI 60 0 0 (Leica Microsystems GmbH, Germany) and the software LAS AF v2.7.9. Images with a resolution of 1024 ×1024 pixels (4.096 pixels/μm) were collected using a 63x/1.40 oil objective. The image files were visualized and extracted in Fiji ( Schindelin et al., 2012 ). Daime (version 2.1) ( Daims, Lücker and Wagner, 2016 ) was used for 3D visualization and digital image analysis of the z-stacks. We used 30 fields of view and 15z stacks for each sample. More detailed information on FISH, CLSM and image analysis can be found in the Supplementary Information.

Biomass sampling, extraction and sequencing
Nucleic acids were extracted based on a method modified from Griffiths et al. ( Griffiths, Whiteley, O'Donnell and Bailey, 20 0 0 ). Duplicate biofilm carriers (n = 3) from every time point were cut into small pieces on a liquid nitrogen bath. Carrier pieces were transferred to 1.5 ml Matrix E lysis tubes (MPbio), and 0.5 ml of both hexadecyltrimethylammonium bromide buffer and phenol:chloroform:isoamylalcohol (25:24:1, pH 6.8) was added. Cells were lysed in a FastPrep machine (MPbio), followed by nucleic acid precipitation with PEG 60 0 0 on ice. Nucleic acid pellets were washed three times with ethanol (70%), and dissolved in 50 μl DEPC treated RNAse free water. DNA quality and quantity was assessed by using agarose gel electrophoresis and a Nanodrop ND-20 0 0c spectrophotometer (Thermo Fisher Scientific, USA). 16S rRNA gene amplicon sequencing was performed by Novogene Ltd. (Hong Kong) on the Illumina NextSeq platform, based on a pairend algorithm (250bp, V3-V4).

Sequence analysis and ecostatistics
Raw sequences were analyzed within the QIIME2 framework ( Bolyen et al., 2019 ). Primer sequences were removed with the cutadapt QIIME2 plugin. After demultiplexing, read pairs were joined, low-quality reads were filtered out and all high quality reads were analyzed with the DEBLUR software to produce amplicon sequence variants (ASVs) based on Illumina Miseq/Hiseq error profiles. Taxonomical assignment of the ASVs was performed within QIIME2 environment with the SILVA 16S V3/V4 classifier (accessed January 2020). After filtering of unclassified and contaminated ASVs, the resulting sequence table consisted of 1038 ASVs. All subsequent analyses were performed on the relative abundances sequence table (rarefied to 15595 reads) to ensure comparability between the samples. Based on the relative abundances, we plotted the top bacterial phyla over time for each sample. Non-metric multidimensional scaling (nMDS) analysis on Bray-Curtis dissimilarity matrices served to visualize patterns of community composition, and PERMANOVA tested for differences among communities. Based on the results of nMDS, and visual observations, we categorized the dataset into colonization (Col), succession (Suc), maturation (Mat). Old carriers were treated as a separate category (Old).
We calculated the Bray-Curtis dissimilarity between adjacent time points and calculated the fraction of dissimilarity contributed by ASVs identified as residents (all ASVs that were present at least in 80% of all samples over all development phases), transients (ASVs that were present in at least 80% of the samples of any two phases and absent during the other), exclusives (ASVs that occurred in > 80% of samples of only one biofilm development phase and were absent during the others) and the ones that could not be categorized (other). To clarify the distribution we divided the transients into early transients (only occurring during Col and Suc) and late transients (only occurring during between Suc and Mat).

Neutral modeling of the stochasticity of the community assembly process
To assess the relative importance of stochasticity and determinism for anammox community assembly, we firstly evaluated the fit of a Sloan neutral community model ( Sloan et al., 2006 ) for each stage. Output plots mainly show: (i) the fit of neutral model, r2, and the migration rate, m; (ii) neutral model predictions and its corresponding 95% confidence intervals (lines) and actual distributions (points) for each successional stage. We then calculated an index-normalized stochasticity ratio (NST) for each stage to verify the results of the Sloan neutral model. Analysis was conducted in R using the package NST ( Ning, Deng, Tiedje and Zhou, 2019 )

Co-occurrence network and C-score modeling of process determinism
To construct co-occurrence networks for microbial community in the anammox biofilm system, Spearman's rank coefficients ( ρ) between all ASVs abundances with occurrence in more than 50% of samples were calculated. Correlations with ρ> 0.8 and Pvalue < 0.05 were considered statistically significant and robust. A series of topological properties including modularity, average clustering coefficient, average shortest path length and network density were calculated and compared across each stage.
The random (R%) and observed (O%) incidences of intra-group and inter-group co-occurrence patterns between microbial entities (i.e., ASVs) were calculated ( Ju and Zhang 2015 ) and a checkerboard score (C-score) model ( https://github.com/GotelliLab/ EcoSimR ) was used as the metric of species segregation ( Stone and Roberts, 1990 ), and its variance (C-scorevar) was used as the metric of both species segregation and aggregation in the microbial community. The higher the observed standardized effect size (SES) values of the C-score, the greater the degree of species segregation, relative to a random distribution. Python and R scripts used for co-occurrence network and C-score analysis are freely available at https://github.com/RichieJu520 . More information on the modeling approach can be found in the SI.

Temporal development of anammox biofilm structure revealed by FISH-CLSM
During the first month of biofilm growth, FISH-CLSM showed small patches of pioneering EUB attached to the carrier material in the void space between the fleece fibers ( Fig. 1 ). While AMX did not contribute to the earliest phase of biofilm development, NOB were already part of the community ( Fig. 1 ). After three months, AMX also started to appear as microcolonies, while existing EUB patches increased in size. After six months, larger and sporadically connected agglomerations of EUB had developed. Numerous small AMX and NOB colonies are attached to, or enclosed within, the EUB patches. After nine months of biofilm growth, AMX clusters displayed a significant increase in volume. EUB patches had grown and combined with each other to form larger agglomerations, while the NOB fraction continued to form additional small colonies. Only after twelve months, NOB patches started to occasionally connect to form larger agglomerations, but never to the extent seen with AMX and EUB ( Fig. 1 ). Finally, in the matured biofilms (Month 24 + ), both AMX and EUB agglomerations had spread even further, now occupying more than 50% of the total carrier volume.
Comparison of bright-field and FISH microscopy revealed that bacteria only grew around, and not within, the fibers of the fleece carrier ( Fig. 1 , Month 24 + ). The available pore space within the carrier appeared to be nearly completely filled with bacterial biomass in the mature and old biofilms.
We used image analysis to determine the average total volume and average number of detected objects from the image data for each time point ( Table 1 ). The average object volume (i.e. average volume of colonies or patches), expressed as the ratio between av- Table 1 Mean volume (per microscopic scan (512 ×512 pixel)) and mean number of objects (per microscopic scan (512 ×512 pixel)) derived from the FISH-CLSM analysis with DAIME. Standard deviation is based on the variation of 30 randomly analyzed pictures for every timepoint.  erage total volume and average number of detected objects, confirmed the qualitative visual observations ( Fig. 2 A). EUB initiated the biofilm formation, but then expanded rather slowly in volume ( Fig. 2 B). Between three and six months, the average object volume increased for all three bacterial consortia. The NOB average object volume remained far lower than for EUB and AMX, and was relatively stable over time, although the overall volume and number of objects increased until nine and twelve months after initiation of the experiment, respectively ( Table 1 ). This suggests a lack of expansion of the existing NOB clusters in parallel with the development of new clusters ( Fig. 2 B). The average object volume for AMX and EUB on the other hand started to increase significantly after six months (t-test, p < 0.005), indicating expanding volume by coalescence of objects ( Fig. 2 B). Our findings for old biofilms also demonstrate that in the additional year (or more) of development the EUB and AMX clusters nearly doubled in size and number of colonies, while NOB clusters remained very stable. All three bacterial groups showed a significant correlation between the mean total volume and the mean number of objects (linear regression; R 2 EUB = 0.960, R 2 NOB = 0.946, R 2 AMX = 0.783). Based on the above observations, we categorized biofilm growth into three distinct biofilm growth phases: The Colonization (Col) phase lasts until three months after carrier addition, when AMX appeared for the first time, but the EUB volume/object ratio is still low. The Succession (Suc) phase includes the dynamic expansion of the EUB and AMX clusters until Month 9. The final Maturation (Mat) phase comprises the remainder of the study period, when the volume/object ratio remains mostly constant, although total biofilm volume still increases.

Temporal dynamics of the biofilm microbial community composition
The development of the biodiversity and composition of the bacterial community was investigated by 16S rRNA gene amplicon sequencing. In total, 1038 ASVs (proxy microbial taxa) were detected in the 102 biofilm samples. After assigning ASVs to the microbial taxonomy, we noted that many of the identified taxa are typical for mainstream anammox in WWTP Lawson et al., 2017 ;Speth et al., 2016 ). Throughout the observation period, the biofilm was dominated by ASVs belonging to the phyla Proteobacteria, Bacteroidetes, Chloroflexi, Actinobacteria, Acidobacteria, Firmicutes, Planctomycetes (i.e., the phylum to which AMX belong) and Nitrospirae (performing nitrite oxidation) ( Fig. 3 A). The microbial community was subject to dynamic changes in the phylum-abundance structure and turnover of representatives of these phyla throughout the one-year experiment. However, the distinct growth phases defined above were also characterized by differences in the dominant microbial taxa and their temporal dynamics. During the Col phase of biofilm development, ASVs that were assigned to the phyla of Proteobacteria (mainly Betaproteobacteriales ) , Bacteroidetes (Chitinophagales and Flavobacteriales), Actinobacteria (Micrococcales) and Firmicutes (Lactobacillales) displayed highest relative abundances. By assessing the we found that all of them are pursuing a heterotrophic lifestyle. Furthermore, members of these bacterial groups were also associated with the formation of suspended aggregates ( Chu et al., 2015 ;Gonzalez-Martinez et al., 2015 ).
Members of the Firmicutes nearly vanished towards the end of the Col phase and were not able to reestablish their prior dominance in any of the later stages. In contrast, ASVs from the phylum Chloroflexi (from the bacteria classes Anaerolineae and Chloroflexia ) increased significantly (Students t-test, p < 0.005) in their relative abundance, reaching up to 20% in the Mat phase, and even higher percentages in the Old samples ( Fig. 3 A). Relative abundance of the Proteobacteria decreased during the Col phase, but they remained the most abundant phylum in all later phases (~40%). Actinobacteria and Bacteroidetes increased during the Col phase but decreased during succession to more or less stable abundance levels. This was also observed for members of the phyla Acidobacteria ( Fig. 3 A). Consistent with the FISH observations, AMX ASVs ( Planctomycetes ) and Nitrospirae displayed very low abundances during the first three months of biofilm growth, but started to increase in numbers during the Suc phase, reaching up to ~8% of the total community. After the Suc phase (nine months) the community composition stabilized with less evident fluctuations in community composition and ASV abundances ( Fig. 3 A). The presence of NOB, despite the lack of oxygen, was very surprising but the overall abundance remained small and we can exclude that substrate competition for NO 2 between AMX and NOB affected process performance. Furthermore, as observed in simultaneous nitritation, anammox and denitrification reactors ( Du et al., 2021 ;Wen et al., 2017 ) heterotrophic denitrification and anammox can occur in parallel. However, denitrification can be excluded as a major competition for the available NO 3 + or NO 2 − in our system as COD concentrations are small and transcription rates of relevant denitrification genes were shown to be very low in comparison to AMX gene transcription ( Niederdorfer et al., 2021 ).
To evaluate how the microbial community structure changed over time, and to verify whether our categorization into different growth phases is supported by community dissimilarity metrics, we performed an ordination analysis on the ASV abundance data (nMDS; Fig. 3 B). We found that community composition of the three growth phases and the Old carriers, were all significantly different (Permanova, p < 0.001, R 2 : 0.5783). Furthermore, community dissimilarities within growth phases were much more pronounced during the Col phase and gradually decreased as biofilms matured. For example, the late stages of Suc and early stages of Mat phase (Week 30-35) did not display any significant differences (Permanova, p > 0.5). As one may expect, the overall community composition on the added carriers appeared to converge towards the Old carriers, but remained significantly different from these (Permanova, p < 0.001). However, even within the Old biofilms we still found significant differences in community composition over time (Permanova, p < 0.001), but without any clear temporal trend.

Successional patterns of individual microbial taxa
A distribution analysis of ASVs occurrence over the distinct growth phases revealed a large number of ASVs that were abundant only during Col phase. Numerous ASVs were only shared between the Suc and the Mat phase and were low in abundance during the Col phase (Supplementary Figure 1). Therefore, we clas-sified ASVs into residents (always present), transients (present in any two phases) and exclusives (present in one phase) depending on their occurrence (or detection) frequencies over the growth phases (see Methods). We found 373 ASVs that were classified as residents. Col and the Old phase displayed the highest number of exclusive ASVs (n = 137 and 52, respectively), Suc and Mat shared the highest number of transient ASVs (n = 211, late transients). Col and Suc shared 41 transient ASVs (early transients), while Col and Mat only shared 17 (re-occurring transients).
To assess the contributions of the defined successional categories to the temporal community turnover, we computed the fraction of the Bray-Curtis dissimilarity explained by each group ( Fig. 4 ). Overall, ASVs belonging to the resident fraction contributed most to the overall community dissimilarity (64 ± 5 %), indicating that they were highly dynamic even though they remained in the biofilm throughout the experimental period. The contribution of the exclusives to the overall community change was highest during Col (14.7 ± 4.5 %), and, unexpectedly, elevated also in old biofilms (11.6 ± 2.9 %) compared to Suc and Mat phases (4.1 ± 1.3 and 4.6 ± 3.6%, respectively). The average contributions of exclusives were significantly different in pairwise comparisons of all phases (Students t-test, p < 0.005). The contribution of transient ASVs increased from 22 ± 3.3 % during Col, to 28 ± 3.1 % in the succession phase, and to 29.7 ± 5% in the Mat phase, respectively. During the Old phase, the contribution of transients to the overall community change decreased to 20.9 ± 3%. ( Fig. 4 ).

Stochasticity versus determinism in community assembly during biofilm development
Given the above observations the dynamics of community assembly clearly changed over the growth phases. We therefore applied a neutral modeling approach to determine the importance of stochasticity in the community assembly. Based on the differences in the neutral model predictions (fitness of the model (r 2 ), immigration rate (m)) and the index-normalized stochasticity ratio NST, this approach led us to further divide the Col phase into three distinct subsections (Col1, Col2, and Col3). We observed a general decrease of r 2 and NST (from 0.75 to 0.59) over the first three months of Col ( Fig. 5 A, B, C; Table 2 ), indicating that the impor-

Deterministic co-occurrence and species segregation patterns
Assuming that both synergistic and antagonistic species interaction play a role in community assemblage, we further explored co-occurrence patterns using network analysis. We built networks from strong and significant correlations of the Spearman's rank coefficient. The networks for each of the growth stages differed substantially in the number of nodes (ASVs), edges (number of connections), size of modules (closely interconnected nodes), modularity (extent to which species interaction are organized into modules or subnetworks in a network), and the clustering coefficient (the degree to which nodes tend to cluster together) ( Fig. 6 ; Table 3 ). For all parameters except modularity, we observed a decrease from Col1 to Col2, whereas all network properties significantly increased again from Col2 to Col3 ( Fig. 6 A, B, C; Table 3 ). The Suc phase displayed lowest values for all network properties with an all-time Fig. 6. Co-occurrence model of AMX biofilm community assembly. A-F Co-occurrence networks colored by modularity class. Each node represents a species with its size proportional to its relative abundance in all samples, and each edge stands for a very strong (Spearman's ρ ≥ 0.8) and significant (FDR-adjusted P-value < 0.05) Spearman's correlation. For more information see Table 3 and 'Material and Methods'. low in nodes (n = 94) and edges (n = 87), but an increase in average modularity from 0.295 (Col3) to 0.795 ( Fig. 6 D). During Mat phase, the number of nodes and edges increased again (n = 213 and 206 respectively). In addition, an increase in modularity was observed ( Fig. 6 F; Table 3 ). Networks for Old carriers that have been in the reactor for over two years displayed a substantial increase in nodes (n = 363) and edges (N = 1027), but a decrease in modularity (0.687) in comparison to the Mat biofilm.
We analyzed the networks for co-occurrence of ASV within and between individual Phyla. We found that in Col1 ASVs tended to not co-occur deterministically, as indicated by the low O/R ratios (Supplementary Table 2). Only ASVs assigned to the phylum of Bacteroidetes appeared to co-occur with other representatives of the same phylum (intra-taxon co-occurrence) during this phase. This changed drastically in Col2 and Suc phases, when the number of significant positive inter-taxa and intra-taxon co-occurrences increased substantially. Especially the Suc phase of biofilm formation displayed a pronounced number of inter-taxa co-occurrences. In stark contrast, during the Col3, Mat, and Old phases, incidences of inter-and intra-taxa co-occurrences occurred less than would be expected by chance (Supplementary Table 2).
Based on the findings from the co-occurrence networks, we applied a null-model approach to make inferences about the ecological species interactions (cooperation and competition) within the microbial community over the growth phases of the biofilm. Here, we found strong consistencies between the C-score and its variance (C var -score) ( Table 4 , Supplementary Figure 2). During the Col phase, SES values from C-score increased from 3.32 to 45.3 while the corresponding C var -score increased from 3.11 to 285, suggesting increasingly less positive species interaction (i.e. cooperation) over time, while species competition increased ( Table 4 , Supplementary Figure 2). SES values and potential species competition decreased for both scores during the Suc phase followed by a significant increase in the Mat phase (SES: 34.6 and 46.8, respectively). Species cooperation increased again in Old carriers in comparison to the Mat phase, with comparable SES values from 17.8 (C-score) and 14.8 (C var -score), respectively ( Table 4 , Supplementary Figure 2).

Discussion
This study was designed to track and model AMX biofilm microbiota development on synthetic carrier material in a controlled engineered ecosystem over the course of a year. Our purpose was to understand in more detail how microbial interaction shape a functioning anammox biofilm. In accordance with our hypotheses i) and ii), we observed distinct developmental stages in biofilm formation. These stages differed in terms of structure ( Fig. 1 , 2 ) and community composition ( Fig. 3 A). In support of hypothesis iii) these stages were further characterized by dynamic shifts between stochastic and deterministic processes ( Fig. 5 ) and changing co-occurrence patterns ( Fig. 6 ). After a year, the new biofilm converged towards community structures that resembled the mature carrier biofilms already present in the reactor at the initiation of the experiment ( Fig. 3 B). In the following sections, we discuss the ecological dynamics at each stage, with Fig. 7 as a conceptual guideline.

Initial colonization phase (Col1)
In support of our hypothesis i) biofilm formation was initiated by a consortium of heterotrophic bacteria with mostly members of Firmicutes, Bacteroidetes, Proteobacteria, Actinobacteria and, surprisingly, NOB ( Fig. 3 A). Using FISH-CLSM ( Fig. 1 ), we observed a stochastic spatial distribution of microcolonies built by pioneering colonizers in the first month. Presumably, detached bacterial cells Table 4 Observed C-scores and C var -scores, mean metric values under null models, standardized effect sizes (SES) for annamox biofilm communities over time. Higher observed SES values of C-score suggest greater degrees of species segregation than would be expected by chance. Higher SES values of C var -score indicate greater degrees of both species segregation and aggregation. All P -values are < 0.0 0 01.  . 7. Conceptual framework of the biofilm microbiota development on synthetic carrier material in an AMX reactor. Starting clockwise from colonization to maturation with the continuous dispersal from the old carriers at the end. Bar charts denote relative contributions of ecological parameters during each observed development stage.
(red: stochastic processes, green: deterministic processes and blue: co-occurrence frequencies). The number of schematized bacteria approaching the biofilm is representative of modelled immigration rates. The black line denotes the oxic/anoxic interface within the biofilm.
from mature biofilms or from the surrounding wastewater nonspecifically and randomly attached to the pristine carrier surfaces. As newly adhered cells are loosely associated with the carriers, they are likely to easily become detached again as they are subjected to hydrodynamic perturbation (e.g., shear forces) in the reactor ( Fig. 1 ). This is clearly reflected in the substantially higher quantity of exclusive ASVs at the very beginning of biofilm formation ( Fig. 4 ). Furthermore, neutral modeling revealed the predominant role of stochastic processes in the microbiota assembly during this initial colonization stage (Col1) ( Fig. 5 A, 7 and Table 2 ), confirming the very dynamic nature of this biofilm growth phase. However, a large fraction of these biofilm-founding taxa established successfully and became part of the biofilm core community. All observed initial colonization patterns during Col1 agree well with the scientific consensus regarding the fundamental steps of biofilm formation Hall-Stoodley et al., 2004;Jackson et al., 2001 ).

Later colonization phases (Col2 & Col3)
During Col2, stochastic processes still dominated the biofilm formation, but the relative contribution of deterministic processes increased due to higher competition ( Fig. 5 B, Table 4 ). We hypothesize that the microscopically observed stagnation in cluster formation ( Fig. 1 , Month 3), directly resulted from increased competition, resulting in a metabolic trade-off . Pioneering populations are forced to invest some of their metabolic capacity into defense mechanisms rather than in fast growth and spatial expansion. The increase in determinism (i.e, competition) accompanied by a considerable decrease in the immigration rate further supports the notion of the community transitioning from an initially reversible to irreversible attachment, resulting in an enhanced community stability. The drastic decrease in exclusive ASVs in Col3 ( Fig. 4 ) supports this result. While pioneering Phyla like Firmicutes slowly vanished, new key players, such as Chloroflexi, emerged ( Fig. 3 A). Based on previous findings on microbial community succession Datta et al., 2016 ) we speculate that these dynamics are modulated by the availability of substrate to microbes within the biofilm. More precisely, during initial biofilm colonization, labile substrates will be consumed first, supporting copiotrophic microbial taxa, such as Firmicutes Nemergut et al., 2016 ). The latter are then replaced by more oligotrophic (slow growing) taxa (e.g. Chloroflexi, Acidobacteria) that are capable of metabolizing the remaining, more recalcitrant, organic carbon pools ( Rui, Peng and Lu, 2009 ). These may eventually outcompete less adapted pioneering consortia . The filamentous members of the Chloroflexi phylum are known to structurally support biofilm formation ( Speirs, Rice, Petrovski and Seviour, 2019 ) and to scavenge secondary metabolites from generalist decomposers ( Kindaichi et al., 2012;Xia et al., 2007 ). The observed dynamics in community turnover are consistent with the results from the dissimilarity analysis ( Fig. 3 B). They also agree with the significantly increased species competition ( Fig. 5 C, 7 ), as well as frequencies of species co-occurrence associations ( Fig. 6 C, 7 ), suggested by the neutral model and network analyses. Overall, the results indicate a substantial increase in the deterministic influences on the microbiota assembly in the later stages of the Col phase.

Succession phase (Suc)
In developing stream biofilms, after successful surface attachment, cells multiply and continue to produce essential EPS matrix components, forming microcolonies and eventually multilayer aggregates ( Hall-Stoodley, Costerton and Stoodley, 2004 ). The observed expansion of the clusters within the carrier material, indicated by the FISH-CLSM imaging ( Fig. 1 , Table 1 (Month 6-9)), clearly resembles this biofilm formation pattern in natural ecosystems. The increasing number of late transient ASVs ( Fig. 4 , Supplementary Figure 1) and the reduced community dissimilarities ( Fig. 3 B) suggest the convergence of the carrier-associated microbiome towards an essentially stable core community. At this point, competition-driven, deterministic dynamics play only a minor role anymore. Indeed, the Suc phase is characterized by a sharp increase in both the importance of stochastic assembly and immigration rates ( Fig. 5 D, 7 and Table 2 ), and by reduced species competition and co-occurrence ( Fig. 6 D and Supplementary Figure 3). A stable core community now enables transitory attachment and immigration events to dominate the temporal dynamics again.
During the Suc phase, AMX started to emerge within the community ( Fig. 1 , 3 A) which confirmed our hypothesis i) that AMX do not participate in the initial biofilm formation but rather grow slowly into the biofilm when substrate conditions are sufficiently conducive. We argue that slow growth rates ( Kuenen, 2008;Mulder et al., 1995 ) and the high sensitivity to inhibition by various environmental factors ( Jin, Yang, Yu and Zheng, 2012 ) make AMX rather weak competitors during initial biofilm formation. However, once a biofilm has established and anoxic/suboxic zones within the deeper layers of the biofilm have formed ( Flemming et al., 2016 ) anaerobic bacteria like AMX increase in biomass and metabolic performance ( Zhang and Okabe, 2020b ).
During natural biofilm development, complex physical and chemical gradients (available nutrients, oxygen, redox gradients) are established ( Veach et al., 2016 ) and create multiple niches that harbor various cells with different metabolic capacities, which in turn creates opportunities for synergistic interactions ( Flemming et al., 2016 ). We observed an increase in inter-taxa co-occurrences between taxonomically distant bacteria (Supplementary Table 2), suggesting a high degree of mutualistic and commensal metabolic cooperation during the Suc phase ( Ju and Zhang, 2015 ).

Mature phase (Mat & Old)
The near complete filling of the pore space with biomass, and the agglomeration of existing clusters during the last stage of the experimental period ( Fig. 1 , 2 ), in combination with the observed stability in community structure ( Fig. 3 , 4 ), confirm that the biofilm has matured into a stable state. During the Mat phase and in Old biofilms, intra-taxa co-occurrences further increased, which suggests established niches and an increase in the coexistence of closely related species ( Röttjers and Faust 2018 ). According to established concepts of biofilm life cycles, bacterial cells disperse from mature biofilm, revert to planktonic growth and start a new cycle of biofilm establishment on new surface patches ( Flemming et al., 2016 ;McDougald et al., 2012 ). In engineered ecosystems, passive separation from the carrier due to shear forces induced by the constant stirring in these systems is the predominant control on biofilm dispersal ( Garny, Horn and Neu, 2008 ), and increases with increasing biofilm thickness ( Telgmann, Horn and Morgenroth, 2004 ). We argue that these dispersal dynamics also create new ecological niches within the biofilm community, as sloughing now creates local patches that are again susceptible to attachment by bacterial cells ( McDougald et al., 2012 ) or growth of already present community members such as observed for Chloroflexi . The observed increase in immigration rate, number of exclusive ASVs and frequency of co-occurrences in the Old biofilm ( Fig. 6 F, 7 ) support this notion. On the other hand, the marked decrease in the community-wide species competition and increased influence of stochastic processes ( Fig. 5 F, Table 2 , Supplementary  Figure 2) on the community assembly might be the direct result of the aforementioned dispersal dynamics.
In summary, our findings show in considerable detail how an AMX biofilm develops on synthetic carrier material within an engineered ecosystem. Stochastic and deterministic processes structure the microbial community, and lead to a multi-phase succession even in the presence of abundant established biofilms in the reactor. The right choice of reactor conditions has proven its potential to facilitate AMX growth within biofilms ( Agrawal et al., 2017 ;Bhattacharjee et al., 2017 ;Tao et al., 2012 ). Our results suggest, however, that the microbial community has to proceed through the different phases of biofilm formation mandated by internal ecological processes to reach the desired core community structure and metabolic capacities of the mature biofilm. Simple changes in reactor conditions alone appear unlikely to significantly accelerate or modulate biofilm formation as the community members themselves, their mutual interactions and the biofilm structure have to develop over time (in our case over a year) and alter the microenvironment they inhabit. This has considerable practical implications. It may prove difficult to considerably speed up the process of establishing functional AMX biofilms, as each phase depends on previous developments, which puts constraints on the startup of new reactors and the regular replacement of carrier material. It seems plausible that a detailed understanding of the community assembly development can be used to optimize the process. On the other hand, given the observed high degree of microbial interactions in AMX reactors Wang et al., 2019 ) and Fig. 6 ) it seems possible that changes in the reactor configuration might alter the community development, which in turn could directly affect the balance of bacterial interactions within the community in unexpected ways.
AMX generally comprise only a small fraction of the microbial community and are, especially under mainstream conditions very susceptible to disturbances ( Joss et al., 2011 ;Laureni et al., 2016 ;Wells et al., 2017 ). However, they are the key players for one of the most promising energy-efficient mechanisms of fixed nitrogen elimination and the effect of changing community interactions on AMX ecology thus needs to be further explored. Encouragingly, our results indicate that the developing biofilm on the new carriers converged towards the community composition and structure observed in already established biofilms. The mature biofilm exhibited considerable long-term stability both in terms of function and microbial composition, in spite of the potential disturbance effects of compositional or seasonal temperature changes of the incoming wastewater. This may point to a stabilizing effect of the complex interactions in the diverse microbial biofilm community.
Future research needs to apply metagenomic sequencing paired with genome-centric approaches, to provide a deeper understanding of the functions involved in the assembly and the exact roles of individual community members, which will in turn help to manipulate their development.

Conclusions
• In the studied pilot-scale mainstream anammox reactor, stochastic and deterministic processes modulated the microbial community structure within carrier-associated biofilms, and led to a multi-phase microbial succession, ultimately leading to the microbial community structure of already established biofilms. • We argue that a microbial community has to proceed through different phases of biofilm formation (constrained by internal ecological processes) to reach the desired stable corecommunity structure and metabolic capacities of the mature biofilm. • It must be assumed that simple changes in reactor conditions alone are unlikely to accelerate, or modulate the biofilm formation, as internal ecological mechanisms seem to play a critical role in structuring the biofilm community over time. • The staged biofilm growth has considerable practical implications for the startup of new reactors and the regular replacement of carrier materials as each growth phase seems to be a prerequisite for the next one and ultimately leading to a stable community structure. • A deeper understanding of the functions involved in the assembly and maintenance of AMX biofilms is clearly necessary to better understand, and potentially manipulate, their development.

Declaration of Competing Interest
None.

Acknowledgements
This work was supported by funding from the Swiss national science foundation Sinergia project ISOMOL: CRSII5_170876 granted to MFL, HB, AJ, JM and the National Key Research and Development Program of China, Grant ID: Project 2018YFE0110500 granted to F.J. The authors would like to thank Joachim Mohn (EMPA) for the helpful scientific discussions during the course of this study. We would like to thank the ScopeM facility and employees of the ETH Hönggerberg for providing access to their infrastructure and support for cryosectioning and microscopic analysis. We would also like to thank the technical staff of the Versuchshalle EAWAG in Dübendorf and the GDC for providing access to and support of bioinformatics analysis performed on the ETH Zurich Euler cluster. Furthermore, we would like to thank the MicEcol group in Kastanienbaum for the fruitful manuscript discussions.

Additional information
Competing interests: The authors declare no competing interest

Data Availability
Raw 16S sequences can be found on the NCBI sequence read archive under the repository number: PRJNA636402 All other data (gene abundance tables as comma-separated tables and R codes) are available from the Eawag Research Data Institutional Collection (Eric) at https://doi.org/10.25678/0 0 045P.

Supplementary materials
Supplementary material associated with this article can be found, in the online version, at doi: 10.1016/j.watres.2021.117225 .