Trophic States Regulate Assembly Processes and Network Structures of Picophytoplankton Communities in Subtropical Coastal Ecosystems

Background: Long-term coastal eutrophication especially in semi-enclosed marine areas is driven by increased amounts of nutrients derived from anthropogenic activities. Given that accelerating nutrients may constitute a strong environmental lter, understanding the diversity, assembly process and co-occurrence pattern of picophytoplankton communities in response to increasing coastal eutrophication is clearly of great importance. Results: We investigated picophytoplankton community changes using rbcL gene amplicon sequencing. The results exhibited that the alpha diversity (ANOVA, p < 0.001) and beta diversity (ANOSIM, p < 0.001) were signicantly different among eutrophic states. Further, phylogenetic based β-nearest taxon distance analyses revealed that stochastic processes mainly provided 69.26% contribution to picophytoplankton community assembly, whereas deterministic processes dominated community assembly in a heavy eutrophic state. Integrated co-occurrence networks modularly responded to eutrophic states and revealed that keystone taxa mainly belonged to the oligo eutrophic group, which may play fundamental roles in network persistence. Importantly, increased environmental disturbances, such as nitrogen and phosphorus nutrients, could alter picophytoplankton community structure and disrupt ecological processes. Conclusion: Stochastic and deterministic processes simultaneously inuenced the assembly of picophytoplankton communities in the subtropical coastal ecosystems. Eutrophic disturbances alert the assembly processes and network structures of picophytoplankton community. Our ndings promote the understanding of fundamental ecological processes along eutrophic gradients in subtropical coastal ecosystems.

Modern molecular methods have revolutionized our knowledge of assessment of detailed picophytoplankton communities and diversity, and 16S rRNA (for pico-cyanobacteria) [7] and 18S rRNA (for pico-eukaryotes) [8] represent useful tools for microbial ecologists to study microbial diversity [9,10].
Besides, functional gene markers, especially rbcL, which encodes the large subunit (L) of ribulose-1,5bisphosphate carboxylase/oxygenase (RubisCO) involved carbon assimilation [11,12], have been extensively used and provide increased resolution for assessments of biodiversity on picophytoplankton community across different coastal environments, including the bay of Bengal [13], Western Iberia [14], Western Mediterranean [15] and Northern South China Sea [9]. Although the general patterns underlying variations in biodiversity have been observed across spatial and temporal scales, the effects of eutrophic conditions that controlling these patterns on picophytoplankton communities remain unclear in a subtropical coastal ecosystem.
Increasing evidence suggests that the capacity of microbial community to resist environmental disturbances is an emerging property arising from the multiple interactions between them, such as competition and/or mutualism [16][17][18]. Intensive interactions in picophytoplankton community undergoing environmental disturbances can markedly determine species co-occurrences patterns. A complex co-occurrence network can reveal intrinsic microbial interactions and modular structures in response to environmental disturbance that may play important roles in governing community assembly [19]. Furthermore, some central or keystone species and environmental factors could be identi ed within an ecological network that may play essential roles in maintaining the structure and governing community assembly. Recently, most studies have focused mainly on co-occurrence networks among bacteria or eukaryotes during picophytoplankton blooms to reveal inter-taxa associations in microbial communities [1,20,21]. To date, there are few studies related to the impact of nutrient level dominance on the intrinsic interactions and structure among picophytoplankton networks along eutrophic gradients.
Moreover, the assembly mechanisms (deterministic and stochastic processes) that are involved in jointly shaping microbial community variation are also intensively in uenced by nutrient level [22][23][24]. A previous study proposed that pH and temperature affect the balance between stochastic and deterministic processes to shape microbial community assembly [25], whereas studies show that nutrient input contributed to a gradual shift from stochastic to more deterministic-based processes in coastal ecosystems [26,27]. Therefore, to advance the understanding of community assembly processes, it is critical to elucidate the balance between stochastic and deterministic processes underlying changing environmental conditions. Additionally, there exists limited knowledge on what underlying nutrient factors lead to the dominance of particular assembly processes for picophytoplankton community in subtropical coastal ecosystems along eutrophic gradients.
As an estuarine environment, the Beibu Gulf is strongly affected by anthropogenic activities, leading to trophic conditions that are more complex and dynamic compared to other ecosystems [28,29]. The addition of nutrients can modify biogeochemical processes and alter microbial communities [30]. In addition, this region exhibits the remarkable nutrient gradients that could provide an ideal simpli ed environment to study the effects of eutrophication on picophytoplankton communities. Here, we explored the diversity and associations of picophytoplankton communities using the rbcL gene based on highthroughput sequencing and co-occurrence network analysis in the subtropical gulf. The goals of this study were to (i) determine how the diversity and composition of the picophytoplankton community vary along eutrophic gradients in coastal ecosystems, (ii) uncover the co-occurrence patterns of picophytoplankton communities, and (iii) explore the underlying nutrient factors regulating community assembly processes along eutrophic gradients.

