1 INTRODUCTION

Urban water bodies are important to human life. Unfortunately, they are suffering serious pollution in recent decades due to anthropogenic activities, posing a threat to the health of aquatic ecosystems, and having changed the physical and chemical environment and microbial community composition (Gücker et al., 2011; Drury et al., 2013). Previous researches have demonstrated that eutrophication by nitrate, phosphate, ammonium, and other organic matter pollutant (permanganate and suspended solids, etc.) had adverse effects on the functions of ecological community and aquatic ecosystem (Toyama et al., 2020; Zhang et al., 2020). As an important part of an aquatic ecosystem, microbial community plays an irreplaceable role in participating in biogeochemical processes and nutrient cycling.

In a natural ecosystem, an individual organism does not exist alone, but coexists with each other to form a network of ecological interactions (Mikhailov et al., 2019). These complex associations ease the influence of biodiversity on ecosystem functions (Duffy et al., 2007; Du et al., 2020). Microorganisms (such as bacteria and microbial eukaryotes) are important contributors to biodiversity as ecological role-players involved in primary production, food web, and biogeochemical cycle (Xue et al., 2018). Community structure is not only affected by interactions between species, but also by environmental fluctuations and regional conditions (Jones et al., 2013; Capo et al., 2017). In addition, the patterns of co-occurrence of bacterial and eukaryotic phytoplankton communities revealed based on next-generation sequencing can predict whether the ecological interactions among environmental variables and various species in aquatic ecosystems are positive or negative (Paver et al., 2013; Bunse et al., 2016), which is important for understanding the structure of microbial communities and provide insights into potential interaction networks (Ma et al., 2016).

Network analysis can determine the keystone species that are necessary to maintain community stability among many species. Vick-Majors et al. (2014) and Martín González et al. (2010) demonstrated that the nodes with higher betweenness centrality values are particularly important in maintaining the connectivity of ecological network and identify them as keystone species. A keystone species is a species that plays a critical and unique role in local natural environment although its abundance may not be high (Paine, 1969). Keystone species can affect the entire microbial community through a series of pathways. For example, they can affect the community structure and function through adjusting some intermediate groups and effect groups. Besides, keystone species are also susceptible to dynamic environments, the disappearance or change of keystone species may cause disturbance to mature communities (Steele et al., 2011).

At present, there have been numerous reports on the temporal and spatial changes and co-occurrence dynamics of microbial communities (Hu et al., 2017; Jiao et al., 2020), but most of them focused on abundant taxa. A microbial community is composed of a few abundant taxa and many rare taxa. Abundant taxa are contributive to the organic matter flux and biomass yield (Pedrós-Alió 2012), whereas certain rare taxa may include more metabolically active microorganisms, which might be considered as keystone species that regulate the function of aquatic environment (Wang et al., 2020). In addition, rare taxa contribute to the stability and persistence of the community, becoming the "seed bank" of microbial and may become dominant under favorable conditions (Shade and Gilbert, 2015). At present, many studies have realize the importance of rare taxa in microbial ecology. For example, Jousset et al. (2017) reported that in terms of taxonomy and functional diversity, rare taxa might increase functional redundancy and enhance the ability of communities against environmental interference, which is critical to the nutrient cycles. Mallon et al. (2015) pointed out that some indigenous rare taxa in the environment could control the invasion of foreign microorganisms (such as pathogens), and indicated that rare microorganisms play an irreplaceable role in the ecosystem. However, due to the limitation in method, the rare biosphere has been largely unexplored. Recently, rapid development of next-generation sequencing and deep sequencing promoted the research on rare biosphere communities, and the ecological importance of rare sub-community attracted widespread attention (Wu et al., 2017).

One of the main purposes of current microbial ecological research is to understand the formation and maintenance mechanism of microbial diversity (including abundance and composition). On this issue, the most mainstream view is that the microbial community structure is influenced by both deterministic processes based on niche theory (including species interactions, environmental selection, etc.) and stochastic processes based on neutral theory (including diffusion limitation, drift, etc.) (Jiao et al., 2020). The neutral theory predicts that community structure and composition are relevant to the geographic distance between samples, which is due to the dispersal restrictions and many species have equal functions in ability to utilize a niche. According to the niche theory, changes in species community composition are related to variations in environmental parameters (Dumbrell et al., 2010), because species with unique characteristics could enable themselves to develop unique habitats for surviving.

The interrelationship networks are complex and can contain the inherent mechanism of microbial interaction in environmental disturbance transmission (Hunt and Ward, 2015). Therefore, identifying the interactions between bacteria and microeukaryotes and the combination of complex ecological communities is critical for the management of biodiversity, especially in the areas of strong human interference and fragile ecosystem. Although many studies are conducted on the interaction between abundant taxa (Ger et al., 2016), few are on how rare taxa affect the community structure of microorganisms in aquatic ecosystem. In addition, knowledge on co-occurrence patterns of rare taxa in urban water bodies remains poor. As the capital of Shanxi Province, Taiyuan is an important coal industrial base in China. Due to the excessive development and consumption of coal resources in the recent decades, the ecological environment in urban Taiyuan has been seriously polluted. At present, the microbial network in the Taiyuan section of the Fenhe River has not been systematically analyzed. Therefore, 16S and 18S rDNA gene sequencing was used to detect the seasonal pattern of bacterial and eukaryotic phytoplankton communities in the study area, and network analysis was used to identify keystone species. Moreover, the species were divided into abundant taxa and rare taxa, from which relationship between keystone taxa and environmental variables are discussed. This research will help understand the interactions and ecological effects of microorganisms and improve our ability to predict the response of microorganisms to environmental changes.

