The study area, Kaiyang, is situated in Guizhou Province, southwest China, in a typical subtropical monsoon climate zone with an annual average temperature of 14.4°C and precipitation primarily concentrated from June to September. There were two karst rivers from west to east, the south river and the north river, with three underground entrances (K4, K8, K9) and three underground exits (K5, K7, and KX) (Figure 1). A total 40 samples were collected. Each 10 samples were collected in January 2017 (Jan-2017), August 2017 (Aug-2017) (Xiang et al., 2020), August 2018 (Aug-2018) and January 2019 (Jan-2019) (Zhang et al., 2021), respectively.

Except for K9 and KX, all sampling sites were in the south river. K1 to K4 were located in the surface river from upstream; K1 was located in a reservoir; K2 and K3 were situated upstream and downstream of a waste water treatment plant (WWTP), respectively; and K4 was situated in a resort. There were agricultural fields and a small amount of livestock activity between K5 and K8. K5 was an outlet of underground river, with evident domestic wastes. K8 was an underground entrance. K6 and K7 were located in the upstream tributaries of K5 and K8, respectively. K6 was a surface river site with domestic wastes, and K7 was the outlet of an underground river. K9 was a sinkhole in the scenic region of the north river, while KX was the confluence of the north and south rivers.

Temperature (T) and pH value of karst river water samples were determined by pH meter (pH30, CLEAN, CA, United States). Electronic conductivity (EC) and dissolved oxygen (DO) were determined by EC meter (CON, CLEAN, CA, United States) and DO meter (CON30, CLEAN, CA, United States), respectively. The total organic carbon (TOC) was determined by TOC analyzer (Shimadzu, Japan). HCO₃⁻ was tested by potentiometric titrator (Metrohm 877 Titrino plus, Swiss), and Cl⁻, SO₄²⁻ were measured by ion chromatograph (DIONEX ICS-900, Sunnyvale, CA, United States). Ca²⁺, Mg²⁺, K⁺, and Na⁺ were measured by ICP-OES (SPECTRO Blue Sop, Germany). And the spectrophotometer UV-1800 (Shimadzu, Japan) was used to determine dissolved inorganic nitrogen (DIN) - NH₄⁺, NO₃⁻ and NO₂⁻. Following our previous studies (Huang H. et al., 2019), auto-solid phase extraction (SPE) instrument (Auto SPE-06C, Reeko Instrument, TX, United States) with an Oasis HLB SPE column (6 mL, 500 mg, Waters, MA, United States) was used to extract antibiotics from water samples of karst river. Samples were collected and immediately chilled in 4°C and transported via overnight express to the laboratory for further test. The sampling process was followed our previous study (Xiang et al., 2020).

The bacteria in each sample were filtered through a 0.22 μm Millipore GSWP membrane for 1.0 l water and then suspended in 10 mL of physiological saline for 6 h. The mentioned suspensions were separated and kept for DNA extraction using a PowerSoil DNA Extraction Kit (MoBio Laboratories, Carlsbad, CA, United States). The extracted DNA was stored at −20°C for further analysis and at −80°C for permanent preservation. Using a NanoDrop ND-2000 spectrophotometer (Thermo Fisher Scientific, Wilmington, DE, United States) and agarose gel electrophoresis (Bio-Rad, Hercules, CA, United States), the quantity and quality of the extracted DNA were established.

The primers 338F and 806R (338F: 5′-barcode-ACTCCTAC GGGAGGCAGCAG-3′ and 806R: 5′-barcode-GGACTACHVGG GTWTCTAAT-3′) were used to amplify the V3-V4 regions of the bacterial 16S ribosomal RNA gene (16S rRNA gene). The PCR amplification method was referred to our previous study (Zhang et al., 2020). Raw fastq files were demultiplexed, then quality-filtered using QIIME (Version 1.17). All clusters were subsampled by the minimum reads (19,551 reads), and operational taxonomic units (OTUs) were clustered with a 97% similarity cutoff by Uparse.¹ Mothur² was used to compute the alpha diversity. The Kyoto Encyclopedia of Genes and Genomes (KEGG) functional profiling was predicted using PICRUSt2.³ The clustered sequence datasets in this study were deposited in the NCBI Genbank with the accession number PRJNA908735 and PRJNA936825.

