Diversity Distribution, Driving Factors and Assembly Mechanisms of Free-Living and Particle-Associated Bacterial Communities at a Subtropical Marginal Sea

Free-living (FL) and particle-associated (PA) bacterioplankton communities play critical roles in biogeochemical cycles in the ocean. However, their community composition, assembly process and functions in the continental shelf and slope regions are poorly understood. Based on 16S rRNA gene amplicon sequencing, we investigated bacterial communities’ driving factors, assembly processes and functional potentials at a subtropical marginal sea. The bacterioplankton community showed specific distribution patterns with respect to lifestyle (free living vs. particle associated), habitat (slope vs. shelf) and depth (surface vs. DCM and Bottom). Salinity and water temperature were the key factors modulating turnover in the FL community, whereas nitrite, silicate and phosphate were the key factors for the PA community. Model analyses revealed that stochastic processes outweighed deterministic processes and had stronger influences on PA than FL. Homogeneous selection (Hos) was more responsible for the assembly and turnover of FL, while drift and dispersal limitation contributed more to the assembly of PA. Importantly, the primary contributor to Hos in PA was Gammaproteobacteria:Others, whereas that in FL was Cyanobacteria:Bin6. Finally, the PICRUSt2 analysis indicated that the potential metabolisms of carbohydrates, cofactors, amino acids, terpenoids, polyketides, lipids and antibiotic resistance were markedly enriched in PA than FL.


Introduction
Microbial communities are highly diverse and play critical roles in biogeochemical cycling cycles across ecosystems, which are fundamental in maintaining climate and ecosystem stability [1][2][3]. In the marine environment, bacteria modulate carbon and nutrient cycling [4] and bacterial community composition is closely related to carbon export [5,6]. Bacterioplankton in the ocean can be classified into two groups based on their lifestyles: Free living (FL) and particle associated (PA) [7]. PA microbes rely on particulate organic matter (POM) as attachment substrate, which constitutes the hotspot microenvironment for the ecological process of marine microorganisms and biogeochemical cycle [8,9]. Indeed, particles of variable sizes, chemical compositions and physical properties conform to the microspatial architecture that structures the microbial environment, which differs from the bulk surrounding seawater in their chemistry and physics [4,10]. There are many different types of particles, such as living or dead protist cells, zooplankton and fish fecal pellets [4]. The patterns of bacteria diversity are known to vary greatly with the size of the particle to Table 1. Criteria to quantify the importance of ecological process accounting for the turnover between communities following Ning et al. (2020).

Homogeneous selection
Consistent selective pressure resulting from consistent environmental conditions is the primary cause of compositional turnover between a pair of local communities <−1.96

Heterogeneous selection
Divergent selective pressure resulting from divergent environmental conditions is the primary cause of compositional turnover between a pair of local communities >+1.96 Homogenizing dispersal High dispersal rates between a pair of local communities is the primary cause of compositional turnover between a pair of local communities <|1.96| <−0.95

Dispersal limitation
Low dispersal rates between a pair of local communities, acting in concert with ecological drift, are the primary cause of compositional turnover between a pair of local communities <|1.96| >+0.95

Drift (and others)
Neither selection nor dispersal is the dominant process driving compositional turnover between a pair of local communities <|1.96| <|0.95| The continental shelf and slope represent a dynamic gradient of terrestrial influences as converging zones of terrestrial currents and ocean basins [45]. There occur intense landsea material exchange, energy transfer and the interaction process between human activities and the natural environment, which play an important role in coupling biogeochemical processes between continents and oceans [46]. Biogeochemical cycles in these regions are closely related to microbial community composition [3]. However, the community assembly mechanisms which underpin microbial community structure in the continental shelf and the continental slope are not well understood.
This study was carried out across the continental shelf and slope in the northern South China Sea, a subtropical marginal sea surrounded by continents, to unravel spatial variability, functions, and processes of free-living and particle-associated bacterial communities. We aim at characterizing the bacterioplankton communities in various pico-to micro-size fractions (free living and particle associated) and describing their variability along horizontal (between the continental shelf and slope) and vertical (from surface waters down into the deep chlorophyll maxima and a deeper layer) gradients in the northern South China Sea. We used high-throughput sequencing to evaluate the diversity and community composition of the bacterioplankton communities and determine the spatial patterns and main drivers of bacterioplankton diversity, community composition, assembly process and potential metabolic functions.

Sample Collection and Environmental Variables
To examine the spatial dynamics of bacterial communities of the costal shelf-slope region, sampling was carried out at the shelf station C6 (117.46 • E, 22.13 • N) and the slope station C9 (117.99 • E, 21.69 • N) in the South China Sea (Figure 1). South China Sea is a marginal sea in the Western Pacific Ocean as it is surrounded by China mainland, Indo-China Peninsula, greater Sunda Islands and the Philippines Islands. Sampling depths included the surface, deep chlorophyll maximum (DCM), which was the depth of 45 m at C6 and 75 m at C9, and for C9 also a deep depth at 150 m (Bottom). The sample codes, locations and study area description are presented in Table S2. For each sample, 40-50 L seawater from more than three CTD rosettes (each with a capacity of 12 L) was Microorganisms 2021, 9,2445 4 of 21 first filtered through a 200 µm bolting cloth to remove large plankton. Two size fractions were obtained by sequential filtration of 200 µm prefiltered water through 3.0 µm and 0.22 µm using 142-mm diameter polycarbonate membranes (Millipore, Billerica, MA, USA) completed within 20 min. Bacteria in the 0.22~3.0 µm size fraction are considered free-living, and those in size fraction 3.0~200 µm are considered particle associated [47]. Finally, the sample-bearing filters were each placed into a 2-mL tube containing 1 mL lysis buffer (0.1 M EDTA and 1% SDS) and stored at liquid nitrogen until DNA extraction. Temperature, salinity and depth were measured using a CTD profiler. A series of nutrients were measured using continuous-flow Technicon AA3 AutoAnalyzer. The detection limits of nitrate + nitrite (NO 3 − + NO 2 − ), nitrite (NO 2 − ), phosphate (PO 3 − ) and silicate (Si 3 O 2 − ) were 0.03 µmol L −1 , 0.04 µmol L −1 , 0.03 µmol L −1 and 0.05 µmol L −1 , respectively.