2 MATERIAL AND METHOD 2.1 Sampling and environmental parameter

Fenhe River is the first tributary of the Huanghe (Yellow) River with 39 471 km2 in area and 710 km in total length. The study area is located in Taiyuan (37°27′N–38°25′N, 111°30′E–113°09′E) (Fig. 1). The regional climate is characterized by a semi-humid continental monsoon climate belonging to a semi-enclosed water body. The mean annual precipitation is 504.8 mm. Besides, it is also the main coal-producing area, industrial area, and urban residential area of Shanxi Province.

For this study, 144 samples in 6 sampling stations were collected in 2019 on March 22, April 18, and May 14 for spring; on June 10, July 17, and August 12 for summer; on September 16, October 11, and November 19 for autumn; and on January 8, December 16 for winter. Three replicates were set for each sample. A 2-L sampler was used to col1ect water samples at the depth of 0.5 m, and then 1 000 mL of water volume was filtered through 0.22-μm polycarbonate filters (Milipore, USA). The > 0.22-μm fraction was used for DNA extraction, and the corresponding water samples were used to determine the physical and chemical parameters. Water temperature and pH were measured by a multi-parameter water quality meter (Hydrolab DS5, Hach Company, USA). Nitrate (NO₃^-), nitrite (NO₂^-), and phosphate (PO₄^3-) were measured as per Wang et al. (2008). The concentration of dissolved organic carbon (DOC) was determined using total organic carbon analyzer (Shimadzu Corporation, Kyoto, Japan). The water and wastewater monitoring analysis method (Zhu et al., 2019) was used to measure chlorophyll-a (Chl-a) concentrations. The significance (One-way ANOVA) of each parameter was calculated by SPSS 26.0.

2.2 DNA extraction, PCR amplification, and products purification

Genomic DNA was extracted from environmental samples using E.Z.N.A^TM 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 the qualified samples were selected for subsequent analysis. The V3–V4 regions 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., 2014). The PCR amplification cycle was: the first round of amplification was 94 ℃ for 3 min, followed by 5 cycles at 94 ℃ for 30 s, 45 ℃ for 20 s, 65 ℃ for 30 s, then followed by 20 cycles at 94 ℃ for 20 s, 55 ℃ for 20 s, 72 ℃ for 30 s, and final extension at 72 ℃ for 5 min. The second round of amplification was 95 ℃ for 3 min, followed by 5 cycles at 94 ℃ for 20 s, 55 ℃ for 20 s, 72 ℃ for 30 s, and a final extension at 72 ℃ for 5 min. Two-round PCR reactions were carried out in 30-μL solution including 15-μL Phusion High-Fidelity PCR Master Mix, 0.2 μmol/L of forward and reverse primers, and 10-ng genomic DNA or PCR products of first-round amplification. Samples were then sequenced on the Illumina MiSeq sequencing platform in the Sangon Biotech Co., Ltd. (Shanghai, China), and paired-end reads were generated by sequencing the forward and reverse strands of each target DNA fragment. The rarefaction curve of each sample in this study tended to be flat, and the coverage rate reached 99.9%, indicating that the depth of sequencing reads was large enough.

2.3 Processing sequence data

The QIIME 2 software package (version 2019.1) was used to process raw sequence data (Bolyen et al., 2019). The DADA2 inference algorithm can perform the primer-free reads to correct sequencing errors and create amplicon sequence variants (ASVs) for microbial communities (Callahan et al., 2016). Therefore, DADA2 was used to perform quality control, denoising, filtering, merging, and removing chimeras on these sequences to generate ASVs. ASVs are more sensitive, specific, and repeatable, and using ASVs can better distinguish ecological patterns than the OTU method (Callahan et al., 2017). All raw sequence data in this study were submitted to the NCBI Sequence Read Archive (SRA) database under the accession number SRP273360.

2.4 Statistical analysis

In this study, the diversity index Simpson and Shannon index for α-diversity analysis was calculated using Mothur software. The significance (ANOVA variance) analysis followed by Tukey's HSD post hoc tests of all indexes and environmental parameters was conducted by SPSS 26.0. Analysis of similarity (ANOSIM) was used to perform differences in microbial communities between groups by "vegan" package for R (version 4.0.3). The classification of rare and abundant species relies on the cutoff level of relative abundance, ASVs with relative abundance greater than 1% are defined as abundant taxa, less than 0.01% is defined as rare taxa (Zhang et al., 2018). Abundant taxa including always abundant taxa (AAT) and conditionally abundant taxa (CAT). Rare taxa including always rare taxa (ART) and conditionally rare taxa (CRT).

