Skip to main content
Download PDF

Growth of microbes in competitive lifestyles promotes increased ARGs in soil microbiota: insights based on genetic traits

Microbiome volume 13, Article number: 8 (2025) Cite this article

Abstract

Background

The widespread selective pressure of antibiotics in the environment has led to the propagation of antibiotic resistance genes (ARGs). However, the mechanisms by which microbes balance population growth with the enrichment of ARGs remain poorly understood. To address this, we employed microcosm cultivation at different antibiotic (i.e., Oxytetracycline, OTC) stresses across the concentrations from the environmental to the clinical. Paired with shot-gun metagenomics analysis and quantification of bacterial growth, trait-based assessment of soil microbiota was applied to reveal the association between key ARG subtypes, representative bacterial taxa, and functional-gene features that drive the growth of ARGs.

Results

Our results illuminate that resistome variation is closely associated with bacterial growth. A non-monotonic change in ARG abundance and richness was observed over a concentration gradient from none to 10 mg/l. Soil microbiota exposed to intermediate OTC concentrations (i.e., 0.1 and 0.5 mg/l) showed greater increases in the total abundance of ARGs. Community compositionally, the growth of representative taxa, i.e., Pseudomonadaceae was considered to boost the increase of ARGs. It has chromosomally carried kinds of multidrug resistance genes such as mexAB-oprM and mexCD-oprJ could mediate the intrinsic resistance to OTC. Streptomycetaceae has shown a better adaptive ability than other microbes at the clinical OTC concentrations. However, it contributed less to the ARGs growth as it represents a stress-tolerant lifestyle that grows slowly and carries fewer ARGs. In terms of community genetic features, the community aggregated traits analysis further indicates the enhancement in traits of resource acquisition and growth yield is driving the increase of ARGs abundance. Moreover, optimizations in energy production and conversion, alongside a streamlining of bypass metabolic pathways, further boost the growth of ARGs in sub-inhibitory antibiotic conditions.

Conclusion

The results of this study suggest that microbes with competitive lifestyles are selected under the stress of environmental sub-inhibitory concentrations of antibiotics and nutrient scarcity. They possess greater substrate utilization capacity and carry more ARGs, due to this they were faster growing and leading to a greater increase in the abundance of ARGs. This study has expanded the application of trait-based assessments in understanding the ecology of ARGs propagation. And the finding illustrated changes in soil resistome are accompanied by the lifestyle switching of the microbiome, which theoretically supports the ARGs control approach based on the principle of species competitive exclusion.

Video Abstract

Introduction

Antibiotic resistance is one of the greatest threats to public health today [1], and the environment is thought to play a critical role in its development and spread [2]. Investigating how microbes balance population growth with the per-microbe enrichment of antibiotic resistance genes (ARGs) is vital, as nutrient scarcity and sub-inhibitory antibiotic stress are the main constraints that influence the increased abundance of ARGs [3, 4].

Microbes navigate trade-offs between reproduction, survival, and competition in conditions with resource limitations and survival stress [5, 6]. These trade-offs could be framed within the trait-based life history strategy (LHS) framework, which was originally developed to elucidate the mechanisms by which organisms adapt to specific environments through trait selection [7, 8]. Advancements in trait-based research methods have currently progressed by microbial community ecologists, as the community aggregate trait (CAT) can be directly captured through high-throughput methods based on genetic abundances [5, 9]. This marks an advantage that facilitates the comparison of disparate communities and the formation of universal ecological hypotheses, from individual-based analysis to CATs analysis [10, 11].

These findings present valuable methodology and hypotheses for identifying core features that soil microbes may take to adapt to antibiotic pressure. Firstly, resistance to antibiotics is not equated with stress tolerance, although the ARG expression serves as a bacterial defense against antibiotic stress. ARGs predominantly encode specialized functions against antibiotics, featuring mechanisms like efflux pumps, antibiotic-inactive enzymes, and protection of drug-binding targets, while stress tolerance traits cover a broader range of defense mechanisms, including the synthesis of the molecular chaperone [12] and sigma factors [13], enhancement of surveillance and damage repairing. Moreover, a nuanced interplay between specific antibiotic resistance and universal stress tolerance strategies is illustrated by the overlap in mechanisms such as reducing permeability, forming resting structures [11], stress sensing, signaling transduction, and the expression of global transcription regulators [14,15,16]. Second, the rising abundance of ARGs is closely tied to the growth yield in terms of bacterial proliferation. Due to this, it is considered that capacities for central carbohydrate metabolism, synthesis of the ribosome, and other precursor biomolecules (e.g., amino acids, fatty acids, nucleotide) are essential for the proliferation of ARGs. Additionally, ribosomal RNA operon copy number (rrn) in bacteria is also reported as a principal trait capturing the rate-yield trade-off [17]. Third, ARGs are essentially DNA fragments whose synthesis depends on the acquisition and conversion of deoxyribose, purines, pyrimidines, and phosphate. Thus, the degradation of complex substrates and the uptake of precursor molecules is crucial for the building of ARGs. Additionally, foraging traits including motility and chemotaxis processes also play an important role in perceiving and utilizing resources, especially in environments with great heterogeneity in resource accessibility [18, 19]. Fourth, as a complement to the Y-A-S frame, the energy-producing traits were emphasized [9, 20]. Indeed, the production and allocation of energy are crucial for the proliferation of ARGs, as it is recognized that bearing ARGs requires additional reproductive costs [21,22,23].

Since antibiotic resistance is deeply intertwined with a range of physiochemical activities in microbes, it cannot be fully captured by a single aspect of features, the core microbial traits that promote ARGs growth were still poorly understood. To this end, the main objective of this study is to reveal the taxonomic and functional-genetic features that support the growth of ARGs, especially in the case of nutritional scarcity. To achieve this objective, changes in complex microbial communities and their resulting impact on ARGs were investigated along a gradient of oxytetracycline (OTC), ranging from no antibiotic presence to a clinically relevant concentration. Subsequently, the key ARGs subtypes that drive the growth of ARGs will be identified and correlated with taxonomic features at the community level and metagenomic assembly genome level. Ultimately, applying CATs analysis to reveal genetic features under conditions of ARGs growth, if present. This study contributes to a broader application of trait-based assessments to understand the causes of antibiotic resistance, leading to a more coherent theory to block the propagation of ARGs in the environmental microbiota.

Materials and methods

Creation of soil suspension and experimental design

To investigate the impact of antibiotic pressures on natural microbial communities, a soil suspension was prepared using bulk soil mixed saline solution (0.85%, w:v = 1:20). The chosen soil with no prior antibiotic exposure or residues was collected from Xixi National Park, Hangzhou, China at a site (120°3′E, 30°16′N), according to our previous studies [24, 25]. It was characterized as a loamy sand texture with pH = 6.8 ± 0.2 (1:5 soil: water mixture), and a carbon to nitrogen (C/N) ratio of 9.44, containing 3.02 g/kg of carbon and 0.32 g/kg of nitrogen.

A concentration gradient of OTC was designed in our previous study [24, 25]. OTC is a broad-spectrum tetracycline antibiotic widely used in human medicine, animal husbandry, and aquaculture [26, 27]. It is one of the main antibiotics detected in farming environments and agricultural soils, effective against both Gram-negative and Gram-positive bacteria [28]. Moreover, its role in driving bacterial resistance has been extensively studied, providing a valuable reference point for broader research [29]. In our previous study, it was found that soil microbial communities exhibited a greater ability to resist OTC compared to other tetracyclines, such as tetracycline, doxycycline, and tigecycline [25]. Building on this observation, the present study continued to focus on OTC to investigate how antibiotic concentrations influence microbial life-history strategies and their subsequent effects on ARG abundance. To obtain more precise concentration-dependent changes, two more concentrations (i.e., 0.5 mg/l and 5 mg/l) were added to the middle of the previous concentration interval of 0.1 mg/l, 1 mg/l, and 10 mg/l, approximately. Unlike premixing 2% LB broth previously, no extra medium was supplied in this study to avoid the potential impacts of exogenous nutrients. Accordingly, the soil suspension was subjected to six groups spiked with designed amounts of OTC, which were antibiotic-free (CK), 0.1 mg/l, 0.5 mg/l, 1 mg/l, 5 mg/l, and 10 mg/l. Every group was replicated on 12 microcosms, of which 3 microcosms were taken out as parallels at every time point of 1, 4, 8, and 24 h, after the initiation of cultivation. All 72 microcosms were set up as 20 ml cultures in 50 ml Mini BioReactor Tubes (Model No. 788211, NEST, China), with vent caps (0.22 µm hydrophobic filters) for gas exchange without contamination. The cultivation condition was kept constantly as shaking at 150 rpm, a temperature of 25 ℃, and in darkness. Three extra microcosms were prepared identically to CK but sampled before cultivation and used as the initial sample. The workflow is illustrated in Fig. S1.1 in Supplementary file 1: S1.

DNA extraction and 16S rRNA gene quantification