Nutrient levels at locations
We measured a number of climatic and nutrient factors to reveal the extent of environmental heterogeneity across sampling sites (Table S1). The environmental parameters of the seawater samples all exhibited characteristic differences and varied signi cantly between the stations sampled from various groups (one way-ANOVA analysis, p < 0.001) (  Fig. 1b), which is consistent with the change in DIP. These results suggested that nutrient levels were markedly different among the locations that were divided into three groups that were based on similar intensities of anthropogenic activity, such as intensive aquaculture. Moreover, other environmental parameters at the time of sampling for all the stations are detailed in Table S1. 2.2 Alpha diversity of the picophytoplankton community across eutrophic states Across all water samples, we obtained 6,299,765 high-quality sequences after normalizing that were clustered into a total of 855 OTUs after quality control at a 90% similarity level (Table S3). Good's coverage ranged from 99.79-100% in all the samples, indicating that the current sequences re ected the actual situation of the majority picophytoplankton community. The richness and diversity of picophytoplankton communities in different groups were compared based on OTU richness and phylogenetic diversity. Picophytoplankton α-diversity, as expressed by the Shannon and Simpson diversity indices, was the highest in the HE and decreased along eutrophic gradients. We also found that α-diversity were signi cantly altered between groups (ANOVA analysis, p < 0.001), and remarkable differences were observed in these two indices between HE and OE (Fig. 2). However, no signi cant differences in these two indices were noted between Chao1 and Observed number indices. Furthermore, the effects of geographical distance and other environmental factors on the alpha diversity (Shannon and Simpson indices) of the picophytoplankton community were compared using contours map and Spearman's correlation, respectively ( Fig. 1c & Table S4). According to the Spearman's correlation analysis, both OTU richness and alpha diversity of the picophytoplankton community were signi cantly negatively (p < 0.01) correlated with salinity (r = -0.46), TOC (r = -0.43), pH (r = -0.42), and COD (r = -0.39) but also positively (p < 0.001) correlated with NH 4 + (r = 0.58), DIP (r = 0.52), TDP (r = 0.5), and DIN (r = 0.42). Furthermore, EI was consistently the best predictor for Shannon diversity with a strong positive signi cant correlation (r = 0.58, p < 0.001) (Table S4).

