Samples were collected in June (summer) and December (winter) 2020 during the day (10 AM to 17 PM) from the Songtao Reservoir of Hainan province, China (Figure 1). In this study, a total of 9 sampling sites were selected, namely BS (upstream, S7-S9), NF (midstream, S1-S3) and FJ (downstream, S4-S6). We used the sterile sample glass bottles (1 L) to collect water samples at a depth of 10–20 cm below the water surface. Six samples were collected at each sampling point (every samples were obtained at intervals of 3 m) and 1 L of water was filtered per sample. The water samples of microbial were transported to a laboratory at 4°C for filtration.

pH, EC (electrical conductivity, mS cm⁻¹), DO (dissolved oxygen, mg L⁻¹), Tem (temperature, °C), TDS (total dissolved solids, g L⁻¹), ORP (oxidation-reduction potential, mV), and Tur (turbidity, NTU) were determined by multi-purpose meters (YSI and turbidity meter) in the field. TN (total nitrogen, mg L⁻¹), TP (total phosphorus, mg L⁻¹), and COD (permanganate index, mg L⁻¹) were determined by spectrophotometric methods (The State Bureau of Quality and Technical Supervision, 1989a,b,c). All chemical analyses were conducted in triplicates. The water samples used were transported to the laboratory at 4°C. Prior to analysis, the water samples were filtered through 0.45 μm porous polycarbonate membranes (Collins, Shanghai, China).

The metal concentrations (Cr, Mn, Ni, Cu, Se, Cd, and Ba) in water were determined in this study by ICP-MS (Thermo scientific iCAP RQ). Firstly, the water samples were filtered through a 0.45 μm membrane. And then nitric acid was added until its concentration reached 2%. The internal standard method was used for the calibration method (Sc, Ge, In, Bi as internal standard elements), and the recovery rate was 98% to 109.8%. The concentration of metal elements in water was expressed by μg L⁻¹.

The water samples were filtered using a 0.22 μm membrane and stored at −80°C. All microbial community analysis was conducted in sextuplicates. The CTAB method was used to extract the total genome DNA from samples. A total of 1% agarose gels were used to monitor DNA concentration and purity. And the sterile water was used to dilute DNA to 1 ng μL⁻¹. We used the specific primer 341F(5′-CCTAYGGGRBGCASCAG-3′) and 806R(5′-GGACTACNNGGGTATCTAAT-3′) to amplified 16S rRNA genes of distinct regions (16S V3-V4). TruSeq^® DNA PCR-Free Sample Preparation Kit (Illumina, USA) was used to generate the sequencing libraries and index codes were added. Then, the Qubit@ 2.0 Fluorometer (Thermo Scientific) was used to assess the library quality. Finally, the library was sequenced on an Illumina NovaSeq platform and 250 bp paired-end reads were generated. The 16S rRNA sequence is deposited at NCBI with the BioProject ID PRJNA933804 and SRA ID SRP422163.

To study the feature table of the operational taxonomic unit (OTU), the initial step involved converting raw data FASTQ files into a format compatible with the QIIME2 system using the qiime tools import program. Subsequently, sequences from each sample were demultiplexed, underwent quality filtering, trimming, and de-noising. Following this, the sequences were merged, and chimeric sequences were identified and eliminated using the QIIME2 dada2 plugin (Behera et al., 2021). We used the classification criterion of similarity ≥ 99% to divide the OTU of (similarity of species level by default) (Bokulich et al., 2018). We used the OTUs, Chao1 index, Shannon index, and Simpson index to assess the microbial diversity within an individual sample. The principal coordinate analysis (PCoA) and analysis of similarities (ANOSIM) were used to visualize the structural variation of microbial community across samples (beta diversity) (Vázquez-Baeza et al., 2013). The redundancy analysis (RDA) and mantel's test were used to assess environmental factors driving microbial community composition. All data were log (X+1) transformed before being analyzed.

