Harnessing phage consortia to mitigate the soil antibiotic resistome by targeting keystone taxa Streptomyces

Abstract

Background

Antimicrobial resistance poses a substantial and growing threat to global health. While antibiotic resistance genes (ARGs) are tracked most closely in clinical settings, their spread remains poorly understood in non-clinical environments. Mitigating the spread of ARGs in non-clinical contexts such as soil could limit their enrichment in food webs.

Results

Multi-omics (involving metagenomics, metatranscriptomics, viromics, and metabolomics) and direct experimentation show that targeting keystone bacterial taxa by phages can limit ARG maintenance and dissemination in natural soil environments. Based on the metagenomic analysis, we first show that phages from activated sludge can regulate soil microbiome composition and function in terms of reducing ARG abundances and changing the bacterial community composition. This effect was mainly driven by a reduction in the abundance and activity of Streptomyces genus, which is well known for encoding both antibiotic resistance and synthesis genes. To validate the significance of this keystone species for the loss of ARGs, we enriched phage consortia specific to Streptomyces and tested their effect on ARG abundances on 48 soil samples collected across China. We observed a consistent reduction in ARG abundances across all soils, confirming that Streptomyces-enriched phages could predictably change the soil microbiome resistome and mitigate the prevalence of ARGs. This study highlights that phages can be used as ecosystem engineers to control the spread of antibiotic resistance in the environment.

Conclusion

Our study demonstrates that some bacterial keystone taxa are critical for ARG maintenance and dissemination in soil microbiomes, and opens new ecological avenues for microbiome modification and resistome control. This study advances our understanding of how metagenomics-informed phage consortia can be used to predictably regulate soil microbiome composition and functioning by targeting keystone bacterial taxa.

Video Abstract

Introduction

Dissemination of antimicrobial resistance in soils and farmlands poses a serious problem for the enrichment of antibiotic resistance genes (ARGs) in trophic networks and food production systems [1, 2]. Specifically, agricultural fields could be critical hotspots for the accumulation of ARGs due to the high input of both antibiotics and resistance genes along with manure and sludge fertilizers [3, 4]. Manure and sludge fertilizers have been recognized as major reservoirs of antibiotic residues and ARGs due to the excretion of antibiotics and resistant bacteria by humans and animals, contributing significantly to the accumulation of ARGs in agricultural soils [4]. While ARGs themselves do not pose health risks, they may cause serious problems when transmitted from commensal to pathogenic bacteria that cohabit with the same soil microbiomes [2]. This process is facilitated by mobile genetic elements (MGEs) that drive horizontal gene transfer in microbial communities [5], enriching genes that confer bacteria selective advantage in soils [6]. Although the mechanisms of ARG spread at the individual strain level are well understood, our predictive understanding of ARG movement within soil microbiomes at the community level remains limited [1]. Specifically, we lack information on how ecological interactions with non-bacterial microbial groups, such as bacteriophages (phages for short), could shape the prevalence and movement of ARGs in soil microbiomes [7].

Phages (i.e., specific bacterial viruses) are the most abundant and genetically diverse biological entities on Earth [8] and are crucial to regulating the diversity and composition of bacterial communities. Clinically, phage therapy has been recognized as a promising approach for targeting multi-drug-resistant bacteria [9], while phage cocktail therapy has demonstrated potential in preventing and managing certain soil-borne bacterial diseases affecting crops [10]. Despite their efficacy in lysing specific hosts, phages typically exhibit a narrow host range, limiting their stability and effectiveness in complex microbiomes [11, 12], such as soil, which harbors highly diverse microbiomes. Additionally, as MGEs, phages can facilitate the spread of ARGs [13]. Despite the attenuating effect of phage lysis on specific drug-resistant bacteria [14, 15], the broader impact of phages on the dynamics of antibiotic-resistant bacteria and their ARGs in complex soil microbiomes remains unclear [16]. Hence, elucidating the effects of phage consortia (i.e., multi-phage cocktails) on antibiotic-resistant bacteria and MGEs is crucial for advancing phage-based interventions in complex environmental microbiomes.

Metagenomic approaches provide powerful tools for identifying potential key species within complex microbiomes, while phage-mediated bacterial lysis enables a subtractive approach to discern their functional roles within these communities. Such metagenome-based phage biocontrol strategies offer a promising avenue for elucidating the roles of key bacterial species in the maintenance and spread of ARGs, as well as for developing ecological and environmentally friendly approaches to microbiome management and antibiotic resistance mitigation. Emerging evidence highlights the presence of keystone taxa within soil microbial communities—low-abundance species that disproportionately influence community composition and functioning [12]. For example, Nitrospira and Gemmatimonas have been identified as keystone taxa for soil microbiome stability due to their specialized roles in nitrogen and phosphorus metabolism [17]. However, the contribution of specific keystone bacterial species to the persistence and spread of ARGs within soil microbiomes remains poorly understood. Moreover, it is unclear whether phage-targeted lysis of such key bacterial taxa could effectively mitigate the spread of ARGs at the community level.

This study aimed to identify the keystone bacterial taxa responsible for maintaining and disseminating ARGs in agricultural soils and to explore the potential of phage-based strategies to mitigate ARG spread. Specifically, we used genome-resolved multi-omics approaches—including metagenomics, viromics, metatranscriptomics, and metabolomics—to investigate how phages isolated from activated sludge influence the structure and function of soil microbial communities, mobilome (all MGEs in a metagenome), and resistome (all ARGs in a metagenome). Municipal sewage sludge was selected because of its high phage diversity and abundance, allowing us to isolate phages that are probably associated with several bacterial pathogens and antibiotic-resistant bacteria. We conducted microcosm experiments using agricultural soils collected from different locations across China to assess the potential of Streptomyces-targeting phages in reducing ARG abundance and reshaping the soil resistome. This work provides a foundation for developing phage-based biotechnological tools to reduce ARG dissemination in complex soil ecosystems under the One Health framework.

Material and methods

Soil sampling and preparation

Soil samples (2.0 kg each) were aseptically collected from two representative cropland sites in southern China, characterized by distinct soil textures and intensive organic amendment practices. The sampling locations included Chongqing (CQ; purple soil, 107.1°N, 29.7°W) and Fujian (FJ; red soil, 119.3°N, 26.0°W), representing alkaline and acidic soils, respectively, that differed in their physicochemical properties (Fig. S1). Samples were taken from the 50 × 50 × 20 cm surface layer after crop harvest (maize), with all sites sampled simultaneously to minimize temporal variation in soil properties. The soils at both sites were collected in their natural state, avoiding anthropogenic disturbances. To focus on phage-bacteria interactions, bacterial and viral communities were separated to eliminate interference from other soil organisms, such as fungi and protozoa. The collected samples were immediately homogenized and sieved by 2-mm meshes after transfer to the lab in Fujian.

Extraction of soil bacterial communities

Indigenous soil bacteria were extracted from sieved fresh soil following a previously described protocol [18], with the following modifications. Approximately 100 g of sieved soil was added to 2 L sterile Erlenmeyer flasks containing 900 mL of phosphate-buffered saline (PBS, 0.01 M, pH 5.5), which maintains normal osmotic pressure and preserves microbial cell integrity. The suspension was homogenized by sonication (47 kHz, five 30 s cycles with 1 min intervals), followed by shaking at 30 °C and 180 rpm for 2 h. After settling at room temperature for 30 min, the suspension was centrifuged at 4500 g for 15 min at 4 °C. The supernatant, containing the extracted bacterial cells, was carefully collected for further analysis. Bacteria suspended in supernatants were purified and concentrated using a Cogent μScale Tangential Flow Filtration (TFF) system (Merck Millipore, MA, USA) with 0.2-μm membrane cassette (Pellicon XL, 50 cm²). This procedure was repeated three times and the retentates obtained from TFF were used as bacterial extracts (Fig. S1a). Finally, the collected soil bacteria were stored at 4 °C overnight prior to using them as inoculants in microcosm experiments the following day (see below).

Extraction and preparation of phages