Composition and beta diversity of the picophytoplankton community across eutrophic states
The dominant microbial phyla included Bacillariophyta and Rhodophyta, accounting for more than 80% of the total sequences. The group related to heavy eutrophication was signi cantly over-dominated in the class of Bangiophyceae, and Coscinodiscophyceae, whereas the medium and oligo eutrophication groups were mostly a liated with Coscinodiscophyceae and Bacillariophyceae ( Figure S1). Nonmetric multidimensional scaling (NMDS) analysis based on Bray-Curtis dissimilarity showed that the picophytoplankton community compositions derived from three groups exhibited a clear separation across the eutrophication gradient, with samples being more dispersed at oligo eutrophic level (Fig. 3a). A higher beta diversity was found in the OE compared with the other two groups. This pattern was further supported by the results of an analysis of similarity (ANOSIM, R = 0.863, p = 0.001), which demonstrated that the Bray-Curtis dissimilarity rank of between group was higher compared with that for each of the three subgroups. These results indicated that the difference in picophytoplankton community composition was more signi cant between groups than within groups (Table S5 & Fig. 3b). Moreover, we also calculated the effects of environmental variables on the beta diversity of picophytoplankton community using Mantel tests ( Table 1). The results suggested that almost all environmental factors were strongly correlated with picophytoplankton community structure, and shifts in picophytoplankton community structure may be due to main environmental factors, such as DIP (r = 0.520, p < 0.001), TDP (r = 0.520, p < 0.001), geographic distance (r = 0.488, p < 0.001), temperature (r = 0.447, p < 0.001), and pH (r = 0.435, p < 0.001) ( Table 1). The results revealed that P nutrient changes maintained stronger effects on the picophytoplankton community composition and exhibited a signi cant positive relationship with differences in Bray-Curtis dissimilarity which indicates that the larger the P nutrient difference, the more dissimilarity between the picophytoplankton community structure. To gain further understanding of the linkages between picophytoplankton community dissimilarity and nutrient changes, we explored linear regressions based on Bray-Curtis distance (Fig. 4). The results showed that geographic distance is an important factor that elicits variation in picophytoplankton community structure which had a strong positive linear relationship with Bray-Curtis dissimilarity (regression, R 2 = 0.412, p < 0.001) (Fig. 4a). DIP and TDP (regression, R 2 = 0.594, p < 0.001) also had signi cant effects on picophytoplankton community structure, and the effect of phosphorus level was stronger than that of geographic distance (Fig. 4).
Overall, these observations strongly suggested that nutrient levels, especially phosphorus, were major drivers in shaping the structure and diversity of the picophytoplankton community. Notably, the dissimilarities among the picophytoplankton communities were also signi cantly correlated with eutrophic states (Fig. 4d), suggesting that the magnitude of changes in picophytoplankton communities is dependent on the selective pressures exerted by increasing eutrophication.
Given that picophytoplankton communities were dramatically altered by increasing eutrophication, we then asked whether picophytoplankton assemblages could be indicative of eutrophic states. Thus, we screened 18 picophytoplankton families for which their relative abundances were signi cantly associated with EI as depicted in the heatmap (Fig. 5). However, the responses of these assemblages to increasing eutrophication were divergent. The majority of picophytoplankton families were signi cantly different as demonstrated nutrient gradient distribution patterns. For example, assemblages a liated with Bacillariophyceae (including the Achnanthidiaceae, Mastogloiaceae, Amphipleuraceae, and Naviculaceae) were signi cantly negatively correlated with eutrophic states, whereas assemblages belonging to Coscinodiscophyceae (Coscinodiscaceae and Lauderiaceae) were signi cantly positively correlated with eutrophic states. Cluster cladogram analysis in the heatmap illustrated that the samples grouped together with three eutrophic states that were clustered separately. These results suggested that the variation in picophytoplankton community varied markedly greater across eutrophic levels. Thus, these sensitive assemblages have the potential to be used as bio-indicators of eutrophic states.

