Prevalence and proliferation of antibiotic resistance genes in the subtropical mangrove wetland ecosystem of South China Sea

Abstract The emerging pollutants antibiotic resistance genes (ARGs) are prevalent in aquatic environments such as estuary. Coastal mangrove ecosystems always serve as natural wetlands for receiving sewage which always carry ARGs. Currently, the research considering ARG distribution in mangrove ecosystems gains more interest. In this work, we investigated the diversity of ARGs in an urban estuary containing mangrove and nonmangrove areas of the South China Sea. A total of 163 ARGs that classified into 22 resistance types and six resistance mechanisms were found. ARG abundance of the samples in the estuary is between 0.144 and 0.203. This is within the general range of Chinese estuaries. The difference analysis showed that abundances of total ARGs, six most abundant ARGs (mtrA, rpoB, rpoC, rpsL, ef‐Tu, and parY), the most abundant resistance types (elfamycin, multidrug, and peptide), and the most abundant resistance mechanism (target alteration) were significantly lower in mangrove sediment than that in nonmangrove sediment (p < 0.05). Network and partial redundancy analysis showed that sediment properties and mobile genetic elements were the most influential factors impacting ARG distribution rather than microbial community. The two factors collectively explain 51.22% of the differences of ARG distribution. Our study indicated that mangrove sediments have the capacity to remove ARGs. This work provides a research paradigm for analysis of ARG prevalence and proliferation in the subtropical marine coastal mangrove ecosystem.


| INTRODUC TI ON
Antibiotics have effectively helped humans fight bacterial infections for many years, but the increase in antibiotic resistance occurrence in clinically relevant microorganisms has rendered ineffective antibiotics. Around the world, antibiotic resistance is rising to a dangerously high level, leading to increased medical costs, bacterial infections, and mortality. The problem of antibiotic abuse and antibiotic resistance in China is quite serious. Antibiotic use per capita in China is more than five times than that of Europe and the United States (Zhang, Ying, Pan, Liu, & Zhao, 2015). In China alone, roughly 24,748 tons of antibiotics from human and animal sources were released in natural waters in 2013 . This exacerbates the problem of antibiotic resistance in aquatic ecosystems.
However, few research focus on the fate of ARGs in mangrove, where is a common, microbe-rich natural wetland.
Mangrove ecosystems are widely distributed in the intertidal zones of tropical and subtropical coasts, and the largest amount (42%) of the world's mangroves is in Asia (Giri et al., 2011). These ecosystems are especially essential for the global carbon cycle, and have important economic value to human society (Bouillon et al., 2008;Brander et al., 2012;Giri et al., 2011). They are composed of highly coordinated mangrove plants, animals, and microorganisms (Andreote et al., 2012;Bai et al., 2013;Holguin, Vazquez, & Bashan, 2001). Microorganisms in mangrove sediments are enriched during the process of mangrove succession (Chen, Zhao, Li, Jian, & Ren, 2016), and play important functions, such as carbon fixation and degradation, nitrogen fixation, sulfur metabolism, and sediment stabilization (Andreote et al., 2012;Bai et al., 2013;Holguin et al., 2001). Mangrove also receives sewage or aquaculture wastewater in coastal regions. The pollutants in water are removed mainly by tidal scour, plant assimilation, or microbial metabolism (Ouyang & Guo, 2016).
Thus, mangroves are often considered as an economical alternative to sewage treatment plants because of their potential in removing and tolerating pollution. However, previous studies have not taken into account their potential in removing ARGs (Leung, Cai, & Tam, 2016;Ouyang & Guo, 2016;Wong, Tam, & Lan, 1997;Wu, Tam, & Wong, 2008). To our knowledge, some literatures reported that novel compounds as potential new antibiotics (Azman, Othman, Velu, Chan, & Lee, 2015) produce from Actinobacteria and antibiotic resistance microbes (Cabral et al., 2016;Ghaderpour et al., 2015;Jalal et al., 2010;Le, Munekage, & Kato, 2005) in the mangrove sediments, but few research study on distribution pattern of ARGs in the mangrove ecosystem. Therefore, in this work, we analyzed the abundance and diversity of ARGs in the mangrove ecosystem located in an urban estuary of a city in the South China Sea. Via shotgun metagenomics approach, one of the appropriate methods to study the ARG diversity of environmental samples such as sediment, soil, water, etc. (Chen et al., 2013;Chen, Yuan, et al., 2016;Zinicola et al., 2015), we investigated the following contents: (a) antibiotic resistance level of the subtropical city estuary, (b) profiles of ARGs in mangrove and nonmangrove areas that both receive the same municipal wastewater, (c) correlations between ARGs and abiotic and biological factors in the sediments, (d) distribution pattern of antibiotic resistance, and (e) impacts of mangrove on the distribution of ARGs. Our work provides a research paradigm and theoretical support for analysis of ARG prevalence and proliferation in the subtropical marine coastal mangrove ecosystem.