Phages were extracted from activated sludge from sewage treatment plants following a previously described protocol, with modifications [19]. Activated sludge collected from Haixiang wastewater treatment plants (WWTPs) located in Fuzhou, China (119.3°N, 26.0°W). Activated sludge flocs result from microbial aggregation under aeration conditions, producing dense microbial communities embedded in extracellular polymeric substances. These flocs harbor high microbial diversity and are rich in phages, including those targeting human-associated bacteria and antibiotic-resistant strains [20,21,22]. Moreover, high phage diversity makes activated sludge a promising reservoir for isolating broad-spectrum phages with potential applications in controlling ARG dissemination in environmental settings. The steps for phage isolation were essentially the same as for bacteria, except that filtrates smaller than 0.22 µm were collected during ultrafiltration. Briefly, phage particles were detached from sludge flocs by sonication (47 kHz, five 30-s cycles with 1-min intervals) and subsequently purified using the same approach as for bacterial extraction, followed by tangential flow filtration (TFF). Sequential 0.2 μm and 100 kDa filters were employed to remove cellular debris and concentrate viral particles [18, 23]. To avoid possible bacterial contamination, the filtration process was repeated three times, and the extracted phages were verified through fluorescent microscopic imaging [23]. The final filtrated phage aliquot was confirmed to be free of living bacteria by incubating phage samples on LB plates for 48 h at 30 °C, resulting in no growth of visible colonies. The phage concentration was measured using fluorescent microscopic imaging and calculated as virus-like particles (VLP) per gram of soil [19]. Concentrated and purified phages (~ 10⁹ VLP/mL) were stored at 4 °C overnight prior to starting the microcosm experiments on the following day. Also, an aliquot of phage sample was autoclaved (121 °C for 20 min) to obtain inactivated phages for the control group treatment in the microcosm experiment (see below). The schematic diagram of bacterial and phage extraction and analytical processes is available in the Supplementary Information (Fig. S2–3).

Experiment 1: testing activated sludge phage community effects on soil microbiomes and resistome

To investigate the effect of sludge phages on the composition and resistome of bacterial communities isolated from both Chongqing and Fujian soils, we performed microcosm experiments using active living phages as the treatment group and inactivated phages as the control group. The homogenized soils were sterilized by autoclaving (121 °C for 12 h) three times with at least 24 h intervals between each sterilization cycle to remove living microbes and other organisms such as protozoa and fungi. The autoclaving was selected due to its accessibility and reproducibility for establishing a sterile microcosm system across multiple replicates. However, this method can alter soil chemistry, including the breakdown of humic and fulvic acids, which might limit extrapolating our findings to natural soils. Both Chongqing and Fujian soils were then distributed to microcosms (500 mL plastic bottles) where every replicate contained 200 g of sterilized soil and 25 mL of bacterial suspension derived from the original unsterilized soil. Subsequently, soil microcosms from Chongqing and Fujian were adjusted to 20% water content using sterile water, consistent with natural soil conditions, and sealed with breathable caps equipped with 0.22 μm filters. The microcosms were incubated at room temperature for 1 week after adding the bacterial suspension, allowing bacterial colonization into the soil. After 1 week, half of the microcosm replicates were spiked with 25 mL of either active or inactivated (control) phage suspensions. While the abundance of bacteria and ARGs in the extracted soil suspension was initially reduced compared to the natural soils, these levels rebounded within 7 days of incubation and stabilized at values comparable to the original soils (Fig. S4). Each of the four treatment groups was replicated 40 times and incubated at room temperature for 60 days (August to October 2023). Although microcosm experiments cannot fully capture the complexity of natural soil ecosystems, they offer a controlled experimental system for investigating the dynamic interactions between phages and bacteria. The four treatments were sampled destructively by randomly removing five replicate microcosms per treatment at the beginning (day 0) and then after 7, 14, 30, and 60 days from the start of the experiment (five sampling time points, resulting in a total of 100 microcosms). All microcosm samples were stored at − 80 °C to extract total DNA or RNA for microbial community analyses as described below.

Analyzing changes in bacterial microbiome, resistome, and metabolome composition

Samples collected at 0, 7, 14, 30, and 60 days were analyzed to investigate the succession of bacterial community composition using 16S rRNA amplicon sequencing. Total genomic DNA was extracted from all the soil sample timepoints using a Fast DNA spin kit (MP Biomedicals, Cleveland) according to the manufacturer’s instructions. Spectrophotometer NanoDrop ND- 2000 (Thermo Fisher Scientific, Wilmington) and gel electrophoresis were used to check DNA concentration and purity. All collected samples (each time point in each group consisting of five replicates) were subjected to 16S rRNA amplicon sequencing targeting the V4-V5 region using a NovaSeq6000 platform (Illumina, PE250 mode), and the raw 16S rRNA gene sequences were processed using QIIME 2 (v2019.7) [24]. The DADA2 pipeline in QIIME2 was adopted to produce amplicon sequence variants (ASVs) [25] followed by taxonomic classification based on the QIIME2 naive Bayes classifier trained on 99% operational taxonomic units available at the SILVA rRNA database (v 138) [26]. Microbial diversity was estimated using alpha diversity (Shannon and observed ASV richness) and community composition using beta diversity index (weighted uniFrac distance) based on the q2-diversity pipeline within QIIME2.

To track the dynamics of ARGs and MGEs during microcosm experiments, qPCR was used to quantify the abundance of 27 widespread ARGs and five common MGEs (Table S1) [27]. The 16S rRNA gene abundances were also quantified to normalize the abundances of ARGs and MGEs with bacterial cell abundances. The qPCR was carried out on a Light Cycler 96 system (Roche, Mannheim, Germany), and relevant details are provided in the Supplementary Information (Text S1). Based on the qPCR results, ARG abundances between phage and control treatments in FJ and CQ soils showed a stable difference after 30 days of incubation. Therefore, 30-day sampling time points were selected for metagenomic and metabolomic analyses (triplicate samples for each treatment).

Samples collected at 30 days were used to characterize the resistome and mobilome using metagenomic sequencing based on Illumina HiSeq X ten platform (2 × 150 paired reads) [28]. The raw sequences with average quality scores less than 30 or lengths smaller than 50 bp were removed using Trimmomatic (v0.39) [29]. The clean reads from metagenomic datasets were analyzed to identify and quantify the ARGs using local pipeline ARGs-OAP (v2.2) as previously described [30]. Briefly, ARG types and subtypes were classified using BLAST against the Structured ARG database (SARG) via the ARG-OAP pipeline. Resistance “Types” refer to the antibiotic classes to which ARGs confer resistance, while “Subtypes” represent individual ARGs. Metagenomic reads were aligned to the SARG database using BLASTX with an e-value cutoff of ≤ 10⁻⁷, 80% identity, and 75% coverage [31]. ARG abundance was calculated using the ARG-OAP pipeline, normalized to bacterial counts with units expressed as “copies of ARGs per cell.” Similarly, MGEs were identified and quantified using the same pipeline, substituting the SARG database with the MGEs database [31] using the same parameters.