The data were normalized except for pH value. Spearman correlation analysis and its significance test among each factor were completed by “corrplot” package.⁴ The ordinary least squares (OLS) analysis among ammonium, nitrate, sulfate, and other environmental factors were calculated by “car” package.⁵ Principal co-ordinates analysis (PCoA) based on Bray-Curtis distance and variance partitioning analysis (VPA) were both performed by “vegan” package.⁶ Venn diagram at genus level was plotted by “vennDiagram” package. The Kruskal-Wallis test (p < 0.05) was completed by “nparcomp” package, and corrected p-values were calculated using Dunn’s test. The microbial community assembly processes were analyzed by “iCAMP” package.⁷ And the heatmap figure was plotted by “pheatmap” package.⁸ All the above data analysis were implemented in R.⁹ The co-occurrence network of bacteria at genus level was visualized with Gephi 0.9.2,¹⁰ and ecological modules were identified using the Louvain algorithm. The saturation index of gypsum (SIg) and saturation index of calcite (SIc) were calculated by PHREEQC (Version 3.6.2). Structural equation model (SEM) was performed by SPSS AMOS (Version 24.0) to discuss the environmental factors effect on ecological module and element metabolic module.

The pH value of karst river water was ranging from 7.70 to 8.82, and the temperature was ranging from 4.0°C to 26.1°C (Supplementary Table 1). DO ranged from 1.88 to 9.93 mg/L, and TOC ranged from 1.09 to 30.63 mg/L. Ca²⁺ and HCO₃⁻ were the major cation and anion due to the abundant carbonate rocks in the study area (Sun et al., 2019). For carbonate rock weathering involving only carbonic acid, the equivalent ratio of [Ca²⁺+Mg²⁺]/[HCO₃⁻] is typically around 1. But the average value of our study was 1.23, while the mean value of [Ca²⁺+Mg²⁺]/[HCO₃⁻ + NO₃⁻ + SO₄²⁻] was 0.93, indicating that the carbonate weathering was influenced by both carbonic acid and sulfuric/nitric acids (Martin, 2017; Sun et al., 2021).

The concentration of NO₃⁻ varied from 0.74 to 45.82 mg/L, with an average value of 18.23 mg/L, and the concentration of SO₄²⁻ varied from 20.14 to 74.25 mg/L. There were severe NH₄⁺ and NO₂⁻ pollution, with NH₄⁺ amounting from 0.02 to 45.35 mg/L and NO₂⁻ amounting from 0.01 to 9.51 mg/L. The concentrations of both NH₄⁺ and NO₃⁻ showed seasonal variations with higher in winter and lower in summer, mainly influenced by the dilution effect of precipitation recharge in rainy season. Meanwhile, the NH₄⁺ were higher in upstream, particularly from K1 to K4, while NO₃⁻ concentrated downstream from K7 to KX, due to nitrification and the NO₃⁻ accumulation for its high solubility (Li X. et al., 2021; Li H. et al., 2021). The negative correlations between NH₄⁺ and NO₃⁻ also supported this hypothesis (Eqs 1, 2):

As shown in Eq. 1, The concentration of NH₄⁺ was positively correlated with antibiotics, as well as its correlation coefficients with macrolides, lincomycin, tetracyclines were 0.78, 0.71 and 0.66 (p < 0.01), respectively. The concentrations of NH₄⁺ and antibiotics showed a similar spatial variation trend, displayed higher upstream. As one of the most prevalent forms of nitrogen pollution in WWTP, NH₄⁺ originated primarily from wastewater plant discharges, the same source of antibiotics.

