Free-living and particle-attached bacterial community composition, assembly processes and determinants across spatiotemporal scales in a macrotidal temperate estuary

Bacteria play an important role in biogeochemical cycles as they transform and remineralize organic matter. Particles are notable hotspots of activity, hosting particle-attached (PA) communities that can differ largely from their free-living (FL) counterparts. However, long-standing questions remain concerning bacterial community assembly processes and driving factors. This study investigated the FL and PA community compositions and determinants within the Aulne estuary and the Bay of Brest coastal waters (France). Our results revealed that the FL and PA community compositions greatly varied with salinity and season, explaining a larger part of the variance than the sampling fraction. Both the FL and PA communities were driven by deterministic assembly processes and impacted by similar factors. The FL-PA dissimilarity varied across space and time. It decreased in the estuarine stations compared to the freshwater and marine ends, and in summer. Interestingly, a significant proportion of the FL and PA communities' β-diversity and dissimilarity was explained by cohesion, measuring the degree of taxa co-occurrence. This suggested the importance of co-occurrence patterns in shaping the FL and PA community compositions. Our results shed light on the factors influencing estuarine bacterial communities and provide a first step toward understanding their biogeochemical impacts.

Sampling strategy. The samplings were carried out at spring tide on February 20 (R02), April 18 (R04), July 18 (R09) and November 14 (R12), 2019, as well as on July 20, 2020 (R14), in the context of a monthly monitoring (Urvoy et al., in prep.). Surface waters were collected at nine stations covering the salinity gradient ( Fig. 1). Station 1 (S1) was a freshwater reference site upstream of the dam and was not subjected to tidal influence. Station 2 (S2) corresponded to a freshwater station (salinity 0) under tidal influence and was collected at a fixed location in front of the dam. Stations 3-8 (S3-S8) were sampled every 5 units of salinity (5 to 30), monitored using a Cond330i thermosalinometer (WTW, Weilheim, Germany). Their location consequently varied depending on tides and river discharges. Finally, Station 9 (S9) was a marine reference sampled at the Service d'Observation en Milieu Littoral (SOMLIT, https:// www. somlit. fr/) station of Saint Anne-du-Portzic (48°21′32.17″ N, 4°33′07.21″ W, salinity range: 33.3-34.7). Sampling was carried out from high to mid-tide. The marine station was first sampled and processed immediately. Sampling was then carried out from S1 to S8 within 4 h. Physicochemical variables. Temperature and salinity were measured in situ using a Cond330i thermosalinometer (WTW). Dissolved inorganic nitrogen (DIN), phosphate, silicate, dissolved organic phosphorus (DOP), suspended particulate matter (SPM), particulate inorganic matter (PIM), total particulate phosphorus (TPP), particulate inorganic phosphorus (PIP), chlorophyll a (Chla), and pheopigments (Pheo) were measured as described in the Supplementary methods. Particulate organic matter (POM) was determined as the difference between SPM and PIM and expressed as a percentage of SPM (%POM). Particulate organic phosphorus (POP) was determined as the difference between TPP and PIP and expressed as a percentage of TPP (%POP). The percentage of Chla (%Chla) was determined as Chla over Chla and Pheo.
DNA sequencing and bioinformatics. The BCC was determined using metabarcoding of the V3-V4 region of the 16S rRNA gene. PA communities were recovered using 3-µm filters (Whatman Nuclepore PC membrane). Several filters were used per station depending on the turbidity to avoid clogging. The 3-µm filtrates were filtrated through 0.2-µm filters (Whatman Nuclepore PC membrane) to recover the FL communities. Filters were flash-frozen in liquid nitrogen and stored at − 80 °C until processing. Blank dry filters were sampled simultaneously and used as contamination controls. DNA was extracted using the NucleoSpin Plant  Figure S1). Linear discriminant analysis (LDA) effect size (LefSe) was performed using a LDA score cutoff of 1.5 and a p-value cutoff of 0.05 (microbiomeMarker 52 , v1.1.1).