Each soil metagenome was assembled individually using SPAdes (v3.13.1) with the parameters “-k 33, 55, 77, 99, 111, 127 –meta” for all samples [32]. Three replicates from one site in the same treatment were merged in the first microcosm experiment and further co-assembled to obtain more contigs using the same parameters. All assembled contigs longer than 2.0 kb in each metagenome were binned individually using Metawrap [33] based on MetaBAT2 [34], MaxBin2 [35], and Concoct [36] using the “Assembly module” with default parameters. Bins were further curated to obtain high-quality genomes using “Bin_refinement module” in Metawrap [33]. All bins from each sample were merged together and were dereplicated using dRep v 2.3.2 [37]. The quality of metagenome-assembled genomes (MAGs) was assessed by CheckM (v1.0.13) [38] and MAGs with higher than 50% completeness and less than 10% contamination were retained for further analyses. The taxonomic classifications of MAGs were based on the Genome Taxonomy Database (GTDB; release 03-RS86) using the GTDB-Tk toolkit (v 0.3.2) with classify workflow [39]. For the ARGs in MAGs, we further assessed their health risk using the bioinformatics tool arg_ranker (https://github.com/caozhichongchong/arg_ranker) based on metagenomic contigs [40]. The health risk of ARGs in MAGs was based on the following three criteria: (1) enrichment in human-associated habitats, (2) gene mobility, and (3) ARG presence in pathogens according to previous studies [40]. ARGs that did not meet the first criterion were assigned to Rank IV (lowest risk). Those that met the first criterion but not the second were assigned to rank III. ARGs that met the first and second criteria, but not the third, were assigned to rank II, and those that met all three criteria were assigned to rank I (highest risk). The enrichment of ARGs in human-associated environments was based on our previous results [41], where we measured the changes in their relative abundance in human-associated habitats (feces of humans and farmed animals) and natural environments (soil, sediment, and seawater), while ARG mobility was determined by the presence of a mobile genetic element associated with the ARGs in same contig. A deep machine learning network approach of DCiPatho was used to identify the potential pathogenic bacteria [42]. To further study the virulence factor of the pathogens, the predicted genes of each MAG were annotated for the virulence factor genes (VFGs) using BLASTP (e-value of 1.0 × 10⁻⁵, > 80% sequence identity) against VFGs databases (https://www.mgc.ac.cn/VFs/download.htm, downloaded on January 2025) [43]. Overall, ARGs ranked as I and II were assigned to high-risk group, while ARGs ranked as III and IV were considered as low-risk group.

The relative abundances of MAGs were quantified using the CoverM pipeline [44] (v0.61, https://github.com/wwood/CoverM) based on the coverage of mapped reads using “genome” mode. Briefly, reads were first mapped to MAGs using “make” command to create BAM files (–percentage_id 0.95 –percentage_aln 0.75). Filtered bam files were then used to generate coverage profiles across samples (–trim-min 0.10 –trim-max 0.90 –min-read-percent-identity 0.95 –min-read-aligned-percent 0.75 -m mean). The coverage of each MAG in the unit of “rpkm” (reads per kilobase of exon per million reads mapped) was directly used for alpha- and beta-diversity analyses (log-transformed matrices were used for Mantel correlations).

Biosynthetic gene clusters (BGCs) were predicted for curated MAGs by using antiSMASH (v6.0) [45]. BiSCAPE was further used to profile the BGCs and cluster them into Big-SCAPE families. The class and number of BGCs were extracted from the antiSMASH output by manual curation, and the number and diversity of predicted BGCs of MAGs were investigated (all the BGC hits within MAGs are listed in Table S2). Moreover, non-targeted metabolomics was conducted for all day 30 soil samples (n = 4 for each treatment group) to understand variation in soil metabolite composition in association with microbial community change [46]. The soil metabolites were extracted by standard procedure [46] and analyzed using one gas chromatograph (Agilent 7890, Santa Clara, USA) coupled with a Pegasus HT time-of-flight mass spectrometer (GC-TOF–MS, LECO, St. Joseph, MI) [46]. As the genus Streptomyces is well recognized for the production of antimicrobials in soils, we specifically tracked the production of four Streptomyces-produced antibiotics (streptomycin, vancomycin, tetracycline, and cycloheximide) using solid phase extraction targeted approach coupled with high-performance liquid chromatography (Agilent 7890B, Agilent, USA) and a tandem mass spectrometer (UPLC ACQUITY, Xevo TQ-S, Waters, USA) [47]. A more detailed description of the analytical methods of antibiotics is included in Supplementary materials (Text S2).

Phage DNA extraction, sequencing, and soil virome analysis

Soil phage particles were extracted and concentrated from the samples through centrifugation and filtration [48]. Briefly, 20 g of soil samples collected on day 30 were mixed with 900 mL of PBS and shaken at 4 °C for 15 min (120 rpm). The soil wash was centrifugated at 5000 × g for 10 min at 4 °C after the supernatant was filtered through 0.22 μm membrane using the tangential flow filtration (TFF) system to separate phages. The filtrates were sequentially concentrated three times using a TFF system with a 100-kDa membrane (Millipore, American). The phage concentrate was treated with 20U DNase I (37 °C, 50 min, Transgen Biotech) to remove naked DNA using a TIANamp Virus DNA/RNA Kit (TIANGEN, Beijing, China). Sequencing libraries were prepared using the DNA Library Prep Kit V2, and then, metagenomic sequencing was performed on the Illumina Novaseq 6000 platform (Illumina, California, USA) to obtain 150 bp paired-end reads. The quality control of raw reads was conducted by Trimmomatic (v0.39, score > 30 and length > 36 bases) [29], and the quality-filtered sequences were individually assembled using Spades using the same parameters as with bacterial metagenomes.

To identify viruses with high confidence, phage contigs larger than 5 kb were recovered from metagenome assemblies and analyzed using a combination of three tools with strict quality thresholds as previously suggested [49]. First, DeepVirFinder (v1.0) [50] was run with loose cutoff values (score 0.7 and p < 0.05) to detect potential phage sequences. Second, these DeepVirFinder-output sequences were used as input files with VirSorter2 (v2.2.1) [51] to identify putative phage sequences with scores ≥ 0.95. The CheckV (v0.9.0) [52] was used for quality assessment to remove potential non-phage sequences after the VirSorter2 analysis. The final phage contig dataset was manually curated and trimmed to remove potential host regions following our previous protocol [53]. The remaining contigs were considered to be phage origin if they met at least one of the three following criteria: (1) contigs contained at least one virus-specific hallmark gene, (2) contigs had VirSorter2 scores ≥ 0.95, and (3) contigs had DeepVirFinder score ≥ 0.9 and p-value < 0.05. The final phage contigs were clustered based on 95% average nucleotide identity and at least 85% coverage using CD-HIT v4.8.1 [54] (parameters: -c 0.95 -aS 0.85), resulting in final virus OTUs (from here on referred as “vOTUs”). The relative vOTU abundances in metagenome datasets were quantified using the CoverM pipeline [44] based on the coverage of mapped reads using the “contig” mode and “rpkm (reads per kilobase of exon per million reads mapped)” calculation method.

Taxonomic assignment of phage contigs was performed using PhaGCN2 and the latest ICTV classification tables [55] against the RefSeq viral database (v216, released in Jan. 2023). The phage lifestyle was predicted using three tools including VIBRANT [56], PhaTYP [57], and manually curated BLAST [58] based on previously described methods [58, 59]. Briefly, VIBRANT [56] and manual BLAST [58] were used to infer temperate lifestyle by identifying viral contigs that contained proteins associated with lysogeny (transposases, integrases, excisionases, resolvases, and recombinases). In case of incompleteness of phage contigs, a machine learning method called PhaTYP [57] was used to predict the phage lifestyle. A virus was inferred as a temperate phage if it was predicted to be lysogenic by any one of these tools. The non-redundant proteins were annotated using Pfam (v32), KEGG (release Sep. 2020) and eggNOG databases (v5.0) with a cutoff of e-value < 0.0001 [19, 44].

To predict associations between viral and bacterial sequences, vOTUs were linked to bacterial MAGs using three computational methods [44, 60]: (1) shared genomic content between viral contigs and bacterial MAGs (bitscore ≥ 50, e-value < 10⁻³, identity ≥ 70%, and matching length ≥ 2500 bp [44, 61]); (2) CRISPR spacers targeting [59]; and (3) matching of vOTU-derived tRNA sequences with bacterial tRNAs (95% coverage and 90% sequence identity, e-value < 1e − 05) [53]. The CRISPR spacers were recovered from bacterial MAGs using CRT (v1.2) with default parameters [62]. Extracted spacer sequences were matched with vOTUs using BLASTn (100% nucleotide identity, 90% coverage, mismatch ≤ 1, and e-value ≤ 10⁻⁵). The tRNA sequences were predicted using tRNAscan (v2.0.9), where the general and bacterial/archaeal models were run sequentially and compared with BLAST after deleting self-hits and duplicates [61]. In case of Streptomyces-enriched phages in the second confirmatory microcosm tests below, we used an integrated machine-learning tool to identify viral hosts based on tools databases and recovered MAGs [63].

Soil RNA extraction, RNA-sequencing, and metatranscriptomics analysis

Soil samples from both microcosm experiments were collected on day 30 and immediately preserved in liquid nitrogen for metatranscriptomic analysis. The total RNA was extracted from soil samples using the RNeasy PowerSoil total RNA kit (Qiagen) following manufacturers’ protocols. RNA concentrations were measured using the Qubit RNA HS assay kit, and RNA integrity was determined using an Agilent 2100 Bioanalyser (Agilent Technologies) [64]. The rRNA was removed with Ribo-minus Transcriptome Isolation Kit (Thermo Fisher, Waltham, USA), and the purified mRNA was prepared for sequencing using the TruSeq stranded mRNA library prep kit (Illumina, California, USA) following the manufacturer’s instructions. The extracted DNA and cDNA were used to construct libraries (300 bp) by DNA Library Prep Kit V2 (Illumina, California, USA). The sequencing was conducted on the NovaSeq6000 platform in PE150 mode, and a more detailed description of RNA extraction is included in Supplementary materials (Text S3).

To analyze the metatranscriptomics data, reads were quality filtered with Trimmomatic (v 0.39, score > 30 and length > 50 bases) [44] followed by the removal of non-coding RNA sequences (tRNA, tmRNA, 5S, 16S, 18S, 23S, and 28S rRNA sequences) using SortMeRNA (v4.3.4) [65]. The remaining total mRNA reads (three replicates per treatment) were mapped back to MAGs and vOTUs to identify active bacterial and viral taxa in terms of average coverage of mapped transcripts using minimap2 [66]. Briefly, metatranscriptomic sequence reads were used as an input with the same parameters as with the metagenomic read mapping, except for choosing “tpm” as the calculation method (transcripts per kilobase per million mapped reads, tpm). All activity was quantified at the level of vOTUs and MAGs, and the genes of viruses and bacteria were deemed as active when the “tpm” values were larger than 0 [53]. To evaluate the transcription of putative BGCs within a MAG, we searched for any transcripts mapped to the same genomic region of predicted biosynthetic genes (identity 95% and 75% coverage).

Experiment 2: Regulating soil resistome using Streptomyces-enriched activated sludge phage community

To verify the influence of Streptomyces-specific phages on ARG abundances, we used Streptomyces-enriched phages to repeat the first soil microcosm experiment using a wider range of farmland soils collected across 16 provinces in Southeast China (16 sites, 3 independent replicate samples per field site, totaling 48 soil samples, details in Table S3). Given the high heterogeneity of soils, variations in soil properties and texture were considered as potential factors influencing microbial community composition. To test the stability and efficacy of phage-induced lysis, additional soil samples were collected. These sampling sites covered the major cropland soil areas with different soil types and textures in Southeast China (Fig. S2 d). To isolate Streptomyces-specific phages, we first isolated eight Streptomyces strains from the CQ and FJ soils using selective glucose yeast malt extract (GYM) media supplemented with 50 μg/mL cycloheximide and nalidixic acid [67]. The identity of Streptomyces isolates was identified based on the 16S rRNA gene (Table S4). Before enrichment, all eight Streptomyces strains were individually cultured in a GYM liquid medium at 30 °C for several days to obtain spore suspensions (1 × 10⁸ CFU/ml for each strain). All strains were then mixed in equal volumes (10 mL each) in 500 ml of liquid GYM medium with the same activated sludge phage extract, which was derived from the first experiment. The enrichment was finished after 1 week of bacteria-phage coculturing at 30 °C, and Streptomyces phages were purified and concentrated using the TFF systems with 100 kDa filters as described previously. The final phage stock (~ 10¹⁰ VLP/mL) was stored at 4 °C for the second microcosm experiment.

To test the effect of Streptomyces-specific phages on the composition of soil resistome and mobilome, all soils were processed and characterized the same way as with CQ and FJ soils used in the first experiment. The soil microcosm experiments were set up the same way as previously except that we used Streptomyces-enriched phages as the active or inactive phage inoculum and that the experiment was run only for 30 days given that the effect of phage treatment on soil resistance stabilized after 30 days of incubation in the first microcosm experiment. Three replicates of each soil under treatment and control were conducted during the microcosm incubation experiments (16 sites * 3 replicates * 2 treatments, a total of 96 samples). The microcosms were sampled only once after 30 days, and changes in microbial community composition (16S rRNA amplicon sequencing and metagenomics), resistome, and mobilome were analyzed as described above. To explore the effect of phage treatments on microbial activity, we randomly selected all replicate soil samples from two sites (Nanchang and Guangzhou) for metatranscriptomics analysis. More detailed information on all soil sites and DNA and RNA sequencing datasets are described in Supplementary materials (Table S5).

Statistical analysis

Statistical analysis was conducted using the R-studio platform (v 4.2.2, https://www.r-project.org/) [68]. The analyses of microbial alpha and beta diversity (including alpha index and PCoA) were conducted using vegan and ggplot2 packages in R. The overall mean differences across treatments and sampling time were analyzed using Student’s t-test, repeated measures ANOVA, and one-way ANOVA followed by multiple comparisons using the Tukey HSD test (p < 0.05). Non-parametric PERMANOVA (Adonis function, 999 permutations) was used to determine the significance of different treatments on the microbiome composition and functional gene differences. Normalized transcript counts per gene cluster per time point were compared between groups (phage treatment and control group) using Student’s t-test.

A partial Mantel test was used to assess the correlation between two multivariate matrices while controlling for the potential effects of resistome using the R package “vegan.” Distances for bacterial community and resistome and mobilome abundances during microcosm experiments were calculated using the Bray–Curtis dissimilarity matrices. Partial Mantel correlations were computed between all pairs of distance matrices for bacterial community and resistome with 999 permutations for each comparison. The false discovery rate was computed using the Benjamini–Hochberg method.

To identify the major taxonomic predictors explaining the resistome dynamics during phage treatment in the second experiment, we compared the contribution of bacterial and phage community composition (based on Bray–Curtis dissimilarities) using default parameters of the Random Forest (RF) analysis [69]. In these RF models, the importance of each microbial predictor (bacterial and phage community composition based on metagenomic data) was measured for the resistome and mobilome dynamics (Euclidean distance dissimilarity). To assess the relative importance of different predictor variables, we compared the percentage increases in MSE (mean squared error), where high MSE percentage values denote for relatively higher contribution by given predictor variables [70]. The significance of the models and cross-validated R² values were assessed with 1000 permutations of the response variable using the “A3” package in R. The analysis was conducted using the “rfPermute” package in R.

Results

Phages isolated from the activated sludge drive changes in soil microbiome composition

We first tested if phage communities isolated from the activated sludge affected the microbiome composition of two agricultural soils collected from Chongqing (CQ) and Fujian (FJ). Based on qPCR, bacterial abundances in both control soil groups inoculated with inactivated phages increased continuously during the 30 days of incubation (Fig. 1a), indicative of successful colonization of soil bacteria. The addition of active phages completely changed this pattern, resulting in clear reductions in bacterial abundances in both soils, indicative of tight top-down control by phages (F_4,45 = 5.5, p = 0.0011 in FJ and F_4,45 = 4.9, p = 0.0022 in CQ; Fig. 1a). In addition to reducing bacterial abundances, phages clearly changed the bacterial community composition (R² = 0.56 and p < 0.001 in FJ; R² = 0.46 and p = 0.001 in CQ, PERMANOVA test) and reduced bacterial richness (p = 0.008, t = 3.5, df = 8 in FJ and p < 0.001, t = 8, df = 8 in CQ) in both soils (16S rRNA amplicon sequencing; Fig. 1b and Fig. S5). Actinobacteriota (now named Actinomycetota) and Firmicutes (now named Bacillota) were the most dominant bacterial taxa in the soil community, accounting for more than 90% of the relative abundance of bacteria. Specifically, phages reduced the relative abundance of Actinobacteriota (dominated by Streptomyces genera) by 13.5% and 66.2% in CQ and FJ soils, respectively (Fig. 1c and Fig. S6), and concomitantly increased the relative abundance of Firmicutes (dominated with Paenibacillus and Alicyclobacillus genera) by 16.4% and 93.1% in CQ and FJ soils.

Fig. 1

The effect of phages on bacterial and viral abundances, diversity, and community composition during the first soil microcosm experiment in CQ and FJ soils. a–d The bacterial (a) and phage (d) abundance between phage (cyan) and control (red) groups in CQ and FJ soils during the microcosm experiment. Bacterial density is based on 16S rRNA gene copy numbers per gram of dry soil (n = 5 biologically independent samples per treatment at each time point). Virus-like particles (VLPs) were counted by epifluorescence microscopy (n = 5 biologically independent samples per treatment at each time point). b Comparison of bacterial community composition between phage and control groups based on distance-based PCoA analysis at a 30-day sampling time point (n = 5 biological replicates per treatment). c The mean relative bacterial abundances at the phyla level between phage and control groups at a 30-day sampling time point. e Changes in viral community composition based on distance-based (Bray–Curtis) PCoA analysis between phage and control groups at a 30-day sampling time point (n = 3 biological replicates per treatment). f The mean relative viral abundances at the family level between phage and control groups at a 30-day sampling time point. The data in panels (b–c) and (e–f) are based on amplicon sequencing and viral metagenomic sequencing, respectively

Overall, phage abundances were higher in soils inoculated with active phages, peaking at day 30 after phage inoculation (p < 0.001, t = − 9.5, df = 8 for CQ and p < 0.001, t = − 6.0, df = 8, Student’s t-test, Fig. 1d). After phage sequence clustering and dereplication (95% identification and over 85% coverage), a total of 2608 viral operational taxonomic units (vOTUs, > 5 kb) were identified (Table S6) that were mainly lytic dsDNA phages (60.0%). Inoculation of active phages had significant effects on phage community composition in both soils (R² = 0.56, p < 0.001 for CQ and R² = 0.36, p = 0.001 for FJ, PERMANOVA test, Fig. 1e) and increased phage richness in the FJ soil (p = 0.03, t = − 2.9, df = 4). While Autographiviridae, Chaseviridae, and Mesyanzhinovviridae were the dominant phage families in both control groups, the relative abundance of the Chaseviridae family significantly increased from 1.9 to 9.7% in CQ (p = 0.003, t = − 5.3, df = 4) and from 43.8% to 51.6% in FJ (p = 0.002, t = 5.6, df = 4) soils inoculated with active phages (Fig. 1f). Together, these results show that phages from the activated sludge imposed strong top-down regulation on bacteria, which was coupled with changes in the composition of both bacterial and phage communities (Fig. S7) likely due to disproportionate lysis of different bacterial taxa.

Phage application decreases the abundance of high risk of ARGs and MGEs

We next used metagenomics to compare the phage effects on the ARG and MGE abundances 30 days after the phage inoculation. The total abundance of ARGs was reduced in both CQ (p < 0.001, t = 17.6, df = 4) and FJ (p = 0.01, t = 4.2, df = 4) soils compared to control groups (Fig. 2a), even though this effect varied between sampling time points (based on qPCR, Fig. S8). Across both control and phage treatment groups, 12 ARG types and 152 ARG subtypes were identified, which included beta-lactam (29.8%), multidrug (29.1%), tetracycline (7.9%), vancomycin (6.6%), macrolide-lincosamide-streptogramin (MLS, 6.6%), and aminoglycoside (5.2%) ARGs (Fig. 2b). Phage treatment decreased the ARG diversity in CQ (p < 0.001, t = 10.7, df = 4) and FJ (p = 0.01, t = 4.01, df = 4) soils (Fig. 2c) and specifically reduced the aminoglycoside, MLS, tetracycline, and vancomycin ARG abundances. Although a few ARG types such as bacitracin and rifamycin experienced some increase in the phage-treated groups, the total abundance of ARGs decreased by 71.9% in the phage-treated groups across both sites (Fig. 2a,b). The enrichment of these ARG subtypes in phage treatments is not clear but could have been due to the enrichment of certain ARG-carrying bacterial taxa, such as Paenibacillus and Alicyclobacillus genera. Overall, the most prominent positive correlations were observed between vancomycin resistance genes and those conferring resistance to beta-lactams, fosfomycin, multidrug agents, and MLS antibiotics. Except for bacitracin, all categories of antibiotic-resistance genes showed significant correlations with MGEs, confirming the biological functions of MGEs in ARG dissemination (Pearson’s correlation coefficient, all p < 0.05; Fig. 2d).

Fig. 2

The effect of phages on the abundance, composition, and diversity of ARGs and MGEs during the first soil microcosm experiment in CQ and FJ soils. a Changes in the abundance of ARGs and MGEs in phage and control groups at a 30-day sampling time point in CQ and FJ soils (n = 5 biological replicates per treatment). d Heatmaps representing fold changes in the abundance of individual ARG and MGE genes at a 30-day sampling time point in the phage relative to the control group. c Changes in the richness of ARGs and MGEs in phage and control groups at a 30-day sampling time point in CQ and FJ soils (n = 5 biological replicates per treatment). d Overall correlations between different ARG and MGE type abundances at a 30-day sampling time points in CQ and FJ soils (n = 5 biological replicates per treatment). Edge widths correspond to the absolute value of the correlation coefficient determined by Mantel’s statistics. Colors indicate correlation types. Solid and dashed lines denote significant and non-significant correlations, respectively. Pairwise comparisons of environmental factors are shown in the matrix with a color gradient denoting Pearson’s correlation coefficient. e The overall correlation between ARG and MGE abundances at a 30-day sampling time point. In panels a and c, box plots encompass 25 th–75 th percentiles, whiskers show the minimum and maximum values, and the midline shows the median (dots present the biologically independent samples, and asterisks denote significant differences (*p < 0.05, **p < 0.01, and **p < 0.001)). The shaded area in panel e represents a 95% confidence interval around the fitted regression line

Similar to ARGs, the MGE abundances (p = 0.012, t = 4.2, df = 4 in CQ and p = 0.0062, t = 5.2, df = 4 in FJ) and diversity (p = 0.002, t = 6.9, df = 4 in CQ and p = 0.038, t = 3.3, df = 4 in FJ) were reduced by active phages in both soils (Fig. 2a–c). Specifically, the total abundances of three dominant types of MGEs (transposase, integrase, and IS91) were decreased in the presence of active phages in FJ (p = 0.041, t = 3.9, df = 4) and CQ soils (p = 0.03, t = 3.6, df = 4, Fig. 2b), with IS91 and transposase abundances being reduced the most (Fig. 2b). Variation in transposase and integrase abundances were positively associated with more than 8 and 5 ARG type abundances, respectively (Mantel’s statistics, all p < 0.05, Fig. 2d). Moreover, the total MGE and ARG abundances were positively correlated (R² = 0.54, p = 0.0039, Fig. 2e), highlighting the importance of horizontal gene transfer for the ARG prevalence. Further, Procrustes analysis showed that changes in bacterial community composition were positively correlated with changes in ARG (M² = 0.71, p = 0.007) and MGE (M² = 0.55, p = 0.012) composition during the microcosm experiments (Fig. S9). Together, these results suggest that phage application changed the soil microbiome resistome and mobilome through effects on bacterial community composition.

Streptomyces is a keystone bacterial taxon that encodes antibiotic synthesis and resistance genes

Metagenomic assembled contigs were binned to gain insights into which bacteria were the most likely carriers of ARGs and how their abundances were affected by the phage application. By using all metagenomic samples in both microcosm experiments, a total of 327 bacterial metagenome-assembled genomes (MAGs) were obtained with medium to high quality (with estimated completeness of ≥ 50% and contamination ≤ 5%). Dereplication reduced the total number of bacterial MAGs to 208 (across the 10 phyla), which varied in their completeness between 50.6% and 98.0% (Fig. 3a and Table S7). The taxonomic composition of MAGs at the phylum level was consistent with the taxonomic composition based on 16S rRNA sequencing (Fig. S10), with Actinobacteriota (45.5% of MAGs), Firmicutes (33.2%), and Proteobacteria (6.2%, now named Pseudomonadota) being the dominant phyla. Based on deep machine learning for pathogen identification [42], 16.3% of ARG-carrying MAGs (n = 34) were identified as potential pathogenic bacteria infecting humans, animals, and plants, and most of them (61.5%) carried virulence-associated genes (Table S8).

Fig. 3

Genomic analysis of ARGs, MGEs, and biosynthetic gene clusters (BGCs) present in 208 MAGs extracted from the first soil microcosm experiment soil samples. a A phylogenetic tree of 208 bacterial MAGs was recovered from metagenomes of all samples based on a concatenated set of 120 conserved bacterial single-copy marker genes. Inner circles are colored by phylum-level taxonomy (see source legend on the left). The different circle colors from inner to outer rings represent the number of ARGs, MGEs, and BGCs in MAGs. The tree scale shows nucleotide substitutions for each site. b The genetic structure of biosynthetic gene clusters (BGC) is associated with known antibiotic production in three Streptomyces MAGs. c Zi-Pi plot exhibiting the distribution of MAGs based on topological roles in a constructed co-occurrence network. Each node represents an MAG and the position of nodes is determined by the scatter plot of within-module connectivity (Zi) and among-module connectivity (Pi). Two Streptomyces taxa were identified as the network hubs among the community in the top right corner. Module hubs were defined with Zi > 2.5 and Pi < 0.62, network hubs were defined with Pi > 0.62 and Zi > 2.5, and connectors were defined with Pi > 0.62 and Zi < 2.5

A total of 495 ARGs were detected across 188 MAGs, and 90.3% (188 of 208 MAGs) could be considered antibiotic-resistant bacteria based on carrying at least one ARG (Table S9). MAGs that carried more than three ARGs mainly belonged to Actinobacteriota, Chloroflexota, Proteobacteria, and Myxococcota phyla (Fig. S11), and at the genus level, Streptomyces (bin_182, bin_183, and bin_184), Protaetiibacter (bin_165), and Spirillospora (bin_180) from the Actinobacteriota phylum could be considered as multidrug-resistant bacterial hosts with eight ARGs per genome on average (Fig. 3a). All MAGs contained at least one MGE (a total of 2965 MGEs across 208 MAGs: 14.2 MGEs per MAG on average, Table S10), which mainly consisted of transposases (73.2%), integrases (21.8%), recombinase (4.8%), and plasmids (0.2%). The MAGs belonging to phyla Actinobacteriota, Firmicutes, and Proteobacteria carried at least two MGEs, and specifically, bin_182 and bin_184 corresponding to Streptomyces, carried a high number of MGEs: 35 and 16 MGEs per MAG, respectively (Fig. 3a). In particular, three Streptomyces MAGs were identified as potential pathogens carrying a large number of mobile ARGs and several virulence-associated genes. Based on the ARG risk ranking (see methods), among the 495 ARGs in the MAGs, the majority of ARGs (91%) were considered low-risk (ranks III and IV), while high-risk ARGs (ranks I and II) only accounted for 9%. Of all human-associated and mobile ARGs, 71% were associated with potential pathogens within the rank I group (Table S9). The presence of ARGs in Streptomyces MAGs (bin_182, bin_183, and bin_184) was also associated with the presence of 20 biosynthetic gene clusters (BGCs) known to encode several antimicrobial compounds, including streptothricin, kanamycin, tetracycline, ficellomycin, mithramycin, and cyphomycin amongst others (Fig. 3b and Table S11). Streptomyces taxon was hence associated with both resistance and biosynthesis of antibiotics.

To explore if the reduction in ARG and MGE abundances by active phages was linked to a reduction in Streptomyces taxon, we explored phage effects on the MAG abundances and composition. As expected, active phages changed the composition of MAG communities in both soils (Fig. S12) and specifically reduced the abundance of Actinobacteriota dominated by the Streptomyces taxon by 50% in both CQ (p = 0.0025, t = 10.1, df = 4) and FJ (p = 0.014, t = 5.5, df = 4) soils (Fig. S13). The relative abundance of the three Streptomyces MAGs reached 14.4% in CQ soil and 18.7% in FJ soil, making them the most dominant taxa in the overall bacterial community. Following phage treatment, their abundance significantly decreased by 1.2-fold in CQ soil and 7.9-fold in FJ soil (p = 0.006, t = 8.8, df = 4 for CQ; p < 0.001, t = 23.4, df = 4 for FJ, Fig. S14). Based on network analysis, Streptomyces genus abundances were positively correlated with 56 genera across 12 phyla (SparCC p < 0.05), indicative of its success in bacterial communities (Fig. S15). Notably, two Streptomyces MAGs were identified as the network hubs (nodes being both a module hub and a connector), indicative of their potential key role in bacterial communities (Fig. 3c). The Random forest model further corroborated that Streptomyces taxa explained the relatively highest proportion of variance in ARG abundances relative to other bacterial genera (Fig. S16). Together, these results suggest that Streptomyces was a likely key taxon for ARG and MGE carriage and antibiotic synthesis whose abundances were significantly reduced by the phage application.

Phages isolated from the activated sludge can control the abundance of Streptomyces-associated ARGs, MGEs, and antibiotic synthesis genes

To identify phages potentially capable of infecting Streptomyces bacterial hosts, computational approaches (see the “Methods” section) were used to link bacterial MAGs with vOTUs (Table S12). Around 4.7% of soil vOTUs (123 of 2608 viral contigs) could be linked with 108 bacterial MAGs as putative hosts, which comprised the relatively most abundant bacterial taxa: Actinobacteriota (50.3% of virus-host pairs) and Firmicutes (26.2% of virus-host pairs, Fig. 4a). Among the 123 vOTUs with putative hosts, 42 vOTUs were predicted to infect Streptomyces (Table S12). These Streptomyces phages consisted of 17 virulent and 25 temperate phages, which represented 3.4% of the total phage abundances (Fig. S17). Phage application significantly increased the proportion of Streptomyces phages by 2.8- and 3.4-fold in both CQ and FJ soils (p < 0.001, t = 33.1, df = 4 for CQ and p < 0.001, t = − 39.6, df = 4 for FJ, Fig. 4b). Moreover, the mean relative Streptomyces MAG abundances decreased from 91.4 to 8.6% in phage treated soils (Fig. 4b, p < 0.001, t = 54.9, df = 4 in CQ and p < 0.001, t = − 10.7, df = 4 in FJ), which was also coupled with reduced transcriptional activity of Streptomyces MAGs (p < 0.001, t = − 10.3, df = 4 in CQ and p < 0.001, t = − 13.7, df = 4 in FJ; i.e., 7.5 fold in CQ and 2.8 fold reduction in FJ, Fig. 4c). Increase in Streptomyces-specific phage abundances correlated negatively with ARG (R² = 0.87, p < 0.001) and MGE abundances (R² = 0.59, p = 0.0019, Fig. S18), suggesting that selective infection of Streptomyces hosts limited the prevalence and horizontal transfer of ARGs.

Fig. 4

Phages are effective at controlling Streptomyces genus abundances, activity, and metabolism, including antibiotic synthesis in the soil. a Predicted virus-host links between vOTUs and bacterial MAGs. The left two panels represent bacterial host taxonomy colored by phylum and genus, and the right two panels show virus families and associated sampling treatments. Gray connecting lines show associations between bacterial hosts (at phylum and family level) and their viruses in the given sampling treatment (shown on the right). b Changes in the relative abundances of Streptomyces MAGs and their phages in phage and control groups at a 30-day sampling time point in CQ and FJ soils. c Changes in the relative transcriptional activity of Streptomyces MAGs and their phages in phage and control groups at a 30-day sampling time point in CQ and FJ soils. d Principal component analysis showing differences in soil metabolomes between phage and control groups in CQ and FJ soils (n = 4 biological replicates per treatment). e Changes in the production of four antibiotics known to be produced by Streptomyces taxa between phage and control groups in CQ and FJ soils. In panels b, c, and e, data show mean ± SD of triplicates per treatment. The asterisks indicate significant differences between phage and control groups based on the Student’s test (*p < 0.05, **p < 0.01, and ***p < 0.001)

To explore if the reduction in ARGs also led to a reduction in the synthesis of antibiotics, we first compared the biosynthetic gene clusters (BGCs) of all MAGs using antiSMASH (v6.0). A total of 940 BGCs were identified across 181 bacterial MAGs (around 5.1 BGCs per MAG, Table S2) that were predicted to be associated with the synthesis of nonribosomal peptides (NRPS; 21.2%), polyketides (PKS; 16.1%), terpenes (17.6%), saccharides, and ribosomally synthesized and post-translationally modified peptides (RiPP; 23.9%, Fig. S19a). Most of the BGCs (44.1%; 415 out of 941) were found within Actinobacteriota phylum (Fig. S19b) that carried more than seven BGCs per MAG, while other phyla contained an average of four BGCs per MAG. Specifically, Streptomyces carried the greatest number of BGCs (mean 24 per genome, Fig. 3a), indicative of their importance in the production of secondary metabolites.

Metabolomic analysis was further employed to link the impact of active phages on soil antimicrobial metabolome at a 30-day sampling time point. PCoA analysis revealed that phages clearly altered the soil metabolome compared to the control group in both CQ and FJ soils (R² = 0.76, p = 0.001 for CQ and R² = 0.54, p = 0.024 for FJ, PERMANOVA test, Fig. 4d and Fig. S20). Furthermore, changes in metabolomes were positively correlated with changes in bacterial community composition (Mantel statistic r = 0.64, p < 0.001; based on beta-dissimilarity of MAGs), suggesting that phage infections drove these changes. As expected, phages were effective at reducing the production of Streptomyces-associated antibiotics, including streptothricin, kanamycin, and vancomycin in both CQ and FJ soils (61.2% and 37.3% mean reductions in CQ and FJ, respectively, Fig. 4e). Together, these results show that in addition to reducing ARG abundances, phages also decreased the synthesis of antibiotic in both soils, which could have further lowered the selective benefit of ARG carriage.

Streptomyces-enriched phages can be used to reduce the ARGs and MGEs in agricultural soils

To confirm the importance of Streptomyces and Streptomyces-specific phages for ARG dissemination in agricultural soils, we enriched Streptomyces phages from the original sludge phage community, and directly tested their effect on resistomes of agricultural soil samples collected from 16 locations across China using a microcosm experiment (Fig. S2 d). The phages enriched using Streptomyces bacteria showed various morphologies (Fig. S21), covering different phage families when visualized using TEM (Figs. 5a). Based on viral metagenomic sequencing, 1059 non-redundant phage contigs were recovered (Table S13). The integrated machine-learning tool [63] revealed that 28% of phages (296 of 1059 viral contigs) were predicted to infect the Actinobacteriota phylum and that 20% of phages (212 of 1059 viral contigs) were associated with the Streptomyces genus (covering 165 species; Table S12). Notably, 8.1% of phages (86 of 1059 viral contigs) were associated with multiple species within the Streptomyces genus, indicative of a broad phage host range (Table S12). Several phage families known for broad bacterial host ranges were identified in the enriched cocktail. These families likely targeted other taxa within Actinobacteriota and Firmicutes, beyond Streptomyces. Additionally, some temperate phages, which were less host-specific and capable of lysogeny, were identified in the cocktail. We found a 45-fold increase of Streptomyces-associated phages in the Streptomyces-enriched phages compared to the original unenriched sludge, indicating substantial enrichment of Streptomyces phages. As 72% (763 out of 1059) of vOTUs in the enriched Streptomyces phage treatment could not be assigned to known bacterial hosts based on the current database, the proportion of Streptomyces-associated phages might be underestimated (Table S12). Also, the significance of enriched, non-Streptomyces-specific phages and their role in shaping bacterial communities, mobilome, and resistome need to be studied in more detail in the future.

Fig. 5

The effect of Streptomyces-enriched phages on ARG and MGE abundances, bacterial community composition, and functioning across 48 agricultural soil samples collected in China. a Transmission electron microscopy (TEM) images of typical phages in Streptomyces enriched-phage suspensions. b Changes in the ARG, MGE, and Streptomyces abundances between phage and control groups averaged over 48 soil samples. c PCoA analysis shows the effect of sampling sites and phage treatment on the ARG and MGE composition based on distance dissimilarity. d Changes in the bacterial community richness (right) and community composition (left) between phage and control groups averaged over 48 soil samples. e The correlation between total ARG and Streptomyces taxa abundances is based on 48 soil samples across China (dots showing means of three replicates per site). f The effect of phages on the abundance (left) and transcriptional (right) activity of Streptomyces-specific phages at the end of the experiment. In panels c and e, shaded clouds represent 95% confidence intervals around the centroids. The asterisks indicate significant differences between phage and control groups based on the Student’s t-test (*p < 0.05, **p < 0.01, and ***p < 0.001)

Across all soil samples, Streptomyces-enriched phages reduced total ARG (p < 0.001, t = 6.9, df = 94), MGE (p = 0.02, t = 2.2, df = 94) and Streptomyces (p < 0.001, t = 3.8, df = 94) abundances, indicative of strong top-down regulation of ARG-carrying bacteria (Fig. 5b). Specifically, the abundance of vancomycin (p < 0.001, t = 8.8, df = 94) and aminoglycoside (p = 0.04, t = 2.0, df = 94) ARGs, and transposase (p = 0.012, t = 2.5, df = 94) and IS91 (p = 0.018, t = 2.4, df = 94) MGEs were significantly reduced by phages (Fig. S22–23). Changes in the dynamics of health risk of ARGs in MAGs in croplands across China (n = 48) were assessed based on a previously established ARG risk framework [40]. Notably, although the number of high-risk ARGs (ranks I and II) in the farm soils represent a relatively small proportion of resistance genes (9%; 45 out of 495), their abundance accounted for 31.4% of all ARGs (Table S15). For example, aadE rank I ARG present in pathogenic Streptomyces bacterium, which confers resistance to aminoglycoside, was detected in all soil samples. After phage treatment, the mean abundances of high-risk and low-risk ARGs decreased by 71.1% and 66.9%, respectively, across all soil samples (Fig. S24). By comparing the removal rate of Streptomyces-encoded ARGs by different phages in both microcosm experiments, Streptomyces-special phages were more efficient than non-enriched phage consortia (83 vs. 57%), possibly due to the higher host specificity of the Streptomyces-enriched phages. Although the magnitude of phage effect on the resistome and mobilome varied between the sampling sites (R² = 0.2961, p < 0.001), the impact of phage treatment was highly significant (R² = 0.5031, p < 0.001, Fig. 5c). While the taxonomic composition of bacterial communities of original soils varied significantly between sampling sites (Fig. S25a), the composition of resistomes did not differ significantly between sampling sites (Fig. S25b). This suggests that the effect of phages on microbial resistance genes was not due to baseline taxonomic differences between soil samples. Furthermore, although soil ARG removal rates were significantly affected by different sampling sites, no significant correlation was observed between the microbial richness and ARG removal rates by phage across all sites (Fig. S26, R² = 0.087, p = 0.56).

Similar to the first experiment, phages reduced the relative abundance of Actinobacteriota (p < 0.001, t = 6.1, df = 94, Fig. S27), leading to reduction in bacterial community richness (p = 0.001, t = 3.1, df = 94) and shifts in overall community composition (R² = 0.29, p = 0.001, Fig. 5e). Specifically, the average abundance of Streptomyces was reduced from 12 to 6.4% after the addition of phage treatment across all sites (Fig. 5b), and the reduction in Streptomyces abundance was associated with the reduction in ARG (R² = 0.55, p < 0.001, Fig. 5f) and MGE (R² = 0.47, p < 0.001) abundances (Fig. S28). Moreover, phage-driven reduction in Streptomyces abundances (p = 0.0018, t = 3.2, df = 94, Fig. S29a) led to more than threefold reduction in the relative activity of Streptomyces (p < 0.001, t = 5.6, df = 10, Fig. S29b) and twofold increase in the transcriptional activity of Streptomyces-specific phages relative to all phages overall (Fig. 5g). Together, these findings suggest that Streptomyces-enriched phages reduced ARGs and MGEs across soils by targeting Streptomyces bacteria that formed a keystone hub for the soil resistome.

Discussion

Understanding how ARGs spread in the complex natural soil environment is important for controlling the rise of antibiotic resistance [71]. Here, we show that Streptomyces plays a critical role in maintaining the diversity of bacterial communities and the spread of ARGs. Specifically, Streptomyces genus formed a crucial ARG-associated hub in soil microbiomes whose abundances could be consistently reduced through the application of unenriched sewage sludge and Streptomyces-enriched phages across different soil samples. Furthermore, the reduction of ARGs was coupled with the reduction in antibiotic synthesis gene abundances and produced antibiotics, which likely further reduced the selective benefit of ARG carriage. Streptomyces is known to play multifunctional ecological roles in antibiotic production, organic matter decomposition, and microbial competition [72]. Thus, their depletion by phages could alter trophic interactions, nutrient availability, or competitive interactions in soil microbiomes. Additionally, phage predation may inadvertently influence non-target organisms through mechanisms such as cross-infection (polyvalence), horizontal gene transfer, or changes in resource availability and niche space following Streptomyces lysis. For example, using phages to selectively remove pathogenic bacterial taxa has been shown to change the community composition and functioning in tomato rhizosphere due to competitive release [73].

We first focused on agricultural soil samples collected from two fields treated by intensively organic waste in China and tested how the phage community isolated from the activated sludge affected the microbiome composition, diversity, and functional gene content specifically focusing on ARGs and MGEs. In addition to clearly controlling bacterial abundances, phages also changed the composition of bacterial communities, which was coupled with clear reductions in ARG and MGE abundances. The reduction in ARGs was linked with the reduction in the abundance of the Streptomyces genus known for carrying antibiotic synthesis and resistance genes. Although previous studies have shown that phage application can impact the native microbial community composition and a few special bacterial strains in soil ecosystems [44, 53, 60, 74, 75], this is the first study demonstrating how phage consortia can be used to modify soil resistome and mobilome at the community level via lysis of keystone bacterial taxa. In line with this, we found that Streptomyces MAGs encoded a relatively high number of ARGs and MGEs among all the MAGs, and for example, three Streptomyces MAGs encoded more than 70 biosynthetic gene clusters, including 20 clusters linked with antimicrobial compounds. Streptomyces are widely distributed in the soil, making up 1–30% of total soil bacterial abundances [76]. Whereas many Streptomyces species have been reported to play important roles in promoting plant growth [77], some Streptomyces species can cause both plant diseases and human infections [78, 79]. Therefore, the role of Streptomyces as an important core group in maintaining soil health requires further validation and should be evaluated from microbial ecological perspectives, including reciprocity, symbiosis, and parasitism.

Reduction in Streptomyces abundances by phages was also linked to reduced production of Streptomyces-encoded antibiotics. This could have had a secondary effect on ARG abundances by reducing the selective advantage of ARG carriage, especially if ARGs were associated with growth or metabolic costs [80]. Random forest model and network analysis further confirmed the importance of Streptomyces genus for soil resistome and mobilome, suggesting that this genus played a keystone role for the presence and movement of ARGs in soil microbiomes. We further verified that Streptomyces-enriched phages had a similar effect in reducing soil resistance genes using a wider range of farmland soils collected across 16 provinces in Southeast China. Although the background of microbial community composition and ARG removal rates varied significantly between soils in different locations, the effect of phages on reducing Streptomyces taxa and the ARGs they carried was consistently significant in all sites compared to the control treatments without phage across all sites. Moreover, no significant correlation was observed between the microbial richness and ARG removal rates by phage across all sites. This suggests that the phage application treatment was not significantly influenced by soil type or microbial community composition, likely due to the strong host specificity of the phages.

Traditionally, phages are used to target specific pathogenic bacteria to selectively reduce their relative abundance, for example, when treating human infections [81] or when bio-controlling plant pathogenic bacteria [10, 82]. While such applications often rely on using phages as cocktails, they most often contain phages that can infect only one certain bacterial species [10, 83]. Here, we expand the phage cocktail approach to much more diverse phage consortia that are less specific but can still predictably change soil microbiome functioning. We also used traditional enriching techniques to increase the relative abundance of Streptomyces-specific phages to show that targeting this bacterial genus was important for reducing ARG abundances across several agricultural soils. The enriched phage cocktails could be further optimized into phage-based biotechnological products that reduce the spread of ARGs in the agricultural environment, offering sustainable solutions for public health within the One Health framework. In particular, based on the previously established framework of different risks of ARGs [40], we investigated the effects of phage consortium treatments on different risk classes of ARGs in soils using metagenomic binning. Whereas the number of high-risk ARGs in soil was relatively small, phage treatment was effective in reducing high-risk ARG abundances in soil samples covering main croplands in China. Streptomyces play important roles in the soil, being significant carriers of mobile ARGs, and secreting large amounts of antimicrobials that can regulate the composition of the bacterial communities. Consistently, the Streptomyces-enriched phage treatments were more effective in removal of high-risk ARGs compared with the non-enriched phage cocktails. While several techniques, such as composting [53], have been developed to reduce ARGs from manure and human waste, we here demonstrate that phages could potentially be used as a biological treatment to remove ARGs from contaminated waste and soils. Such manipulations however rely on the activity of lytic phages—an effect which could be complicated by the presence of lysogenic phages capable of promoting ARG dissemination via horizontal gene transfer [84]. Hence, a targeted approach where lytic phages specific to important ARG-carrying bacterial taxa are used is likely to be more successful.

While reducing the prevalence of ARGs is desirable, concomitant loss of antibiotic synthesis genes could have broader implications for soil functionality [85]. Future studies should evaluate whether phage-mediated lysis of Streptomyces taxa affects plant growth and health. Notably, soil suppressiveness (i.e., the soil’s ability to deter plant pathogens) often relies on the production of antimicrobial compounds [44, 68, 69], with the Streptomyces genus playing a pivotal role in synthesizing pathogen-suppressing antimicrobials [86,87,88]. Using phages that target antimicrobials-producing bacteria could hence make agricultural soils more prone to pathogen invasions, as has been shown recently in the case of bacterial wilt caused by R. solanacearum [82]. Yet, the use of phages that target ARG-carrying bacteria should be less problematic in waste treatment biotechnology where phages can be applied for example during composting or directly inoculated onto sand filters during biological water treatment.

The present study investigated the potential effects of phage-mediated modification of soil microbiomes based on extracted bacterial communities using microcosm experiments. While such experiments may not fully capture the natural complexity and heterogeneity of native soil microbiomes, they offer reproducibility and high experimental control to address causal hypotheses. To move beyond one soil type and bacterial community, we conducted phage treatment experiments across 48 soils derived from 16 provinces throughout China. These experiments confirmed that Streptomyces-enriched phages consistently reduced ARG levels in different background soils with different microbial communities. While the use of a diverse phage consortium offers advantages such as a broad host range and improved robustness against resistance development [89], applying diverse phage consortia into the soil could have unexpected off-target effects, resulting in unwanted phage infections with multiple non-target microbial taxa. This could lead to unpredictable shifts in bacterial community composition, functioning, and interaction networks, which should be studied in more detail to ensure the safety of phage applications. For instance, phages with overlapping host ranges may alter bacterial competition dynamics, while phage-mediated horizontal gene transfer could inadvertently mobilize ARGs or virulence factors. Furthermore, the heterogeneous spatial structure of soil may influence phage infection efficiency, resulting in uneven ARG suppression in space and across different bacterial taxa. These considerations highlight the need for targeted phage design that considers potential ecological risks for the surrounding microbiota and quantifies the phage survival and degradation in the soil after bacterial infections and lysis. Moreover, several limitations of phage therapy should be considered, such as the high host specificity leading to a narrow range of lysed hosts, the potential evolution of phage resistance associated with long-term use, challenges related to the stability and preservation of phage preparations, and the risk of some phages carrying virulence genes.

Conclusions

Our multi-omics analysis combined with phage biocontrol experiments revealed that Streptomyces are the keystone taxa for ARG carriage and dissemination in agricultural microbiomes and that phages can be used to regulate the soil microbiome resistome and mobilome by reducing the relative abundance of Streptomyces taxa. Metagenomic-informed phage biocontrol design could hence be used to effectively control desired bacterial taxa to trigger community-wide changes in the composition and functioning of other microbiomes beyond agricultural soils. Our approach offers a novel approach for the identification and verification of keystone taxa in soil microbiomes and expands the scope of phage biocontrol from bacterial species to functional gene level in complex environments.