| Sampling sites
The sampling location was the estuary of Fengjiajiang River in Beihai City, Guangxi Province, China (21°24′43.43″N, 109°9′50.98″E), which is surrounded by mangrove and a beach on two sides ( Figure 1) bank is a beach with no plants, and the sediment type is sandy soil.
The northern bank is covered by mangroves, most of which are identified as Avicennia marina; its sediment type is sandy or sandy loam.
Outside the mangrove is a sandy beach. Samples were collected in the mangrove area from inner sites A and D, center site B, and outer site C. Samples were also collected from site E outside the mangrove area. Samples from the nonmangrove area were X, Y, and Z, which were located from the inner to the outer site. Figure 1 shows the geographic map of the sampling sites.

| Sampling collection
All the sediment samples were collected on August 7, 2017. Sterile polyvinyl chloride tubes with a diameter of 3 cm and a height of 15 cm were used to collect samples from the sampling sites, each with a quadrat of 20 × 20 m. Five sediment samples were collected from each site and placed in a sterile bag inside a box filled with ice immediately after sufficient mixing. All samples were stored at −80°C in the same day. The data output from each DNA sample was more than 5 Gb.

| DNA extraction and highthroughput sequencing
The metagenomic sequencing data were deposited in the NCBI SRA database under the BioProject PRJNA472209.

| Sediment property analysis
Sediment properties, such as pH, total organic carbon (TOC), and metal elements (Cu, Zn, Hg, Cd, and metalloid As), were determined using

