Organic fertilization co-selects genetically linked antibiotic and metal(loid) resistance genes in global soil microbiome

Soils are major reservoirs of microbial diversity and hubs for antibiotic resistance genes (ARGs), posing risks of transfer to pathogens. Agricultural soils receive substantial anthropogenic inputs via fertilization, especially organic fertilizers (e.g., manure) that can introduce both antibiotics (tetracyclines, sulfonamides, quinolones, macrolides) and metal(loid)s (Cu, Zn, As, Cr, Hg, Mn). Prior studies often inferred co-occurrence of ARGs and metal(loid) resistance genes (MRGs) from abundance correlations, but direct evidence for genetic linkage (co-localization on the same genetic elements such as plasmids) in complex soil communities was lacking. While mining-impacted environments show ARG–MRG associations and some isolate genomes exhibit ARG–MRG linkages, it remained unclear whether organic fertilization co-selects genetically linked ARG–MRG pairs at a global scale in agricultural soils. This study tests the hypothesis that organic fertilization promotes genetic linkage and co-selection of ARGs and MRGs in soil microbiomes, evaluates transcriptional co-regulation under arsenic stress, identifies host genomes carrying co-resistant determinants, and maps global risk of antibiotic and metal(loid) co-resistant bacteria (AMRB).

Previous work showed increased antibiotics and metal(loid) inputs to soils via manure fertilization, with long-term manure use elevating soil metal concentrations. Studies in mining-impacted environments observed frequent co-occurrence and genetic associations between multidrug ARGs and multimetal MRGs, suggesting metal-driven co-selection. Some bacterial isolates showed specific ARG–MRG linkages (e.g., Zn with beta-lactam resistance; bacitracin with polymyxin; Cd with aminoglycoside resistance). However, most agricultural soil studies relied on abundance correlations without demonstrating genetic linkage within complex communities. The present work leverages expanded metagenomic datasets to directly assess genetic co-localization of ARGs and MRGs, their regulation under As stress, and associated health risks.

