Emerging antimicrobial resistance is considered by the World Health Organization as a global public health crisis that must be managed with the utmost urgency^[1]. Part of the problem is the contamination of aquatic ecosystems with the antibiotics that we use for the treatment or prevention of bacterial infection, as well as in agriculture and aquaculture^[2-3], which is attributed to the treatment of wastewater often being poor or even non-existent^[4]. The input of antibiotics may challenge the bacterial communities^[5] and accelerate the appearance and spread of antibiotic resistance genes (ARGs)^[6]. ARGs can be transferred horizontally between different microbial communities, increasing the risk of ARG dissemination^[7]. Therefore, ARGs are considered emerging contaminants and raise increasing public concern^[8], especially with respect to food safety as they have been detected in bacteria inside edible aquatic species^[9].

ARG contamination has been investigated in various aquatic environments, such as lakes^[10-11], sewage^[12-13], aquacultures^{[3, 14]}, oceans^[15], and even drinking water^[16-17]. Several past publications indicated that urban water systems around the world are highly contaminated by ARGs^[18-20]. However, few studies have examined the occurrence and abundance of ARGs in urban recreational waters, despite the latter being perceived as aesthetically and ecologically beneficial elements of urban environments^[21]. Moreover, pathogenic bacteria and ARGs may find their way to recreational waters because of the input of sewage discharge and animal excrement^[22]. The intensity of anthropogenic activities has been shown to affect the abundance of ARGs^[6]. Recreational water contaminated with ARGs may pose a risk for public health, as citizens often come in contact with it through swimming and fishing. Thus, it is necessary to understand the occurrence of ARGs in urban recreational lakes.

To date, studies of ARGs in aquatic environments, including recreational waters, have been using traditional approaches of molecular biology, such as PCR which only allowed them to focus on a specific gene^{[3, 23-24]}. The bias deriving from the very selection process of the genes to amplify, as well as limitations related to how well the designed primers work, shed doubt on whether these studies reflect the real composition of ARG and microbes (especially uncultured microbes) in recreational waters^[25]. In recent years, the advent of next-generation sequencing and the development of bioinformatics gave rise to metagenomics, a powerful tool for exploring the complex microbial diversities and creating broad-spectrum profiles of ARGs in samples from various aquatic environments^[26-28].

This study aimed to utilize a high-throughput sequencing-based metagenomic approach to characterize the composition and abundance of the microbes and ARGs in Jiushan Lake, a recreational lake located in the city of Wenzhou, and estimate the potential public health risk. Moreover, because the function of urban recreational lake waters as potential reservoirs of ARGs and bacteria is subject to the temporal variation of environmental conditions, we chose to take the temporal factor into account. This work should provide useful information for a better understanding of the ecosystem structure and function of urban recreational lakes.

1 Materials and Methods 1.1 Main reagents and equipment

PowerWater DNA Isolation Kit, Qiagen company; Portable pH meter, Mettler Toledo company; Portable dissolved oxygen meter, YSI company; Polycarbonate filters, Millipore company; Microspectrophotometer, DeNovix company; GF/F glass fiber filters, GE Healthcare company; Illumina HiSeq platform, Illumina company.

1.2 Site information and sampling method

Wenzhou, a Chinese city located in southeastern Zhejiang province, is close to the East China Sea and has abundant water resources. Jiushan Lake is a shallow body of water situated in the central city district that serves as a very popular natural swimming pool for locals. Water samples were collected at two different time points, in August and December 2018, representing summer and winter, respectively. Three different sampling sites (locations with frequent citizen activities) were selected, and three surface water samples (to be exact, from a depth of about 0.5 m) were collected for each season. The samples were stored in sterilized Nalgene 1-liter polycarbonate bottles and transported to the laboratory in an icebox within 1 hour from collection. Each sample was divided into two parts, one for physicochemical analysis and the other for analysis of microbial genomes.

1.3 Physicochemical parameters