To enhance the representativeness of the sequencing result, equal volumes (i.e., 15 ml) of three parallels per sample were pooled together for DNA extraction. Then the pooled sample was centrifuged at 3200 × g, 4 ℃ for 10 min, and the supernatant and cell pellets were collected separately. Cell pellets were subjected to DNA extraction with a DNeasy PowerSoil Pro Kit (Cat. No. 47014, QIAGEN, Germany). Extracted DNA samples, supernatants, and remaining samples were frozen at − 20 ℃ until further utilization.

The amount of 16S rRNA gene per sample was quantified with a real-time qPCR approach [24]. Briefly, 2 µl DNA samples were measured with a 25-µl array which was prepared as the introduction of manufacture (TB Green Premix Ex Taq II, RR820A, TaKaRa, Japan). Triplicates were measured for each sample. qPCR was performed on CFX96™ Real-Time System (Bio-RAD, USA), and initial start with 95 ℃ for 2 min, and 2-step amplification was achieved by 40 cycles of denaturation 5 s at 95 ℃ and annealing 30 s at 60 ℃ [30]. The amount of the 16S rRNA gene per sample, as copies/µl DNA, was converted to copies/ml culture by the assumption that DNA from all suspensions was equally and completely extracted. The quantification result of 16 rRNA genes per sample is shown in Figure S2.1

Bacterial 16S rRNA amplicon sequencing and community analysis

V3–V4 region of the bacterial 16S rRNA genes was amplified and sequenced at Guangdong Magigene Biotechnology Co., Ltd. (Guangdong, China) on the Illumina Nova 6000 platform using a 250-bp paired-end protocol. The primer pair 338F (5′-ACTCCTACGGGAGGCAGCA-3′) and 806 R (5′-GGACTACHVGGGTWTCTAAT-3′) were used which were well employed to profile bacterial identities [31, 32]. Raw sequencing data was cleaned by using Fastp [33]. Afterward, amplicon sequence variants (ASVs) were generated with QIIME2 software by using the DADA2 algorithm. After removing ASVs that only appeared in one sample, ASVs were annotated against the SILVA database [34] (Version 138). In addition, rrn numbers were calculated per each ASV with the web service of rrnDB [35]. After rarefaction, a normalized abundance table with 72,260 reads per sample, 6752 ASVs in total, and rrn number per ASV was generated (Supplementary file 2: Table S1). In this study, the bacterial community composition was primarily analyzed at the family level, due to its superior resolution in characterizing community composition compared to other levels (Fig. S3.1 in Supplementary file 1: S3).

To identify the representative bacterial families in communities exposed to different OTC concentrations. A ternary diagram was generated with “ggtern” package. While samples were separated into three groups (G1, G2, and G3) according to inter-community compositional dissimilarities. Then the representativeness of a bacterial population across groups depends on its proportion, yielding three scenarios: (1) predominant in one group with a proportion > 50% indicating a single group specificity; (2) common in two groups with a proportion range of 25–50% meanwhile the third < 25%, suggesting its shared specificity in first two groups; (3) presence of 25–50% in all groups, denoting non-specificity to any group.

A correlation network was constructed to identify bacteria facilitating the increase of ARGs abundances, while samples from six groups at 24 h were used. Pearson correlations between key ARG subtypes and bacterial families were determined upon their changes in quantity by using Hmisc package and the p-value was adjusted with the Benjamini and Hochberg method [36]. The positive correlations between ARGs and bacteria were compatible with our aim (i.e. simultaneous increase or decrease), and only significant correlations were kept for the network construction (r > 0.8, p < 0.01). The network was visualized in Gephi (version 0.10.1) [37], in which modularity was then defined for every node (representing an ARG or a bacterial family). The entire network was unfolded into smaller modules that are internally cohesive yet externally distinct, as illustrated in Supplementary file 1: Fig. S8.1.

Metagenomic sequencing and genome assembly

Metagenomic sequencing was performed on 13 samples (one for initial, six for 8 h, and six for 24 h) since only a few changes per community in the structural composition and 16S rRNA copies were observed before 8 h. Gene library preparation and shotgun sequencing were ref to [38, 39]. Raw sequencing data of 1.08 × 1010 bases and 7.21 × 107 reads on average were obtained (n = 13). After removing adapters and poor quality reads (quality score < 20 and length < 50 base pair) with Trimmomatic (Ver. 0.36) [40], Clean Data were generated for subsequent use, there were 9.30 × 109 bases and 6.27 × 107 reads on average, with Q20 = 100% and Q30 > 99.69%. Details about the raw data and data cleaning per sample were given in Supplementary file 2: Table S2.

Ref to [39] the metagenomic assembly genome (MAG) was constructed upon clean data by de novo assembly. Scaffolds were assembled from short reads with MEGAHIT (v1.0.6, parameters: k-min 35, k-max 115, k-step 20) [41]. Then strictly dissociated scaffolds at unknown bases to clean chimera, and contiguous sequences with length > 500 bases without unknown bases were used as scaftigs. The summary of assembled scaftigs per sample is listed in Supplementary file 2: Table S3. The scaftigs underwent metagenomic binning utilizing MetaBAT2, employing a sensitive parameter setting as described by ref [42]. The quality of the resulting MAGs was assessed with CheckM [43] and the outcomes are shown in Supplementary file 2: Table S4. Scaftigs were also employed to predict unigenes (i.e., non-redundant gene clusters, Supplementary file 2: Sheet S5). These unigenes were subsequently mapped with functional databases, to acquire genetic traits and their abundance per sample to calculate CATs. Details about the processing of unigene prediction and functional database mapping were given in Supplementary file 1: S4.

Identification and quantification of ARGs with shot-gun data