DNA Extraction and Sequencing
Total DNA extraction was carried out using a CTAB protocol combine DNA Clean & Concentrator kit (Zymo Research Corp, Irvine, USA) as previo [48]. DNA concentrations were determined using a NanoDrop ND-2000 Sp eter (Thermo Fisher Scientific, Carlsbad, CA, USA). The bacterial 16S rRN amplified using the primer pair 515F (5-GTGCCAGCMGCCGCGGTAA-3) GGACTACHVGGGTWTCTAAT-3) [49]. The polymerase chain reaction sisted of an initial denaturation at 98 °C for 30 s; followed by 30 cycles of d 98 °C for 10 s, annealing at 50 °C for 30 s and elongation at 72 °C for 60 elongation step at 72 °C for 5 min. The PCR products were purified using A Pure XP beads and eluted in the Elution buffer. The purified amplicons checked using the Agilent 2100 Bioanalyzer. The validated amplicon sam quenced on MiSeq (Illumina, San Diego, CA, USA) in a barcoded multiplex more than 20 Mbp data output of each sample. The generated raw sequen deposited at the Sequence Read Archive of the NCBI, BioProject number PR

DNA Extraction and Sequencing
Total DNA extraction was carried out using a CTAB protocol combined with Zymo DNA Clean & Concentrator kit (Zymo Research Corp, Irvine, CA, USA) as previously reported [48]. DNA concentrations were determined using a NanoDrop ND-2000 Spectrophotometer (Thermo Fisher Scientific, Carlsbad, CA, USA). The bacterial 16S rRNA genes were amplified using the primer pair 515F (5-GTGCCAGCMGCCGCGGTAA-3) and 806R (5-GGACTACHVGGGTWTCTAAT-3) [49]. The polymerase chain reaction protocol consisted of an initial denaturation at 98 • C for 30 s; followed by 30 cycles of denaturation at 98 • C for 10 s, annealing at 50 • C for 30 s and elongation at 72 • C for 60 s; and a final elongation step at 72 • C for 5 min. The PCR products were purified using Agencourt AMPure XP beads and eluted in the Elution buffer. The purified amplicons were quality-checked using the Agilent 2100 Bioanalyzer. The validated amplicon samples were sequenced on MiSeq (Illumina, San Diego, CA, USA) in a barcoded multiplex fashion, with more than 20 Mbp data output of each sample. The generated raw sequence data were deposited at the Sequence Read Archive of the NCBI, BioProject number PRJNA782430.

DNA Sequence Analysis and Function Inference
We used the QIIME2 pipeline (version 2020.2) to process the 16S rRNA gene [50]. The demultiplex sequences were uniformly trimmed to 250 bp (forward and reverse), and then were denoised to infer Amplicon Sequence Variants (ASVs) that were 462 bp in length using the DADA2 plugin with default settings [51]. Taxonomy was assigned to representative sequences using the Silva 138 Naive Bayes 515F/806R taxonomy classifier. All ASVs affiliated to archaea, chloroplasts and mitochondria, as well as singletons, were removed from the dataset. A phylogenetic tree was generated from representative sequences by aligning sequence fragments via MAFFT, masking ambiguous alignments and inferring a tree using the FastTree algorithm [52]. A rooted tree was created by putting the root at the midpoint of the farthest tips in the tree. To make samples comparable, the feature table was rarefied to a depth of 9000 sequences per sample. We estimated beta diversity using both the weighted and unweighted UniFrac distance in QIIME 2. In addition, PICRUSt2 was used to predict metagenomic functions based on the normalized feature tables [53].

Quantification of Community Assembly Processes
To analyze the relative importance of various ecological processes, variation partitioning, null model and neutral community model (NCM) were performed. The variation partitioning approach was used to evaluate the relative influence of environmental factors, spatial variabilities and their coinfluences on community composition variations. To further elucidate community assembly processes, null model analysis was conducted. The governing processes were determined from the models of variable selection, homogeneous selection, dispersal limitation, homogenizing dispersal and undominated using weighted beta Net Relatedness Index (βNRI) in combination with Raup-Crick metric (RC) [38]. As shown in Table 1, for βNRI, values < −1.96 indicate homogeneous selection, whereas values > 1.96 indicate heterogeneous selection. The RC metric was used to partition the remaining parts with |βNTI| ≤ 1.96. Here, RC < −0.95 represents homogenous dispersal, RC > 0.95 represents dispersal limitation and |RC| ≤ 0.95 represents ecological drift (undominated). Null models were constructed using 999 randomizations, as described by Ning et al. [38].
Moreover, NCM [54] was also used to assess the potential contribution of stochastic processes to the assembly of bacterial communities by modeling and predicting the relationship between ASVs relative abundance and their occurrence frequency across the wider meta-community. The influences of drift (birth-death immigration process) and stochastic dispersal processes can be incorporated in this neutral model. The R 2 values represent the goodness of fit to the NCM, and larger values of R 2 indicate stochastic dispersal and ecological drifts are more important than selection in community assembly by the neutral model [54,55]. Model calculations were performed using the minpack.lm and HMisc packages in R [55,56].

Statistical Analyses
All subsequent analyses were mainly carried out in R (v4.0.2). To investigate the significance of changes in microbial community composition between a priori groups of samples (grouped as shown in Figure 2b and Table S2), non-metric multidimensional scaling (NMDS) ordination analysis of similarity (ANOSIM) tests based on the Bray-Curtis distance were performed using the "vegan" package in R v4.0.2. Single-factor ANOVA and Tukey's HSD test were performed to determine the significance of differences between each sample group with respect to environmental parameters, alpha indices and the relative abundances of specific bacterial taxonomic groups. The alpha diversity indices of bacterial communities, including Pielou, faith_pd, Shannon and Observed_ASVs, were calculated using PAST v3.0 [57]. Two-way analysis of variance (ANOVA) was used to test the significance of season, location and their interaction on environmental factors and alpha diversity indices using PAST v3.0 [57]. The relative abundance of dominant bacteria at different taxonomic levels was envisioned by "ggplot2 package in R v4.0.2. Mantel test was per-Microorganisms 2021, 9, 2445 6 of 21 formed to evaluate the linkages between community structure of the planktonic groups and environmental parameters, which were visualized using "vegan" package. Bacterioplankton with significantly different abundances in different size-fractions were recognized by STAMP based on Welch's t-tests with an adjusted false discovery rate (FDR) [58].
Microorganisms 2021, 9, x FOR PEER REVIEW groups and environmental parameters, which were visualized using "ve Bacterioplankton with significantly different abundances in different sizerecognized by STAMP based on Welch's t-tests with an adjusted false (FDR) [58].