Water temperature and pH were measured on-site using an pH meter. Dissolved oxygen (DO) concentration was also measured on-site with a YSI 550 portable dissolved oxygen meter. After filtration through GF/F glass fiber filters, Chlorophyll-a (Chl-a) was extracted from the filters with pure methanol for 24 h, centrifuged at 3 000 r/min for 10 min, and then spectrophotometrically quantified^[29]. Soluble reactive phosphorus (SRP) was measured using the phosphorus molybdenum blue method, ammonium nitrogen (N-NH₄) with the sodium hypobromite oxidation method, nitrate-nitrogen (N-NO₃) with the Zinc cadmium reduction-azo colorimetric method, and nitrite-nitrogen (N-NO₂) using the diazo colorimetric method^[30].

1.4 DNA extraction and shotgun metagenomic sequencing

Half of each water sample, i.e., 500 mL, was pre-filtered through 5 μm pore-size polycarbonate filters to remove impurities. Then, each filtrate was filtered again using a 0.22 μm membrane filter to capture and concentrate microbes. The filters were cut into strips using sterile scissors, put into 5 mL sterile tubes, and stored at -20 ℃ until DNA extraction. DNA was extracted using the PowerWater DNA Isolation Kit according to the manufacturer's instructions. The DNA quality was assessed using a microspectrophotometer. All samples had OD_260/280 ratios between 1.8 and 2.0. DNA degradation degree and potential contamination were electrophoretically examined on 1% agarose gels. Qualified DNA samples were sent to Novogene (Beijing, China) for Illumina shotgun high-throughput sequencing using the PE 150 (Paired-end sequencing, 150-bp reads) sequencing strategy. A total amount of 1 μg DNA per sample was used as input material for library construction with an insert size of 350 bp, followed by high-throughput sequencing on an Illumina HiSeq platform.

1.5 Bioinformatics

The raw sequencing data were pre-processed for quality control using Readfq V8 (public domain) as follows: a) Removal of reads that contained more than 26.7% (40/150) low-quality bases (the default quality threshold value was set at 38); b) Removal of reads with more than 6.7% (10/150) undetermined bases (N), and; c) Removal of reads that shared an overlap with the adapter sequence longer than 15 bp. Q20, Q30, and GC-content were also calculated. The sequencing output of each water sample was about 12 Gb with more than 85% clean-read bases, i.e., having quality scores higher than 30. The base sequencing error rate was 0.001. A total of 505 097 533 clean reads were generated from 6 samples with an average of 84 182 922 reads per sample. The filtered high-quality sequences were used for subsequent analysis.

To investigate the microbial composition of the water, genes were aligned to the NR database (Version: 2018-01-02) of NCBI using DIAMOND V0.9.9 software (free domain)^[31] with an E-value≤10^-5. The abundance of a taxonomic group was calculated by summing the abundance of genes annotated to a feature. The identification of human pathogenic bacteria referred to the list of human pathogenic bacteria (51 species) provided in a previous study^[32]. All values are averages from three different sequencing runs.

To explore the diversity and abundance of ARGs, genes were aligned to the CARD database using Resistance Gene Identifier (RGI) software^[33-35] for annotation and BLASTp comparison (E-value≤10^-30). The alignment results were used to calculate both the total and the relative abundance of ARGs. The proportions of different ARG sequences in "total detected ARG sequences" and "total metagenome sequences" were defined as " relative percentage" (%) and "relative abundance" (ppm (parts per million), one hit in one million aligned sequencing genes), respectively^[36-37]. All values are averages from three different sequencing runs.

1.6 Statistical analysis

Diversity index (Shannon, Inverse Simpson) was calculated based on the genera profiles. Histograms were performed at OriginPro 8.5. Redundancy analysis (RDA) was performed to investigate the correlation between the top 10 most abundant genera and environmental factors. Canonical correspondence analysis (CCA) was used to investigate the relationships between environmental factors, human pathogenic bacteria, and the distribution of ARGs. RDA and CCA were conducted in CANOCO V5 (Biometris, Wageningen, the Netherlands) for windows. One-way ANOVA (LSD test), performed with IBM SPSS Statistics V22 (IBM, Armonk, NY, USA), was used to assess the statistical differences between data sets. The significance level was set at P=0.05.

