Skip to main content
Download PDF

Natural selenium stress influences the changes of antibiotic resistome in seleniferous forest soils

Environmental Microbiome volume 17, Article number: 26 (2022) Cite this article

Abstract

Background

Metal(loid)s can promote the spread and enrichment of antibiotic resistance genes (ARGs) in the environment through a co-selection effect. However, it remains unclear whether exposure of microorganisms to varying concentrations of selenium (Se), an essential but potentially deleterious metal(loid) to living organisms, can influence the migration and distribution of ARGs in forest soils.

Results

Precisely 235 ARGs conferring resistance to seven classes of antibiotics were detected along a Se gradient (0.06–20.65 mg kg−1) across 24 forest soils. (flor)/(chlor)/(am)phenicol resistance genes were the most abundant in all samples. The total abundance of ARGs first increased and then decreased with an elevated available Se content threshold of 0.034 mg kg−1 (P = 2E−05). A structural equation model revealed that the dominant mechanism through which Se indirectly influences the vertical migration of ARGs is by regulating the abundance of the bacterial community. In addition, the methylation of Se (mediated by tehB) and the repairing of DNA damages (mediated by ruvB and recG) were the dominant mechanisms involved in Se resistance in the forest soils. The co-occurrence network analysis revealed a significant correlated cluster between Se-resistance genes, MGEs and ARGs, suggesting the co-transfer potential. Lelliottia amnigena YTB01 isolated from the soil was able to tolerate 50 μg mL−1 ampicillin and 1000 mg kg−1 sodium selenite, and harbored both Se resistant genes and ARGs in the genome.

Conclusions

Our study demonstrated that the spread and enrichment of ARGs are enhanced under moderate Se pressure but inhibited under severe Se pressure in the forest soil (threshold at 0.034 mg kg−1 available Se content). The data generated in this pilot study points to the potential health risk associated with Se contamination and its associated influence on ARGs distribution in soil.

Introduction

Antibiotic resistance and/or multidrug resistance of clinical pathogens has become a serious global health problem, which has been aggravated by the increasing overuse of antibiotics [1,2,3]. Overuse of antibiotics poses selective pressure on bacterial communities and facilitates the dissemination of antibiotic resistance genes (ARGs) in the environment [4]. Indeed, ARGs can be horizontally transferred (commonly termed horizontal gene transfer, HGT) between pathogenic and non-pathogenic bacteria through mobile genetic elements (MGEs) [5]. Metal(loid)s, widely distributed and persistent contaminants of the soil, differ from antibiotics but can serve as a co-selective pressure for antibiotic-resistant bacteria and ARGs [6]. Indeed, several studies have demonstrated that under favorable conditions, metal(loid)s can enhance antibiotic resistance in various soil environments, including multi-metal(loid)s contaminated swine farms soil [7], copper contaminated soil [8], arsenic enriched pig manure [9], as well as arsenic and cadmium polluted paddy soil [10].

Selenium (Se) is an essential nutrient element for many microorganisms and animals, including humans [11]. It exists as different species in the soil, including highly soluble Se oxyanions (selenate and selenite), less mobile elemental Se and reduced selenides [12]. In organisms that need Se, selenocysteine is involved in the formation of selenoproteins and plays an important role in detoxifying and scavenging reactive oxygen species (ROS) [13, 14]. Also, a moderate concentration of Se is beneficial to the growth of bacteria, while a high concentration of Se causes deleterious effects [15]. When bacteria were exposed to toxic concentrations of selenite, intracellular selenite acted as a pro-oxidant, leading to the accumulation of ROS, which in turn caused damage to proteins and DNA [16, 17]. To resist the toxicity of selenite, a series of Se resistance genes, including superoxide dismutase (sodA and sodB) [18], DNA damage repair enzymes (recG and ruvB) [19], and tellurite-resistance genes (tehB) [20] coordinate the detoxification and repair process. Experimental evidence has shown that Se-enriched bacteria have a stronger ability to scavenge ROS and are more resistant to antibiotics than non-Se-enriched bacteria [21]. In addition, analyses of massive complete genome collections found close genetic linkages and high co-transfer potentials between Se resistance genes and bacitracin as well as between kasugamycin and β-lactam resistance genes [22]. Therefore, we hypothesize that Se pressure may affect the distribution of ARGs in soils.

In humans, Se has many important biological functions, including anti-oxidation and anti-cancer activity [23]. Studies have shown that some populations living in soil Se deficient areas, including China, Europe and New Zealand, lack adequate Se intake and, therefore, are at risk of Se deficiency [24]. Developing Se-enrich agricultural foods have been recognized as an effective strategy to argument the Se intake of humans [25, 26]. The application of selenite fertilizers in agriculture or the addition of selenite feeds in animal husbandry is gaining popularity [27]. However, overuse of selenite fertilizers may contaminate the soil. Since metal(loid)s have been found to drive the spread of ARGs in soil ecosystems [7,8,9,10], it remains unclear whether the exposure of microorganisms to Se pressure may promote the enrichment and migration of ARGs in soils.

Forests are an important component of terrestrial ecosystems globally and a huge ARGs reservoir for exchange with clinical pathogens [28]. The Enshi Autonomous Prefecture of Hubei province, China, is a natural Se-rich area where the average soil Se content reaches 9.6 mg kg−1 (the world average Se content is 0.2 mg kg−1) [29]. The Yutangba area, located in the north of Enshi Autonomous Prefecture, is a typical high-Se basin, with soil Se content ranging from 0.41 to 42.3 mg kg−1 [29]. Seleniferous forest soils with low anthropogenic influences are ideal areas to investigate the relationship between Se exposure and the distribution of ARGs as well as the underlying mechanisms. Therefore, we aimed to investigate the patterns of ARGs and bacterial community along a Se gradient by analyzing soils from 24 forest sites in Enshi. Furthermore, we used a culture-dependent approach to investigate the association between Se-enrichment and antibiotic resistance in bacteria. Lastly, we investigated the genetic mechanisms of the association between Se-enrichment and antibiotic resistance by analyzing the whole genome of Se-enriched bacteria.

Methods

Study area and soil sampling

Soil samples were collected in the forest of Enshi, an autonomous prefecture located in the southwest mountain area of Hubei Province, China (see Fig. 1 and Additional file 7: Table S1 for geographic information on the area). The sampling area is a subtropical monsoon climate with an annual average temperature of 16.2 °C and annual average precipitation of 1600 mm [30]. Based on previous studies which fully investigated the soil Se content in our sampling region, we selected 24 sampling sites around the Se-rich Yutangba area (Fig. 1). To the best of our knowledge, animal manure or fertilizers have never been applied to the soils in these sampling sites. Also, prior to sampling, these sites have experienced very minimal human impacts. In each site, three replicate samples were obtained randomly from a 20 m × 20 m plot. After removing the litter layer, samples were aseptically collected from the topsoil (10 cm depth) using a sterile spade. Three soil sub-samples were collected from 1 m × 1 m sub-plot and mixed to form a composite sample. Composite samples were immediately placed on ice and transported to the laboratory within 24 h of collection. Soil sub-samples for the isolation of culturable bacteria were stored at 4 °C, and soil sub-samples for the extraction of total microbial community DNA were freeze-dried and stored at − 80 °C. The soil sub-samples for the analyses of chemical properties were air-dried, crumbled to pass through a 2-mm sieve and thoroughly homogenized using standard methods.

Fig. 1
figure 1

Map of the study area. Twenty-four sampling sites around the high Se content area (red point) in the forest of Enshi Autonomous Prefecture of Hubei province, China, were sampled

Full size image

Chemical properties analyses and Se determination

The soil water content was determined gravimetrically by drying the soil samples at 105 °C for 24 h. Soil pH was determined with a soil to water ratio of 1: 2.5 using a pre-calibrated FiveEasy Plus™ pH-meter (Mettler Toledo, Zurich, Switzerland). The total soil carbon (TSC) and total soil nitrogen (TSN) were determined on an automated VarioMax CN Element Analyzer (Elementar, Hanau, Germany). Soil total phosphorus was measured with Mo-Sb colorimetric method using a spectrophotometer (Shimadzu, Kyoto, Japan). Soil total potassium was digested with hydrofluoric acid and perchloric acid and subsequently measured using a flame photometer (Shanghai Jing Ke Instrument Co., Ltd, Shanghai, China).

Se contents in soil were determined using HG-AFS (coupling of hydride generation to atomic fluorescence spectrometer), a sensitive method for the determination of Se [31]. Preprocessing of samples and determination of total and available Se contents were conducted according to the method of García et al. [31] with slight modifications (Additional file 8: Text S1). Se content in the samples was determined by extrapolation to a standard curve of Se concentrations analyzed using a double-channel atomic fluorescence spectrometry (AFS-930, Beijing Titan Instruments Co., Ltd, Beijing, China). The recovery rate of soil quality standard GBW07428(GSS-14) was 90–110%.