The normalized abundance of ARG was determined by mapping the clean data against the SARG v2.0 database which comprises 12,307 unique ARG sequences [44]. Then the ARG-like sequences per sample were identified by using the ARGs-OAP automatic pipeline (https://github.com/xinehc/ARGs_OAP), with performance parameters set according to the recommended guidelines (alignment length > 75 nucleotides, e-value < 10−7, identity > 80%). The abundance of ARGs at the genotype level, as well as their corresponding antibiotic categories (i.e., ARG type) and assigned gene affiliations (i.e., ARG subtype), were generated. The normalized abundance of ARGs was calculated using three incorporated approaches, such as parts per million (ppm), copies per 16S rRNA gene copy (short as copies/16S in following), and copies per cell (the outcomes presented as Tables S1–S4 in Supplementary file 4). The relationship between bacterial quantity and the relative abundance of ARGs per sample is shown in Fig. S5.1 in Supplementary file 1. To indicate the absolute changes in the abundance of ARGs, the normalized abundance of ARGs (per 16S rRNA gene) was further weighted by the quantification of the 16S rRNA gene (results documented in Supplementary file 2: Sheet S10).

SARG database including both acquired and intrinsic ARGs. This would facilitate the distinction between intrinsic housekeeping ARGs and those that have been mobilized [45, 46], providing a clearer understanding of their respective roles and potential contributions to the observed increase in abundance.

Determination of key ARG subtypes that respond to antibiotic stress

Collectively, the key ARG subtypes were defined as fast-growing subtypes that also significantly contribute to the increase of total ARGs abundance and resistome variation. The fast-growing ARGs were termed in this study, as the ARG type or subtype (1) with a growth rate that exceeds bacterial growth and (2) with an absolute abundance greater than it in antibiotic-free control. The period growth rate (8 to 24 h, \({\mu }_{8-24}\)) was calculated to indicate the growth of ARGs or 16S rRNA genes (see Eq. S6.1 in Supplementary file 1: S6.1). Principal coordinate analysis (PCOA) was employed to discern the differences in resistome among samples, while random forest (RF) analysis was utilized to identify ARG subtypes that significantly contribute to the increase of total ARGs abundance and resistome variation, respectively. Collectively the key ARG subtypes were screened out with a custom R script (Supplementary file 5) and more details can be seen in Supplementary file 1: S6.

Determination of bacterial CATs and the global trait space

To evaluate CATs and allocate them to varied LHS [5, 47], three functional databases (i.e. KEGG, eggNOG, and CAZy) were mapped, and data were aggregated at the (1) KEGG pathways, (2) eggNOG functional categories, and (3) CAZy modules, outcomes were shown in Supplementary file 2: Tables S6, S7, and S8, respectively. A fourth database was defined with bacterial traits that have been evaluated in life-history theory [5]. Where seven traits were calculated per sample as relative abundances of genes associated with (1) sigma factor, (2) exopolysaccharides (EPS), (3) siderophores, (4) molecular chaperones, as well as the amount of (5) 16S rRNA gene copies, (6) average rrn, and (7) sum of relative abundance of ARGs (Supplementary file 2: Table S9). The full list of CATs used in this study is given in Supplementary file 3.

Then the multivariate axes that best explain the global variation in CATs were identified according to ref [5]. A multiple co-inertia analysis (MCOA) was performed with the R package “ade4” [48]. It is an exploratory analysis that leveraged together 121 CATs from 4 databases (7 defined LHTs, 11 eggNOG categories, 96 KEGG pathways, and 6 CAZy modules). By this, the co-relationships between the different databases were identified and a covariance optimization criterion to show a common structure of the information shared by multi-table multivariate was generated, as a global trait space.

Sample coordinates on the first and second dimensions of the MCOA were extracted (i.e., MCOA1 and MCOA2) and used to represent community samples in the global trait space (Fig. 4). Random forest (RF) models were then used to identify the important predictors for these coordinates. The RF analysis was performed with “randomForest” package and the significance of the importance metrics of RF was calculated with “rfPermute” package (setting 500 decision trees, 1000 times random permutations). Estimate importance per predictor with significance (p < 0.05) was plotted for MOCA1 and MCOA2, as Figures S7.3 and S7.4 of Supplementary file 1.

Analyses of the phylogenetic and functional features of MAGs

The good-quality MAGs (completeness > 60%, contamination < 10%) were selected for further taxonomic assignment and functional gene annotation. Phylogenetic analyses of these MAGs were performed using the GTDB-TK [49] (version 2.1.1) with the taxonomic classification (Supplementary file 2: Sheet S11) visualized using iTOL. The sequencing coverage (depth) of MAGs per sample was documented in Supplementary file 2: Table S12. Functional genes, i.e. ARGs and KEGG Orthology (KO) terms, presented in high-quality MAG were identified with the DIAMOND software (e-value < 10−3, pident > 0.6), and outcomes were compiled in Supplementary file 2: Sheet S13 and S14, respectively. Enriched KEGG pathways were estimated through profiled KO terms using the clusterProfiler package in R [50], which aimed at elucidating the metabolic features of the high-quality MAGs.

Results

Non-monotonic changes in ARGs along the OTC concentration gradient

A holistic perspective of changes in ARGs along the OTC concentration gradient was taken at first. The absolute quantities of total ARGs abundance (Fig. 1A), as well as the normalized abundance of ARGs per 16S rRNA gene (Fig. 1B), were plotted. Compared to the initial conditions, markable growth of ARGs was only observed between 8 and 24 h with spiked OTC concentrations not exceeding 1 mg/l. It should be noted that in the antibiotic-free CK, the absolute counts of ARGs also raised to 6.38 × 107 copies/ml, and the relative abundance increased to 0.71 ARG copies/16S rRNA gene copy, which is about 16.8 times and 4.7 times the initial, respectively. Moreover, the peak of relative enrichment was observed at 0.1 mg/l, with 0.82 ARG copies/16S copy, a 16.4% increase over CK. Meanwhile, due to a greater number of measured 16S rRNA genes (Fig. S2.1), the highest absolute abundance of ARGs was observed at 0.5 mg/l (7.19 × 107 copies/ml), which marked a 12.7% rise compared to CK.

Fig. 1
figure 1

Changes in antibiotic resistance genes along the antibiotic concentration gradient. Total abundances of ARGs in samples are plotted as A absolute quantification and B relative quantification (per 16S rRNA genes); The height of the shaded area indicates the abundance of ARGs in the initial sample, while the empty bar and filled bar represent samples taken at 8 h and 24 h, respectively. C The richness of ARGs was counted as the number of unique genotypes of ARGs per sample. D The evenness of ARGs was calculated using the Pielou index, which considers both the number and abundance of unique genotypes of ARGs per sample. The horizontal reference line indicates the level of richness and evenness in the initial sample. E Resistome per sample profiled with the absolute abundance of ARG types, corresponding to antibiotic categories. TRG is short for transcriptional regulator gene; MLS is short for macrolides, lincosamides and streptogramines

Full size image

In terms of α-diversity, richness (i.e., counts of unique ARG per sample, Fig. 1C) generally showed an increasing trend while evenness decreased (i.e., calculated with Pielou index, Fig. 1D). Numerically, after 24 h of cultivation, variations in both diversity indices were less pronounced at 5 mg/l (< 35%) and 10 mg/l (< 10%) compared to other conditions, which was in line with the observed changes in quantities of ARGs. The distinction lay in that the richness of ARGs across the antibiotic concentration gradient displayed an unimodal pattern (Fig. 1C), reaching its maximum at 0.5 mg/l with 1157 unique genotypes, which is 75 more than CK, and 2.09 times the initial counts. Conversely, the evenness slightly oscillated in the range from CK and 1 mg/l, presenting values around 0.707 ± 0.01, and thus showing a basin-like pattern (Fig. 1D). The diversity results indicated that abundances of certain ARGs increased faster than general, leading to their dominance within the community and consequently reducing evenness.

To further investigate which kind of ARGs drove the changes in ARG quantity and α-diversity, 2090 ARG genotypes were classified into 226 ARG subtypes (i.e., named functional genes) across 21 categories of antibiotics (Supplementary file 4). The top abundant 15 antibiotic categories of ARGs are displayed (Fig. 1E), and the remaining six categories together accounted for less than 1% quantity per sample. Intuitively, ARGs related to multidrug and bacitracin were identified as playing a pivotal role in the growth of the total abundance of ARGs, which could be underpinned by their interval growth rates. Absolute abundances of multidrug- and bacitracin-resistant genes grew faster than those of total ARGs, at OTC additions lower than 1 mg/l (µ8-24, Table S6.1 in Supplementary file 1). Interestingly, the growth rate of tetracycline-resistant genes was not faster than that of total ARGs in any scenario, although microbes were exposed to OTC.

At the subtype level, with finer resolution, we identified 127 fast-growing ARG subtypes, 46 important ARG subtypes for absolute abundance of ARGs and 58 important ARG subtypes for the structural variation of resistome (listed in Table S6.1, Supplementary file 1: S6) and 19 key ARG subtypes were screened out by cross-referencing as an intersection of the three lists (Table 1). Although the key ARGs were mainly dominated by functional ARGs, a large number of acquired ARGs were shown in the list of fast-growing and important ARGs. For instance, ARGs such as tet40, tet43, tetG, tetP, tetR, tetV, tetX2, tetX3, and tetX6, were highly relevant to the OTC resistance.

Table 1 Key subtypes driving ARG abundance growth. Corresponding antibiotic categories, resistance mechanisms, and abilities to resist OTC are noted accordingly
Full size table

Key ARGs represented 8.41% counts of all 226 ARG subtypes, yet accounting for 46.31% of the total abundance of ARGs. Corresponding to the specificity of the used antibiotic, ten of the key ARGs have been identified to mediate resistance to OTC. Moreover, most of the key subtypes (ca. 84.2%) are related to multidrug efflux pumps, whose encoding genes are typically located on the bacterial chromosomes and phylogenetically conserved [51]. This result indicated that changes in bacterial quantity and taxonomy might promote the observed increase in ARG abundance.

Representative bacterial taxa at different OTC concentrations

To capture the taxonomic profile relating to changes in resistome, the compositional features of communities at different antibiotic concentrations were investigated (Fig. 2A) and then identified representative bacterial populations regarding their relative abundance across groups (Fig. 2B).

Fig. 2
figure 2

Changes in bacterial community composition along the OTC concentration gradient. A Bacterial communities in 24-h samples are profiled based on the relative abundance of families. Abundant bacterial families (i.e., with average relative abundance > 1%) are shown in color otherwise are gathered as “others”, “uncultured”, and “unclassified”. The six samples are grouped into three clusters (G1, G2, and G3) based on similarities in community composition, as evaluated using the Bray–Curtis index. B representative bacteria per each group are plotted on a ternary diagram. Each dot represents a bacteria family, and its abundance value is adjusted relative to the total abundance across all three groups, i.e., adjusted abundances per family in the three groups always sum to 100%. A dot located in triangles near corners indicates a single group specificity; located in the center triangle denotes non-specificity to any group; located in the remaining three small triangles indicates shared specificity in two groups. C networks constructed upon significant positive correlations between key ARG subtypes and bacterial families. Nodes in green are key ARG subtypes, in orange are bacterial families that have positive correlations with ARGs, and in grey are other families. Key ARG subtypes and abundant families are labeled. The layout of networks is identical to Fig. S8.1 in Supplementary file 1 where five modules (#M1-5) with more than two nodes are defined

Full size image

According to the compositional similarities (Fig. 2A), six samples from 24 h were gathered into three groups, i.e., G1, G2, and G3 for the low, mediate, and high OTC stress, respectively. In G1, the representative bacterial populations were Bacillaceae, Micrococcaceae, and Rhizobiaceae (Fig. 2B). The prevalence of Bacillaceae in G1 was highlighted with an average proportion of 10.1%, which sharply contrasted to its lower observations at 0.6% and 0.3%, for G2 and G3 respectively. In G2, Pseudomonadaceae, Burkholderiaceae, and Microbacteriaceae were representative, with Pseudomonadaceae having a higher average abundance, accounting for 25.6% in G2 and 15.8% across all samples. Streptomycetaceae is the most representative species in G3 (52.7%). Furthermore, four bacterial families, Planococcaceae, Intrasporangiaceae, Nocardioidaceae, and Xanthomonadaceae, were noticed as common populations shared by G1 and G2 (i.e., G1/2 region, in Fig. 2B).

Subsequently, networks showcasing positive correlations between key ARGs and bacteria were constructed to associate enriched bacteria per treatment with ARGs growth (Fig. 2C). The representative families well exhibited across modules (i.e., #M1-2), in line with the ARGs growing (Fig. 2C). In particular, within module #M1 only representative bacteria in G1 were included, where Micrococcaceae linked two external ARGs unassociated with OTC resistance, indicating these species contribute to bacterial growth with minimal influence on ARGs. Module #M2 encompassed all 19 key resistance genes, playing a central role in the increase of ARGs. This module extensively connected ARGs with common families shared by G1 and G2, accounting for 78.1% of all bacteria-ARG links, though over half absent OTC resistance (Group G1/2, Supplementary file 2: Table S19). In addition, three representative bacterial families for G2 occupied another end of #M2 highly connected with OTC resistance genes (with an average OTC resistance gene linkage rate of 84.9%), which contributes to the growth of community size as well as to specific ARGs under antibiotic selection pressures. Phylogenetic analysis of the contigs carrying ARGs could underpin this finding (Supplementary file 1: Fig. S9.1). The remaining modules (#M3-5) did not identify correlations with increases in key ARG abundances, and the included microbes (i.e. Streptomycetaceae, Vicinamibacteraceae, and Xanthobacteraceae) might correspond to the minimal growth of resistance genes under high selection pressure.

To further illustrate the differences in life-history strategies among representative microbes identified under varying OTC concentrations, MAGs were employed to investigate the relationships between phylogenetic classification, ARG profiles, and metabolic traits at the assembled genome level. Twenty-eight high-quality MAGs were obtained from de novo binning. These MAGs were characterized based on their taxonomic features and ARG contents, along with their variations in abundance across samples (Fig. 3). MAG1 and MAG11 were identified to be genera Pseudomonas and Neobacillus, respectively, those were phylogenetically belonging to the representative families, Pseudomonadaceae for G2 and Bacillaceae for G1 (Fig. 3A and Fig. 2B). And quantitatively, the sequence depth of these two MAGs across different OTC concentrations (Fig. 3B) align with the relative abundance trends of the representative microbes (Fig. S3.2 in Supplementary file 1). Thus, the enrichment of KEGG pathways for these two MAGs was performed to provide a direct comparison of their phylogenetic and metabolic features (Fig. 4C). MAG1 has relatively greater value in genome size (4.83 Mbp), ARG content (10 ARGs), and mapped KO terms (2273 KOs) than MAG11 (3.56 Mbp, 2 ARGs, and 889 KOs, respectively), but the enriched gene pathways of MAG11 suggested a stronger proliferation ability (Fig. 3C). Particularly in carbon assimilation, ribosome synthesis, and DNA replication, as well as a bigger rrn which are closely related to the fast bacterial proliferation. In contrast, the enrichment in the oxidative phosphorylation pathway of MAG1 stood out and indicated a stronger energy production capacity. Additionally, enriched gene pathways of MAG1 also show stronger capacities in xenobiotics biodegradation and metabolism, bacterial secretion, and biofilm formation (Table S16), indicating advantages in metabolizing complex substrates.

Fig. 3
figure 3

Taxonomy, ARGs carried, and metabolic features of metagenome-assembled genomes. A 28 metagenome-assembled genomes (MAGs) with binning completeness above 60%. Their taxonomic classification at levels of phylum (inside circle) and genus (outside circle) are shown as colored squares, while other level classifications are listed in Supplementary file 2: Table S11. And the number of ARGs carried by each MAGs is indicated by the length of the internal red bars. B Normalized sequence depth per MAG. Based on the changing pattern of sequencing depth along with the antibiotic concentration increase, all MAGs were divided into four groups, i.e., Cluster1: decreasing; Cluster2: decreasing follow increase; Cluster3: almost unchanged; Cluster4: increasing gradually. C Enriched KEGG pathways and their enrich folds on MAG1 (orange) and MAG11 (green)

Full size image
Fig. 4
figure 4

Life history strategies of bacterial communities under antibiotic treatment. Two-dimensional trait space from an MCOA depicting community aggregated traits across OTC-treated bacterial communities, with traits inferred from enriched genes in bacteria metagenomes. Dots represent the positions of 13 sampled bacterial communities used in this study along these two dimensions. In the trait lists (Supplementary file 3: Table S1), letters in brackets represent how life history traits (high growth yield, resource acquisition, stress tolerance) have been associated with classical theoretical works [5, 10, 11, 20]

Full size image

Predominant CATs at different OTC concentrations

To explore microbial features associated with an increased abundance of ARGs, genetic traits were explored at the community level per OTC concentration. The multiple co-inertia analysis (MCOA) was performed accordingly [5]. As shown in Fig. 4, two dimensions captured most of the overall variation in metagenomic CATs. The first axis MCOA1 emerged as the most dominant, capturing 92.19% of the co-inertia variance, and the second axis MCOA2 served as a complementary dimension, accounting for 6.42%. A triangle of sample distributions dominated by these two dimensions emerged, by which the feature of CATs was illustrated.

Samples with relatively high ARG abundances (i.e., CK and 0.1–1 mg/l OTC at 24 h) were found to cluster at the upper end of the MCOA1. To elucidate this pattern, the quantity of ARGs for each sample was plotted against its coordinates within the trait space (Fig. S7.2 in Supplementary file 1). A significant positive correlation was observed between the ARG abundance and projections on the MCOA1 (Pearson’s r = 0.98, p-value < 0.05), but not with the MCOA2 axis (Pearson’s r = 0.14). This indicates that CATs which significantly contribute to the rise in MCOA1 values are also key determinants of the growth in ARGs abundance.

RF models were used to identify CATs that significantly contributed to co-inertia values of MCOA1 and MCOA2 (Figs. S4.1 and S4.2, respectively). At the upper end of MCOA1, microbial communities with higher ARG abundances demonstrated enhanced bacterial growth and resource acquisition capabilities (Fig. S4.1). This was presented by the higher abundance of genes related to energy production and conversion, central carbon metabolism, and carbon fixation pathways. Additionally, genes involved in exopolysaccharide (EPS) production, signal transduction, membrane, and DNA repair were also enriched, indicating a stress-tolerant capacity to withstand environmental stress. In contrast, samples in the lower end of the MCOA1 (Fig. 4) showed a decrease in these genetic traits but an increase in defense mechanisms, mRNA surveillance, and molecular chaperone synthesis. Overall, the first trait dimension captured variation associated with resource acquisition and productivity which supported the growth of bacteria and ARGs.

Microbes differentiated along the dimension MCOA2, particularly at the negative zone of MCOA1, reflecting the adaptation of communities from native soil to cultivated conditions (Fig. 4). At the low MCOA2 region, enhanced environmental sensing and tolerance were evident, marked by a high abundance of genes for sigma factors and defense mechanisms. As MCOA2 increased, genes related to cell cycling regulation, lipid metabolism, and the pentose phosphate pathway became more abundant, enhancing the capacities for reduction power and intermediary metabolite production. Moreover, genes related to EPS synthesis, glycosyltransferases, and glycosidases increased, indicating a high capacity for sugar uptake and utilization.

Focusing on samples with a relatively higher abundance of ARGs, and using CK as a reference with a 5% change threshold, we identified genetic traits with increased and decreased numbers to clarify traits mediating simultaneous increases in ARGs and bacteria (Table S15 in Supplementary file 2). Under the condition of the highest ARGs abundance (i.e., 0.5 mg/l OTC), further increases in gene abundance related to traits of energy production, damage repair, biofilm formation, and ribosome synthesis were found, while streamlining genetic traits involved in complex sugar and lipid metabolism, active transport, motility and chemotaxis. This demonstrates that under conditions of antibiotic pressure and resource scarcity, the enhancement of energy production and utilizing efficiency leads to further growth in both microbes and ARGs, compared to the intrinsic resistance expansion detected in the antibiotic-free condition.

Discussion

Collectively, the main results obtained from this study are summarized in Table 2, providing a comprehensive view.

Table 2 Summary of all observations in this study
Full size table

Our findings illustrate the amplification of ARGs by microbial communities is a nuanced response to the antibiotic pressure, promoted by both the relative enrichment of ARGs and bacterial growth. At this point, both the relative and absolute abundances of ARGs display a non-monotonic growth trend across the antibiotic gradient; however, due to variations in bacterial quantities, their peak values manifest at different concentrations (Fig. 1). Unlike the clinical focus, which primarily concentrates on the antibiotic resistance capabilities of specific pathogens [52, 53], the environmental viewpoint, guided by the One Health approach, emphasizes understanding the overall presence of ARGs within a region and the associated risks [4, 54]. To this end, concentrating exclusively on the changes in relative enrichment would result in biased perceptions of the true situation. Quantification through the use of tagged gene spiking is acknowledged as a more accurate approach, though it is more demanding [55]. The method of weighting by 16S rRNA gene copies to characterize changes in the absolute abundance of ARGs is also deemed appropriate in this study, given that DNA recoveries from samples were expected to be essentially identical.

It is also important to clarify that the abundance of ARGs does not equate to antibiotic resistance capacity. The mere presence of ARGs does not necessarily indicate active expression. Intrinsic functional ARGs, chromosomally carried by the host microbes may remain inactivated or unexpressed due to insufficient external stimuli, which contrasts with the expression efficiency of acquired ARGs carried on plasmids or integrons [4]. To serve different research purposes, various ARG databases have been developed, such as ResFinder, which primarily focuses on acquired resistance [46], ResFinderFG for functional ARGs [45], and the CARD database combines both [56], however, no single database is suitable for all purposes. In this study, we employed the SARG database, which encompasses both intrinsic and acquired resistance genes, as its automated pipeline for short-read analysis is more aligned with the need for absolute quantification of ARG abundance.

With the absolute quantification of ARGs, an interesting phenomenon was observed under intermediate antibiotic stress (i.e., 0.5 mg/l), there was a simultaneous increase in bacterial quantity and the relative abundance of ARGs compared with CK (Table 2). The explanation for a similar phenomenon is that low dose antibiotics induce a hormesis effect, meaning bacteria overcompensate for the inhibition of antibiotics [57]. In addition, antibiotics at sub-lethal levels are known for accelerating the emergence and spread of antibiotic-resistant bacteria, which act as signaling molecules in promoting biofilm formation and bacterial virulence [3]. These perspectives have been preliminarily validated using bacterial populations and simplified microbial communities [58, 59]. However, their applicability to complex microbial communities remains uncertain, as fluctuations in the abundance of ARGs may result from taxonomic shifts rather than subinhibitory antibiotic pressure [4]. By synthesizing findings from representative bacteria and CATs classification (Table 2), it could be proposed that this further increase in both bacteria and ARGs is the consequence of shifts in the LHS driven by microbial communities.

The simultaneous growth is not due to opportunist or passive stress tolerance; instead, it signifies a proactive succession of microbes to antibiotic selection pressure. From the perspective of representative bacteria, the transition from Bacillaceae in G1 to Pseudomonadaceae in G2; indicates shifts from the ruderal to a competitive lifestyle [60]. Generally, bacteria from Bacillaceae were considered to have small genome size, high maximum growth rate, and great ability to germinate from spores within a short period [61], ensuring its opportunistic growth under antibiotic-free conditions; in contrast, Pseudomonadaceae is featured by their high catabolic diversity, siderophore production and capacity to form biofilms, enhancing their survival competitiveness in stressed environments. The Pseudomonadaceae family is also known for its substantial intrinsic resistance to tetracyclines and toxins [14, 62]. This resistance is primarily mediated by the expression of a diverse array of functional ARGs that encode multidrug efflux pumps, which actively expel antibiotics from the bacterial cell [51], including mexAB-oprM and mexEF, as detected on MAG1. These characteristics are consistent with our findings obtained from MAGs (Fig. 3; Tables S11, S13–14, S16). Additionally, enrichment of genes in the oxidative phosphorylation pathway of MAG1 (belonging Pseudomonadaceae) indicates an advantage in energy production trait (Fig. 3C). From the CATs analysis, the end with higher quantities of bacteria and ARGs demonstrates stronger substrate acquisition capabilities and rapid proliferation traits; build on this, enhanced energy production and substrate utilization efficiency emerge as primary genetic traits that support the further increase in the abundance of ARGs, which is reflected by the strengthening of central carbon metabolism and the streamlining of low competitiveness function genes. Moreover, key ARGs encoding tetracycline resistance (Table 1), linking with the representative bacteria in G2, become prominent in the positive correlation networks (Fig. 2C), which demonstrates the active adaptation of the community to increased antibiotic selection pressure.

In this study, a genetic trait-based assessment approach was utilized to reveal altered LHSs under various antibiotic selection pressures. It further extends the application of this approach in the development of universal ecological theories. We mainly followed the Y-A-S theoretical framework [10], while taking into account the effects of energy production [20], motility, and chemotaxis on opportunists [11], and using MCOA as an integrating statistic analytical tool to describe aggregated genetic traits at the community level [5]. Although the development of these theories and methods was not specifically aimed at antibiotic resistance, they greatly enhance our understanding of how microbial communities balance proliferation with responses to environmental stresses. To obtain a more coherent theoretical framework for the antibiotic resistance issue, the focus could be placed on traits in terms of accessibility and mobility of ARGs [4], and the distinction between resistance and tolerance [63], in subsequent studies. Which could be complemented by the use of a genome-informed trait-based dynamic energy budget model [9].

Further, the limits of the experimental design of this study should be stated here, which focuses more on the effects of antibiotics in resource-limited conditions to mimic the typical nature. However, under conditions influenced by human activities the dominant LHSs may vary with the substantial increase in available resources [64, 65]. In comparison with a previous study with the same soil source and experimental conditions, the maximum abundance of ARGs was observed at a higher concentration of 10 mg/l OTC with extra LB medium addition (even only 2%) [25]. This is mainly due to the increased abundance of specific efflux pumps for tetracyclines (tet39) and inactivated enzymes such as tetX. However, the extra fitness cost of carrying these two ARGs might be unsustainable for microbial communities under resource-scarce conditions. Moreover, the relatively short experimental duration imposes some limitations on the generalizability of the conclusions. While differences in microbial community composition and ARG abundance across OTC treatments were observed within 24 h, it remains uncertain whether these differences would intensify or diminish with extended reaction times. This uncertainty is particularly relevant as the selective pressure exerted by OTC may gradually decrease over time due to its slow degradation. Additionally, exploring multiple classes of antibiotics represents a crucial avenue for future research, despite the added complexity introduced by potential antagonistic or synergistic interactions between antibiotics. Notably, investigating the combined effects of tetracycline and sulfonamide antibiotics could yield meaningful insights, especially when aiming to replicate conditions commonly found in farming environments [66].

The scarcity of resources prompted another consideration in the process of controlling the spread of ARGs. That is blocking the propagation of ARGs based on the principle of competitive exclusion [67], in particular, artificially enhancing competition between alternative less-ARG populations and more-ARG populations under limited resource conditions. At the population level, it was well predicted based on the coexistence theory that a fluctuating resource environment can effectively affect the competitiveness of resistant bacterial strains [23]. Gut symbionts protect against pathogens through nutrient blocking is another experimental evidence with a more complex microbial context [68]. Returning to this study, the replacement of Bacillaceae by Pseudomonadaceae with increasing concentrations of antibiotics and ARGs caught our attention. Species from these two families both ubiquitously occur in natural environments and are often co-isolated from soil samples [60, 61, 69], indicating comparable adaptive capacity. Indeed, interactions between members of these two families range from mutualism to competition depending on cultural conditions, however, were mainly reported as either competition or amensalism [61]. Most of the Bacillaceae members have no pathogenic potential and are relatively sensitive to antibiotics [70, 71], and Bacillus subtilis, in particular, is commonly used as an intestinal probiotic and as a biocontrol bacterium for the protection of human, animal, and plant health [72,73,74]. Thus, the use of Bacillus subtilis spore preparations to deplete limited growth resources might be a viable method of controlling ARGs spread under sub-inhibitory antibiotic stress. Meanwhile, it would be an attempt to use the traits-based life history theory to guide microbial practice.

Conclusions

This study uncovers the capacity of soil microbial communities to grow ARGs across gradients of antibiotic concentrations, shedding light on the correlating bacterial taxa and their functional genetic features. Results of the analysis of the community aggregated traits demonstrated that a competitive lifestyle facilitating the growth in ARGs by enhanced energy production and resource utilization, particularly under conditions of nutrient scarcity met the sub-inhibitory of antibiotics.

Data availability

All the original sequences are accessible at CNCB-NGDC (https://ngdc.cncb.ac.cn/) with a BioProject NO. PRJCA024366. Shotgun metagenome sequences [GSA: CRA016954]; 16S rRNA amplicon sequences [GSA: CRA017053]. Other data used in the analyses are provided in the text and supporting materials.

References

  1. Murray CJL, Ikuta KS, Sharara F, Swetschinski L, Robles Aguilar G, Gray A, et al. Global burden of bacterial antimicrobial resistance in 2019: a systematic analysis. Lancet. 2022;399(10325):629–55.

    Article  Google Scholar 

  2. UNEP. Bracing for Superbugs: Strengthening environmental action in the One Health response to antimicrobial resistance. Geneva; 2023. https://www.unep.org/resources/superbugs/environmental-action.

  3. Andersson DI, Hughes D. Microbiological effects of sublethal levels of antibiotics. Nat Rev Microbiol. 2014;12(7):465–78.

    Article  PubMed  Google Scholar 

  4. Larsson DGJ, Flach C-F. Antibiotic resistance in the environment. Nat Rev Microbiol. 2022;20:257–69.

    Article  PubMed  Google Scholar 

  5. Piton G, Allison SD, Bahram M, Hildebrand F, Martiny JBH, Treseder KK, et al. Life history strategies of soil bacterial communities across global terrestrial biomes. Nat Microbiol. 2023;8(11):2093–102.

    Article  PubMed  Google Scholar 

  6. Yao X, Wang J, Hu B. How methanotrophs respond to pH: a review of ecophysiology. Front Microbiol. 2023;13:1034164.

    Article  PubMed  PubMed Central  Google Scholar 

  7. Pianka ER. On r- and K-selection. Am Nat. 1970;104(940):592–7.

    Article  Google Scholar 

  8. Grime JP. Evidence for the existence of three primary strategies in plants and its relevance to ecological and evolutionary theory. Am Nat. 1977;111(982):1169–94.

    Article  Google Scholar 

  9. Marschmann GL, Tang J, Zhalnina K, Karaoz U, Cho H, Le B, et al. Predictions of rhizosphere microbiome dynamics with a genome-informed and trait-based energy budget model. Nat Microbiol. 2024;9(2):421–33.

    Article  PubMed  PubMed Central  Google Scholar 

  10. Malik AA, Martiny JBH, Brodie EL, Martiny AC, Treseder KK, Allison SD. Defining trait-based microbial strategies with consequences for soil carbon cycling under climate change. ISME J. 2020;14(1):1–9.

    Article  PubMed  Google Scholar 

  11. Wood JL, Malik AA, Greening C, Green PT, McGeoch M, Franks AE. Rethinking CSR theory to incorporate microbial metabolic diversity and foraging traits. ISME J. 2023;17(11):1793–7.

    Article  PubMed  PubMed Central  Google Scholar 

  12. Caporaso JG, Lauber CL, Walters WA, Berg-Lyons D, Lozupone CA, Turnbaugh PJ. Global patterns of 16s rrna diversity at a depth of millions of sequences per sample. Proc National Academy Sci. 2011;108(supplement_1):4516–22.

    Article  Google Scholar 

  13. Mitchell AM, Silhavy TJ. Envelope stress responses: balancing damage repair and toxicity. Nat Rev Microbiol. 2019;17(7):417–28.

    Article  PubMed  PubMed Central  Google Scholar 

  14. Du D, Wang-Kan X, Neuberger A, van Veen HW, Pos KM, Piddock LJV, et al. Multidrug efflux pumps: Structure, function and regulation. Nat Rev Microbiol. 2018;16(9):523–39.

    Article  PubMed  Google Scholar 

  15. Janssens A, Nguyen VS, Cecil AJ, Van der Verren SE, Timmerman E, Deghelt M, et al. Slyb encapsulates outer membrane proteins in stress-induced lipid nanodomains. Nature. 2024;626(7999):617–25.

    Article  PubMed  Google Scholar 

  16. Hao Z, Lou H, Zhu R, Zhu J, Zhang D, Zhao BS, et al. The multiple antibiotic resistance regulator marr is a copper sensor in Escherichia coli. Nat Chem Biol. 2014;10(1):21–8.

    Article  PubMed  Google Scholar 

  17. Roller BRK, Stoddard SF, Schmidt TM. Exploiting rrna operon copy number to investigate bacterial reproductive strategies. Nat Microbiol. 2016;1(11):16160.

    Article  PubMed  PubMed Central  Google Scholar 

  18. Wadhwa N, Berg HC. Bacterial motility: machinery and mechanisms. Nat Rev Microbiol. 2022;20(3):161–73.

    Article  PubMed  Google Scholar 

  19. Keegstra JM, Carrara F, Stocker R. The ecological roles of bacterial chemotaxis. Nat Rev Microbiol. 2022;20(8):491–504.

    Article  PubMed  Google Scholar 

  20. Karaoz U, Brodie EL. Microtrait: a toolset for a trait-based representation of microbial genomes. Front Bioinform. 2022;2:918853.

    Article  PubMed  PubMed Central  Google Scholar 

  21. Rajer F, Sandegren L. The role of antibiotic resistance genes in the fitness cost of multiresistance plasmids. Bio. 2022;13(1):e03552–21.

    Google Scholar 

  22. Andersson DI, Hughes D. Antibiotic resistance and its cost: Is it possible to reverse resistance? Nat Rev Microbiol. 2010;8(4):260–71.

    Article  PubMed  Google Scholar 

  23. Letten AD, Hall AR, Levine JM. Using ecological coexistence theory to understand antibiotic resistance and microbial competition. Nat Ecol Evol. 2021;5(4):431–41.

    Article  PubMed  Google Scholar 

  24. Liu Z, Jin Y, Yu Z, Liu Z, Zhang B, Chi T, et al. Vertical migration and dissipation of oxytetracycline induces the recoverable shift in microbial community and antibiotic resistance. Sci Total Environ. 2023;905:167162.

    Article  PubMed  Google Scholar 

  25. Tang H, Liu Z, Hu B, Zhu L. D-ring modifications of tetracyclines determine their ability to induce resistance genes in the environment. Environ Sci Technol. 2023;58(2):1338–48.

    Article  PubMed  Google Scholar 

  26. Bacanli M, Basaran N. Importance of antibiotic residues in animal food. Food Chem Toxicol. 2019;125:462–6.

    Article  PubMed  Google Scholar 

  27. Scaria J, Anupama KV, Nidheesh PV. Tetracyclines in the environment: An overview on the occurrence, fate, toxicity, detection, removal methods, and sludge management. Sci Total Environ. 2021;771:145291.

    Article  PubMed  Google Scholar 

  28. Qiao M, Ying GG, Singer AC, Zhu YG. Review of antibiotic resistance in China and its environment. Environ Int. 2018;110:160–72.

    Article  PubMed  Google Scholar 

  29. Thaker M, Spanogiannopoulos P, Wright GD. The tetracycline resistome. Cell Mol Life Sci. 2010;67(3):419–31.

    Article  PubMed  Google Scholar 

  30. Yao X, Wang J, He M, Liu Z, Zhao Y, Li Y, et al. Methane-dependent complete denitrification by a single Methylomirabilis bacterium. Nat Microbiol. 2024;9(2):464–76.

    Article  PubMed  Google Scholar 

  31. Liu Z, Cichocki N, Hübschmann T, Süring C, Ofiteru ID, Sloan WT, et al. Neutral mechanisms and niche differentiation in steady-state insular microbial communities revealed by single cell analysis. Environ Microbiol. 2019;21(1):164–81.

    Article  PubMed  Google Scholar 

  32. Li S, Abdulkadir Nu, Schattenberg F, Nunes da Rocha U, Grimm V, Müller S, et al. Stabilizing microbial communities by looped mass transfer. Proc Natl Acad Sci USA. 2022;119(17):e2117814119.

  33. Chen S, Zhou Y, Chen Y, Gu J. Fastp: An ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–90.

    Article  PubMed  PubMed Central  Google Scholar 

  34. Quast C, Pruesse E, Yilmaz P, Gerken J, Schweer T, Yarza P, et al. The SILVA ribosomal RNA gene database project: Improved data processing and web-based tools. Nucleic Acids Res. 2012;41(D1):D590–6.

    Article  PubMed  PubMed Central  Google Scholar 

  35. Stoddard SF, Smith BJ, Hein R, Roller BRK, Schmidt TM. Rrndb: improved tools for interpreting rRNA gene abundance in bacteria and archaea and a new foundation for future development. Nucleic Acids Res. 2014;43(D1):D593–8.

    Article  PubMed  PubMed Central  Google Scholar 

  36. Benjamini Y, Hochberg Y. Controlling the false discovery rate: a practical and powerful approach to multiple testing. J R Stat Soc Ser B-Stat Methodol. 1995;57(1):289–300.

    Article  Google Scholar 

  37. Bastian M, Heymann S, Jacomy M. Gephi: An Open Source Software for Exploring and Manipulating Networks. Proceedings of the International AAAI Conference on Web and Social Media. 2009;3(1):361–2.

    Article  Google Scholar 

  38. Liu Z, Zhao Y, Zhang B, Wang J, Zhu L, Hu B. Deterministic effect of pH on shaping soil resistome revealed by metagenomic analysis. Environ Sci Technol. 2023;57(2):985–96.

    Article  PubMed  Google Scholar 

  39. Zhao Y, Liu Z, Zhang B, Cai J, Yao X, Zhang M, et al. Inter-bacterial mutualism promoted by public goods in a system characterized by deterministic temperature variation. Nat Commun. 2023;14(1):5394.

    Article  PubMed  PubMed Central  Google Scholar 

  40. Bolger AM, Lohse M, Usadel B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics. 2014;30(15):2114–20.

    Article  PubMed  PubMed Central  Google Scholar 

  41. Li D, Liu C-M, Luo R, Sadakane K, Lam T-W. Megahit: An ultra-fast single-node solution for large and complex metagenomics assembly via succinct de Bruijn graph. Bioinformatics. 2015;31(10):1674–6.

    Article  PubMed  Google Scholar 

  42. Kang DD, Li F, Kirton E, Thomas A, Egan R, An H, et al. Metabat 2: An adaptive binning algorithm for robust and efficient genome reconstruction from metagenome assemblies. PeerJ. 2019;7:e7359.

    Article  PubMed  PubMed Central  Google Scholar 

  43. Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW. CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes. Genome Res. 2015;25(7):1043–55.

    Article  PubMed  PubMed Central  Google Scholar 

  44. Yin X, Jiang X-T, Chai B, Li L, Yang Y, Cole JR, et al. ARGs-OAP v2.0 with an expanded SARG database and Hidden Markov Models for enhancement characterization and quantification of antibiotic resistance genes in environmental metagenomes. Bioinformatics. 2018;34(13):2263–70.

    Article  PubMed  Google Scholar 

  45. Gschwind R, UgarcinaPerovic S, Weiss M, Petitjean M, Lao J, Coelho Luis P, et al. ResFinderFG v2.0: a database of antibiotic resistance genes obtained by functional metagenomics. Nucleic Acids Res. 2023;51(W1):W493–500.

    Article  PubMed  PubMed Central  Google Scholar 

  46. Bortolaia V, Kaas RS, Ruppe E, Roberts MC, Schwarz S, Cattoir V, et al. ResFinder 4.0 for predictions of phenotypes from genotypes. J Antimicrob Chemother. 2020;75(12):3491–500.

    Article  PubMed  PubMed Central  Google Scholar 

  47. Bahram M, Hildebrand F, Forslund SK, Anderson JL, Soudzilovskaia NA, Bodegom PM, et al. Structure and function of the global topsoil microbiome. Nature. 2018;560(7717):233–7.

    Article  PubMed  Google Scholar 

  48. Dray S, Dufour A-B. The ade4 package: Implementing the duality diagram for ecologists. J Stat Softw. 2007;22(4):1–20.

    Article  Google Scholar 

  49. Parks DH, Chuvochina M, Rinke C, Mussig AJ, Chaumeil P-A, Hugenholtz P. GTDB: an ongoing census of bacterial and archaeal diversity through a phylogenetically consistent, rank normalized and complete genome-based taxonomy. Nucleic Acids Res. 2021;50(D1):D785–94.

    Article  PubMed Central  Google Scholar 

  50. Wu T, Hu E, Xu S, Chen M, Guo P, Dai Z, et al. clusterProfiler 4.0: a universal enrichment tool for interpreting omics data. Innovation. 2021;2(3):100141.

    PubMed  PubMed Central  Google Scholar 

  51. Henderson PJF, Maher C, Elbourne LDH, Eijkelkamp BA, Paulsen IT, Hassan KA. Physiological functions of bacterial “multidrug” efflux pumps. Chem Commun. 2021;121(9):5417–78.

    Google Scholar 

  52. Willyard C. The drug-resistant bacteria that pose the greatest health threats. Nature. 2017;543(7643):15-.

    Article  PubMed  Google Scholar 

  53. Bakkeren E, Huisman JS, Fattinger SA, Hausmann A, Furter M, Egli A, et al. Salmonella persisters promote the spread of antibiotic resistance plasmids in the gut. Nature. 2019;573(7773):276–80.

    Article  PubMed  PubMed Central  Google Scholar 

  54. Zhang A-N, Gaston JM, Dai CL, Zhao S, Poyet M, Groussin M, et al. An omics-based framework for assessing the health risk of antimicrobial resistance genes. Nat Commun. 2021;12(1):4765.

    Article  PubMed  PubMed Central  Google Scholar 

  55. Yin X, Chen X, Jiang X-T, Yang Y, Li B, Shum MH-H, et al. Toward a universal unit for quantification of antibiotic resistance genes in environmental samples. Environ Sci Technol. 2023;57(26):9713–21.

  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. 2019;48(D1):D517–25.

    PubMed Central  Google Scholar 

  57. Iavicoli I, Fontana L, Agathokleous E, Santocono C, Russo F, Vetrani I, et al. Hormetic dose responses induced by antibiotics in bacteria: A phantom menace to be thoroughly evaluated to address the environmental risk and tackle the antibiotic resistance phenomenon. Sci Total Environ. 2021;798:149255.

    Article  PubMed  Google Scholar 

  58. Kumar A, Saha SK, Banerjee P, Prasad K, Sengupta TK. Antibiotic-induced biofilm formations in Pseudomonas aeruginosa strains KPW.1-s1 and HRW.1-s3 are associated with increased production of eDNA and exoproteins, increased ROS generation, and increased cell surface hydrophobicity. Curr Microbiol. 2023;81(1):11.

    Article  PubMed  Google Scholar 

  59. Chatterjee A, Willett JLE, Dunny GM, Duerkop BA. Phage infection and sub-lethal antibiotic exposure mediate Enterococcus faecalis type VII secretion system dependent inhibition of bystander bacteria. PLOS Genet. 2021;17(1):e1009204.

    Article  PubMed  PubMed Central  Google Scholar 

  60. Fierer N. Embracing the unknown: disentangling the complexities of the soil microbiome. Nat Rev Microbiol. 2017;15(10):579–90.

    Article  PubMed  Google Scholar 

  61. Lyng M, Kovács ÁT. Frenemies of the soil: Bacillus and Pseudomonas interspecies interactions. Trends Microbiol. 2023;31(8):845–57.

    Article  PubMed  Google Scholar 

  62. Langendonk RF, Neill DR, Fothergill JL. The building blocks of antimicrobial resistance in Pseudomonas aeruginosa: Implications for current resistance-breaking therapies. Front Cell Infect Microbiol. 2021;11:665759.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Brauner A, Fridman O, Gefen O, Balaban NQ. Distinguishing between resistance, tolerance and persistence to antibiotic treatment. Nat Rev Microbiol. 2016;14(5):320–30.

    Article  PubMed  Google Scholar 

  64. Fu Y, Wang F, Xiang L, Harindintwali JD, Elsner M, Amelung W, et al. Combating antibiotic resistance in the human-impacted environment with carbon-based materials: Applications and challenges. Crit Rev Environ Sci Technol. 2024;54(9):699–721.

    Article  Google Scholar 

  65. Tiedje JM, Fu Y, Mei Z, Schäffer A, Dou Q, Amelung W, et al. Antibiotic resistance genes in food production systems support One Health Opinions. Curr Opin Environ Sci Health. 2023;34:100492.

    Article  Google Scholar 

  66. Chen Z, Zhang W, Peng A, Shen Y, Jin X, Stedtfeld RD, et al. Bacterial community assembly and antibiotic resistance genes in soils exposed to antibiotics at environmentally relevant concentrations. Environ Microbiol. 2023;25(8):1439–50.

    Article  PubMed  Google Scholar 

  67. Gause GF. Experimental analysis of Vito Volterra’s mathematical theory of the struggle for existence. Science. 1934;79(2036):16–7.

    Article  PubMed  Google Scholar 

  68. Spragge F, Bakkeren E, Jahn MT, B. N. Araujo E, Pearson CF, Wang X, et al. Microbiome diversity protects against pathogens by nutrient blocking. Science. 2023;382(6676):eadj3502.

  69. Banerjee S, Schlaeppi K, van der Heijden MGA. Keystone taxa as drivers of microbiome structure and functioning. Nat Rev Microbiol. 2018;16(9):567–76.

    Article  PubMed  Google Scholar 

  70. Harirchi S, Sar T, Ramezani M, Aliyu H, Etemadifar Z, Nojoumi SA, et al. Bacillales: from taxonomy to biotechnological and industrial perspectives. Microorganisms. 2022;10(12):2355.

    Article  PubMed  PubMed Central  Google Scholar 

  71. Mandic-Mulec I, Stefanic P, Elsas JDv. Ecology of Bacillaceae. Microbiol Spectr 2015;3(2):https://doiorg.publicaciones.saludcastillayleon.es/10.1128/microbiolspec.tbs-0017-2013.

  72. Tachibana L, Telli GS, Dias DdC, Gonçalves GS, Guimarães MC, Ishikawa CM, et al. Bacillus subtilis and Bacillus licheniformis in diets for nile tilapia (Oreochromis niloticus): Effects on growth performance, gut microbiota modulation and innate immunology. Aquac Res. 2021;52(4):1630–42.

  73. Mahapatra S, Yadav R, Ramakrishna W. Bacillus subtilis impact on plant growth, soil health and environment: Dr Jekyll and Mr Hyde. J Appl Microbiol. 2022;132(5):3543–62.

    Article  PubMed  Google Scholar 

  74. Arnaouteli S, Bamford NC, Stanley-Wall NR, Kovács ÁT. Bacillus subtilis biofilm formation and social interactions. Nat Rev Microbiol. 2021;19(9):600–14.

    Article  PubMed  Google Scholar 

Download references

Acknowledgements

The authors gratefully acknowledge the Zhejiang University's Convergence Innovative Research Program of Environmental Hi-Tech and Ecological Civilization. And acknowledge Peng Zhang from Magigene Biotechnology Ltd. (Guangdong, China) for providing technical support.

Funding

This work was supported by the National Natural Science Foundation of China [grant numbers 22193061, 22206166]; National Key Research and Development Program of China [grant number 2020YFC1806903]; Open Project of State Key Laboratory of Urban Water Resource and Environment, Harbin Institute and Technology [grant number No. ES202118]; China Postdoctoral Science Foundation awards the fellowship [2023M733056]; Zhejiang Province Ecological Environment Research and Results Promotion Project [2022HT0025].

Author information

Authors and Affiliations

  1. Key Laboratory of Environment Remediation and Ecological Health, Ministry of Education, College of Environmental and Resource Sciences, Zhejiang University, Hangzhou, 310058, China

    Zishu Liu, Lizhong Zhu & Baolan Hu

  2. College of Environmental and Resource Sciences, Zhejiang University, Hangzhou, 310058, China

    Zishu Liu, Xiangwu Yao, Chengyi Chen, Yuxiang Zhao, Chifei Dong, Lingtao Sun, Junxian Zhao, Lizhong Zhu & Baolan Hu

  3. Key Laboratory of Water Pollution Control and Environmental Safety of Zhejiang Province, Hangzhou, 310058, China

    Zishu Liu & Baolan Hu

  4. Hangzhou Ecological and Environmental Monitoring Center, Hangzhou, 310007, China

    Baofeng Zhang

  5. School of Medical Technology and Information Engineering, Zhejiang Chinese Medical University, Hangzhou, 310053, China

    Zhendi Yu & Dongqing Cheng

Authors
  1. Zishu Liu
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Xiangwu Yao
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Chengyi Chen
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Yuxiang Zhao
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Chifei Dong
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Lingtao Sun
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Junxian Zhao
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Baofeng Zhang
    View author publications

    You can also search for this author inPubMed Google Scholar

  9. Zhendi Yu
    View author publications

    You can also search for this author inPubMed Google Scholar

  10. Dongqing Cheng
    View author publications

    You can also search for this author inPubMed Google Scholar

  11. Lizhong Zhu
    View author publications

    You can also search for this author inPubMed Google Scholar

  12. Baolan Hu
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

ZL: make the conception and experimental design, bioinformatic analysis, prepared Figs. 1, 2, 3 and 4 and manuscript writing. XY: conception. CC: was a major contributor in community cultivation and sampling. CD: Phylogenetic analysis of contigs carrying ARGs. LS, JZ and BZ: experimental work. YZ: assembly-based analyses. ZY: qPCR analysis. DC, LZ and BH: supervision, manuscript writing, and editing. All authors read and approved the final manuscript.

Corresponding author

Correspondence to Baolan Hu.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare 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

40168_2024_2005_MOESM1_ESM.pdf

Supplementary Material 1: S1. Experimental settings. Figure S1.1: Schematic diagram of soil suspension preparation and microcosmic culture array. S2: Quantification of 16S rRNA gene with real-time qPCR. Table S2.1: Primer sequences and the receipt of the reaction array for the real-time qPCR quantification of 16S rRNA gene. Figure S2.1: Quantification result of 16S rRNA genes per sample. S3: Analysis of the taxonomic composition of microbial communities. Figure S3.1. Profile of bacterial communities at levels of phylum, class, family and genus. Only the most abundant 30 taxa per level were coloured and the others were summed as‘others’, ‘unclassified’ and ‘uncultured’ affiliations were shown separately. Figure S3.2 16S rRNA gene weighted abundance per family at 24 hours across different OTC concentrations. S4: Unigene prediction and functional database mapping. Table S4.1: Summary of the ORF prediction. S5: Abundance of ARG per antibiotic per sample. Figure S5.1: Relationship between bacterial quantities and relative abundance of ARGs per sample. The number of 16S rRNA gene copies was plotted against the relative abundance of ARGs, which normalized as A) ARGs per 16S; B) ARGs per million reads; and C) ARGs per cell. S6 Screening for fast-growing ARGs and key ARG subtypes. Table S6.1: Period growth rate µ8-24 (/hour) per ARG types in different OTC concentrations. Figure S6.1: Screening results of fast-growing ARGs for A) ARG types and B) ARG subtypes. Each dot on the plot represents an ARG type or subtype, the colour and size indicate its corresponding antibiotic category and absolute abundance, respectively. Figure S6.2: Important ARG subtypes for the changes of absolute abundance of ARGs. R package ‘rfPermute’ was used to calculate the IncMSE value and significance, with settings of 500 decision trees and 1000 times random permutations. The identity of ARG subtypes was shown as y-axis labels, while the colour and size of the point indicate its corresponding antibiotic category and absolute abundance, respectively. Figure S6.3: Principal coordinates analysis of resistome per sample with the relative abundance of ARGs at subtype level. Each dot represents a sample, colours indicate the concentration of added OTC and the point size indicates the sampling time. PCoA1 and PCoA2 together explain 87.82% of the total variation between samples. Figure S6.4: Important ARG subtypes for the changes of resistome per sample. R package ‘rfPermute’ was used to calculate the IncMSE value and significance, with settings of 500 decision trees and 1000 times random permutations. The identity of ARG subtypes was shown as y-axis labels, while the colour and size of the point indicate its corresponding antibiotic category and relative abundance, respectively. Table S6.1: Identified fast-growing and important ARG subtypes in this study. Bold text indicates ARGs that are present in all. Figure S6.5: Scheme of key ARG subtypes screening. A) 19 key ARG subtypes were screened out by cross-referencing as an intersection of the three lists (Tab. S6.1); B) distribution of 19 key ARG subtypes corresponding to antibiotic categories. Identities of 19 key ARG subtypes are shown in the main text Table 1. S7 Multitable co-inertia analysis for community aggregated traits. Figure S7.1. Stress plot representing the % of variation of the global dataset captured by each dimension of the MCOA. Table S7.1: Coordinates of samples on MCOA1 and MCOA2. Figure S7.2. linearly regression between the sum absolute abundance of ARGs per sample and its coordinates on MCOA1 and MCOA2. Figure S7.3: CATs contributions to the dimension MOCA1. Based on the RF model analysis, only CATs significantly (p < 0.05) contributed to the changes in coordinates reported. Bar colours indicate the LHS affiliation in line with the framework of YAS theory, and the direction indicates the positive or negative associations between the CATs and the MCOA dimensions. Figure S7.4: CATs contributions to the dimension MOCA2. Based on the RF model analysis, only CATs significantly (p < 0.05) contributed to the changes in coordinates reported. Bar colours indicate the LHS affiliation in line with the framework of YAS theory, and the direction indicates the positive or negative associations between the CATs and the MCOA dimensions. S8 Construction and modular analysis of positive correlation networks. Figure S8.1: Positive correlation networks and modules. The correlation networks contain 109 nodes and 1669 edges (i.e. significant correlations between key ARG subtypes and bacterial families, only positive correlations with Pearson coefficient > 0.8 and p < 0.01 are plotted). Modules are defined with the algorithm of Blondel et al. (2008) and five modules with more than two nodes are obtained (circled by colour dash). S9 Phylogeny analysis of assembled contigs that carrying ARGs. Table S9.1: Number of contigs carrying ARGs with phylogenetic information per sample. Figure S9.1: Phylogenetic information (at the family level) of contigs carrying ARGs in correspondence with their ARGs. From left to right, columns are OTC treatment conditions, family-level classification, and the types of ARGs carried.