Data are available at the NCBI Sequence Read Archive under project No. PRJNA613875.

2 Results 2.1 Environmental parameters

The physical and chemical properties of the water samples are shown in Table 1. The pH and dissolved oxygen (DO) values did not display significant seasonal differences between the two seasons. The temperature was 34 ℃ in the summer and 18 ℃ in the winter. Chlorophyll-a concentration was 37.6 μg/L and 17.1 μg/L in the summer and winter, respectively, i.e., the summer concentration was 2.2 times higher than the winter concentration. With respect to nutrients, SRP was 1.5 times higher in the summer than that in the winter water samples. In contrast, the concentration of dissolved inorganic nitrogen, i.e., the sum of N-NO₂^-, N-NH₄⁺, and N-NO₃^-, in the summer samples was only 30% of the one in winter samples.

2.2 Microbial abundance and composition

At the kingdom level, in summer water samples, 54% of genes were aligned to Bacteria, 9% of genes to Viruses, 0.4% to Archaea, and 0.04% to Eukaryota, while about 36% of the genes could not be annotated using the NR database. In winter samples, 64% of genes were aligned to Bacteria, 5.9% to Viruses, 0.3% to Archaea, and 0.04% to Eukaryota, while 29.41% of the genes were not annotated. Notably, the relative abundance of Bacteria in the winter was higher than that in the summer.

This study identified 148 phyla in the summer and 152 phyla in the winter samples (including Bacteria, Eukaryota, and Archaea), respectively. As seen in Figure 1, the phyla with the highest relative abundance greater than 1% in both seasons were Proteobacteria, Actinobacteria, Bacteroidetes, Cyanobacteria, and Verrucomicrobia. Several important observations can be made. Proteobacteria and Actinobacteria dominated samples from both seasons, together accounting for nearly 35% and 46% of the microbial community in summer and winter samples, respectively. Next in order of abundance in the summer water samples were Cyanobacteria and Bacteroidetes, while the corresponding positions in winter samples were occupied by Bacteroidetes and Verrucomicrobia. Proteobacteria, Actinobacteria, and Bacteroidetes had significantly higher abundance in winter samples than that in summer samples (P < 0.05), while both Cyanobacteria and Verrucomicrobia were more abundant in summer than those in winter (P < 0.01).

At the genus level, bacteria with high relative abundance in both seasons were those belonging to the genera Limnohabitans (a genus of Comamonadaceae family, Burkholderiales order, Betaproteobacteria class, Proteobacteria phylum), Synechococcus (a genus of Chroococcales family, Cyanobacteria class, Cyanobacteria phylum), Candidatus Planktophila (a genus of Candidatus Nanopelagicales order, Actinobacteria class, Actinobacteria phylum), Polynucleobacter (a genus of Burkholderiaceae family, Burkholderiales order, Betaproteobacteria class, Proteobacteria phylum), and Candidatus Methylopumilus (a genus of Methylophilaceae family, Nitrosomonadales order, Betaproteobacteria class, Proteobacteria phylum) (Figure 2). Synechococcus was the most common genus in summer samples (2.78%) and its abundance was significantly higher in summer samples than that in winter (P < 0.01), while Limnohabitans was significantly higher in winter than that in summer (P < 0.05). RDA was performed to investigate the correlation between the top 10 most abundant genera and environmental factors. The results (Figure 3) demonstrated that NH₄⁺-N, NO₂^--N, and temperature were found to be significantly correlated with the distribution of the main planktonic bacteria (P < 0.05).

In the diversity analysis, the trend of Shannon and Simpson index was consistent, and the results showed that a high α diversity was observed in summer than that in winter (Table 2). Of particular interest was the presence among the identified species in this study associated with human opportunistic infections, such as Pseudomonas aeruginosa, Escherichia coli, Chlamydia trachomatis, Clostridium botulinum, and Bacillus cereus. Amongst, P. aeruginosa (a species of Pseudomonas genus, Pseudomonadaceae family, Pseudomonadales order, Gammaproteobacteria class, Proteobacteria phylum) (0.013 3%/0.020 8%) had the highest abundance in all samples.

