Seasonal Water Level Fluctuation and Concomitant Change of Nutrients Shift Microeukaryotic Communities in a Shallow Lake

: Seasonal water level fluctuations (WLFs) impose dramatic influences on lake ecosystems. The influences of WLFs have been well studied for many lake biotas but the microeukaryotic community remains one of the least-explored features. This study employed high-throughput 18S rRNA gene sequencing to investigate the spatiotemporal patterns of microeukaryotic communities in the dry and wet seasons with concomitant change of nutrients in Poyang Lake, which experiences huge seasonal WLFs. The results showed that the dry season and wet season had distinct microeukaryotic community compositions and structures. In the dry season, Ciliophora (13.86–40.98%) and Cryptomonas (3.69–18.64%) were the dominant taxa, and the relative abundance of these taxa were signiﬁcant higher in the dry season than wet season. Ochrophyta (6.88–45.67%) and Chlorophyta (6.31–22.10%) was the dominant taxa of microeukaryotic communities in the wet season. The seasonal variation of microeukaryotic communities was strongly correlated to seasonal nutrient variations. Microeukaryotic communities responded signiﬁcantly to dissolved organic carbon, total nitrogen, nitrate, and soluble reactive phosphorus in the dry season, and correlated to nitrate and total phosphorus in the wet season. The microeukaryotic community showed di ﬀ erent modular structures in two seasons, and nutrient variations were the key factors inﬂuencing seasonal variations of the modular structures. Moreover, microeukaryotic community networks based on di ﬀ erent seasons indicated that the microeukaryotic community co-occurrence patterns were not constant but varied largely associating with the nitrogen and phosphorus variations under the e ﬀ ects of WLFs. Our results are important for understanding how microeukaryotic communities respond to nutrient variation under seasonal water level ﬂuctuation.


Introduction
Hydrological regimes are key drivers of freshwater ecosystem structures and functions [1]. In response to hydrologic imbalance, water level fluctuations (WLFs) are a natural hydrological regime and an inherent feature of lake ecosystems, controlling ecological structures and processes in lacustrine ecosystems [2][3][4][5]. WLFs are more important for lakes, because even small WLFs can result in enormous WLFs have complex effects on lake ecosystems, and it is crucial to understand the dynamic of microeukaryotic communities experiencing excessive seasonal WLFs for maintaining the health of lake ecosystems under extreme climate. As the largest freshwater lake in China, Poyang Lake is a seasonal shallow lake experiencing dramatic WLFs, some low-lying areas be inundated in the whole year, and higher areas are inundated for a few weeks or months during the wet season [47,48]. The maximum inundation area in the wet season is more than 14 times of the minimum inundation in the dry season [49]. Here, using high-throughput 18S rRNA gene sequencing, we investigated the spatiotemporal patterns of microeukaryotic communities in the dry and wet seasons with special emphases on the effects of nutrients. The main objectives are to (1) compare the microeukaryotic communities in the dry and wet seasons; (2) assess the response patterns of microeukaryotic communities to nitrogen and phosphorus variations; and (3) uncover the co-occurrence patterns of microeukaryotic communities in different seasons. We hypothesized that WLFs and concomitant nutrient variations are the dominant drivers for community compositions and individual interaction of microeukaryotes. This study can enrich our knowledge of the influences of WLFs and coupled seasonal dynamics of lake ecosystems on microeukaryotic communities.