40168_2024_2005_MOESM2_ESM.xlsx

Supplementary Material 2: Table S1: ASV table. Table S2: Shotgun metagenomic raw data cleaning. Table S3: The summary of assembled scaftigs. Table S4: Summary of metagenome-assembled genomes (MAGs). Table S5: Summary of predicted unigenes per sample. Table S6: Gene abundance per KEGG pathways. Table S7: Gene abundance per eggNOG catalogue. Table S8: Gene abundance per CAZy module. Table S9: Abundance per defined trait in this study. Table S10: Absolute abundance of ARGs (copies/ml). Table S11. Taxanomic classification of MAGs (GTDB taxanomy, Release 207). Table S12: Sequencing depth of MAGs.Tab. S12 sequencing depth of MAGs. Table S13: ARG detected on MAGs. Table S14: Mapped KO terms per each MAG. Table S15: Enhanced and diminished traits compared to CK. Table S16: Enrichment of KEGG pathways in MAG 1 and MAG 11. Table S17: Relative abudnance and 16S weighted abundance per bacterial family. Table S18: Correlation analysis between the CATs values of samples and the corresponding MCOA coordinate values. Table S19: Number of significant correlations between representative bacteria and key ARG subtypes. Table S20: Fast-growing ARG subtypes per each OTC treatment.

Supplementary Material 3: Table S1. The full list of CATs used in this study and their LHS affiliation.

40168_2024_2005_MOESM4_ESM.xlsx

Supplementary Material 4: Table S1 summary of ARGs-OAP analysis results. Table S2: Normalized abundance of ARGs (ppm). Table S3: Normalized abundance of ARGs to predicted 16S rRNA gene copies. Table S4: Normalized abundance of ARGs to predicted cell number.

Supplementary Material 5.

Supplementary Material 6.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. 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-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Z., Yao, X., Chen, C. et al. Growth of microbes in competitive lifestyles promotes increased ARGs in soil microbiota: insights based on genetic traits. Microbiome 13, 8 (2025). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40168-024-02005-6

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s40168-024-02005-6

Keywords

Microbiome

ISSN: 2049-2618

Contact us