The co-occurrence network of each season was established through Molecular Ecological Network Analysis (MENA) based on the relative abundance of ASVs (Deng et al., 2012), by which the Pearson Correlation Coefficients (PCC) were computed for each ASV pair, and then the statistical significance of PCC values was calculated in permutation test. The similarity threshold of microbial communities was determined according to the method of random matrix theory. The Benjamini-Hochberg method was used to standardize the correlation coefficients to correct the significance P-value obtained from the original hypothesis test, and finally the corrected P-value was used to keep the relevant ASV of P < 0.05 to construct the correlation network. Moreover, we used the Network Analyzer tool to obtain network topology parameters such as the characteristic path length, the number of connections, the number of nodes, the aggregation coefficient, the network density, and the average connectivity of the network. Edges were set between pairs of ASVs, for which the PCC was significant. Then, the networks were visualized with Gephi (Version 0.9.2) (Heymann, 2014) and Cytoscape (version 3.2) (Shannon et al., 2003). Meanwhile, 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 (Version 4.0.3) (Lieberman et al., 2005).

In addiiton, 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. To determine the topological properties of each node in the network, the nodes were divided into four categories through their among-module connectivity (P_i) and within-module connectivity (Z_i): network hubs are nodes that are highly connected to other nodes in the entire network (Z_i > 2.5; P_i > 0.6), module hubs are nodes that are highly connected to other nodes in the network modules (Z_i > 2.5; P_i < 0.6), connectors are nodes that connect different modules (Z_i < 2.5; P_i > 0.6), and Peripherals are nodes that have less connectivity with other nodes (Z_i < 2.5; P_i < 0.6) (Deng et al., 2012).

Ordination analysis was performed by Canoco 5.0. Detrended correspondence analysis (DCA) was used to determine whether unimodal or linear ordination methods should be applied. Results show that the relationship between the bacteria and environmental variables was only applicable to the unimodal model (canonical correspondence analysis, CCA), and the relationship between the eukaryotic phytoplankton and the environmental variables applied to both redundancy analysis (RDA) and CCA. The Monte Carlo permutation tests (499 permutations) were applied to identify the environmental variables with significant effects (Sarker et al., 2013).

3 RESULT 3.1 Physical, chemical, and biotic characteristics of water quality

The characteristics of physical, chemical, and biotic are shown in Table 1. The lowest values of nitrate, nitrite, phosphate, Chl a, and DOC were detected in winter. Nitrate concentration in spring was generally higher than those in other seasons. Nitrite and phosphate concentrations in summer ranged 0.083–0.170 mg/L, 0.128–0.212 mg/L, respectively, which were also higher than those in other seasons. Water temperature was highest in summer, with an average of 27.4 ℃, and post hoc HSD tests indicated there was no significant difference between spring and autumn. The pH value in winter was higher than in other seasons. The level of Chl a in autumn was much greater than that in spring and winter (P < 0.01). The maximum average DOC concentration was detected in spring. One-way analysis of variance showed that all parameters were significantly different (P < 0.05), and the DOC concentration in spring was significantly higher than that in winter (P < 0.05).

3.2 Sequence quality and α-diversity indexes

The rarefaction curve obtained by measurement of the primers in this study tended to be flat, indicating that the amount of sequencing data was suitable. In addition, the coverage rate can also reflect whether the sequencing results represent the true situation of the sample. The actual test results of the primers all reached more than 99%, which reflected the real situation of the community composition in the tested sample. In this study, information of abundant taxa and rare taxa is shown in Table 2. Singletons were deleted before further analysis. The results show that the abundant taxa for bacteria accounted for only 0.23%–1.11% of the total sequences, while the number of rare taxa was more than 93% of total ASV number. Abundant taxa for phytoplankton accounted for 0.63%–14.75%, while rare taxa accounted for more than 63.87%. It can be seen that the number of rare taxa was much greater than those of abundant taxa.

The average number of high-quality sequences was the highest in autumn, 53 756 16S rDNA sequences and 23 690 18S rDNA sequences. However, the number of ASVs was higher in spring, with 327 for bacteria and 128 for eukaryotic phytoplankton (Supplementary Table S1). Sequences of 16S rDNA had an average length of 420 bp. Sequences of 18S rDNA had a mean length of 400 bp. At the ASV level of phytoplankton communities, ANOSIM was used to test the temporal and spatial similarity, and gave the results of R²=0.92, P < 0.01 (between seasons) and R²=0.04, P > 0.07 (between sites), which demonstrated that a significant difference existed across four seasons. Similarly, at the ASV level of bacteria, ANOSIM gave the results of R²=0.89, P < 0.01 (between seasons) and R²=0.11, P > 0.05 (between sites). The α-diversity indexes of bacteria and phytoplankton are shown in Supplementary Table S2. ANOVA test showed significant differences in diversity between seasons, indicating that when environmental factors change in season, and the species composition and abundance also change significantly. From the perspective of bacteria or phytoplankton, post hoc HSD tests indicated that the biodiversity in spring was significantly higher than that in autumn (P < 0.05). The positive correlation between nitrate and DOC and Shannon index indicated that the structure and composition of the microbial community in this study mainly respond to the fluctuation of carbon and nitrogen.

3.3 Seasonal variation of bacterial and eukaryotic phytoplankton community