2.3 Diversity and abundance of ARGs

A total of 449 ARGs were detected in this lake, 370 in summer and 383 in winter water samples, respectively. More than half of the ARGs (304 ARGs) were shared, 66 were summer-exclusive, and 79 winter-exclusive (Figure 4). The relative abundance of ARGs in the summer samples (359.02±7.0 ppm) was higher than that in the winter samples (318.94±11.6 ppm) (Figure 5A), with the difference being statistically significant (P < 0.05). The ARGs with the relative percentage were shown in Figure 5B. The predominant resistance gene in the summer samples was MCR-1.2 (polymyxin E resistance gene) (4.66%), followed by LsaC (aminoglycoside resistance gene) (4.46%), APH6-Ia (aminoglycoside resistance gene) (3.16%), and LpeB (macrolide resistance gene) (2.40%). In winter samples, BcI (β-lactam resistance gene) (7.35%) was the predominant ARG, followed by LRA-13 (cephalosporin resistance gene) (3.73%), acrD (aminoglycoside resistance gene) (3.53%), and SHV-180 (β-lactam resistance gene) (2.60%). The canonical correspondence analysis (CCA) (Figure 6) indicated that the environmental variables (NO₃^--N, NO₂^--N, and SRP) and human pathogenic bacteria (P. aeruginosa, E. coli, and C. botulinum) showed significant relationships with the distribution of ARGs (Monte Carlo test, P < 0.05).

Figure 7 shows the percentage of resistance mechanisms represented in the ARGs repertoire of summer and winter water samples. In both seasons, three resistance mechanisms were dominant, accounting for more than 90%: antibiotic efflux, antibiotic inactivation, and antibiotic target alteration. However, their relative representation depended on the season. Specifically, antibiotic inactivation had a much higher representation in winter than in summer, whereas the percentages of the other two major mechanisms showed the opposite trend. In both seasons, antibiotic efflux and inactivation were the largest groups, with a collective percentage of 73.94% in summer and 80.37% in winter.

3 Discussion

In this study, we used Illumina high-throughput shotgun metagenomic sequencing to investigate the seasonal composition and abundance of microbial communities and ARGs in Jiushan Lake. Proteobacteria and Actinobacteria were the predominant phyla in all water samples, which was consistent with previous reports on other freshwater environments, such as general freshwater lakes^[38-39] and freshwater aquacultures^[40-41]. At the genus level, the relative abundance of Synechococcus spp. was much higher in summer than that in winter, while the opposite was true for Limnohabitans spp. and Candidatus Planktophila spp. This seasonal difference in bacterial community composition may be caused by the seasonal variation in physical and chemical properties such as nitrogen availability, temperature, pH, and total phosphorus^[42-44]. For example, inorganic nutrients are often the main limiting factor for microbial growth^[45]. Working with samples from five lakes with different trophic status and humic content, Lindström showed that the quantity and community structure of bacteria vary in different habitats and nutrient conditions^[46]. Synechococcus spp. are photoautotroph organisms belonging to Cyanobacteria^[47], the result of RDA demonstrated that Synechococcus was positively correlated with temperature, thus the higher summer temperature may provide favorable conditions for its growth^[48]. On the other hand, Limnohabitans spp. belongs to βetaproteobacteria, which are nitrogen cycle-related organisms with a wide distribution in freshwater systems^[42]. Total nitrogen plays an important role in the growth and reproduction of this genus^[49]. Consistent with these reports, our study shows a positive correlation between NH₄⁺-N, NO₂^--N and Limnohabitans (Figure 3).