Antibiotic residue quantitation

Two target tetracyclines, namely tetracycline (TCN) and doxycycline hydrochloride (DOX); two (fluoro)quinolones, namely ofloxacin (OLF) and enrofloxacin (ENR); two β-lactam, namely amoxicillin (AMX) and clavulanic acid (CA); the aminoglycoside kanamycin (KAN); and the erythromycin (ENY), a macrolide lincosamide streptogramin_b (MLSB), were quantified in this study. The antibiotic residues in soils were extracted using solid-phase extraction [32]. Thereafter, the concentrations of antibiotics were determined using liquid chromatography-tandem mass spectrometry (LC-MS/MS) as previously described [33].

Determination of the profile and abundances of ARGs and MGEs through high-throughput quantitative PCR (HT-qPCR)

High-throughput quantitative PCR (HT-qPCR) was used to determine the profile and abundance of ARGs and MGEs in soil samples. The profile of MGEs was targeted in order to measure the potential of horizontal gene transfer of ARGs. Firstly, microbial DNA was extracted from 0.25 g soil sample with the MoBio PowerSoil DNA isolation kit (Qiagen, Hilden, Germany) following the manufacturer's instruction. The quantity and quality of the extracted DNA were determined using a NanoDrop ND2000c spectrophotometer (NanoDrop Technologies, Wilmington, DE, USA). The DNA was kept at − 20 °C until analyses.

Next, the profile of ARGs and MGEs were measured using HT-qPCR array on a CFX384™ Real-Time system (Bio-Rad Laboratories, Hercules, CA, USA) [34]. The HT-qPCR array contains 364 pairs of primers which can simultaneously target 328 ARGs including seven main classes of antibiotics (aminoglycosides, β-lactam, erythromycin, (flor)/(chlor)/(am)phenicol, MLSB, tetracycline, and vancomycin), and series of multidrug/efflux genes; 34 MGEs, including transposon-transposase genes, integron-integrase genes and plasmids; and 2 replicate 16S rRNA genes (housekeeping genes). Each 10 μL qPCR reaction contained 5 μL of Premix Ex Taq™ (TaKaRa Biotechnology, Dalian, China), 0.5 μL of each primer (10 mM), 1 μL template DNA, and nuclease-free water. The thermal cycling conditions were: 95 °C for 10 min and then 40 cycles of denaturation at 95 °C for 30 s and annealing at 60 °C for 30 s. Each reaction was conducted in triplicate and contained negative controls. The qPCR amplification efficiency ranged from 80 to 120%. The relative abundance of ARGs and MGEs was determined after dividing target gene abundance by the bacterial 16S rRNA gene copy number [35].

The quantitative PCR of the 16S rRNA gene was conducted on a CFX96™ Real-Time system (Bio-Rad Laboratories) by standard curve method [36]. The primer set was Bact1369F/Prok1492R with the ProbeTM1389F [37]. The reaction wells contained 12.5 μL Premix Ex Taq™ (Takara Biotechnology, Dalian, China), 10 μM of each primer, 5 ng μL−1 template DNA, 1 μL ProbeTM1389F and nuclease-free water to a final reaction volume of 25 μL. Information on the plasmid DNA used to generate the standard curve and the amplification procedure are described in a previous study [38].

Absolute quantification of Se resistance genes

