The mesocosm experiment was constructed to simulate the aquatic ecosystem. Surface water (1–10 cm) and surface sediments (1–10 cm) were collected from lakes at the Guangdong Ocean University (Zhanjiang, China), removing larger stones and other debris, and divided into transparent plastic boxes (50 × 40 × 30 cm). A plastic box contained 30 L of lake water and 10 cm of sediment. Keep at ambient temperature (25 ± 3°C) for 3 days until the aquatic mesocosm systems have stabilized. Florfenicol solution was added so that the concentration of florfenicol in the water column was 0, 0.1, 1, 10, and 100 mg/L. Water samples (numbered W0 to W4) were collected on days 0, 7, 30, and 60 (D0, D7, D30, and D60). Among them, the 0-day water samples (D0) were collected before the addition of florfenicol. The experiment had 5 groups with 3 replicates each. Three replicate samples from each group were mixed and immediately split into two parts for the determination of physicochemical parameters and DNA extraction.

Several experiments were performed to assess the physiochemical characteristics of the samples. The pH and conductivity of the water samples were measured using a pH meter and conductivity meter, and nitrate nitrogen (colorimetric method) and ammonium nitrogen (indigo method) were estimated using an ultraviolet-visible spectrophotometer. Each sample was analyzed in triplicate.

According to the manufacturer's instructions, genomic DNA was extracted from water samples using the water DNA kit (Omega Bio-Tek, Norcross, GA, USA). The concentration and purity of the extracted DNA were assessed using a nanodrop UV-vis spectrophotometer. Three DNA of samples were extracted from each group for analysis. Based on previous work, the SYBR-primer method was employed with the CFX Connect Real-Time System instrument (BIO-RAD) to determine the abundance of the nirS gene. The denitrifying bacterial nirS gene was amplified using the primer pair CD3AF/R3CDR. The qPCR reaction system was: 10 μL SYBR, 1 μL DNA template, 1 μL upstream and downstream primers (10 mol/L) each, and 8 μL ddH₂O. The qPCR reaction program was: 94°C, 5 min; 94°C for 30 s, 58°C for 30 s, 72°C for 30 s, 39 cycles.

The DNA samples from each of the three replicates were combined into the mixed sample and sent to Jinweizhi (Suzhou, China) for Illumina MiSeq paired-end sequencing of the nirS gene amplicon regions. The nirS gene was amplified for sequencing analysis using the same primer pairs as qPCR. The high-throughput sequencing raw data was optimized using QIIME 1.9.1 by splicing overlapping regions at the ends of the sequences, removing sequences that were shorter than 200 bp, and removing chimeric sequences to produce valid data. An operational taxonomic unit (OTU) cluster analysis was carried out based on 97% similarity. The species taxonomic annotation was completed using the NCBI database. The 20 water samples were divided into different groups according to the sampling time and florfenicol addition concentration, and then OTU clustering analysis was performed. According to the OTU cluster analysis results, a Venn diagram was drawn. Based on OTU analysis results are obtained, using the method of random sampling sample sequences is flat, calculate Shannon, Chao1 alpha diversity index, community species abundance and diversity. Through the (UN) weighted unifrac analysis and comparison between samples whether there are significant differences in the microbial community. PCoA is display beta diversity visualization, PCoA is based on the distance between the matrix Brary - Curtis. Through the weighted clustering hierarchy and the group average method construct UPGMA (Unweighted pair group method with arithmetic mean) clustering tree. And diversity index was examined using QIIME 1.9.1 software and R language. The water properties and nirS gene abundance were used as independent variables, and denitrifying bacteria relative abundance was used as the response variable. By establishing a multiple linear regression model, the response variables were explained by the independent variables, and the correlation between them was explored.

Excel 2016 was used to analyze the nirS gene abundance, R 3.3.1 was used to analyze alpha diversity, and the correlation between environmental factors and community structure was examined using CANOCO 5.0 software. The two-factor variance analysis was applied to test the effects of florfenicol concentration, time, and their interactions on water properties and nirS gene abundance. Statistical analyses were performed with IBM SPSS Statistics 24. All significant differences were at P ≤ 0.05.

As shown in Table 1, florfenicol concentration and sampling time both affected clearly water properties. A higher effect on water properties was observed for sampling time than florfenicol addition concentrations. This trend was not broken by the lower florfenicol additions, but by the higher florfenicol additions. The pH of all water samples was alkaline and weakly alkaline, with a significant decrease in pH on day 7 and a gradual increase in water pH on days 30 and 60. The treatment groups with higher florfenicol additions showed a larger variation amplitude of pH. Higher concentrations of florfenicol caused an accumulation of nitrate nitrogen in the water compared to the control group. In addition, 100 mg/L florfenicol caused a marked accumulation of ammonium nitrogen and a marked rise in water conductivity at day 30. According to the statistical analysis results, concentrations of florfenicol and time significantly affected the water properties in an aquatic microcosm model (P < 0.05).

As shown in Figure 1, the absolute abundance of nirS genes in the water bodies was 2.53 × 10⁴-2.02 × 10⁷ copies/μL. The nirS denitrifying community showed two response patterns to the florfenicol concentration in the water body when stressed by florfenicol. On Day 7, after the addition of florfenicol, the abundance of the nirS gene was 1–2 orders of magnitude higher. Moreover, the abundance of nirS gene increased more with increasing florfenicol addition concentration. On 30 and 60 days, the abundance of nirS gene decreased continuously over time, with a slower decrease in 10 mg/L florfenicol treatment groups. However, with 100 mg/L florfenicol, the abundance of nirS gene kept increasing over time until 60 days, when the abundance of nirS gene reached maximum value and increased by 3 orders of magnitude compared with no florfenicol addition. According to the statistical analysis results, concentrations of florfenicol and sampling time significantly affected the abundance of nirS genes in an aquatic microcosm model (P < 0.05).

