Global profiling of antibiotic resistomes in maize rhizospheres

The spreading of antimicrobial resistance (AMR) in crops and food products represents a global concern. In this study, we conducted a survey of resistomes in maize rhizosphere from Michigan, California, the Netherlands, and South Africa, and investigated potential associations with host bacteria and soil management practices in the crop field. For comparison, relative abundance of antibiotic resistance genes (ARGs) is normalized to the size of individual metagenomes. Michigan maize rhizosphere metagenomes showed the highest abundance and diversity of ARGs, with the detection of blaTEM-116, blaACT-4/-6, and FosA2, exhibiting high similarity (≥ 99.0%) to those in animal and human pathogens. This was probably related to the decade-long application of manure/composted manure from antibiotic-treated animals. Moreover, RbpA, vanRO, mtrA, and dfrB were prevalently found across most studied regions, implying their intrinsic origins. Further analysis revealed that RbpA, vanRO, and mtrA are mainly harbored by native Actinobacteria with low mobility since mobile genetic elements were rarely found in their flanking regions. Notably, a group of dfrB genes are adjacent to the recombination binding sites (attC), which together constitute mobile gene cassettes, promoting the transmission from soil bacteria to human pathogens. These results suggest that maize rhizosphere resistomes can be distinctive and affected by many factors, particularly those relevant to agricultural practices.