Co-occurrence patterns of the picophytoplankton community and environmental variables
In addition to simple composition and diversity, co-occurrence network analysis can potentially provide a deep and unique perspective on picophytoplankton interactions and ecological assembly rules that are constructed based on Spearman's rank correlation (Fig. 6). The picophytoplankton co-occurrence networks consisted of 240 nodes linked by 4432 edges (Table S6 & Fig. 6c). For the picophytoplankton co-occurrence network, Bangiaceae (15.42%), Thalassionsiraceae (14.17%), Naviculaceae (9.58%), Skeletonemataceae (6.67%), Lithodesmiaceae (6.67%), Bacillariaceae (5.42%), and Chaetocerotaceae (5.42%) mainly occupied the nodes (Fig. 6b). In the picophytoplankton co-occurrence network, 30 OTUs with high degrees and eigenvector centrality were identi ed as keystone species, and these species were signi cantly preferred given that almost all were distributed mainly in the OE group ( Figure S1). In addition, the top 10 of keystone OTUs mainly belonged to Naviculaceae, Lithodesmiaceae, Mastogloiaceae, Cymatosiraceae, Rhizosoleniaceae, Amphipleuraceae, and Catenulaceae ( Figure S1). Given that OTUs from the same module co-occurred more frequently than nodes in other modules, cooccurrence network analysis also shows that nodes from the same family or class may play an important role in determining modular structures. Consequently, the picophytoplankton network was clearly divided into four major modules, of which modules I, II, III, and IV accounted for 34.58%, 43.33%, 16.37%, and 6.67% of the whole network, respectively (Fig. 6a). Module I was primarily occupied by Coscinodiscophyceae and Bacillariophyceae (including Thalassionsiraceae, Skeletonemataceae, and Naviculaceae), whereas module II was mainly composed of Bangiophyceae. Ternary plot analysis suggested that these modules were speci c to a particular eutrophic group ( Figure S2). For example, modules I and III were speci c to the OE group, whereas modules II and IV were speci c to the HE group and ME group. Furthermore, the observed modularity, clustering coe cient and average path length were all greater than those of their respective Erdös-Réyni random networks, indicating the picophytoplankton network was non-random with "small-world" properties and modular structure (Table S6). For the integrated networks, the variation of correlations was also signi cantly correlated with the Euclidean distance of environmental factors. Salinity and pH were mainly attached to module I, which was associated with oligo eutrophication, whereas nutrients clustered into module II were associated with heavy eutrophication. DIP, TDP, and EI exhibited higher degrees and eigenvector centrality compared to other environmental factors, indicating that nutrient levels had stronger relationships with picophytoplankton OTUs compared to other environmental factors in the Beibu Gulf (Fig. 6c).

Assembly processes of picophytoplankton communities along eutrophic gradients
To discriminate between the deterministic and stochastic processes in shaping picophytoplankton community along eutrophic gradients, we compared the different assembly processes of picophytoplankton community in all samples and the three groups using within-community (nearesttaxon index [NTI]) and between community (βNTI) assessments (Fig. 7). For the overall group, stochastic processes, which encompassed dispersal limitation (10%), homogenizing dispersal (36.34%), and ecological drift (22.92%), contributed a larger fraction (69.26%) to the community assembly, indicating that stochastic processes dominated community assembly. In contrast, deterministic assembly processes (80%) mainly in uenced the picophytoplankton community in the HE group. In different groups, the relative contribution of stochastic processes varied from 20% in HE and 87.59% in ME to 99.23% in OE. We also found the NTI values decreased linearly in relation to increasing DIN and EI, indicating that the community assembly tended to strengthen signi cantly with increasing eutrophication (Fig. 8a & 8c). To infer alterations in the deterministic/stochastic assembly processes along nutrients gradients, we examined the relationships between βNTI and environmental variables using Mantel test analyses (Table 3). Pairwise comparisons of βNTI values for the picophytoplankton communities were mainly signi cantly positively correlated to changes in DIN (r = 0.665, p < 0.001), NO 3 − (r = 0.635, p < 0.001), salinity (r = 0.520, p < 0.001), EI (r = 0.520, p < 0.001), NH 4 + (r = 0.469, p < 0.001) and temperature (r = 0.300, p < 0.001), indicating that increases in their divergence were associated with a shift in picophytoplankton community assembly processes from stochasticity to heterogeneous selection (Table 2). Additionally, the deterministic processes that shaped the community assembly in the HE group tended to strengthen signi cantly (Mantel-High, EI, r = 0.966, p < 0.001) with increasing coastal eutrophication. Overall, the regression models showed that most of the βNTI values were between − 2 and 2, indicating a dominant role of stochasticity in the entire picophytoplankton community (Fig. 8b & 8d) that was signi cantly affected by nitrogen concentration, salinity, and temperature.   Increasing evidence has shown that disturbance could select species exhibiting speci c biological traits and markedly alter the composition of picophytoplankton community [31,32]. A high relative abundance of Rhodophyta was observed and widely spread in the heavy eutrophication Maowei Sea, whereas Bacillariophyta was signi cantly over-dominated in the Qinzhou Bay and Beibu Gulf open sea as eutrophic states decreased, which is consistent with previous studies [33,34]. Bacillariophyta exhibited high adaptability for oligotrophic conditions as to be dominant in nutrient-de cient areas [35]. This assertion is supported by our previous nding that nutrient variables exhibited clear spatially structured patterns. Additionally, the variation in the relative abundance of 18 picoplankton families closely correlated with eutrophic states, indicating that sensitive assemblages could be considered as bioindicators for evaluating eutrophic states (Fig. 5). For example, the Bangiaceae family can e ciently cope with nutrients from nearshore eutrophic coastal areas for both absorption and growth, which is consistent with their ability to better adapt to nutrient-enriched conditions [36]. In contrast, the high relative abundance of most oligotrophic assemblages has been extensively detected and was negatively correlated with EI (Fig. 5). In general, changes in the patterns of these sensitive assemblages were signi cantly associated with eutrophic states, which is consistent with their ecological strategies.

