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 .