Heavy Metal-Induced Co-Selection for Antibiotic Resistance in Terrestrial Subsurface Soils


 Background: Terrestrial surface ecosystems are important sinks for antibiotic resistance genes (ARGs) due to the continuous discharge of contaminants from human-impacted ecosystems. However, factors determining the abundance and resistance types of ARGs in terrestrial subsurface soils remain largely unknown. In this study, we investigated the abundance and diversity of ARGs, and their correlations with metal resistance genes (MRGs), mobile genetic elements (MGEs), bacteria, and heavy metals in subsurface soils in a global scale using high throughput quantitative PCR and metagenomic sequencing approaches. Results: Abundant and diverse ARGs were detected with high spatial heterogeneity among the sampling sites. Vertically, there was no significant difference in the ARG profiles between the aquifer and non-aquifer soils. Heavy metals were the key factors shaping ARG patterns in soils with high heavy metal contents, while they induced no significant effect in low contents. Moreover, heavy metals could trigger the proliferation of antibiotic resistance by increasing MGE abundance or influencing bacterial communities. Metagenomic analysis also revealed the widespread co-occurrence of ARGs and MRGs, with heavy metals possibly aggravating the co-selection of ARGs and MRGs in soils with high heavy metal contents. Conclusions: This study highlighted the heavy metal-induced co-selection for ARGs and MRGs and revealed the occurrence of ARG pollution in terrestrial subsurface soils.


