2.1 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 significantly between the stations sampled from various groups (one way-ANOVA analysis, p < 0.001) (Table S2). The DIN was comparable between the Maowei Sea (0.564 ± 0.091 mg/L) and Qinzhou Bay (0.160 ± 0.057 mg/L). In addition, the EI dropped sharply (range from 0.55 to 25.48) and tended to decrease from the estuary to the open sea area. The lowest value was noted at the Beibu Gulf open sea (1.189 ± 0.322), and the highest value was observed at the Maowei Sea (16.380 ± 5.152) (Table S1 & 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 reflected 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 significantly 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 significant 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 significantly 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 NH4+ (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 significant correlation (r = 0.58, p < 0.001) (Table S4).
2.3 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 significantly over-dominated in the class of Bangiophyceae, and Coscinodiscophyceae, whereas the medium and oligo eutrophication groups were mostly affiliated 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 significant 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 significant 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, R2 = 0.412, p < 0.001) (Fig. 4a). DIP and TDP (regression, R2 = 0.594, p < 0.001) also had significant 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 significantly 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 significantly 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 significantly different as demonstrated nutrient gradient distribution patterns. For example, assemblages affiliated with Bacillariophyceae (including the Achnanthidiaceae, Mastogloiaceae, Amphipleuraceae, and Naviculaceae) were significantly negatively correlated with eutrophic states, whereas assemblages belonging to Coscinodiscophyceae (Coscinodiscaceae and Lauderiaceae) were significantly 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.
2.4 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 identified as keystone species, and these species were significantly 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, co-occurrence 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 specific to a particular eutrophic group (Figure S2). For example, modules I and III were specific to the OE group, whereas modules II and IV were specific to the HE group and ME group. Furthermore, the observed modularity, clustering coefficient 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 significantly 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).
2.5 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 (nearest-taxon 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 influenced 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 significantly 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 significantly positively correlated to changes in DIN (r = 0.665, p < 0.001), NO3− (r = 0.635, p < 0.001), salinity (r = 0.520, p < 0.001), EI (r = 0.520, p < 0.001), NH4+ (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 significantly (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 significantly affected by nitrogen concentration, salinity, and temperature.
Table 1
Mantel test for β-diversity. The mantel test (Pearson correlations) showed the correlations between β-diversity of picophytoplankton communities (Bray-Curtis dissimilarity) and environmental factors (Euclidean distance) with 999 permutations. (*: p < 0.05, **: p < 0.01, ***: p < 0.001)
Variables
|
All
|
High
|
Medium
|
Oligo
|
Geo-distance
|
0.488***
|
0.629***
|
0.413***
|
0.397***
|
Env-factors
|
0.603***
|
0.622***
|
0.629***
|
0.475***
|
Temp
|
0.447***
|
0.133*
|
0.149**
|
0.038
|
pH
|
0.435***
|
0.187*
|
0.369***
|
0.170**
|
Salinity
|
0.403***
|
0.619***
|
0.324***
|
0.359***
|
NO2−
|
0.509***
|
0.770***
|
0.260***
|
0.542***
|
NO3−
|
0.437***
|
0.690***
|
0.450***
|
0.288***
|
NH4+
|
0.324***
|
0.552***
|
0.400***
|
0.608***
|
Chl a
|
0.261***
|
0.581***
|
0.302***
|
0.400***
|
DIN
|
0.370***
|
0.568***
|
0.363***
|
0.269***
|
TDN
|
0.288***
|
0.391***
|
0.250***
|
0.030
|
DIP
|
0.520***
|
0.759***
|
0.417***
|
0.304***
|
TDP
|
0.520***
|
0.760***
|
0.431***
|
0.390***
|
TOC
|
0.132***
|
0.458***
|
0.372***
|
0.258***
|
COD
|
0.166***
|
0.548***
|
0.450***
|
0.199**
|
EI
|
0.371***
|
0.368**
|
0.552***
|
0.492***
|
Table 2
Mantel test for βNTI. The mantel test (Pearson correlations) showed the correlations between βNTI and environmental factors (Euclidean distance) with 999 permutations. (*: p < 0.05, **: p < 0.01, ***: p < 0.001)
Variables
|
All
|
High
|
Medium
|
Oligo
|
Geo-distance
|
0.041
|
0.594***
|
0.079
|
-0.070
|
Env-factors
|
0.516***
|
0.952***
|
0.332***
|
-0.001***
|
Temp
|
0.300***
|
-0.016
|
-0.128
|
0.098
|
pH
|
-0.221
|
0.658***
|
0.570***
|
0.044
|
Salinity
|
0.520***
|
0.378**
|
0.530***
|
0.016
|
NO2−
|
0.068*
|
0.166
|
0.544***
|
0.043
|
NO3−
|
0.635***
|
0.249*
|
0.135*
|
-0.095
|
NH4+
|
0.312***
|
0.867***
|
0.323***
|
0.016
|
Chl a
|
0.053
|
0.656***
|
0.426***
|
-0.011
|
DIN
|
0.665***
|
0.721***
|
0.490***
|
-0.112
|
TDN
|
-0.008
|
0.971***
|
0.530***
|
0.244***
|
DIP
|
-0.164
|
0.126
|
0.475***
|
0.057
|
TDP
|
-0.172
|
0.083
|
0.396***
|
0.056
|
TOC
|
0.030
|
0.361**
|
0.314***
|
0.016
|
COD
|
0.028
|
0.556***
|
0.471***
|
-0.082
|
EI
|
0.469***
|
0.966***
|
-0.172
|
0.160**
|