Spatial Variations in Alpha and Beta Diversity
A total of 45 water bacterioplankton communities at different depths nental shelf and slope were analyzed based on the bacterial 16S rRNA ge

Spatial Variations in Alpha and Beta Diversity
A total of 45 water bacterioplankton communities at different depths from the continental shelf and slope were analyzed based on the bacterial 16S rRNA genes (Table S2). Overall, 770,050 tags were obtained in all the samples, 9026 to 19,617 per sample (on average 17,112 tags per sample). After quality control and rarify, a total of 405,000 highquality tags (average of 9000 tags per sample) were retrieved for taxonomic annotation of the bacterioplankton.
As shown in Figure 2a and Figure S2, clear differences were observed between the particle-associated bacterial community (PA, >3 µm) and the free-living community (FL, 0.22~3.0 µm) in alpha diversity. The Shannon-Wiener and Pielou indices of PA were lower than that of FL, while Faith_pd and species richness (observed_ASVs) of PA was significantly greater relative to FL. For FL, the Faith-pd and species richness increased with increasing depth (S < D < B), whereas the Shannon and Pielou indices were the lowest at the surface and the highest at DCM. For PA, the Shannon, Pielou and Faith-id indices decreased with increasing depths, a reverse trend relative to FL ( Figure S2), whereas the species richness had a similar trend to FL ( Figure S2). The two-way ANOVAs with factors for habitat and depth interaction for the alpha diversity indices of bacterioplankton communities suggested depth was a major influencing factor for the Shannon, Pielou and Observed_ASVs, whereas location (station) was mainly responsible for the Faith_pd (Table S3).
NMDS ordination revealed that the bacterioplankton communities were grouped by size fraction (free living vs. particle associated), habitat (shelf vs. slope) and water depth (surface vs. DCM vs. Bottom) (ANOSIM R-value: 0.52 vs. 0.39 vs. 0.30, respectively) ( Figure S3). Meanwhile, the variability of the bacterioplankton communities related to particle size (i.e., size-fractions) was represented by the factor "size-fraction", the coastal to oceanic variability of the communities by the factor "station" and their vertical variability by the factor "depth". The PERMANOVA test revealed statistically significant differences in bacterioplankton community structure due to all three factors (p < 0.001) (Table S1), yet the variability due to "size-fraction" (19%) and "station" (14%) was much higher than that explained by "depth" (6.9%). In addition, the community composition of the smaller size fraction was more variable between depths and stations (i.e., higher beta diversity) than those of the largest size fraction (Table S4).
More importantly, the samples of bacterioplankton formed ten groups that were significantly different for all the microbial communities in an NMDS ordination space built from ASVs-based Bray-Curtis distances (p = 0.001), and the global R values for the entire bacterioplankton among these groups were 0.7176 ( Figure 2b and Table S5). However, SHSF vs. SHDF and SHSP, SLBF vs. SHSP and SHDP, SLDP vs. SLBP had no significant differences (p > 0.05) (Table S5). Overall, the mean Bray-Curtis distance between samples belonging to FL was lower than the mean distance between samples belonging to PA ( Figure S2 and Table S5).