Present partial least squares structural equation modeling (PLS-SEM) was used to evaluate the direct and indirect relationships between microbial OTUs and season (summer and winter), region (upstream, midstream, and downstream), physicochemical factors (Tem, DO, C, TDS, pH, ORP, Tur), nutrients (TN, TP, COD), and metals (Cr, Mn, Ni, Cu, Se, Cd, and Ba). The season and region were given assignments (1: summer, 2: winter, and 1: upstream, 2: midstream, 3: downstream, respectively). SmartPls 4.0 was used to construct the model. Firstly, we constructed a structural model consisting of six main latent variables, namely Season, Region, Physicochemical factors, Nutrients, Metals, and Microbial, and established relationships among the latent variables. Secondly, all observed variables were then imported and the model was fitted. Finally, the model was optimized by eliminating the indicators with low reliability and validity through the fitting results. The results of reliability and validity was shown in Supplementary Table S3. The detailed information on the constructive model referenced by Wen and Li (2019).

One-way ANOVA was used for the inter-group difference with P < 0.05 of significant difference. The significant difference in data was detected by IBM SPSS statistics 20. The R language program and GraphPad Prism 8.0 were used to draw heat maps, bar plots, and other figures. The figure of sampling sites was drawn by Arcgis 10.7 and Google Earth.

According to the sequencing results, there were significant variation in the microbial community across different periods and regions (Figure 2). A total of 24,686 OTUs were identified in all samples. Notably, distinct OTUs were identified during summer (9,391, 38.04%) and winter (14,250, 57.73%), with only a minimal overlap (1,045, 4.23%) between the two seasons (Figure 2B). Additionally, variations in OTUs number were observed in different regions (Figure 2A). The largest number of OTUs was recorded in NF (11,015, 44.62%), followed by FJ (9,250, 37.47%), and the lowest in BS (8,970, 36.34%). The number of OTUs shared by the three regions was 1,559 (6.32%). The unique OTUs number of NF, FJ, and BS were 8,680 (35.16%), 6,624 (6.83%), and 6,392 (25.89%), respectively.

As shown in Figures 3A2, A3, A4, the Chao1 index, Shannon index, and Simpson index in winter exhibited significantly higher values compared to summer (P < 0.01). Regarding spatial distribution characteristics, while the number of OTUs and the Chao1 index remained relatively stable, while the alpha diversity indexes in BS were significantly higher than the other two regions (P < 0.001, Figures 3B1, B2, B3, B4).

In Figure 4, we used the Bray-Curtis distance to study the similarity of all samples in PCoA. It is possible to calculate the quantitative characteristics of the different species composition of a biological sample species by the Bray-Curtis distance. The microbial community composition in the data was explained by PC1 (50.31%) and PC2 (12.30%) (Figure 4). Notably, the distribution of microbial communities at the same location varied significantly across different seasons. Meanwhile, there was a significant spatial difference between BS and FJ in the microbial community composition distribution. Overall, the seasonal variation in microbial community was larger than the spatial variation. ANOSIM results (Supplementary Figure S1) further confirmed that there were significant spatial and temporal differences in the microbial composition of the Songtao Reservoir.

The environmental factors exhibited significant variations over time, with pH, EC, DO, Tem, Ni, Se, and Cd registering significantly higher levels during summer compared to winter. Conversely, ORP, TN, TP, and Mn exhibited significantly higher levels during winter than summer. Furthermore, in terms of spatial distribution, the BS region exhibited peak values solely for TN, TP, COD, and TDS (Supplementary Table S1). Moreover, Pearson's correlation showed that there were correlations between environmental factors except for the Ba, TN, and Tur (Figure 8). Significant correlations between the Tem and the Se and Ni were found, especially. Pearson's correlation coefficients were < 0.90 among Se, Ni, and Tem; among Cu, Se, and Ni and between pH and Se.

In this study, RDA was used to study the impact of environmental factors (pH, Tem, DO, TN, TP, Cr, Mn, etc, Supplementary Table S2) on microbial community composition (Figure 9). The explanation amount of RDA1 and RDA2 were 61.70% and 21.95%, respectively. The highest correlation coefficient with RDA1 was TP (r = 0.9999), followed by TN (r = 0.9949), DO (r = −0.9933), Cd (r = −0.9845), Se (r = −0.9206), Tem (r = −0.9135), Ni (r = −0.8931), Cr (r = –0.8669), pH (r = −0.8567), Mn (r = 0.8496), Cu (r = −0.8474), and ORP (r = 0.8244). Meanwhile, Tur was correlated with RDA2 (r = 0.9951), followed by TDS (r = 0.8786), Ba (r = 0.8636), and EC (r = 0.8377, Supplementary Table S2).