- Global metagenome collection: Retrieved 511 agricultural topsoil (depth <20 cm) metagenomes (Illumina paired-end, >1 Gb/sample) from SRA, ENA, DDBJ, and NGDC (search date 2022-10). Plant-associated (rhizosphere/rhizoplane) samples were excluded. - Grouping by fertilization: 109 samples identified as organic fertilized (OF) and 109 as non-organic fertilized (NOF) from metadata. For 293 with missing fertilization info, a random forest (RF) classifier (F1=0.97) trained on genus-level abundances of the labeled samples assigned OF/NOF labels, yielding 227 NOF and 284 OF samples. - Processing: Reads trimmed with fastp; assemblies per sample with MEGAHIT (min contig 500 bp). ORFs predicted with Prodigal (-p meta). MAGs recovered and refined with MetaWRAP; quality assessed by CheckM; dereplicated with dRep; taxonomy via GTDB-Tk; novelty via MiGA. MAG abundances computed by mapping with BLAT and normalized to coverage per genome equivalents (CPG; MicrobeCensus). - Identification of ARGs/MRGs/MGEs/VFGs and genetic linkage: ORFs annotated against SARG (ARGs), a custom MRG database (BacMet, copper resistance proteins, UniProt As genes), a custom MGE database (ICEberg, INTEGRALL, ISfinder, nanoMGEs, Transposon Registry), and VFDB using DIAMOND. Criteria: ARG/MGE/VFG: ≥60% identity, ≥70% coverage, e≤1e−5 (reads-based screening used ≥90% identity, ≥25 aa alignment, e≤1e−7); MRG: ≥50% identity, ≥70% coverage, e≤1e−5. Contigs >1000 bp classified as plasmid/chromosome by PlasFlow. ARG–MRG-carrying contigs (AMCCs) were contigs with ≥1 ARG and ≥1 MRG. AMCC abundance per sample computed as the mean of ORF-level RPKG across ARG/MRG-like ORFs on AMCCs, representing copies per cell. Coexistence networks (ARG/MRG subtypes) visualized in Gephi. - Co-occurrence networks and diversity: Taxonomy classified by Kraken2/Bracken; genus-level co-occurrence networks constructed (RMT-based thresholds; RMThresh), robustnes assessed by random removals; alpha/beta diversity computed with R packages. - Arsenic stress validation: Collected 12 paddy soils from four Chinese sites. Classified into relatively low (avg 9.32 mg/kg As) and high As (avg 19.49 mg/kg As) groups (threshold 15 mg/kg). Performed metagenomes and metatranscriptomes (rRNA removal with SortMeRNA; multitrim). Assembled AMCCs carrying vanS-arsM, macB-wtpC, and mdtB-mdtC. Mapped meta-omics reads to AMCCs (BLAT ≥90% identity, ≥70% coverage). Computed AMCC abundance (metagenomes) and transcriptional activity (metatranscriptome RPKM normalized by metagenome abundance). - Drivers of coexistence: Partial least squares structural equation modeling (PLS-SEM) tested effects of fertilizer type, soil properties (pH, organic carbon content and density, total nitrogen), climate (mean annual temperature, annual precipitation), socioeconomics (population, GDP, HDI), community diversity (Shannon, richness), and elemental metabolism gene abundances (C, N, P, S, K). Model quality assessed by DG.rho, loadings, AVE, GOF. Random forest identified predictor importance (increase in MSE). - MAG analysis: Identified MAGs carrying AMCCs (completeness ≥70%, contamination ≤5%). Screened for MGEs and VFGs in MAGs, assessed taxonomy and novelty, quantified class-level composition and abundances (CPG) in NOF vs OF. - Machine learning risk mapping: Compiled 34 geospatial features (e.g., EarthStat yield, FAO livestock density, HDI). Discretized AMRB abundance (host MAGs carrying AMCCs with MGEs/VFGs) into six levels by K-means. Split data 80/20; trained ANN, KNN, SVM, XGBoost, RF with 10-fold cross-validation, hyperparameter tuning, and feature selection. Selected best RF model (23 features, n_estimators=50, max_depth=27, min_samples_split=7, min_samples_leaf=1; F1=0.722 on test). Predicted global AMRB risk on croplands (Google Earth Engine coordinates) at 0.083° resolution; visualized in Python.

