Succession of microbial community composition and secondary metabolism during marine biofilm development

Abstract In nature, secondary metabolites mediate interactions between microorganisms residing in complex microbial communities. However, the degree to which community dynamics can be linked to secondary metabolite potential remains largely unknown. In this study, we address the relationship between community succession and secondary metabolism variation. We used 16S and 18S rRNA gene and adenylation domain amplicon sequencing, genome-resolved metagenomics, and untargeted metabolomics to track the taxons, biosynthetic gene clusters, and metabolome dynamics in situ of microorganisms during marine biofilm succession over 113 days. Two phases were identified during the community succession, with a clear shift around Day 29, where the alkaloid secondary metabolites, pseudanes, were also detected. The microbial secondary metabolite potential changed between the phases, and only a few community members, including Myxococotta spp., were responsible for the majority of the biosynthetic gene cluster potential in the early succession phase. In the late phase, bryozoans and benthic copepods were detected, and the microbial nonribosomal peptide potential drastically decreased in association with a reduction in the relative abundance of the prolific secondary metabolite producers. Conclusively, this study provides evidence that the early succession of the marine biofilm community favors prokaryotes with high nonribosomal peptide synthetase potential. In contrast, the late succession is dominated by multicellular eukaryotes and a reduction in bacterial nonribosomal peptide synthetase potential.


Introduction
A central focus of microbial ecology is to determine the forces that govern structural and functional changes of microbial communities, including assembly and succession dynamics found in microbial biofilms [1][2][3][4][5].However, it is only recently that microbial chemical ecologists have begun to acknowledge the many roles microbial secondary metabolites may play in microbial communities [6].In fact, the term "keystone metabolites" was recently used to describe the potential power of these complex molecules in not only microbial communities but also whole ecosystems [7].Thus, to determine the full functionalities of microbial communities and their many different roles, it is essential to understand if and how microbial secondary metabolites shape their development.This implies that a temporal scale must be included in such studies.
Microbial secondary metabolites play a role in several microbial interactions and are important modulators of microbial behavior [8].They can mediate interference competition, virulence, and nutrient scavenging and are produced by a multitude of microorganisms [9][10][11].Much of our understanding of secondary metabolites and microbial community processes has been gained from pure culture or in vitro model systems.While providing essential knowledge about single-handed mechanisms that govern communities (i.e.colonization patterns, interspecies interactions, and phenotypic traits), such studies are constrained to narrow phylogenetic groups and/or performed under conditions far away from the ecological context [12].In contrast, large-scale field studies are still limited, probably not only due to the lack of experimental control in complex natural ecosystems but also due to technology still not allowing the dissection of molecular mechanisms directly in nature.However, the understanding obtained from model organisms on molecular mechanisms of secondary metabolites and the increasing quality of computational tools has allowed us to use genomic traits to predict the biosynthetic gene cluster (BGC) potential, allowing predictions on how microorganisms cope with environmental conditions [13].Most in situ studies on secondary metabolites have focused on profiling the BGC potential of different environments for bioprospecting purposes [14][15][16][17].Efforts within this research area recently revealed that secondary metabolites co-occur in a habitat-specific manner, thus facilitating the distinction of environments on Earth [18].However, these large-scale omics studies only provided "snapshots" of the microbial communities and secondary metabolite potential, hence neglecting the importance of microbial community intrinsic processes and successional changes.One may assume that f luctuations in environmental conditions and community composition will lead to changes in the presence of secondary metabolites.To the best of our knowledge, only one study has examined the changes in the chemical composition of a microbial community through ecological succession [19].Therefore, studies examining the genetic and metabolic response of the microbial communities to different successional phases are scarce.
Here, we analyzed the genetic potential for microbial secondary metabolite production during microbial community assembly in situ and succession in a marine biofilm community.We used a combination of targeted amplicon and genome-resolved metagenomics.We further combined the metagenomic results with metabolomic analyses to track functional dynamics during the succession of a marine community.Besides the fundamental biological understanding, such studies of temporal dynamics can guide more applied bioprospecting, indicating not only the sites but also the times at which the discovery of novel compounds is most likely at a given site.