Figure 1. (A)
Map of the Bay of Brest and the Aulne estuary; sampling dates and their associated tidal coefficient, river flow and mean water temperature. Freshwater (S1, S2) and marine water (S9) stations were collected at fixed locations, while the other stations were variable depending on the salinity. The locations of the April sampling (R04) stations are indicated as an example. Modified from data.shom.fr. (B) Sampling dates and their associated tidal coefficient, river flow and mean water temperature. n: number of sampled stations. n: number of sampled stations. Flow data were retrieved from http:// www. hydro. eaufr ance. fr/ at the Châteaulin station (J3821820). www.nature.com/scientificreports/ Species richness, Shannon index, and Faith's phylogenetic diversity (PD) α-diversity indices were computed (phyloseq and PhyloMeasures 53 , v2.1) on rarefied counts (phyloseq, rarefaction depth = 26,700 reads per sample, rngseed = 999), and differences were tested using the Wilcoxon test (ggpubr). Rarefaction resulted in the removal of two samples (PA communities in S9 of April 2019 and S9 of July 2020) ( Figure S1).
Two β-diversity indices were computed on the read count table transformed to relative abundances (total sum scaling) 54 . Bray-Curtis dissimilarity (phyloseq) considers the ASV abundance and ranges from 0 (no dissimilarity) to 1 (complete dissimilarity). The abundance-weighted β-mean nearest taxon distance (βMNTD) (microeco) considers the mean pairwise phylogenetic distance between each ASV in one community and its closest relative in the second community, weighted by the ASV abundance. β-diversity was visualized using principal coordinate analysis (PCoA, ape 55 , v5.5) applying the Cailliez correction for negative eigenvalues 56 . Permutational multivariate analysis of variance (PERMANOVA) was performed, and variance homogeneity was checked (adonis2 and betadisper test with the Cailliez correction, respectively, vegan package).
The community assembly processes were assessed in the framework developed by Stegen et al. 57,58 based on phylogenetic null modeling. Briefly, the abundance-weighted βMNTD was calculated for all pairwise comparisons of samples. The null βMNTD distribution was quantified by randomly shuffling the ASV name and abundance across the tips of the phylogenetic tree (999 permutations, microeco). The β-nearest taxon index (βNTI) represents the number of standard deviations that βMNTD departs from the mean of the null distribution. A βNTI > 2 or < − 2 indicates that the community assembly is governed by deterministic processes (variable and homogeneous selection, respectively). A βNTI between − 2 and 2 indicates the predominance of stochastic processes. For this range of βNTIs, the Bray-Curtis-based Raup-Crick score (RC Bray ) was estimated (999 permutations, microeco). An RC Bray < − 0.95 suggests the predominance of homogenizing dispersal, an RC Bray > 0.95 suggests dispersal limitation, and an RC Bray between − 0.95 and 0.95 can be interpreted as drift or undominated mechanism 57 .
Cohesion, a measure of the degree of connectivity of a community, was calculated for all samples as proposed by Herren and McMahon 28 . Briefly, the pairwise Pearson correlation of all ASVs was computed and corrected by the "expected" correlation generated from null modeling. The average positive and negative corrected correlations were calculated for each ASV (connectedness metric). The positive and negative cohesion values were computed for each sample by summing the abundance-weighted positive and negative connectedness metrics, thus ranging between − 1 to 0 and 0 to 1, respectively. The absolute value of cohesion increases when the abundances of highly connected taxa increase, reflecting the degree of correlation or connectivity within a community 28 . The cohesion metrics were computed using the authors' script (https:// github. com/ cherr en8/ Cohes ion; persistence cutoff = 0.1; number of iterations = 200; shuffle algorithm = taxa shuffle). As per author guidelines, low-prevalence ASVs (comprising fewer than 150 reads across all samples) were removed beforehand, removing most ASVs (1959 or 12.3% remaining ASVs) but retaining most reads (92.1% remaining reads).
The effect of cohesion and abiotic variables on the composition of the FL and PA communities was assessed using distance-based redundancy analysis (dbRDA) and variance partitioning on the Bray-Curtis dissimilarity (with the Cailliez correction, vegan). The absolute value of negative cohesion was used so that the amount of negative correlations increased with the absolute value of negative cohesion. All predictors were z-score standardized, and multicollinearity was checked using the variance inflation factor (VIF, usdm 59 , v1.1-18). DIN, silicate, and %Chla were thus removed (final VIF < 10). The explained variance was adjusted and the model and axis significance was tested (anova.cca, vegan, 1000 permutations). The effect of cohesion and abiotic variables on the FL-PA pairwise dissimilarity was estimated using redundancy analysis (RDA, vegan) on Bray-Curtis dissimilarity and βMNTD, as described for dbRDA. DIN, silicate, %Chla, and PA negative cohesion were removed (final VIF < 10). Table S1 and plotted in Figure S2. The temperature varied from 6.3 to 23.4 °C throughout the samplings. In February and November, the river end (S1) was colder than the marine end (S9), while the opposite trend was observed on the other sampling dates. DIN and silicate followed the theoretical dilution curve along the salinity gradient (respective Spearman correlation of − 0.93 and − 0.97 with salinity, Figures S2, S3), and both decreased during summer throughout the estuary. In contrast, phosphate did not follow the theoretical dilution curve but increased in S3-S8 compared to S1 and S9, and was higher in July 2019 than on the other sampling dates ( Figure S2). DOP did not follow a clear pattern along the estuary and increased in summer (July 2019, 2020) compared to the other seasons. SPM increased between salinity 5 and 10 (S3-S4), which corresponded to the maximum turbidity zone (MTZ). SPM was especially important in February and April 2019, reaching 197 and 112 mg L −1 , respectively. The MTZ (S3-S4) and downstream estuarine stations (S5-S8) were characterized by a lower %POM and %POP, consistent with the presence of less labile particles compared to the freshwater and marine stations. Chla and %Chla increased in summer (Figs. 2, S2), concomitantly with the increase in temperature and decrease in turbidity that favored light penetration and phytoplankton growth. Consequently, the samplings were differentiated based on salinity and sampling season, with a clear difference between the summer samplings (July 2019, 2020) and the other samplings (Fig. 2).
β-Diversity. The Bray-Curtis-based (taxonomic composition) PCoA showed a marked spatial pattern with a separation between the freshwater (S1-S2) and marine samples (S9) along the 1st axis (16.9% of variance) www.nature.com/scientificreports/ (Fig. 5A, left). There was also a clear seasonal pattern with the 2nd axis (14.6% of variance) separating the summer samples (July 2019, 2020) from the other seasons. The 3 rd axis (9.3% of variance) opposed the end-stations (S1-S2, S9) to the intermediate estuarine stations (S3-S8) (Fig. 5A, right). The communities also differed The Betadisper test showed no difference in dispersion for the sampling time and sampling fraction (p > 0.05), supporting that the PERMANOVA results came from a difference in centroid location. Significant heterogeneity in dispersion was found for the sampling stations (p < 0.05); however, the PCoA plot clearly supported a difference in both centroid location and dispersion.
Bacterial community cohesion. The positive cohesion, reflecting the degree of positive correlations, was overall higher in the PA communities than in the FL communities (Wilcoxon, p < 0.05), especially at the intermediate stations (S3-S8) in February, April and November (Fig. 6). The FL positive cohesion decreased from S1 to S9 (range: 0.19-0.29). The PA positive cohesion was more variable (range: 0.13-0.34). It first increased in the MTZ (S3) (especially in February, April and November 2019), then decreased from S3 to S9, with a sharper decrease at the marine station (S9). Both FL and PA positive cohesion increased in the summer samplings, Community assembly processes. All FL and PA community assemblies were largely dominated by deterministic homogeneous selection (βNTI < − 2, Figure S7), with a mean contribution of 96% to assembly processes. In February, April and November 2019, PA communities were more significantly impacted by homogeneous selection than FL communities, as the mean βNTIs were significantly lower (Wilcoxon, p < 0.005). In contrast, in the summer group, the mean PA community's βNTI was either not significantly different (July 2019) or significantly higher (July 2020) than the mean FL-βNTIs. Stochastic processes also contributed to the FL community assembly by 16 Biotic and abiotic factors driving the FL and PA community compositions. The dbRDA performed on the Bray-Curtis dissimilarity explained most of the variance for the FL (R 2 ajd = 63%, p < 0.001) and PA communities (R 2 ajd = 57%, p < 0.001) (Fig. 7). For both communities, the first two axes clearly separated the sampling stations along the salinity gradient (FL community: 1st axis, 23.1% of variance; PA community: 2nd axis, 12.3% of variance) (Fig. 7A, B). They also separated the summer samplings (July 2019, 2020) from the other  . The summer samplings were correlated with seasonal factors such as temperature, DOP or Chla but also with positive and negative cohesion. Variance partitioning showed that salinity, nutrients and organic matter, and cohesion explained similar amounts of variance on their own (FL: 11%, 9% and 11%, respectively; PA: 6%, 7% and 10%, respectively) (Fig. 7C, D). Temperature alone explained low amounts of variance (FL: 2%; PA: 1%), but an important part was shared with cohesion and/or organic matter and nutrients, reflecting their complex relationships.
Biotic and abiotic factors driving the FL-PA dissimilarity. The RDA (R 2 ajd = 76%, p < 0.001) performed on the pairwise FL-PA Bray-Curtis dissimilarity and βMNTD showed that the increased dissimilarity in the end-stations (S1, S9) was correlated with markers of particle lability (%POM, %POP), while the decreased dissimilarity in the estuarine stations was correlated with phosphate, PA positive cohesion and FL negative cohesion (Fig. 8). The decreased phylogenetic pairwise βMNTD in summer was correlated with seasonal variables (e.g., temperature, DOP, Chla) but also with FL positive cohesion.