SO₄²⁻ is mainly derived from both natural (dissolution of soluble sulfate, oxidation of sulfide minerals, and atmospheric precipitation) and anthropogenic inputs(agricultural fertilizers, domestic sewage, industrial wastewater, and mine waste water)(Li et al., 2018). In the study area, the [Mg²⁺]/[Ca²⁺] values ranged from 0.39 to 1.14 with average 0.65, which indicated that the chemical composition is mainly influenced by low-magnesium minerals such as calcite and gypsum. Furthermore, the SIg values were − 2.43 ~ −1.75 while SIc were − 0.79 ~ 1.28, implied that the dissolution of gypsum were the main source of sulfate. Additional, SO₄²⁻ were positively correlated with antibiotics (R = 0.66, p<0.01), which implied that SO₄²⁻ also originated from WWTP, which was consistent with its concentration spatial distribution of higher in upstream.

A total of 28 antibiotics were detected, with total concentrations ranging from 9.13 to 1411.33 ng/L (Supplementary Table 2). The mean concentrations of four sampling times from high to low were Jan-2019 (average value of 514.15 ng/L), Jan-2017 (472.60 ng/L), Aug-2017 (392.44 ng/L), and Aug-2018 (199.37 ng/L), with a seasonal trend of higher in winter and lower in summer (Figure 2), due to the dilution effect of precipitation recharge, which was consistent with our previous study (Huang F. et al., 2019; Huang H. et al., 2019). The maximum value for each time was frequently recorded at K2 and K3, suggesting that WWTP effluent discharge was the major source of antibiotics in karst river, which indicated the anthropogenic inputs. Antibiotic demonstrated a spatial changes that fluctuates and decreases along the stream, which was affected by the adsorption of soil and sediment and the degradation of bacteria (Zhao et al., 2012). Moreover, antibiotic concentrations in tributaries K6 and K7 were lower than those in the mainstream.

Quinolones (QNs), lincomycin (LIN), and macrolides (MLs) were the major antibiotics, accounting on average for 43.45, 34.58, and 13.99% of total concentration, respectively, which were consistent with their frequent use in southwest China (Sun et al., 2015; Zeng et al., 2019). QNs ranged from 3.10 to 994.13 ng/L (Figure 2A) with a 95% detection frequency. QNs were persistent pollutants due to their high chemical stability and prolonged half-life (Andriole, 1999). Ofloxacin and norfloxacin, as the maximum content subclass of QNs in most WWTP (Dinh et al., 2017), also dominated in karst river, ranging from 4.60 to 486.00 ng/L and 4.09 to 204.00 ng/L, respectively. The detection frequency of LIN was 95%, and the concentration was 4.77–227.53 ng/L (Figure 2B). The detection frequency of MLs was 88%, with concentration ranging from 0.83 to 205.13 ng/L (Figure 2C). Erythromycin and roxithromycin were the main subclasses of MLs. While antibiotics were transferred underground by continuous precipitation recharge in summer, their natural degradation was weakened due to weak light and then accumulated (Liu et al., 2021). However, the high connectivity of karst aquifers enabled rapid hydration and a more quickly response to contaminants.

The concentrations of tetracyclines (TCs), chloramphenicols (CAP), and sulfonamides (SAs) were ranged from 3.43 to 177.80 ng/L (Figure 2D), 1.70 to 76.20 ng/L (Figure 2E), and 0.65 to 33.56 ng/L (Figure 2F), respectively. The detection rate of TCs was 60%, with tetracycline and oxytetracycline being the most commonly used antibiotics. The detection frequency of CAP was 53%. Both TCs and CAP were higher in winter. The detection rate of SAs was 75%, with sulfapyridine being the most common subclass. The maximum content and the most subclasses of SAs were observed in Aug-2018. Only samples from Aug-2018 contained sulfadimethoxine and sulfadimidin, both of which were widely applied in veterinary medicine. Increased livestock vaccination in study area in 2018 may result in an increase use of veterinary drugs.¹¹ As a notable veterinary drug, enrofloxacin was also detected only in Aug-2018.