Extremely worrisome is our detection of species typically associated with human opportunistic infections, such as P. aeruginosa, E. coli, C. trachomatis, C. botulinum, and B. cereus, as these findings suggest a clear health hazard to citizens coming into contact with the water of the lake during recreational activities. P. aeruginosa was the most abundant potential pathogen in all water samples. It has been reported that the infection of human external ear and skin in recreational waters is related to the presence of P. aeruginosa^[50] and it has attracted more attention because of the multidrug resistance of the widely specific multi-drug efflux system^[51]. B. cereus primarily causes food poisoning and severe eye infections^[52], Resistance of B. cereus to beta-lactam and carbapenem has also been reported^[53]. Drug-resistant pathogens can cause intestinal problems in swimmers by swallowing in lake water, and opportunistic pathogens found in recreational waters can pose a considerable potential threat to public health.

With respect to ARGs, our study is consistent with previous reports according to which urban aquatic environments, such as rivers and lakes (including recreational waters), may serve as ARG reservoirs^{[23, 26]}, as we were able to detect hundreds of ARGs in Jiushan Lake. Most of the detected ARGs endowed resistance to aminoglycoside, β-lactam antibiotics, or macrolide. Moreover, we observed a significant difference in the abundance and composition of ARGs between summer and winter water samples. Human pathogenic bacteria (P. aeruginosa, E. coli, and C. botulinum) showed significant relationships with the distribution of ARGs, implying the threatening risk of antibiotic resistance to the public health. The relative abundance of ARGs was higher in summer than that in winter. With respect to the cause, a seasonal study of ARGs in a peri-urban river denoted temperature as a potential factor driving the shift of ARG profiles^[54]. However, other factors may also be at play, such as the increased human presence in the lake during summer.

MCR-1.2 was the predominant ARG in summer water samples. It is a novel variant of MCR-1 (mobilized colistin resistance gene -1), a plasmid-mediated resistance gene for polymyxin E (also called colistin, hence the gene's name)^[55-56]. Polymyxin E is a polypeptide antibiotic produced by Bacillus polymyxavar and is figuratively called the "last line of defense" against gram-negative bacterial infection^[57-58]. To date, MCR-1.2 has been found in Klebsiella pneumonia and E. coli^{[55, 59]}, both of which were detected in this study. As ARGs can be transferred between different bacteria via plasmids during conjugation, MCR-1.2 will ultimately find its way to opportunistic human pathogens, at which point it will pose a hard challenge to clinical medication and a significant threat to human health. The detection of MCR-1.2 in Jiushan Lake, whose waters come in contact with a large number of swimming enthusiasts, should arouse the public's concern.

BcI was the preponderant ARG in winter water samples. It codes for Bacillus cereus beta-lactamase, which is a zinc metallo-beta-lactamase that hydrolyzes a large number of penicillins and cephalosporins^[60]. β-lactam antibiotics are the most widely used antibiotics^[61], which might explain why ARGs against them are widely distributed in aquatic environments. An example is the detection of the other β-lactamase-coding gene, bla_-TEM, in recreational lakes in parks of Beijing^[23].

The predominant resistance mechanisms detected in Jiushan Lake were antibiotic efflux and antibiotic inactivation, which were also found to be the chief mechanism in paddy soils from South China^[27] and the receiving surface water of a wastewater treatment plant^[62]. Antibiotic inactivation is mainly associated with resistance to common anthropogenic antibiotics, such as aminoglycosides, β-lactams, and macrolides^[36]. Efflux pumps can transfer not only antibiotics but also heavy metals and other toxins from the cell to the external environment and are involved in various processes such as detoxification of metabolic intermediates, virulence, and signal trafficking^[63]. Thus, antibiotic extrusion by efflux pumps is the most effective method to resist the different pressures of diverse environments^[64-65].

4 Conclusion

This study provides a description of the microbial community and ARG composition in an urban recreational lake as assessed by metagenomic analysis. The microbial diversity and relative abundance of ARG were higher in summer than in winter. MCR-1.2 and BcI were the predominant ARGs in the summer and winter water samples, respectively. An expanded study with samples from several months should help to understand the composition and dynamic change of bacteria and ARGs more clearly. More importantly, as the high abundance of observed ARGs may pose a public health risk, our results suggest that more attention should be paid to the ARG contamination of urban recreational waters.