Drivers of alpha and beta diversity in picophytoplankton communities
In this study, we used Spearman's correlation to investigate the effects of environmental variables on the picophytoplankton community and found that Shannon and Simpson diversity exhibited signi cant positive relationships with nitrogen and phosphorus concentration but negative relationships with pH and salinity (Table S4). There are two possible explanations for such signi cant effects. One possible explanation for this positive effect could be attributed to the fact that the accumulation of nutrients increased the nutritional pressure to select tolerance of many species, thereby increasing the community alpha diversity. The high diversity in the HE group may increase the functional redundancy of the picophytoplankton community and further provide biological buffering capacity to resist environmental changes [37]. Another possible explanation for this negative effect was that the increasing salinity prevented taxa from adapting to the enhanced osmotic pressure, so these taxa may die or become inactive [38], thus reducing alpha diversity.
Consistently, we found that increasing eutrophic states drove gradually stronger changes in picoplankton community composition and exhibited a clear separation from oligotrophic-dominated to eutrophicdominated assemblages (Fig. 3). A higher beta diversity was observed in the OE compared with the other two groups, and the OE also exhibited the lowest alpha diversity (Figs. 2 & 3). This nding may be attributed to the fact that many taxa of Bacillariophyta prefer the relatively low turbidity and low nutrient concentrations in the seawater [39]. As expected, community dissimilarity between different groups exhibited signi cantly greater variances, indicating greater variances within the communities that exhibited lower alpha diversity. Interestingly, moderate eutrophication caused the selectivity and stability of more convergent picophytoplankton communities, thereby resulting in reduced beta diversity (Fig. 3). Additionally, the HE also exhibited lower community dissimilarity, which may be due to aquaculture activities; for example, shell sh farming not only controlled the size of picophytoplankton populations but also the growth of certain algal species through their selective ltration [40]. Thus, a reduction in the impact of aquaculture activity is an effective method to control the occurrence of microalgae blooms in such a semi-closed and intensively eutrophic gulf. As previous studies reported [9], relatively high temperature stimulated both photosynthesis and respiration, thereby impacting the biodiversity of phytoplankton at the sea surface given its signi cant effects on the metabolic rate of picophytoplankton [41]. We observed that the variation in picophytoplankton community composition was related to a direct in uence of temperature, salinity, nitrogen and phosphorus, which displayed stronger signi cant correlations based on mantel test analysis (Table 1). A previous study also found that picophytoplankton community diversity dramatically decreased with the addition of DIN and DIP [42]. In other words, high nutrient input could reduce beta diversity driven by phytoplankton competition. In addition, we observed a scale-dependent distance-decay relationship between community similarity and geographic distance in picophytoplankton communities (Fig. 4), which is similar to observations in bacterial [43] and fungal communities [44]. However, our results were in agreement with a previous study [45] demonstrating that environmental variables were more important than geographic distance on a regional scale ( Fig. 4 and Table 1).

