DNA Extraction, PCR Amplification, and Products Purification
Genomic DNA was extracted from environmental samples using E.Z.N.ATM Mag-Bind DNA Kit (OMEGA) as per manufacturer’s instructions. Agarose gel was used to detect the purity and concentration of DNA extracts, and then selected the qualified samples for subsequent analysis. The V3-V4 region of 16S rDNA gene was amplified using primers 341F (5’-CCTACGGGNGGCWGCAG-3’) and 805R (5’-GACTACHVGGGTATCTAATCC-3’) (Herlemann et al., 2011); the V4 region of 18S rDNA gene was amplified using primers V4F (5’-GGCAAGTCTGGTGCCAG-3’) and V4R (5’-ACGGTATCTRATCRTCTTCG-3’) (Sun et al., 2015). The PCR amplification cycle was: the first round of amplification was 94°C for 3 min, followed by 5 cycles at 94°C for 30 s, 45°C for 20 s, 65°C for 30 s, then followed by 20 cycles at 94°C for 20 s, 55°C for 20 s, 72°C for 30 s, and a final extension at 72°C for 5 min; the second round of amplification was 95°C for 3 min, followed by 5 cycles at 94°C for 20 s, 55°C for 20 s, 72°C for 30 s, and a final extension at 72°C for 5 min. Two round PCR reactions were carried out in 30 μL reactions, including 15 μL 2×Taq master Mix, 10 μM of forward and reverse primers, and 20 ng genomic DNA or PCR products of first round amplification, respectively.
Sequencing and Data Analysis
The gene amplicons were sequenced on the Illumina Miseq platform, effective tag readings were gained by trimming the terminal bases, removing primer sequences and barcodes, and filtering the low-quality bases. The non-specific amplified sequences were removed using Usearch, chimeras were detected by Uchime algorithm (version 4.2.40) (Edgar, Haas, Clemente, Christopher, & Rob, 2011). Subsequently, we then blastn aligned the representative sequence in the database with the chimeric-removed sequence. Alignment results below the threshold were considered to be the sequence outside the target region and eliminated this part of the sequence (Haas et al., 2011). Then using the open reference OTU to cluster the qualified sequences into the species-level OTU (Operational Taxonomic Units), and the similarity cut-off value of using UCLUST in QIIME was 97%. The classification of rare and abundant species relies on the cutoff level of relative abundance, set rare taxa to 0.1 or 0.01% and abundant taxa to 1% (Zhang et al., 2018). In our study, we classified OTUs with a relative abundance greater than 1% as abundant taxa and less than 0.01% as rare taxa. In order to gain the species classification information corresponding to each OTU, the OTU representative sequence was aligned with the Silva database using the RDP Classifier algorithm with an 80% confidence threshold, and annotate species information at various levels (Wang, Garrity, Tiedje, & Cole, 2007). All raw sequences data in this study were deposited into the NCBI Sequence Read Archive (SRA) database under the accession number of SRP273360.
The co-occurrence network of each season was established through Molecular Ecological Network Analysis (MENA) based on the relative abundance of OTUs (Deng et al., 2012; Zhou et al., 2010), which computed the Pearson Correlation Coefficients (PCC) for each OTU pair, and then calculated the statistical significance of PCC values by a permutation test. The similarity threshold of microbial communities was determined according to the method of random matrix theory. And the correlation is significant (P < 0.01) when the correlation between the two OTUs is greater than the similarity threshold. Edges were set between pairs of OTUs which the PCC was significant. Then, the networks were visualized with Gephi and Cytoscape (Assenov, Ramirez, Schelhorn, Lengauer, & Albrecht, 2008). At the same time, in order to compare the difference between the molecular ecological network and the corresponding random network, the topological properties based on Erdós–Réyni random network were calculated using R language.
Then, we used the greedy modularity optimization method to characterize the network modularity of each network established in each season (Shi et al., 2016). Modules greater than 0.4 were used as the threshold of the detection module. In order to determine the topological properties of each node in the network, the nodes were divided four categories through their among-module connectivity (Pi ) and within-module connectivity (Zi ): Network hubs (Zi > 2.4 and Pi > 0.6), Module hubs (Zi > 2.5 and Pi< 0.6), Connectors (Zi < 2.5 and Pi> 0.6), and Peripherals (Zi < 2.5 andPi > 0.6) (Deng et al., 2012; Zhou et al., 2010).
PICRUSt is a bioinformatics tool that predicts metagenome gene functional content using 16S rDNA. PICRUSt uses the existing gene content annotations in the Integrated Microbial Genome (IMG) database (Markowitz et al., 2012) v3.5 (2,590 genomes) as well as the 16S copy number and functional classification scheme of the reference microbial genome to classify the expected metagenome content. And the current PICRUSt software (version 1.1.0) supports two types of functional prediction. We applied the more popular KEGG Orthologs in this study (Kanehisa, Goto, Sato, Furumichi, & Tanabe, 2012) for subsequent analysis.
Results
Physical, chemical, and biotic characteristics of water quality
The nutrient concentrations, physical properties, and Chl a concentration in temporal were shown in Table 1. The lowest values of nitrate, nitrite, phosphate, Chl a, and DOC were detected in winter. The concentrations of nitrate in spring were generally higher than in other seasons. The nitrite and phosphate concentrations in summer ranged from 0.083-0.170 mg L-1, 0.128-0.212 mg L-1, respectively, which were also higher than those in other seasons. The level of Chl a in autumn was far in excess of the other seasons. The maximum of average DOC concentration was detected in spring.
Seasonal succession of bacterial and eukaryotic phytoplankton communities
The bacterial community composition at the genera level was shown in Fig. 2. Our results demonstrated that Proteobacteria (12.35-92.4%) was the most abundant bacterial phyla, followed by Firmicutes (0.9-8.19%), Bacteroidetes (2.83-7.11%), Planctomycetes (0.94-5.27%), Actinobacteria (0.49-4.64%), Verrucomicrobia (0.14-1.78%), and Cyanobacteria (0.04-2.1%). Cyanobacteria were represented by Synechocystis andPlanktothrix . At higher taxonomic resolution, the most enriched genus based on relative abundance include Comamonas ,Straphylococcus , Sphingomonas , Acinetobacter ,Brevundimonas , Vogesella , Citrobacter , andClostridium in spring; Acinetobacter ,Chryseobacterium , Comamonas , Exiguobacterium ,Pseudomonas , and Sphingomona in summer;Acinetobacter , GpIIa , Chryseobacterium ,Pseudomonas , Comamonas , Stenotrophomonas ,Flavobacterium , Exiguobacterium , and Sphingobium in autumn; Pseudomonas , Flavobacterium , Blastomonas ,Undibacterium , Clostridium , and Janthinobacteriumin winter.
The variations in eukaryotic phytoplankton community composition at the genera level were shown in Fig. 3. The identified sequences were assigned that eukaryotic phytoplankton phylum was composed of Chlorophyta (75.4-93.67%), Bacillariophyta (4.13-17.23%), and Ochrophyta (0.19-43.27%). Results demonstrated that Chlorophyta were the major dominate phyla in Fenhe River. The peak of Ochrophyta occurred in winter. Desmodesmus was the dominant genus in the phylum of Chlorophyta, followed by Mychonastes ,Nephrochlamys , Coelastrum , Discostella , andChlorella , and in the phylum of Bacillariophyta Cyclotellawas the main genus in spring. Mychonastes , Desmodesmus ,Pseudopediastrum , Nephrochlamys , and Golenkiniawere the most enriched genus in summer. The community was characterized by green algae Mychonastes , Scenedesmus ,Polyedriopsis , Neodesmus , and diatoms Discostellain autumn. The dominant genera were Monoraphidium ,Parietochloris , Chlorella , Nannochloropsis , andLessonia in winter.
Co-occurrence networks ofbacterial communities
The interaction between different plankton plays a vital role in shaping the distribution pattern of microorganisms. They do not exist in the natural environment in isolation, but form a complex ecological interaction network. Aquatic ecosystems are composed of networks connecting different organisms of different sizes, and are maintained through species-to-species interactions (including symbiosis, competition, and parasitism) (Zhu, Hong, Zada, Hu, & Wang, 2018). Therefore, co-occurrence patterns of each season were established to explore the succession of OTU associations, which could not only reflect the potential interactions between microbes in complex communities, but also reflect ecological processes such as historical effects, cooperation, and habitat filtering.
The topological properties of bacterial networks over time were presented in Table S1. Overall, the bacterial network in spring was larger and more complex (Fig. 4). The proportion of positive interactions in spring, summer, and winter accounted for 100%. In the co-occurrence networks, the positive interaction was mainly considered as cooperation (Ju, Xia, Guo, Wang, & Zhang, 2014). Although in the autumn network, in addition to the positive interactions accounted for 82%, the negative interactions also accounted for 18%. However, the numbers of positively correlated connections between nodes were always higher than the negatively correlated connections, indicating that more mutual cooperation can make the microbial communities better resistant to changes in the external environment, thereby maintaining the stability of microbial network (Xue et al., 2018). Compare some network topology parameters of empirical network and the Erdös-Réyni random network, including the modularity index, the average path length, and the average clustering coefficient. The results demonstrated that the parameters of the empirical network were significantly higher than random network, which reflected that the network structure was not randomly distributed. The modularity of all networks was higher than 0.4, demonstrating that networks had the modular structure composed of closely connected nodes and forming a ”small world” topology.
Among the nodes of the bacterial networks, Proteobacteria occupied the largest number of OTUs in all seasons. And Proteobacteria occurred in almost all of the sub-network regardless of season (Fig. 4). These results demonstrated that some Proteobacteria members can adapt to different ecological environments, and this adaptation could also explain the great abundance of Proteobacteria in all sites in each season. Therefore, bacteria of Proteobacteria always maintained a higher abundance compared with other phyla. The network is divided into Network hubs, Module hubs, Connectors, and Peripherals according to the among-module connectivity (Pi ) and within-module connectivity (Zi ). OTUs that are divided into Network hubs, Module hubs, and Connectors are regarded as keystone species (Fig. S1 & S2). These nodes can keep the network at the lowest hierarchical structure under the premise of completing its ecological functions (Fan, Weisenhorn, Gilbert, & Chu, 2018).
Networks hubs were only detected in spring and autumn. In the spring network, one Network hubs identified originated from Proteobacteria (Pseudomonas ); two Module hubs belonging to Planctomycetes (Pirellula ) and Proteobacteria (Sphingomonas ); five OTUs were divided into Connectors, which were two Proteobacteria OTUs, two Bacteroidetes OTUs, and one Planctomycetes OTU (Gemmata ), respectively. In the summer network, one OTU was divided into Module hubs, belonging to Verrucomicrobia at the phylum level; six OTUs were classified as Connectors, the genera that could be identified as keystone taxa were mainly Luteimonas , Gemmata ,Mycobacterium , Stenotrophobacter , and Sphingopyxis . In the autumn network, the Bosea genus was divided into Network hubs; two out of three Module hubs were Proteobacteria (Agrobacterium and Hyphomicrobium ) and one was Bacteroidetes (Chryseobacterium ); four OTUs were classified as Connectors, genera that could be identified as key taxa includedPseudomonas , Acidovarax , and Rhodoplanes . In the winter network, three Module hubs identified originated from Proteobacteria and Planctomycetes; two Connectors belonging to Verrucomicrobia and Proteobacteria at the phylum level.
Co-occurrence networks of eukaryotic phytoplankton communities
In the network topological properties of eukaryotic phytoplankton, the average path coefficient, average path length, and modularity index in the molecular ecological network were greater than the corresponding values of the random network, indicating that the network constructed in this study conformed to the characteristics of scale-free, small world, and modular, which can be used for the subsequent study of phytoplankton interaction. The number of nodes and connections in the spring network was significantly higher than in other seasons, with a larger and more complex network, followed by autumn (Fig. 5). The greater complexity of the microbial network usually contributes to the stability of the microbial communities. It can be seen that high and low temperatures are not conducive to the stability of eukaryotic phytoplankton communities. In summer and autumn, the network modularity was relatively higher, and the connections between microorganisms were closer. In the summer and winter networks, positive correlation accounted for 47% and 43%, negative correlation accounted for 53% and 57%, respectively. On the contrary, in spring and autumn, the algal communities almost completely tended to co-express, that is, the positive relationship accounted for a higher proportion (spring: 100%, autumn: 94%). In the spring network, the Module hubs identified originated from Neodesmus ; Connectors between network came from Ochrophyta (Monodus ) and Chlorophyta (Tetradesmus ). In the summer network, Desmodesmus genus was divided into Network hubs; two OTUs were divided into the Module hubs, both of which were Chlorophyta, but only appraisal forMychonastes ; twelve Connectors were identified as keystone taxa were mainly Dangeardinia , Chlamydomonas ,Lagerheimia , Choricystis , Monoraphidium , andTetradesmus . In the autumn network, two Module hubs were Ochrophyta (Monodus ) and Chlorophyta; five Connectors mainly identified originated from Scenedesmus , Characiopsis , andCyclotella . In the winter network, one Bacillariophyta OTU (Amphora ) was divided into Network hubs; nine OTUs were classified as Connectors, all of which were Chlorophyta at the phylum level, the identified genera including Pseudomuriella ,Messastrum , Mychonastes , and Monoraphidium .