1 Introduction
The persistence and spread of antimicrobial resistance (AMR) in the environment, is considered a global silent pandemic and one of the leading health issues facing the world (Mendelson et al., 2022). It is estimated that AMR leads to approximately 700,000 deaths annually and the number is projected to rise to 10 million deaths per year in 2050 (Joshua et al., 2022). Additionally, by 2050, AMR could cost $100 trillion in lost economic output (Burki, 2018). Given the use of antibiotics in aquaculture, fish farms have been suggested as potential sources and sinks for antimicrobial resistance (Watts et al., 2017). Addressing AMR associated with aquaculture practices requires a comprehensive approach, including setting standards for antibiotic use, promoting responsible and prudent use of antibiotics, and adopting alternative disease management strategies and biosecurity to prevent disease outbreaks and reduce reliance on these drugs (Bhat and Altinok, 2023).
One of the strategies commonly used in fish farming to decrease the use of antibiotics is treatment with low dose disinfectants (Rico and Van den Brink, 2014). These disinfectants can help decrease bacterial pathogens on the fish or in the water column (Hushangi and Hosseini Shekarabi, 2018), thereby lowering the risk of systemic infections and minimizing antibiotic dependence. Topical disinfectants can be administered in aquaculture systems continuously, as low dose indefinite baths, or at higher dose as short-term dips (Mota et al., 2022). Regardless of the strategy used, the process rarely eliminates all the bacteria in the system. The remaining bacteria have the potential to regrow and drive the microbial community (Pei et al., 2019; Huang et al., 2023).
It has been shown that bacteria that remerge after disinfection may be more resistant not only to the disinfectant used, but also to antibiotics (Schweizer, 2001; McBain et al., 2002; Sharma et al., 2016). The mechanisms for this selection have still to be elucidated, but specific gene mutations, changes in cell envelope permeability, increased expression of efflux pump, and increased ability to form biofilms have been hypothesized (McDonnell and Russell, 1999; Poole, 2002; Rozman et al., 2021). These findings highlight the importance of considering the potential impacts of disinfection on ARGs particularly in aquaculture where disinfection is used frequently. The objective of this study was to evaluate the effect of low dose ozone on the relative abundance of ARGs in pond water using shotgun metagenomic sequencing analysis.
2 Materials and methods
2.1 Experimental design
The experiment was conducted in August 2022 at the Au Tau Fisheries Office, Agriculture, Fisheries and Conservation Department in Yuen Long, Hong Kong SAR, PR China. The pond water we used had not been treated with any antibiotics and the last H2O2 (7 mg/L) treatment for algae control was done over a year prior to our study. The original aim of the study was to assess the effect of different ozone delivery systems on the microbial communities. The detailed methods and results of the initial study are published elsewhere (Huang et al., 2023). We used the sequence data from this prior study to subsequently assess the effect of ozone on the relative abundance of ARGs. In brief, pond water was pumped into twelve 75L tanks. There were 4 groups in total, two ozone treatment groups, ozone macrobubbles (O3MB) and ozone nanobubbles (O3NB), and two control groups, air macrobubbles (AirMB) and air nanobubbles (AirNB). All groups were triplicated.
To generate ozone nanobubbles, a nanobubbler (aQua+075MO, AquaPro Solutions Pte Ltd., Singapore) was connected to an ozone generator (DNO-15G, Dino Purification Co., Ltd, China) with an incoming oxygen stream flow rate of approximately 0.3 LPM. Water from the tanks was circulated into the nanobubbler and subsequently reintroduced back into the tanks. A dissolved ozone meter (DOZ-30, Dino Purification Co., Ltd, China) was used to monitor the ozone concentration meanwhile, and the system was stopped when the desired ozone concentration of 0.15 mg/L was attained in the three O3NB tanks. the same ozone concentration of 0.15 mg/L was achieved in O3MB group by attaching the same ozone generator to an air stone to deliver gaseous ozone macrobubbles to the tanks. For the AirMB and AirNB group, air pumps (Dazs model AP-528, Hong Kong) with a flow rate of 1.2 L/min were attached to air stones and the nanobubbler, respectively to provide microbubble and nanobubble aeration, and these were run for the same amount of time as the pumps used in the O3MB and O3NB treatment group, respectively.
Water samples (400 mL) were collected before and 24 h after the different treatments, and ozone was not detected at either of these time points (0 mg/L). Water was filtered, DNA was extracted using QIAamp Fast DNA Stool Mini Kit (Qiagen™), and then submitted to Novogene for 150 bp pair-end metagenomic shotgun sequencing. A sequencing depth of 6 Gb (Giga base pairs) was applied for each sample on Illumina HiSeq platform.
2.2 Metagenomic data processing
After the sequencing, raw reads were cleaned using Readfq (V8).1 The clean data were assembled using MEGAHIT software (v1.0.4-beta) (Li et al., 2015). Fragments shorter than 500 bp in the contigs were removed from the assembly. MetaGeneMark (v2.10) (Gemayel et al., 2022) was then used to predict genes on contigs, and the information with a length less than 100 nt in the prediction results was filtered out. We used the CD-HIT software (v4.5.8)2 with the default setting to create a unique initial gene catalog. Clean data of each sample was aligned to this initial gene catalog to determine the number of reads that matched the genes in our catalog using Bowtie2 (Bowtie2.2.4). Genes with reads ≤2 in each sample were filtered out. The filtered data was used for identification of ARGs and subsequent analysis. Based on the number of reads aligned and the length of the genes, the relative abundance of each gene in each sample was calculated by the following formula (Cotillard et al., 2013; Le Chatelier et al., 2013; Zeller et al., 2014), in which r was the number of gene reads, and L is the length of gene.
G
k
=
r
k
L
k
⋅
1
∑
i
=
1
n
r
i
L
i
Notes: “k” represents a certain gene and “i” was a value from 1 to n, used to cycle through all genes.
2.3 Bacterial taxonomy and ARGs analysis
The presence of acquired antimicrobial resistance genes was determined using the ResFinder tool (Florensa et al., 2012). Because we did not know the bacterial origin of most of our genes it was sometimes difficult to classify genes as acquired or intrinsic in this database as in cases such as (blaPAO, blaSGM and oleB) the nature of the gene could differ depending on the bacterial origin. In these cases, we classified the gene as acquired because it has been described in the literature as acquired in at least one bacterial species. To find the intrinsic ARGs belonging to the following 4 families, ATP binding cassette (ABC), resistance nodulation cell division (RND), small multidrug resistance (SMR), and major facilitator super (MFS) family, we used the CARD-RGI tool (McArthur et al., 2013). The threshold values for inclusion in our study were, e-value < 1e−9, identity >70%, and alignment length >100 bp. The Comprehensive Antibiotic Resistance Database (CARD) was also used to determine the mechanism of action and the antibiotic class that the genes conferred to.
DIAMOND (Version 0.9.9.110) was utilized to search genes with the e-value threshold of 1e−5 against sequences of bacteria, fungi, archaea, and viruses obtained from the NR database of the National Center for Biotechnology Information (NCBI). We used the LCA algorithm to ensure accurate species annotation. Finally, we identified whether the ten most abundant bacterial genera, which dominated post treatment had either intrinsic or acquired AMR genes associated with efflux pump proteins.
2.4 Statistical analysis
All data, including specific acquired ARGs, as well as their associated antibiotic classes, and resistance mechanisms, were graphically represented using box plots and stacked bar plots. The data were expressed as mean ± the standard deviation. To assess significant changes in relative abundance between pre- and 24-h post-treatments within treatment groups, we performed paired t-tests after log-transforming the data when the model residuals did not meet the assumptions of normality and homogeneity of variance. A p-value threshold of less than 0.05 was used to determine statistical significance. We also evaluated the difference in the relative abundance of intrinsic resistance genes for efflux protein families, RND, ABC, SMR and MFS, between pre and 24 h-post treatment in the four treatment groups using the same statistical methods. All statistical analysis were done in R Statistical Software (v4.1.2) (R Core Team, 2021).
A dendrogram was used to visualize the hierarchical cluster analysis, which was performed using Ward’s clustering method with squared Euclidean distance to determine if the acquired ARG profiles between treatments were different (Ward, 1963).
Additionally, we investigated whether acquired or intrinsic ARGs co-existed within the same genera of bacteria for overall and post-ozone treated samples by using a Sankey chart, which linked the relationships between these ARGs and the bacterial genera, while highlighting the interplay of resistance mechanisms. This was done using the “networkD3” package in the R Statistical Software (v4.1.2) (R Core Team, 2021).
4 Discussion
Our shotgun metagenomic analysis of pond water treated with low dose ozone (0.15 mg/L) indicated that a one-time treatment, irrespective of how this gas was delivered (nanobubbles or macrobubbles), played a role in influencing the composition and relative abundance of antimicrobial resistant genes (ARGs) 24 h after treatments. Our findings suggested an increase in intrinsic efflux and acquired ARGs. Five out of 33 acquired ARGs detected in our samples increased significantly in relative abundance 24 h after a low dose ozone treatment. These 5 genes were associated with antibiotic efflux, antibiotic activation, and antibiotic target protection mechanisms; however, the efflux pump genes appeared to be more prominent in the ten most abundant bacterial genera post treatment. Other researchers have also reported similar ARG changes in bacterial communities exposed to different disinfectants, such as chlorine, hydrogen peroxide, quaternary ammonium, peracetic acid, chlorhexidine diacetate and trisodium phosphate (Su and D’Souza, 2012; Capita et al., 2013; Biswal et al., 2014; Molina-González et al., 2014; Gadea et al., 2016). Interestingly, we found an increase in relative abundance for a few specific non-efflux ARGs after ozone treatments, but closer inspection revealed that they usually co-existed in bacterial genera associated with intrinsic or acquired efflux ARGs.
The increased relative abundance of both intrinsic efflux and acquired ARGs associated with an efflux pump mechanism after treatment with ozone was somewhat expected. Many of these genes play a crucial role in the survival of bacteria, particularly Gram-negative bacteria, by expelling toxic substances, including antibiotics and other compounds such as those generated during oxidative stress (Nelson, 2002; Huwaitat et al., 2016; Masi et al., 2017). In fact, many researchers have found both acquired efflux pump genes such as qepA, tet(G) and oqxB, and the intrinsic efflux genes in the RND family to be up-regulated in response to different stressors, suggesting they play an important role in bacterial survival (Shi et al., 2013; Tandukar et al., 2013; Xu et al., 2017; Gu et al., 2020; Hu et al., 2023; McCarlie et al., 2023). Therefore, it is not surprising that we detected an increase in the relative abundance of genes associated with the efflux pump mechanism after exposure to low dose ozone, which is associated with oxidative stress (Pell et al., 1997).
It is important to better understand the selective pressure that ozone exerts on bacterial communities at different doses used in aquaculture particularly low dose treatments that may not kill all bacteria (Iakovides et al., 2019, 2021; Huang et al., 2023). Our study suggests at lower ozone concentrations, which are used to treat fish, bacteria are not completely eliminated (Huang et al., 2023), and the surviving bacteria may be more resistant to ozone. This shift in population could have additional consequences on other types of treatments, particularly antibiotic therapies as the genes suspected of conferring ozone resistance may also bestow resistance to antibiotics.
Besides an increase in ARGs associated with efflux pumps, we also observed a significant increase in the relative abundance of a few non-efflux genes either involved in antibiotic target protection or antibiotic inactivation after our ozone treatments. Previous studies also found an increase in antibiotic specific genes in wastewater after treatment of with 0.9 ± 0.1 g ozone per 1 g DOC. In one study, there was an increase in the abundance of vanA and blaVIM genes, encoding vancomycin and carbapenem resistance (Alexander et al., 2016). Another study found after an ozone treatment, the surviving R. ornithinolytica bacterial population had an increase in the blaCTX–M–62 gene (Azuma et al., 2022). It is possible genes that conferred ozone resistance in our study (i.e., efflux pump genes) co-occurred with AMR specific (i.e., non-efflux pump) genes, which could therefore lead to their co-selection under ozone exposure. Co-selection of genes under different environment stressors has been reported by others (Ding et al., 2019). In our study, there was evidence for this hypothesis as we found several non-efflux pump ARG often co-occurring within the same genera of bacteria as acquired efflux ARGs (i.e., OqxB) or the intrinsic RND efflux genes (i.e., adeF). In fact, all but two of the ten most abundant bacteria that survived and colonized our samples after an ozone treatment had either intrinsic RND efflux genes or acquired efflux ARGs. Interestingly, these two genera Sphingomonas and Citrobacter, are known to form biofilms. This protective mechanism protects bacteria against disinfectants (Donlan and Costerton, 2002; Gulati and Ghosh, 2017) and may explain their survival post ozone treatment.
Notably, a few ARGs such as dfrA1, floR and ant(3″)-Ia, which increased in relative abundance after 24 h after ozone treatments are known to be located on mobile genetic elements (MGEs) (Zhang et al., 2017; Sánchez-Osuna et al., 2020; Raza et al., 2022; Zhu et al., 2022). If the increased ARGs observed in our study post treatment were present within MGEs, the enrichment of those ARGs could lead to nonpathogenic microbial species acting as an ecological reservoir for pathogenic bacteria (Peterson and Kaur, 2018). Bacteria can acquire new ARGs through horizontal gene transfer from plasmids, transposons, or other MGEs originating from other bacteria or the environment (Jian et al., 2021).
This study had several limitations that should be acknowledged. Firstly, our results were based on only one pond. Different pond water sources could be included in future studies to increase the external validity of the findings. Secondly, we used shotgun metagenomic analysis which only provided us with short-read sequencing. This is particularly problematic for determining the origin of genes within the mixed bacterial community. Despite attempts to address this limitation by employing read assembly techniques (Li et al., 2022), we could not determine the source of many ARGs (i.e., bacterial genomic or MGEs). Another limitation is that the classifiers we used to determine whether genes were intrinsic or acquired were not consistent across all bacterial genera. There were some genes (i.e., blaPAO, blaSGM, oleB) that are likely intrinsic for some bacteria and considered acquired for others, therefore, it was difficult for us to categorize these as we did not know the bacterial origin of all the genes. Despite this limitation, of the genes we classified as acquired (n = 33), there were only a few which could have been intrinsic depending on the origin (Li et al., 2022). Finally, it is important to note that our reporting of relative changes in certain genes between samples, was achieved by normalizing the reads to the total sequence reads, therefore we could only assess relative abundance. This measure does not necessarily translate to phenotypic antibiotic resistance or absolute abundance. To investigate the absolute abundance of ARGs, quantitative polymerase chain reaction (qPCR) analyses would be required, and this was beyond the scope of this study.