Co-occurrence networks reveal picophytoplankton community dynamics
We constructed co-occurrence networks to explore interspecies interactions that drove the picophytoplankton community responses to eutrophic disturbances that correlated with ecosystem stability (Fig. 6). In this study, the inter-domain co-occurrence network exhibited "small-world" features, non-randomly connected properties, and modular structure ( Fig. 6 & Table S6). The degree and betweenness centrality value described the level of connectedness between OTUs, which provides information about the importance of OTUs to network connectivity. OTUs act as important keystone species in this study that mainly belonged to the OE group, indicating that these taxa may play irreplaceable roles in maintaining picophytoplankton community structure ( Figure S1). Modularity may re ect synergistic relations, competitive interactions, and niche differentiation, inducing non-random patterns of interaction and the complexity of ecological networks [46]. The association of picophytoplankton competition and succession played a signi cant role in altering community structure when ecosystems suffered nutrient variations [47]. Most modules in our network exhibited distinct eutrophic variation, potentially identifying different groups that performed different functions [48].
Speci cally, OTUs which were mainly attached to Rhodophyta in module II were speci c to the HE. The majority of OTUs in module IV primarily occurred in the ME, while OTUs that were mostly a liated with Bacillariophyta in modules I and III were speci c to the OE (Fig. 6b & S2). Interestingly, the correlations associated with Chlorophyta were greater in these modules (Table S6). Furthermore, network variation was signi cantly correlated with Euclidean distance of nutrients, which may play an important role in the composition and connection of corresponding module. Evidence indicates that increasing eutrophication stress enhances the strength of interactions among species [49], whereas a densely competitive relationship may exist in the oligotrophic area. These results showed that the complexity of interspecies interaction varied with eutrophic levels because niche conservatism led to the coexistence of closely related picophytoplankton taxa. Therefore, the complexity of interspecies interactions was dependent on eutrophic states.

Stochastic and deterministic processes shape the picophytoplankton community
It is now generally accepted that deterministic and stochastic processes jointly control the diversity of microbial communities, but uncertainty remains about their relative importance in disturbance intensity [50]. In this study, deterministic processes of environmental selection play a more important role than stochastic processes in shaping picophytoplankton community assembly at the regional scale of heavy eutrophication (HE), while stochastic processes such as dispersal and drift play a predominant role in structuring picophytoplankton communities on large scales (Fig. 7). This divergence was attributed to the notion that picophytoplankton assemblages are resilient to short-term eutrophication disturbances, whereas continuous disturbances induced environmental selection pressure that become dominant with increasing eutrophic states [51]. Accordingly, although stochastic processes were therefore dominant in governing picophytoplankton community assembly in moderate, oligo eutrophication areas and the entire Beibu Gulf, the underlying deterministic processes depend on the severity and duration of eutrophic disturbances [52]. By quantifying the effects of an individual environmental factor, we found that nitrogen concentration (DIN and TDN) was the most important factor associated with environmental selection, especially in the HE ( Table 2). Although nutrients exhibited strong selective pressure on the picophytoplankton community, it should be noted that NTI decreased with increasing DIN and EI, indicating that picophytoplankton community assembly was less phylogenetically clustered in increasing eutrophic states (Fig. 8). These ndings were consistent with many studies demonstrating that environmental ltering is a potentially essential driver of variation in picophytoplankton community structure at small spatial scales, whereas dispersal limitation was not strong [53]. In summary, both deterministic and stochastic processes contribute to picophytoplankton community structure and assembly in the Beibu Gulf. Taken together, our study provides valuable insights into how community assembly processes shape biogeographic patterns of picophytoplankton communities in the Beibu Gulf.

Conclusion
Picophytoplankton community is considered to be one of the indicators of environmental change and ecosystem stability owing to their fast and strong responses to environmental disturbances such as eutrophication in coastal ecosystems. In this study, we found that picophytoplankton communities exhibited different and complex compositions, diversity and assembly processes in response to environmental changes under eutrophic states. The phyla/class of Rhodophyta, Bangiophyceae, and Coscinodiscophyceae were signi cantly over-dominated in the HE related to heavy eutrophic state, whereas the ME and OE groups related to medium and oligo eutrophic states, respectively, were mostly a liated with Bacillariophyta, Coscinodiscophyceae, and Bacillariophyceae. Both stochastic and deterministic processes simultaneously in uenced the assembly of picophytoplankton communities along eutrophic states. Specially, stochastic processes dominantly govern governing picophytoplankton community assembly, whereas the underlying deterministic processes depend on the severity and duration of eutrophic disturbances. In addition, co-occurrence network analyses revealed that the interspecies interaction of picophytoplankton community varied based on different eutrophic states. Our results provide strong evidence for understanding the general pattern of picophytoplankton community diversity, assembly and interaction in response to increasing eutrophication.