The bacterial community composition at the genus level is shown in Fig. 2. Our results demonstrate that Proteobacteria was the most abundant bacterial phylum in all seasons. Besides, Bacteroidetes were also the dominant phylum in spring. Dominant bacteria in summer samples included Actinobacteria and Firmicutes, and autumn samples included Bacteroidetes and cyanobacteria. The relative abundances of Planctomycetes and Verrucomicrobia in winter were higher than those in other seasons. At higher taxonomic resolution, the most enriched genus based on relative abundance included Proteobacteria (Acinetobacter, Brevundimonas, Comamonas, Citrobacter), Bacteroidetes (Cloacibacterium) in spring, of which Acinetobacter and Comamonas were the abundant taxa, Brevundimonas and Citrobacter the rare taxa; Proteobacteria (Acinetobacter, Comamonas, Pseudomonas, and Sphingomonas), Bacteroidetes (Chryseobacterium), Firmicutes (Exiguobacterium) in summer, of which Acinetobacter, Pseudomonas, and Sphingomonas were the abundant taxa, Comamonas the rare taxa. The relative abundance of Acinetobacter was the highest in all samples in autumn, accounting for 80%–90%, only Acinetobacter, Delftia, and Klebsiella were the abundant taxa, all other bacteria were the rare taxa. The dominant genus in winter differed greatly from that in other seasons, mainly including Proteobacteria (Janthinobacterium, Pseudomonas, and Undibacterium), Bacteroidetes (Flavobacterium), Firmicutes (Clostridium) in winter, of which only Janthinobacterium was the abundant taxa, the others were the rare taxa. In addition to the genera listed in Fig. 2, the rare subcommunity also included 190 genera in spring, 138 genera in summer, 102 genera in autumn, and 141 genera in winter.

The variations in eukaryotic phytoplankton community composition at genus level are shown in Fig. 3. The eukaryotic phytoplankton was composed of Chlorophyta (56.87%–98.46%), Bacillariophyta (0.10%–35.99%), Cryptophyta (0.05%–34.27%), and Ochrophyta (0.06%–13.74%). Results demonstrate that Chlorophyta was the major dominant phylum in Fenhe River. The peak of Cryptophyta occurred in winter, and the abundance of Ochrophyta was highest in autumn. At genus level, Desmodesmus was the dominant genus in the phylum of Chlorophyta, followed by Mychonastes of Chlorophyta and Cyclotella of Bacillariophyta in spring, of which Desmodesmus and Mychonastes were the abundant taxa, Cyclotella was the rare taxa. In summer, Mychonastes, Golenkinia, and Pseudopediastrum were the abundant taxa, while Desmodesmus, Messastrum, and Monoraphidium with higher relative abundance were the rare taxa. The community was characterized by Nannochloropsis and Tetraëdron of Chlorophyta, as well as diatoms Discostella and Stephanodiscus in autumn, of which Stephanodiscus was the rare taxa. The dominant genera including green algae Mychonastes and Monoraphidium in winter, all phytoplankton were the rare taxa except for Mychonastes. In addition to the genera listed in Fig. 3, the rare subcommunity also included 40 genera in spring, 49 genera in summer, 50 genera in autumn, and 33 genera in winter.

3.4 Co-occurrence networks of bacterial communities

The interaction between different planktons plays a vital role in shaping the distribution pattern of microorganisms. 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 et al., 2018). Therefore, co-occurrence patterns of each season were dipicted to explore the associations between ASVs, which reflect the potential interactions between microbes in complex communities and ecological processes as well such as historical effects, cooperation, and habitat filtering.

The topological properties of bacterial networks over time are presented in Supplementary Table S3. We compared the network topology parameters of the empirical network with the Erdös-Réyni random network, including the modularity index, the average path length, and the average clustering coefficient. The results demonstrate that the parameters of the empirical network were significantly higher than the random network, reflecting that the network structure was not randomly distributed. The modularity of all networks was higher than 0.4, demonstrating that networks had a modular structure composed of closely connected nodes and forming a "small world" topology. Overall, the proportion of negative interactions accounted for 68%–100% of the potential correlations observed in four seasons (Fig. 4), indicating the co-exclusion (negative interactions, green lines) among taxa is greater than co-occurrence (positive interactions, red lines).

Among the nodes of bacterial networks, Proteobacteria occupied the largest number of ASVs in all seasons, with 43 ASVs in spring, 32 ASVs in summer, 22 ASVs in autumn, and 53 ASVs in winter. In addition, Proteobacteria occurred in almost all of the sub-network regardless of season (Fig. 4). These results demonstrate 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. The network is divided into network hubs, module hubs, connectors, and peripherals according to the amongmodule connectivity (P_i) and within-module connectivity (Z_i). ASVs that are divided into network hubs, module hubs, and connectors are regarded as keystone species (Supplementary Figs.S1–S2). These nodes can keep the network at the lowest hierarchical structure under the premise of completing its ecological functions (Fan et al., 2018). In order to determine the taxonomic composition of each network, we drew a phylogenetic diagram of the main modules.