Results
We developed an in situ system to study microbial community development in a marine biofilm on solid surfaces (BioElements) in a coastal seawater environment.These BioElements were enclosed in nine separated plastic cages; thus, BioElements were taken from each of the nine plastic cages at each time point (Fig. 1A and Supplementary Fig. S1A).We chose to investigate marine microbial biofilms because the emergent properties of microbial biofilms, such as structural rigidity, facilitate microbial interactions and, therefore, have the potential to facilitate microbial secondary metabolite production [20,21].This in situ system allowed us to link the bacterial secondary metabolism dynamics and the community composition over spatially distinct but temporally synchronized BioElements over four months from June to October.At discrete time intervals, we collected BioElements from each cage (one for chlorophyll-a, one for f low cytometry, one for metabolomics, and one for DNA) to assess the temporal dynamics of the biofilm community composition, diversity, BGC profiles, and metabolome.

The temporal abundance of the biofilm microbial community
We quantified the number of microbial cells (cell size <100 μm) detached from the BioElements using f low cytometry (Fig. 1B) to characterize the total microbial abundance and succession dynamics in the microbial communities.Across nine BioElement replicates, the dynamics roughly followed a logistic growth model where the abundance was positively correlated with temperature (spearman rho = 0.37, P-value <.001).In particular, at day 3, the microbial community cell abundance was estimated to mean ± standard deviation (SD) log 5.1 ± 0.3 cells/BioElement, whereafter it increased rapidly until day 29, where it peaked with log 5.9 ± 0.2 cells/BioElement.Subsequently, the microbial cell abundance significantly decreased (Dunnett's test; P-values <.01) and remained below the cell abundance observed from the earlier time points.Based on the temporal patterns in cell abundances, we identified three discrete phases of community succession.In the first phase (Early; day = 3-23), the microbial community steadily increased in absolute abundance.In the second phase (Peak; day = 29), the microbial community "peaked" in cell abundance.The third phase (Late; day = 44-113) was characterized by a decline in microbial cell abundance from day 29 to 44, followed by variable, but relatively low abundances that reached a minimum at day 99 at log 5.25 ± 0.4 cells/BioElement.
We measured chlorophyll-a concentrations as a biomass measure over time to account for photosynthetic organisms on the BioElements (Fig. 1B).The highest chlorophyll-a concentration at 1.8 ± 03 μg/BioElement was measured 3 days after immersion.It was followed by a steep decrease until day 10, suggesting the rapid loss of primary colonizing photosynthetic microorganisms.Over the next two months (days 10-71), the mean chlorophylla concentration increased, resulting in a significant difference between day 10 and day 71 (Wilcoxon test: P-value <.01), thus indicating the increase of photosynthetic biomass over time.However, the SD increased over time, indicating increased variation among replicates (Fig. 1B).

Taxonomic composition patterns during community succession
To assess the prokaryotic and eukaryotic community compositions and diversities over time, we used amplicon sequencing of the 16S and 18S rRNA genes to generate amplicon sequencing variants (ASVs).From the immersion of the BioElements, the prokaryotic ASV richness increased and reached a plateau at day 29 with 1477 ± 215 ASVs (Fig. 1C).The richness remained within this level in the remaining time with the highest richness reached at day 113 with 1928 ± 367 ASVs (Fig. 1C).In contrast, the eukaryotic richness decreased significantly from Days 23 to 44 (Tukey's test: P-value <.01, Fig. 1C).During this decrease, colonization of the bryozoan Conopeum seurati (order: Cheilostomatida; visual identification based on its morphological features, Supplementary Fig. S1B and C) was observed by the naked eye (Supplementary Fig. S1B).The bryozoan colonization was also ref lected in the eukaryotic ASV composition, where a few ASVs, namely, ASV1/ASV3 (class: Copepoda, order: Harpacticoid) and ASV2 (Phylum: Bryozoa, order: Cheilostomatida), increased rapidly in relative abundance after day 29 (Fig. 1D).This trend could also be observed by decreased Shannon diversity after day 29, due to a skewed evenness (Supplementary Fig. S2A).This decline in the eukaryotic diversity and the arrival of the bryozoans along with the copepods indicate a successional shift between Days 29 and 44.This successional shift in the eukaryotic fraction co-occurred with the decrease in microbial cell abundance (<100 μm).This cooccurrence could suggest a shift from "bottom-up" to "top-down" control as previously described [3].
A total of 14.3% of the prokaryotic and 11.3% of the eukaryotic compositional variation were described by the three previously defined community phases (Permutational multivariate analysis of variance [PERMANOVA], P-values <.001), and this was confirmed by the clustering in the principal coordinate analysis (PCoA) (Fig. 1E).Within the early and peak phase for the prokaryotic and eukaryotic communities, the communities clustered clearly according to discrete time points, in contrast to the late phase, where the communities converged (Supplementary Fig. S3).Additionally, the eukaryotic community composition variances within the peak and late phase were significantly larger compared to the early phase (beta-dispersion; P-value <.001; Tukey's t-test; P-values <.0001), while the dispersion within the two phases and the peak for the prokaryotic community compositions was similar (beta-dispersion; P-value = .474).This could suggest that the eukaryotic communities were more susceptible to stochasticity (e.g.temporal and/or environmental changes).Thus, the communities became more distinct from each other over time compared to the prokaryotic communities.Stochastic f luctuations throughout time and genetic variation (drift and diversification) are inherent properties of time-dependent community studies that commonly lead to divergent or convergent communities [22][23][24].
Since the biofilm communities originated from spatially distinct replicates (maximum ∼3 m apart, Supplementary Fig. S1A), spatial location (replicate) could be responsible for the variation observed in community compositions.However, replicates only explained a minor part of the eukaryotic community composition (R2 = 0.01443, P-value = .031)and none of the prokaryotic community compositions (P-value = .091).On the other hand, a significant association was found between the prokaryotic and the eukaryotic composition supported by the Procrustes analysis (M2 = 0.303, P-value = .001).

The early succession phase of the biofilm community possesses a higher nonribosomal peptide potential
To assess temporal changes in the biosynthetic potential of the microbial biofilm communities over time, we used an ampliconbased approach targeting the adenylation domains (ADs) in the nonribosomal peptide synthetase (NRPS) genes followed by clustering into 99% operational biosynthetic units (OBUs) [14,25].
Overall, the biosynthetic potential was higher in the biofilm communities compared to the surrounding seawater both during the early and peak phase (Linear Mixed Model [LMM] and estimated marginal means [EMMs] = 110, P-values <.001, Supplementary Fig. S2B).In the late phase, the biosynthetic potential of the biofilm communities had a higher variability (Tukey multiple comparisons of means; P-values <.001) than in the peak and early phase.This ref lected the bacterial compositional variance across time supported by Procrustes analysis (M2 = 0.359, P-value = .001)and the clustering of the AD OBU composition in the PCoA ordination (Fig. 1E).Interestingly, the AD OBU richness pattern in the developing biofilm did not follow that of the prokaryotic richness (Fig. 1C).This was further supported by the fact that the AD OBU richness was not correlated with the prokaryotic richness (Spearman correlation, rho: 0.02, P-value = .83,Supplementary Fig. S2C).
We hypothesized that the NRPS BGC potential would positively correlate with prokaryotic richness.However, since we did not observe this correlation (Supplementary Fig. S2C), we speculated that the increasing AD OBU richness within the early to peak phase could be explained by a few prokaryotic taxa having a high NRPS potential, and that these were reduced or eliminated from the community during the late phase, hence explaining the significant drop in the AD OBU richness between days 29 and 44-99 in the late phase (Dunnett's test; P-values <.03).Indeed, several ASVs belonging to the Proteobacteria, Bdellovibrionota, Planctomycetota, and Verrucomicrobiota disappeared between peak and late phase (days [44][45][46][47][48][49][50][51][52][53][54][55][56][57], where especially ASVs belonging to the phylum Myxococcota, Actinobacteriota, Spirochaetota, and Fusobacteriota were a pronounced part of the peak phase (Supplementary Fig. S4).Several members within these phyla, especially Actinobacteria and Myxobacteria, have a large biosynthetic potential [26,27].However, the direct linkage of AD OBUs to prokaryotic ASVs is fundamentally not feasible.

Few taxa were responsible for the majority of the nonribosomal peptide synthetase potential
To determine if particular taxa in the biofilm community could be linked to the NRPS biosynthetic potential, shotgun metagenomes were generated to link functionality to taxonomy.Shotgunsequenced metagenomes were generated using coassembly from eight samples originating from eight different days (days 10, 15, 23, 29, 44, 71, 85, 99, and 113), resulting in 200 Gb of raw sequence data.The co-assembly consisting of the nine samples resulted in 444 232 contigs, with an N50 of 5295 bp and more than 2 124 856 open reading frames (ORFs).Saturation curves of ORFs found in the coassembly, particularly the samples from the early phase, were f lattening (Supplementary Fig. S5A), which we suggest is related to low biodiversity in the early phase leading to higher coverage per bacterial genome rather than differences in total sequencing output.To estimate the completion of metagenomeassembled genomes (MAGs), we used the method of singlecopy core genes (SCGs) following the Anvi'o pipeline (https:// anvio.org/).The abundance of 71 SCGs indicated the potential to assemble 340 MAGs, which matched well with the 306 that were recovered (Supplementary Fig. S5B).
A total of 17.5% of the 2 209 unique AD OBUs were mapped to the 306 MAGs with identity >95% and query length > 200 bases, whereas 34.4% of the AD OBUs could be mapped to the full coassembly (Supplementary Fig. S5C).We considered if the unmapped AD OBUs (∼65.5%) were low-abundant OBUs and therefore not represented in the MAGs due to low sequencing effort.However, the relative abundance of AD OBUs was equally distributed across both high-and low-abundant AD OBUs (Supplementary Fig. S5C), suggesting that an unknown fraction of the AD OBUs could be artificial variants generated during the inference of the OBUs or due to Polymerase chain reaction (PCR) biases.Consequently, the amplicon-based method may overestimate the true complexity of AD domains in the marine biofilm (Supplementary Fig. S5C).
To identify MAGs with genes encoding AD domains, we targeted the AMP-binding domains homologous to AD domains [28].In the coassembly, 5 696 ORFs were annotated as AMP-binding domains, and of these, 3655 AMP-binding domains were allocated among the curated MAGs, with the largest proportion found in the peak phase (Day 29; Fig. 2B), corresponding well with the AD amplicon-predicted richness (Fig. 2C).The majority of MAGs had between 2 and 20 AMP-binding domains, yet we found one MAG (JH21_MAG_00132) belonging to an unknown phylum, with 80% completeness, a genome size of 5.65 Mbases, and an N50 of 6217 bases, that contained a total of 123 AMP-binding domains (Fig. 2A and B).This MAG reached 5.41 × Q2-Q3 mean coverage in the peak phase, together with four other MAGs classified as Myxococcota (JH21_MAG_00014, JH21_MAG_00071, JH21_MAG_00065, JH21_MAG_00094), which contained more than 20 AMP-binding domains per MAG but were completely absent after the peak phase (day 29; Fig. 2A).In conclusion, this suggested that specific prolific NRPS producers were disappearing in the late phase.

High abundance of secondary metabolite genes and stress genes characterize the early-phase biofilm
Next, we analyzed the broader BGC potential in the MAGs.The analysis confirmed that the high number of AD domains found in JH21_MAG_00132 and the myxobacteria (JH21_MAG_00014, JH21_MAG_00071, JH21_MAG_00065, JH21_MAG_00094) was associated with NRPS BGCs (Fig. 2C).Besides this, any of these MAGs also had a high genomic potential to produce RiPPs, terpenes, and PKS-like compounds that were significantly allocated to the early/peak phase of the succession in contrast to MAGs with low abundance in the late phase (Fig. 2C).
Besides the BGC potential, the functional genomic potential varied between the early/peak and late phase within specific Clusters of Orthologous Groups (COG) categories (Fig. 3).Differential abundance analysis of the coassembly data revealed that genes encoding defense mechanisms, such as emrA and ercB encoding multidrug eff lux pump [29,30] and salY encoding ABCtype antimicrobial peptide transport system [31], were significantly more abundant in the early-peak phase (adjusted P-value <.05) (Fig. 3).Accordingly, genes related to secondary metabolisms were significantly enriched in the early-peak phase.Specifically, bstA/dinB, encoding Bacillithiol/mycothiol S-transferase [32,33], and hutI, encoding a part of the imidazolonepropionase or a related amidohydrolase biosynthesis [34], were enriched in this period.Lastly, prokaryotic stress defense genes, such as genes encoding outer membrane receptor protein for iron transport (CirA) [35] and genes associated with the universal stress protein (UspA) [36], were also significantly more enriched in the early-peak phase.

The temporal metabolic landscape and detection of secondary metabolites
To complement our genomic prediction of secondary metabolism patterns, we aimed at generating a metabolic landscape of the bacterial communities from Days 10 to 113 using liquid chromatography-tandem mass spectrometry (LC-MS/MS).To   ease the analysis of this complex data, feature-based molecular networking (MN) was used, allowing for features with similar fragmentation patterns to be linked and form clusters or molecular families.Inherent in this analysis is the idea that structurally similar metabolites will fragment in similar fashions, allowing for researchers to pull more out of complex data through relating features in a more automated process.From a total of 795 features detected in the metabolome, 45, 36, and 476 features were unique to the early (n = 23), peak (n = 8), and late (n = 50) phases, respectively ( Supplementary Fig. S6A).A total of 34.36% of the metabolomic feature variance could be described by the first and second components of the PCoA ordination (Fig. 4A).Due to higher beta-dispersion between samples in the late phase compared to the early phase (P-value = .002),the significant effect described by the successional phases on metabolome composition (PERMANOVA; R2 = 0.062, P-value = .002)may be partly attributed to group heterogeneity as well.
Differential analysis revealed that 14 features, including one predicted polyketide (m/z = 300.2169),dominated the early/peak phase compared to 11 features, including one predicted NRP (m/z = 494.3458) in a molecular family with two predicted alkaloids (m/z = 468.33055,m/z = 496.3629),dominating the late phase (Fig. 4B).
Based on the AD amplicon-predicted alpha diversity, we hypothesized that we would find more putative amino acids and peptides during the early/peak phase compared to the late phase.However, the metabolic richness of the putative tripeptides or larger peptides (m/z > 300) were similar across all timepoints (Dunnett's test BH: P-value >.05; Supplementary Fig. S6B) and was not comparable with the amplicon-predicted NRPS richness.Thus, the predicted genetic NRPS potential was either not produced, below detection level [37] or further metabolized [8,38].

Discussion
Microbial secondary metabolites are expected to be important biotic modulators of microbial community populations, especially in spatially structured environments [6][7][8].Despite the acknowledged importance of secondary metabolites in community interactions, little is known about their dynamics and link to microbial community intrinsic processes.Previous studies have described the microbial community succession of marine biofilm [2,3], but have not addressed the potential role of secondary metabolites.In other studies, the secondary metabolite potential in a microbial community has been characterized but only at a single time point [14,15,18,47].Here, we bridge this gap and present the first comprehensive temporal investigation of the interlinkage between secondary metabolite dynamics and metataxonomic composition during community succession in a marine microbial biofilm in situ.We show that the microbial secondary metabolite potential changes during marine biofilm succession and that only a few BGC enriched community members drive the BGC succession pattern.
The in situ biofilm system was used to assess the general succession pattern for the marine microbial biofilm communities.The early succession (<29 days) revealed increasing numbers of microbial cells (<100 μm) until the arrival of multicellular eukaryotes after Day 29.Chlorophyll a concentrations highlighted the prevalence of photosynthetic microorganisms as primary colonizers on the BioElements.This could be explained by the inert surface (plastic) of the BioElements, which does not degrade in the same way as surfaces with organic carbon (e.g.wood).Degradable surfaces typically attract mixotrophs and heterotrophic organisms as primary colonizers, whereas early communities on inert surfaces are primarily composed of autotrophs [3].Subsequently, the microbial abundance and the eukaryotic richness decreased, which could be attributed to the establishment of sessile bryozoans and an increase in benthic copepoda [48,49].Combined, these observations suggest a successional shift around Day 29, which was also observed in the increased variability of chlorophyll a densities and divergence of the eukaryotic community compositions after this time point.
Successional shifts have been described in other microbial systems [1][2][3][50][51][52].They have been suggested to be triggered by both abiotic [53,54] and biotic factors [52,55], resulting in "bottom-up" and "top-down" selective pressures [3].In this study, the successional shift resulted in a significant decrease in the microbial BGC potential, particularly the genomic AD amplicon richness.However, what specific ecological drivers that were causing this succession pattern was not elucidated in our study.It is noteworthy that we observed elevated levels of genes (COGannotated) encoding stress-related functions at the same time, several of which were related to antimicrobial defense mechanisms [29][30][31][32][33][34][35][36].Higher levels of stress-related genes can serve as an indicator of a stressful environment [43,47,48]; thus, these genes may aid in rapid colonization, as surface colonizers require initial acquisition of space and nutrients.However, it is important to note that the higher levels of stress-related genes and defense-encoding genes identified in this analysis could be attributed to both the secondary metabolite producer itself and non-producing strains.We suggest that secondary metabolite producers in the early succession could contribute to the formation of stressful environments (e.g.increased competition and microbial warfare) or assist in coping with abiotic stress-induced conditions (e.g.nutrient limitations).In contrast, the late succession may favor host-associated microbiomes, which provide a relatively less stressful environment, perhaps because it shields the community from abiotic stress factors [56].
Metagenomic analyses revealed that a high portion of the genomic NRPS potential was concentrated in a few bacterial MAGs, supporting previous findings where it was found that a few taxa were highly enriched with BGCs [15].These belonged to an unidentified phylum and species within the Myxococcota class, which suggests that these organisms may be responsible for inf lating the AD potential during the peak phase.Myxobacteria have a huge arsenal of BGCs [57].They can display complex behaviors such as "pack hunting" predation and fruiting body formation that require a high level of cell-to-cell communication and intercellular coordination mediated by their secondary metabolites [58].Evidently, microbial secondary metabolism may be vital for the establishment and survival of bacteria in the environment [9,10,21,55], and secondary metabolite producers, found in marine biofilm, may even serve as initiators for the settlement of invertebrates [59].However, the secondary metabolites may have less impact on the assembly of the invertebrate host-associated dominated microbiomes.
Deciphering the ecological role and function of secondary metabolites has remained a serious challenge, in part due to the in situ detection of the metabolites in the environment as well as low structural annotation rates in untargeted metabolomics [18].Using tandem mass spectrometry, we were able to observe microbial secondary metabolites in situ and track changes in the secondary metabolite feature composition during the community succession.The pattern of secondary metabolites was similar in both succession phases, including NRPS-related features, which contrasts with the trend observed in the sequence-based predicted NRPS potential.Furthermore, the low dispersion of the metabolomes indicated by the PCoA suggested a more conserved microbiome during the early succession compared to the late communities.This observation was further supported by the high levels of stress-related genes in the early phase.
Despite methodological advances, very few metabolites can be detected in the marine environment, and even fewer can be annotated [60].Despite this, we could detect a suite of pseudanes at the successional shift point (Day 29).Pseudanes have multiple biological functions and interactions, including quorum sensing in Pseudomonas [61], antimicrobial activity [45], and acting as an interkingdom signal between Pseudoalteromonas sp. and microalgae [62].The presence of Pseudane IV (also named HHQ) has previously been suggested to be structuring marine microbial communities [6,62].Since we observe Pseudane IV in the marine community, it is plausible that they contribute to community structuring of both prokaryotes and eukaryotes [62].
Although single pseudane compounds have previously been detected in the environment, such as in tunicate extracts (see Supplementary Table S3), our findings present the first evidence of multiple pseudanes being simultaneously detected in situ.However, it is noteworthy that the pseudanes were solely detected in one of the nine replicates, and only at one time point, which could indicate either limitations in chemical detection, or that microbial secondary metabolite production might only occur within a brief time window.Nonetheless, these findings highlight the importance of spatiotemporal resolution when studying microbial assembly dynamics and the secondary metabolite production during biofilm succession.
Our results support similar studies, in that our present metabolite annotation tools are limited when investigating environmental samples, especially from the marine environment [63].Furthermore, bulk chemical extractions may miss low concentration metabolites, which nevertheless could be important for the microbial community dynamics [64].Yet, the addition of metabolomes in microbial community ecology is needed to get a comprehensive and holistic perspective on essential factors governing microbial communities.We therefore argue that future studies should aim to implement metabolomes to strengthen the analysis and enhance the cheminformatic tools.

Conclusion and perspectives
This study has demonstrated that the microbial secondary metabolite potential in marine biofilm communities undergoes succession driven by only a few community members that carry a high genomic potential.This signifies that integration of secondary metabolite investigations into the field of microbial community ecology can aid to a better understanding of intrinsic community processes, such as community succession.Further work should consider a higher spatiotemporal resolution to elucidate the local and temporal effects specifically related to secondary metabolites.

Experimental design, immersion site, and environmental parameters
RK BioElements light (0.93 g/cm 3 ; RK Bioelements, Skive, Denmark) were used as solid surface elements to study the biofilm development.In total, 350 BioElements were added to each of the nine (biological replicates, n = 9) 0.5-L Nalgene ® bottles perforated with multiple holes (diameter ∼1 cm) to allow an exchange of water, thus functioned as cages enclosing the BioElements.The cages were immersed in Jyllinge Harbour, Denmark (55.744920N, 12.094893 E) in June 2021.The attached weights (100 g) allowed the immersion of the cages at 10-20-cm depth.BioElements from cages were sampled on days 4, 7, 14, 28, 42, 56, 70, 84, 98, 112, and 126.Salinity was measured on each sampling day, and water temperatures were measured continuously using a HOBO Pendant ® MX Temperature/Light data logger (Onset Computer Corporation, Pocasset, MA).

Microbial community abundances and chlorophyll a levels
At each sampling time point, two BioElements from each biological replicate (n = 2 × 9) were removed to estimate microbial community abundances by f low cytometry (Supplementary section 01.01) and to estimate the chlorophyll-a concentration (see Supplementary Section 01.02).

DNA extraction
At each sampling point, a third BioElement was removed from each biological replicate (n = 9) in parallel with seawater samples and used for DNA extraction (see Supplementary Section 01.03).

Amplicon sequencing of the 16S rRNA, 18S rRNA, and adenylation domain gene regions
The 16S rRNA V3-V4 region, 18S rRNA V9 region, and NRP ADs were amplified by the primer pairs tagged with octameric barcodes, cleaned, and pooled in equimolar ratios as previously described [25].Primer information is listed in Supplementary Table S1.The amplification was carried out in 75-μl reactions.Further information for each region has been specified in the supplementary information (see Supplementary Section 01.04).All amplicons were cleaned and pooled in equimolar ratios.All sequencing was done by Novogene, Cambridge (UK), on the Illumina NovaSeq 6000 250PE (16S V3-V4 and AD amplicons) or 150PE (18S V9) platform.
ASVs with <10 reads and 16S ASVs with a length below 380 bp were filtered out.Phylogenetic trees were constructed by aligning sequences with Mafft (v.7.490) [67], using a gap penalty of one, and adjusting sequence direction if necessary.FastTree created an approximately maximum-likelihood phylogenetic tree using generalized time-reversible [68].The tree was rooted by automatically picking an outgroup with the longest branch.AD ASVs were clustered into OBUs (99%) using the function TreeLine (method = complete, cutoff = 0.01) from DECIPHER (v2.26.0) [69].

Bacterial amplicon profiling and compositional analysis
A linear mixed-effects (LMEs) model of the relationship between time and environment (fixed effects) on the alpha diversities was conducted using lme4 (v4_1.1-29)[70], with biological replicates as the random effect.P-values were obtained by the Anova function from car (v1.0-9) and post hoc multiple pairwise comparisons using emmeans (EMMs) (v1.7.5) [71], with the inclusion of P-value adjustment (Bonferroni).ASV tables were normalized using total sum scaling (TSS) to 100 000 reads per sample for compositional analysis, using Bray-Curtis dissimilarities.Nonmetric dimensional scaling (nMDS) ordination plots were used to visualize the results.The significance of terms was determined by a PERMANOVA test using adonis2 and betadisper from vegan (v2.6-2) [72].ANCOM-BC (v1.6.2) [73], with default settings, was used to identify 16S ASVs that were significantly different between Days 23-29 and Days 44-57.

Metagenomic sequencing
A subset of the samples was used for shotgun sequencing.The extracted DNA was shipped to Novogene (Cambridge, UK) postextraction for DNA fragmentation by sonication.Library preparation was carried out using Novogene NGS DNA Library Prep Set (Cat No. PT004).Further quality control included qPCR for quantification.Size distribution was detected using Bioanalyzer 2100 (Agilent), and quantified libraries were pooled and sequenced on a NovaSeq 6000 (Illumina) with a 150-bp paired-end strategy.

Metabolome extraction and data acquisition
A BioElement from each biological replicate was used for MeOH extraction.Samples were sonicated in a formic acid (0.1%) MeOH solution for 4 min at 25 • C (28-kHz, 2 × 150 W sonication bath, Delta 220, Deltasonic, Meaux, France).Subsequently, samples were vortexed at max speed for 30 s, evaporated, eluted in 100 μl MeOH, and analyzed on an LC-MS system.LC was performed on an Agilent Infinity 1290 UHPLC system, where 20-μl extract was injected onto an Agilent Poroshell 120 phenyl-C6column (2.1 × 150 mm, 1.9 μm) at 40 • C using CH3CN and H2O, both containing 20 mM formic acid.Initially, a linear gradient of 10% CH3CN/H2O to 100% CH3CN over 10 min was employed, followed by isocratic elution of 100% CH3CN for 2 min.Then, the gradient was returned to 10% CH3CN/H2O in 0.1 min and, finally, an isocratic condition of 10% CH3CN/H2O for 1.9 min, all at a f low rate of 0.35 min/ml.HRMS data were recorded in positive ionization on an Agilent 6545 QTOF MS equipped with an Agilent Dual Jet Stream electrospray ion source with a drying gas temperature of 250 • C, drying gas f low of 8 min/L, sheath gas temperature of 300 • C, and sheath gas f low of 12 min/L.The capillary voltage was 4000 V, and the nozzle voltage was 500 V.The HRMS data were processed and analyzed using Agilent MassHunter Qualitative Analysis B.07.00.HPLC-grade and LC-MS-grade solvents (VWR Chemicals) were used for extractions and LC-MS, respectively.

Metabolomic analysis
For the alpha-and beta-analysis, features with m/z < 300 [79], carbohydrates, fatty acids (NPC pathway predicted), and lignans (NPC superclass predicted) were removed to limit metabolome to potential secondary metabolite production.The data were centered log-transformed after zero-handling with zCompositions pseudo-counts as described in Bates et al. [89] between the early/peak phase (grouped), and the late phase was carried out using the Kruskal-Wallis rank sum test on the non-transformed LC-MC "peak area" across the samples.

Figure 1 .
Figure 1.Microbial community dynamics during marine biofilm succession; (A) experimental setup; nine plastic cages containing plastic elements (BioElements) were submerged in coastal seawater and sampled over four months; at each sampling time point, one plastic element was collected, one from each of the nine plastic cages, to represent the microbial biofilm community at the given time point; illustration is created with Biorender.com; (B) the upper plot shows the absolute abundance of microbial cells (cell size <100 μm) in the collected microbial biofilms; mean and SD are shown as a solid line and error bars, respectively; lower plot shows the surrounding seawater temperature and the mean chlorophyll-a concentration over time; ribbons represent SD; N = 9. (C) 16S, 18S, and AD observed richness over time in the biofilm communities; mean and SDs are shown; each point corresponds to one sample and is colored by the respective phase defined in (B); N = 2-9; (D) relative abundance of top three eukaryotic ASVs per time point in the eukaryotic communities; N = 2-9; (E) community compositions visualized as Bray-Curtis distances in PCoA ordinations; samples are connected to the centroid of the corresponding phase; from left: 16S (n = early: 27, peak: 9, late: 52), 18S (n = 31, 9, 53), and ADs (n = 34, 9, 55); ellipses represent 95% confidence intervals.

Figure 2 .
Figure 2. MAGs found in the marine microbial biofilm community across time; (A) mean coverage abundance of curated MAGs from the coassembly, constituting a subset of samples from 9 days (Days 10, 15, 23, 29, 44, 71, 85, 99, and 113); color indicates phylum classification; N = 1; (B) bubble plot with distribution of AMP-binding domains with total counts of annotated MAGs across time; each bubble represents one MAG colored by phylum annotation; size relates to mean coverage abundance at each time point; the horizontal histogram indicates the sum of AMP-binding domains per time point, and the vertical histogram indicates the distribution of AMP-binding domains per MAG; (C) phylogeny of the recovered MAGs; inner layer: phylum classification; second layer: heatmap indicating the mean coverage abundance of the MAGs across time; outer layer: antiSMASH predicted BGC counts colored by BiG-SCAPE class.

Figure 3 .
Figure 3. Functional analysis of stress-related genes; differential abundance of functional stress-related genes; genes with negative Log 2 fold change were more abundant in the early-peak phase (n = 4), and genes with positive Log 2 fold change were more abundant in the late phase (n = 5); each point represents a functional gene, colored by COG category.

Figure 4 .
Figure 4. Ordination and differential analysis of the metabolomic LC-MS features; (A) metabolome composition of features >300 m/z visualized as Euclidean distances in PCoA ordinations; ellipses represent 95% confidence intervals; each sample is connected to the centroid of the corresponding phase (early, peak, or late); N = 23, 8, and 50; (B) differential analysis of metabolic features dominating the early/peak (n = 29) vs. late phase (n = 49) (Kruskal-Wallis rank sum test; comparisons with adjusted BH P-values <.05 are shown); the LC-MS peak area of each feature in the early/peak (left side; circles) vs. late phase (triangles) colored by their predicted NPC pathway; error bars represent SDs, and a circle/triangle represents the mean peak area of each feature.