Sampling Sites and Environmental Parameters
A total of 90 water samples were collected across three prede ned sampling areas encompassing Maowei Sea estuary (H1-H4 sites), adjoining Qinzhou Bay (M1-M6 sites) and Beibu Gulf open sea (O1-O8 sites), with 5 duplicates per site. All samples were collected within two weeks to negate the effects of seasonality on picophytoplankton community during October, 2018. Each sample was collected at a depth 5 m using a sterilized water-sampling device and immediately transferred to the laboratory within 24 h for further processing. The 10 L seawater of each sample was pre ltered using 10-µm pore-size polycarbonate, followed by 2-µm pore-size lters, and subsequently ltered through 0.2-µm pore-sized polycarbonate membranes to obtain the < 2-µm ('pico') fraction (Millipore Corporation, Billerica, MA) using a vacuum pump. Then, the lters were immediately stored in sterile 50 mL centrifuge tubes at -20 °C until DNA extraction. At the time of sampling in each station, a portable meter (556 MPS; YSI, USA) was used to measure seawater salinity, pH and temperature. Additionally, other environmental parameters including dissolved inorganic nutrients (NO 2 − , NO 3 − , NH 4 + , and PO 4 3− ) measured by spectrophotometric and colorimetric methods [54], the total organic carbon (TOC) and chemical oxygen demand (COD) measured by the standard method [55], the total dissolved nitrogen (TDN) and total dissolved phosphorus (TDP) determined by the Cu-Cd column reduction and spectrophotometric phosphomolybdate blue method, respectively [56], and Chl a measured by a uorescence spectrophotometer [57] were recorded. Dissolved inorganic nitrogen (DIN) and dissolved inorganic phosphorus (DIP) were composed of NO 2 − , NO 3 − , NH 4 + and PO 4 3− , respectively. The eutrophication index (EI) was calculated as: EI = DIN × DIP × COD × 10 5 /4500 [58]. Eutrophic states were categorized as follows: oligo eutrophication ranged from 1.0 < EI < 3.0; medium eutrophication ranged from 3.0 < EI < 9.0; and heavy eutrophication ranged from EI > 9.0 [58]. We determined the that trait distribution of eutrophication showed a decreasing trend from the Maowei Sea to Beibu Gulf open sea. Consequently, sampling sites were divided into three groups: heavy eutrophication (HE), medium eutrophication (ME) and oligo eutrophication (OE).

DNA Extraction, PCR Ampli cation and Sequencing
Environmental DNA extraction was performed using a DNeasy Power Water Kit (TIANGEN, China) through 0.2-µm pore size polycarbonate membranes according to the manufacturer's protocols. DNA yield and purity were evaluated using a Nanodrop-2000 Spectrophotometer (Thermo Scienti c, USA). The DNA sample was preserved at − 80 °C. Form ID rbcL gene fragment (550 bp) was ampli ed from environmental DNA for all the stations using previously published rbcL primers (F:5′-GATGATGARAAYATTAACTC-3′ and R:5′-ATTTGDCCACAGTGDATACCA-3′) [59].

Statistical Analysis
Group differences among the various environmental variables and α-diversity index were tested with a permutation test (number of permutations = 10,000) using the One-way analysis of variance (ANOVA) in SPSS Statistics v24.0 [60]. The α-diversity of each sample between different groups was determined by calculating Shannon, Simpson, Chao1 and Observed numbers index and calculated using 'vegan' R package [61]. Spearman's rank correlation between the α-diversity index and environmental factors was calculated using 'psych' R package [62] with P-value adjusted by the FDR correction. The community composition was analyzed to investigate the variation of relative abundance and visualized using the 'pheatmap' R package [63]. To investigate successional patterns of picophytoplankton communities between groups, the β-diversity was undertaken using Bray-Curtis dissimilarity index and visualized by Nonmetric multidimensional scaling (NMDS). The analysis of similarities (ANOSIM) was performed to analysis similarities variation in community using the 'vegan' R package. Correlations between picophytoplankton community and constraining variables were tested using Mantel tests between Bray-Curtis distance matrices of community composition and Euclidean distance matrices of environmental factors. Linear regressions were performed to determine the relationships between the Bray-Curtis similarity of picophytoplankton community and geographical distance, environmental factors, respectively.