Network hubs were detected in spring and autumn only (Table 3). In the spring network, one network hub identified originated from Proteobacteria (Sphingobium) with max betweenness centrality. In addition, nine ASVs were divided into connectors, which were seven Proteobacteria ASVs, one Bacteroidetes ASV (Flavobacterium), and one Firmicutes ASV (Clostridium), respectively. Moreover, three sub-modules were highly correlated. They belong to the rare taxa except ASV36 and ASV51. In addition to Proteobacteria, Verrucomicrobia and Planctomycetes were also predominate in sub-networks. In the summer network, two ASVs were divided into module hubs, belonging to Proteobacteria (Salinivibrio and Sphingopyxis), of which Salinivibrio affiliated to rare taxa had the highest betweenness centrality; three ASVs were classified as connectors, the genera that could be identified as keystone taxa were mainly Chryseobacterium, Polaromonas, and Mycobacterium, of which Polaromonas and Mycobacterium were rare taxa. The summer network was relatively diverse, the main components of which were Proteobacteria, Actinobacteria, and Bacteroidetes. The autumn network had fewer nodes and links, the only Acinetobacter genus was divided into network hub and six ASVs were classified as connectors. Genera that could be identified as key taxa included Acinetobacter, Salinivibrio, Enterobacter, and Legionella. In the winter network, one module hub identified was originated from Proteobacteria. And five ASVs were divided into connectors, belonging to Bacteroidetes and Proteobacteria at phylum level. The winter network had the highest modularity. Although the negative correlation accounts for the majority of the entire network, the sub-modules composed of ASV65 (Pseudomonas) and ASV569 (Pirellula) exhibited a high positive correlation (Fig. 4).

The interactions between ASVs within phylum or high-level taxa (Supplementary Fig.S3) indicated that taxonomically closely associated bacteria were also ecologically closely related, which in turn also reflected their synergistic relationships or common niche preferences. Apart from the co-occurrence patterns between phyla, such non-random patterns also occurred between different phyla. A typical example was that members of the phylum Proteobacteria showed significant correlations with Bacteroidetes, Planctomycetes, and Actinobacteria.

3.5 Co-occurrence networks of eukaryotic phytoplankton communities

In the network topological properties of eukaryotic phytoplankton (Supplementary Table S3), 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. In spring and summer, the network modularity was relatively higher, and the connections between microorganisms were closer (Fig. 5). In the summer and autumn networks, positive correlations accounted for 79% and 59%, negative correlations accounted for 21% and 41%, respectively. On the contrary, in winter, the algal communities almost completely tended to co-exclude, that is, the negative relationships accounted for a higher proportion (92%).

Keystone species identified by topological roles of eukaryotic phytoplankton networks are shown in Table 4. In the spring network, one network hub belonging to Scenedesmus (rare taxa) of Chlorophyta was identified. The module hubs that originated from Mychonastes with the maximal betweenness centrality were identified and the connectors came from Ochrophyta (Nannochloropsis) and Chlorophyta (Choricystis, Pseudoschroederia, Monoraphidium). The spring network presented a co-occurrence pattern between Chlorophyta and Bacillariophyta. ASV16 (Mychonastes) was an important hub connecting all modules, forming a sub-network with Cyclotella, Pseudopediastrum, and Nannochloris. In the summer network, Pseudoschroederia genus was divided into several module hubs; three ASVs could be divided into connectors, of which two were Chlorophyta and one belonging to Bacillariophyta, the identified genera including Mychonastes, Nannochloris, and Cyclotella. As the abundant taxa, ASV39 formed a sub-module with Messastrum, Rotundella, Monoraphidium, Pseudopediastrum, and Pseudoschroederia, then indirectly connected with other taxa through ASV95 (Pseudoschroederia). In the autumn network, one module hub was based on Cyclotella (Bacillariophyta), including four connectors mainly identified as being originated from Scenedesmus, Monoraphidium, Mychonastes, and Follicularia. Ochrophyta was abundant in autumn and closely related to Chlorophyta in the entire network. In the winter network, six ASVs were divided into several network hubs, of which four were Chlorophyta at the phylum level; and the identified genera included Desmodesmus, Mychonastes, ASV454 Pseudopediastrum, and Monoraphidium, and two ASVs were affiliated to Cryptomonas of Cryptophyta. Compared with other seasons, more Cryptophyta (Rhodomonas, Cryptomonas, Rhinomonas, Teleaulax) appeared in the winter network, and these taxa also played an important role in the sub-network. In addition, the co-expression patterns appeared in sub-modules composed of ASV331 (Mychonastes) only.

The less relevant co-occurrence pattern between ASVs within the phylum was observed (Supplementary Fig.S4), which is consistent with previous observation of taxa over-dispersion (no obvious clustering and correlation across taxa) in plant or animal natural communities (such as sedge community and mammalian community) (Ju et al., 2014). This not only confirmed that over-dispersion of phylogeny may be a common feature of all biological communities from microorganism to macro-plants and animals (Horner-Devine and Bohannan, 2006), but also showed the important influence of negative interaction (such as competition) on community assembly.

3.6 Relationship between keystone species and environmental factors

The information of keystone species in bacterial co-occurrence is shown in Table 3. Most of the keystone species detected were affiliated to the rare taxa, and only a small proportion belonging to the conditionally rare and abundant taxa (CRAT). CCA was performed to identify the effects of water quality physicochemical variables on these keystone taxa (Fig. 6a). The first two axes explained 35.33% and 31.02% of the cumulative variance of the relationship of species-environmental variables, and the eigenvalues were 0.960 3 and 0.843 1, respectively (Supplementary Table S4). The species-environment correlation values of axis 1 and axis 2 were 0.998 8 and 0.979 1, respectively. The Monte Carlo test showed that water temperature, nitrite, and DOC had significant effects on the bacterial community (P < 0.05). Flavobacterium, Salinivibrio, Mycobacterium, Chryseobacterium, and Sphingopyxis were affected by nitrite. Polaromonas showed significantly positively related to water temperature and phosphate. The most genus of Proteobacteria in spring showed a significant positive correlation with DOC.

