Insights into antibiotic resistomes from metagenome-assembled genomes and gene catalogs of soil microbiota across environments

Sample collection, DNA extraction, and sequencing

Soil samples were collected from four Chinese provinces: Sichuan (25), Yunnan (28), Guizhou (28), and Jiangsu (30). These regions were chosen to represent diverse geographic and ecological settings, including mountainous and plateau areas in the southwest and agricultural plains in the east. All four provinces are characterized by intensive human activities and livestock production, providing representative environments for studying human-animal-soil interactions. Sampling across these regions allowed the capture of ecological and agricultural variation relevant to soil microbial communities and ARGs. In total, 111 soil samples were obtained. With detailed metadata including location, altitude, and GPS coordinates recorded for each site (Table S1). At each sampling site, surface soil (0–1 cm depth) was collected according to a five-point sampling protocol using a sterile soil shovel. The five sub-samples obtained from each point were then combined to form a single composite sample. Subsequently, the composite soil was immediately passed through a 2-mm sieve to remove plant debris and gravel. Approximately 1–3 g of the composite soil was collected for each sample. To preserve sample integrity, all samples were frozen on dry ice after collection and subsequently stored at −80 °C. Total genomic DNA was extracted from soil samples using the Mag-Bind^® Soil DNA Kit (Omega Bio-tek, Norcross, GA, USA) following the manufacturer’s instructions. DNA concentration and purity were determined using a NanoDrop 2000 spectrophotometer, and integrity was assessed by 1% agarose gel electrophoresis.