- Diversity and community changes: OF soils had lower microbial Shannon diversity but more complex and robust co-occurrence networks than NOF (p<0.001). OF soils showed significantly higher richness of ARGs, risk ARGs, and MRGs: 60.3 vs 44.1 (ARGs), 5.5 vs 2.3 (risk ARGs), and 84.7 vs 69.1 (MRGs). In OF soils, risk ARG abundance correlated positively with MRGs (R²=0.10, p<0.001) and MGEs (R²=0.19, p<0.001), not observed in NOF. - Genetic linkage (AMCCs): AMCC abundance was higher in OF than NOF soils (0.44 vs 0.11 copies per cell; p<0.001). AMCCs with MGEs and VFGs were more frequent in OF: MGEs 4.59% vs 0.86%; VFGs 2.72% vs 0.05%; both 1.14% vs 0.98%. AMCC types were ~3× more diverse in OF than NOF (63 vs 22 types). Multidrug-Zn was dominant but reduced in proportion in OF (NOF: ~71%; OF: ~45%), while OF harbored many additional ARG–MRG combinations (e.g., multidrug with Al, As, Au, Cd, Cr, Cu, Fe, Se; plus kasugamycin, vancomycin, bacitracin, tetracycline, fosmidomycin, rifamycin, beta-lactam, quinolone linkages). A recurrent plasmid gene cluster carrying As resistance (acr3, arsC2, arsH, arsR3), transposons (Tn6183, Tn6302), and sulfonamide resistance (sul2) was identified across OF soils, evidencing ARG–MRG co-location. - Arsenic stress experiment: In 12 Chinese paddy soils, AMCCs (vanS-arsM; macB-wtpC; mdtB-mdtC) had higher abundances under high As vs low As (0.52 vs 0.13 copies per cell) and higher transcriptional activity (0.44 vs 0.02). Metatranscriptomic read recruitment showed even expression across ARGs and MRGs within AMCCs, indicating co-regulation under As stress. - Drivers of coexistence: PLS-SEM (GOF=0.55) explained 41% of variance in AMCC abundance. Fertilizer type, soil properties, climate (MAT, AP), socioeconomics (HDI), and biofactors (ARG/MRG abundances, potassium metabolism) significantly contributed. Fertilizer effects were largely indirect via soil properties and ARG/MRG abundances. Climate had a stronger positive total effect in NOF (0.53) than OF (-0.14), while ARG/MRG abundances were predominant drivers in OF. - Host genomes (MAGs): 95/617 (≈15%) high-quality MAGs carried AMCCs, all containing MGEs and VFGs, indicating HGT potential and pathogenicity risk. Many represented novel taxa (3 potential new genera, 72 new species, 17 new subspecies). AMRB abundance was higher in OF vs NOF (0.011 vs 0.005 CPG). OF AMRB were dominated by Gammaproteobacteria (44%), Bacteroidia (20%), Actinomycetia (15%); NOF by Acidobacteriae (41%), Nitrospiria (18%), Gammaproteobacteria (16%). Pathogenic genera (Enterobacter, Pseudomonas) had high VFG counts. - Global risk mapping: A tuned RF model achieved F1=0.722 on the test set. Global predictions at ~0.083° resolution indicated higher AMRB risk (levels 4–6) in central North America, Eastern Europe, Western Asia, and Northeast China compared to other regions (levels 1–3).

The study provides direct evidence that organic fertilization co-selects genetically linked ARGs and MRGs in agricultural soils worldwide. Compared to NOF, OF soils exhibit increased AMCC abundance and diversity, enriched MGEs and VFGs on AMCCs, and broader ARG–MRG linkage types, supporting both cross-resistance (e.g., mdtABC with Zn) and co-resistance mechanisms (multidrug with multiple metals). Metatranscriptomic validation under arsenic stress demonstrates coordinated expression of co-located ARGs and MRGs, confirming functional co-selection in situ. Community network complexity and enriched MGEs in OF soils suggest enhanced potential for horizontal gene transfer, facilitating ARG spread within microbial communities. SEM and RF analyses indicate that fertilization practices influence soil properties and ARG/MRG abundances, which in turn drive ARG–MRG coexistence, with climate and socioeconomics modulating these pathways. The identification of numerous AMRB MAGs, including potential pathogens harboring MGEs and VFGs, underscores dissemination risks. The machine learning-based global map highlights geographic hotspots where agricultural practices and environmental conditions likely elevate co-resistance risks, aligning with regions of intensive agriculture and human activity. Overall, findings address the research question by demonstrating genetic linkage and co-selection of ARGs and MRGs associated with organic fertilization, emphasizing implications for public health and agricultural management.

This work demonstrates that organic fertilization globally promotes the genetic linkage and co-selection of antibiotic and metal(loid) resistance genes in agricultural soils. OF soils harbor higher abundances and a wider range of AMCC types, increased MGE/VFG association, and active co-expression under arsenic stress. Drivers include fertilization-mediated changes in soil properties and ARG/MRG abundances, with climate and socioeconomic factors exerting direct and indirect effects. The study identifies diverse AMRB host genomes, many novel, and constructs a machine learning-based global risk map highlighting hotspots in central North America, Eastern Europe, Western Asia, and Northeast China. These results underscore the need to monitor and manage organic fertilization practices, reduce inputs of antibiotics and metals, and develop interventions to curb ARG–MRG co-selection and dissemination. Future research should refine causal pathways, integrate temporal datasets, validate transfer potential experimentally, and evaluate management strategies that mitigate co-resistance while maintaining soil health.