The eukaryotic phytoplankton in the freshwater ecosystem responded rapidly and strongly to environmental disturbances and was considered the most important indicator to environmental changes and ecosystem status. Half of the keystone species detected in the eukaryotic phytoplankton co-occurrence patterns were affiliated to the rare taxa, and the rest belonged to the abundant taxa and CRAT (Table 4), indicating that rare taxa were crucial in the aquatic ecosystem. Therefore, RDA was performed to identify the effects of water quality physicochemical variables on these keystone taxa (Fig. 6b). The first two axes explained 71.43% and 22.00% of the cumulative variance of the relationship of species-environmental variables, and the eigenvalues were 0.454 8 and 0.140 0, respectively (Supplementary Table S4). The species-environment correlation values of axis 1 and axis 2 were 0.874 7 and 0.698 6, respectively. The Monte Carlo test showed that water temperature, nitrate, and phosphate had significant effects on the eukaryotic phytoplankton (P < 0.05). Water temperature was positively correlated with the abundance of Mychonastes, while negatively correlated with those of Cryptomonas, Desmodesmus, and Pseudopediastrum. DOC, nitrogen, and phosphorus (especially nitrate and phosphate) were positively related to Nannochlorisand Pseudoschroederia, whereas negatively related to Cyclotella and Scenedesmus. Follicularia and Monoraphidium showed a negative correlation with pH.

3.7 Co-occurrence patterns between bacterial communities and phytoplankton taxa

Interrelations between bacteria and microalgae are multifaceted and complicated; for example, bacteria can rely on photosynthetic phytoplankton to obtain the organic carbon needed to maintain their growth (Falkowski et al., 2008). In turn, phytoplankton depends on bacteria to remineralize organic matter into inorganic substitutes, ultimately supporting the growth of algae (Worden et al., 2015). The significant interaction between phytoplankton and bacteria has become widely examined in models of, and research on, the microbial loop, particularly in aquatic ecosystems (e.g., Azam et al., 1983; Williams and Ducklow, 2019). Therefore, the research of phycosphere and bacterial communities is important because they control the metabolic interactions of algal-bacteria in the microenvironment. We studied the non-random, co-occurrence networks, and important inter-taxa relationships using the network analysis approach (Fig. 7).

In Proteobacteria taxa, Bosea in spring, Brevundimonas in autumn, and Pseudomonas in winter could promote the growth of Desmodesmus, while Enterobacter and Ochrobactrum were mutually exclusive with Desmodesmus. Flavobacterium was positively correlated with most green algae, whereas it had a negative correlation with Chlorotetraëron and Goniomonas. The node of green algae Radiococcus had the highest connectivity in autumn, with 7 positively correlated nodes and 3 negatively nodes. As an abundant bacterium in winter, Janthinobacterium was only helpful to the growth of green algae Nannochloropsis. Comamonas belonging to Comamonadaceae family had a relatively higher abundance in spring and summer, which was considered as denitrifying polyphosphate accumulating microorganism (Calderer et al., 2014), which showed a positive relationship with Follicularia and Polyedriopsis of green algae.

Spring cyanobacteria had only a negative correlation with Sphingobium, and probably promoted the growth of the remaining bacteria (mainly Proteobacteria), Bacillariophyta (Amphora), and green algae (Golenkinia, Phacotus, and Nannochloris). Autumn cyanobacteria were more mutually exclusive with Bacteroidetes; specifically, Synechococcus was negatively related to Flavobacterium. Winter Mesorhizobium was also the prevalent (negative) group among the Synechococcus-associated bacteria.

Discostella is a relatively abundant genus in autumn, the co-existence of Discostella and bacterial Proteobacteria (especially Reyranella and Bosea) was more conducive to their growth compared to the individual population. The blooms of the small centric diatom Stephanodiscus were regularly found every year in lakes and reservoirs under low-temperature conditions (Kang et al., 2007). We found that Stephanodiscus was positively correlated with bacterial Caulobacter and Cyanobacterium, whereas negatively related to most green algae in winter. The relationships between Cyclotella and bacterial communities did not appear to be one-sided, Cyclotella constituted symbiotic relationships with Sphingopyxis, Novosphingobium, and Pseudoxanthomonas. Besides, Cyclotella also repelled each other with Chlamydomonas and Chloroidium of green algae and Cryptomonas of Cryptophyta. The genus Amphora is one of the larger genera belonging to the family Naviculaceae and widespread including freshwater, brackish water, and marine habitats (Nagumo, 2003). Amphora was positively related to cyanobacteria and most Proteobacteria (Acinetobacter, Hyphomicrobium, and Bradyrhizobium), and negatively related to Choricystis.

4 DISCUSSION 4.1 The ecological effects of keystone species (rare taxa)

Our results show that Proteobacteria, Actinobacteria, and Bacteroidetes were the dominant phyla of freshwater habitats, which is consistent with previous studies (e.g., Šimek et al., 2005). At the genus level, the seasonality showed a certain trend of change, because some genera dominated in one season while tended to be rare or non-existent in other seasons (Fig. 2). The reason might be related to the availability of nutrients and metabolites produced by the rare taxa, which are essential for the reproduction of the abundant taxa (Nyirabuhoro et al., 2020). Another reason might depend on the rarity threshold, below which the dormant groups cannot grow even under favorable conditions (Ruiz-González et al., 2015).