A total of 3,955,799 optimum sequences were obtained in 40 samples (Supplementary Table 3). After subsampling by the minimum sequences, 5,242 OTUs were obtained. The species richness in Aug-2017 was the lowest. Nonetheless, the Shannon index had no significant difference at four sampling times (Supplementary Table 3).

Proteobacteria (relative abundance ranged from 36.03 to 61.15%), Bacteroidetes (18.97 to 38.03%), Cyanobacteria (6.14 to 13.85%), and five other phyla with relative abundance greater than 5% were the dominant phyla. Flavobacterium, Arcobacter, norank_c_Cyanobacteria, Pseudomonas and Acinetobacter were prevalent in a total 945 genera, with the relative abundance of 12.68, 5.78, 5.18, 4.41 and 4.35%, respectively.

PC1 and PC2 explained 22.65% of the variance in community structure together (Figure 3A). The bacterial communities of Jan-2017 and Jan-2019 were similar. The bacterial community structures of Aug-2017 and Aug-2018 were divided from each other, which apparently differed from samples taken in winter.

The relative abundance of 275 shared genera varied from 89.04 to 96.27% (Figure 3B). Most of them showed obvious seasonal characteristics. Seven genera were mainly distributed in winter, including Flavobacterium, Pseudomonas and hgcI_clade. In previous study, we observed strong correlations between ARGs and Flavobacterium, Pseudomonas (Zhang et al., 2021), which were more likely to survive in winter with higher antibiotic concentrations. Six genera were enriched in summer, including Dechloromonas, Hydrogenophaga and Cloacibacterium. The relative abundance of unclassified_f__Comamonadaceae gradually increased while Limnohabitans showed no noticeable trend in this study. Nine out of the nineteen dominant genera (relative abundance above 1%) (Supplementary Figure 1) exhibited statistically significant differences (p < 0.05) (Figure 3C). Flavobacterium, Pseudomonas, hgcI_clade and Albirhodobacter were considerably enriched in winter with average relative abundances varying from 8.15 to 28.06%, 6.26% to 8.08, 3.53 to 4.54% and 1.49 to 1.62%, respectively. Acinetobacter, Malikia, Dechloromonas and Hydrogenophaga were highly expressed in Aug-2017 with average relative abundances of 12.37, 11.90, 5.24 and 1.94%, respectively. And the relative abundance of norank_c_Cyanobacteria (7.15%) was the highest in Aug-2018.

Compared with shared genera, the richness and abundance of specific genera were lower. Forty seven genera only presented in winter, with relative abundances of 0.23 and 0.24% in Jan-2017 and Jan-2019, including Candidatus_Azambacteria, Empedobacter, etc. The winter genera were mainly distributed at K4 and K6 in Jan-2017 and at K6 in Jan-2019. As a sampling site in the tributary, there might be a relatively higher flow at K6 that introduce a great number of exotic species in winter. Forty six genera were only found in Aug-2017, with an average relative abundance of 1.13%, and distributed mainly at K2 and K3, such as Macromonas, Candidatus_Aquiluna. Forty nine genera such as Sporacetigenium and Desulfitobacterium were only found in Aug-2018, with an average relative abundance of 1.73% and primarily enriched at K8.

Five clusters (MOD0-MOD4) were obtained through co-occurrence network analysis of genera with relative abundance above 0.1% (Figure 4A). Twenty one genera were included in MOD0, primarily Flavobacterium, Malikia, Dechloromonas, etc., each of which exhibited distinct seasonal differences. Thus, there was no discernible seasonal trend in MOD0 but a peak value in Jan-2019 (Figure 4C), where Flavobacterium concentrated. Twenty three genera of MOD1 were mainly included Limnohabitans and Planktothrix, which concentrated in summer. Twenty five genera were included in MOD2, such as Arcobactor, Pseudomonas, Cyanobacteria. MOD3 was consisted of 17 genera, such as hgcI_clade and Hydrogenophaga, etc., which were abundant in winter. Prosthecobacter, Nitromonas and Nitrospira were dominant in MOD4 among the total 24 genera.