After the quality optimization of the original data, the number of valid sequences obtained in water samples was 63,042–186,410. After OTU cluster analysis, the number of OTUs in each sample ranged from 285 to 2,607, with an average of 1,135 OTUs per sample, of which the W4D7 group had the fewest OTUs. The results of the cluster analysis showed that sampling time had an obvious influence on OTU composition in the water column without florfenicol. As shown in Figure 2, the 0-day samples had 512 common OTUs; the four samples without added florfenicol but collected at different times had 82 common OTUs. 1 mg/L florfenicol could markedly reduce the number of common OTUs. Florfenicol concentrations affect its degradation rate and produce different effects on microbial communities in the water body.

The number of valid sequences and OTUs in each sample was obtained from the OTU analysis results. The random sampling method was used to extract a different number of sequences from each sample and calculate the number of OTUs after extraction. The dilution curve between OTU and the number of sample sequences was constructed by using the number of valid sequences extracted from the samples as the abscissa and the number of observed OTUs as the ordinate. As presented in Figure 3, with the increased number of extracted sequences, the number of detected OTUs increases and gradually tends to be flat, indicating that the sequencing is saturated, the sampling is essentially reasonable and the most abundant microbiome is included in the sequencing results. This provided a more realistic representation of the microbial community structure.

The denitrifying microbial community diversity indices mainly include the Shannon and Simpson indices, and the richness indices mainly include the Chao1 and ACE indices. Larger indices represent a greater total species number and higher community diversity. As shown in Table 2, denitrifying microbial community richness and diversity showed a decreasing trend in the water body as the concentration of florfenicol increased; with the change of time, the community richness and diversity were first suppressed and then gradually recovered. Its recovery rate was bounded by the florfenicol concentration. As shown in Figure 4, the Chao1 and Shannon indices of the denitrifying microbial community decreased with increased florfenicol concentrations. The most inhibition of denitrifying microbial community richness and diversity index was observed in the 100 mg/L florfenicol treatment group on day 7.

The denitrifying microbial community composition of each sample was statistical under different species classification levels, as shown in Figure 5. The identified species were classified into 7 phylum, 12 classes, 2 orders, 23 families, 52 genera, and 106 species. The relative abundance of each species in the denitrifying communities was calculated, with the relative abundance of unclassified species ranging from 4.57% to 87.83%. At the phylum level, the relative abundance of Proteobacteria ranged from 8.43% to 95.37% and mainly included Betaproteobacteria, Alphaproteobacteria, and Gammaproteobacteria, which were the absolute dominant bacteria. The relative abundance of Betaproteobacteria was the highest and the relative abundance of Gammaproteobacteria was the lowest, but the relative abundance of Gammaproteobacteria in the W4D60 treatment group was 91.22%. However, the relative abundance of Gammaproteobacteria was 91.22% in 60-day water sample with 100 mg/L florfenicol. At 7 days, the relative abundances of Gemmobacter and Comamonas were increased by lower florfenicol concentrations while higher florfenicol concentrations decreased the abundance of both. However, Comamonas was more sensitive to florfenicol than Gemmobacter, and 1 mg/L florfenicol reduced its abundance. The relative abundance of Pseudomonas and Pelomonas was < 0.2% on day 0, but they showed different responses to the florfenicol concentration on day 60. The relative abundance of Pelomonas was suppressed in groups with florfenicol added, except for 0.1 mg/L, where the relative abundance of Pelomonas reached 71.79%. The relative abundance of Acidovorax was not markedly changed at 7 days and markedly suppressed at 30 days, whereas higher florfenicol concentrations increased the relative abundance of Acidovorax at 60 days and returned to the 0-day level. In addition, the relative abundance of Azovibrio showed a steep increase to 29.84% in the 30-day water sample with 10 mg/L florfenicol. Higher florfenicol concentrations also caused enrichment of some unclassified species, for example, the relative abundance of OTU5, OTU13 and OTU14 could reach 41.83%, 25.77% and 37.55% respectively.

PCoA (Principal Co-ordinates Analysis) and UPGMA (Unweighted pair group method with arithmetic mean) clustering methods were used to express the similarity and differences degree between different samples. As shown in Figure 6A, the contribution of PC1 and PC2 principal coordinates to sample differences were 22.11% and 13.96% respectively in PCoA. The 0-day samples were completely separated from others samples collected at other times, indicating that sampling time had a greater effect on the degree of variation between samples. Samples with higher florfenicol concentrations added were closer together, where they had a higher similarity of denitrifying microbial communities. The similarity of denitrifying microbial communities in samples with lower florfenicol concentrations added was mainly influenced by time. Samples with similar denitrifying microbial community structure were clustered together in the clustering analysis. As shown in Figure 6B, 0-day samples were clustered together; 30- and 60-day samples with lower florfenicol concentrations were clustered with control (W0D30, W1D30, W0D60, W1D60 and W2D60); samples with higher florfenicol concentrations (W3D7, W4D7, W4D30W, W3D30, and W3D60) were clustered together.

The top 10 denitrifying bacteria with higher relative abundance in each sample were selected, and redundancy analysis (RAD) was used to investigate its correlation with water properties and nirS gene abundance. The results are shown in Figure 7, beta proteobacterium and Pseudomonas were strongly correlated with nirS gene abundance and showed a positive correlation, while Gemmobacter, Pelomonas and Paucibacter showed a clearly negative correction. Pseudomonas, beta proteobacterium, Azovibrio and Dechloromonas showed positive correlations with NO 3- and NH 4+; while Gemmobacter, Comamonas and Paucibacter showed negative correlations. Acidovorax was sensitive to variations in pH and showed a significant negative correlation.