To further determine the contribution and role of abundant and rare taxa, ASVs were classified according to their relative abundance. However, this might cause certain species could not be divided into abundant taxa despite they belonged to a phylum with higher relative abundance. Regardless of this point, we still found that rare taxa were more important in the seasonal dynamics of microbial community because rare taxa contributed more than 93% of bacteria and 63% of eukaryotic phytoplankton, while the contribution rate of abundant taxa was less than 2%. The species rarity of rare taxa may be due to its very narrow environmental niche and is only abundant in a few habitats (Jousset et al., 2017). Besides, rare taxa are generally more sensitive to environmental fluctuations and likely to become extinct because of their low abundance (Zhou et al., 2020). The changing environment is also more likely to influence rare taxa. Therefore, most of the diversity composed of rare taxa may play a critical role in the microbial community composition and dynamics of the Fenhe River. Wang et al. (2020) pointed out that CRAT could switch between abundant and rare taxa under certain environmental conditions and play a key role in microbial community dynamics. Considering that CRAT also accounted for a large proportion in the keystone species of eukaryotic phytoplankton communities, it is of great significance for subsequent research to know the contribution of CRAT to co-occurrence patterns of microbial community.

Moreover, rare taxa are important to the co-occurrence dynamics of microbial community. The interaction between microbial taxa was analyzed through the topological characteristics of co-occurrence network. The co-occurrence network in this study showed the same non-random features and modular structure as the previous studies (Xue et al., 2018; Wang et al., 2020). The number of negative interactions in the seasonal co-occurrence network of bacteria was greater than those of positive interactions, indicating that bacterial taxa tend to compete each other. Among them, most of the negative interactions were observed within rare taxa (spring: 93, summer: 53, autumn: 50, winter: 132), while the interactions between rare and non-rare taxa exist in spring only, indicating that the competition within rare taxa was stronger than that between rare and non-rare taxa. For eukaryotic phytoplankton, the positive interactions within rare taxa (spring: 74, summer: 31, autumn: 6, winter: 15) and between rare and non-rare taxa (only summer: 22, autumn: 45) were higher than the negative interactions, indicating that the cooperation within rare taxa and between rare and non-rare taxa was stronger than that within non-rare taxa, which is beneficial to species isolation and modularization of the network (Hu et al., 2017) and maintain the stability of ecosystems. For example, the rare Symbiodinium taxa can significantly improve the stability of host-symbiont communities in a changing environment (Ziegler et al., 2018). Besides, the cooperation among these taxa also can contribute to the adaptability of microbial communities in complex environments, because the interaction networks can offer a buffer for resisting environmental disturbance (Konopka et al., 2015).

Rare microorganisms maintain a huge functional gene library, which can enhance indirectly the function of abundant microorganisms (Jousset et al., 2017). As shown in Figs. 4–5, most sub-modules were composed mainly of rare taxa, showing the key role of rare taxa in community assembly. Correspondingly, variations in rare taxa also may lead to drastic changes in entire networks and sub-modules (Guimerà and Nunes Amaral, 2005). Degree-based keystone species in co-occurrence networks are often regarded as species that play an irreplaceable role in the microbial communities (Hu et al., 2017). The disappearance or change of these keystone species can lead to the disintegration of modules and networks, even sub-network (Shi et al., 2016). In this study, there were 24 keystone species belonged to rare taxa, which indicated that rare taxa may have a great impact on other communities (Wang et al., 2020). Thus, we can conclude that rare taxa were more important than abundant taxa in maintaining the composition and structure of entire microbial community.

Nitrogen and phosphorus are key elements for the synthesis of protein and nucleic acid, and they play an important role in the growth, development, and division of cells. Many bacteria populations were significantly correlated with inorganic nitrogen concentrations, such as Salinivibrio and Sphingopyxis of Proteobacteria, Chryseobacterium, Mycobacterium, indicating these bacteria groups have potential competition to respond to or use the variations of inorganic nitrogen in water. This is similar to those Proteobacteria populations that could actively uptake nitrate in coastal surface waters (Morando and Capone, 2016). Besides, the relative abundance of several bacteria populations was significantly related to organic carbon concentrations, including Enterobacter, Aeromonas, and Acinetobacter belonging to Proteobacteria. These groups could respond to the organic carbon availability and utilize them as a carbon source. In addition, Polaromonas was easier to grow in warmer water.

Different nitrogen and phosphorus concentrations also had a greater impact on the growth and physiological and biochemical reactions of algae. It was observed that Nannochloris and Pseudoschroederia had a positive correlation with nitrate, DOC, and phosphate. However, Shen et al. (2015) showed that compared with organic carbon sources, supply of carbon dioxide under suitable nitrogen and phosphorus nutrient conditions increased the growth rate of microalgae.

4.2 Network associations between bacterial communities and phytoplankton taxa

Although co-occurrence networks can predict the direct associations between taxa, microorganisms may also show positive or negative correlation based on indirect reasons such as environmental preferences (Weiss et al., 2016). The extent to which closely related microbial groups may coexist is still a topic worthy of discussion (Mayfield and Levine, 2010). Previous studies on aquatic microbial communities have demonstrated that closely related taxa have consistent temporal dynamics and similar niches (Andersson et al., 2010).