Discussion
Characterizing FL and PA bacterial communities and their determinants in estuaries is a first step toward understanding their biogeochemical impact on the organic matter and nutrients exported to coastal ecosystems. Understanding the differences between FL and PA communities is especially interesting since they contribute differently to ecosystem functioning. This study depicted, for the first time, the diversity, composition and driving factors of the FL and PA bacterial communities within the macrotidal Aulne estuary and the adjacent Bay of Brest coastal waters. FL and PA assembly processes were dominated by homogeneous selection. Ecological null modeling showed that FL and PA bacterial community assemblies were largely driven by homogeneous selection, indicating that the biotic and abiotic environmental conditions led to bacterial communities that were more similar than what would be expected by random changes. Few studies have quantified assembly processes for FL and PA communities in coastal ecosystems; however, homogeneous selection was an important driver in other land-sea transition areas [60][61][62] . The predominance of deterministic processes has been suggested to lead to communities being well adapted to their environment, thus potentially having a stronger biogeochemical impact than communities assembled through stochastic processes 63 . This result is consistent with the idea of estuaries being natural bioreactors where organic matter is heavily transformed by estuarine bacteria 32 . FL and PA communities followed pronounced spatiotemporal patterns driven by similar factors. The β-diversity of the FL and PA communities differed largely across space and time and was driven by similar variables with comparable relative importance. Since the sampling stations and sampling time explained www.nature.com/scientificreports/ most of the variance in the FL and PA communities' β-diversity ( Fig. 5 and associated PERMANOVA), the salinity and seasonal variables (e.g., temperature, Chla, DOP) explained a large part of the bacterial community variation (Fig. 7). Salinity is logically an important driver of bacterial communities along fresh to marine water gradients such as fjords, estuaries and coastal margins 16,19,22,23,64,65 . Salinity reflects the mixing of freshwater and marine water masses and directly impacts the BCC by selecting taxa capable of growing at a specific saline concentration. Similar to other studies 16,22,64 , the freshwater communities shifted from a predominance of Actinobacteriota (e.g., Frankiales) and Betaproteobacteria (e.g., Burkholderiales) to a majority of Alphaproteobacteria (e.g., Rhodobacterales, SAR11), Gammaproteobacteria (e.g., Alteromonadales, Oceanospirillales, Vibrionales) and Bacteroidetes (e.g., Flavobacteriales) in the estuarine and marine samples. However, the BCC in the estuarine stations (S3-S8) did not result only from the mixing of water masses: these stations hosted native FL and PA bacterial communities that differed from both the freshwater (S1-S2) and the marine (S9) stations (Fig. 5A, right panels). A native estuarine community can occur only if the doubling time of the bacteria within the estuary is smaller than the flushing time of the water. This is the case in the Aulne estuary, where the bacterial generation times estimated from bacterial production were extremely fast for all sampling times (11-55 h and 1-12 h for the FL and PA communities, respectively, for S3-S8; Urvoy et al., in prep.) compared to the residence time of the water masses (3-30 days 38 ). This reinforces the idea that the Aulne estuary is a biogeochemical bioreactor hosting adapted and active specific communities. Among seasonal variables, water temperature is an important driver of microbial communities and can directly influence BCC and metabolic activity 66,67 . However, temperature on its own explained low amounts of variance, as it was correlated with other variables, such as Chla, %Chla, DOP, and phosphate (Fig. 7), which can select for specific bacterial taxa. For instance, Flavobacteria (Bacteroidota), members of the Roseobacter clade  ) and Gammaproteobacteria (e.g., Alteromonadaceae) have been consistently linked with algal bloom occurrences 68 . In our study, the communities in summer largely differed from the communities in the other seasons. They were less diverse, enriched in autotrophic Synechoccocales and groups usually involved in carbohydrate degradation (Chitinophagales, Verrucomicrobiales or Planctomycetales) 69,70 . The latter can be linked the high Chla and %Chla values measured in summer within the estuary, suggesting the presence of phytoplankton likely impacting both DOM and POM pools. In addition to variables related to space and time, cohesion explained an important fraction of the variance on its own as well as shared with organic matter, temperature or salinity (Fig. 7). This suggested the importance of cohesion in explaining BCC patterns and its link with environmental variables. Negative cohesion reflects the importance of negative co-occurrences, which can result from differences in niches or antagonist interactions such as competition (i.e., the taxa compete for the same resource), among others 71,72 . In contrast, positive cohesion reflects the magnitude of positive co-occurrences, which can result from beneficial interactions (e.g., division of labor to exploit the same resource) or niche partitioning 71,72 . As such, even though they are not appropriate for implying accurate interactions 72 , co-occurrence patterns can drive community evolution by creating feedback loops that can stabilize or differentiate the communities 71,72 . Consequently, cohesion was found to be an important factor explaining BCC in other studies 28,73,74 . In our study, the summer samples contained higher positive and negative cohesion levels than did the other sampling dates, possibly reflecting their importance in structuring the community from winter to summer.