Spatial Variation of Bacterioplankton Community
Gammaproteobacteria (39.02%) was the dominant subphylum in all the microbial communities, followed by Alphaproteobacteria (23.14%), Bacteroidota (11.64%), Cyanobacteria (9.60%), Actinobacteria (3.20%), Bdellovibrionota (1.83%), SAR324_clade (Marine_group_B) (1.45%), Marinimicrobia_(SAR406_clade) (1.43%) and Firmicutes (1.06%) (Figure 3a). We further analyzed the compositions of the FL and PA communities. Significantly higher abundances of Marinimicrobia_(SAR406_clade), Actinobacteria, SAR324_clade (Marine_group_B), Alphaproteobacteria and Cyanobacteria were found in the FL community than in the PA community (Welch's t-test, p < 0.05, Figure 3b). By contrast, the abundances of Gammaproteobacteria, Planctomycetota, Firmicutes, Verrucomicrobiota and Desulfobacterota were significantly higher in the PA than those in the FL community (Welch's t-test, p < 0.05, Figure 3b).  Besides, we also found similar phylum/subphylum distribution patterns acro depth between the FL and the PA communities ( Figure S4). For example, Cyano decreased quickly with increasing water depth. By contrast, the relative abundanc rinimicrobia_(SAR406_clade), SAR324_clade (Marine_group_B), Bdellovibrion Gammaproteobacteria showed continuous increases with increasing depth. The abundance of Firmicutes, Bacteroidota, Gammaproteobacteria, Bdellovibrionota phaproteobacteria were stable across all bacterioplankton communities, with no cant difference from each other ( Figure S4). At the family level, the significantl abundances of Clade_I, SAR86_clade, Actinomarinaceae, Marinimicrobia_SAR4 and SAR324_calde (Marine_group_B) were found in the FL community than in community (Welch's t-test, p < 0.05, Figure S5). The relative abundances of Actino ceae, Marinimicrobia_SAR406_clade and SAR324_clade(Marine_group_B) in the munity showed continuous increases with increasing depth. However, The relati dance of Clade_I showed continuous decreases with increasing depth in the FL Besides, we also found similar phylum/subphylum distribution patterns across water depth between the FL and the PA communities ( Figure S4). For example, Cyanobacteria decreased quickly with increasing water depth. By contrast, the relative abundance of Marinimicrobia_(SAR406_clade), SAR324_clade (Marine_group_B), Bdellovibrionota and Gammaproteobacteria showed continuous increases with increasing depth. The relative abundance of Firmicutes, Bacteroidota, Gammaproteobacteria, Bdellovibrionota and Alphaproteobacteria were stable across all bacterioplankton communities, with no significant difference from each other ( Figure S4). At the family level, the significantly higher abundances of Clade_I, SAR86_clade, Actinomarinaceae, Marinimicrobia_SAR406_clade and SAR324_calde (Marine_group_B) were found in the FL community than in the PA community (Welch's t-test, p < 0.05, Figure S5). The relative abundances of Actinomarinaceae, Marinimicrobia_SAR406_clade and SAR324_clade(Marine_group_B) in the FL community showed continuous increases with increasing depth. However, The relative abundance of Clade_I showed continuous decreases with increasing depth in the FL and PA communities (Welch's t-test, p < 0.05, Figure S5).

Environmental and Spatial Factors Influencing Bacterioplankton Community
The comparisons of environmental factors related to size-fractions and location are shown in Figure S1. Seawater environment factors included nitrite, phosphate, silicate, nitrate, temperature and salinity. The depths we sampled at each station presented pronounced vertical physicochemical variation. Water temperature decreased with increasing depth, from 30.62 • C at the surface to 16.95 • C at the 150 m layer. Salinity and the concentrations of phosphate and nitrate showed a gradual increase with the lowest at surface water and the highest at the 150 m layer, while the concentration of nitrite fluctuated, with an overall decreasing trend, along the whole water column ( Figure S1).
To explore factors influencing the community structure, the mantel test and VPA (Variance Partitioning Analysis) between community composition and environmental-spatial factors were analyzed. Results showed that salinity, temperature, nitrate, phosphate depth and silicate were significantly correlated with the bacterioplankton composition for all samples (  nitrate, temperature and salinity. The depths we sampled at each station presented pronounced vertical physicochemical variation. Water temperature decreased with increasing depth, from 30.62 °C at the surface to 16.95 °C at the 150 m layer. Salinity and the concentrations of phosphate and nitrate showed a gradual increase with the lowest at surface water and the highest at the 150 m layer, while the concentration of nitrite fluctuated, with an overall decreasing trend, along the whole water column ( Figure S1). To explore factors influencing the community structure, the mantel test and VPA (Variance Partitioning Analysis) between community composition and environmentalspatial factors were analyzed. Results showed that salinity, temperature, nitrate, phosphate depth and silicate were significantly correlated with the bacterioplankton composition for all samples (  Based on VPA, the conditional effects of environmental and spatial factors on the FL community were 0.18% and 17.39%, respectively ( Figure 5a). Meanwhile, temperature and salinity were identified as the most important environmental factors shaping the FL community, with conditional effects of 7.21% and 6.15%, respectively (p < 0.05, Figure 5c). The conditional effect of environmental factors on the PA community was 10.72%, whereas spatial factors had no effect (Figure 5b). Moreover, phosphate, silicate and nitrite had significant effects on the PA community, and the conditional effects were 5.99%, 12.47% and 5.23%, respectively (p < 0.05, Figure 5d). Based on VPA, the conditional effects of environmental and spatial factors on the FL community were 0.18% and 17.39%, respectively (Figure 5a). Meanwhile, temperature and salinity were identified as the most important environmental factors shaping the FL community, with conditional effects of 7.21% and 6.15%, respectively (p < 0.05, Figure 5c). The conditional effect of environmental factors on the PA community was 10.72%, whereas spatial factors had no effect (Figure 5b). Moreover, phosphate, silicate and nitrite had significant effects on the PA community, and the conditional effects were 5.99%, 12.47% and 5.23%, respectively (p < 0.05, Figure 5d). When no numbers are shown, this matrix had no explanatory wer (zero or negative). Asterisks represent statistical significance: * p < 0.05, as estimated using canonical correspondence analysis. FL community (a,c); PA community (b,d).

Community Assembly Mechanisms of Bacterioplankton
To explore mechanisms underpinning the observed spatial patterns, the relative roles of niche and neutral processes in community assembly were analyzed. The Sloan neutral model was fitted to the data set to investigate the contribution of stochastic processes to free-living and particle-associated bacterioplankton community assembly. The neutral model explained a large fraction of the relationship between the occurrence frequency of ASVs and their relative abundance variations ( Figure 6), with 64.6%, 62.9% and 60.2% of explained variance for all samples, the free-living and particle-associated communities. Furthermore, the R 2 value ranged from 0.517 to 0.655, with the lowest in the group SHDP and the highest in SLSF ( Figure 6, Table S6). In the same water layer, the neutral interpretation provided a better fit for the continental slope communities than for the continental shelf communities ( Figure 6, Table S6), indicating that the influence of stochastic processes increased in the continental slope compared with the continental shelf. The majority of the ASVs were well predicted by the model, with values from 87.6% to 94.6% ( Figure 6). Concurrently, small variations of the immigration rate m were observed in the whole community, the FL community and the PA community. The estimated migration rate (m value), a measure of the influence of dispersal on community composition, was 0.011, 0.014 and 0.015 for the whole community, the FL community and the PA community, respectively ( Figure 6). The estimated migration rates tended to be higher in the PA than the FL, suggesting that the dispersal limitation had a stronger effect on the PA community than the FL community.

Community Assembly Mechanisms of Bacterioplankton
To explore mechanisms underpinning the observed spatial patterns, the relative roles of niche and neutral processes in community assembly were analyzed. The Sloan neutral model was fitted to the data set to investigate the contribution of stochastic processes to freeliving and particle-associated bacterioplankton community assembly. The neutral model explained a large fraction of the relationship between the occurrence frequency of ASVs and their relative abundance variations ( Figure 6), with 64.6%, 62.9% and 60.2% of explained variance for all samples, the free-living and particle-associated communities. Furthermore, the R 2 value ranged from 0.517 to 0.655, with the lowest in the group SHDP and the highest in SLSF (Figure 6, Table S6). In the same water layer, the neutral interpretation provided a better fit for the continental slope communities than for the continental shelf communities ( Figure 6, Table S6), indicating that the influence of stochastic processes increased in the continental slope compared with the continental shelf. The majority of the ASVs were well predicted by the model, with values from 87.6% to 94.6% ( Figure 6). Concurrently, small variations of the immigration rate m were observed in the whole community, the FL community and the PA community. The estimated migration rate (m value), a measure of the influence of dispersal on community composition, was 0.011, 0.014 and 0.015 for the whole community, the FL community and the PA community, respectively ( Figure 6). The estimated migration rates tended to be higher in the PA than the FL, suggesting that the dispersal limitation had a stronger effect on the PA community than the FL community.
Microorganisms 2021, 9, x FOR PEER REVIEW 11 of 21 Figure 6. Fit of the neutral model for microbial community. The predicted occurrence frequencies for (a-c) bacterioplankton communities representing all, free-living and particle-associated samples, respectively. ASVs that occur more frequently than predicted by the model are shown in red while those that occur less frequently than predicted are shown in orange. Dashed lines represent 95% confidence intervals around the model prediction (black line).

Effects of Lifestyle on Bacterioplankton Community Assembly
The null model based on iCAMP analysis was applied to quantify the contribution of various ecological processes to FL and PA bacterioplankton community assembly. Furthermore, the relative contribution of the observed taxa divided into different groups ('bins') to various ecological processes was also analyzed. The null model results revealed that drift (and others), homogeneous selection and dispersal limitation were more dominant mechanisms shaping the bacterioplankton communities, with an average relative importance of 49.09~61.57%, 13.63~33.61% and 12.68~21.72% (Figure 7a,b), respectively. In contrast, heterogeneous selection and homogenizing dispersal only represented a small percentage of variation in the bacterioplankton communities, with an average relative importance of 0.32%~0.45% and 2.76%~4.17% (Figure 7a,b), respectively. The relative importance of dominant processes between FL and PA bacterioplankton was significantly different (Figure 7c-e). Since other processes had relatively low estimated relative importance (<4.17%), we primarily focused on the effects of homogeneous selection, dispersal limitation and drift on FL and PA bacterioplankton in subsequent analyses.
Overall, the PA community had higher relative importance of drift and dispersal limitation and lower relative importance of homogeneous selection, when compared with the FL community. Significant variations in different groups were observed (Figure 7c-e). The PA community showed significantly higher ratio of drift (Effect size: Cohen's d = −2.19~−0.41, p < 0.05) and dispersal limitation (effect size: Cohen's d = −6.56~0.94, p < 0.05), but lower ratio of homogeneous selection (effect size: Cohen's d = 3.67~5.79, p < 0.05) than the FL community (Figure 7c-e). Drift (and others) contributed less to the community assembly with increasing water depth, however, the opposite trend was found in the PA community at the continental slope station (Figure 7e and Table S7). These results indicate that drift became more dominant in bacterioplankton community assembly with increasing water depth.
Dispersal limitation, which causes divergence in community composition because of a limited exchange of microbes [39], decreased quickly with increasing water depth in the PA community. In contrast, the relative contribute in the FL community continuously increased with increasing depth (Figure 7c and Table S7). Homogeneous selection had a relatively higher importance in the FL community as compared to the PA community. Regardless of whether the free-living or particle-associated community, homogeneous selection increased with increasing water depth in the continental shelf, but an opposite trend was noticed in the continental slope.

Stochastic vs. Deterministic Bacterial Assembly
Based on the principle of the null models employed by iCAMP, the fractions of dispersal limitation, homogenizing dispersal and drift are largely considered stochastic [31]. Fit of the neutral model for microbial community. The predicted occurrence frequencies for (a-c) bacterioplankton communities representing all, free-living and particle-associated samples, respectively. ASVs that occur more frequently than predicted by the model are shown in red while those that occur less frequently than predicted are shown in orange. Dashed lines represent 95% confidence intervals around the model prediction (black line).

Effects of Lifestyle on Bacterioplankton Community Assembly
The null model based on iCAMP analysis was applied to quantify the contribution of various ecological processes to FL and PA bacterioplankton community assembly. Furthermore, the relative contribution of the observed taxa divided into different groups ('bins') to various ecological processes was also analyzed. The null model results revealed that drift (and others), homogeneous selection and dispersal limitation were more dominant mechanisms shaping the bacterioplankton communities, with an average relative importance of 49.09~61.57%, 13.63~33.61% and 12.68~21.72% (Figure 7a,b), respectively. In contrast, heterogeneous selection and homogenizing dispersal only represented a small percentage of variation in the bacterioplankton communities, with an average relative importance of 0.32%~0.45% and 2.76%~4.17% (Figure 7a,b), respectively. The relative importance of dominant processes between FL and PA bacterioplankton was significantly different (Figure 7c-e). Since other processes had relatively low estimated relative importance (<4.17%), we primarily focused on the effects of homogeneous selection, dispersal limitation and drift on FL and PA bacterioplankton in subsequent analyses.
Overall, the PA community had higher relative importance of drift and dispersal limitation and lower relative importance of homogeneous selection, when compared with the FL community. Significant variations in different groups were observed (Figure 7c-e). The PA community showed significantly higher ratio of drift (Effect size: Cohen's d = −2.19~−0.41, p < 0.05) and dispersal limitation (effect size: Cohen's d = −6.56~0.94, p < 0.05), but lower ratio of homogeneous selection (effect size: Cohen's d = 3.67~5.79, p < 0.05) than the FL community (Figure 7c-e). Drift (and others) contributed less to the community assembly with increasing water depth, however, the opposite trend was found in the PA community at the continental slope station (Figure 7e and Table S7). These results indicate that drift became more dominant in bacterioplankton community assembly with increasing water depth.
Thus, the sum of their estimated relative importance can be used to estimate the stochasticity of community assembly. Based on iCAMP results, there showed a significantly higher ratio of stochastic processes in the particle-associated community in comparison with the free-living community (86.04% vs. 66.02%, Cohen's d = 5.772, p < 0.001, Figure  7f), suggesting the bacterioplankton community assembly was even more stochastic for the PA community than the FL community.  Dispersal limitation, which causes divergence in community composition because of a limited exchange of microbes [39], decreased quickly with increasing water depth in the PA community. In contrast, the relative contribute in the FL community continuously increased with increasing depth (Figure 7c and Table S7). Homogeneous selection had a relatively higher importance in the FL community as compared to the PA community. Regardless of whether the free-living or particle-associated community, homogeneous selection increased with increasing water depth in the continental shelf, but an opposite trend was noticed in the continental slope.

Stochastic vs. Deterministic Bacterial Assembly
Based on the principle of the null models employed by iCAMP, the fractions of dispersal limitation, homogenizing dispersal and drift are largely considered stochastic [31]. Thus, the sum of their estimated relative importance can be used to estimate the stochasticity of community assembly. Based on iCAMP results, there showed a significantly higher ratio of stochastic processes in the particle-associated community in comparison with the free-living community (86.04% vs. 66.02%, Cohen's d = 5.772, p < 0.001, Figure 7f), suggesting the bacterioplankton community assembly was even more stochastic for the PA community than the FL community.

Assembly Mechanisms across Different Phylogenetic Groups
In contrast to other approaches, iCAMP allows us to infer information on the relative importance of different ecological processes in individual lineages (bins). Based on the iCAMP analysis, the observed 2515 ASVs were divided into 94 phylogenetic bins. Relative abundances of different phyla in the FL community and PA community, Gammaproteobacteria:Others (20.48% vs. 19.09%) was the dominant phylum/class in the all microbial communities, followed by Cyanobacteria:Bin6 (13.80% vs. 6.16%), Gammaproteobacteria:Bin91 (11.41% vs. 26.30%), Bacterioidota (11.40% vs. 11.85%), Alphaproteobacteria:Others (10.79% vs. 14.55%) ( Figure S6). We further analyzed their relative contributions to homogeneous selection (HoS), dispersal limitation (DL) and drift (DR) for the FL and PA communities ( Figure S6). Higher relative contributions to HoS in the PA community than the FL community, including Alphaproteobacteria:Bin63, Gammaproteobacteria:Bin91, Gammaproteobacteria:Others, Alphaproteobacteria:Others, Planctomycetota, Bacteroidota and Actinobacteriota: Others. By contrast, SAR324_clade (Marine_group_B), Cyanobacteria:Bin6, Alphaproteobacteria:Bin69, Actinobacteriota:Bin30 had more contributions to HoS in the FL community. Besides, we also found that these phyla in the FL community had opposite contributions to Drift and dispersal limitation compared to homogeneous selection compared with the PA community. For example, Gammaproteobacteria:Others had more contributions to drift and dispersal limitation in the FL community than the PA community, whereas Cyanobacteria:Bin6 had more contributions to homogeneous selection in the FL subcommunity than the PA community ( Figure S6).

Predicting Metabolic Functions That Might Underlie Changes between the FL and PA Community Structure
Metagenome inference from the normalized feature tables was performed using PICRUSt2. The weighted nearest sequenced taxon index (NSTI) values ranged between 0.08-0.19 (average 0.12) for the PA communities and 0.1~0.16 (average 0.12) for the FL lineages, indicating that these ASVs were representative of the reference genome database used and could generate reasonable functional predictions. In this study, we inferred the metagenome of the communities using PICRUSt2, and searched for specific genes and metabolic pathways predicted to be differentially abundant between the PA and FL communities ( Figure 8). As shown in Figure 8, more pathways were enriched in the PA community compared with the FL community, consistent with the generally larger genomes and richer metabolic capacities of particle-associated bacteria [59]. The enriched pathways in the PA community were related to secondary metabolites, such as steroid biosynthesis, secondary bile acid biosynthesis, sesquiterpenoid biosynthesis, flavonoid biosynthesis, isoflavonoid biosynthesis, arachidonic acid metabolism and betalain biosynthesis. The additional enriched pathways in the PA community were for Staphylococcus aureus infection, bacterial invasion of epithelial cells and endocytosis. Regarding nutrient metabolism, pathways we found included cyanoamino acid metabolism, which plays a vital role in the nitrogen metabolism of marine microbe [60]. Pathways related to antibiotic resistance were also significantly enriched in the PA community, such as beta-Lactam resistance. Only five pathways were enriched in the FL community compared with the PA community, consistent with their typically smaller genomes. These pathways were related to N-glycan biosynthesis, brassinosteroid biosynthesis, mRNA surveillance pathway, basal transcription factors and the proteasome.

General Structure of Bacterioplankton Community Structure across Size Fraction, Water Depth and Station
As shown in Figure 2a,b and Table S1, clear differences were observed between the particle-associated bacterioplankton community (PA, >3 μm) and the free-living community (FL, 0.22~3 μm). Indeed, the Bray Curtis dissimilarity was highest between the FL and PA communities, followed by the difference with station, with the lowest variability observed between depth (Figures 2b and S3, Table S1). Thus, lifestyle seems to be the most important driver of differential microbial community composition in the continental shelf and slope. As shown in Figure S2, the FL community was more diverse (Shannon indices) and even (Pielou indices) than the PA, although Shannon indices were not significantly different. These results are consistent with previous studies from major ocean gyres [61] and the offshore western Mediterranean [62]. In contrast, the PA community was more abundant (observed_ASVs) than the FL counterparts. We should caution that these differences could, in principle, be due to more 16S rRNA genes per genome in the PA bacteria compared with the FL ones [12].
The microbial community structure varied between size fractions. The FL and PA communities were enriched in different taxa and differed significantly (Figure 3b). For example, Gammaproteobacteria, Planctomycetota, Firmicutes, Verucomicrobiota were enriched in the PA, which is consistent with previous studies, such as those in the Atlantic [27], North Pacific Subtropical Gyre [63], Baltic [21] and Mediterranean Seas [11]. In contrast, Marinimicrobia_(SAR406_clade), Actinobacteriota, SAR324_clade (Ma-rine_group_B), Alphaproteobacteria were enriched in the FL. It is possible that taxa enriched in PA or FL perform different ecological functions. These results agree with the findings in previous studies [12,63]. In addition, the change of the community with sampling depth was clear. Cyanobacteria and Bacteroidetes were relatively more abundant in the surface, whereas Gammaproteobacteria, Marinimicrobia and Bdellovibrionota were relatively more abundant in deeper samples. A similar pattern has been documented previously in the Mediterranean stratified and mixed water columns [64].  Table S1, clear differences were observed between the particle-associated bacterioplankton community (PA, >3 µm) and the free-living community (FL, 0.22~3 µm). Indeed, the Bray Curtis dissimilarity was highest between the FL and PA communities, followed by the difference with station, with the lowest variability observed between depth (Figure 2b and Figure S3, Table S1). Thus, lifestyle seems to be the most important driver of differential microbial community composition in the continental shelf and slope. As shown in Figure S2, the FL community was more diverse (Shannon indices) and even (Pielou indices) than the PA, although Shannon indices were not significantly different. These results are consistent with previous studies from major ocean gyres [61] and the offshore western Mediterranean [62]. In contrast, the PA community was more abundant (observed_ASVs) than the FL counterparts. We should caution that these differences could, in principle, be due to more 16S rRNA genes per genome in the PA bacteria compared with the FL ones [12].
The microbial community structure varied between size fractions. The FL and PA communities were enriched in different taxa and differed significantly (Figure 3b). For example, Gammaproteobacteria, Planctomycetota, Firmicutes, Verucomicrobiota were enriched in the PA, which is consistent with previous studies, such as those in the Atlantic [27], North Pacific Subtropical Gyre [63], Baltic [21] and Mediterranean Seas [11]. In contrast, Marinimicrobia_(SAR406_clade), Actinobacteriota, SAR324_clade (Marine_group_B), Alphaproteobacteria were enriched in the FL. It is possible that taxa enriched in PA or FL perform different ecological functions. These results agree with the findings in previous studies [12,63]. In addition, the change of the community with sampling depth was clear. Cyanobacteria and Bacteroidetes were relatively more abundant in the surface, whereas Gammaproteobacteria, Marinimicrobia and Bdellovibrionota were relatively more abundant in deeper samples. A similar pattern has been documented previously in the Mediterranean stratified and mixed water columns [64].

Different Environmental Factors Affecting the FL and PA Bacterioplankton Communities
Macronutrient profiles showed that from the surface down to the bottom, various environmental factors changed markedly with depth ( Figure S1). These results suggest that the water column was stratified [10,62]. Regarding environmental factors influencing FL and PA bacterioplankton composition and abundance, our study revealed that physical variables (e.g., salinity and temperature), nutrient conditions (e.g., nitrate, nitrite, silicate and phosphate) and spatial factors together contributed to the shifts in bacterioplankton communities in the continental shelf and slope regions (Figure 5a,b). However, the main factors driving the FL and PA community structures were fundamentally different. The variation partitioning analysis (VPA) clearly showed the largest contribution to the variation in the particle-associated bacterioplankton was from environmental factors, whereas the largest contribution to the variation in the FL community was spatial factor. It is possible that the PA community may in general be more metabolically plastic or more responsive to local conditions than the FL [12]. As previously reported, in comparison to the FL, bacterioplankton that dwell on particles exhibit a copiotrophic rather than oligotrophic lifestyle, and are better suited to nutrient-rich and variable marine environments rather than the relatively stable oligotrophic seawater environments [12,46].
In this study, temperature and salinity appeared to play a key role in shaping the FL bacterioplankton community (Figure 5c); conversely, phosphate, nitrite and silicate can explain more community variation in the PA (Figure 5d). The results of our study indicate that FL bacteria are more susceptible to environmental changes (in this case temperature and salinity) compared to PA bacteria, because FL cells are more directly exposed to the surrounding water [65,66]. This can be expected as particles, to which PA bacteria are attached, may act as a "buffer" or micro-island, thereby, PA bacteria are impacted by macronutrients determined by the local environment [65]. Nevertheless, the unexplained contributing part to the variation in the bacterioplankton community was more than 50%, suggesting that other environmental factors, which were not measured in this study, are also important drivers of the bacterioplankton community. For example, a previous study indicated that the combination of the matrices of pigments and depth had significant explanatory power on FL community [67]. Consequently, our understanding of the relative importance of all processes governing community assembly remains incomplete.

Stochastic Processes Dominated Bacterioplankton Community Assembly
Because of the limitation of variation partitioning analysis for inferring ecological processes [68], it is important to interrogate further the ecological mechanisms governing the assembly of microbial communities in community ecology [31,38,54]. In the current study, we used two different conceptual frameworks respectively based on neutral mode and null model, which have been considered useful for investigations of the assembly processes of bacterial communities in diverse environments [38,69,70]. Our results revealed that environmental drift (stochastic process) overwhelmed homogeneous selection (deterministic process) in the bacterioplankton community assembly regardless of freeliving or particle-associated fractions. It agrees with some recent studies that focused on bacterioplankton in estuarine and coastal environments [25]. The null model based on iCAMP showed that homogeneous selection was more responsible for the assembly and turnover of the FL community, while drift and dispersal limitation contributed more to the community assembly in the PA community (Figure 7). That may be because of a higher habitat homogeneity and hydrologic connectivity in the FL community than in the PA community. Homogeneous selection led to the observed reduction of species richness and more similar community composition under consistent environmental conditions [21,71], which is consistent with the lower species richness and Pielou indices in the FL community than the PA community (Supporting Information Figure S2). It may be because of the effect of different light, oxygen and nutrient conditions in size-fractionated samples (free living and particle associated) on the bacterioplankton or due to biotic interactions (e.g., competition, facilitation, mutualism and predation) [68,72]. Stochastic processes, such as dispersal and drift, are important for controlling the assembly of microbial communities in aquatic environments, including rivers, lakes and marine ecosystems [73], where dispersal is relatively unrestricted. In the current study, drift and dispersal limitation contributed more to the community assembly in the particle-associated fractions. It is possible that hydrologic mixing increases dispersal-related processes and ecological drift in open fluidic systems [74,75].
Since the null model approach could more closely reflect the actual stochastic population dynamics, we also performed a neutral community modeling analysis [31,68]. Our finding suggests that neutrally distributed ASVs contributed a large proportion of community richness and thus that community assembly was predominantly driven by neutral processes (Figure 6). The value of the neutral model parameter R 2 was slightly higher in the FL community than in the PA community, indicating that the influences of stochastic processes were stronger for the FL. These observations might be explained by greater habitat homogeneity and hydrologic connectivity in the FL than in the PA, as indicated by the difference of the Bray-Curtis distance between the two fractions ( Figure 2; Supporting Information Figure S3). Moreover, the result of immigration ratio showed that the PA was more dispersal limited than the FL, which was consistent with the result of null model analysis (Figures 6 and 7a,b). However, the neutral interpretation provided a better fit for the continental slope samples than for the continental shelf samples in the FL, just the reverse in the PA, indicating that the influence of stochastic processes increased in the PA and decreased in the FL from the shelf to the slope. These observations might be explained by the decreased contribution of land-derived particles and more variable hydrology in the environment farther away from the coast at a marginal sea, for example, seawater density and temperature, which were the most important environmental modulators of the balance between stochastic and deterministic assembly processes [29].
These two statistical frameworks provide complementary information and have their advantages and limitations. Many previous studies have pointed out that aquatic microorganisms may frequently alternate between a free-living and a particle-associated stage under environmental stress [76,77]. Whether this alternation occurs for the bacterioplankton communities in the South China Sea remains to be investigated in the future.

Predicted Metabolic Potential in the FL and PA Communities
Understanding the functional profiles of bacterial communities is of great importance because it may shed light on ecosystem processes and community assembly mechanisms. In our present study, functional predictions using PICRUSt2 revealed that the bacteria in the subtropical marginal sea in China were involved in many diverse pathways (Figure 8), most of which were the metabolism of carbohydrates, cofactors, amino acids, terpenoids, polyketides, lipids and biosynthesis of other secondary metabolites. However, some specific metabolic pathways were enriched in the PA community. For example, the PA community significantly enriched in genes mediating antibiotic resistance than in the FL community [78]. Compared to the PA community, proteasome was more abundant in the FL community. Studies have shown that proteasome plays a vital role in mycobacterial and actinobacterial stress response pathways in hostile conditions [79]. How this potential function fits the environment being studied remains unclear currently.

Conclusions
In this study, bacterioplankton communities at the continental shelf and slope stations in a subtropical marginal sea were analyzed from 45 samples based on 16S rRNA gene high-throughput sequencing. Both bacterioplankton communities and environmental variables exhibited significant spatial variations. More diverse bacteria were found in the FL community compared to the PA community. Moreover, some members of Marinimicro-bia_(SAR406_clade), Actinobacteriota, SAR324_clade(Marine_group_B), Alphaproteobacteria and Cyanobacteria were more abundantly represented in the FL commuities, whereas Gammaproteobacteria, Planctomycetota, Firmicutes, Verucomicrobiota and Desulfobacterota were enriched in the PA community. Moreover, our results demonstrated that salinity and water temperature were the key factors modulating turnover in the FL bacterioplankton, whereas nitrite, silicate and phosphate were the key factors modulating turnover in the PA counterparts. Finally, our findings suggested that deterministic and stochastic processes simultaneously affected the assembly of bacterioplankton communities, in which the stochastic processes exerted stronger influences on the FL bacterioplankton community, whereas deterministic processes had stronger effects on the PA community. Meanwhile, we found that homogeneous selection was more responsible for the assembly and turnover of the FL community, while drift and dispersal limitation contributed more to the community assembly in the PA community. Importantly, the major contributor to HoS in the PA community was Gammaproteobacteria:Others, whereas in the FL community was Cyanobacteria:Bin6. Addtionally, the major contributors to drift and dispersal limitation in both the PA community and FL community were Gammaproteobacteria:Bin91 and Gammaproteobacteria:Others, respectively. In addition, the predicted metagenomic analysis (PICRUSt2) indicated that the potential metabolism of carbohydrates, cofactors, amino acids, terpenoids, polyketides, lipids, biosynthesis of other secondary metabolites and antibiotic resistance were also markedly more enriched in the PA than the FL community. The findings of differential influencing factors and assembly processes between the FL and PA bacterioplankton communities provide a foundation for further studies to understand the dynamics and ecological niches of these communities in the continental shelf and slope regions.
Supplementary Materials: The following are available online at https://www.mdpi.com/article/10 .3390/microorganisms9122445/s1, Figure S1. An overview of environmental and prokaryotic data at each depth in both the C6 and C9 sampled stations, Figure S2. Distribution of average alpha-diversity at each depth, size-fraction for both the C6 and C9 sampled stations, Figure S3. nMDS ordinations representing the Bray-Curtis distance between bacteiroplankton communities, Figure S4. Relative abundance of top nine bacterial taxa in the Phylum/class level, Figure S5. Relative abundance of top nine bacterial taxa at the family level, Figure S6. Ecological processes controlling major phylogenetic bins, Table S1. Permutational multivariate analysis of variance (PERMANOVA) examining the effects of the factors station, depth and size-fraction on bacterioplankton community structure, Table S2. Group information of all samples, Table S3. Results of two-way ANOVAs with factors for habitat and depth interaction for the alpha diversity indices of bacterioplankton communities, Table S4. Summary of the coefficient of determination (R 2 ) from the Permutational multivariate analysis of variance (PERMANOVA) examining the effects of the factors station and depth for each of the two size-fractions, Table S5. ANOSIM (analysis of similarity) statistics of bacterioplankton community composition between samples grouped following Figure 2, Table S6 Fit of the neutral model for bacterioplankton community, Table S7 Summary of the contribution of the ecological processes that determine community assembly of bacterioplankton Continental Shelf-Slope Region in northern South China Sea.  Institutional Review Board Statement: Not applicable. Informed consent was not deemed necessary since the samples included in this study were used after collection for routine bacterial biodiversity survey in the field.

Informed Consent Statement: Not applicable.
Data Availability Statement: All data generated or analysed during this study are included in this article and its supplementary material files. Further enquiries can be directed to the corresponding author (ED).