The co-occurrence network of major microorganisms helped to reveal functional group structure and potential interactions between bacteria (Curtis et al., 2002; Cardona et al., 2016). There were 165 positive and 68 negative edges among five MODs (Figure 4B), with a modularity of 0.41, suggesting a strongly associated community among the major communities which were mainly positive synergies. MOD0, MOD3 and MOD4 were strongly connected with other MODs with 100, 98 and 65 positive edges, respectively. There were more negative edges between MOD1 and MOD4, indicating a negatively competitive relationship among genera. Genera involved in carbon, nitrogen and sulfur metabolism were widely distributed in each module. Flavobacterium, Dechloromonas, Massilia in MOD0 participated in denitrification, nitrogen reduction and sulfate reduction (Tanikawa et al., 2018; Chang et al., 2019). The diversity of metabolic function makes it most closely related to other MODs. MOD3 were primarily denitrifying bacteria, such as hgcI_clade (Herrmann et al., 2015; Luo et al., 2018). Conversely, ecological modules with different functions often exhibit a tendency to compete for limiting nutrients (Tate, 2020). This competitive relationship may also exist between MOD1 and MOD4. MOD1 was dominated by Limnohabitans, Planktothrix, and other genera with both carbon and nitrogen metabolism function, whereas the MOD4 was ruled by nitrifying bacteria such as Prosthecobacter, Nitrosomonas, and Nitrospira. In addition, Hydrogenophaga and Limnohabitans in MOD1 and MOD4 participate in sulfur metabolism (Herrmann et al., 2015; Luo et al., 2018). Conjecturally, bacterial communities participated in nitrogen and sulfur cycles were potentially affected by the pollution of antibiotics and nitrogen in karst river.

The bacterial community assembly processes were consisted of deterministic processes (heterogeneous selection and homogeneous selection) and stochastic processes (dispersal limitation, homogenizing dispersal and drift) (Stegen et al., 2013). Deterministic processes represented the influence of environmental factors such as physical and chemical indicators on microbial community, while stochastic processes were related to spatial variation. In general, deterministic processes were dominant in small spatial scale while stochastic processes dominated in a larger spatial scale (Dini-Andreote et al., 2015). However, frequent anthropogenic impacts and intense variation in karst hydrodynamic conditions produced high variance and increased the influence of stochastic processes.

In Figure 5, stochastic processes had more advantage in Aug-2017 (56.96%) and Aug-2018 (50.66%), and primarily affected by dispersal limitation and drift effects. Stochastic processes played a major role in microbial community assembly with less environmental pressure (Dini-Andreote et al., 2015). Lower antibiotic concentrations in summer reduced the survival pressure for less resistant bacteria, which increased their relative abundance. The type of antibiotics in Aug-2017 and Aug-2018 were remarkable different, and their additive, synergistic and antagonistic effects led to greater differences in community structure (De Liguoro et al., 2018; Liang et al., 2022). The rapid surface-underground exchange and the drastic environmental changes in karst river resulted in short adaptation and rapid turnover in summer, which benefited distinguishable bacterial communities as well. Deterministic processes enhanced in winter, and homogeneous selection was predominant, with 49.51 and 68.14% in Jan-2017 and Jan-2019, respectively, which was caused by the survival pressure of high concentration of antibiotics and DIN. Therefore, we suggested that the environmental selection caused by the multiple pressure of antibiotics and DIN created a more stable bacterial community with antibiotic resistance in winter. Consequently, despite the two-year sampling interval between Jan-2017 and Jan-2019, the structure of bacterial community showed a highly similarity, while the stronger survival pressure from pollutions led to a decline in the relative abundance of specific genera. The proportion of dispersal limitation in the stochastic processes decreased gradually, while the proportion of homogenizing dispersal increased, indicating that the composition of microbial community tended to be more similar over time.