Study Area and Sample Collection
Poyang Lake, the largest freshwater lake in China, located in the lower reach of Yangtze River ( Figure 1). Poyang Lake is one of the two lakes retaining free connection with Yangtze River. As a shallow seasonal lake, Poyang Lake is a water-carrying and throughput lake and its water level was mainly controlled by Yangtze River and its five tributaries [50,51]. Poyang Lake watershed belongs to the climate region of East Asian Monsoon. The seasonal floods of its tributaries and Yangtze River lead to large seasonal fluctuations of the water level and surface area of Poyang Lake [6,52]. The seasonal variation curve of water level was unimodal ( Figure S1), in every year, the water level raised regularly from January to July, and then declined until December [45,53]. In the dry season, the lake flows from south to north and finally discharges to the Yangtze River. Many connected and disconnected sub-lakes are formed and the lake area shrinks to less than 1000 km 2 due to the drop of the water level. In the wet season, however, the floods from Yangtze River flow reversely to Poyang Lake [54], expanding the lake area over 4000 km 2 [53,55]. The spectacular water level fluctuation results in profound seasonal dynamics in biological and physicochemical properties of the lake ecosystem.
Water 2020, 12, x FOR PEER REVIEW 4 of 19 Agency, Washington, DC, USA, https://www.epa.gov/cwa-methods), total nitrogen (TN) (EPA 300.0), and total phosphorus (TP) (EPA 365.3) was analyzed using the ion chromatography method and the ascorbate acid colorimetric method after oxidation. Filtered water sample by glass fiber filters (Whatman, Maidstone, UK) was used to test nitrate (NO3 − ) with ion chromatography method (EPA 300.0), ammonium (NH4 + ) with indophenol colorimetric method (EPA 350.1), soluble reactive phosphorus (SRP) with ascorbate acid colorimetric method (EPA 365. 3), and dissolved organic carbon (DOC) with a Shimadzu TOC Analyzer (TOC-VCPH, Shimadzu Scientific Instruments, Columbia, Maryland). The information of the nutrients was summarized in Table S2.   During the dry season (10-14 January) and wet season (2-5 August) 2017, we collected 13 and 10 samples, respectively (GPS coordinates of the sampling sites were shown in Table S1). The water level was 16.5 + 0.1 m and 9.7 + 0.1 m during the sampling in the wet and dry seasons, respectively (measured at the Xingzi Hydrological Station, Figure 1). In general, the water level of Poyang Lake reaches the lowest in December or January and the highest in July or August, which causes some wetlands to be seasonally exposed. Due to the water surface narrowing and severe sand excavation in the northern part of the lake in the dry season, this area was not sampled during the dry season. In the wet season, we distributed the sampling sites across the whole lake, including the northern part of the lake. Moreover, the lake environments are more heterogeneous under low water level during the dry season, while more homogeneous with strong ecological drifts under frequent water exchange during the wet season. Thus, there were more sampling sites in the dry season than wet season. Although the sample sites did not coincide exactly in the dry and wet seasons, these samples were representative to explore the influences of extreme water level fluctuations in Poyang Lake within a year. In order to avoid the influences of storms on wet season sampling, the samples were collected before the storm. There was also no rainfall before and during dry season sampling.
Water samples at the depth of 0.5 m were collected at each sampling site. Sterilized polycarbonate membrane filter (0.2-µm, Whatman, Maidstone, UK) was used to filter 200-mL subsample. The filter was then stored in a sterilized centrifuge tube and immediately frozen in liquid nitrogen in the field. In the lab, filter samples were stored at −80 • C until DNA extraction. Another 500-mL subsample was collected for chemical analyses and transported to the laboratory at 4 • C. According to the Clean Water Act Analytical Methods (United States Environmental Protection Agency, Washington, DC, USA, https://www.epa.gov/cwa-methods), total nitrogen (TN) (EPA 300.0), and total phosphorus (TP) (EPA 365.3) was analyzed using the ion chromatography method and the ascorbate acid colorimetric method after oxidation. Filtered water sample by glass fiber filters (Whatman, Maidstone, UK) was used to test nitrate (NO 3 − ) with ion chromatography method (EPA 300.0), ammonium (NH 4 + ) with indophenol colorimetric method (EPA 350.1), soluble reactive phosphorus (SRP) with ascorbate acid colorimetric method (EPA 365.3), and dissolved organic carbon (DOC) with a Shimadzu TOC Analyzer (TOC-VCPH, Shimadzu Scientific Instruments, Columbia, Maryland). The information of the nutrients was summarized in Table S2.

DNA Extraction, PCR, and Sequencing
DNA was extracted from the filter samples using the TIANGEN-DP336 DNA Isolation Kit (TIANGEN-Biotech, Beijing, China) according to the instruction. DNA extracts were quantified using a Qubit 3.0 Fluorometer (Life Technologies, Darmstadt, Germany). A total of 50-100 ng DNA was used to generate amplicons. The 18S rDNA hypervariable regions of V9 were amplified using forward primers 1380F-CCCTGCCHTTTGTACACAC and reverse primer 1510R-CCTTCYGCAGGTTCACCTAC [56]. The forward primers were barcoded, and the barcodes were designed with considering the balanced guanine-cytosine (GC) content, minimal homopolymer runs, and no self-complementarity of more than two bases to reduce internal hairpin propensity [57]. The barcodes were prepended to DNA libraries for identification of individuals in pooled biomolecule populations. Each plankton DNA sample was PCR-amplified in triplicate. The 30 µL PCR mixture contained 15 µL of Phusion High-Fidelity PCR Master Mix (New Engliand Biolabs, Beverly, MA, USA), 0.2 µM of each primer, and 10 ng of sample DNA. The PCR reaction was conducted on a thermal cycler (ABI GeneAmp ® 9700, Foster City, CA, USA) using the following program: 3 min initial denaturation at 94 • C, 30 s denaturation at 94 • C with 24 cycles, 90 s annealing at 57 • C, 10 s extension at 72 • C, and finally 10 min extension at 72 • C. Amplified DNA were verified on 2% agarose gels (biowest agarose, Spain). The triplicate PCR products for each sample were pooled in equal quantity for purification using the Gel Extraction Kit (Qiagen, Hilden, Germany) and quantified using a Qubit 3.0 Fluorometer (Life Technologies, Germany). According to the manufacturer's instructions, DNA libraries were multiplexed and loaded on an Illumina MiSeq platform (Illumina, San Diego, CA, USA) at GENEWIZ, Inc. (Suzhou, China).

Analyses
Raw sequences were analyzed using the software package QIIME (version 1.9.1) [58]. The forward and reverse reads were joined and assigned to samples based on barcode and truncated by cutting off the barcode and primer sequence. Joined sequences were quality filtered based on following criteria: sequence length < 200 bp, no ambiguous bases, and mean quality score ≥ 20. The chimeric sequences were detected using the UCHIME algorithm and removed before downstream analyses. The effective sequences were clustered into operational taxonomic units (OTUs) at a threshold of 97% similarity level against the Silva 132 database [59]. The raw sequence data was available at the National Center for Biotechnology Information (PRJNA436872).
Relative abundance of OTUs was used to define abundant and rare taxa based on some thresholds [30,32]. Following recent studies [31,60], the OTUs were classified into four categories. Abundant taxa (AT) were defined as the OTUs with relative abundance ≥ 1% in all samples and the OTUs with relative abundance ≥ 1% in at least one sample and ≥ 0.01% in all samples. Rare taxa (RT) were defined as the OTUs with relative abundance < 0.01% in all samples and the OTUs with relative abundance < 0.01% in at least one sample and < 1% in all samples. Moderate taxa (MT) were defined as the OTUs with relative abundance between 0.01% and 1% in all samples. Conditionally rare and abundant taxa (CRAT) were defined as the OTUs with relative abundance < 0.01% and ≥ 1% in at least one sample.
Seasonal differences of microeukaryotic communities were assessed using principal coordinates analysis (PCoA) as well as analysis of variance using distance matrices (ADONIS) and multi-response permutation procedure analysis (MRPP) based on Bray-Curtis distance (Vegan package 2.5-3 in R). Moreover, to further explain the patterns of beta-diversity [61], we calculated Levin's niche width index [62] for each OTUs (spaa package 0.2.1 in R). Redundancy analyses (RDA) were generated using Vegan package (version 2.5-3) in R to reveal the association of the microeukaryotic communities with nutrient factors. The goodness of fit for each nutrient factor was estimated by applying the envfit function with 999 permutations. Variance partitioning analysis (VPA) was performed using Vegan package (version 2.5-3) in R to determine the relative contributions of C-factor (DOC), N-factor (TN, NO 3 − , and NH 4 + ), P-factor (TP and SRP), and the interactions of these factors (C × N, C × P, N × P, and C × N × P) [4,63]. t-tests were used to compare the differences of various parameters between different seasons. Network analyses were conducted to reveal co-occurrence patterns of the microeukaryotic communities [64,65]. To reduce the complexity of the data sets, OTUs with more than five sequences were retained for the construction of networks. All pairwise Spearman's rank correlations between those OTUs were calculated. Only robust (R > 0.8 or R < −0.8) and statistically significant correlations (p-values < 0.05, p-values were adjusted using FDR method) were considered [66]. R package igraph 1.2.4 [67] was used for network visualization, topological, and node/edge metrics (numbers of nodes, average clustering coefficient, network heterogeneity, and average number of neighbors), and modular analysis. Nodes represent the OTUs and edges represent Spearman's correlations relationships. The clustering coefficient quantified nodes clustering in the networks. Network heterogeneity quantified the diversity of connections between nodes in networks even with different topologies [68]. The average number of neighbors was calculated to represent the average connectivity of a node in the network [69]. The modular structure was analyzed and the modularity value > 0.4 suggests that the network is modular [70]. Mantel tests were used to analyze the relationship between environmental variables and modules in different seasons (PAST 3.18). The topological roles of each OTUs were determined according to the values of within-module connectivity (Zi) and among-module connectivity (Pi) [71]. Module hubs (i.e., highly linked species within their own module) with the value of Zi > 2.5 and Pi < 0.62, and connectors are linking different modules in a co-occurrence network have the connectivity value of Zi < 2.5 and Pi > 0.62 [71]. Both deterministic and stochastic processes simultaneously affect the community assembly of eukaryotic plankton, and random patterns (e.g., ecological drifts) are particularly pronounced for the rare taxa. Thus, 999 Erdős-Rényi random networks with the same number of nodes and edges as the real networks were generated in the igraph package 1.2.4 [67] to explore randomness of real networks. Z-tests were used to compare the topological parameters between real and random networks. t-tests were used to compare the topological parameters of two real networks. All the statistical analyses were carried out in R 3.4.1 [66].

Community Composition and Alpha Diversity
Overall, 1,146,504 high-quality sequences and 2985 OTUs at 97% nucleotide similarity level were yield from the sequencing of 18S rRNA genes. In the dry season, 2213 OTUs were detected, within which 32 OTUs were abundant representing 48.1% of all sequences and 2095 OTUs were rare representing 28.3% of all sequences (Figures 2 and 3). In wet season, 2132 OTU were detected, within which 32 OTUs were abundant and representing 43.5% of all sequences and 1999 OTUs were rare representing 29.6% of all sequences (Figures 2 and 3). Venn diagram shown that 1360 OTUs were shared by dry and wet seasons, within which only three were abundant in both dry and wet seasons ( Figure 2a). Abundant taxa had significantly higher niche width values than rare taxa ( Figure S2). The number of observed OTUs, Chao 1, ACE (abundance based coverage estimator), and phylogenetic diversity were higher in the wet season than in the dry season, while Simpson and Shannon index were not significantly different between two seasons ( Figure 4).   In the dry and wet seasons, SAR, Cryptophyceae, Opisthokonta, Archaeplastida, and Excavata were major super groups (relative abundance > 1%, Figure S3). The wet season had significantly higher Archaeplastida and Excavata than dry season ( Figure S3). For the dominant lineages (relative abundance > 1%) of the whole communities (all taxa), the relative abundances of Cryptomonas and Ciliophora were significant higher in the dry season, while Chlorophyta, Discicristata, Fungi, and Ochrophyta were significantly higher in the wet season ( Figure 3). Most of the lineages in abundant and rare subcommunities were also significantly different between seasons (Figure 3). For abundant subcommunities (abundant taxa), there were significantly higher numbers of Chlorophyta, Fungi, Ochrophyta, and Protalveolata in the wet season, which similar with the community composition of the whole communities. However, the composition of rare subcommunities (rare taxa) was quite different, Chlorophyta, Discicristata, and Dinoflagellata were significantly higher in the wet season, while Bicosoecida, Ochrophyta, and Protaleolata were higher in the dry season ( Figure 3). Moreover, principal coordinates analyses (PCoA) revealed that eukaryotic communities (all, abundant, and rare) were strongly clustered by seasons (Figure 2b). The results of ADONIS and MRPP further suggested clear seasonal patterns of eukaryotic communities in the dry and wet seasons (Figure 2c). Besides of the sample D04, there were no spatial variation of microeukaryotic community structure within the dry season (Figure 5a). In the wet season, the main microeukaryotic community structure showed spatial variation between the northern part of the lake samples (W01-W04) and the other lake samples (W05-W10) (Figure 5b).

Relationships between Community and Environmental Variables
Redundancy analyses (RDA) indicated that the variation of eukaryotic communities across dry and wet seasons were significantly affected by the variations of nutrient factors, including DOC, TN, − +

Relationships between Community and Environmental Variables
Redundancy analyses (RDA) indicated that the variation of eukaryotic communities across dry and wet seasons were significantly affected by the variations of nutrient factors, including DOC, TN, NO 3 − , NH 4 + , TP, and SRP (p < 0.05, Figure 6a). Variance partitioning analyses (VPA) revealed that nutrient factors explained 53.5%, 61.2%, and 53.7% of eukaryotic composition variations across dry and wet seasons for the overall communities, abundant subcommunities, and rare subcommunities, respectively ( Figure 6b). P and N factors, as well as the interaction of C, N, and P factors (C × N × P) had high contributions to the overall explanation (Figure 6b).  In the dry season, the results of RDA showed that taxonomic dissimilarities were significantly associated with DOC, TN, NO 3 − , and SRP (Figure 7a). Rare subcommunities also associated with NH 4 + . In the wet season, however, taxonomic dissimilarities were only associated with NO 3 − and TP ( Figure 7b). Abundant and rare subcommunities as well as overall communities showed similar patterns to nutrients, suggesting strong seasonal patterns and nutrient influences on eukaryotic communities. variables (arrows) across different seasons in terms of overall communities as well as abundant and rare subcommunities. All nutrient variables had goodness of fit at the significant level p < 0.05 by envfit function (999 permutations). (b) Variance partition analysis (VPA) determined the relative contributions of C-factor (including DOC), N-factor (including TN, NO3 − , and NH4 + ), P-factor (including TP and SRP), and the interactions between two or three of these factors (C × N, C × P, N × P, and C × N × P).

Co-Occurrence Network
Co-occurrence networks of eukaryotes were built based on spearman correlation relationships, and several topological parameters of the networks were calculated to reveal community assembly rules and interactions in different seasons ( Figure 8 and Table 1). The dry season networks consisted of 574 nodes (i.e., OTUs) by 2263 edges (i.e., significant interactions), with only 11.1% of the correlations were negative (Figure 8 and Table 1). The networks in the wet season consisted of 682 nodes by 5151 edges, with 38.9% of the correlations were negative ( Figure 8 and Table 1). The results suggested stronger competing relationships in the wet season than in the dry season. The networks exhibited a scale-free characteristic with the degree distribution follows a power law (p < 0.001, Figure S4), indicating that the networks were non-random. Through comparing these topological parameters (modularity, average clustering coefficient, network heterogeneity, and average number of neighbors), we can get that real network in both dry and wet seasons was the complexity of ecological networks rather than the random networks. The strong interactions of microeukaryotic communities in real networks suggesting that these no-random networks in both dry and wet seasons had 'small world' properties and modular structure (modularity > 0.4, Table 1). The dry season network was clearly parsed into five major modules (Figure 8b and Table S4). Modules A, B, C, D, and E accounted for 24.7%, 16.7%, 16.6%, 8.9%, and 8.5% of the whole network. According to the within and between module connectivity (Zi and Pi, respectively), there was one module hub and four connectors in the dry season network ( Figure S5 and Table S3). The wet season network was also clearly parsed into five major modules (Figure 8b and Figure S4). Modules A, B, C, D, and E accounted for 28.2%, 23.5%, 12.9%, 10.0%, and 5.3% of the whole network. There were 14 module hubs and 4 connectors in the wet season network ( Figure S6 and Table S3). Water 2020, 12, x FOR PEER REVIEW 11 of 19

Discussion
WLFs strongly impair the community structures and decrease biodiversity of many lake biotasincluding bacterioplankton, periphyton, macroinvertebrates, macrophytes, and fish [14,16,[18][19][20][21]especially due to hydrological regime disruption, water quality deterioration, as well as habitat loss at very low water level [19,20,72,73]. Microeukaryotic communities in PYL were also influenced by seasonal WLFs. In Poyang Lake, the dry season and wet season had distinct microeukaryotic communities. The relative abundances of Cryptomonas (3.69-18.64%) and Ciliophora (13.86-40.98%) were significant higher in the dry season, while Chlorophyta (6.31-22.10%), Fungi (1.62-6.81%), and Ochrophyta (6.88-45.67%) were significantly higher in the wet season (Figure 3). Cryptomonas and Ciliophora favor growing in nutrient enrichment environment [29], thus, the higher concentration of N (TN, NO 3 − , NH 4 + ) and SPR in the-dry season might promote the growth of Cryptomonas and Ciliophora. Temperature, organic matter and P content were the key factor of Chlorophyta, Fungi, and Ochrophyta [29,42], moreover, temperature and P could simultaneous influencing these taxa. In the present study, the high temperature and high content of DOC during the wet season could provide suitable environmental for growing of these taxa. There were higher observed OTUs, Chao 1, and ACE in the wet season than that in the dry season ( Figure 4). Previous studies found that increase of water temperature would favor the growth of microeukaryotic species [64]. Moreover, seasonal variations in microeukaryotic communities could be driven by the variation in chemical, physical, and biological parameters in their habitats [44]. Thus, higher nutrient concentrations related to low water lever would promote the growth of dominant eukaryote communities which caused lower biodiversity [31,44,64]. Microeukaryotes can be considered to be an indicator for environment change and freshwater ecosystem health, due to their rapid and strong responses to environmental disturbance [39]. More importantly, abundant subcommunities and rare subcommunities also had clear seasonal pattern (Figure 3), and these different taxa did not respond equally to the seasonal changes and linked water level fluctuations. Abundant taxa showed greater differences in community composition between two seasons, and the shared OTUs of abundant taxa between the dry season and wet season was lower than that of rare ones (Figure 2). Previous study also indicated that rare plankton communities could determine the changes of eukaryotic plankton communities in disturbed environment [31]. This result indicated that rare taxa include more ecologically redundant of the communities, and could provide biological buffering capacity to withstand environmental changes (e.g., seasonal WFLs in PYL) as previous observed in freshwater ecosystems [31,35].
Many studies suggest that nutrients are potential drivers of aquatic organism communities in freshwater ecosystems [16,74,75]. Our study also found nutrient variations had high explanations not only for the seasonal variations of microeukaryotic communities but also for the seasonal pattern of abundant subcommunities. The response of bacterial communities on WLFs found that DOC, NO 3 − , and SRP was the key factors influencing bacterial communities in the dry season, and TP was the dominant factor in the wet season [4]. Compared with bacterial communities, we can find that the response of microeukaryotic communities to nutrients in different seasons was the same as that of bacterial communities. DOC, TN, NO 3 − , and SRP had strong effects on microeukaryotic communities in the dry season, while NO 3 − and TP were the key drivers on microeukaryotic community structure in the wet season (Figures 6 and 7). WLFs could control nutrient distribution during a hydrological year in lake ecosystems [12,13]. Water quality usually has fast responses to WLFs with high nutrients concentrations and deterioration in the dry season due to low water capacity, low biological nutrient consumption, as well as high nutrient release from sediment [12,[76][77][78]. Previous studies in Poyang Lake revealed that water quality was controlled by WLFs with high nitrogen concentration in the dry season [77,79]. Moreover, DOC and phosphorus were also related to hydrological regimes [4,53,77], because there were high autochthonous and allochthonous sources contributing to the DOC pool in the lake ecosystem during the rainy seasons [80]. Stormwater also contributes high loading of TP with a high proportion of particulate phosphorus during the wet season [81][82][83][84]. The seasonal pattern of nutrients might cause seasonal responses of microeukaryotic communities to nutrients. Thus, seasonal nutrient variables in PYL were crucial drivers not only of bacterial communities, but also of microeukaryotic community dynamic patterns during seasonal WLFs. Aquatic organism communities are typically complex and compose of taxa with strong interactions. For highly complex microbial communities, co-occurrence patterns are ubiquitous and particularly important in understanding community assembly rules and potential biotic interactions [85,86]. The dry season network had higher average path length, clustering coefficient, and centralization betweenness, suggesting that microeukaryotic communities in the dry season constituted a more connected and complex co-occurrence network than in the wet season. In general, communities with low diversity, tight co-occurrence interactions, and high complexity have lower stability and are more susceptible to disturbance [87][88][89], which suggested that microeukaryotic communities in the dry season were more sensitive and vulnerable to water level perturbations in PYL.
Microeukaryotes in both dry and wet seasons exhibited modular structures (Table 1). Modularity is a characteristic of most large complex systems in which highly interconnected species are grouped into a module of a biotic network [70,90]. Environmental variations can drive the pattern of the co-occurrence networks by affecting a module [31,91]. There were many modules parsed in the networks (Figure 8). Some modules showed similar responses to environmental variables as the overall communities, while some modules only had very weak or no significant relationships with environmental variables. The results suggested that WLFs and concomitant change of nutrients had different influences on different groups of taxa. For example, module A, B, C, and D in the dry season were dominated by Ochrophyta that were mostly affiliated to rare taxa (Figures 3 and 8). DOC, nitrogen, and phosphorus were important nutrient factors which could promote the growth of Ochrophyta [42]. In the wet season, module B, C, and E were dominated by Chlorophyta (Figures 3 and 8), which could be strongly driven by high temperature, high light, and high organic matter [43,92]. Like the overall microeukaryotic communities, most major modules in the dry season were closely correlated with nutrients, for example, high concentration of available nitrogen and SRP driven the growth of Cryptomonas and Ciliophora [33,93,94]. However, in the wet season, most major modules only had weak or no correlations with nutrients. The results suggest that nutrients had stronger effects on microeukaryotic communities in the dry season at low water level.
With respect to modularity, different taxa played different roles in the networks. In our study, the dry season network had 1 module hub and 4 connectors, while the wet season network had 14 module hubs and 4 connectors. Connector nodes are responsible for inter-module fluxes, and when these modules are poorly connected or not connected at all to each other, the elimination of these nodes would probably have a large impact on the global structure of fluxes in the network [71]. Moreover, although the pathways in which module hubs are involved within the module, elimination of these nodes may also have a large impact on the networks especially when a network is mainly controlled by a large module [71]. Thus, module hubs and connectors were the keystone taxa which are crucial in maintain network structures and ecosystems stability. In PYL, Chlorophyta, Ochrophyta, and Ciliophora were the module hubs and connectors in the microeukaryotic networks, which were usually acted as the dominant phylum in many other lake and reservoir ecosystems [31,95]. Moreover, most of these module hubs and connectors were rare OTUs in PYL, suggesting that rare taxa can act as keystone taxa in the microeukaryotic networks, and the disappearance of these taxa could lead to the breaking apart of modules and networks [31]. Thus, less module hubs in the dry season network further suggested that microeukaryotic communities in the dry season are more vulnerable to environment perturbation.
In the present study, we mainly explored the influence of intra-annual WLFs on plankton communities, but inter-annual WLFs will also be important to elucidate the long-term effects of WFLs in the further study. In addition, although in the current study, we mainly focus on the effects of nutrient variations on microeukaryotic communities under the effect of WFLs, more environmental variables-such as pH, temperature, turbidity, and organic pollution-are also important in the future research related to WFLs.

Conclusions
Excessive WLFs impose dramatic influences on lake ecosystems and have been regarded as a pervasive threat globally mainly due to anthropogenic landscape modification reducing wetlands and building dams, which has reduced the natural water buffering capacity of the landscape. Our results demonstrated that seasonal environmental changes and WLFs were correlated to changes in microeukaryotic communities in Poyang Lake. The dry season and wet season had distinct microeukaryotic community compositions and structures. Moreover, the dry season had lower alpha-diversity than the wet season. Seasonal community patterns were strongly controlled by seasonal nutrient patterns caused partly by WLFs. In different seasons, microeukaryotic communities responded differently to nutrient variations. Nutrients had stronger correlation to microeukaryotic communities in the dry season. Less module hubs and tight co-occurrence interactions in the dry season network suggested that microeukaryotic communities in the dry season are more sensitive and vulnerable to environment perturbation. Network analyses further revealed different assembly mechanisms of microeukaryotic communities in different seasons, which supported that cooccurrence patterns of microeukaryotic communities were not static, but they changed over time. Thereby highlighting that the correlations of microeukaryotic communities strongly depend on the environmental fluctuation caused by WLFs in Poyang Lake. Our results add knowledge to the understanding of seasonal dynamics of lake ecosystems in response to seasonal WLFs and provide insight into assembly mechanisms of microeukaryotic communities.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2073-4441/12/9/2317/s1. Figure S1. Average water depth of Poyang Lake from January to December. Figure S2. Niche width value of microeukaryotes in different subcommunities in (a) dry season and (b) wet season. The top, middle, and bottom lines of the box indicate the 75th quartile, median, and 25th quartile, respectively. The whiskers indicate the minimum and maximum. The black dots indicate outliers. Different letters above each box indicate significant differences (ANOVA, p < 0.05). Figure S3. Chord diagram showing the relative abundances of microeukaryotic super groups in dry season and wet season. Statistical significance between dry season and wet season was assessed by t-test and indicated by asterisks (** represents p < 0.01). Figure S4. The node degree distributions of real co-occurrence networks (colored) and random networks (grey) of microeukaryotic communities in (a) the dry season and (b) wet season. Figure S5. Composition of modules. Pie charts showing the proportion of nodes number in (a) axonomic groups and (b) subcommunities. Figure S6. Zi-Pi plot indicates the distribution of OTUs based on their topological roles. Each dot represents an OTU colored by (a) taxonomic groups and (b). The topological roles were determined according to the values of within-module connectivity (Zi = 2.5) and among-module connectivity (Pi = 0.62). The detail of the module hubs and connectors are shown in Table S3. Table S1. GPS coordinates of the sample sites in the dry and wet seasons of Poyang Lake. Table S2. Nutrient concentrations (mg/L) in the dry and wet seasons of Poyang Lake. The minimum (Min), maximum (Max), and mean values, as well as standard deviation (SD) are shown. Statistical significance between dry and wet seasons was assessed by t-test. Table S3. List of module hubs and connectors in co-occurrence networks according to the connectivity of each node. Table S4. Mantel test between environmental variables and modules in different seasons. Spearman correlations were calculated. * represents p < 0.05 and ** represents p < 0.01.