The RDA results showed the environmental factors correlated with the seven dominants microbial at the phylum level including TP, Se, Tem, etc (Figure 9). Specifically, Actinobacteria was positively correlated with TN, TP, and Mn. Proteobacteria was negatively correlated with most environmental factors except DO and Cd. In addition, there was a positively correlated between Cyanobacteria, Bacteroidetes, and Verrucomicrobia with EC, TDS, Tur, and Ba (Figure 9). The mantel's test result (Figure 8) showed that the environmental factors include Tem, ORP, TP, and Cd were significantly correlated with the relative abundance of archaea and bacteria community (P < 0.05). And, except for the Ba, all environmental factors were highly significantly correlated with the bacteria (P < 0.01).

The results of the structural model (Figure 10) revealed that the factor with the highest influence on the OTUs was the metals (path coefficient = −0.959, R² = 0.843, P < 0.001), followed by the physicochemical factors (path coefficient = 0.812, R² = 0.935, P < 0.01), and the lowest was nutrients (path coefficient = 0.292, R² = 0.621, P < 0.01). In addition, there were significant (P < 0.001) path coefficients for spatial and physicochemical factors, nutrients, and metals. Although the direct effect of season on the OTUs was not significant, we found that season changed the OTUs by affecting physicochemical factors, metals, and nutrients. Spatial variation had the same significant indirect effects on the OTUs (Supplementary Table S4).

Our results indicated that Proteobacteria, Actinobacteria, Cyanobacteria, and Bacteroidets accounted for 90% of the total microbial community at the phylum level (Figure 5). This may be related to Proteobacteria, Actinobacteria and Bacteroidetes are all believed to degrade complex organic macromolecules (González et al., 2000; Kirchman, 2002; Zaitlin and Watson, 2006). Previous study showed that Proteobacteria and Actinobacteria were the dominant microbial community in surface water of reservoir which account for more than 50% of the total bacteria (Yu et al., 2014). Earlier research also found that Proteobacteria was a general component in freshwater microbial community, accounting for 60%−70% of all microbial community in rivers, lakes, and reservoir ecosystems (Newton et al., 2011; Parfenova et al., 2013). Similarly, Actinobacteria is also a relatively large group in freshwater (Newton et al., 2011; Ji et al., 2019; Niu et al., 2019). Some studies, for instance, found that Actinobacteria existed in almost all lakes, reservoirs, and rivers, which is related to its involvement in the decomposition of organic matter (Parfenova et al., 2013; Ghai et al., 2014). Bacteroidetes were the most diverse group of bacteria living in mesotrophic and eutrophic water environment (Ligi et al., 2014). Therefore, it is essential for Actinobacteria, Proteobacteria, and Bacteroidetes in determining the microbial composition of water bodies, although the main microbial community differs in freshwater (Jiang et al., 2021). Our results showed that the dominant populations of all samples were similar to other rivers, lakes, and reservoirs in China (Yu et al., 2014; Jiang et al., 2021). However, the composition of microbial communities in this study differed from that of eutrophic rivers and reservoirs, which the dominant microbial were Actinobacteria, Proteobacteria and Bacteroidets (Debroas et al., 2009).

In this study, the relative abundance of Cyanobacteria was high in Songtao Reservoir (Figure 5) than in the subtropical region (Guo et al., 2021; Zhu and An, 2021). Previous studies showed that Cyanobacteria are the simplest and most primitive green autotrophic bacteria, widely distributed in freshwater. Under the influence of global climate change, the abundance of Cyanobacteria in lakes is gradually expanding (Hou et al., 2022). Lürling et al. (2018) found that the Cyanobacteria had a competitive advantage at elevated temperatures. This study was sampled in a tropical region, which may explain the relatively high abundance of cyanobacteria. And, excessive Cyanobacteria are the main feature of blooms (Brooks et al., 2016). Therefore, more attention should be paid to monitoring and controlling Cyanobacteria in the water environment in summer in Songtao Reservoir.

The results indicated that there were significant temporal and spatial differences in microbial community composition across Songtao Reservoir (Figure 6). Specifically, both diversity (Figures 3, 4) and OTUs number (Figure 2) of microbial community showed significant differences among the two seasons and three regions, the PLS-SEM result also indicated that season and region have significant indirect and direct effects on microbial community composition (Figure 10).