Co-occurrence network analysis
Only OTUs in at least ve samples with more than 20 sequences in picophytoplankton community were selected to construct co-occurrence networks, as previous study [21]. To identify correlation between OTUs with trait distributions, spearman's rank correlations were calculated as being positively or negatively correlated with a certain trait using "Hmisc" R package [64]. Only robust (|r| > 0.6) and statistically signi cant (FDR-adjusted P-value < 0.05) correlations were selected for network analysis. A GML format network le was generated using the "igraph" R package and visualized using Gephi 0.9.2 [65,66].. Networks are composed of nodes, which correspond to individual taxon, connected by links (or edges), that represent signi cant co-occurrence relationships. Nodes with high degree (> 65) are recognized as potential keystone species of the state in co-occurrence networks. The network-level (average node degree, clustering coe cient, average path length, and modularity) topological features of a network was calculated. Meanwhile, the real networks were compared with 1000 Erdös-Réyni random networks, which had the identical number of nodes and edges as the real networks, were generated in the "igraph" R package [67]. We also tested whether environmental factors related to its modularization in the network using correlation tests in R.

Analysis of picophytoplankton community assembly processes
To assess the dominance of stochastic and deterministic assembly processes along the eutrophic gradients, we compared all possible pairwise comparisons of β-nearest taxon index (βNTI) within certain group and chemical properties ranges with 999 randomizations in the "picante" R package. If βNTI values are > 2 or < -2, the dominance of deterministic processes in shaping the community composition, whereas βNTI values are between − 2 and 2, community assembly was dominated by stochastic processes. To further evaluate the variation in community assembly processes along nutrient gradients, nearest-taxon index (NTI) and βNTI values were regressed against Euclidean distance matrices of nutrition factors. We        Co-occurrence patterns of picophytoplankton communities at the family level. Node colors indicate the different modularity classes (a) and different taxa (b). A connection represents a strong (Spearman's |p| > 0.60) and signi cant (FDR-adjusted P < 0.05) correlation between picophytoplankton community in modules and environmental factors. The size of each node is proportional to the eigenvector centrality.  The percent of turnover in picophytoplankton community assembly. Picophytoplankton community assembly was governed primarily by various deterministic, including homogeneous and variable selection, and stochastic processes, including dispersal limitations, homogenizing dispersal, and ecology draft, in all groups and three groups.

Figure 7
The percent of turnover in picophytoplankton community assembly. Picophytoplankton community assembly was governed primarily by various deterministic, including homogeneous and variable selection, and stochastic processes, including dispersal limitations, homogenizing dispersal, and ecology draft, in all groups and three groups.

Figure 8
The relationship among NTI, DIN, and EI for picophytoplankton community (a). The relationship between βNTI and change in DIN and EI which was calculated by Euclidean distance (b). Linear regressions models and associated correlation coe cients are provided on each panel. Horizontal dashed black lines indicate upper and lower signi cance thresholds at βNTI = +2 and −2, respectively. Boxplots within the panels showed the variation in NTI and βNTI in three groups. Boxplots that do not share a letter are signi cantly different (p < 0.05; multiple comparisons with Kruskal-Wallis).

Figure 8
The relationship among NTI, DIN, and EI for picophytoplankton community (a). The relationship between βNTI and change in DIN and EI which was calculated by Euclidean distance (b). Linear regressions models and associated correlation coe cients are provided on each panel. Horizontal dashed black lines indicate upper and lower signi cance thresholds at βNTI = +2 and −2, respectively. Boxplots within the panels showed the variation in NTI and βNTI in three groups. Boxplots that do not share a letter are signi cantly different (p < 0.05; multiple comparisons with Kruskal-Wallis).