Heterogeneity of interactions of microbial communities in regions of Taihu Lake with different nutrient loadings: A network analysis

To investigate the differences in the interactions of microbial communities in two regions in Taihu Lake with different nutrient loadings [Meiliang Bay (MLB) and Xukou Bay (XKB)], water samples were collected and both intra- and inter-kingdom microbial community interactions were examined with network analysis. It is demonstrated that all of the bacterioplankton, microeukaryotes and inter-kingdom communities networks in Taihu Lake were non-random. For the networks of bacterioplankton and inter-kingdom community in XKB, higher clustering coefficient and average degree but lower average path length indexes were observed, indicating the nodes in XKB were more clustered and closely connected with plenty edges than those of MLB. The bacterioplankton and inter-kingdom networks were considerably larger and more complex with more module hubs and connectors in XKB compared with those of MLB, whereas the microeukaryotes networks were comparable and had no module hubs or connectors in the two lake zones. The phyla of Acidobacteria, Cyanobacteria and Planctomycetes maintained greater cooperation with other phyla in XKB, rather than competition. The relationships between microbial communities and environmental factors in MLB were weaker. Compared with the microbial community networks of XKB, less modules in networks of MLB were significantly correlated with total phosphorous and total nitrogen.

. Topological properties of the empirical species-species networks of microbial communities in lake zones with different nutrient loading levels and an associated random network. Random networks were generated by rewiring all of the links with the same numbers of nodes and edges to the corresponding empirical network. The numbers in parentheses indicate the standard deviation (SD) of topological properties of 1000 random networks. MLB, lake zone with high nutrient loading; XKB, lake zone with low nutrient loading. a Significant difference (P < 0.001) between the empirical network and the random network (Z-test). b Significant difference (P < 0.001) between network indexes for the two lake zones (Student t-test). c Numbers in parenthesis represent the number of microeukaryotes-microeukaryotes edge in the inter-kingdom network.
SCiEntifiC REPORts | (2018) 8:8890 | DOI: 10.1038/s41598-018-27172-z For bacterioplankton, microeukaryotes and inter-kingdom communities, the networks in MLB and XKB had higher CC than those in random networks, indicating that there are more highly interconnected (clustered) nodes in the observed networks than in random networks. A small-world network means that most nodes can be reached from every other node by a small number of hops or steps 38 . The clustering nodes in all networks (Fig. 1,  Supplemental Fig. S2) and the remarkably high average degree (AD) ( Table 1) suggested that these networks have 'small-world' properties 38,39 , especially for the microbial community networks in XKB. Although identical thresholds were used to define the networks, the sizes of the networks in the two regions were different (Table 1). The network properties for bacterioplankton in MLB were closer to those of random networks (Table 1). Our previous study found that stochastic processes played non-negligible roles in controlling the assembly of both bacterioplankton and microeukaryotes communities in MLB 14 . The higher nutrient loadings in MLB may account for the similarity to random networks in this region.
As shown in Figs 1 and S2, the network in XKB was significantly larger and more complex than that in MLB for bacterioplankton and inter-kingdom, but the networks in the two regions were comparable for microeukaryotes. Significant differences in all of the network indexes were observed between MLB and XKB for the bacterioplankton as well as the inter-kingdom networks (Student's t-test, P < 0.001) (Table 1). However, there were no significant differences in MD or CC (P > 0.05) between the microeukaryotes networks in the two regions. The greater complexity of the network in XKB (Fig. 1) may be due to the relatively lower nutrient loading level, which would lead to relatively stronger niche selection due to the greater competition for resources and less-diverse resources 6 . It has been suggested that bacterioplankton and microeukaryotes have similar cellular-mineral-environmental constraints 40 and it is well known that bacterioplankton and microeukaryotes are very sensitive to environmental variations 14,35 . Although the lower nutrient loading and stronger environmental filtering effects in XKB may lead to stronger niche selection, microeukaryotes are better able to adapt to environmental perturbations and compete with each other 40 , which contributes to the comparable microeukaryotes networks in MLB and XKB as well as the correlations between microeukaryotes and environmental variables (Supplemental Fig. S3).  (Table 1), and thus had modular structures 17,41,42 . Therefore, the species-species association networks were divided according to modules and obvious differences were found for the networks of bacterioplankton ( Fig. 1) and inter-kingdom (Supplemental Fig. S2) in MLB and XKB.

Co-occurrence/co
The network for XKB included more interactions (edges) ( Table 1), and the APL and ND of the network for bacterioplankton and inter-kingdom were significant lower in XKB than in MLB (Table 1). However, the average degree (AD), which is the most robust measure of network topology along with CC, of the network for bacterioplankton and inter-kingdom were significant higher for the networks in XKB 43 . In general, there were higher CC and AD but lower APL and ND indexes for the network of XKB, which indicated that the nodes in the network of XKB were more clustered and closely connected with plenty edges for each point, thus the network in XKB was more complex than that for MLB 43 . It is also observed that the interaction of inter-kingdom network was markedly stronger in XKB (195 positive correlations and 157 negative correlations) than those of MLB (39 positive correlations and 12 negative correlations) (Supplemental Table S1). However, the species-species association networks for microeukaryotes communities were not as obviously different as those for bacterioplankton between MLB (Fig. 1c) and XKB (Fig. 1d).
Environmental conditions would affect the interactions among microbial communities. It has been demonstrated that species that share similar ecological niches may exhibit competition when resources are scarce, but may show positive interactions under resource-rich conditions 3,44 . In ecological systems, coexistence is supported by niche processes like environmental filtering [45][46][47] . Species pairs that co-occur may share similar ecological characteristics 8,9,48 , which can be used to infer life-history strategies 49,50 . Glöckner et al. 51 found that bacterioplankton and microeukaryotes showed significant differences in the abundance and relationships among phyla under different nutrient loading conditions. The favorable nutrient state in MLB could weaken niche selection by providing more resources and reducing competition among species, which may explain the weak correlation and simple network 28,52 in MLB. However, the microeukaryotes networks in the two regions were comparable. It has been demonstrated that, in pelagic ecosystems, microeukaryotes are the primary consumers of phytoplankton, heterotrophic bacteria and archaea 34,53 . In addition, they serve as important trophic links for transferring carbon between the microbial food web and the metazoan food web 34 . Furthermore, uncommonly delineated microeukaryotes species may contribute to the difference between networks for bacterioplankton, inter-kingdom community and microeukaryotes in the two lake zones 14,54 . Intra-and inter-phyla co-occurrence/co-exclusion patterns for bacterioplankton. The O/R ratios (O: observed incidence of the co-occurrence of two taxa; R: random incidence of the co-occurrence of two taxa) for bacterioplankton were calculated to determine if OTUs from the same phylum/different phyla tended to exhibit co-occurrence or co-exclusion. Few significant O/R ratios in microeukaryotes networks were observed and the overwhelming majority of significant O/R ratios in inter-kingdom network were the co-occurrence/co-exclusion patterns for bacterioplankton, thus neither of them was shown. An O/R value > 1 means that the observed incidence of co-occurrence (O) of two taxa was higher than that expected at random. As shown in Supplemental Table S2, the O/R values were almost all significantly higher than 1 in both XKB and MLB for positive interaction, indicating a very strong co-occurrence pattern for intra-phylum OTUs, especially in MLB (O/R values almost all > 3). However, the O/R values for negative intra-phylum interaction were almost 0 or NA, and not significant in both XKB and MLB. Co-occurrence reflects commonly preferred conditions or cooperative behaviors 50 ; the higher nutrient loading in MLB provided more suitable environmental conditions for these phyla to live and co-occur. It has been demonstrated that strong ecological Only correlations between species that were statistically significant (P < 0.01, Q-value < 0.05) and strong (r ≥ 0.9) were shown. Red solid line means positive correlation and black lines mean negative correlation. Different microbial phyla were represented with different colors and the number on each node means the number of OTUs clustered at 97% similarity. The circles consist of some nodes mean modules. The O/R values for inter-phyla interactions in XKB were more significant than those in MLB for both positive and negative interactions. For example, Chloroflexi showed significant and strong co-occurrence relationships (O/R > 3, P < 0.01) with many phyla (Acidobacteria, Actinobacteria, Cyanobacteria, Gemmatimonadetes, Planctomycetes, Betaproteobacteria and unclassified) in the bacterioplankton community of XKB, but showed no significant co-occurrence patterns in MLB. Other phyla, such as Acidobacteria (with Actinobacteria, Chloroflexi, Cyanobacteria, Gemmatimonadetes and Planctomycetes), Cyanobacteria (with Acidobacteria, Bacteroidetes, Betaproteobacteria, Chloroflexi, Gammaproteobacteria, Gemmatimonadetes, Planctomycetes and unclassified) and Planctomycetes (with Acidobacteria Actinobacteria, Bacteroidetes, Betaproteobacteria, Chloroflexi, Cyanobacteria, Firmicutes and Gemmatimonadetes and unclassified), showed a similar difference in co-occurrence patterns as Chloroflexi in the two lake zones.
The phyla of Acidobacteria, Cyanobacteria and Planctomycetes maintained strong relationships with other phyla in XKB (Supplemental Table S2), which may reflect greater cooperation with other phyla, rather than competition. Therefore, once in a suitable environment, bacterioplankton will live, and intra-phyla co-occurrence will usually occur. On the other hand, an unsuitable environment would dramatically increase inter-phyla relationships. Thus, there were more relationships, such as symbiosis or competition, among species in XKB than in MLB, indicating that the networks between bacterioplankton may change with different nutrient loading levels.
Detection of the topological roles of nodes for bacterioplankton and inter-kingdom network. The topological roles of the OTUs identified in these four networks are shown as a Zi-Pi plot (Fig. 2).
The detection of the topological roles for microeukaryotes was omitted because there was neither module hub nor connector in the microeukaryotes network. It was observed that most of the OTUs (98.8% and 97.7% for MLB and XKB, bacterioplankton, respectively; 100% and 99.38% for MLB and XKB, inter-kingdom, respectively) were peripherals, with most of their links inside their modules. Furthermore, among these peripherals, most OTUs (91.7% and 67.8% for MLB and XKB, bacterioplankton, respectively; 92.71% and 82.08% for MLB and XKB, inter-kingdom, respectively) had no links outside their own modules (Pi = 0). For the network of bacterioplankton, 4 (1.2%) OTUs were identified as module hubs network of MLB, 3 (1.4%) and 5 (0.9%) OTUs were identified as module hubs and connectors in the network of XKB respectively ( Fig. 2a and Table 2). For the network of inter-kingdom, 2 (0.41%) and 1 (0.21%) OTUs were identified as module hubs and connectors in the network of XKB respectively ( Fig. 2c and Table 2), no connectors or module hubs were observed in the network of MLB for inter-kingdom ( Fig. 2c and Table 2). It is worth mentioning that OTUE00117 (phylum-Archaeplastida) was identified as module hubs in the network of XKB for inter-kingdom (Table 2), showing the important component of microeukaryotes in the whole ecosystem 35,36 . It is reported that the Archaeplastida are a major group of microeukaryotes, which may be the reason why it was the module hubs of the network 57 . Furthermore, no network hubs were observed in all of the networks in two lake zones (Fig. 2). Table 2 shows that all of the module hubs were from different modules (BM1, BM2, BM3 and BM6 for bacterioplankton; X3, X4, X12 for inter-kingdom) in the network of MLB for bacterioplankton and both lake zones for inter-kingdom. However, for bacterioplankton in XKB, module hubs were mostly from module BX3 while the connectors were mostly derived from module BX2 and were classified into different phyla (Planctomycetes, Actinobacteria and unclassified). The other two connectors from modules BX7 and BX4 were classified as Gammaproteobacteria and Cyanobacteria (Table 2). Furthermore, Table 2 shows that all of the module hubs and connectors were from different phyla in the specified region (Bacteroidetes, Alphaproteobacteria, Actinobacteria and Gemmatimonadetes for MLB, Actinobacteria, Alphaproteobacteria and Betaproteobacteria for XKB) and Alphaproteobacteria and Actinobacteria appeared to be module hubs in both MLB and XKB. Members of the Actinobacteria and Alphaproteobacteria were widely present in both MLB and XKB (Fig. 1). Hub OTUs in the co-occurrence network mostly belonged to Actinobacteria and Alphaproteobacteria ( Table 2), suggesting that some OTUs of Actinobacteria and Alphaproteobacteria play important roles in the networks in the two regions. Associations among bacterioplankton are usually established by a cluster of multiple highly interacting species with similar ecological niches or cooperation 28 . Actinobacteria play an important role in the decomposition of organic materials and the production of secondary metabolites with very diverse activities 28,58 , which may result in their influential status. It has been reported that Alphaproteobacteria isolates can either promote or inhibit the growth of coexisting blooming Cyanobacteria in freshwater lakes, which implies a strong functional interaction 3 . In other studies, Alphaproteobacteria were prominent members of modules at all time points and co-occurred with Actinobacteria and other phyla 59 , which may explain the module-hub role of Alphaproteobacteria. Module hubs in the MLB bacterioplankton community network ensured that species within a module were linked more tightly, which explains why the MLB network had much higher modularity than XKB (Table 1). Similarly, connectors in the XKB network made much tighter interactions among modules than those in MLB, which also confirmed that the network structures in the two regions were different.

Relationships between modules and environmental factors.
To explore the effects of environmental factors on species-species association networks, environmental factors were added to these networks (Supplemental Figs S3 and S4, in which red lines signify positive correlations and black lines signify negative correlations). Compared with the environmental variables in MLB, those in XKB had greater effects on these modules (Supplemental Figs S3 and S4). Most nodes in modules BX1, BX2 and BX4 for bacterioplankton (Supplemental Fig. S3) and modules X1 and X2 for inter-kingdom (Supplemental Fig. S4) network were positively correlated with total phosphorous (TP) and total nitrogen (TN) (Supplemental Table S3).
An eigengene analysis 60 was performed for modules that were positively or negatively associated with environmental variables to quantitatively describe the relationships between modules and environmental variables. The eigengene network analysis was omitted for microeukaryotes and inter-kingdom since the module of microeukaryotes and inter-kingdom (MLB) had only a very few significant correlations with environmental variables (Supplemental Table S3). The results of the eigengene analysis on the modules of bacterioplankton networks are shown in Supplemental Fig. S5. Figure 3 shows the relationships between eigengenes in the MLB and XKB networks and environmental variables for bacterioplankton, and represents the first evidence that environmental factors have positive or negative effects on particular bacterioplankton modules. In MLB, module BM7 had significant positive correlations with TN and TP, module BM8 had a significant positive correlation with NO 2 , module BM9 had a significant positive correlation with NH 4 + -N and a negative correlation with pH, modules BM11 and BM13 had significant negative correlations with TN and TP, and module BM14 had a significant positive correlation with dissolved organic carbon (DOC) (Fig. 3a). In XKB, modules BX1, BX2, BX4, BX5 and BX10 had significant positive correlations with environmental variables including TN, TP and pH, whereas modules BX6 and BX7 had significant negative correlations with these variables (Fig. 3b). The results regarding the correlations between the module eigengenes and environmental variables in the two regions for bacterioplankton and microeukaryotes networks are shown in Supplemental Table S3.
In our previous study, both nutrient variables (TN, TP) and pH significantly affected the compositions of both the bacterioplankton and microeukaryotes communities in XKB (P < 0.05), whereas environmental factors were not significantly related to the composition of the microbial communities in MLB, except for a weak correlation between DOC and the microeukaryotes community 14 . Strom 61 found that the functional traits of microorganisms are products of multiple populations within these communities rather than those of a single population. All of the results of our previous study are fairly consistent with those shown in Fig. 3, and proved that the modules are composed of bacterial clusters with similar ecological niches. The relationships between OTUs and ecological relatedness are complicated, and depend on individual microbial groups, as well as environmental conditions 28,61,62 . For instance, nutrient loading had a very strong effect on the modules in both MLB and XKB. Meanwhile, environmental factors affected each module differently (Fig. 3). The node composition is substantially different among different modules in the two regions, since they have different nutrient loadings. The different nutrient loadings suggest different extents of environment stress, which would influence the composition and turnover of the microbial communities 6 . Therefore, it is reasonable that the   relationships between the microbial communities and environmental factors in MLB are weaker than those in XKB, since niche selection in the former was weakened by reduced competition.

Conclusions
The results of the present study demonstrated that the network structure and co-occurrence patterns were significantly different between the MLB and XKB regions of Taihu Lake for bacterioplankton, microeukaryotes and inter-kingdom community. The properties of the obtained networks were significantly different from those of random networks, indicating that the assembly of microbial communities in these lake zones was non-random. The region with lower nutrient loading (XKB), and stronger environmental filtering effects, maintained a higher complexity for the whole network and more complex co-occurrence pattern compared with those in MLB for bacterioplankton and inter-kingdom community. It is also observed that the inter-kingdom interactions were stronger in XKB than those of MLB, whereas the networks for the two regions were comparable for microeukaryotes. Non-random co-occurrence of taxonomically related bacterioplankton was also observed, and OTUs from the same phylum tended to co-occur in both lake zones. The relationships between microbial communities and environmental factors in MLB were weaker than those in XKB for bacterioplankton, microeukaryotes and inter-kingdom community, and environmental factors affected each module differently. This study was limited in that it considered only two lake zones; a wider range of study areas will be needed to determine the impacts of environmental factors on the interactions among microbial communities under distinct nutrient conditions.

Materials and Methods
Sample collection, Illumina second-generation sequencing and analysis. Ten water samples were collected from two regions of Taihu Lake (Meiliang Bay and Xukou Bay, which represent different nutrient loadings) in October 2015, respectively (Supplemental Table S4). Physicochemical characteristics of the water samples (including the salinity, temperature, oxidation reduction potential (ORP), pH, and conductivity) were measured in situ using a calibrated multifunction water-quality sonde (YSI 6600, Yellow Springs, OH, USA -N) were measured in laboratory as described by Zhao et al. 14 . The sample used in the present study were the same as our previous study and the methods of DNA extraction, amplification, pyrosequencing and data analysis have been described 14 . The total OTU richness was 3247 and 2059 at 97% similarity cutoff for all the rarefied samples for the 16S rRNA and 18S rRNA, respectively. The raw reads were deposited into the NCBI Sequence Read Archive (SRA) database (Accession Number: SRP090623).
The computational procedures for network construction. All samples were divided into two groups [Meiliang Bay (MLB) and Xukou Bay (XKB)] representing higher (MLB) and lower (XKB) nutrient loading levels. To improve the network reliability, only OTUs that appeared in at least 8 samples in each group were considered 63 . The relative proportions of sequence numbers were used for the following correlation analysis, since the sequence numbers of individual OTUs significantly varied among the samples 58 . In each group, a correlation matrix was constructed according to the relative abundance of the OTUs in each sample. Both the correlation matrix (R matrix) and the significance matrix (P matrix) were calculated using the 'Hmsic' package in R by calculating all possible pairwise Spearman rank correlations between all OTUs 64 . Only robust (Spearman correlation coefficient ≥0.9 (or ≤−0.9)) and statistically significant (P < 0.01) correlations were considered 65 . The correlation approach was justified by the analysis for the sampling effectuated according to Weiss et al. 66 and then improved by Q-value as described below. The possibility of obtaining false results was reduced by calculating a Q-value, which represented the fraction of false positives or negatives if a given pair was identified to be significant (Q-values < 0.05), using the 'qvalue' package in R 7 . The node degree (i.e., number of edges connected to the node) was plotted against the probability P(k) that a node would have that degree in the network. Three methods (Power law, Exponential law and Truncated law) of power law-fitting of the degree distribution in networks in the two regions were applied 65 . The existence of meaningful, non-random associations in the networks of the two regions was demonstrated by the structural similarity among these ecological networks, in comparison to a Gaussian connectivity distribution predicted by an expectation of randomness.
Network characterization. The resulting correlation matrix was transformed into a Cytoscape dataset in R. Cytoscape v2.8.2 was then used 67 for network visualization and topological analysis. Other information regarding nodes (OTUs), taxonomy, module, edge, weights, and positive and negative correlations, was also imported into Cytoscape.
Each network was separated into modules using the fast greedy modularity optimization 68 . Various indexes, including modularity (MD), clustering coefficient (CC), average path length (APL), network diameter (ND), average degree (AD) and graph density (GD), were used to describe the properties and the overall topologies or structures of the networks. Most of these parameters were calculated using the 'igraph' packages in R 69 . Random networks were also generated using the 'igraph' packages in R. For each network in this study, 1000 random networks were generated, and all of the network indexes were calculated individually. The average and standard deviation for each index of all of the random networks were then obtained. The statistical Z-test was used to test the significance of differences between the indexes of the observed and random networks. Student's t test was used to test the significance of differences between the indexes of the networks in MLB and XKB for bacterioplankton, microeukaryotes and inter-kingdom communities.
Patterns of co-occurrence and co-exclusion. An R script was developed to check the observed (O) and random (R) incidences of the microbial patterns of co-occurrence and co-exclusion. The O/R ratio has been used as a benchmark for checking non-random assembly patterns in complex bacterial communities 17 . Here, we calculated the observed incidence of the co-occurrence (O) of two taxa as the relative percentage of the number SCiEntifiC REPORts | (2018) 8:8890 | DOI:10.1038/s41598-018-27172-z of observed edges between them, whereas the random incidence of co-occurrence (R) was calculated as the mean value of the observed incidence of co-occurrence for 1000 random networks. Hence, the degree of disagreement between the O and R incidences of co-occurrence may be used as a benchmark for exploring non-random assembly patterns in complex microbial communities 51 . In this study, the observed (O) and mean value of the random (R) incidences of co-occurrence and their significance levels were calculated according to Zhao et al. 3 . The R code for calculating the O/R has been attached in the supplementary information (Supplemental R code 1).

Topological roles of individual nodes. Visualization of the topological roles of individual nodes reveals
the effects of the nutrient loading level on key microbial populations. Topologically, different OTUs (nodes) play distinct roles in the network 70 . The topological roles of different OTUs can be described by two parameters: within-module connectivity (Zi), which describes how well a node is connected to other nodes within its own module, and connectivity among modules (Pi), which reflects how well a node connects to different modules 69,71 . Zi and Pi are calculated as described by Guimera and Amaral 71 . The R code for calculating the Pi and Zi has been attached in the supplementary information (Supplemental R code 2). According to the simplified classification used in networks 72 , the nodes in a network are divided into four subcategories: (i) peripheral nodes (Zi ≤ 2.5, Pi ≤ 0.62), which have low Z and P values (i.e., they have only a few links and almost always to species within their modules); (ii) connectors (Zi ≤ 2.5, Pi > 0.62), which have a low Z but a high P value (i.e., these nodes are highly linked to several modules); (iii) module hubs (Zi > 2.5, Pi ≤ 0.62), which have a high Z but a low P value (i.e., they are highly connected to many species in their own modules); and (iv) network hubs (Zi > 2.5, Pi > 0.62), which have high Z and P values (i.e., they act as both module hubs and connectors) 58,71 . Relationships between modules and environmental variables in the networks of bacterioplankton and microeukaryotes in two lake zones. To investigate the relationships between the distribution of nodes in networks and environmental variables, environmental variables were integrated into the networks. Only correlations between environmental variables and species that were statistically significant (P < 0.05) and strong (r ≥ 0.6 or r ≤ −0.6) were considered. In addition, to quantitatively describe the relationships between modules and environmental variables, an eigengene network analysis was performed. In this approach, each module was decomposed into a single representative abundance profile called the module eigengene. The molecular ecological network consisted of many nodes and edges. It was difficult to retrieve information intuitively, but the network could be simplified using various methods, such as module partitioning. Modules can be treated as single units for biologically motivated data reduction 73 . First, all of the nodes in module i were selected, and their eigengene values were calculated using the 'WGCNA' packages in R 74 . Second, the Spearman correlations were calculated between each eigengene and the environmental variables 3 . The calculations were performed as described by Zhao et al. 3 .