| Bioinformatic analysis
Based on previous studies, adapter sequences were removed using SeqPrep v1.33. Moreover, sequences <100 bp and those with quality F I G U R E 1 Geographic map of sampling sites. The red labels are the sampling sites that located at the north bank of the estuary, the A, B, C, and D are mangrove sampling sites, and the E is a sandy beach site outside the mangrove. The yellow labels are the south bank sites, all of them are sandy beach sites Z <20 as well as reads containing an N base were removed using Sickle v1.2; then, clean reads were created (Yang, Sun, Yang, & Li, 2014).
Function blastx (set as --evalue 1e-5 -query-cover 75 --id 90 -k 1 --compress 1) in DIAMOND v0.9.14 (Buchfink, Xie, & Huson, 2015) was executed to align the clean reads against the Comprehensive Antibiotic Resistance Database (CARD) (Jia et al., 2017) to annotate the ARG reads. ARGs were classified according to Antibiotic Resistance Ontology terms "antibiotic molecule" and "mechanism of antibiotic resistance" of CARD, which are called "resistance type" and "resistance mechanism" in this study, respectively. The multiple resistance genes were classified into three resistance types: (a) sulfonamide and sulfone; (b) macrolide-lincosamide-streptogramin; and (c) multidrug, other multiple resistance genes. The mobile genetic elements (MGEs) integrase and insertion sequence common region (ISCR) were annotated using the same method as that of ARGs, and the annotation of plasmids was according to a previous study (Chen et al., 2013).
The abundance of ARGs or MGEs was indicated as the ratio of normalized ARGs (or MGEs) to normalized 16S rRNA according to the method performed by Li et al. (2015). Briefly, the abundance of ARGs or MGEs was calculated as the following equation: where N Gene is the reads number of the ARG or MGE; L GeneRef is the reference sequence length of the corresponding ARG or MGE; N 16S rRNA is the reads number of the 16S rRNA; L 16S Ref is the average length of the 16S sequence in the Silva database; n is the number of the ARG or MGE genotype; L reads is the sequence length of the Illumina reads (150 nt) in this study.
Before taxonomic classification, unigenes were annotated according to previous work  with two modifications: we used MEGANHIT v1.1 and SOAPaligner v2.21 with parameters "--presets meta-large" and "-u, −2, -m 200," respectively; we removed the genes that mapped reads less than or equal to 2 in the step "Gene prediction from scaftigs." For taxonomic classification, the total unigenes were aligned against NCBI NR database (V20151204) using DIAMOND v0.9.14 (set as blastp, -e 1e-5). The sequences with e-value not more than the minimum e-value × 10 were selected for the taxonomic annotation by LCA algorithm from MEGAN5 (Huson, Mitra, Ruscheweyh, Weber, & Schuster, 2011). Then, by combining the unigene abundance information, the abundance of all taxonomic levels in each sample can be calculated (Huson et al., 2011;Karlsson et al., 2013).