Introduction
According to the Centers for Disease Control (CDC) and World Health Organization (WHO), antimicrobial resistance (AMR) is considered one of the greatest challenges for human health (CDC 2019;WHO 2014). It is estimated that 1.27 million deaths in 2019 were directly caused by AMRrelated infections across the world (Murray et al. 2022). Unfortunately, if the situation is left unsolved, AMR can cause as many as 10 million deaths per year, making it the number 1 killer by 2050 (O'Neill 2016). The spreading of AMR and its transmission into human pathogens through food production and consumption have significantly contributed to this global threat (Arias and Murray 2009;Durso and Cook 2014). Furthermore, the extensive use of antibiotics in animal husbandry and the subsequent application of manure as agricultural soil amendments can worsen the issue of AMR in agriculture.
Soil is known as a major reservoir of natural antibiotics and intrinsic antibiotic-resistant bacteria (ARB) (Allen et al. 2010). For instance, native Actinobacteria can produce a variety of antibiotics, such as beta-lactam, tetracycline, and glycopeptide, which thus promote the detection of associated antibiotic resistance genes (ARGs) (e.g., bla, tetM, and vanA) even in pristine environments (D'Costa et al. 2011). When manure is used as soil fertilizer, it may not only introduce extrinsic ARGs and ARB, but also affect the abundance of ARGs and ARB that are native to the soil. The rich organic matters in manure may promote the proliferation of both intrinsic and extrinsic ARB (Xie et al. 2018). At the same time, antibiotic residuals in manure can elicit selective pressure that enriches ARGs, promote their transmission, and even stimulate gene mutations to evolve new AMR mechanisms (Xie et al. 2018). A recent study revealed that tetracycline remained high at 342.49 mg/kg in chicken litters and long-term application of this contaminated manure for over two years significantly increased tetracycline resistance Communicated by Sunita Varjani. genes (e.g., tet(X)/(X2), tet(X3), and tet(X4)) by approximately 40-fold in crop soil as compared with the control (p < 0.01, 2-tailed t test) (He et al. 2021). Furthermore, longterm application of manures can increase the resilience of resistomes by shuffling extrinsic ARGs through horizontal gene transfer into native bacteria with better fitness and viability in the environment (Xie et al. 2018). Recent reports indicated both intrinsic (e.g., aacC and blaSFO) and extrinsic ARGs (e.g., InuB and ermC) can persist in soils for years after manure application is ceased without re-amendments Xie et al. 2016).
Maize (Zea mays) is one of the most important food crops in the world. The annual production of maize was 1137 million metric tons in 2019 (Erenstein et al. 2022). Maize is widely grown in 160 countries throughout the continents of America, Europe, Oceania, Africa, and Asia (Ranum et al. 2014). As a multipurpose crop, maize can be used for human food products, animal feed, and biofuel (Watson 2017). Maize rhizosphere can serve as "hotspots" of ARGs and ARB since root exudates, and nutrients can promote microbial proliferation and accelerate gene swapping through horizontal transfer (Berg et al. 2005;Chen et al. 2019). Using qPCR, a recent study demonstrated the use of manures from sulfadiazine-treated pigs as soil fertilizer significantly increased the abundance of sul1 genes in maize rhizosphere by approximately tenfold greater than the control (Chen et al. 2018). Li et al. found a significant elevation in the abundances of five β-lactamase genes and two tetracycline-and vancomycin-resistant genes in maize rhizosphere than those in bulk soil via DNA microarray (pair t test, p < 0.05) (Li et al. 2014). The use of these molecular techniques (e.g., qPCR and microarray) is limited by the availability of ARG primers, inability to link to their hosts, and other issues, such as false negatives caused by PCR inhibitors and false positives caused by nonspecific amplification (Li et al. 2015;Wind et al. 2021). Furthermore, there is a knowledge gap regarding the profiling of AMR derived from maize agriculture and its potential risks to human and animal health.
With the advance in high-throughput sequencing, metagenomics enables the profiling of a broad spectrum of ARGs, identification of their hosts, and discovery of novel ARGs, without sacrificing the quantitative power (Luby et al. 2016). In this study, we conducted a global survey of resistomes in maize rhizosphere and profiled the diversity and abundance of ARGs in 52 metagenomes of maize rhizosphere samples from North America, Europe, and Africa. We focused on ARGs that are of high abundance and/or frequency in maize rhizosphere and investigated their potential origins (i.e., intrinsic vs extrinsic) and transmission ability. We also assessed their correlations with host bacteria and soil management practice in the crop field. Certain ARGs and their hosts are highlighted given their potential risks of human exposure via maize product consumption and transmission to human pathogens. This in-depth investigation provides the molecular insights into ARG origins and dissemination potentials in maize rhizosphere, calling for sustainable management frameworks to control and mitigate the spread of AMR in agroecosystems.

Collection and assembling of global metagenomic datasets of maize rhizosphere samples
We identified 52 metagenomes of maize rhizosphere available in the databases of the Joint Genome Institute (JGI, https:// genome. jgi. doe. gov/ portal/) and the National Center for Biotechnology Information (NCBI) Sequence Read Archive (SRA) (https:// www. ncbi. nlm. nih. gov/ sra). All metagenomes were obtained by the Illumina sequencing-bysynthesis technology, resulting in the assembly-based dataset sizes of 0.4 -7.6 Gb. Sample and sequencing information, including the sampling time and location and sequencing reads and post-assembly dataset sizes, is summarized in Table S1. Soil management at the maize fields where samples were collected is outlined in Table S2. Among them, 27 metagenomes for maize rhizosphere samples collected in Michigan (n = 15) and California (n = 12) were assembled and downloaded directly from JGI. (Nordberg et al. 2014) For the other 25 metagenomes in South Africa (n = 21) and the Netherlands (n = 4), only raw sequencing reads were available from NCBI SRA (Sayers et al. 2020). Trimmed metagenomic reads were then assembled into contigs using SPAdes v3.13 (Prjibelski et al. 2020) with the metaSPAdes option and default parameter settings.

Bioinformatics analysis
ARG-carrying contigs (ACCs) were screened by searching against the ResFinder database (Bortolaia et al. 2020) and the Comprehensive Antibiotic Resistance Database (CARD) (Alcock et al. 2020) using Blastn with an identity of ≥ 70% and a query coverage of ≥ 60% within the Pan Resistome Analysis Pipeline (PRAP) (He et al. 2020). Redundant ARGs were manually removed in the combined data set. Then, sequences of ARGs were extracted from ACCs using the seqtk toolkit. ARGs were categorized into types and subtypes based on the antibiotics they confer resistance to and the associated molecular mechanisms, respectively. For example, all ARGs conferring resistance to beta-lactam are considered as a singular ARG type, whereas blaOXA is one of the beta-lactam resistance subtypes. For quantitative comparison among different metagenomics datasets, the relative abundance of ARGs was calculated by normalizing the number of reads mapped to the identified ARGs (based on average coverage) to million sequencing reads (rpm) in individual metagenomes (Yang et al. 2013).
To further investigate the potential mobility of ARGs, three online tools (i.e., ISfinder, ICE, and Integrall) were used to predict insertion sequences, integrative and conjugative elements, and integrons in ACCs, respectively. An E-value cut-off was set at 10 −5 for all three tools. For plasmid prediction, ACCs were searched against the NCBI plasmid database PLSDB (retrieved on March 4, 2021) with an E-value cut-off of 10 −5 and a minimum similarity of 95%.
For taxonomic assignments, ACCs were annotated against the NCBI non-redundant (NR) protein database (retrieved on March 10, 2021) using DIAMOND v.2.0.8 at an E-value cut-off of 10 −5 (Buchfink et al. 2015). A customized database of human and animal pathogenic bacteria, designated as the pathogenic database, was constructed to search for pathogens that contain the closest sequence homologues to ARGs obtained from metagenomic datasets of maize rhizosphere. The flow chart for identifying the closest human pathogens of ACCs from the maize rhizosphere metagenome is shown in Figure S1. Briefly, for each selected ARG subtype (e.g., dfrB) of high abundance (> 0.020 rpm) or prevalence (n ≥ 15), a list of accession numbers of human and animal pathogen-associated contigs was obtained from the human and pathogen databases on the CARD website and NCBI's Microbial Browser for Identification of Genetic and Genomic Elements (https:// www. ncbi. nlm. nih. gov/ patho gens/ micro bigge/#). Then, nucleotide sequences of these contigs were downloaded as FASTA files from the NCBI GenBank (http:// www. ncbi. nlm. nih. gov/) using the Entrez Programming Utilities (E-utilities) and used to construct the local pathogenic database. Nucleotide sequences of highly abundant ARG subtypes obtained from metagenomic datasets of maize rhizosphere were searched against the local pathogenic database using BLAST v.2.9.0 + at an E-value cut-off of 10 −5 . The BLAST hit outputs were subsequently filtered to identify the closest pathogens at a stringent nucleotide identity threshold of 99% (Forsberg et al. 2012). The ACCs carrying these filtered ARGs were further searched against the genomic data (i.e., GENBANK files) of the closest pathogens retrieved from NCBI GenBank to identify homologous flanking regions that embrace the ARGs.

Metagenome-assembled genome analysis
The average metagenome assembly size of maize rhizosphere samples from Michigan and California was 4.2 ± 1.2 Gb, which is sufficient to assemble genomes of high quality. However, samples from South Africa and the Netherlands have low assembly sizes (0.78 ± 0.59 Gb), precluding the recovery of genomes (Table S1). Considering the higher diversity of ARGs detected in the Michigan rhizosphere than that in the California rhizosphere sample, we focused on genomes recovered from the metagenomes in Michigan rhizospheric samples.
A total of 236 metagenome-assembled genomes (MAGs) from Michigan rhizosphere samples were downloaded from the JGI portal system. Among them, 159 MAGs were selected since they have a quality of ≥ 50% (estimated as "completeness − 5 × contamination"). (Zhao et al. 2020) ARGs present in MAGs were identified using the PRAP pipeline with the same parameters as aforementioned for ACC determination. The relative abundance of ARG-carrying MAGs was calculated based on the number of reads mapped to MAGs (based on the average coverage) normalized to the total number of sequencing reads. (Pasolli et al. 2019;Zhao et al. 2020) The host of individual ARG is determined by the phylogeny of its carrying MAG.

Network analysis
The ARG co-occurrence patterns in 17 Michigan samples were assessed using the Spearman's rank correlation with a rho coefficient ρ ≥ 0.6 and a significant FDR corrected p-value ≤ 0.01 (Ju and Zhang 2015). RStudio v1.2.5033 was used to perform the statistical and network analyses. Gephi v.0.9.2 was used for the network visualization (Bastian et al. 2009).

Statistical analysis
Principle coordinate analysis (PCoA) was processed to visualize dissimilarities among resistomes using the Bray-Curtis distance. Pairwise permutational multivariate analysis of variance (pairwise perMANOVA) was used to evaluate the significance of the difference between resistomes (Anderson 2014). FDR corrected p values according to the Benjamini-Hochberg Procedure (a False Discovery Rate multiple test correction) (Benjamini and Hochberg 1995) were estimated for pairwise perMANOVA tests.

Higher diversity and abundance of ARGs detected in Michigan maize rhizosphere
Across four studied regions, a total of 653 ACCs were detected in 38 out of 52 maize rhizosphere metagenomes. The total ARG abundance varied greatly from 0.09 ~ 4.50 rpm, and the detected ARGs can be categorized into 18 types and 68 subtypes ( Fig. 1a and b). Michigan maize rhizosphere metagenomes showed the greatest abundance and diversity of ARGs. As shown in Fig. 1a, the averaged relative abundance of ARGs in Michigan maize rhizosphere metagenomes was 1.52 rpm, 2.6 -7.1 times greater (p < 0.01, Mann-Whitney U test) than those in the other three regions (0.58 rpm for California, 0.26 rpm for the Netherlands, and 0.21 rpm for South Africa). Michigan maize rhizosphere metagenomes also showed a greater diversity of ARG subtypes than those in the other three regions (54 vs 4 -24 ARG subtypes) (Fig. 2). The most dominant ARG subtypes in Michigan maize rhizosphere confer resistance to beta-lactam, elfamycin, fosfomycin, fluoroquinolone, triclosan, and multidrug. These antibiotics are often used for animal growth promotion and disease prevention (Li et al. 2015). The extensive use of these antibiotics may entail selective pressure to enrich ARGs in animal fecal microbiomes and further on soil and rhizosphere microbial communities once manure is applied for fertilization.
Notably, blaOXA (an averaged relative abundance of 0.21 rpm, n = 9), blaACT (0.19 rpm, n = 6), and blaTEM (0.20 rpm, n = 3) were the top three abundant ARGs, together accounting for approximately 30% of the total ARGs in Michigan maize rhizosphere metagenomes on average (Fig. 2). These ARGs encode beta-lactamases that can inactivate a broad spectrum of commonly used betalactams, particularly those belonging to the cephalosporin and penam classes, by hydrolyzing the β-lactam ring (Bush 2018). Frequent detections of blaOXA and blaTEM have been reported in the phyllosphere (i.e., the above-ground Fig. 1 Abundance (a) and diversity (b) of ARGs in maize rhizospheres across the 4 different regions. Significant differences between Michigan and the other three regions were tested using two-way Mann-Whitney U test (**, p < 0.01)

Fig. 2
Ln-scale relative abundance of ARG subtypes in maize rhizospheres across the 4 different regions. Four putative intrinsic ARGs are highlighted in red boxes surface of plants) of a variety of vegetables, such as tomatoes, carrots, and cucumbers when they were fertilized with manure (Cerqueira et al. 2019;Knapp et al. 2010;Udikovic-Kolic et al. 2014). Another study also detected a high abundance of blaOXA and blaTEM (10 −4 ~ 10 −1 ARG copies per 16S rRNA gene copy) in almost all 12 samples collected from animal feces and wastewater from livestock farms (Li et al. 2015). BlaACT was found in greenhouse soils amended annually with commercial organic manure for 5 years . Our results concur with others, revealing the dominance of beta-lactam resistance genes in the maize rhizosphere and other environments, which become concerning as they may pose serious threats to public health. At lower abundance or frequency, blaBJP (0.022 rpm), blaJOHN (0.022 rpm), blaPEDO (0.043 rpm), EF-Tu (0.13 rpm), FosA2 (0.021 rpm), UhpA (0.022 rpm), emrR (0.021 rpm), and fabI (0.041 rpm) were also detected in Michigan maize rhizosphere metagenomes, but rarely in those from other regions (Fig. 2).
As depicted in the PCoA diagram (Fig. S2), resistomes of Michigan samples formed a distinct cluster from the other three regions along the primary coordinate PC1. This was supported by the pairwise perMANOVA test (FDR corrected p < 0.05) ( Table S3). The data above demonstrated a higher diversity and abundance of ARGs in the Michigan maize rhizosphere. This can also be supported by a higher proportion of antibiotic inactivating genes (see SI) and unique detection of blaTEM-116, FosA2, and blaACT-4/-6 that are of great relevance (≥ 99.0% identity) with animal and human pathogens (see discussion below). All Michigan samples were collected between 2014 and 2016 from the Kellogg Biological Station where manure or composted manure from animals treated with antibiotics, such as ceftiofur and oxytetracycline, were used for fertilization for over 10 years before 2008 (Guo et al. 2021;Hussain et al. 2020;Munir and Xagoraraki 2011). This previous manure application at the Kellogg Biological Station may be associated with the high abundance and diversity of ARGs in the Michigan maize rhizosphere metagenomes. Nonetheless, many other factors, such as soil properties (e.g., nutrients), environmental characteristics (e.g., meteorology), and farming practices (e.g., harvest time), differ greatly among maize fields from 4 regions we studied and they may also affect rhizospheric metagenomes and associated resistomes in addition to previous exposure to manure applications.

Prevalence of RbpA, vanRO, mtrA, and dfrB in global resistomes of maize rhizosphere
In the global resistomes of maize rhizosphere, RbpA and vanRO were the most abundant and prevalent ARGs, which were detected in 34 and 32 samples of the total 57 metagenomes from all four different regions, respectively (Fig. 2). The other two highly abundant ARGs, mtrA and dfrB, were identified in 30 and 27 samples from Michigan, California, and the Netherlands, but not in samples from South Africa. The relative abundances of RbpA, vanRO, mtrA, and dfrB were in the range of 0.02 -0.08 rpm, 0.09 -0.19 rpm, 0.01 -0.04 rpm, and 0.02 -0.07 rpm accounting toward 5.1 -14.4%, 12.8 -96.0%, 1.9 -12.0%, and 2.8 -25.1% of the total ARGs in different maize rhizosphere samples, respectively. RbpA is an actinobacterial transcription factor that binds directly to the principal sigma subunit of RNA polymerases promoting the transcription and precluding the interference of rifamycin, which thus confers rifamycin resistance (Newell et al. 2006). VanRO is a component of the vanO gene cluster that confers resistance to glycopeptide (Gudeta et al. 2014). MtrA is a transcriptional activator of the gene cluster MtrCDE that encodes a multidrug efflux pump conferring resistance to penam and macrolide (Rouquette et al. 1999). DfrB encodes plasmidborne type II dihydrofolate reductases, which can reduce the inhibitory efforts of trimethoprim in the DNA replication (Toulouse et al. 2017).
In addition, RbpA, vanRO, and mtrA have frequently been detected in environments without previous anthropogenic interruptions. For instance, a recent study detected RbpA from Actinobacteria in a biofilm community on the wall of a pristine cave in Western New Guinea (Turrini et al. 2020). VanRO has been detected widely in a variety of soils, including pristine Tibet soils and farm soils (0.01 ~ 0.22 vanRO copies/16S rRNA gene copies) (Li et al. 2020;Wind et al. 2021;Zhang et al. 2020). Furthermore, mtrA was observed in tropical agricultural soil which was used to cultivate mostly maize and vegetables and has no history of manure use or other known AMR sources for over 15 years (Salam et al. 2021).
A majority of bacterial hosts of RbpA, vanRO, and mtrA belonged to the phylum Actinobacteria that are likely native to maize rhizosphere (Li et al. 2014). We found more than 90.0% of RbpA and mtrA were carried by Mycobacterium and Mycolicibacterium. VanRO was also predominantly identified in members of Actinobacteria (94.6% in total), mostly Mycobacterium and Streptomyces. Figure S8 showed Mycobacterium and Mycolicibacterium were dominant resistant taxa detected in maize rhizosphere samples in Michigan, California, and the Netherlands.
As key members in soil, Mycobacterium spp. are known for their roles in nutrient transformation, such as autotrophic carbon fixation and phosphorus utilization (Deng et al. 2021;Li et al. 2014). In addition, Mycobacterium frequently co-exist with Streptomyces, which are the producers of a wide variety of commonly used antimicrobials, such as aminoglycoside, glycopeptide, and rifamycin (Peterson and Kaur 2018). The co-existence with Streptomyces is known to lead to the evolution of resistance mechanisms and intrinsic ARGs in Mycobacterium (Peterson and Kaur 2018).
Furthermore, mobile genetic elements were rarely detected in contigs that carry RbpA, vanRO, and mtrA, and these three ARGs were reported to be typically located on the chromosome Zhao et al. 2020). This information suggested their low mobility and transmission potential in the soil microbial community or beyond. Collectively, these findings indicated that RbpA, vanRO, and mtrA were likely intrinsic ARGs prevailing in maize rhizosphere and have a low potential of being transmitted to human pathogens.

Evolutionary analysis of dfrB and its clinical concern
DfrB was the other ARG widely detected in the maize rhizosphere metagenomes. We conducted the alignment of protein sequences encoded by dfrB genes detected in maize rhizosphere metagenomes with reference proteins available in CARD. According to CARD, dfrB gene family contains 8 reference genes, including dfrB1, dfrB2, dfrB3, dfrB4, dfrB5, dfrB6, dfrB7/dfrB8 (equivalent to each other), and dfrB9. Results in Figs. S3 and S4 showed several essential amino acid residues, including Lys32, Gln67, Ile68, and Tyr69, are highly conserved. These amino acid residues can be involved in ligand binding and catalysis for the dfrB family (Alonso and Gready 2006). Presence of conversed amino acids validated the annotation of these dfrB genes detected in the maize rhizosphere.
The prevalence of dfrB genes in maize rhizosphere metagenomes suggested they are probably intrinsic. In addition to trimethoprim inactivation, dfrB-coding dihydrofolate reductases were predicted to participate in certain unknown metabolic processes, since the occurrence of dfrB was found prior to the first clinical use of trimethoprim as antibiotics (Alonso and Gready 2006). This is in good agreement with a recent study that reported dfrB3 as a soil-intrinsic ARG since it was one of the most abundant ARGs (accounting for 30 to 90% of resistomes) in six different farm soils before manure application. In addition, this study also found the abundance of dfrB3 genes decreased by fourfold 4 days after the application of manure probably as their hosts were outcompeted by extrinsic ARB (Macedo et al. 2021).
Interestingly, the combination of phylogenetic analysis and ARG-flanking region analysis revealed a portion of these detected dfrB genes are potentially mobile and of clinical concern. Phylogenetic analysis based on their amino acid sequences revealed two distinct groups for dfrB genes in maize rhizosphere metagenomes (Fig. 3). Particularly, 17 out of 48 dfrB genes formed a distinctive clade (group I) from all known dfrB reference genes, while the remaining dfrB genes (group II) clustered with dfrB2, dfrB3, and dfrB9 reference genes. Many group II dfrB genes were clustered with those associated with human pathogens, such as dfrB3 (WP_172688537.1) in Pseudomonas aeruginosa (NG_047746.1), dfrB2 (WP_015060228.1) in Burkholderia dolosa (NZ_CM002277.1), and dfrB3 (WP_010890144) in Klebsiella pneumoniae (NZ_WUKZ01000186.1).
Furthermore, as highlighted in Fig. 3, 9 group II dfrB gene-carrying contigs also contain a recombination binding site (attC), which together constitute a mobile gene cassette. The attC sequences in these group II dfrB-carrying contigs showed high similarities (73.9 ~ 93.6% of nucleotide sequence identity) to two attC reference sequences that are located downstream of the integron class 1 (IntI1) in Pseudomonas aeruginosa strain 60,503 plasmid p60503-DIM (MN208063) and human pathogen Klebsiella pneumoniae (JADAYE010000036) (Fig. S5).
DfrB was frequently found in a mobile gene cassette with an attC recombination site located within integrons or exists as a free circular non-replicative molecule (Alonso and Gready 2006). Thus, dfrB was often detected in a broad spectrum of hosts within the phylum of Proteobacteria (e.g., Escherichia, Pseudomonas, and Staphylococcus) in human and animal samples, wastewater samples, and field soil samples regardless of manure amendment (Alonso and Gready 2006;Zhang et al. 2020). Therefore, it is likely that the group II dfrB genes from maize rhizosphere are highly mobile and transmissible as promoted by the adjacent MGEs. They can be of health concern considering the possibility of transferring the group II dfrB genes from maize rhizosphere to human pathogens through food production and consumption. The group I dfrB genes, on the other hand, show less mobility and horizontal gene transfer is less likely to occur.
It is worth noting that all three blaTEM-116 genes were embedded within transposons Tn3 located in the fragments with the length between 0.9 and 2.2 kbp. These blaTEM-116-carrying fragments share > 99.0% identity with the sequences detected in several pathogenic bacteria, such as Acinetobacter baumannii, Salmonella enterica, and Staphylococcus aureus as depicted by blue bars in Fig. 4. Similarly, ARGs that belong to the blaTEM family are often located within Tn3 transposons, promoting transmission of blaTEM and other ARGs to bacteria, some of which are urgent and serious pathogens, such as vancomycin-resistant Enterococci (VRE) and methicillin-resistant Staphylococcus aureus (MRSA) (Partridge et al. 2018).
Glutathione-S-transferase encoded by FosA2 can cleave the epoxide ring of fosfomycin, consequently reducing the susceptibility of the bacterial host to this antibiotic (Falagas et al. 2019). In the Michigan maize rhizosphere metagenomes, 3 contigs were found to carry FosA2 genes. These FosA2 genes showed > 99.5% similarity to the reference genes in human and animal pathogens, E. asburiae and E. ludwigii (Table S4). Furthermore, the FosA2 flanking regions (1.1 ~ 58.4 kbp) overlap with genomes of both pathogens with high identities of > 98.3%. For example, Figure S7 depicted that the contig Ga0207674_10000029 recovered from the Michigan maize rhizosphere metagenomes and E. asburiae shared a region of ~ 58.4 kbp, both of which encompass an MGE (i.e., a prophage integrase IntA) and virulent factors, including hcp (type VI secretion system tube protein) and tsr (methyl-accepting chemotaxis protein). Despite the discovery of FosA2 flanked with the integrase IntA in only one contig from maize rhizosphere, this ARG is frequently identified on plasmids harbored by the bacterial family Enterobacteriaceae (Falagas et al. 2019), suggesting their mobility in the environment matrices and clinical settings.
Two contigs were found to carry blaACT-4 and blaACT-6, showing ≥ 99.0% identity to their reference genes in pathogenic Enterobacter asburiae. For instance, as shown in Figure S6, the contig Ga0207674_10001472 is highly similar (99.2%) to the ~ 29.6 kbp-long flanking region in the genome of E. asburiae (LZCS01000012). This shared region also contained LysR (a family transcriptional regulator of betalactamase), sugE (a quaternary ammonium compoundresistance protein), and a virulent factor GroEL (a heat shock protein involved in intestinal colonization), forming a putative resistance island with multiple ARGs and genes that promote the dissemination of antimicrobial resistance.
Moreover, in this study, we also investigated the host bacteria for these three ARGs. Acinetobacter spp. were Fig. 3 Unrooted phylogenetic tree of dfrB gene sequences identified across the maize rhizospheres from 4 regions and reference genes selected from NCBI and CARD. The phylogenetic tree was constructed using a maximum likelihood algorithm with 500 bootstrap replicates assigned as the hosts of bla-TEM116, while Enterobacter spp. were identified as hosts of blaACT-4, blaACT-6, and FosA2. Members of Acinetobacter and Enterobacter were reported associated with manure application. A previous study reported that long-term manure amendment for 5 years can enrich Acinetobacter and Enterobacter in soil samples while they were absent in none amended controls, suggesting these species were transmitted from manure to soil samples . It is likely that the long-term application of manure at Michigan fields can facilitate the possible transmission of extrinsic ARGs (e.g.,blaTEM116,and FosA2) and their hosts (i.e., Acinetobacter and Enterobacter species) from manure to maize rhizosphere microbiomes. Collectively, the high similarity of blaACT-4/-6, blaTEM-116, and FosA2 with those found in human and animal pathogens (≥ 99.0%) and the presence of MGEs in ARG-flanking regions implied that they were probably extrinsic ARGs which might be transferred from manure. However, the potential enrichment of indigenous bacteria, such as Acinetobacter and Enterobacter carrying those ARGs by rich organic sources provided from manure cannot be ruled out.

Multidrug-resistant bacteria in Michigan maize rhizosphere
Among the 159 MAGs recovered from Michigan maize rhizosphere metagenomes, 9 drug-resistant MAGs were identified (Table S5). Among them, MAG 3300013105_29 and 3300026116_28 annotated as Klebsiella planticola and Enterobacter sp., respectively, contain multiple ARGs. In particular, the MAG of Klebsiella planticola harbors 7 ARGs, conferring resistance to elfamycin (EF-Tu), fosfomycin (Ptsl), triclosan (fabI), and multidrug (CRP, H-NS, marA, and soxS). The Enterobacter MAG carries 15 ARGs encoding resistance to beta-lactam (blaACT-6), fluoroquinolone (parE), fosfomycin (FosA2, ptsI, and UhpA), triclosan (fabI and gyrA), and multidrug (acrA, cpxA, CRP, emrR, HN-S, marA, oqxB, and ramA). It is known that many Enterobacter species, such as E. cloacae and E. asburiae species, are multidrug-resistant and responsible for causing nosocomial infections, spanning from respiratory infections to bacteremia (Gou et al. 2020). In this regard, we screened the presence of ARGs in the genome database of E. cloacae and E. asburiae that were isolated from soils and plants. All eight genomes of E. cloacae and E. asburiae available on JGI revealed the presence of four ARGs, including blaACT for beta-lactam, rpoB for rifamycin, acrB for multidrug, and parC for fluoroquinolone. Note that acrA and acrB, a component of acrA-acrB-TolC multidrug efflux system (Burse et al. 2004), were also identified in the Enterobacter MAG recovered from Michigan rhizosphere metagenomes. Moreover, Enterobacter species are known to contribute significantly to the dissemination of ARGs, such as blaTEM, blaCTX, and blaKPC genes (Chavda et al. 2016;Davin-Regli et al. 2019). Hence, extra attention is needed to be drawn to Enterobacter and other species that carry multiple ARGs, since they can adapt to the selection pressure of antibiotics and promote the acquisition and dissemination of ARGs in the maize rhizosphere.

Co-occurrence patterns of ARGs in Michigan maize rhizosphere metagenomes
The network analysis revealed co-occurrence patterns of ARGs in Michigan maize rhizosphere metagenomes. The network shown in Fig. 5 consists of 17 nodes (ARG subtypes) and 50 correlation links (average degree of 5.88 and average shortest path length of 1.67). Two ARG clusters I and II were distinctively formed (Fig. 5, Table S6). The multidrug resistance genes (acrA, adeR, cpxA, CRP, H-NS, marA, and ramA), which accounted for 81.6% of the total ARGs in Cluster I, were connected with beta-lactam (blaACT ), elfamycin (EF-Tu), fluoroquinolone (emrR and parE), fosfomycin (FosA2, PtsI, and UhpA), and/or triclosan (fabI) resistance genes. The other cluster (II) Fig. 4 Comparison between three blaTEM-116-carrying contigs in Michigan maize rhizosphere metagenomes and those identified in the closest human pathogenic isolates. The blue bars depict genetic fragments of high similarity (> 99%), and the red lines indicate Tn3 fragments consisted of a connection of only the beta-lactam resistance gene (blaOXA) and the multidrug resistance gene (soxS). Interestingly, in both clusters, ARGs from different ARG types co-existed more commonly than those from the same ARG types.
Network analysis also revealed multidrug resistance efflux pumps often co-occurred with the other resistant mechanisms in the maize rhizosphere. Multidrug-resistance efflux pumps confer resistance to not only different types of antibiotics but also natural substances produced by the host plants, thus enhancing environmental adaptation and interaction of bacterial species with the plants (Piddock 2006). For example, acrA detected in Cluster I is the periplasmic portion of the acrA-acrB-TolC multidrug efflux complex, which confers resistance to antibiotics including cephalosporin, fluoroquinolone, glycylcycline, penam, phenicol antibiotic, rifamycin, tetracycline, and triclosan, as well as plant-derived antibacterial toxins, such as flavonoids, isoprenoids, and alkaloids (Burse et al. 2004). Therefore, hosts of arcA genes can survive and thrive under diverse environmental stresses as reported by Burse et al.'s study. Escherichia coli KAM3 transformants carrying acrA-acrB genes from Erwinia amylovora exhibited eightfold greater resistance (MIC of 500 µg/L vs MIC of 62.5 µg/L) to phloretin (i.e., an antibacterial agent produced by apples) than the control (Burse et al. 2004). Moreover, the marA and ramA tend to co-exist with acrA as all three genes are grouped in Cluster I. The co-existence of these three ARGs is plausible since marA and ramA can act as positive regulators of acrA-acrB and thus upregulate the activity of the efflux pump system in response to appropriate antibiotics (Holden and Webber 2020).

Limitations and environmental implications
Overall, we investigated a broad spectrum of antibiotic resistance in maize rhizosphere from 4 regions in 3 continents across the world, including Michigan, California, South Africa and the Netherlands. Since this study is based on metagenomic sequence information from publicly available databases, it is limited by the fact that many variables spanning from maize cultivation, sample preparation, to sequencing depth cannot be well controlled. However, following standardized data processing and normalization protocols, it is feasible to extract ARGs and their carrying fragments from different metagenomes for the phylogenetic and diversity comparison. With the increasing availability of relevant metagenomes in maize rhizosphere through variable-controlled environments and intercontinental sampling, more conclusive correlations can be drawn between ARG profiling and agricultural practices and other impacting factors.
Taken together, we found 4 ARGs, including vanRO, RbpA, mtrA, and dfrB, prevalently present across all maize rhizosphere samples. Among them, vanRO, RbpA, and mtrA are likely to stem from native Actinobacteria, posing low clinical concern. In contrast, dfrB can be troublesome given its transmission potential to some human pathogens, such as Pseudomonas aeruginosa and Klebsiella pneumoniae. To gain a better understanding of the evolution and transmission of dfrB genes, future studies are needed to investigate genomic and physiological characteristics of dfrB-carrying isolates from soil and plant rhizosphere, as well as the horizontal gene transfer potentials of this concerning ARG.
Furthermore, regions like South Africa and the Netherlands showed low detection of a similar collection of ARGs, which are mostly intrinsic based on our investigation. In contrast, Michigan samples exhibited a much higher diversity and abundance of ARGs, many of which were of human/ animal origins and/or located adjacent to MGEs. These concerning ARGs may be associated with the long history of manure fertilization, despite that many factors (e.g., meteorology and soil properties) and agricultural practices can affect the resistomes in maize rhizosphere. The presence of these resistance determinants can increase the likelihood of their entry into the crop products and subsequently promote the spreading of their resistance genes to human pathogens through direct contact with the agricultural soils or consumption of contaminated maize products.

Fig. 5
The co-occurrence network of ARG subtypes in Michigan maize rhizosphere metagenomes. A connection represents a significant correlation (Spearman's rank correlation with rho coefficient ρ ≥ 0.6 and significant FDR corrected p-value ≤ 0.01). The nodes were colored based on ARG types and weighted according to their relative abundances (rpm). The thickness of the connecting lines was proportional to a correlation coefficient