DNA samples were fragmented to an average size of approximately 350 bp using a Covaris M220 instrument (Gene Company Limited, China) for the preparation of paired-end sequencing libraries. Library construction was carried out with the NEXTFLEX Rapid DNA-Seq Kit (Bio Scientific, Austin, TX, USA). Adapters comprising the full complement of sequencing primer binding sites were ligated to the blunt-ended DNA fragments. Paired-end sequencing was conducted on an Illumina NovaSeq 6000 platform (Illumina Inc., San Diego, CA, USA) at Majorbio Bio-Pharm Technology Co., Ltd. (Shanghai, China), utilizing the NovaSeq 6000 S4 Reagent Kit v1.5 (300 cycles) in accordance with the manufacturer’s protocols (https://www.illumina.com/).

Metagenomic assembly

Quality control and sequence assembly were performed using the Majorbio Cloud Platform (http://www.majorbio.com). Adapter sequences and low-quality reads (length <50 bp or quality value <20) were removed using fastp (version 0.23.0) (Chen et al., 2018). Host-derived reads were removed by aligning against the human (GRCh38.p13), chicken (GRCg7b), and pig (Sscrofa11.1) genomes using BWA (http://bio-bwa.sourceforge.net, version 0.7.9a) (Li & Durbin, 2009). Clean reads were assembled using MEGAHIT (version 1.1.2) (Li et al., 2015).

Metagenomic binning and MAG quality control

Three tools, Metabat2 (version 2.12.1) (Kang et al., 2015), MaxBin2 (version 2.2.5) (Alneberg et al., 2014), and CONCOCT (version 0.5.0) (Wu, Simmons & Singer, 2016), were used for binning. The resulting bins were consolidated using DAS Tools (version 1.1.0) (Sieber et al., 2018). RefineM (version 0.0.24) (Parks et al., 2017) was used to filter contigs based on genomic features such as GC content, tetranucleotide signatures, coverage, and taxonomy. MAG completeness and contamination were assessed with CheckM (version 1.0.12) (Parks et al., 2015) using lineage-specific marker genes. Only MAGs with ≥ 50% completeness and <10% contamination were retained (Parks et al., 2017).

Pairwise comparisons were performed using Mash (Ondov et al., 2016) and dereplication was conducted with dRep (version 3.4.2) (Olm et al., 2017), retaining the highest-quality MAG for each cluster at ≥ 99% average nucleotide identity (ANI). MAG coverage was calculated using CoverM (version 0.6.1). Taxonomy was assigned using Genome Taxonomy Database Toolkit (GTDB-Tk) (version 2.3.0) (Parks et al., 2018) based on 120 universal single-copy marker genes from the Genome Taxonomy Database (GTDB).

Gene prediction and functional annotation

Gene prediction for all MAGs was performed using Prodigal (version 2.6.3) with the -p meta parameter (Hyatt et al., 2010). Only genes with a nucleic-acid length ≥ 100 bp were retained and subsequently translated into amino-acid sequences. The annotated protein sequences were automatically aligned against multiple databases using Diamond (Version 0.8.35, E-value ≤ 1e−5), including Kyoto Encyclopedia of Genes and Genomes (http://www.genome.jp/kegg/), COG (https://www.ncbi.nlm.nih.gov/research/cog-project/), CAZy (http://www.cazy.org/), CARD (https://card.mcmaster.ca/), and MGE (https://github.com/KatariinaParnanen/MobileGeneticElementDatabase).

Data analysis

In this study, violin plots, clustered heatmaps, pie charts, bar plots, correlation scatter plots, and polar bar plots were generated using the Microbioinformatics platform (https://www.bioinformatics.com.cn). The MAG species Sankey diagram, correlation heatmap, and circular phylogenetic tree were generated using Python (v2.7.10). UpSet Venn diagrams, dual-matrix correlation heatmaps, KEGG histograms, COG annotation classification statistics, antibiotic resistance gene prediction classification statistics, and gene potential mobility analysis charts were produced using the Majorbio platform (http://www.majorbio.com). Additionally, UpSet Venn diagrams and interactive heatmap bar plots were generated with the assistance of the Paisainuo Gene Cloud platform (https://www.bioinformatics.com.cn).

To investigate the relationships between species distribution and environmental variables (including altitude, latitude, and longitude), as well as the associations between mobile genetic element MGE types and ARG categories, Spearman’s correlation coefficient was employed for analysis. Furthermore, Spearman’s rank correlation test was used to evaluate the correlation between the abundance of ARGs and MGEs. Spearman’s correlation calculation was employed, and significance was denoted by p-values, with * indicating p < 0.05, ** indicating p < 0.01, and *** indicating p < 0.001.

Assembly of 1,136 microbial genomes from soil

We constructed a metagenomic assembly gene catalog using sequencing data from 111 soil samples across four Chinese provinces (Fig. 1). High-throughput sequencing generated 1.276 Tb of raw data, averaging 11.50 Gb per sample. After stringent quality filtering, 1.258 Tb of high-quality data were retained, reflecting a 98.59% retention rate and an average effective sequencing depth of 11.33 Gb per sample (Table S1).

Dereplication at an average nucleotide identity (ANI) threshold of ≤ 99% and quality filtering yielded 1,136 MAGs that met our predefined standards for medium-quality assemblies (≥ 50% completeness and <10% contamination; Table S2). Among them, 176 MAGs met the MIMAG standard for high-quality MAGs defined by the Genome Standards Consortium (≥ 90% completeness and ≤ 5% contamination) (Bowers et al., 2017), while the remaining 960 were classified as medium-quality. The average completeness across all MAGs was 74.19%, with a mean contamination of 3.40% (Table S2).

To examine the genome characteristics of these soil-derived MAGs, we used metrics including total genome length, contig number, predicted coding sequences (CDS), GC content, N50, tRNA genes, rRNAs, and repetitive elements (Fig. S1A). The total genome length of all MAGs combined was 3.17 Gb. On average, each MAG contained 725 contigs, ranging from 1 to 3,135, indicating substantial variation in final assembly continuity among MAGs. We predicted 3,548,612 CDS in total, with a mean of 3,123 CDS per MAG, offering a rich resource for functional annotation.

The average GC content was 58.75%, ranging from 28.83% to 75.07%, reflecting the diversity characteristics of the soil microbial genomes. It is worth noting that the upper limit of 75.07% GC content comes from a Cellulosimicrobium funkei genome. Compared with the NCBI database, this value has reached the highest level of GC content reported for this species so far.

The overall average N50 was 12,271 bp. Notably, according to the MIMAG criteria (≥ 90% completeness and ≤ 5% contamination) (Bowers et al., 2017), the 176 high-quality MAGs also exhibited greater assembly continuity, with an N50 of up to 7,575,282 bp—considerably higher than the average of 6,364,488 bp observed for medium-quality MAGs (Table S3). In addition, we identified 36,249 tRNA genes, 684 rRNA genes, and 74,838 repetitive elements across the MAGs (Fig. S1A; Tables S3–S5).

Microbial classification and composition

To classify the 1,136 MAGs, we aligned their sequences with the GTDB. The taxonomic distribution improved at higher ranks: 1,115 MAGs were assigned to a class, 1,047 to an order, 891 to a family, 223 to a genus, and 144 to a species (Table S6). All MAGs were taxonomically annotated using GTDB (Fig. S1B), which revealed 23 bacterial phyla (1,103 MAGs) and three archaeal phyla (33 MAGs). The most abundant bacterial groups included Actinomycetota (324 MAGs), Pseudomonadota (252 MAGs), and Bacteroidota (119 MAGs). Among archaeal MAGs, 20 out of 33 were unclassified at the species level. These archaeal genomes were distributed across Thermoproteota (31 MAGs), Thermoplasmatota (1 MAG), and Methanobacteriota (1 MAG). The average relative abundance of archaeal phyla is approximately 7.12% (Table S6).

Using MAG-based Sankey diagrams, we analyzed species abundance across samples at multiple taxonomic levels, visualizing patterns of community structure and regional distribution. At the phylum level, Actinomycetota dominated the top 30 most abundant MAGs. As the taxonomic resolution increased to genus and species, the Sankey diagram displayed more complex branching. Some species showed regional enrichment; for instance, MAG442 (Lactobacillus amylovorus) was significantly enriched in Yunnan Province (Fig. 2A). Additionally, correlations were observed between MAGs abundance and environmental variables such as altitude, latitude, and longitude. For example, Dormibacterota abundance was positively correlated with altitude, while Armatimonadota abundance was positively associated with both latitude and longitude (Fig. 2B).

Further analysis revealed patterns of microbial community composition across provinces. Among the four regions, 58 MAGs were shared, including those from Pseudomonadota (18 MAGs), Actinomycetota (17 MAGs), Thermoproteota (14 MAGs), Cyanobacteriota (6 MAGs), Chloroflexota (2 MAGs), and Nitrospirota (1 MAG) (Fig. S2). We also identified 651 unique MAGs: 261 from Sichuan, 216 from Yunnan, 27 from Guizhou, and 147 from Jiangsu (Fig. 2C). Bacterial taxa found in more than 90% of samples from a region were defined as core bacteria. In Guizhou, five phyla met this threshold, indicating their wide distribution and potential ecological importance (Table S7).

To explore regional differences, the top 10 MAGs with the highest relative abundance were selected from the shared and endemic MAGs in each province in this study. We compared the top 10 most abundant MAGs, both shared and region-specific, within each province. At the genus level, shared MAGs in Yunnan, Guizhou, and Jiangsu showed greater similarity compared to those from Sichuan. Genera such as Micrococcus (MAG1095), Arthrobacter (MAG845), and VBCG01 (MAG104) were dominant in Yunnan, Guizhou, and Jiangsu, while Chroococcidiopsis (MAG929), Arthrobacter (MAG802), and Paracoccus (MAG292) were more abundant in Sichuan, indicating a distinct MAG composition in that region (Fig. S3).

Genus-level analysis of region-specific MAGs also demonstrated clear provincial differentiation. Dominant genera included VAYN01, UBA4720, and JACDBZ01 in Sichuan; o_Gloeomargaritales; g_Unclassified in Guizhou; Sphingomonas_I, 40CM-4-68-19, and CADCTB01 in Yunnan; and Sphingomicrobium and VAYN01 in Jiangsu (Table S6).

Functional annotations of MAGs

The assembled MAGs were functionally annotated using three major databases: COG, KEGG, and CAZy. We found that 42.69% (2,643,751), 30.8% (2,346,501), and 3.47% (127,544) of the predicted proteins had at least one COG, KEGG, and CAZy function, respectively (Fig. 3A). According to the COG classification, the MAGs were categorized into four main groups: Cellular Processes and Signaling, Information Storage and Processing, Metabolism, and Poorly Characterized functions. These groups include functional descriptions of 25 specific COG types, with Metabolism being the most prevalent (Fig. 3A; Table S8).

KEGG annotations revealed that protein functions were grouped into six major metabolic systems. Approximately 77.52% of the annotated genes were involved in metabolic functions (Fig. 3A), particularly those related to carbohydrate and amino acid metabolism (Fig. 3B). These results were consistent with the COG data, reinforcing the conclusion that soil microbes predominantly perform metabolic roles. The abundance of carbohydrate- and amino acid-related pathways is likely due to the soil environment’s characteristics, such as abundant humus, plant root exudates, and microbial residues (Abdelrahman et al., 2016; Chantigny, Olk & Angers, 2025). These characteristics provide ample substrates for these metabolic processes and drive soil microorganisms to adapt to such conditions.

To further analyze genes related to carbohydrate metabolism, we used the CAZy database (Fig. 3A). A total of 127,544 CAZyme-encoding genes were identified, classified into seven CAZyme categories. Glycosyl transferases (GT) were the most abundant, followed by glycoside hydrolases (GH) and carbohydrate esterases (CE). These genes were unevenly distributed across MAGs. Notably, GH and GT were especially enriched in Actinomycetota and Pseudomonadota (Fig. 3C). Among all MAGs, MAG975 (g_Actinocrispum; s_Unclassifiedand) contained the highest number of functionally annotated genes across all three databases: 6,860 in COG, 5,775 in KEGG, and 517 in CAZy (Table S9).

To assess functional differences among region-specific MAGs, we analyzed the top 10 relatively abundant unique MAGs from each province. Genes related to cell wall/membrane/envelope biogenesis and amino acid transport/metabolism were predominant in MAGs from Sichuan, while MAGs from Guizhou, Yunnan, and Jiangsu were mainly enriched in genes associated with general function prediction, signal transduction, and amino acid metabolism (Fig. S4; Table S8). At the genus level, KEGG analysis showed that UBA4720 (Sichuan), Ferruginibacter (Guizhou), Sphingomonas_I (Yunnan), and Sphingomicrobium (Jiangsu) had relatively high proportions of annotated genes. Functional pathways common across these genera included global and overview maps, carbohydrate metabolism, and amino acid metabolism (Fig. S5A; Table S1). Similarly, CAZy analysis indicated that genera such as UBA4720, Ferruginibacter, Sphingomonas, and Flavisolibacter had relatively high proportions of CAZyme-encoding genes (Fig. S5B).

Antibiotic resistance gene profiling

To investigate antibiotic resistance in the soil microbiome, ARGs in 1,136 MAGs from 111 soil samples were annotated using the CARD database. A total of 35 drug resistance classes, comprising 162 distinct ARG types, were identified (Table S10). Among all phyla, Actinomycetota and Pseudomonadota harbored the greatest numbers of ARGs. ARGs were also detected in archaeal phyla such as Thermoproteota and Methanobacteriota (Fig. 4A). Among these, the phylum Methanobacteriota was found to carry 34 ARGs, encompassing 14 distinct ARG classes (Table S10). Multidrug, peptide, and glycopeptide ARGs were the most widely distributed in the soil environment (Fig. 4B).

Genera QHVT01, Sphingobacterium, and UBA5704 contained the greatest number of ARGs per genome among taxa represented by at least five MAGs (Fig. 4C), suggesting their potential as ARG reservoirs. The widespread presence of ARGs raises concerns about environmental dissemination through food chains or soil-water-plant systems, potentially threatening ecological and human health.

MAG975 (g_Actinocrispum; s_Unclassified) and MAG635 (g_JADGHX01; s_Unclassified), found only in Yunnan Province, harbored the highest number of ARGs according to CARD predictions (Table S9), indicating possible high resistance potential and unique ecological adaptability. All 1,136 MAGs contained at least five ARGs (Table S10), further underscoring the pervasive nature of resistance genes in soil microbial communities.

Of particular note, one Escherichia coli strain (MAG817) was identified with 81 distinct ARG types spanning 27 drug resistance classes (Table S11). These included resistance to multidrug, peptide, glycopeptide, tetracycline, macrolide, aminoglycoside, and fluoroquinolone antibiotics. Given the pathogenic nature of E. coli, this strain could potentially be a drug-resistant superbug.

ARGs were also observed in several probiotic species, including Arthrobacter oxydans (MAG82, MAG174, MAG670, MAG845), Glutamicibacter arilaitensis (MAG79, MAG161, MAG537, MAG825), Pseudomonas_E helleri (MAG467), Priestia megaterium (MAG750), and Lactobacillus reuteri (MAG783) (Table S10).

We also compared ARG abundance among soil-specific MAGs from the four provinces. Sichuan, Yunnan, and Jiangsu soils contained significantly more ARGs than Guizhou, with aminoglycoside and phosphonic acid resistance genes being dominant across all regions. In terms of ARG diversity, Sichuan had twice as many antibiotic classes (10) as Guizhou (5), while Yunnan exhibited the greatest diversity with 12 types. Jiangsu followed closely with nine types (Fig. 4D).

Mobile genetic elements

MGEs are instrumental in the horizontal transmission of ARGs among microbial cells. Understanding their distribution and the relationship between MGEs and ARGs is essential for assessing antibiotic resistance dynamics. A total of 46,662 MGE-related genes were identified across the 1,136 MAGs by aligning gene catalog protein sequences with the MGE Database. These genes were grouped into 68 distinct MGEs and further classified into five types: transposase, integrase, recombinase, conjugative transfer protein, and transposon. Among these, transposase genes were the most abundant, suggesting their major contribution to the dissemination of ARGs (Table S12).

MGEs can carry and promote the spread of ARGs through horizontal gene transfer (HGT). Thus, the presence and distribution of ARGs and MGEs can serve as indicators of the degree of antibiotic resistance contamination. To explore regional mechanisms of ARG transmission, we analyzed MGE content within the unique MAGs from the four provinces. Results revealed clear regional differences: Sichuan demonstrated the highest diversity (46), followed by Yunnan (36), Jiangsu (27), and Guizhou (23). Regarding province-specific MGE types, Sichuan had 12, Yunnan 5, Guizhou 2, and Jiangsu 1 (Fig. 5A). Among all regions, transposase, XerD (site-specific tyrosine), and Tn3 were the dominant classes (Fig. 5B). Further analysis of gene mobility across the 1,136 MAGs showed that among contig regions shorter than 5 kb, ARGs most frequently appeared adjacent to transposase genes. The genome of MAG975 (g_Actinocrispum; s_Unclassified) contained the highest number of such co-located regions. Among these, the bcrA gene (n = 13) was most commonly associated with transposase genes within the same contigs. This suggests a likely mechanism for the horizontal transfer of bcrA mediated by transposases in this genome (Table S12).

Previous studies have highlighted E. coli as a potential pathogen across various environments. To evaluate the risk of ARG transfer at the strain level, we examined the distribution of MGEs in the genome of MAG817 (E. coli). Transposase genes were located adjacent to multiple ARGs within contigs shorter than five kb. For example, a transposase gene was closely linked to the efrB gene (Fig. 5C). Given that MAG817 (E. coli) was detected in eight out of 111 samples and has a relatively high abundance, this E. coli strain may represent a heightened risk for ARG transmission through MGE-mediated horizontal transfer.

To quantify the relationship between MGEs and ARGs, we performed Spearman and linear regression analyses. A significant positive correlation was observed (y = 0.2349x + 0.8653, R² = 0.5229, p < 0.001), indicating that higher MGE abundance may facilitate ARG proliferation (Fig. 6A). This correlation was also supported by a heatmap (Fig. 6B). Notably, when delving into the individual components of the mobilome, we found a more complex relationship. While transposase genes exhibited a strong positive correlation with peptide ARGs (Fig. 6C), intriguingly, the abundance of TraG (a key ATPase gene essential for bacterial conjugation) was negatively correlated with the total ARG abundance (Fig. 6B).