FL-PA dissimilarity differed across space and time and was driven by various factors. Even
though the sampling fraction explained less variance than the sampling station and date, the FL and PA communities were largely dissimilar (Fig. 5). The pairwise FL-PA Bray-Curtis dissimilarity (taxonomic composition) and βMNTD (phylogenetic composition) varied spatially and temporally. Most strikingly, both FL-PA taxonomic and phylogenetic dissimilarities decreased in the intermediate estuarine stations compared to the end-stations for all sampling dates. The pairwise βMNTD also showed that FL-PA dissimilarity significantly decreased in summer (July 2019, 2020) compared to the other seasons throughout the estuary, which was concomitant with a decrease in PA α-diversity. These patterns can be explained by the superposition of several factors (Fig. 8). The higher FL-PA dissimilarity in the freshwater and marine stations (S1-S2, S9) was correlated www.nature.com/scientificreports/ with the presence of more labile POM compared to the inner estuarine stations, as suggested by higher %POM and %POP. Other studies have also shown that the marine station (S9) POM pool was be largely dominated by phytoplankton 75 , while the presence of numerous weirs favored phytoplankton development in the freshwater station (S1) 76 . The colonization of labile POM was suggested to provide a larger advantage than the colonization of refractory terrestrial POM commonly found in estuaries, thus increasing the FL-PA dissimilarity 15,16 . Conversely, the estuarine stations contained less labile POM, possibly explaining the lower FL-PA dissimilarity. Estuary are known to contain refractory POM that originates from terrestrial run-off, which is further transformed by physicochemical (e.g., flocculation, photoalteration) and microbial processes, favored by numerous sedimentation/resuspension cycles 33,34 . In addition, these estuarine stations were characterized by higher PA positive cohesion and FL negative cohesion levels. The increased PA positive cohesion could suggest that colonizing bacteria cooperated to degrade complex, refractory estuarine POM. Interestingly, SPM did not seem to explain the FL-PA dissimilarity (Fig. 8). However, they may contribute to the FL-PA dissimilarity in February, April and November 2019, when strong hydrodynamic forcings may resuspend benthic communities. This would be consistent with the enrichment of anaerobic bacterial groups usually found in marine sediments (Anaerolineales, Desulfobulbales, Desulfobacteriales) on these sampling dates.
The decrease in phylogenetic FL-PA dissimilarity in summer was logically correlated with seasonal variables (temperature, Chla, DOP) and FL positive cohesion. POM and DOM are expected to be more labile in summer since autochthonous processes (e.g., primary production) increase compared to allochthonous inputs (e.g., refractory terrestrial matter run-off), which is suggested by increased Chla and %Chla in summer (July 2019, 2020). Consistently, a different study led in the Aulne estuary has reported a decreased particulate organic carbon to nitrogen ratio in summer (6.6 ± 0.3) compared to winter, where high values (8.6 ± 0.1) were reflective of terrestrial inputs 77 . In particular, algal components were shown to account for 65-80% of POM in summer, compared with 5% in the other seasons 77 . As such, these results seem conflicting: on the one hand, more labile POM was linked to a higher taxonomic and phylogenetic dissimilarity in the end-stations; on the other hand, FL-PA phylogenetic dissimilarity decreased in summer, which most likely contains more labile organic matter. A possible explanation is that the autochthonous labile DOM and POM produced during summer have similar molecular compositions (e.g., phytoplankton components such as polysaccharides or proteins). This may lead to the selection of closely related taxa able to assimilate similar substrates in both fractions and thereby increase the communities' positive cohesion (e.g., increase in positive co-occurrence driven by similar niches). Simultaneously, PA communities are known to solubilize POM into DOM through their extracellular hydrolysis enzymes 6,78 . The resulting DOM can be assimilated by FL bacteria, likely shaping their BCC and favoring co-occurrence patterns between the two compartments. It is possible that these DOM releases increase in summer due to a higher bacterial metabolism, especially within the estuary (Urvoy et al., in prep.). This would likely participate in the homogenization of substrate pools and reduce the FL-PA dissimilarity. The decreased river flow in summer also likely favors the exchange of bacteria and organic matter between both fractions by increasing particles and water residence time.

Conclusion
This study depicted the FL and PA bacterial communities within a temperate macrotidal estuary, the Aulne estuary, for the first time, a step toward understanding their biogeochemical impact. It showed the importance of deterministic processes in structuring both the FL and PA communities along the salinity gradient. Their compositions were driven by salinity and seasonal variables as well as by cohesion, highlighting the need to integrate co-occurrence patterns in such studies. Interestingly, our study showed that FL-PA dissimilarity varied across space and time, and that taxonomic and phylogenetic dissimilarity painted complementary pictures. Given the large taxonomic differences between the FL and PA communities, it is likely that they differ in their metabolic capacities and possibly differentially affect the inputs to the Bay of Brest, although future studies are needed to investigate this hypothesis further.

Data availability
The sequencing data generated during this study are available on online repositories under the accession number PRJNA825348.