Variations in microbial community were accompanied by alterations in community function (Zhong et al., 2018). The results of KEGG functional profiling revealed that the average relative abundance of 16 functional genes exceeded 1%, which were mainly sub-functions of metabolism, environmental information processing, and genetic information processing and so on (Figure 6). Among them, the average relative abundance of energy metabolism ranged from 5.58 to 5.79% and was higher in summer, possibly due to nutrient limitation (Cline and Zak, 2015; Jiang et al., 2022). In winter, the copious nutrients conditions made nutrients more accessible to microorganisms, and hence the corresponding metabolism gene was less prevalent. However, nutrient availability was restricted by dilution effects in summer, and microorganisms maintained their growth by increasing the metabolic function gene (Jiang et al., 2022).

From the results of function profiling, 13 carbon fixation pathways, 5 nitrogen metabolism pathways, and 2 sulfur metabolism pathways were annotated (Figure 7; Supplementary Table 4). The carbon fixation pathways mainly included 7 carbon fixation pathways in prokaryotes (average relative abundance ranged from 5.98 to 6.58%) and 8 pathways in photosynthetic organisms (from 3.01 to 3.26%), such as reductive citrate cycle (M00173), dicarboxylate-hydroxybutyrate cycle (M00374), and reductive pentose phosphate cycle (Calvin cycle) (M00165), etc. The nitrogen metabolism pathways included nitrogen fixation (M00175), nitrification (M00528), denitrification (M00529), and 2 nitrate reduction to ammonium pathways (M00530 and M00531), with average relative abundances from 0.43 to 0.81%. The sulfur metabolism pathways included 2 sulfate reduction pathways (M00176 and M00596), with relative abundance ranging from 0.55 to 0.66%. But all the metabolism pathways above showed no clear seasonal variations.

SEM showed that MOD0, MOD3 and MOD4 contributed significantly to nitrogen and sulfur metabolism (Figure 8A). MOD0 was mainly composed of denitrification and nitrate reduction bacteria, so the denitrification was significantly influenced directly by MOD0, with path coefficient value of 0.95, and the nitrate reduction was not significantly influenced directly by MOD0. MOD3 also directly contributed significantly to nitrate reduction and denitrification, with coefficients of 0.84 and 0.44, respectively. Meanwhile, MOD4 played an important role in sulfate reduction.

On the other hand, the microbial community function was indirectly affected by environmental factors through certain ecological modules. Consistent with the result of VPA (Figure 8B), the antibiotics and DIN had significant impacts on microbial communities (with 5.34 and 0.53% explanation, respectively). MOD3 and MOD4 were resistant and promoted by the antibiotics. Trichococcus and norank_f__Leptotrichiaceae in MOD3 and Acidovorax and norank_p_Saccharibacteria in MOD4 carried resistance genes of tetracycline and multidrugs (Zhang et al., 2021), providing a stronger survival advantage in winter with higher concentration of antibiotics. DIN showed negative influence to MOD0, MOD3 and MOD4. The high nitrate content would produce significant inhibition or toxicity to bacteria (Cua and Stein, 2011), as well as to sulfate-reducing bacteria (He et al., 2010; Zhou and Xing, 2021). Besides, the high nitrate suppression to sulfate-reducing bacteria was diminished by the denitrifier (Zhou and Xing, 2021), which also explained the strong correlation between MOD4 and MOD0.

Denitrification, nitrate reduction and sulfur reduction were indirectly promoted by antibiotics via MOD3 and MOD4, whereas they were inhibited by DIN. As the substrate of nitrogen metabolism, DIN also indirectly promoted the process of nitrification, denitrification and nitrate reduction. SEM results showed that the total effect value of environmental factors on nitrification and nitrate reduction was −2.53 and 20.39, respectively. We suggested that the antibiotics and DIN inhibited nitrification, but strongly promoted the nitrate reduction, which lead a further accumulation of NH₄⁺. Although the effects of antibiotics and inorganic nitrogen pollution on microbial ecological modules and their functions in karst rivers was indicated, the influence relationship of different ecological groups on environmental factors and the possible environmental utility was identified, the influence on metabolic process and the fluxes of nitrogen and sulfur still requires further researches such as long-term monitoring, isotope tests and other experimental methods.