| Statistical analysis
Bar charts were generated using Excel 2016 (Microsoft, USA). Venn diagrams were created online (http://bioin forma tics.psb.ugent.be/ webto ols/Venn/). Difference between two groups was analyzed by SPSS v23.0 (IBM), statistical test is independent-samples T test, Confidence Interval Percentage is 95%. were chosen for the RDA and pRDA similar to a previous study . These axes account for 87.28%-99.90% of the variance in these factors. Network analysis was performed using Gephi v0.9.2 (Bastian, Heymann, & Jacomy, 2009). At the class level, z-scores represent the standardized relative abundance of taxa were calculated by the following equation: where SG is the relative abundance of individual taxa in a sample; TGA is the average relative abundance of the individual taxa in total samples; SD is the standard deviation of the individual taxa in total samples.

| Abundance and diversity of ARGs, resistance type, and resistance mechanism
The abundance of the ARGs was described as the number of ARGs which normalized by 16S rRNA according to previous studies . The abundances of ARGs in the samples were between 0.144 (sample B) and 0.203 (sample X), which are mangrove site and nonmangrove site, respectively (Figure 2a).
A total of 163 different ARGs were annotated (Table S1 and  (Table S1) according to drug class. The abundances of the major resistance types were 0.703, 0.312, 0.186, 0.097, and 0.030 according to elfamycin, multidrug, fluoroquinolone, aminocoumarin, and peptide, respectively, and the abundances of these five resistance types were 96.39% of the total abundance ( Figure 2c and Table S2).
The following six of the seven resistance mechanisms in CARD were detected: antibiotic efflux, antibiotic inactivation, antibiotic target alteration, antibiotic target protection, antibiotic target replacement, and reduced permeability to antibiotics (Figure 2d). In addition, the term "unknown mechanism" was utilized for the ARG tet34 (oxytetracycline resistance gene) (Table S1). Target alteration, efflux, inactivation, and target protection existed in all samples (Table S3) as major resistance mechanisms with 99.47% of the total abundance ( Figure 2d).

| Antibiotic resistance differences between mangrove and nonmangrove sediments
Excluding sample E, we divided all the samples into two groups: M1 and S1. Group M1 contained samples A, B, C, and D from the mangrove sediments of the north bank of the estuary. As a control group, S1 contained samples X, Y, and Z from the sandy beach (nonmangrove) on the south bank of the estuary.
Abundances of elfamycin, multidrug, and peptide-resistant ARGs showed significant differences (p < 0.05) between the two groups, and they are higher in S1 (Figure 4). The Venn diagram analysis ( Figure   A1b) shows that group M1 and S1 shared most resistance types.
Group M1 contained all the five resistance mechanisms that found in S1, and had unique unknown mechanism (Table S3 and Figure   A1c).

| Sediment properties
To study the effects of sediment properties on antibiotic resistance, the pH, TOC, heavy metal (Cu, Zn, Hg, Cd) and metalloid As contents of the samples were determined (Table S4). The average pH and TOC of M1 group were both higher than S1 group, and the difference of TOC was significant (p < 0.05). Zn was the highest metal (metalloid) with a content up to 52.9 mg/kg in the sample X. The contents of As, Hg, and Cd were significantly higher in M1 group than that in S1 group (p < 0.05). Figure 5 shows the significant correlations (p < 0.01) between sediment properties and antibiotic resistance variables (ARVs) consisting of the abundances and diversity of ARGs, resistance types and resistance mechanisms. The significant negative correlations between sediment properties and ARVs were found by correlations analysis (p < 0.01).
Metal Cd was the most frequent factor related to ARVs. It negatively correlated to total ARG abundance, abundance of ARGs rpoB, gyrA, and ef-Tu, abundance of resistance type elfamycin, peptide and aminoglycoside, and abundance of resistance mechanism target alteration. The pH and TOC showed significant (p < 0.01) correlation with ARG and resistance type of ARVs. Similar correlations were found in Cu, As, and Hg. No significant (p < 0.01) correlation was found between Zn and ARVs.

| MGES
The integrase genes, ISCR, and plasmids that are potentially related to HGT were annotated, as shown in Table 1 In the MGEs, intIA and plasmid showed significant difference (p < 0.05) between M1 and S1. The average abundance of plasmid in the group S1 was significantly greater (p < 0.05) than group M1, while intIA was different. Figure 5 shows the significant correlations (p < 0.01) between MGEs and ARVs. The abundance of total integrase genes correlated to resistance type sulfonamide/sulfone, ISCR1 and plasmid mainly showed positive correlation with the ARVs, while intIA was showed only negative correlation with the ARVs. No significant (p < 0.01) correlation were found between intI and ARVs.

| Microbial community
To research the effect of microbial community on the distribution of ARGs, a total of 93 classes were identified in the samples.
Groups M1 Groups S1 Sample E F I G U R E 4 The relative abundance of antibiotic resistance genes (ARGs), resistance types and resistance mechanisms which was significant difference between group M1 and S1. And the more abundant variable is marked as 100% To investigate the major change of microbial community between M1 and S1, we chose the most abundant 35 classes from all the 93 classes that accounted for more than 98.96% of total abundance in each sample and represented the dominant classes. Figure 6 shows the abundances of the most abundant 35 classes, which belong to 15 phyla. About half of the classes were more abundant in samples of M1 than in S1, while the others were different.

| Correlation among ARGS, sediment properties and MGES, as well as microbial community
The results of pRDA showed that the microbial community, sediment properties, and MGEs explained 79.95% of the differences in the distribution of ARGs (Figure 7). Sediment properties explained 23.07% alone, MGEs explained 2.87% alone, and they collectively explained 51.22%, while the microbial community contributed 11.50% alone.
For sediment properties and microbial community, the two major factors contributed to the ARG distribution alone, the Mantel test showed a matrix-level significant correlation between the ARGs

| Impact of mangroves as natural wetland on ARGs
The estuarine water flow is complex, and the sediments are periodically soaked by tides, making it difficult to research straightforwardly the impact of mangrove area on ARGs through the analysis of ARG abundance in influents and effluents, unlike in sewage treatment plants and constructed wetlands (Chen et al., 2015;Yang, Li, Zou, Fang, & Zhang, 2014). However, we achieved this aim by comparing F I G U R E 5 Significant correlation between (a) sediment properties and antibiotic resistance variables (ARVs) and (b) mobile genetic elements (MGEs) and ARVs. a RT, resistance type; b RM, resistance mechanism. *p < 0.05; **p < 0.01; ***p < 0.001 ARVs between mangrove and nonmangrove sediments, which were considered to receive identical pollutions in the same estuary.
For the abundance, the most abundant ARGs and the total ARGs, the most abundant resistance types, and the most resistance mechanism in mangrove sediment were significantly lower than that in nonmangrove sediment in this study (p < 0.05) (Figure 4). The results indicated that mangrove ecosystems have the capacity to remove ARGs as a natural wetland. Recently research indicates that antibiotics used in aquaculture were commonly occurred in mangroves, and the antibiotic pollutions can be mitigated with mangrove vegetation to some extent . The lack of selective pressure from antibiotics may contribute to the ARGs reducing, because antibiotic resistance always carries a fitness cost (Vogwill & MacLean, 2015).
Biodegradation is considered to be one of the reasons for removing ARGs in the integrated constructed wetland (Chen et al., 2015), and anaerobic digestion can effectively remove ARGs during swine manure treatment . And we found that the sediment properties were much different between the M1 and S1 groups, such as TOC, As, Hg, and Cd were significantly lower in M1 than that in S1 (Table S4). The sediment properties (pH, TOC, Cu, As, Hg, and Cd) were significantly correlated to ARGs such as mdtB ( (Chen et al., 2015;Fang et al., 2017), while others only observed a slight effect (Berglund et al., 2014;Hsu, Mitsch, Martin, & Lee, 2017).
Regardless, our work may serve as a precedent for ARG research in the mangrove ecosystem.

| Correlations between antibiotic resistance and abiotic factors, as well as biological factors
To further understand the potential factors that influence the antibiotic resistance, we analyzed the relationship between ARVs and abiotic factors (properties of sediment), as well as biological factors (microbial community and MGEs) in the samples.
Metals such as lead, copper, zinc, cadmium, and metalloid arsenic are often used in agriculture and aquaculture (Singer, Shaw, Rhodes, & Hart, 2016). Since resistance to antibiotics and metals can share certain mechanisms (Singer et al., 2016), the phenomenon called coselection between them has been reported, and the stress of widespread metals in the environment may promote this phenomenon (Li, Xia, & Zhang, 2017). Especially, Cu and Zn that were reported to have an extensive positive correlation with antibiotic resistance (Ji et al., 2012;Wang et al., 2016;Zhu et al., 2013) in the samples from livestock farms. However, in the total metal or metalloid elements detected in this study, only As and Hg positively correlated to abundance of arnA.
While the others showed negative correlation to the different ARVs, except zinc has no correlation with ARVs ( Figure 5). The previous study showed that the occurrence of coselection dependent on whether the metal concentration reaches the minimum coselective concentration of an environment (Seiler & Berendonk, 2012 , 1995). And the concentrations of Cu and Zn are much lower than the samples from livestock farms in the above study (Ji et al., 2012;Wang et al., 2016;Zhu et al., 2013). In our study, sediment TOC was found to negatively correlate to gene rpoB, resistance type aminoglycoside and amount of resistance types ( Figure 5), and pH was positively correlated to abundance of gene mdtB and ARG diversity, and negatively correlated to abundance of gene arnA. Correlations between ARGs and pH or TOC in different soil researches are variant (Knapp et al., 2011;Wang et al., 2016), it seems that TOC and pH affect ARGs depends on the specific environmental background.
In this study, results of pRDA showed that the sediment properties and MGEs were the main causes for the differences in the distribution of ARGs, and the sediment properties is the most influential factor for ARG distribution rather than microbial community.
We also analyzed the cooccurrence pattern among ARGs, MGEs, sediment properties, and microbial community using network analysis (Spearman's r > 0.7, p < 0.05). The network was divided into F I G U R E 7 Venn analysis for the contribution of microbial community, sediment properties, and mobile genetic elements to variation of antibiotic resistance gene two major modules and one mini module after modularity analysis ( Figure A2). MGEs including plasmid, ISCR1 and intI1 located at the key nodes of the module one, indicating that these MGEs play an essential role in the module. And the sediment properties such as Hg, TOC, Cd, and As were located at the key nodes of the module two, these sediment factors may affect the ARG distribution in the module. This is similar to the results of studies on Chinese estuaries at a continental scale (Zhu et al., 2017). However, another research on sewage treatment process showed that bacterial community and MGEs are the main factors impacting the fate of ARGs . This suggested that the most influential factors affecting the ARG distribution are different in the natural environment and biological treatment processing, but the MGEs have a substantial impact on ARGs. The network analysis also confirmed that a variety of ARGs and MGEs connected to the same class, this indicated the possibility of the ARGs migrating in the modules through HGT, which may lead to the emergence of antibiotic-resistant strains.

| Profiles of antibiotic resistance in the subtropical urban estuary adjacent to the south china sea
It was reported that ARG abundance range from 0.066 to 0.543 among 18 Chinese estuaries (Zhu et al., 2017). The average ARG abundance in our study is 0.170, which is in the same order of magnitude as a previous survey (Zhu et al., 2017). It seems that this abundance range is the current situation of Chinese estuaries. However, it is much lower than the ARG abundances of samples from urban sewage treatment plants influents and wastewater from livestock farm (Zhu et al., 2013).
Twenty-two antibiotic resistance types were found in this work.
The most abundant three resistance types in each sample are elfamycin, multidrug, and fluoroquinolone ( Figure 2). Elfamycin resistance in the study was mainly contributed by the most abundant ARG ef-Tu which is a prokaryotic elongation factor gene that becomes ARG due to mutation (Cappellano, Monti, Sosio, Donadio, & Sarubbi, 1997;Prezioso, Brown, & Goldberg, 2017). Multidrug resistance is mainly contributed by the second highest abundance ARG rpoB and vast efflux genes. The rpoB is also an evolved ARG from RNA polymerase gene, were widely found in pathogenic mycobacteria (Andre et al., 2017). The majority of efflux pumps are encoded by chromosomes (Blanco et al., 2016) and considered as a general metabolic mechanism (Borges-Walmsley, Mckeegan, & Walmsley, 2003). This coincided with a point of view that lots of ARGs evolved from natively functional genes of the bacterium (Wright, 2007). Fluoroquinolone is the most used antibiotic in China and the fluoroquinolone resistant genes were prevalent in China (Xu et al., 2009;Zhang et al., 2015;Zhang, Lu, Wang, Pang, & Zhao, 2014), similar situation was found in several relatively abundant resistance types such as β-lactam and tetracycline ZHANG et al., 2017ZHANG et al., , 2015Zhu et al., 2017).

| CON CLUS ION
In this study, profiles of ARGs in a mangrove area at an urban estuary in the South China Sea were investigated. The abundance of ARGs was significantly lower in mangrove sediment than that in nonmangrove sediment (p < 0.05), suggesting the ARGs-removed capacity of the subtropical mangrove wetland ecosystems. And sediment properties together with MGEs were the most important factors impacting the distribution of ARGs rather than microbial community.

CO N FLI C T O F I NTE R E S T S
All authors declare no competing interests.

E TH I C S S TATEM ENT
The study was carried out under the stipulation of ethics ap-

DATA ACCE SS I B I LIT Y
All the metagenomic sequencing data analyzed in this work were deposited in the NCBI SRA database under the BioProject PRJNA472209.