The Mychonastes are distributed worldwide and have a wide range of habitats, including streams and large still waters, for example, Baringo and Victoria Lakes in Kenya, Erken Lake in Sweden, and Stechlin Lake in Germany, etc. (Krienitz et al., 2011). However, as the most abundant genus in spring, summer, and winter in our study, Mychonastes was first reported by Li et al. (2013) as a new record genus of freshwater Chlorophyceae isolated from Dianchi Lake. Mychonastes is a taxon that can naturally adapt to the desert climate as its optimal temperature is 30 ℃, which is significantly related to its abundance. Synechococcus is a common prokaryotic picocyanobacteria, which had a relatively high abundance in summer. Liu et al. (2019) studied the growth and interspecies competition of Mychonastes and Synechococcus under different N: P ratios. Although we have not found the competitive relationship between Mychonastes and Synechococcus, Mychonastes had a significant correlation with other cyanobacteria (Cyanobacterium) in autumn. Flavobacterium is considered as a bacterium that can promote plant growth, which was reported to exist during the algal blooms and can improve the biological phosphorus removal ability, and participate in the cleavage of polysaccharides (Park et al., 2007). In our study, based on positive correlation, Flavobacterium probably promoted the growth of most green algae. The node of green algae Radiococcus had the highest connectivity in autumn. The stability of ecosystem function lies on the interaction between species (Zhou et al., 2010). Lower connectivity in the community will lead to higher functional stability of the system (for example, a non-scale network), because the entire network module is less affected by a node loss (Yang et al., 2017).

Cyanobacteria are the only prokaryotes with oxygen-producing photosynthesis. They exist in different niches and are important participants in the global carbon and nitrogen cycle. The planktonic cyanobacterium Synechococcus is ubiquitous in the ocean and fresh waters and plays an important role in total carbon sequestration on a global scale in oligotrophic aquatic environments (Callieri et al., 2012). We observed that Synechococcus were negatively related to other bacteria, the negative correlations in co-occurrence patterns might imply the predation or competition among taxa. The presence of correlations between populations increases the probability that these populations interact. However, the relationship of these populations may change when different populations are influenced by the same factor. Therefore, not all algae-bacteria connections can be identified by correlation methods, and correlations may not necessarily come from the connections between populations (Paver et al., 2013), which may also be indirectly affected by the external environment.

Numerous pieces of evidence have been demonstrated that a link occurs between diatoms and microorganisms (Behringer et al., 2018); however, whether they cooperate or exclude each other still lacks relevant knowledge. In our results, we found significant correlations between diatom communities and bacteria communities. In addition, these relationships were greatly subjected to seasonal variations. There is considerable evidence that the combination of diatoms and bacteria is beneficial to the growth of algae (Grossart et al., 2005) and bacteria (Stanish et al., 2013). The symbiotic relationships between Bacillariophyta and phylogenetically related bacteria verified that cell yield and growth rate of bacteria increased when grown in mixed culture, because the organic material secreted by diatoms was absorbed by and would support the growth of bacteria (Jolley and Jones, 1977).

The growth stages of plankton throughout the year include mainly: winter (transition stage, December–January), spring (accumulation stage: March–April; climax transition stage: May), summer (water bloom outbreak stage: June–August, depletion stage: September–November). Different microorganisms may play the same ecological role over time, which may explain the unique keystone taxa detected in different ecological networks. Therefore, our results provide evidence for the profound influence of physico-chemical conditions on the seasonal dynamics of microbial community composition.

5 CONCLUSION

The eukaryotic phytoplankton and bacterial networks constructed in our research conformed to the properties of having non-randomly connections, modular structure, and "small world" In co-occurrence patterns, the positive interactions between bacteria and eukaryotic phytoplankton ASVs may be evidence of mutual influence (cooperation), while the negative correlations may indicate predation or competition. Through its within-module connectivity (Z_i) and among-module connectivity (P_i), we identified the keystone species that can maintain the stability of microbial networks. Most of the keystone species were defined as rare taxa, a few abundant taxa were still observed although rare taxa accounted for a large proportion. In the bacterial networks, ASVs of Proteobacteria (Aeromonas, Devosia, Salinivibrio, Sphingopyxis, and Acinetobacter) and Acidobacteria (Mycobacterium) can play an important role in the stability of ecosystem; and in eukaryotic phytoplankton networks, ASVs of Chlorophyta (Mychonastes, Pseudoschroederia, Follicularia, and Monoraphidium), Bacillariophyta (Cyclotella), and Cryptophyta (Cryptomonas) were significantly related to water temperature and inorganic nitrogen. With the disappearance of the "keystone species" in the network center, the aggregation structure will suffer major damage, which may bring huge biogeochemical consequences. The abundance of these keystone species was significantly related to many ecosystem functions and environmental changes, highlighting the key ecological role of rare species, and providing a new ecological significance for the seasonal co-occurrence patterns of bacteria and eukaryotic phytoplankton in urban aquatic ecosystems.

6 DATA AVAILABILITY STATEMENT

The datasets generated and/or analyzed during the current study are available in the Sequence Read Archive (SRA) of NCBI repository, under accession number SRP273360.