Introduction
The increasing spread and accumulation of antibiotic resistance genes (ARGs) in the environment due to unregulated use pose risks to human/animal health and challenge life-saving antibiotic therapies, which have raised global concerns [1][2][3][4]. Antibiotic resistance is prevalent in many natural environments including soil [5][6], surface water [7], groundwater [8], air [9], even in the gut microbiota [10][11]. The dispersal of ARGs in these environments presently is mainly associated with the improper disposal of environmental wastes containing contaminants as products of anthropogenic activities [5,[12][13].
Groundwater is a vital public resource used in a myriad of things including agricultural irrigation, industry, and drinking water. However, the groundwater has also become an important gene pool for ARGs recently [7,[14][15]. Many studies have revealed the prevalence of groundwater ARGs under contaminated locations, such as livestock facilities [16] and waste land lls [8]. Terrestrial surface water is another important sink for the ARGs [7,15]. Mutual replenishment process of surface water and groundwater may accelerate ARG spread in both surface and subsurface ecosystems. In addition, anthropogenic contaminants in surface waters can seep into the deep soil and groundwater [17], continuously exerting selection pressures on ARGs in subsurface systems. Therefore, it can be hypothesized that there may be no signi cant difference in ARG pro les between the subsurface aquifer and non-aquifer soils.
In the environment, determinants of antibiotic resistance triggering the occurrence and proliferation of ARGs are complex [12,[18][19]. Among these, heavy metal contamination is one of the most important factors in uencing ARGs [5,12,20]. Heavy metals can co-select for ARGs and metal resistance genes (MRGs) using several mechanisms, such as co-resistance [21]. The potential risks correlated with ARGs and MRGs have been emphasized in both environmental and medical implications, with their cooccurrence especially being more prevalent in human pathogens [22]. Moreover, it is of particular concern that heavy metals can serve as strong and long-term selective pressures for ARGs and MRGs [23][24].
The environmental conditions, such as pH, oxygen level, and temperature [25][26][27], can also affect heavy metal toxicity by impacting metal ions valence state and thus their bioavailability [28][29]. The toxicity of heavy metals also depends largely on its concentration. Toxic heavy metals can further urge microbial communities to develop the resistance [28,[30][31][32]. However, the environmental situation in subsurface soils is relatively anoxic, which might impact heavy metal toxicity and ARG attenuation [33][34][35]. Presently, many studies have been made to explore antibiotic resistance in a variety of environments but most mainly focused on surface ecosystems [13,[36][37]. Knowledge about the distribution and resistance determiners of ARGs in subsurface vertical soils on the other hand remain limited. This study then aims to i) assess the vertical occurrence and distribution of ARGs in the terrestrial subsurface aquifer and nonaquifer soils, and ii) explore the associations between ARG prevalence and environmental parameters, particularly with heavy metals.
The ARGs and MGEs in the soil samples showed strong spatial heterogeneity with abundances ranging from 1.9 × 10 4 to 8.1 × 10 6 copies/g, 428 to 4.7 × 10 6 copies/g, respectively. Among the samples, Zurich had the highest absolute abundance (3.7 ± 2.9 × 10 5 copies/g), relative abundance (0.1 ± 0.01 copies/16S rRNA gene copy) and number (up to 151 subtypes) of ARGs. Bremen in contrast had the least number of subtypes with only 28 subtypes of ARGs and also lowest absolute (1.9 ± 0. 9 × 10 5 copies/g) and relative (0.02 ± 0.008 copies/16S rRNA gene copy) abundance ( Fig. 1c; Additional le 1: Figure S3). There were signi cant correlations between ARGs and MGEs (Additional le 1: Figures S4 and S5, Table   S1), indicating potential horizontal ARG transfer via MEGs. Vertically, the ARG and MGE pro les in the aquifer soils were similar with those in the non-aquifer soils (Additional le 1: Figure S6).

Co-occurrence patterns between ARGs and bacterial taxa
High throughput sequencing of the bacterial 16S rRNA gene revealed that the with most prevalent phyla among the samples was Chloro exi (16.5%), followed by Proteobacteria (16.4%), and Actinobacteria (15.8%) (Fig. 2). NMDS plots further showed a similar pattern in bacterial community composition between the aquifer and non-aquifer soils ( Fig. 2a; Additional le 1: Figure S7). However, soils with high heavy metal contents had bacterial communities that were further separated by heavy metal types (Fig.   2c). A close relationship among bacteria and ARGs was revealed by the strong (r > 0.8) and signi cant (P < 0.05) co-occurrence patterns in the network (Additional le 1: Figure S8). Proteobacteria, Chloro exi, and Actinobacteria were connected with diverse ARGs mainly demonstrated potential aminoglycoside, multidrug, and β-lactamase resistance. Procrustes analysis also indicated the highly signi cant associations between ARG pro les and the bacterial communities (M 2 = 0.8, r = 0.4, P = 0.002, 999 permutations), as con rmed by Mantel test (Spearman's r = 0.3, P = 0.005) and Spearman's correlation analysis (p < 0.05).

Factors in uencing ARGs
Variation partitioning analysis was performed to assess the lone and interactive effects of heavy metals, bacteria communities, MGEs, and physico-chemical factors on ARG pro les (Fig. 3). In soils with low heavy metal concentrations, bacterial communities had the highest variances (40.0%) for ARG composition, which were signi cantly higher than that of heavy metals (19.6%), MGEs (20.1%), and physico-chemical factors (0.4%) (Fig. 3b). However, in soils with high heavy metal concentrations, the main in uencing factors shifted to the heavy metals with the most explanation of 28.6%. The signi cant associations between heavy metals and ARG pro les were further con rmed by Procrustes analysis (M 2 = 0.805, r = 0.441, P = 0.02, 999 permutations). Higher heavy metals concentration signi cantly in uenced ARG pro les (Additional le 1: Table S2). Speci cally, As, Co, and Ni demonstrated the strongest correlations with ARGs con rmed by RDA and spearman's correlation ( Fig. 3a; Additional le 1: Figure  S9). In summary, heavy metals seemed to have strongly in uenced the proliferation of ARGs in the subsurface soils.

Effects of heavy metals on the overall patterns of ARGs
To determine the effects of heavy metals on ARGs, path analysis was conducted based on structural equation model (SEM) (Fig. 4). In soils with low heavy metal contents, heavy metal toxicity showed a low in uence on ARG abundance and diversity. Furthermore, spearman's correlation analysis showed no signi cant relationship between heavy metals and ARGs (Additional le 1: Table S3). However, in soils with high heavy metal contents, path analysis showed that soil heavy metal toxicity had a signi cant effect on ARG abundance and diversity (Fig. 4). Moreover, heavy metal toxicity obviously changed the bacterial community composition (λ = 0.595, p < 0.001) and MGEs (λ = 0.276, p < 0.05), subsequently positively in uencing the ARGs. This con rmed the positive relationship between heavy metals and ARGs via bacteria or MGEs. Correlation analysis further showed signi cant positive associations among metal toxicity, bacterial community composition, MGEs, and the diversity and abundance ARGs (Additional le 1: Tables S2, S4, and S5).
Signi cant and positive co-occurrence patterns (r = 0.6, p < 0.05) between heavy metals and ARGs were assessed in soils with high heavy metal contents (Additional le 1: Figure S10). Heavy metal As, Pb, Co, Cr, and Ni all correlated with ARGs conferring resistance to different antibiotics including multidrug, aminoglycosides, tetracycline and vancomycin (Additional le 1: Table S6). The correlations between heavy metal contaminants and ARGs suggested a signi cant degree of heavy metal selection via linkage with other factors.

Metagenomic insights into the correlations between ARGs and MRGs
Since heavy metals were the main factors in uencing ARGs, the co-occurrence of ARGs and MRGs in soils with high heavy metal content should be expected. Shotgun metagenomic analysis was used to con rm the co-selection of ARGs and MRGs induced by heavy metals (Fig. 5). A total of 571 ARGs and 302 MRGs were detected with relatively high associations (R 2 > 0.9). Consistent with HT-qPCR results, the most abundant ARGs were the multidrug resistance genes. The detected MRGs contained the resistance to heavy metals such as As, Cu, Zn, Pb, of which, the multi-metal resistance genes were the most abundant (1.1 ± 0.02×10 -2 ), followed by Cu (0.5 ± 0.02×10 -2 ) and Zn (0.3 ± 0.03×10 -2 ) resistance genes. In the MRG-metal network, Cu, Zn, Pb, Cr closely connected with genes confer resistance to different heavy metals (Additional le 1: Figure S11), indicating the strong in uence from heavy metals on MRGs.
Based on ARG and MRG co-occurrence pro les, multidrug resistance genes (n = 82) showed the most associations with MRGs, followed by resistance genes belonging to the macrolide antibiotic (n = 52) and aminoglycoside (n = 45). The abundance and number of genes conferred with multi-resistance to both antibiotic and metals were signi cantly higher than that of single resistance (Fig. 5a). Most of the ARG-MRG subtypes were detected in both soils with high and low heavy metal contents. However, the unique genes were more diverse and abundant (p = 0.048, Fig. 5c) and the signatures of ARG and MRG cooccurrence were more frequent in soils with higher heavy metal content. ARG abundance in soils increased with high heavy metal contents, much greater than in low heavy metal contents (Fig. 5b). Moreover, heavy metals Cu (p = 0.03, r = 0.929), Pb (p < 0.001, r = 1), and Zn (p = 0.014, r = 0.857) were signi cantly positively correlated with gene abundance that has resistance to both antibiotic and heavy metals, as con rmed by RDA analysis (Additional le 1: Figure S12).

Discussion
This study provided good insights into the heavy metal-driven co-selection of ARGs and MRGs in terrestrial subsurface soils. The ARGs and MRGs occurred abundantly and diversely in subsurface soils. Heavy metals could aggravate ARG pollution resulting in the high abundance of ARG-MRG co-occurrence in soils. Moreover, they were the dominant factors shaping ARG patterns competed over bacteria, MGEs, and physico-chemical parameters. These results highlighted the key role of heavy metals in disseminating ARGs in terrestrial subsurface soils.

ARG pollution in terrestrial subsurface soils
The abundance and number of ARGs were up to 10 6 copies/g and 151 respectively, higher than those reported in soils receiving swine and dairy manures [38] and sewage sludge [39]. Our results also revealed the ARG pollution in subsurface deep soils. Some ARG subtypes, such as sul1, sul2, and ermF, can easily permeate into deep soils via rainwater [40]. Beyond this, the environmental conditions of deep soils are advantageous anymore to ARG attenuation [33][34][35], leading to the ARG accumulation in the terrestrial subsurface soils. Furthermore, no apparent differences in ARG pro les were observed between aquifer and non-aquifer soils, consistent with our hypothesis. Mutual replenishment process of surface water and groundwater could have resulted in the mixing of ARGs between the aquifer and non-aquifer soils.

Heavy metals induced co-selection for ARGs
Metals at low concentrations are important for the growth of microorganisms [28][29], but high concentrations of which could be toxic to the environment [28,41]. The heavy metals signi cantly induced positive effects on ARG abundance and diversity, which were only detected when they occurred in high contents. Such effects largely depended on heavy metal concentrations as a result of the positively close linkages between resistance gene abundance and metal concentrations [42][43][44]. Moreover, heavy metals can trigger and catalyze the co-selection of antibiotic resistance [42,[45][46]. They can upregulate ARG subtypes by decreasing microbial susceptibility to antibiotics, ultimately resulting in the enhancement of ARGs [21,47].
Correlations in the metagenome between ARGs and MRGs further provided insights on the heavy metal induced co-selection of ARGs. In this study, diverse and abundant co-occurrence of ARGs and MRGs were observe, such as the aminoglycoside-Cu resistance. The co-occurrence of ARGs and MRGs in soils with high heavy metal levels were signi cantly higher than those with low contents. The genetic link between ARGs and MRGs has been con rmed [22]. The increase in MRGs induced by heavy metal pressure could facilitate the proliferation of ARGs. Moreover, coupling of genes for multidrug and multi-metals resistance were abundant, suggesting that many more types of antibiotic and metal resistance could be subjected to heavy metal selection pressure [5,12,20,23].
A signi cant effect of heavy metal was observed for both abundance and diversity of ARGs when induced by changes in MGE abundance. The higher abundance of MGEs and their strong positive associations with heavy metals were detected in soils with high heavy metal contents. The increased abundance of MGEs may enhance horizontal gene transfer of ARGs [13,[48][49]. Also, heavy metals tend to favor the antibiotic resistance by promoting the spread of MGEs via co-selection mechanism [5,21,[50][51]. Moreover, heavy metals could signi cantly in uence ARGs by shaping the bacterial communities. It is vital for bacteria to develop the resistant systems under selection pressures of heavy metal contaminations [30][31][32]. The toxic heavy metals drove the shift in the abundance and composition of some bacteria, which could conversely impact antibiotic resistance.
Previous studies reported that metal pollution is a key selector of ARGs [44,46,52], which was also con rmed in this study. The dominant factors contributing to the variation of ARG pro les shifted from bacteria to heavy metal contaminations. Such observations could be due to the signi cant heavy metal co-selection via MGEs, resulting in higher horizontal genes transfer and swift acquisition of ARGs in soils with high heavy metal contents. The remarkable responses to heavy metal contamination of ARGs found in this study has caused an alarm on the threats and risks for possible ARG spread in terrestrial subsurface ecosystems.

Conclusions
This study comprehensively investigated the occurrence of terrestrial subsurface ARGs and highlighted the considerable impact of heavy metals in the emergence and proliferation of antibiotic resistance combined with bacteria and MGEs. Once the subsurface soils are polluted by heavy metals, heavy metalinduced co-selection effects greatly accelerated the occurrence and propagation of ARGs and MRGs. This co-selection mechanism can cause compound pollution and pose risks in terrestrial subsurface environments. Thus, it is critical and necessary to control heavy metal contaminations in order to reduce the continuous spread of ARGs in terrestrial subsurface soils.  (Fig. 1a). These sampling sites covered different climatic conditions, soil types, geological zones and were demonstrated to be affected by distinct environmental factors (Additional le 1: Table S7). Three parallel quadrats were dug to collect subsamples at each sampling site. For each quadrat, ve parallel sample pits were drilled. The soil samples collected from the same depth in all ve pits were thoroughly mixed to generate a composite sample. All three replicate samples were analyzed separately and averaged to represent site conditions. A portion of the sample was separated and used for chemical analyses, another subsample was stored at -20 °C for further analysis.

Physico-chemical analyses
Soil ammonium, nitrite, and nitrate were extracted from 6 g soil added with 30 ml KCl solution (2 mol L -1 ) and measured by ow analyzers (Germany, SEAL, AA3). Soil moisture content were detected by oven-drying 10 g fresh soil at 105 °C until they reached constant weight. The pH value was determined by utilizing a pH Analyzer (Mettler Toledo, USA) in a supernatant after preparing the soil and water mixture (1:2.5). Total organic matter was determined using air-dried soil based on the LOI550 method [53]. Soil total nitrogen, total carbon and total sulfur were determined using an elemental analyzer system (Vario EL III, Elementar). Soil metal concentrations were measured by X-ray uorescence (XRF). Each parameter measured was determined in triplicate samples.
To provide a normalized contamination level of metals, metal toxicity index (TI) was expressed and calculated for each sample (Fig. 1b). TI was calculated using the formula in previous studies [54][55]: Here is the content of heavy metal i in the sample (mg/kg), is the half-maximal effective concentration for heavy metal i referred to the study of Welp [54]. Details about the heavy metal concentrations and types of each sample were in the Supplementary Information. DNA extraction and High-throughput quantitative PCR About 0.25 g freeze-dried and sieved soil was used to extract soil DNA using the PowerSoil kit (MO BIO, Carlsbad, CA, U.S.A.) referring to the instructions. The extracted genomic DNA were checked by 1% agarose gel and its concentration was determined using Nano Drop ND-2000 UV spectrophotometry (Thermo Fisher Scienti c Inc., USA).
To better reveal the antibiotic resistance in terrestrial subsurface ecosystems, high-throughput quantitative PCR was performed targeting the 283 ARGs, 12 MGEs, and one 16S rRNA gene, utilizing the WaferGen Smart Chip Real-time PCR system [13] (Additional le 1: Table S8). The 283 primer sets targeted ARGs included the major antibiotic resistance classi cations. The major integrons and transposases, such as intI-1, cintI-1, tnpA, were also included. All primer sets included no-template negative controls. Each quantitative PCR had technical triplicates and the e ciencies were all in the range of 1.7-2.3 and an R 2 above 0.98.
The standard curve method of quanti cation was used to determine the absolute 16S rRNA copy numbers in Roche 480 system. For standard calculation, a plasmid standard containing a cloned and sequenced 16S rRNA gene fragment was used to generate calibration curves. Each qPCR was carried out in triplicates with negative controls. The relative copy numbers of ARGs and MGEs generated by the HT-qPCR were transformed into absolute copy numbers by normalization, using the absolute 16S rRNA gene copy number.

Bacterial 16S rRNA gene and metagenomic shotgun sequencing and analysis
The V3-V4 region of the 16S rRNA gene was ampli ed with 338F and 806R prime set in an Illumina Miseq platform (Majorbio, China) to characterize bacterial communities. Paired end reads were demultiplexed, quality-ltered using QIIME [56] and were grouped into Operational Taxonomic Units (OTUs) using a 97% identity threshold with UPARSE. The taxonomic classi cation of each OTU was analyzed using RDP Classi er (http://rdp.cme.msu.edu/) against the Silva (SSU115) with con dence threshold of 70% [57].
Metagenomic shotgun library was constructed using TruSeqTM DNA Sample Prep Kit (Illumina, San Diego, CA, USA) and sequenced using an Illumina HiSeq4000 platform (Illumina Inc., San Diego, CA, USA). To obtain the high-quality sequences, the adapter sequences and low-quality reads were ltered from raw sequencing reads using fastp. After quality control, the clean sequences were assembled using Megahit [58]. The open reading frames (ORFs) of assembled contigs were predicted using MetaGene [59] with nucleic acid length ≥ 100 bp selected. The non-redundant gene sequences were constructed using CD-HIT [60].

Identi cation of ARGs and MRGs
The non-redundant gene sequences constructed from metagenomic data were used to identify ARGs and MRGs. The CARD [61] and BacMet [62] databases were respectively used for the annotations of ARGs and MRGs using BLASTP search [63] with an e-value cutoff ≤ 10 -5 . The relative abundance of each gene was determined using SOAPaligner [64] with 95% identity. The relative abundance of gene i was calculated by the following formula [65]: Here represents the relative abundance of gene i in a sample. represents the reads of gene i in a sample.
represents the total reads of all the genes in a sample.

Statistical analysis
Spearman's rank correlation and one-way ANOVA analyses were done in SPSS 10.0. Non-metric multidimensional scaling (NMDS) and variation partitioning were performed using R 3.6.2. Procrustes and Mantel tests were conducted in the R environment using the 'vegan' package [13]. Redundancy analysis (RDA) was performed based on forward selection in CANOCO 5.0. Network maps based on Spearman's rank correlations were built with the psych package in R 3.6.2. Only a correlation between two items with the r > 0.8 and the p-value < 0.01 was considered statistically robust and used to construct network in Gephi [66]. Correlation heatmaps were generated using vegan, ggcor, and dplyr packages in R 3.6.2. Linear regression was analysis and generated by Origin 9.0. A p-value < 0.05 indicated statistical signi cance. Circular stacked barplot were generated in R 3.6.2 using tidyverse and virdis packages. Upset plot were analyzed and generated in R 3.6.2 using UpSetR and yyplot packeages [67]. All bar graphs, box plots, and pie charts were generated by Origin.
Path analysis based on structural equation model (SEM) was used to assess the direct and indirect effects of soil metal contamination on ARG pro les by SPSS AMOS [5,45]. Theoretical assumptions were made to establish the model, such that (i) heavy metal could directly in uence the MGE abundance, bacterial diversity, ARG abundance, and ARG diversity; (ii) heavy metal could indirectly in uence ARG abundance and diversity via MGEs and bacteria; (iii) bacteria could indirectly in uence ARG abundance and diversity via MGEs. The ARG abundance and diversity, MGE abundance and bacterial diversity were used to t the model based on maximum-likelihood estimation method. The model t and its overall goodness-of-t were indicated by several parameters, namely (i) non-signi cant Chi square value (p > 0.05); (ii) low RMSEA value (< 0.08); (iii) high values of goodness-of-t indexes such as CFI, RFI, IFI, and TLI (> 0.9).

Declarations
Ethics approval and consent to participate Not applicable.

Consent for publication
Not applicable.

Availability of data and materials
All the sequence data from 16S rRNA gene amplicon and metagenome shotgun sequencing were submitted to NCBI SRA under BioProject accession number PRJNA421677.

Competing interests
The authors declare that they have no competing interests.  The global distribution of antibiotic resistance genes (ARGs) in subsurface vertical soils. a) Sampling sites of the vertical soils. Names with red and blue color respectively represent sampling sites with high and low heavy metal concentrations (The detailed information about the concentrations and types of heavy metals was in Supplementary Information). b) Box plots for heavy metals of each sampling site.
Whisker shows range from the 5th to 95th percentile. The inset box plots of As, Co, and Ni represent      Metagenomic insights into correlations between ARGs and MRGs. a) Circular stacked barplot showing the numbers of genes conferring resistance to both antibiotic and metals in microbiome. The top six of the microbial communities were displayed. "Proteo, Acido, Actino, Gemma, Chloro, Nitro" were the abbreviation of "Proteobacteria, Acidobacteria, Actinobacteria, Gemmatimonadetes, Chloro exi, Nitrospirae". A list of ARG-MRG were shown accounting for over 60%. b) ARG and MRG co-occurrence pro les in high and low heavy metal concentration samples. The inset depicts the MRG numbers encountered within 20. c) Upset plot showing the overlap among the samples. The transverse bar chart