Se resistance genes, including recG, ruvB, sodA, sodB and tehB were retrieved from BacMet (http://bacmet.biomedicine.gu.se/). The genes recG and ruvB are ATP-dependent DNA helicases involved in repairing DNA damage caused by selenite [19], while sodA and sodB are Mn- and Fe-dependent superoxide dismutases, respectively [18]. The products of sodA and sodB are involved in eliminating ROS caused by excessive intracellular selenium. The tehB gene is a S-adenosyl-L-methionine-dependent methyltransferase involved in the methylation of Se and responsible for Se detoxification [20]. In addition, based on published research, we targeted eleven Se metabolism genes, including two methyltransferase genes mmtA and ubiE [39], two selenocysteine synthesis genes SelA and SelB [40], two groups of selenate reductase: serA, serB, serC and serD (group 1) [41] and srdA, srdB, srdC (group 2) [42]. Libraries for all the above sixteen genes were constructed via PCR cloning and generation of recombinant T-Vectors harboring the 16 genes. Thereafter, these plasmid DNAs were used to generate standard curves and amplification profiles for the absolute quantification of the abundance of the sixteen genes in soil samples on a CFX96™ Real-Time system (Bio-Rad). The qPCR was a 10 µL reaction volume and included 5 µL Premix Ex Taq™ (Takara Biotechnology, Dalian, China), 0.5 µL each primer, 0.5 µL template DNA, and nuclease-free water. The primer sequences and PCR conditions are provided in Additional file 7: Table S2. Each PCR reaction was performed in triplicate and included negative controls. The amplification efficiency was valid (ranged from 80 to 120%). The relative abundance of Se resistant genes was determined after dividing target gene abundance by the bacterial 16S rRNA gene copies.

Bacterial 16S rRNA gene sequence-based diversity analysis

The bacterial community composition was determined by amplifying the hypervariable V4 region of the bacterial 16S rRNA gene with primer pair 515F and 907R [43]. PCR reactions were performed on GeneAmp 9700 thermal cycler (Applied Biosystems, Waltham, MA, USA) and conducted in triplicate. Each PCR reaction contained 5 μL of Premix Ex Taq™ (TaKaRa Biotechnology, Dalian, China), 0.5 μL of each primer, 10 ng template DNA and nuclease-free water in total 10 μL volume. The thermocycler condition was 94 °C for 3 min, 27 cycles of 95 °C for 30 s, 55 °C for 30 s and 72 °C for 45 s at, with a final extension at 72 °C for 10 min. The PCR products were purified using the AxyPrep DNA Gel Extraction Kit (Axygen Biosciences, Union City, CA, USA). The sequencing was performed on an Illumina MiSeq system (Illumina, San Diego, CA, USA) at Majorbio Bio-Pharm Technology, Shanghai, China.

Raw sequence reads were processed using the Quantitative Insights Into Microbial Ecology (QIIME, version 1.9) with the suggested operational procedure [44, 45]. Usearch v10.0 [46] was used to merge quality-filtered forward and reverse reads as well as eliminate chimeric sequences. Thereafter, the sequence reads were clustered into operational taxonomic units (OTUs) at 97% 16S rRNA gene similarity, and representative sequences aligned against the Silva reference database (Release 132) for taxonomic assignment (https://www.arb-silva.de/). All 16S rRNA gene sequences have been uploaded to Bio-Med Big Data Center (https://www.biosino.org/node/index) with the accession number OEP001738.

Isolation of Se-resistant bacteria and growth kinetics

To investigate the association between a high concentration of Se and resistance genes in bacteria, we selected two soil samples that were respectively collected from sites with the highest and lowest soil selenium contents. The soil samples were serially diluted (up to 10−3) with normal saline (0.85% NaCl) and inoculated onto tryptic soy broth (TSB) (1:200), supplemented with 50 μg mL−1 ampicillin and 0, 1, 10, 100 or 1000 mg kg−1 sodium selenite. The ampicillin concentration and the gradient of sodium selenite concentration were determined according to preliminary experiments. That is, the concentrations of 50 μg mL−1 ampicillin and 1000 mg kg−1 sodium selenite inhibited the growth of most sensitive bacteria as determined by enumeration of colony forming unit. Each ampicillin-sodium selenite combination was set up in triplicates. Cultures were incubated at 28 °C for 48 h with gently shaking at 50 rpm. Bacterial growth was monitored every 6 h by optical density measurements at a wavelength of 600 nm. Subsequently, for enumeration of viable bacteria, cultures in the middle and early stages of the logarithmic growth phase were selected and spread plated on tryptic soy agar (TSA) plates supplemented with 3000 mg kg−1 sodium selenite. The cultures were then inoculated in TSB supplemented with a series of sodium selenite (0, 1, 10, 100 or 1,000 mg kg−1) and 50 μg mL−1 ampicillin to generate a growth curve. Furthermore, to isolate Se-resistant bacteria, three single clones were picked up from each of the above TSA plates and then cultured in TSB. Finally, sodium selenite susceptibility assays were performed on the pure cultures by growing them in TSB amended with 50 μg mL−1 ampicillin and 0, 500, 1000, 3000, 6000, or 9000 mg kg−1 sodium selenite at 28 °C with shaking for 48 h. Thereafter, the optical density was determined at 600 nm with a spectrophotometer.

Genome sequencing, assembly and gene annotation

The DNA of the pure microbial cultures was extracted using the QIAamp DNA Microbiome Kit (Qiagen, Hilden, Germany) by following the manufacturer’s instructions. Thereafter, the V4 region of the 16S rRNA gene was sequenced using Sanger sequencing [43]. The complete genome sequencing of the isolates was performed on the PacBio RS II platform (Pacific Biosciences of California, Inc, Menlo Park, CA, USA) and Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, USA) at the Beijing Genomics Institute (BGI, Shenzhen, China).

The sequences of the V4 region of the 16S rRNA gene were manually checked for quality and edited using Chromas software, and aligned against the NCBI GenBank using BLAST to identify the bacterial species. The complete genome sequence data was assembled using multiple software, including Falcon v0.3.0, proovread v2.12 and Celera Assembler v8.3 as previously reported [47,48,49]. The genome gene prediction was performed using glimmer3 (http://www.cbcb.umd.edu/software/glimmer/) [50, 51]. The genome gene was annotated through the NCBI PGAP using BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Seven databases, including Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.genome.jp/kegg/), Clusters of Orthologous Groups (COG, https://www.ncbi.nlm.nih.gov/research/cog-project/), Non-Redundant Protein Database (NR) [52], Swiss-Prot [53], TrEMBL [53], Gene Ontology (GO) [54], and EggNOG v5.0 [55] were used for general functional annotation. Virulence factors and drug resistance genes were identified using the Virulence Factors of Pathogenic Bacteria database (VFDB, http://www.mgc.ac.cn/VFs/main.htm) and The Comprehensive Antibiotic Resistance Genes Database (CARD, https://card.mcmaster.ca/) [56], respectively.

To obtain insights into whether co-regulation between Se metabolism and antibiotic resistance exists, the whole genome sequence information was further mined for the presence of secondary metabolite biosynthetic gene clusters (BGCs) using antiSMASH v6.0 [57]. The complete genome sequences are available in the Bio-Med Big Data Center (https://www.biosino.org/node/index) under the accession number OEP001741.

Statistical analysis

Pearson correlation and normality test were performed in SPSS 21 (IBM, Armonk, NY, USA). Five models, namely linear, quadratic, cubic, segmented, and stegmented regression, were used to fit the relationship between available Se content and the detected number of ARGs or relative abundance of ARGs with the “chngpt” package in R (v.3.4.1) (https://www.r-project.org/) [58]. Redundancy analysis (RDA) was performed to investigate the relationship of soil chemical properties, bacterial phylum, relative abundance of MGEs and Se resistance genes on the pattern of ARGs using the “vegan” package of R [59]. Structural equation modeling (SEM) was used to examine the direct and indirect effect of edaphic factors, Se resistance genes, bacterial abundance and MGEs on the pattern of ARGs. SEM was constructed and analyzed using the “lavaan” [60] and “piecewiseSEM” packages [61] of R, respectively. For the SEM, the first principal coordinate of the relative abundance matrix of Se resistant genes, MGEs, ARGs and phylum-level abundance of bacterial community were extracted with the “vegan” package of R. The goodness of fit was tested by “fitMeasures” function in “lavaan” package. Co-occurrence network analysis was performed to explore the interactive correlation among the relative abundance of ARGs, MGEs, and Se resistance genes (r > 0.8, P < 0.05). The pairwise correlations were calculated by Pearson correlation analysis using the "psych" package [62] in R and then visualized and further annotated using Gephi v0.9.2 [63].

Results

Se contents and major soil characteristics

Among all samples, soil total Se content ranged from 0.06 to 20.65 mg kg−1, with a mean value of 3.40 mg kg−1 (Additional file 7: Table S1), which was higher than the world mean total Se value of 0.2 mg kg−1 [29]. The highest total Se content sample, with a mean content of 20.65 mg kg−1, was collected from TZY2 site. Soil available Se content ranged from 0.003 to 0.14 mg kg−1 with a mean of 0.03 mg kg−1. Pearson correlation analysis showed that pH and available Se content were significantly positively associated (r = 0.61, P < 0.05) (Additional file 7: Table S3). In addition, available Se and total Se contents were significantly positively correlated (r = 0.90, P < 0.05). Given that available Se is the bioavailable form of Se in the soil [12], it was selected for further analyses instead of total Se content.

Diversity and abundance of ARGs and MGEs

A total of 264 genes, including 235 ARGs and 29 MGEs, were detected in all samples (Additional file 1: Fig. S1). The numbers of detected ARGs in each sample ranged from 7 to 84 (Fig. 2a). The detected ARGs confer resistance to seven classes of antibiotics, including aminoglycoside, β-lactam, (flor)/(chlor)/(am)phenicol, sulfonamide, tetracycline, MLSB, and vancomycin, and also codes for a series of multidrug/efflux proteins. The detected number of ARGs coding for multidrug/efflux pumps was the highest, accounting for 18.9% of the total numbers of ARGs. ARGs for β-lactam (17.8%) and tetracycline (17.2%) were the next most abundant after multidrug/efflux pumps. The relative abundance of ARGs detected in these seleniferous soil samples ranged from 2.30E−06 to 9.25E−03 copies per 16S rRNA gene copy (Fig. 2b). Among them, multidrug resistance genes were detected in the highest abundance (63.2% of total ARGs) in all samples. Among the 235 ARGs, 11 highly abundant genes were detected in as much as 79.2% (19/24) of samples, including five multidrug/efflux genes (oprJ, mexF, acrR-03, acrA-04, acrA-05), two β-lactam resistance genes (oprD, fox5), one MLSB resistant gene (oleC), one aminoglycoside resistant gene (aadA1), one tetracycline (tetPB-01), and one other type of ARGs (sat4). MexF, an important multidrug efflux gene, had the highest abundance with a mean relative abundance of 5E−03 copies per 16S rRNA copy across all samples. The normalized copy numbers of MGEs ranged from 9.79E−06 to 1.32E−01 copies per bacterial cell. The intI1 gene encoding class 1 integron-integrase was detected in only 3 of the 24 samples, with a mean of 5.52E−05 copies per bacterial cell. In addition, Pearson correlation analysis showed that the total abundance of MGEs was not correlated with the total abundance of ARGs in each sample (data not shown).

Fig. 2
figure 2

ARGs abundance and correlation with Se content. a The number of detected ARGs. b The relative abundance of ARGs. c Correlation between available Se content and the detected number of ARGs. d Correlation between available Se content and the relative abundance of ARGs (d). In a, c, ARGs conferring resistance to seven classes of antibiotics or multidrug/efflux were shown in columns of different colors. Others refer to ARGs are not included in the seven classes of antibiotics. Blue solid lines are segmented according to the threshold regression model and show the significant relationship below and above the threshold (b, d)

Full size image

Concentration of antibiotics residues

To investigate whether the ARGs pattern was affected by antibiotics residues in soil, the concentrations of eight major classes of antibiotics residues were analyzed. As shown in Additional file 7: Table S4, the concentrations of amoxicillin (β-lactam class) and tetracycline (tetracycline class) were between 0 and 12.42 μg kg−1 and between 0 and 0.45 μg kg−1, respectively. The concentrations of the other six antibiotics were below the detection limit in all samples. Furthermore, the Pearson correlation analysis showed that the relationship between the concentration of detected antibiotics and the diversity or abundance of ARGs was not significant (P > 0.05).

16S rRNA gene quantification and bacterial diversity

The absolute copy number of 16S rRNA genes (for the total abundance of the bacterial community) ranged from 4.00E+09 to 6.26E+10 copies g−1 soil. Based on a second-order polynomial fit, we found that the absolute copy number of the 16S rRNA genes exhibited a unimodal relationship with available Se content (R2 = 0.312, P < 0.001) (Fig. 3a) and peaked at 0.08 mg kg−1 available Se content.

Fig. 3
figure 3

The correlation between available Se content and the absolute copy number of genes. a 16S rRNA. b ruvB. c tehB. d recG. Significant quadratic correlations were well fitted between available Se content and the abundance of 16S rRNA and tehB. Significant positive correlations were found between available Se content and the abundance of ruvB and recG

Full size image

From high-throughput amplicon sequence analysis of the V4 hypervariable region of the 16S rRNA gene, a total of 1,486,320 high-quality sequences (50,006 to 71,070 per sample) were obtained and clustered into 1,178,016 OTUs (38,072 to 58,547 for each sample; mean = 49,533). Proteobacteria, Acidobacteria, Actinobacteria and Planctomycetes were the dominant phyla in all the samples, accounting for 11.76% to 33.39% of total bacterial 16S rRNA gene sequences (Additional file 2: Fig. S2). At the class taxonomic level, Acidobacteria, Alphaproteobacteria, and Actinobacteria were the dominant bacteria, accounting for 22.94%, 17.79%, 12.58% of the total 16S rRNA gene sequences, respectively.

Quantitative analysis of sixteen Se resistance genes

To investigate how Se resistance genes impact Se distribution and ARGs transfer in these forest soil under Se selective pressure, we constructed libraries of sixteen Se resistance genes and then performed absolute quantitative PCR to determine their abundance in soils across a Se content gradient. Amplifications were obtained for only 12 of the 16 genes, and these were further analyzed. The absolute copy number of Se resistance genes varied from 1.33E+02 copies g−1 soil to 5.40E+08 copies g−1 soil (Additional file 7: Table S5). Of the 12 genes, ruvB was the most abundant, with a mean copy number of 1.93E+08 per gram of soil. Interestingly, the relative abundance of ruvB, tehB, and recG were significantly correlated with available Se content (Fig. 3b–d). Specifically, the relative abundance of tehB exhibited a significant quadratic relationship with available Se content (R2 = 0.31, P = 0.02) while the copy numbers of ruvB and recG exhibited a positive correlation with available Se content (ruvB: R2 = 0.24, P < 0.01; recG: R2 = 0.23, P < 0.01). More specifically, a threshold for available Se content existed at 0.08 mg kg−1, at which point, the relative abundance of tehB exhibited a unimodal relationship. The abundance of tehB first increased with elevated available Se ranging from 0.003 to 0.008 mg kg−1 and then decreased with available Se ranging from 0.08 to 0.14 mg kg−1.

Network analysis

The co-occurrence pattern between ARGs, Se resistance genes, and MGEs was explored by network analysis (Fig. 4 and Additional file 3: Fig. S3). The modularity index of the co-occurrence network is 0.797, indicating a high modularity of the network. The network contained six modules, with tightly connected nodes within each module and few cross-connections between different modules. Module I exhibited a close and independent cluster, which consisted of five Se resistance genes, three ARGs, one transposon, and one resistant plasmid. The other five modules were composed of ARGs, Se resistance genes, or MGEs (Fig. 4 and Additional file 3: Fig. S3).

Fig. 4
figure 4

Co-occurrence network analysis showing the correlation between ARGs, MGEs, and Se resistance genes (Spearman' r > 0.8, P < 0.05). Node size indicates connectivity, and the color of nodes is grouped by resistance type

Full size image

Correlation between environmental factors, Se resistance genes, MGEs, bacterial community and ARGs

The minimum Akaike information criterion (AIC) value showed that the stegmented regression model was the best fit, and that the thresholds for available Se content were 0.015 mg kg−1 (available Se versus the detected number of ARGs, P = 8E−04) (Fig. 2c) and 0.034 mg kg−1 (available Se versus the abundance of ARGs, P = 2E−05) (Fig. 2d). The number of detected ARGs increased with increasing available Se content from 0.003 to 0.015 mg kg−1 (R2 = 0.38, P = 0.04) and decreased with available Se content from 0.015 to 0.14 mg kg−1 (R2 = 0.43, P = 0.01). The abundance of ARGs first increased with elevated available Se content ranging from 0.003 to 0.034 mg kg−1 (R2 = 0.67, P < 0.001) and then decreased with available Se content ranging from 0.034 to 0.14 mg kg−1 (R2 = 0.47, P = 0.03).

Redundancy analysis (RDA) was used to examine the impact of environmental factors, Se resistance genes, MGEs, and bacterial community on the patterns of ARGs (Fig. 5). The RDA model explained 54% of the total variance of ARGs. Environmental factors, including the available Se content and C/N ratio, significantly affected the patterns of ARGs in the first and second canonical axes. The relative abundance of Se resistance genes, including ruvB and srdB, also significantly affected the patterns of ARGs. Actinobacteria and Chloroflexi were the two bacterial phyla that had significant effects on the pattern of ARGs. A few MGEs including intI2-02, tnpA-05/IS26 and tnpA-07/ISEcp1 had significant (P < 0.05) effects on the composition of ARGs.

Fig. 5
figure 5

Redundancy analysis (RDA) of the significant soil chemical properties, bacterial phylum, MGEs, and Se resistance genes on the pattern of ARGs. C/N refers to the ratio of total soil carbon and total soil nitrogen. ASe refers to available Se content

Full size image

Structural equation modeling (SEM) was established to integrate direct and indirect influences of variables on ARGs patterns (Fig. 6a, b). A Chi-square P value of 0.702 was obtained for the SEM, and it indicated that the model described the actual relationship among the variables (P > 0.05). The goodness of fit test showed that root mean square error of approximation (RMSEA) was 0.075, and adjusted goodness-of-fit index (AGFI) was 0.96. As shown in Fig. 6a, the overall variables explained 41% of the total variance of ARGs pattern. The model suggested that of the three edaphic variables (available Se, pH and C/N), available Se, with the highest standardized total effect of 0.299 (compared to 0.045) and -0.058 for pH and C/N, respectively), was the primary driving factor of the ARGs patterns (Fig. 6b). Furthermore, available Se content had a significant direct effect on ARGs patterns and also affected ARGs pattern through the indirect effects of Se resistance genes and bacterial abundance.

Fig. 6
figure 6

Structural equation models (SEM) showing relationships between chemical and biological parameters. a Direct and indirect effect of soil properties, selenium resistance genes, bacterial abundance, and MGEs on the pattern of ARGs. b The standardized total effects on ARGs R2 values are the explained proportion of variance. Significance levels of each path are indicated: *P < 0.05, **P < 0.01, ***P < 0.001. The hypothetical models fit the data well and were examined by CMIN/DF = 1.475, RMSEA = 0.075, AGFI = 0.96

Full size image

Isolation and cultivation of Se resistant and antibiotic-resistant bacteria

To further verify the correlation between Se content and ARGs, we used culture-based sodium selenite treatment to select Se resistant and antibiotic-resistant bacteria. The growth of total culturable bacteria was inhibited by sodium selenite (Fig. 7a, b). By comparison, cultures originating from the highest Se content site had a stronger ability to tolerate sodium selenite than those from the lowest Se content site. In addition, Se resistant and antibiotic-resistant bacteria were isolated on TSA plates amended with 50 μg mL−1 ampicillin and 3000 mg kg−1 sodium selenite from the highest and the lowest selenium content sites.

Fig. 7
figure 7

Effects of sodium selenite on the growth of total cultured bacteria originating from high Se content site (a) and low Se content site (b) over 48 h. Effects of different concentrations of sodium selenite on the growth of pure cultures after cultured for 48 h (c). Each value is the mean ± SD of three replicates. The number of putative ARGs in the genome of P. putida (low selenite resistance) and L. amnigena YTB01 (high selenite resistance) (d). The colonies of P. putida and L. amnigena YTB01 cultured for 48 h on TSB supplemented with 3000 mg kg−1 sodium selenite (e)

Full size image

Subsequently, three colonies were selected from each of the two plates and were subsequently identified (using phylogenetic analysis of 16S rRNA gene sequences) as two bacteria, namely Lelliottia amnigena YTB01 (from highest selenium content site) and Pseudomonas putida (from lowest selenium content site). The sodium selenite susceptibility assay in TSB showed that sodium selenite had a dual effect on the growth of Lelliottia amnigena YTB01 and Pseudomonas putida. Specifically, the growth of L. amnigena YTB01 was promoted by sodium selenite concentration of up to 9000 mg kg−1. In contrast, the concentration of sodium selenite below 500 mg kg−1 promoted the growth of Pseudomonas putida, while sodium selenite concentrations higher than 500 mg kg−1 inhibited the growth of Pseudomonas putida (Fig. 7c). The two isolates grew on TSA supplemented with 3000 mg kg−1 sodium selenite (Fig. 7e). The colonies of Lelliottia amnigena YTB01 on TSA were bigger in size and much more in numbers than those of Pseudomonas putida. These results demonstrated that the growth of total culturable bacteria was inhibited by sodium selenite and suggested that the two strains differed in their abilities to tolerate sodium selenite.

The genome analyses of Se-resistant bacteria

Next, we sequenced the whole genome of the Se resistant Lelliottia amnigena YTB01 and Pseudomonas putida in order to investigate the co-occurrence pattern of Se resistance genes and ARGs. The genome of Lelliottia amnigena YTB01 is 4,941,692 bp in size and has a G+C content of 54.43%, including 4605 coding genes and 141 non-coding RNA (ncRNA) (Additional file 4: Fig. S4a). The genome annotation revealed various putative proteins involved in Se resistance and antibiotic resistance (Additional file 5: Material 1), including multidrug resistance protein EmrD, fosfomycin resistance protein FosA, tellurite resistance proteins TehA/TehB, and selenium resistant proteins RuvA/RuvB. The genome of Pseudomonas putida is 6,287,193 bp in size and has a G+C content of 62.34%, including 5,624 coding genes and 139 ncRNA (Additional file 4: Fig. S4b and Additional file 6: Material 2). Importantly, the number of ARGs and Se resistance genes in the genome of Lelliottia amnigena YTB01 were 38 and 11, respectively, which were more than the 25 and 9, respectively, in the genome of P. putida. This observation is striking given that the complete genome of Lelliottia amnigena YTB01 isolated from the highest Se content site is smaller than the genome of Pseudomonas putida isolated from the lowest Se content site (Fig. 7d). Furthermore, in the genome of Lelliottia amnigena YTB01, a hypothetical antibiotic resistance RND efflux gene was identified, and an RNA polymerase was predicted as a regulatory gene (Fig. 8a). In the genome of Pseudomonas putida, two gene clusters both harboring resistance and regulator genes were predicted. In one cluster, a multidrug efflux RND transporter was predicted as a resistance gene downstream, followed by a DNA-binding response regulator (Fig. 8b). In another cluster, a macrolide ABC transporter was predicted as a resistance gene while five regulatory genes were detected (Fig. 8c).

Fig. 8
figure 8

Schematic representation of the important structures related to the regulation of antibiotics resistance genes expression in two typical secondary metabolism BGCs extracted from the genomes of the two isolates. a In the genome of Lelliottia amnigena isolation from the highest Se content site. b and c In the genome of Pseudomonas putida isolation from the lowest Se content site

Full size image

Discussion

Soil Se induces the dual effects of ARGs pattern

Soil is a vast reservoir of ARGs, and the migration of ARGs in the natural environment might pose a public health threat [64, 65]. In our study, we detected a high diversity and abundance of ARGs in seleniferous forest soils and found a significant dual effect between available Se and the diversity or the abundance of ARGs with available Se thresholds of 0.015 mg kg−1 and 0.034 mg kg−1. These results indicated that the ARGs pattern are different below and above the thresholds. Interestingly, the two thresholds are not consistent, suggesting the diversity and abundance of ARGs respond differently to Se stress. Considering that the abundance of ARGs is more relevant for transferring ARGs between non-pathogenic and pathogenic bacteria than the diversity of ARGs, the threshold of 0.034 mg kg−1 available Se is further discussed in this study. The role of metal(loid)s as co-selective pressures driving ARGs distribution has been widely demonstrated [66,67,68]. The present study showed that Se can induce a dual effect on ARGs pattern, unlike the positive correlation between metal(loid)s and ARGs reported in several previous studies [10, 69]. A possible explanation is the varied co-occurrence mechanisms between different metal(loid)s and ARGs. In addition, the concentration of the only two detected antibiotics is quite low and is not correlated with the abundance of ARGs, indicating the weak influence of anthropogenic activities in this area. Therefore, we speculate that antibiotics residues are not the drivers of ARGs changes in this area. Altogether, the data regarding an association between available Se and ARGs demonstrate that the diversity and abundance of ARGs might be increased in moderate Se content soils but decreased in high Se content soils. Given that the available Se content of most soils in the world is below the threshold of 0.034 mg kg−1, the elevated soil Se content caused by Se-enrich agricultural farming can potentially lead to the enrichment of ARGs.

Bacterial community is an important contributor to the variance of ARGs

The bacterial community is the carrier of ARGs and is widely considered to be an important contributor to the variance of ARGs [70]. In this study, we found that available Se had a dual effect on the abundance of the bacterial community, and this effect was consistent with the correlation between available Se content and the abundance of ARGs. This observation indicates that bacterial communities may play a distinct and important role in influencing the ARGs pattern, which was further corroborated by SEM. The explanation for the observed dual effect may be that moderate content of Se is essential and beneficial for bacterial proliferation [71], which can lead to an increase in bacterial abundance, thereby enriching their associated ARGs. Previous evidence from culturable bacteria also suggested that Se-enriched bacteria are more resistant to antibiotics than those not treated with Se [21]. However, excessive Se, especially soluble and toxic selenite, can inhibit the growth of bacteria [13, 71] and thus reduce the corresponding ARGs’ abundance.

Further evidence showed that the dual effect of available Se on bacterial abundance was comprehensively influenced by ruvB, recG and tehB. Hence, it was more likely that the toxicity of Se to bacterial cells will enhance with the increase of available Se concentration, and the abundance of DNA damage repair genes such as ruvB and recG will raise [19]. When the available Se concentration is lower than the threshold (0.08 mg kg−1), the abundance of tehB gene will increase and involves Se methylation. The methylation is also able to detoxify Se [20], which contributes to cell proliferation, thus increasing bacterial abundance. On the contrary, when the available Se is higher than the threshold, the abundance of tehB decreases, which weakens Se-detoxication via Se methylation. Interestingly, complete genome analyses verified the presence of both DNA repair (ruvB and recG) and Se methylation (tehB) genes in the genome of both Lelliottia amnigena YTB01 and Pseudomonas putida.

In this study, HGT may not be the main mechanism for the transfer of ARGs in the forest soils. First, Intl1, a gene responsible for HGT or acquisition of ARGs in environmental microorganisms [28], had a low frequency and abundance in this study. Interestingly, this is consistent with previous reports that the abundance of intl1 is relatively low in environments with minimal anthropogenic impact [5, 72]. Second, a weak correlation between the total abundance of MGEs and ARGs suggests that HGT is not the main mechanism contributing to ARG transfer. Rather, the fact that Se may indirectly mediate the migration of ARGs by regulating the abundance of bacterial communities, indicated that vertical genetic transfer is the dominant mechanism of ARG migration in forest soils.

Co-selection of Se resistance genes and ARGs

An increasing number of studies have demonstrated that metal(loid)s contamination may promote the spread of ARGs through co-selection effect [8,9,10]. Importantly, ruvB (involved in the repairing of DNA damage) and srdB (selenate reductase subunit B) had a significant effect on the pattern of ARGs, suggesting that Se pressure may influence the co-selection of Se resistance genes and ARGs. The co-selection effect is due to the coupling of resistance mechanisms against metal(loid)s and antibiotics [73, 74]. Co-occurrence network analysis revealed a close correlated cluster among ARGs, Se resistance genes, and MGEs, suggesting that ARGs and Se resistance genes may co-occur on the same genetic element or is co-present within the same cell.

Another potential mechanism for Se detoxication is the reduction of Se oxyanions to elemental Se. For example, bacteria can reduce soluble and higher toxic Se (IV) and Se (VI) oxyanions to less toxic Se (0) and bind the elemental Se to the selenate reductase A protein to form a Se nanosphere, which is subsequently exported to the extracellular environment [75, 76]. This was also observed in the culturable bacteria experiment where the two isolates transformed colorless sodium selenite to red elemental Se. The active efflux of the elemental Se may be accompanied by the efflux of antibiotics from the cell to the extracellular environment. This potential co-efflux of Se and antibiotics from the cell is considered a cross-resistance mechanism [77]. However, further studies are required to investigate the hypothesis.

Using culture-based methods, we confirmed that the growth of total cultured bacteria was inhibited by sodium selenite. This observation was in contrast to the observed correlation between available Se and the abundance of the bacterial community (Fig. 3a). The reason might be that the cultured bacteria only accounted for less than 1% of the total bacterial community, and 16S rRNA gene analysis cannot differentiate between dead and alive bacteria. Moreover, the total bacterial cultures from the higher soil Se content site had a stronger ability to tolerate sodium selenite than those from the lower soil Se content site, and this observation was consistent with the results of pure cultures. In addition, the numbers of putative ARGs on the genome of Lelliottia amnigena YTB01 (high Se tolerance) were much more than those of Pseudomonas putida (low Se tolerance), which indicated that Se may have cooperative associations with bacterial ARGs. Altogether, these observations supported that Se influences the pattern of ARGs in Se-rich forest soil. However, compared with the environmental concentration, the concentration of sodium selenite used to select culturable bacteria may be too high; this constitutes a limitation of this study, and further studies are required to explore lower concentration of sodium selenite in this regard.

The genomic mining of both Lelliottia amnigena YTB01 and Pseudomonas putida indicated that some multidrug efflux RND transporter and macrolide ABC transporter harbored dual functions as transporter protein and resistance genes, regulated by RNA polymerase and DNA-binding response regulator, respectively. Altogether, these gene repertoire confirms the insights obtained from the metagenomic on the mechanism of Se detoxification in the sampled soils. Furthermore, it suggests possible co-selective or co-regulatory resistance mechanisms for transporting both Se and antibiotics to the extracellular space. In summary, high Se pressure may dominantly drive the co-selection of Se resistance and antibiotic resistance through cross-resistance and co-regulatory resistance), thereby influencing ARGs diversity and abundance.

Conclusion

In this study, by investigating soils with a wide range of Se contents, we provided evidence that the diversity and abundance of ARGs are enhanced under moderate Se content but inhibited under severe Se pressure. In addition, the methylation of Se (mediated by tehB) and the repair of DNA damages (mediated by ruvB and recG) are the dominant mechanisms of Se resistance in the study area, and such Se resistance may influence ARGs pattern. Bacterial communities also play an important role in mediating the dual effects of soil Se on the distribution of ARGs. Furthermore, the cultured high Se-tolerant bacteria possessed more ARGs on the genome than the low Se-tolerant bacteria. Remarkably, the world average soil Se content is lower than the threshold of 0.034 mg kg−1; therefore, elevated soil Se content may facilitate the spread of antibiotic resistance.

Availability of data and materials

The datasets generated during the current study are available in the NODE repository (https://www.biosino.org/node/index), including all 16S rRNA gene sequences (Accession No. OEP001738) and the complete genome sequences of two isolates (Accession No. OEP001741).

Abbreviations

ARGs:

Antibiotic resistance genes

BGCs:

Biosynthetic gene clusters

HGT:

Horizontal gene transfer

MGEs:

Mobile genetic elements

RDA:

Redundancy analysis

SEM:

Structural equation modeling

VGT:

Vertical gene transfer

ncRNA:

Non-coding RNA

MLSB:

Macrolide lincosamide streptogramin_b

OTUs:

Operational taxonomic units

SOC:

Soil organic carbon

TN:

Total nitrogen

TCN:

Tetracycline

DOX:

Doxycycline hydrochloride

OLF:

Ofloxacin

ENR:

Enrofloxacin

AMX:

Amoxicillin

CA:

Clavulanic acid

KAN:

Kanamycin

ENY:

Erythromycin

References

  1. Zowawi HM, Harris PN, Roberts MJ, Tambyah PA, Schembri MA, Pezzani MD, et al. The emerging threat of multidrug-resistant Gram-negative bacteria in urology. Nat Rev Urol. 2015;1210:570–84.

    Article  CAS  Google Scholar 

  2. Pamer EG. Resurrecting the intestinal microbiota to combat antibiotic-resistant pathogens. Science. 2016;3526285:535.

    Article  CAS  Google Scholar 

  3. Fresia P, Antelo V, Salazar C, Gimenez M, D’Alessandro B, Afshinnekoo E, et al. Urban metagenomics uncover antibiotic resistance reservoirs in coastal beach and sewage waters. Microbiome. 2019;7:35.

    Article  PubMed  PubMed Central  Google Scholar 

  4. Pruden A, Pei RT, Storteboom H, Carlson KH. Antibiotic resistance genes as emerging contaminants: studies in northern Colorado. Environ Sci Technol. 2006;4023:7445–50.

    Article  CAS  Google Scholar 

  5. Aminov RI. Horizontal gene exchange in environmental microbiota. Front Microbiol. 2011;2:158.

    Article  PubMed  PubMed Central  Google Scholar 

  6. Seiler C, Berendonk TU. Heavy metal driven co-selection of antibiotic resistance in soil and water bodies impacted by agriculture and aquaculture. Front Microbiol. 2012;369:399.

    Google Scholar 

  7. Zhu YG, Johnson TA, Su JQ, Qiao M, Guo GX, Stedtfeld RD, et al. Diverse and abundant antibiotic resistance genes in Chinese swine farms. Proc Natl Acad Sci USA. 2013;1109:3435–40.

    Article  Google Scholar 

  8. Berg J, Thorsen MK, Holm PE, Jensen J, Nybroe O, Brandt KK. Cu exposure under field conditions coselects for antibiotic resistance as determined by a novel cultivation-independent bacterial community tolerance assay. Environ Sci Technol. 2010;4422:8724–8.

    Article  CAS  Google Scholar 

  9. Sun W, Qian X, Gu J, Wang XJ, Zhang L, Guo AY. Mechanisms and effects of arsanilic acid on antibiotic resistance genes and microbial communities during pig manure digestion. Bioresour Technol. 2017;234:217–23.

    Article  CAS  PubMed  Google Scholar 

  10. Zhao X, Shen J-P, Zhang L-M, Du S, Hu H-W, He J-Z. Arsenic and cadmium as predominant factors shaping the distribution patterns of antibiotic resistance genes in polluted paddy soils. J Hazard Mater. 2020;389:121838.

    Article  CAS  PubMed  Google Scholar 

  11. Combs GF, Scott ML. Nutritional interrelationships of vitamin-E and selenium. Bioscience. 1977;277:467–73.

    Article  Google Scholar 

  12. Nancharaiah YV, Lens PN. Ecology and biotechnology of selenium-respiring bacteria. Microbiol Mol Biol Rev. 2015;791:61–80.

    Article  Google Scholar 

  13. Heider J, Bock A. Selenium metabolism in microorganisms. Adv Microb Physiol. 1993;35:71–109.

    Article  CAS  PubMed  Google Scholar 

  14. Bulteau AL, Chavatte L. Update on selenoprotein biosynthesis. Antioxid Redox Sign. 2015;2310:775–94.

    Article  CAS  Google Scholar 

  15. Letavayova L, Vlasakova D, Vlckova V, Brozmanova J, Chovanec M. Rad52 has a role in the repair of sodium selenite-induced DNA damage in Saccharomyces cerevisiae. Mutat Res-Gen Tox En. 2008;6522:198–203.

    Article  CAS  Google Scholar 

  16. Schiavon M, Ertani A, Parrasia S, Vecchia FD. Selenium accumulation and metabolism in algae. Aquat Toxicol. 2017;189:1–8.

    Article  CAS  PubMed  Google Scholar 

  17. Van Hoewyk D, Takahashi H, Inoue E, Hess A, Tamaoki M, Pilon-Smits EAH. Transcriptome analyses give insights into selenium-stress responses and selenium tolerance mechanisms in Arabidopsis. Physiol Plantarum. 2008;1322:236–53.

    Google Scholar 

  18. Bebien M, Lagniel G, Garin J, Touati D, Vermeglio A, Labarre J. Involvement of superoxide dismutases in the response of Escherichia coli to selenium oxides. J Bacteriol. 2002;1846:1556–64.

    Article  CAS  Google Scholar 

  19. Miranda AT, Gonzalez MV, Gonzalez G, Vargas E, Campos-Garcia J, Cervantes C. Involvement of DNA helicases in chromate resistance by Pseudomonas aeruginosa PAO1. Mutat Res-Fund Mol Mech. 2005;5781–2:202–9.

    Article  CAS  Google Scholar 

  20. Liu M, Turner RJ, Winstone TL, Saetre A, Brenzinger MD, Jickling G, Tari LW. Escherichia coli TehB requires S-adenosylmethionine as a cofactor to mediate tellurite resistance. J Bacteriol. 2000;182:6509–13.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  21. Li Z, Tan J, Shao L, Dong X, Ye RD, Chen D. Selenium-mediated protection in reversing the sensitivity of bacterium to the bactericidal antibiotics. J Trace Elem Med Biol. 2017;41:23–31.

    Article  CAS  PubMed  Google Scholar 

  22. Li L-G, Xia Y, Zhang T. Co-occurrence of antibiotic and metal resistance genes revealed in complete genome collection. ISME J. 2017;113:651–62.

    Article  CAS  Google Scholar 

  23. Kieliszek M. Selenium-fascinating microelement, properties and sources in food. Molecules. 2019;24:1298.

    Article  CAS  PubMed Central  Google Scholar 

  24. Rayman MP. Selenium and human health. Lancet. 2012;379:1256–68.

    Article  CAS  PubMed  Google Scholar 

  25. Zhang L, Song H, Guo Y, Fan B, Huang Y, Mao X, et al. Benefit-risk assessment of dietary selenium and its associated metals intake in China (2017–2019): is current selenium-rich agro-food safe enough? J Hazard Mater. 2020;398:123224.

    Article  CAS  PubMed  Google Scholar 

  26. Bermingham EN, Hesketh JE, Sinclair BR, Koolaard JP, Roy NC. Selenium-enriched foods are more effective at increasing glutathione peroxidase (GPx) activity compared with selenomethionine: a meta-analysis. Nutrients. 2014;610:4002–31.

    Article  CAS  Google Scholar 

  27. Phiri FP, Ander EL, Bailey EH, Chilima B, Chilimba ADC, Gondwe J, et al. The risk of selenium deficiency in Malawi is large and varies over multiple spatial scales. Sci Rep-UK. 2019;9:6566.

    Article  CAS  Google Scholar 

  28. Forsberg KJ, Reyes A, Wang B, Selleck EM, Sommer MOA, Dantas G. The shared antibiotic resistome of soil bacteria and human pathogens. Science. 2012;337:1107–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  29. Zhu J, Wang N, Li S, Li L, Su H, Liu C. Distribution and transport of selenium in Yutangba, China: impact of human activities. Sci Total Environ. 2008;3922–3:252–61.

    Article  CAS  Google Scholar 

  30. Yu T, Hou WL, Hou QY, Ma WJ, Xia XQ, Li YT, et al. Safe utilization and zoning on natural selenium-rich land resources: a case study of the typical area in Enshi County. China Environ Geochem Hlth. 2020;429:2803–18.

    Article  CAS  Google Scholar 

  31. García JB, Krachler M, Chen B, Shotyk W. Improved determination of selenium in plant and peat samples using hydride generation-atomic fluorescence spectrometry (HG-AFS). Anal Chim Acta. 2005;5342:255–61.

    Article  CAS  Google Scholar 

  32. Wang FH, Qiao M, Su JQ, Chen Z, Zhou X, Zhu YG. High throughput profiling of antibiotic resistance genes in urban park soils with reclaimed water irrigation. Environ Sci Technol. 2014;4816:9079–85.

    Article  CAS  Google Scholar 

  33. Yuan X, Qiang Z, Ben W, Zhu B, Qu J. Distribution, mass load and environmental impact of multiple-class pharmaceuticals in conventional and upgraded municipal wastewater treatment plants in East China. Environ Sci-Proc Imp. 2015;173:596–605.

    Google Scholar 

  34. Muziasari WI, Pitkanen LK, Sorum H, Stedtfeld RD, Tiedje JM, Virta M. The resistome of farmed fish feces contributes to the enrichment of antibiotic resistance genes in sediments below baltic sea fish farms. Front Microbiol. 2017;7:2137.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Schimel J, Balser TC, Wallenstein M. Microbial stress-response physiology and its implications for ecosystem function. Ecology. 2007;886:1386.

    Article  Google Scholar 

  36. Han XM, Hu HW, Shi XZ, Wang JT, Han LL, Chen D, et al. Impacts of reclaimed water irrigation on soil antibiotic resistome in urban parks of Victoria, Australia. Environ Pollut. 2016;211:48–57.

    Article  CAS  PubMed  Google Scholar 

  37. Suzuki MT, Taylor LT, DeLong EF. Quantitative analysis of small-subunit rRNA genes in mixed microbial populations via 5 ’-nuclease assays. Appl Environ Microb. 2000;6611:4605–14.

    Article  Google Scholar 

  38. He JZ, Shen JP, Zhang LM, Zhu YG, Zheng YM, Xu MG, et al. Quantitative analyses of the abundance and composition of ammonia-oxidizing bacteria and ammonia-oxidizing archaea of a Chinese upland red soil under long-term fertilization practices. Environ Microbiol. 2007;9:2364–74.

    Article  CAS  PubMed  Google Scholar 

  39. Ranjard L, Prigent-Combaret C, Favre-Bonte S, Monnez C, Nazaret S, Cournoyer B. Characterization of a novel selenium methyltransferase from freshwater bacteria showing strong similarities with the calicheamicin methyltransferase. Biochim Biophys Acta-Gene Struct Expr. 2004;16791:80–5.

    Article  CAS  Google Scholar 

  40. Driscoll DM, Copeland PR. Mechanism and regulation of selenoprotein synthesis. Annu Rev Nutr. 2003;23:17–40.

    Article  CAS  PubMed  Google Scholar 

  41. Krafft T, Bowen A, Theis F, Macy JM. Cloning and sequencing of the genes encoding the periplasmic-cytochrome B-containing selenate reductase of Thauera selenatis. DNA Seq J DNA Seq Mapp. 2000;106:365–77.

    Article  Google Scholar 

  42. Kuroda M, Yamashita M, Miwa E, Imao K, Fujimoto N, Ono H, et al. Molecular cloning and characterization of the srdBCA operon, encoding the respiratory selenate reductase complex, from the selenate-reducing Bacterium Bacillus selenatarsenatis SF-1. J Bacteriol. 2011;1939:2141–8.

    Article  CAS  Google Scholar 

  43. Zhou JZ, Wu LY, Deng Y, Zhi XY, Jiang YH, Tu QC, et al. Reproducibility and quantitation of amplicon sequencing-based detection. ISME J. 2011;58:1303–13.

    Article  CAS  Google Scholar 

  44. Caporaso JG, Kuczynski J, Stombaugh J, Bittinger K, Bushman FD, Costello EK, et al. QIIME allows analysis of high-throughput community sequencing data. Nat Methods. 2010;75:335–6.

    Article  CAS  Google Scholar 

  45. Washburne AD, Morton JT, Sanders J, McDonald D, Zhu QY, Oliverio AM, et al. Methods for phylogenetic analysis of microbiome data. Nat Microbiol. 2018;36:652–61.

    Article  CAS  Google Scholar 

  46. Edgar RC. Search and clustering orders of magnitude faster than BLAST. Bioinformatics. 2010;26(19):2460–1. https://doi.org/10.1093/bioinformatics/btq461.

    Article  CAS  PubMed  Google Scholar 

  47. Chin CS, Peluso P, Sedlazeck FJ, Nattestad M, Concepcion GT, Clum A. Phased diploid genome assembly with single-molecule real-time sequencing. Nat Methods. 2016;13(12):1050–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  48. Hackl T, Hedrich R, Schultz J, Förster F. Proovread: large-scale high-accuracy PacBio correction through iterative short read consensus. Bioinformatics. 2014;30(21):3004–11.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  49. Denisov G, Walenz B, Halpern AL, Miller J, Axelrod N, Levy S. Consensus generation andvariant detection by Celera Assembler. Bioinformatics. 2008;24(8):1035–40.

    Article  CAS  PubMed  Google Scholar 

  50. Delcher AL, Harmon D, Kasif S, White O, Salzberg SL. Improved microbial gene identification with GLIMMER. Nucleic Acids Res. 1999;2723:4636–41.

    Article  Google Scholar 

  51. Delcher AL, Bratke KA, Powers EC, Salzberg SL. Identifying bacterial genes and endosymbiont DNA with Glimmer. Bioinformatics. 2007;236:673–9.

    Article  Google Scholar 

  52. Bagheri H, Dyer R, Severin AJ, Rajan H. Comprehensive Analysis of Non Redundant Protein Database. Off Biotechnol Publ. 2020;9.

  53. Bairoch A, Apweiler R. The SWISS-PROT protein sequence data bank and its supplement TrEMBL in 1999. Nucleic Acids Res. 1999;27:49–54.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Ashburner M, Ball CA, Blake JA, Botstein D, Butler H, Cherry JM, et al. Gene ontology: tool for the unification of biology. Nat Genet. 2000;251:25–9.

    Article  CAS  Google Scholar 

  55. Cepas JH, Szklarczyk D, Heller D, Plaza AH, Forslund SK, Cook H, et al. eggNOG 5.0: a hierarchical, functionally and phylogenetically annotated orthology resource based on 5090 organisms and 2502 viruses. Nucleic Acids Res. 2019;47(D1):D309–14.

    Article  CAS  Google Scholar 

  56. Alcock BP, Raphenya AR, Lau TTY, Tsang KK, Bouchard M, Edalatmand A, et al. CARD 2020: antibiotic resistome surveillance with the comprehensive antibiotic resistance database. Nucleic Acids Res. 2020;48(D1):D517–25.

    CAS  PubMed  Google Scholar 

  57. Blin K, Shaw S, Steinke K, Villebro R, Ziemert N, Lee SY, et al. antiSMASH 5.0: updates to the secondary metabolite genome mining pipeline. Nucleic Acids Res. 2019;47:81–7.

    Article  CAS  Google Scholar 

  58. Fong Y, Huang Y, Gilbert PB, Permar SR. chngpt: threshold regression model estimation and inference. BMC Bioinform. 2017;18:454.

    Article  Google Scholar 

  59. Oksanen J, Blanchet FG, Kindt R, Legendre P, O'Hara RG, Simpson G, et al. vegan: community ecology package 2.5–7. Ordination methods, diversity analysis and other functions for community and vegetation ecologists. 2020.

  60. Rosseel Y. lavaan: an R package for structural equation modeling. J Stat Softw. 2012;48(2):1–36.

    Article  Google Scholar 

  61. Lefcheck JS. piecewiseSEM: piecewise structural equation modeling in R for ecology, evolution, and systematics. Methods Ecol Evol. 2015;7(5):573–9.

    Article  Google Scholar 

  62. Revelle W. psych: procedures for personality and psychological research. Evanston: Northwestern University; 2018.

    Google Scholar 

  63. Bastian M, Heymann S, Jacomy M. Gephi: an open source software for exploring and manipulating networks. In: International AAAI conference on weblogs and social media. 2009.

  64. Vaz-Moreira I, Nunes OC, Manaia CM. Bacterial diversity and antibiotic resistance in water habitats: searching the links with the human microbiome. Fems Microbiol Rev. 2014;384:761–78.

    Article  CAS  Google Scholar 

  65. Li Y, Cao W, Liang S, Yamasaki S, Chen X, Shi L, et al. Metagenomic characterization of bacterial community and antibiotic resistance genes in representative ready-to-eat food in southern China. Sci Rep-UK. 2020;101:15175–15175.

    Article  CAS  Google Scholar 

  66. Baker-Austin C, Wright MS, Stepanauskas R, McArthur JV. Co-selection of antibiotic and metal resistance. Trends Microbiol. 2006;144:176–82.

    Article  CAS  Google Scholar 

  67. Gullberg E, Albrecht LM, Karlsson C, Sandegren L, Andersson DI. Selection of a multidrug resistance plasmid by sublethal levels of antibiotics and heavy metals. MBio. 2014;5:e01918-e2014.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  68. Seiler C, Berendonk TU. Heavy metal driven co-selection of antibiotic resistance in soil and water bodies impacted by agriculture and aquaculture. Front Microbiol. 2012;3:399.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Hu H-W, Wang J-T, Li J, Shi X-Z, Ma Y-B, Chen D, et al. Long-term nickel contamination increases the occurrence of antibiotic resistance genes in agricultural soils. Environ Sci Technol. 2017;512:790–800.

    Article  CAS  Google Scholar 

  70. Forsberg KJ, Patel S, Gibson MK, Lauber CL, Knight R, Fierer N, et al. Bacterial phylogeny structures soil resistomes across habitats. Nature. 2014;509:612-+.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Stolz JE, Basu P, Santini JM, Oremland RS. Arsenic and selenium in microbial metabolism. Annu Rev Microbiol. 2006;60:107–30.

    Article  CAS  PubMed  Google Scholar 

  72. Gillings MR, Gaze WH, Pruden A, Smalla K, Tiedje JM, Zhu YG. Using the class 1 integron-integrase gene as a proxy for anthropogenic pollution. ISME J. 2015;96:1269–79.

    Article  CAS  Google Scholar 

  73. Hasman H, Aarestrup FM. tcrB, a gene conferring transferable copper resistance in Enterococcus faecium: occurrence, transferability, and linkage to macrolide and glycopeptide resistance. Antimicrob Agents Chemother. 2002;465:1410–6.

    Article  CAS  Google Scholar 

  74. Hasman H, Aarestrup FM. Relationship between copper, glycopeptide, and macrolide resistance among Enterococcus faecium strains isolated from pigs in Denmark between 1997 and 2003. Antimicrob Agents Chemother. 2005;491:454–6.

    Article  CAS  Google Scholar 

  75. Stolz JF, Basu P, Oremland RS. Microbial transformation of elements: the case of arsenic and selenium. Int Microbiol. 2002;54:201.

    Google Scholar 

  76. Rosenfeld CE, Kenyon JA, James BR, Santelli CM. Selenium (IV, VI) reduction and tolerance by fungi in an oxic environment. Geobiology. 2017;15:441–52.

    Article  CAS  PubMed  Google Scholar 

  77. Dhanjal S, Singh AK, Cameotra SS. Global gene expression analysis of bacterial stress response to elevated concentrations of toxic metalloids—selenium and arsenic. Geomicrobiol J. 2014;316:480–92.

    Article  CAS  Google Scholar 

Download references

Acknowledgements

We gratefully acknowledge Dr Hongzhao Yuan for her technical support of the instrumental operation, and express our gratitude to EditSprings (https://www.edits-prings.com/) for the expert linguistic services provided.

Funding

This study was financially supported by the Cooperative Innovation Center for Sustainable Pig Production of Huazhong Agricultural University and National Natural Science Foundation of China (Grant No. 41771336).

Author information

Author notes

  1. Fang-Fang Wang, Guo-Ping Liu and Fan Zhang contributed equally to this work.

Authors and Affiliations

  1. College of Animal Science, Yangtze University, Jingzhou, 434025, China

    Fang-Fang Wang, Guo-Ping Liu, Zong-Ming Li & Xiao-Lin Yang

  2. State Key Laboratory of Urban and Regional Ecology, Research Center for Eco-Environmental Sciences, Chinese Academy of Sciences, Beijing, 100085, China

    Fang-Fang Wang

  3. University of the Chinese Academy of Sciences, Beijing, 100049, China

    Fang-Fang Wang & Zong-Ming Li

  4. College of Horticulture and Gardening, Yangtze University, Jingzhou, 434025, Hubei, China

    Fan Zhang & Chao-Dong Yang

  5. Key Laboratory of Agro-Ecological Processes in the Subtropical Region, Institute of Subtropical Agriculture, Chinese Academy of Sciences, Changsha, 410125, China

    Zong-Ming Li & Jian-Lin Shen

  6. Faculty of Veterinary and Agricultural Sciences, The University of Melbourne, Parkville, VIC, 3010, Australia

    Ji-Zheng He

  7. Department of Botany and Plant Sciences, University of California, Riverside, CA, 92521-0124, USA

    B. Larry Li

  8. College of Veterinary Medicine, Hunan Agricultural University, Changsha, 410128, China

    Guo-Ping Liu & Jian-Guo Zeng

Authors
  1. Fang-Fang Wang
    View author publications

    You can also search for this author in PubMed Google Scholar

  2. Guo-Ping Liu
    View author publications

    You can also search for this author in PubMed Google Scholar

  3. Fan Zhang
    View author publications

    You can also search for this author in PubMed Google Scholar

  4. Zong-Ming Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  5. Xiao-Lin Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  6. Chao-Dong Yang
    View author publications

    You can also search for this author in PubMed Google Scholar

  7. Jian-Lin Shen
    View author publications

    You can also search for this author in PubMed Google Scholar

  8. Ji-Zheng He
    View author publications

    You can also search for this author in PubMed Google Scholar

  9. B. Larry Li
    View author publications

    You can also search for this author in PubMed Google Scholar

  10. Jian-Guo Zeng
    View author publications

    You can also search for this author in PubMed Google Scholar

Contributions

GPL, JZH, BLL and FFW designed the research, FFW, FZ, XLY, CDY and JLS participate in field work, FFW and FZ lead the lab work; FFW, GPL, ZML performed data analysis; FFW, GPL, XLY, CDY, JZH, BLL, JLS and ZJG wrote the manuscript in close consultation from all authors.

Corresponding author

Correspondence to Guo-Ping Liu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher's Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig. S1

. Heatmap of the relative abundance of ARGs. The column is labeled with the content of available Se content of each sampling from low to high. The row is the ARGs detected in each sample and grouped according to antibiotic resistance type.

Additional file 2: Fig. S2

. Relative abundance of the bacterial community at the phylum level.

Additional file 3: Fig. S3

. Co-occurrence network analysis showing the correlation between ARGs, MGEs, and Se resistance genes (R > 0.8, P < 0.05). Node size indicates connectivity, and the color of nodes is grouped by modularity class.

Additional file 4: Fig. S4

. The genome circles graphs of Lelliottia amnigena YTB01 (a) and Pseudomonas putida (b). From outer to inner: 1: Genome Size; 2: Forward Strand Gene, colored according to the cluster of orthologous groups (COG) classification; 3: Reverse Strand Gene, colored according to the cluster of orthologous groups (COG) classification; 4: Forward Strand ncRNA; 5: Reverse Strand ncRNA; 6: repeat; 7: GC; 8: GC-SKEW.

Additional file 5: Material 1

. The genome annotation of various putative proteins involved in Se resistance and antibiotic resistance of Lelliottia amnigena YTB01.

Additional file 6: Material 2

. The genome annotation of various putative proteins involved in Se resistance and antibiotic resistance of Pseudomonas putida.

Additional file 7: Table S1

. Sampling site locations and soil chemical properties. Table S2. Primer sequences and quantitative PCR reaction conditions of 16 selenium resistance genes. Table S3. Pearson correlation between soil chemical properties. “*” significant correlations at P < 0.05; “**” Significant correlations at P < 0.001. Table S4. The detected concentration of antibiotic residues (µg kg−1). Table S5. The absolute copy numbers of selenium resistance genes (copies g−1).

Additional file 8: Text S1

. Determination of total and available Se contents of soil.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Wang, FF., Liu, GP., Zhang, F. et al. Natural selenium stress influences the changes of antibiotic resistome in seleniferous forest soils. Environmental Microbiome 17, 26 (2022). https://doi.org/10.1186/s40793-022-00419-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s40793-022-00419-z

Keywords

Environmental Microbiome

ISSN: 2524-6372