ARTICLE https://doi.org/10.1038/s41467-022-29283-8 OPEN Assessment of global health risk of antibiotic resistance genes Zhenyan Zhang1,6, Qi Zhang1,6, Tingzhang Wang2,6, Nuohan Xu1, Tao Lu 1, Wenjie Hong2, Josep Penuelas 3,4, Michael Gillings 5, Meixia Wang2, Wenwen Gao2 & Haifeng Qian 1 ✉ 1234567890():,; Antibiotic resistance genes (ARGs) have accelerated microbial threats to human health in the last decade. Many genes can confer resistance, but evaluating the relative health risks of ARGs is complex. Factors such as the abundance, propensity for lateral transmission and ability of ARGs to be expressed in pathogens are all important. Here, an analysis at the metagenomic level from various habitats (6 types of habitats, 4572 samples) detects 2561 ARGs that collectively conferred resistance to 24 classes of antibiotics. We quantitatively evaluate the health risk to humans, defined as the risk that ARGs will confound the clinical treatment for pathogens, of these 2561 ARGs by integrating human accessibility, mobility, pathogenicity and clinical availability. Our results demonstrate that 23.78% of the ARGs pose a health risk, especially those which confer multidrug resistance. We also calculate the antibiotic resistance risks of all samples in four main habitats, and with machine learning, successfully map the antibiotic resistance threats in global marine habitats with over 75% accuracy. Our novel method for quantitatively surveilling the health risk of ARGs will help to manage one of the most important threats to human and animal health. 1 College of Environment, Zhejiang University of Technology, 310032 Hangzhou, P. R. China. 2 Key Laboratory of Microbial Technology and Bioinformatics of Zhejiang Province, Hangzhou 310012, P. R. China. 3 CSIC, Global Ecology Unit CREAF-CSIC-UAB, Bellaterra, 08193 Barcelona, Catalonia, Spain. 4 CREAF, Campus Universitat Autònoma de Barcelona, Cerdanyola del Vallès, 08193 Barcelona, Catalonia, Spain. 5 ARC Centre of Excellence in Synthetic Biology, Faculty of Science and Engineering, Macquarie University, Sydney, NSW 2109, Australia. 6These authors contributed equally: Zhenyan Zhang, Qi Zhang, Tingzhang Wang. ✉email: hfqian@zjut.edu.cn NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 A ntibiotic resistance is an increasing global threat to human the Comprehensive Antibiotic Research Database (CARD). Of health and to the clinical treatment of disease1. Antibiotic these, 2401 were genes conferring resistance to only one drug resistance genes (ARGs) have been detected in the last class, and 160 conferred resistances to multiple drug classes decade in all environments, including natural2–4, engineered5–8, (Supplementary Data 2). Twenty-five ARGs were found in more and clinical9,10 habitats. Anthropogenic activities, including the than 75% samples, however, the frequency of most ARGs (2313/ clinical use of antibiotics, are widely regarded as the main drivers 2561) were <10% (Supplementary Fig. 1a). On the other hand, of the dissemination of ARGs6,11,12. nearly half of 2561 ARGs were commonly shared by diverse ARGs, however, do not arise from current human activities. habitats (Supplementary Fig. 1b), especially genes conferring They existed prior to the antibiotic era, having been detected in resistance to widely used antibiotics26 like aminoglycosides, tet- permafrost13 and human paleofeces14. Anthropogenic activities racyclines, and beta-lactams (Supplementary Fig. 1c). These instead drive the selection of genes from environmental and results implied that anthropogenic activities, like the use of cellular sources, and these can be subsequently co-opted to confer antibiotics, were critical for dissemination of ARGs globally. antibiotic resistance. These genes originally had a range of We examined the abundance and composition of ARGs in environmental functions, such as liberating phosphorus from diverse sub-habitats at a global scale. The human-associated phosphate15, or encoding efflux pumps. Diverse genes can habitats, including the digestive system and skin, had the highest potentially confer resistance to antibiotics, so we need to deter- abundances of ARGs (Fig. 1b). Built environments, mainly mine whether genes bioinformatically identified as ARGs pose a including urban subways, also had considerable abundances of health risk to human health, and to do so from multiple ARGs, confirming these as hotspots of ARGs7. Genes conferring perspectives16–18. resistance to tetracyclines and aminoglycosides, two widely used The first consideration is the potential for the transmission of antibiotics in the clinic26, dominated in digestive system and skin, ARGs from the environment to bacteria in humans16, defined respectively, while ARGs conferring multidrug resistance had a here as “human accessibility”. ARGs can transfer from the high proportion in built environments (Fig. 1c). Although environment to humans, via their bacterial hosts, and then have a geographic factors like latitude were reported to influence the negative impact on human health6,19. Mass sampling and the abundance of the ARGs4, confirmed in this study (Supplementary development of metagenome-assembled genomes (MAGs) have Fig. 2), anthropogenic activities are also critical for dissemination allowed us to understand the distribution of ARGs and their hosts of ARGs. To further determine the impacts of anthropogenic in specific habitats4–7,9. Danko et al. (2021) assembled the first activities on the dissemination and abundance of ARGs, we atlas of urban metagenomics using 4728 metagenomic samples calculated population densities at each sample site. from mass-transit systems in 60 cities7. Understanding the dis- tribution and dissemination of ARGs in global habitats from the Effect of anthropogenic activities on the dissemination of perspective of global health is important, especially from envir- ARGs. Sampling sites were clustered into two groups according to onmental compartments to humans20. The second consideration their general population density, one with high-intensity activity is the potential for the transmission of ARGs from environmental (>58 people/km2, the global average population density27), and bacteria to pathogens16–18, here defined as “mobility” and the other with low population density. Regions of high-intensity “human pathogenicity” of the ARGs. From this perspective, only activity had significantly higher total abundances of ARGs and ARGs in pathogenic hosts that can infect humans pose an ele- genes conferring resistance to specific classes (p < 0.001, two- vated risk to human health16. ARGs frequently move by hor- tailed Welch’s t-test; Fig. 1d and Supplementary Fig 3). A total of izontal gene transfer (HGT) from nonpathogens to pathogens, 671 ARGs were specifically detected in environments with high- and this activity has been a major driver in the evolution of intensity activity (Fig. 1e). For ARGs shared between high- and pathogens resistant to antibiotics12,21,22. The third perspective is low-intensity environments, 715 were significantly more abun- the current clinical use of antibiotics17, here defined as the dant in high-intensity environments (p < 0.05, two-tailed Welch’s “clinical availability” of the ARGs. The use of antibiotics has t-test; Fig. 1f) and most of these conferred beta-lactam and increased globally in recent years23, and some new antibiotics multidrug resistance (Supplementary Fig. 4). have been developed for clinical use24,25. In contrast, some Thirty-four ARGs were specific to low-intensity environments, antibiotics are now rarely used. Considering clinical relevance in mainly annotated as beta-lactam resistance genes (adjust p > 0.05, evaluating the risk of ARGs is thus necessary. two-tailed Welch’s t-test; Fig. 1e and Supplementary Fig. 4). The In this work, based on the three perspectives outlined above, distributions of 1102 ARGs were not significantly influenced by we use 4572 metagenomic datasets to reveal the distribution and anthropogenic activities (adjust p > 0.05, two-tailed Welch’s t-test; dissemination of 2561 ARGs as well as their hosts in global Supplementary Data 3). These genes are probably not resistance habitats. In the next step, we carry out a framework to quanti- genes that will affect human health and likely perform different tatively evaluate the health risk of each ARG and sample by functions for their bacterial hosts in the natural environment13,15,17. considering the four indicators (human accessibility, mobility, We identified potential ecological functions of ARGs in biogeochem- human pathogenicity, and clinical availability). There are 23.78% ical cycling by annotating and mapping all ARG-like reads with genes of the ARGs pose a health risk, especially those which confer associated with the cycling of carbon, nitrogen, phosphorus and multidrug resistance. Finally, we predict and map the antibiotic sulfur (Supplementary Fig. 5). There were 43 genes initially identified resistance threats in marine habitats using machine learning. as ARGs that clearly perform biological functions in addition to antibiotic resistance. Results indicated that different ARGs exhibit the different level of correlation to the anthropogenic activities, which Results will influence the health risk of ARGs on human lives. We then Global patterns of ARG distribution. We used a set of 4572 quantitatively evaluated the health risk of each ARGs considering the metagenomic samples to illustrate the global patterns of ARG four indicators (human accessibility, mobility, human pathogenicity, distribution (Supplementary Data 1). These samples were col- and clinical availability) in the following sections. lected from six types of habitats: air, aquatic, terrestrial, engi- neered, humans and other hosts (Fig. 1a and Supplementary Data 1). From these samples, we identified a total of 2561 ARGs Human accessibility of ARGs. We first examined the ARGs that that conferred resistance to 24 drug classes of antibiotics based on were shared between humans and the other three main habitats to 2 NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 ARTICLE a b c Sub-habitat Human Digestive system Skin Engineered Built Lab enrichment Wastewater Bioreactor Aquatic Freshwater Marine Thermal springs Non-marine Saline and Alkaline Terrestrial Soil Deep subsurface Habitat Air (9) Human (1020) Other host (276) 2000 1500 1000 500 0 50 100 Aquatic (1697) Terrestrial (396) Engineered (1174) Abundance (RPKM/sample) Relative abundance of ARG classes (%) ARG classes Multidrug e Beta-lactams d g p=7.60e−168 Aminoglycoside 20000 Tetracycline High Glycopeptide Low 34 1846 671 Terrestrial Peptide Diaminopyrimidine Phenicol 15000 Low High 1193 Fluoroquinolone f 8.00 MLS Abundance (RPKM) Elfamycin Log 2 (Average abundance) 6.00 Fosfomycin 10000 1460 Human Others 4.00 715 genes 1223 2.00 5000 0.00 3.50 2.50 1.50 29 genes 0 0.00 Aquatic Group Low High Built Fig. 1 Distribution patterns of antibiotic resistance genes (ARGs) globally. a Geographic distribution of samples with ARG abundance in various habitats. Each point indicates one sampling location, rounded to the nearest degree, with point size reflecting the number of samples, and point color indicating the habitat. Samples with the unclassified location was showed as longitude = 0 and latitude = 0. b Abundance of ARGs in each sub-habitat. Digestive system of human, mainly including the fecal samples, had the highest abundances of ARGs. c Composition of antibiotic resistome in each sub-habitat. Only sub- habitats containing at least 20 samples are shown. d High-intensity human activities significantly promoted the abundance of ARGs. Each dot represents one sample (n = 1643 and 2309 samples for Low- and High-group, respectively). The p value represents the statistical significance (two-tailed Welch’s t- test). e Number of ARGs specific or shared in the areas with low- or high-intensity human activities. There were 671 ARGs specifically detected in high- intensity human activities environment. f The abundance of 715 and 29 ARGs significantly increased in high- and low-intensity human activities environment, respectively (adjust p < 0.05, two-tailed Welch’s t-test; Supplementary Data 3). g ARGs shared between the human-associated and three main habitats. Number in the circles represents the number of shared ARGs. investigate the accessibility of ARGs to the human microbiota. As Distribution of ARG hosts and mobile genetic elements expected, built environments had the most ARGs shared with (MGEs) in different habitats. From the 4572 metagenomic human habitats (1460), with terrestrial (1193), and aquatic (1223) samples we used, 18,465 metagenome-assembled genomes environments having fewer ARGs (Fig. 1g). Most of these ARGs (MAGs) have been recovered by Nayfach et al.28, which greatly were annotated as conferring resistance to multidrug and beta- aided our work in determining the global distribution of ARG lactams. We then determined the average abundance and pre- hosts. However, directly predicting the host of ARGs based on valence of each ARG in the human habitats and calculated human MAGs is challenging and could lead to misleading results because accessibility (see the “Methods” section; Supplementary Data 4). MAGs are composite genomes and do not represent an actual Only 1714 of 2561 ARGs were detected in the human habitats, genome of a single microbe in the community29. Here, we tried to most of them with average abundances <50 reads per kilobase per improve the accuracy of host identification of ARGs by imple- million (RPKM) per sample and prevalence <10% (Supplemen- menting strict quality criteria. We only considered ARGs in tary Fig. 6 and Supplementary Data 4). The gene tetQ, conferring contigs longer than 10 kb and made sure that the taxonomic resistance to tetracycline antibiotics, had the highest human affiliation of any genes found in those ARG-containing contigs accessibility (Supplementary Fig. 6). These results indicated that agreed with the overall taxonomy of the MAG (Fig. 2a). Subse- the human accessibility of ARGs were variable and only a fraction quently, 7555 MAGs were identified as being host of ARGs of ARGs exhibited high accessibility to humans and posed (Supplementary Data 5). The hosts of individual ARGs differed potential risk. significantly in different habitats (Fig. 2b and Supplementary NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 a ARG-contig 3 b c abundance (Bacteria A, > 4 10 kb) ARG-contig G-c-con 4 (Bacteria ria A, < 10 kb) Aquatic Terrestrial 2 kb) 0 ARG-contig 2 (Bacteria A, > 10 kb) MAG 5 > 10 r2=0.2011, p<2.2e-16 ARG (Bacte G-contig -2 (Bacteria A) 3 richness a C, ria AR 2 Built Human 1 ARG-ccontig kb) , > 10 1 (Bacteria B Other Pyhlum Proteobacteria Alpha- Beta- Delta- Gamma- Other- r2=0.1892, p<2.2e-16 0 Bacteroidetes Firmicutes Verrucomicrobia Chlorobi Acidobacteria 2 3 4 5 MGE d 5000 e f g 6000 Single 6000 Non-pathogen p=3.32e-11 Mutiple Non-pathogen Pathogen Abundance of MGE per sample Pathogen p<2.2e-16 p<2.2e-16 Number of Genome 4000 Number of host 4000 464 81 57 4000 3000 Low (13.4%) 2000 2000 2000 1000 High 449 281 232 (29.2%) 0 0 0 Low High Low High Low High Density h 27013 completed genome i k 10124 600 16889 Non-pathogen Pathogen 400 Mobility j IS Integrase 200 Transposase <5kb <5kb Integrase/Transposase MGE1 MGE2 0 ARG 0 1000 2000 3000 4000 0.00 0.25 0.50 0.75 1.00Density Number of MGEs Human pathogenicity Fig. 2 Human pathogenicity and mobility of ARGs. a The host of ARGs were determined by only considering ARGs in contigs longer than 10 kb and making sure that the taxonomic affiliation of those ARG-containing contigs agree with the overall taxonomy of the MAG. b Composition of antibiotic resistance host at the phylum level (and classes in Proteobacteria) showed the different distributions of ARG hosts in four main habitats. c MGEs exhibited a significantly positive correlation with ARGs in abundance and richness. Each point represents one sample (n = 4572 samples). The value of abundance and number of ARGs were showed in log10. Result of linear regression are shown as r2 and p value, evaluated by F-statistic (one-sided). d Abundance of MGEs in high-intensity human activities areas were significantly higher than in low-intensity human activities areas. n = 1643 and 2309 biologically independent samples for Low- and High-group, respectively. Significance between two groups was evaluated by two-tailed Welch’s t-test. Data were showed as mean ± standard error of mean (SEM). e Number of genomes containing multiple ARGs increased in high-intensity human activities areas, compared to the low-intensity human activities areas (Fisher’s exact test, one-sided). f Shared ARGs in pathogens or non-pathogen increased in high-intensity human activities areas, compared to low. The percentage of shared ARGs in the two areas are shown. g Pathogenic hosts of ARGs increased in high-intensity human activities areas, compared to low (Fisher’s exact test, one-sided). h We collected 27,013 completed genome from NCBI RefSeq database for determining the human pathogenicity and mobility of ARGs. Five kb upstream and downstream of the ARGs detected in all completed genome for annotating the MGEs. i Among the 27,013 completed genome, 16,889 were recognized as pathogens’ genome. j We totally identified 4612 MGEs from completed genomes, and most of them were referred to the transposase. k The human pathogenicity of ARGs exhibited a bimodal distribution, while most of ARGs (2266/2561) showed mobility <10. 4 NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 ARTICLE Fig. 7a). Hosts of most ARGs (89.61%) were specific to one method determined a potential mobility of ARGs, which was habitat, with only few hosts shared across two or three habitats critical for risk assessment. Most of the evaluated ARGs (2265/ (Supplementary Fig. 7b). ARG hosts in the built and human- 2561) were carried by less than 10 different MGEs (Mobility <10) associated habitats were less diverse than the hosts in more nat- (Fig. 2k). ural habitats (Fig. 2b), perhaps a consequence of selective pres- We also determined the potential human pathogenicity of an sures from anthropogenic activities. ARG based on the proportion of pathogens that carried it to The habitats contained distinct ARG hosts, suggesting that they evaluate the health risk of clinical ARGs (see the “Methods” directly or indirectly selected these hosts. Mobile genetic elements section; Supplementary Data 8). The human pathogenicity of (MGEs), including insertion sequences (ISs), integrons and ARGs exhibited a clearly bimodal distribution (Fig. 2k), indicating transposons, are all capable of horizontally transferring ARGs that most ARGs had specific genetic contexts and were exclusively from one bacterium to another in association with plasmids and carried by non-pathogens only (human pathogenicity = 0) or phages21,30. Abundance and distribution of MGEs can be shaped pathogens only (human pathogenicity = 1). by environmental selection31–33. Thus, we determined the composition of MGEs in each sample from the metagenomic reads to discover if they influence the distribution of ARG hosts Health risk assessment of ARGs and samples. In summary, we in diverse habitats. MGEs were significantly and positively analyzed the characteristics of ARGs to determine their human correlated with ARGs in abundance and richness, respectively accessibility, mobility, and human pathogenicity based on their (linear regression: p < 0.001) (Fig. 2c). Few MGEs (3.80%) were potential to move from the environment to humans and to drive shared by the four habitats, similar to the distribution of hosts the evolution of pathogens resistant to antibiotics. We further carried ARG (Supplementary Fig. 8). The unique environmental determined the clinical relevance of ARGs by systematically conditions in each habitat type could influence the distribution of evaluating the risk of ARGs to human health. Data for the global MGEs31–33 and thus determine the distinct compositions of the use of antibiotics were collected37, which indicated that penam hosts that contained ARGs in each habitat. The abundance of (55.25%) and cephalosporins (13.07%), two beta-lactam anti- MGEs per sample was significantly higher in areas with high-, as biotics, were used the most (Supplementary Fig. 9). The total use opposed to low-intensity anthropogenic activities (p < 0.001, two- of antibiotics for each ARG was calculated as clinical availability tailed Welch’s t-test; Fig. 2d). This could lead to the insertion of (see the “Methods” section and Supplementary Data 9 for details). more ARGs into genomes (Fig. 2e), which then drives an increase Genes conferring resistance to clinically available antibiotics were in ARGs shared between non-pathogens and pathogens (Fig. 2f) a high proportion of the 2561 ARGs we detected. The multidrug while also increasing the range of potentially pathogenic hosts of resistance gene tolC, for example, can confer resistance to nearly ARGs (Fig. 2g). These results indicate that MGEs likely all common antibiotics and thus had the highest clinical avail- contributed to the HGT of ARGs between hosts, and significantly ability (Supplementary Data 9). from non-pathogens to pathogens. We evaluated the overall health risk for each ARG using the four calculated metrics: human accessibility, mobility, human pathogenicity and clinical availability, as determined above Mobility and human pathogenicity of ARGs. Considering that (Fig. 3a). All, except clinical availability, only covered about half MGEs can be either not binned, or incorrectly binned into the of the ARGs (Fig. 3b). We calculated the risk index (RI) as wrong MAGs34,35, we collected 27,013 completed genomes in RI = HA × MO × HP × CA (also see the “Methods” section). This NCBI RefSeq database36 to determine the mobility and human formula was quite strict, and only 23.78% of the 2561 ARGs with pathogenicity of ARGs (see the “Methods” section; Fig. 2h). All all four indicators were predicted as a health risk (Fig. 3c; these completed genomes were sequenced by whole-genome Supplementary Data 10). Most high-risk ARGs were multidrug sequencing, the accurate and standard approach for discovering resistance determinants (Fig. 3d). However, this formula is still MGEs35, and 16,889 of them were recognized as pathogens reasonable. For example, the ARGs with high clinical availability (Fig. 2i and Supplementary Data 6). For determining the mobility but with no MGEs may only be genes intrinsic to specific of ARGs, we extracted 5 kb upstream and downstream of the hosts17,38 and cannot transfer between hosts. We divided the ARGs detected in all completed genome for annotating the MGEs ARGs with an RI > 0 into four categories based on the rank of (Fig. 2h). We only considered the ISs, integrases and transposases their RI for further examination: Q1 (top 25%), Q2 (50–75%), Q3 in this step. Some sequences attributed to the function of plas- (25–50%), and Q4 (bottom 25%). Interestingly, most ARGs mids and phages, but that did not directly affect the mobility of conferring multidrug resistance or resistance to commonly used ARGs were excluded. We used such conservative approach antibiotics, such as tetracycline belonged to Q1, and ARGs mainly because it is difficult to identify phages and plasmids at conferring resistance to rarely used antibiotics such as glycopep- the gene level. Some genetic elements close to ARGs may be tides belonged to Q4 (Fig. 3e). These results confirmed that the involved in the function of plasmids and phages, however, they use of antibiotics increased the risk of ARGs to human health and cannot contribute to the HGT of ARGs and result in false potentially caused the failure of clinical treatments of infection38. positives21. In total, there were 4612 MGEs identified from To validate the performance of our assessment method, we completed genomes, and most of them (3061/4612) were trans- determined ARGs in 568 hospital pathogenic MAGs from posase (Fig. 2j). another catalog of human gut microbiota39 (Supplementary It is now clearer than ever that MGEs were greatly responsible Data 11) as well as a subset of completed genomes from to the dissemination of ARGs and used for determining the pathogens. The method used for assigning ARGs in MAGs was mobility of ARGs in the previous studies, which assessed the the same as Fig. 2a. Results clearly showed that the number of health risk of ARGs qualitatively17,18. In the present study, for ARGs belonging to Q1 were significantly higher than other ranks quantitative analysis, the mobility of ARGs was defined as the per genome (Fig. 3f), confirming that these ARGs were highly number of associated MGEs detected (see “Methods” section; dangerous to human health and complicated the clinical Supplementary Data 7). It should be noted that it is almost treatment of disease. On the contrary, ARGs belonging to Q3 impossible to measure the absolute value of the mobility of ARG, and Q4 were seldomly carried by genomes of pathogens. These which can be changed with the genetic contexts in specific results testified to the validity of our work in health risk species, because of the fitness costs in HGT17. However, our assessment. NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 23 d Average risk index .7 e a b 3000 8% 1e+6 1e+4 1e+2 1e+0 ARG classes c >0 109 Multidrug Number of ARGs, excluding zero 2560 Human accessibility 0 36 MLS 45 Tetracycline Mobility 2000 Risk index 1714 37 Fluoroquinolone 79 Aminoglycoside 1231 172 Beta-lactams Health Risk on Human 80.00 1000 907 25 Peptide Proportion (%) 60.00 12 Fosfomycin Human pathogenicity 40.00 23 Phenicol 22 Diaminopyrimidine Clinic availability 0 20.00 HA MO HP CA 27 Glycopeptide 0.00 Q1 Q2 Q3 Q4 f g 4 10 Habitat MAG Risk rank Completed genome Risk rank Number of ARGs per genome Aquatic a a Q1 8 Q1 3 Terrestrial Q2 Q2 Built Q3 6 Q3 Human 2 Q4 Q4 4 Average RI b 1e+9 1 b 2 2e+9 c d c e 0 0 3e+9 Fig. 3 Health risk evaluation of ARGs. a We evaluated the health risk of human for each ARG with four indicators, including human accessibility (HA), mobility (MO), human pathogenicity (HP), and clinical availability (CA). b Number of ARGs after excluding the zero value in each indicator for risk index calculation. c Only 23.78% of all evaluated ARGs exhibited risk index (RI) > 0. RI of most ARGs was zero because of the strict formula we used. d The number and average RI of ARGs in each class. Only classes that had more than 10 ARGs with RI > 0 were shown. Most of the ARGs with RI > 0 were referred to multidrug resistance with the highest average RI. e Composition of ARGs with RI > 0 in different classes. Most of the ARGs which resisted commonly used antibiotics belonged to Q1, while most of ARGs which resist rarely used antibiotics belonged to Q4. f Number of ARGs per pathogenic genome, which assembled by hospital fecal metagenomes (n = 568 MAGs) or from the completed genome dataset (n = 15,596 MAGs). Pathogens carried much more ARGs belonged to the Q1, compared with other ranks. Different letters represented the significant difference by Kruskal–Wallis H-test with the pairwise comparisons. Data were showed as mean ± standard error of mean (SEM). g Global ARG risk map of four main habitats. Antibiotic resistance risk was detected all around the world, even in the polar region. Human-associated habitats posed the highest risk of antibiotic resistance than other habitats. The average RI of each sampling site was calculated by the combination of abundance and RI of ARGs and showed as the size of points. Habitats were showed as color. Samples with the unclassified location was showed as longitude = 0 and latitude = 0. We evaluated the risk of antibiotic resistance for each sample in machine learning for better performance and avoiding the using a combination of abundance and RI of ARGs (see the overfitting of prediction model43. “Methods” section) to model global surveillance, based on a data The different methods discretized the dataset in different results set of 4005 metagenomic samples from four main habitats (Supplementary Figs. 12 and 13). Equal frequency resulted in a well- (aquatic, terrestrial, built, and human-associated) around the distributed dataset, however, it failed to clearly distinguish the world. The risk of antibiotic resistance was global, even in the samples in ranks 1–5. On the contrary, equal width clearly polar region (Fig. 3g). Human-associated habitats posed the differentiated the samples in each rank, but nearly all the samples highest risk of antibiotic resistance, and health risk was also were grouped as rank 1 with only one sample in some ranks. k- influenced by anthropogenic activities (Supplementary Fig. S10), means algorithm, the most known and used clustering method40, as expected. balanced the sample number (not strictly even but better than the original dataset and equal width) and dissimilarity in each rank. We further determined the prediction performance of machine learning Global mapping of the antibiotic resistance threats in marine based on the dataset discretized by k-means and equal frequency, and environments. After evaluating the health risk of each sample, we the former exhibited much higher accuracy rate than the latter wanted to map the antibiotic resistance threats around the world, (Fig. 4b). Results of 10-fold cross-validation also showed the using machine learning (Fig. 4a). We used 712 samples from prediction performance of machine learning with the k-means marine habitats to establish the prediction model (Supplementary method (Supplementary Fig. 14). We then chose the best model in Data 12). The distribution of the risks for marine samples were 10-fold cross-validation with the accuracy rate = 76.06% for further uneven (Supplementary Fig. 11), so for better prediction accu- analysis. This model classified each risk rank well (Fig. 4c). Latitude, racy, the dataset was discretized by three unsupervised methods which has been confirmed to significantly influence the abundance of (k-means40, equal width, and equal frequency41), and the samples ARGs (Supplementary Fig. 2), was the most important predictor in were then divided into 10 ranks according to their risks (rank 10 the prediction model (Fig. 4d). In the meantime, climate change for the highest risk and rank 1 for the lowest risk). Seventeen stressors including ultraviolet radiation changes (UV), sea level rise anthropogenic drivers in marine habitats provided by Halpern (SLR), surface temperature rise (STR), and ocean acidification (OA) et al.42 as well as the latitude were chosen as 18 factors (Sup- were also critical for the prediction model. plementary Data 12) that influenced the risk ranks. The random With the prediction results from machine learning, global forest algorithm combined with 10-fold cross-validation was used mapping of antibiotic resistance risk in marine habitats with each 6 NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 ARTICLE a b d 0.8 Accuracy rate 76.06% Latitude Discretization UV 0.7 k-means SLR equal width k-means equal frequency DNL equal frequency 0.6 OP 1 2 3 4 5 6 7 8 9 10 Marine sample 10 Risk ranks 18 Factors Round STR (712) c 1.0 PL DNH 0.8 OA 1 2 10 Shipping True positive rate SLN 0.6 PH R1 AUC = 0.91 IP All data R2 AUC = 0.88 R3 AUC = 0.91 ... ... 0.4 R4 AUC = 0.91 FP R5 AUC = 0.94 PP R6 AUC = 0.98 Training 0.2 R7 AUC = 0.99 OR Machine learning R8 AUC = 0.99 Fishing Test R9 AUC = 0.99 (Random Forest) Average AUC = 0.96 R10 AUC = 1.00 IS 0.0 0.0 0.2 0.4 0.6 0.8 1.0 0.00 0.05 0.10 0.15 10-fold cross-validation False positive rate Importance e ARG rank 10 9 8 7 6 5 4 3 2 1 Fig. 4 Global mapping of antibiotic resistance threats in marine habitats. a Machine learning were trained by 712 samples from marine habitats and used to predict the antibiotic resistance threats in global marine habitats. b Accuracy rate of machine learning with different discretization methods. K-means exhibited higher accuracy rate than equal frequency and the best model (accuracy rate = 76.06%) were chosen for the further prediction. c The ROC plots confirmed the high performance of the best model in classification of risk ranks. d Latitude as well as the climate change stressors exhibited the high importance in predicting the antibiotic resistance risk. The full name of each indicator can be found in Supplementary Data 12. e The map of ARG risk in marine habitats with prediction results by machine learning, which was drawn by ArcGIS in 20′ × 20′ resolution. pixel (20′ × 20′ resolution) were done (Fig. 4e). Note that the risk marine areas, with 712 samples. Thus, for more comprehensive rank here only represents the relative risk (antibiotic resistance and accurate mapping of antibiotic resistance threats in marine risk in rank 10 areas is higher than rank 1), areas in rank 1 did environments, we need more high-quality metagenomic datasets. not mean no risk. The global distribution of antibiotic resistance On the other hand, we could not map the antibiotic resistance is related to the prediction of the abundance of ARGs in each threats in terrestrial and human-associated habitats, mainly country and territory by Hendriksen et al.5. For example, the because of the uneven distribution of samples and the lack of ARG risks in the marine areas close to Brazil and Africa were metadata on original soil properties. higher than in the marine areas close to eastern USA. This result was consistent with the higher ARG abundance of urban sewage Discussion in Brazil and Africa, compared to the USA5. For comparison We annotated 4572 metagenome samples with CARD to identify between oceans, the antibiotic resistance threats in Pacific and ARGs and obtain an unprecedented view of their global ARG Atlantic Ocean were the higher than others. The antibiotic distribution. A total of 2561 ARGs covering various habitats were resistance risk of marine areas close to the Antarctic Pole were used to link their abundance to the intensity of anthropogenic higher than the areas near the Arctic Pole. This map provided an activities. The dissemination of ARGs across habitats and HGT of overview of the antibiotic resistance threat in global marine, ARGs from non-pathogenic hosts to pathogens increased with which was important for surveillance. However, it is still limited the intensity of anthropogenic activities, associated with an by the number of samples, for example, samples in high-risk increase in the abundance of MGEs. These results were pre- ranks were less than lower ones, and there was low coverage of dictable, but the global distribution of the main hosts of ARGs NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 was surprising. By examining different habitats, we demonstrated of ARGs. We downloaded 4572 metagenomic samples for determining the global that particular ARGs were found in distinct host taxa in different distribution of ARGs and MGEs across diverse habitats from European Nucleotide Archive49. We fortunately have the permission to download the 18,465 MAGs environments. Sequencing the 16S rRNA gene combined with constructed by Nayfach et al. 28 from the IMG/M portal50,51. For more accurately high-throughput qPCR is a common approach to identify determining the mobility and human pathogenicity of ARGs, 27,013 completed microorganisms that are ARG hosts by correlation analysis2,3,44. genomes in NCBI RefSeq database36 were collected. To validate the performance of Metagenomic sequencing can more accurately identify hosts that our assessment method, we determined ARGs in 568 hospital pathogenic MAGs of have ARGs, although studies are often limited to a single another catalog39 from European Nucleotide Archive49. The pathogenic genomes in this study were defined by comparison of their taxonomical information to the habitat4,7 and metagenomic approaches can cause false positives. A-to-Z database which is continually updated with clinically relevant pathogens Here we examined all available MAGs from multiple habitats by and aligned with information from the Robert Koch Institute in Germany and the applying stricter quality criteria, which improved the accuracy of WHO52. Detailed information for these datasets is provided in Supplementary host identification of ARGs. Data 1, 5, 6 and 11, including the ARGs annotation, sampling location, habitat information and pathogenic identification. IBM® Aspera Connect (v4.1.1) was used The global distribution of ARGs and their hosts, mainly based for downloading all these data. on their number and abundance is important, but insufficient for evaluating the health risk of antibiotic resistance16–18. By sys- tematically considering their accessibility, potential to contribute ARG and MGE annotation and abundance calculation at the metagenomic level. The raw data of metagenomic samples were qualified by FastQC (v0.11.5; to pathogenicity, mobility between hosts and relevance to clinical https://github.com/s-andrews/FastQC), and then trimmed and quality-filtered treatment, we improved determinations of the health risk for the using Trimmomatic53 (v0.36). ARGs were annotated with the CARD using reads 2561 ARGs detected in the metagenomic analysis. Only 23.78% of by their recommended tool, RGI54 (v5.1.1), with default parameters for metage- the 2561 ARGs were determined as “risk ARGs” with an RI > 0 nomic reads. BWA55 (v0.7.13) was used for mapping reads to ARGs in each because of our stringent, but reasonable, methods. We confirmed sample, and the unmapped reads were removed using Samtools56 (v1.3.1). The number of mapped reads of ARGs in each sample were counted using a script, that the use of antibiotics was tightly linked to the health risk of which is available at GitHub (see “Code availability” below). ARG abundance was the ARGs, but some ARGs were exceptions. For example, ARGs calculated as RPKM with the number of reads and gene lengths. ARGs were conferring resistance to beta-lactams were common, but only 172 manually reclassified based on the drugs to which they confer resistance. ARGs of 1631 ARGs in this class had a significant risk index. ARGs referring to penam, cephalosporin, carbapenem, cephamycin, penem and mono- bactam were grouped into the beta-lactam class. ARGs referring to macrolides, conferring resistance to glycopeptide antibiotics (such as vanco- lincosamides and streptogramins were grouped into the MLS class. ARGs referring mycin, a last line of defense antibiotic) were mainly in the lowest to more than one drug class were grouped into the multidrug class. risk Q4 category. For these ARGs, we need information about The annotation of MGEs in each sample from the metagenomic reads was whether they co-exist with multidrug-resistant ARGs in the same performed to discover if they could influence the distribution of ARG hosts in diverse habitats. Reads of ISs were annotated with the ISfinder database57 using pathogens30 before the target antibiotics are used. Vaccines would BWA55 (v0.7.13). Reads were also annotated with the nucleotide sequence referred then be more efficient and less risky than antibiotics45. to the NCBI Reference Sequence of transposases and integrases, which were Global prediction or surveillance of antibiotic resistance has clustered using CD-HIT58 (v4.7) with the threshold ≥90%, using BWA55 (v0.7.13). already begun5,46,47. Hendriksen et al. precisely predicted the The abundances of the MGEs were calculated in the same manner as for the ARGs. global abundance of ARGs using World Bank variables, mainly concerning sanitation and health5, but abundance did not directly Potential ecological functions of ARGs in biogeochemical cycling. We compiled represent health risk. In this study, we provided a feasible a comprehensive database of genes involved in the cycling of carbon, nitrogen, quantitative method for global mapping of antibiotic resistance phosphorus, and sulfur (CNPS), including SCycDB59, NCycDB60, and QMEC61. threats for each pixel (20′ × 20′) in marine habitats with machine Potential ecological functions of ARG-like reads were annotated with the CNPS database we created by BWA55 (v0.7.13). learning. However, there is still an immediate need for a global cooperative system for reporting ARGs and their hosts, especially in clinics. The Global Antimicrobial Resistance Surveillance Determination of the ARG hosts in MAGs. The MAGs and ARG contigs were System of the WHO provides the opportunity for global coop- taxonomically assigned using Kraken2 (v2.1.2) with the default parameter. After that, we tried to improve the accuracy of host identification of ARGs by imple- eration on surveilling antimicrobial resistance47. Most impor- menting stricter quality criteria: only considering ARGs in contigs longer than tantly, these data should be collected under the supervision of 10 kb and making sure that the taxonomic affiliation of any genes found in those governments and shared on public platforms such as the World ARG contigs agreed with the overall taxonomy of the MAG (Fig. 2a). Health Organization (WHO) with detailed metadata, including clear coordinates of sampling sites, physicochemical properties of ARGs and MGEs annotation in completed genomes. ARGs in completed gen- samples (for water or soils), age, gender, race and dietary habits omes were annotated with CARD by their recommended tool, RGI54 (v5.1.1), with (for humans), and so on. default parameters for genomes. We then extracted 5 kb upstream and downstream Our study provides a novel and straightforward method for of the ARGs detected in all completed genomes for annotating the MGEs, because quantification and standardization of the antibiotic resistance such close proximity of MGEs and ARGs is more likely to induce HGT21. For genomes, ISs were annotated with the ISfinder database57 using BLASTN (v0.7.13; threats all around the world that will improve decision-making in e-value ≤ 10−10, identity ≥ 80%, coverage ≥ 80%), while the transposases and inte- clinics and public health management. Predicting future threats, grases were annotated with the NCBI Reference Sequence, which were clustered as demonstrated by COVID-19, is difficult. The most serious using CD-HIT58 (v4.7) with the threshold ≥ 90%, at an e-value ≤ 10−10 with a threats, however, are likely those we already recognize48. ARGs minimum amino acid identity of 60% over 60% query coverage using Diamond62 (v0.9.36.137). and their hosts are unquestionably a serious global threat, cur- rently causing about 700,000 deaths a year around the world, and this could increase to 10 million deaths a year by 205048. We used Collection of population data and global antibiotic use. The population density a systematic and targeted surveillance of ARGs to evaluate the of each sampling location was collected from the SEDAC Population Estimation antibiotic resistance threats around the world, which provides Service63, which provides population data in a specific area using longitude and latitude. We only considered samples from aquatic, terrestrial, engineered and some early warning. Comprehensive systems for surveilling the human-associated habitats with clear location coordinates. All samples were risk of antibiotic resistance under global cooperation should be a clustered by population density into two groups regardless of habitat: one was an priority. area with high-intensity anthropogenic activities with >58 people/km2, and the other area was summarized as low-intensity activity. This threshold was chosen because it is the global average population density calculated according to the Methods global population and land area27. Data for global antibiotic use were collected Collecting the datasets of metagenome, completed genomes, and MAGs. In from an online website37 containing data on antibiotic use in hospitals and from this study, we used several datasets for comprehensively evaluating the health risk retail outlets in 76 countries. 8 NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 ARTICLE Determining the health risk of ARGs to humans. We emphasized the impor- Statistical analysis and visualization. Significant differences were identified using tance of determining the health risk of ARGs to humans from multiple perspectives multiple methods, mentioned in the main text and legends of figures. Two-tailed rather than basing risk only on abundance. We defined four indicators from these Welch’s t-tests were performed in Excel Analysis Tools (Microsoft Corporation, perspectives: “human accessibility” (HA), “mobility” (MO), “human pathogenicity” Redmond, USA), F-statistic and Fisher’s exact test were performed by the basic (HP), and “clinical availability” (CA). statistic package in R (v3.6.3), and Kruskal–Wallis H-test with the pairwise com- HA represents the ability of ARGs to transfer from the environment to the parisons were performed in IBM® SPSS® Statistics (v20.0.0). For multiple com- bacterial groups in humans16 and is calculated as parisons, p values were adjusted by FDR in R (v3.6.3). The geographic distribution of the samples and the global RI map of diverse habitats was constructed using the HA ¼ Average abundancehuman ´ Prevalencehuman ð1Þ ggmap65 and ggplot266 packages in R (v3.6.3). All Venn diagrams were generated where Average abundancehuman and Prevalencehuman represent the average using EVenn (http://www.ehbio.com/test/venn; v1.0). The regression analysis, abundance and prevalence of ARGs in human-associated habitats, respectively. point plots, density plots and box plots visualized using the ggplot2 package in R MO represents the ability of ARGs to transfer between hosts by HGT. We (v3.6.3). All the schematic diagrams or elements in this paper were drawn using determined the MO of ARGs as the number of the related MGEs in completed BioRender (https://app.biorender.com/; v1.0) with full publishing rights. Heatmaps genomes, which could potentially transfer them between genomes or to plasmids. were constructed using TBtools67 (v1.082). Other plots were constructed using HP specifically represents the ability of ARGs to transfer from nonpathogenic GraphPad Prism (v7.00). hosts to pathogenic hosts, leading to the evolution of pathogens resistant to antibiotics and the failure to control clinical infections17. HP was calculated for Reporting summary. Further information on research design is available in the Nature each ARG as Research Reporting Summary linked to this article. HP ¼ Numberpathogenic =Numberall ð2Þ where Numberpathogenic is the number of pathogenic ARG hosts and Numberall is Data availability the total number of hosts containing ARGs. All raw data used in this study are available in NCBI RefSeq database, IMG/M portal, and CA represents the clinical availability of ARGs, so we indicated that their health European Nucleotide Archive. Information for all metadata used in this study as well as risk was higher than others if the ARGs conferred resistance to the most commonly the important data for analysis are provided in Supplementary Data 1, 5, 6 and 11, used antibiotics. We then calculated CA for each ARG as including the accession-codes. Supplementary Data contained the critical supplementary information in this study are also publicly available online (https://doi.org/10.6084/ CA ¼ ∑ni¼1 Usage of antibiotici ð3Þ m9.figshare.19189652). The raw data underlying figures are provided as Source data where n is the number of antibiotics to which the ARG conferred resistance. For which can be obtained in public repository (https://doi.org/10.6084/ example, according to the classification in CARD, tetC confers resistance to m9.figshare.19189568). The nucleotide sequences in our CNPS database can be available macrolide, fluoroquinolone, aminoglycoside, carbapenem, cephalosporin in https://doi.org/10.6084/m9.figshare.19193489. Source data are provided with (cephamycin), glycylcycline, penam, tetracycline, peptide, aminocoumarin, this paper. rifamycin, phenicol, triclosan, and penem antibiotics. So, the number of antibiotics to which tetC conferred resistance is 14, and the clinical availability of tetC was the sum of the consumption of these 14 antibiotics. Code availability All the scripts and codes for gene annotation, abundance calculation, machine learning, Finally, the risk index (RI) of the ARGs to human health was calculated as statistical analysis and visualization used in this study are available online at https:// RI ¼ HA ´ MO ´ HP ´ CA ð4Þ github.com/ZhenyanZhang/ARG-global68. We also calculated RI for each sample as RIsample ¼ ∑ni¼1 Abundancei ´ RIi ð5Þ Received: 30 September 2021; Accepted: 2 March 2022; where Abundancei is the abundance of ARGi in the sample, and RIi is the RI of ARGi calculated using Eq. (4). To validate the performance of our assessment method, the ARG annotation from hospital pathogenic MAGs were performed as described above. We calculated References the number of ARGs carried by per pathogenic MAG with different risk ranks 1. Allen, H. K. et al. Call of the wild: antibiotic resistance genes in natural (Q1–Q4), respectively. We also compared the number of ARGs in each risk rank environments. Nat. Rev. Microbiol. 8, 251–259 (2010). detected in the pathogenic completed genomes. 2. Zhang, Q. et al. Cyanobacterial blooms contribute to the diversity of antibiotic-resistance genes in aquatic ecosystems. Commun. Biol. 3, 737 Global mapping for antibiotic resistance threats in marine environments. (2020). Seventeen data sets of anthropogenic drivers in marine habitats were collected from 3. Zhu, Y.-G. et al. Continental-scale pollution of estuaries with antibiotic the database provided by Halpern et al.42 and downloaded from https:// resistance genes. Nat. Microbiol. 2, 16270 (2017). knb.ecoinformatics.org/view/doi:10.5063/F1S180FS in August 2021 and used for 4. Bahram, M. et al. Structure and function of the global topsoil microbiome. prediction the antibiotic resistance threats in marine habitats. In the meantime, Nature 560, 233–237 (2018). latitude was chosen as a geographic factor because it has been reported to influence 5. Hendriksen, R. S. et al. Global monitoring of antimicrobial resistance based on the abundance of the ARGs4, as confirmed in this study. The 712 samples from metagenomics analyses of urban sewage. Nat. Commun. 10, 1124 (2019). marine habitats were divided into 10 ranks according to their risk after dis- 6. Kang, K. et al. The environmental exposures and inner- and intercity traffic cretization. The 18 selected indicators associated with the 712 marine samples were flows of the metro system may contribute to the skin microbiome and collected for machine learning using random forest algorithm combined with the resistome. Cell Rep. 24, 1190–1202.e1195 (2018). ten-fold cross-validation. Random forest was an accurate algorithm by using 7. Danko, D. et al. A global metagenomic map of urban microbiomes and bootstrap sample: each tree was built by about 2/3 samples of all data and model antimicrobial resistance. Cell 184, 3376–3393.e3317 (2021). performance was validated by the remaining out-of-bag data64. To further ensure 8. Pärnänen, K. M. M. et al. Antibiotic resistance in European wastewater the performance of random forest and avoid overfitting, we used 10-fold cross- treatment plants mirrors the pattern of clinical antibiotic resistance validation. The original dataset was randomly partitioned into 10 folds. In each prevalence. Sci. Adv. 5, eaau9124 (2019). round, nine folds were used to train the model as the training set, and the 9. Hu, Y. et al. Metagenome-wide analysis of antibiotic resistance genes in a large remaining data evaluated model as the test set43. The results of ten-fold cross- cohort of human gut microbiota. Nat. Commun. 4, 2151 (2013). validation were evaluated by confusion matrix. The final prediction model was 10. Chng, K. R. et al. Cartography of opportunistic pathogens and antibiotic constructed based on the k-means discretization methods because it exhibited the resistance genes in a tertiary hospital environment. Nat. Med. 26, 941–951 best performance. We further evaluated the performance of machine learning by (2020). receiver operator characteristic curve (ROC) plots, which constructed in each rank. The importance of each indicator was also determined by machine learning to 11. Bougnom, B. P. & Piddock, L. J. V. Wastewater for urban agriculture: A estimate the most critical factors influencing the antibiotic resistance threats in significant factor in dissemination of antibiotic resistance. Environ. Sci. marine habitats. We also collected the information of 18 indicators in 380,887 sites Technol. 51, 5863–5864 (2017). with 20′ × 20′ resolution in marine areas for prediction. After calculating risk for 12. Groussin, M. et al. Elevated rates of horizontal gene transfer in the each site, we converted all data to the Mollweide projection with a WGS84 datum industrialized human microbiome. Cell 184, 2053–2067.e2018 (2021). for showing the global antibiotic resistance threats in marine habitats by ArcGIS 13. D’Costa, V. M. et al. Antibiotic resistance is ancient. Nature 477, 457–461 (v10.8). All the scripts for machine learning, including building the prediction (2011). model, plotting ROC, determining the importance of indicators, and predicting the 14. Wibowo, M. C. et al. Reconstruction of ancient microbial genomes from the risk of unobserved marine sites, can be found in the GitHub (see the “Code human gut. Nature 594, 234–239 (2021). availability” section). Machine learning were performed in Python (v3.9) using 15. McRose, D. L. & Newman, D. K. Redox-active antibiotics enhance phosphorus PyCharm Community Edition (version: 2021.2.2). bioavailability. Science 371, 1033 (2021). NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 16. Manaia, C. M. Assessing the risk of antibiotic resistance transmission from the 47. World Health Organization. Global Antimicrobial Resistance Surveillance environment to humans: non-direct proportionality between abundance and System—Manual for Early Implementation https://www.who.int/publications/ risk. Trends Microbiol. 25, 173–181 (2017). i/item/9789241549400 (2017). 17. Martínez, J. L., Coque, T. M. & Baquero, F. What is a resistance gene? Ranking 48. May, M. Tomorrow’s biggest microbial threats. Nat. Med. 27, 358–359 (2021). risk in resistomes. Nat. Rev. Microbiol. 13, 116–123 (2015). 49. European Nucleotide Archive. http://www.ebi.ac.uk/ena EMBL-EBI (2021) 18. Zhang, A.-N. et al. An omics-based framework for assessing the health risk of 50. Chen, I. M. A. et al. The IMG/M data management and analysis system v.6.0: antimicrobial resistance genes. Nat. Commun. 12, 4765 (2021). new tools and advanced capabilities. Nucleic Acids Res. 49, D751–D763 19. Chen, Q.-L., Cui, H.-L., Su, J.-Q., Penuelas, J. & Zhu, Y.-G. Antibiotic (2021). resistomes in plant microbiomes. Trends Plant Sci. 24, 530–541 (2019). 51. Mukherjee, S. et al. Genomes OnLine Database (GOLD) v.8: overview and 20. Hernando-Amado, S., Coque, T. M., Baquero, F. & Martínez, J. L. Defining updates. Nucleic Acids Res. 49, D723–D733 (2021). and combating antibiotic resistance from One Health and Global Health 52. HARTMANN Science Centre. Relevant Pathogens from A–Z https:// perspectives. Nat. Microbiol. 4, 1432–1442 (2019). www.bode-science-center.com/center/relevant-pathogens-from-a-z.html 21. Ellabaan, M. M. H., Munck, C., Porse, A., Imamovic, L. & Sommer, M. O. A. (2021). Forecasting the dissemination of antibiotic resistance genes across bacterial 53. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for genomes. Nat. Commun. 12, 2435 (2021). Illumina sequence data. Bioinformatics 30, 2114–2120 (2014). 22. Jiang, X. et al. Dissemination of antibiotic resistance genes from antibiotic 54. Jia, B. et al. CARD 2017: expansion and model-centric curation of the producers to pathogens. Nat. Commun. 8, 15784 (2017). comprehensive antibiotic resistance database. Nucleic Acids Res. 45, 23. Klein, E. Y. et al. Global increase and geographic convergence in antibiotic D566–D573 (2017). consumption between 2000 and 2015. Proc. Natl Acad. Sci. USA 115, E3463 55. Li, H. & Durbin, R. Fast and accurate short read alignment with (2018). Burrows–Wheeler transform. Bioinformatics 25, 1754–1760 (2009). 24. Batchelder, H. R. et al. Development of a penem antibiotic against 56. Li, H. et al. The sequence Alignment/Map format and SAMtools. Mycobacteroides abscessus. Commun. Biol. 3, 741 (2020). Bioinformatics 25, 2078–2079 (2009). 25. Imai, Y. et al. A new antibiotic selectively kills Gram-negative pathogens. 57. Siguier, P., Perochon, J., Lestrade, L., Mahillon, J. & Chandler, M. ISfinder: the Nature 576, 459–464 (2019). reference centre for bacterial insertion sequences. Nucleic acids Res. 34, 26. World Health Organization. WHO Report on Surveillance of Antibiotic D32–D36 (2006). Consumption: 2016–2018 Early Implementationhttps://www.who.int/ 58. Fu, L., Niu, B., Zhu, Z., Wu, S. & Li, W. CD-HIT: accelerated for clustering the publications/i/item/9789241514880 (2018). next-generation sequencing data. Bioinformatics 28, 3150–3152 (2012). 27. Wikipedia. List of Countries and Dependencies by Population Density https:// 59. Yu, X. et al. SCycDB: a curated functional gene database for metagenomic en.wikipedia.org/wiki/List_of_countries_and_dependencies_by_population_ profiling of sulphur cycling pathways. Mol. Ecol. Resour. 21, 924–940 (2021). density (2021). 60. Tu, Q., Lin, L., Cheng, L., Deng, Y. & He, Z. NCycDB: a curated integrative 28. Nayfach, S. et al. A genomic catalog of Earth’s microbiomes. Nat. Biotechnol. database for fast and accurate metagenomic profiling of nitrogen cycling 39, 499–509 (2021). genes. Bioinformatics 35, 1040–1048 (2019). 29. Rice, E. W., Wang, P., Smith, A. L. & Stadler, L. B. Determining hosts of 61. Zheng, B., Zhu, Y., Sardans, J., Peñuelas, J. & Su, J. QMEC: a tool for high- antibiotic resistance genes: a review of methodological advances. Environ. Sci. throughput quantitative assessment of microbial functional potential in C, N, Technol. Lett. 7, 282–291 (2020). P, and S biogeochemical cycling. Sci. China Life Sci. 61, 1451–1462 (2018). 30. Botelho, J. & Schulenburg, H. The role of integrative and conjugative elements 62. Buchfink, B., Reuter, K. & Drost, H.-G. Sensitive protein alignments at tree-of- in antibiotic resistance evolution. Trends Microbiol. 29, 8–18 (2021). life scale using DIAMOND. Nat. Methods 18, 366–368 (2021). 31. Shaw, L. P. et al. Niche and local geography shape the pangenome of wastewater- 63. SEDAC. Population Estimation Service https://sedac.ciesin.columbia.edu/ and livestock-associated Enterobacteriaceae. Sci. Adv. 7, eabe3868 (2021). mapping/popest/pes-v3/ (2021). 32. Bethke, J. H. et al. Environmental and genetic determinants of plasmid 64. Bylander, T. Estimating generalization error on two-class datasets using out- mobility in pathogenic Escherichia coli. Sci. Adv. 6, eaax3173 (2020). of-bag estimates. Mach. Learn. 48, 287–297 (2002). 33. Brito, I. L. et al. Mobile genes in the human microbiome are structured from 65. Kahle, D. & Wickham, H. ggmap: spatial visualization with ggplot2. R. J. 5, global to individual scales. Nature 535, 435–439 (2016). 144–161 (2013). 34. Maguire, F. et al. Metagenome-assembled genome binning methods with short 66. Villanueva, R. A. M. & Chen, Z. J. ggplot2: elegant graphics for data analysis reads disproportionately fail for plasmids and genomic Islands. Microb. (2nd ed.). Meas.: Interdiscip. Res. Perspect. 17, 160–167 (2019). Genom. 6, mgen000436 (2020). 67. Chen, C. et al. TBtools: an integrative toolkit developed for interactive analyses 35. Brito, I. L. Examining horizontal gene transfer in microbial communities. Nat. of big biological data. Mol. Plant 13, 1194–1202 (2020). Rev. Microbiol. 19, 442–453 (2021). 68. Zhang, Z. Y. et al. Assessment of global health risk of antibiotic resistance 36. NCBI-FTP-Server. ftp://ftp.ncbi.nlm.nih.gov/blast/db/16SMicrobial.tar.gz. genes. Github https://doi.org/10.5281/zenodo.6125695 (2022). NCBI (2021) 37. The Center for Disease, Dynamics Economics & Policy. ResistanceMap: Antibiotic Resistance https://resistancemap.cddep.org/AntibioticUse.php Acknowledgements (2021). We appreciate Nayfach et al. (2020)28 for their work on the genomic catalog of Earth’s 38. Bungau, S., Tit, D. M., Behl, T., Aleya, L. & Zaha, D. C. Aspects of excessive microbiomes and its publication as a public resource, and Almeida (2021)40 for their antibiotic consumption and environmental influences correlated with the work on the human gut genome and its publication as a public resource. We also occurrence of resistance to antimicrobial agents. Curr. Opin. Environ. Sci. acknowledge all the PIs for uploading the metagenomic samples as precious public Health 19, 100224 (2021). resources. All the PIs’ name were listed in Supplementary Data 1. H.F.Q. was financially 39. Almeida, A. et al. A unified catalog of 204,938 reference genomes from the supported by the National Natural Science Foundation of China (21976161 and human gut microbiome. Nat. Biotechnol. 39, 105–114 (2021). 21777144). T.L. was financially supported by the National Natural Science Foundation of 40. Sinaga, K. P. & Yang, M. Unsupervised K-Means clustering algorithm. IEEE China (41907210). Josep Penuelas was financially supported by the Fundación Ramon Access 8, 80716–80727 (2020). Areces grant ELEMENTAL-CLIMATE, the Spanish Government grant PID2019- 41. Helskyaho, H., Yu, J. & Yu, K. in Machine Learning for Oracle Database 110521GB-I00, and the Catalan Government grants AGAUR-2020PANDE00117 and Professionals: Deploying Model-Driven Applications and Automation Pipelines SGR 2017-1005. (eds Helskyaho, H., Yu, J. & Yu, K.) 39–95 (Apress, 2021). 42. Halpern, B. S. et al. Spatial and temporal changes in cumulative human impacts on the world’s ocean. Nat. Commun. 6, 7615 (2015). Author contributions 43. Ban, Z. et al. Machine learning predicts the functional composition of the Z.Y.Z. and Q.Z. designed the study with guidance from H.F.Q. Z.Y.Z. wrote the first draft protein corona and the cellular recognition of nanoparticles. Proc. Natl Acad. of the manuscript and H.F.Q., J.P., and M.G. contributed substantially to revisions. Sci. USA 117, 10492 (2020). W.J.H., M.X.W., and W.W.G. contributed to the all functional annotation of 44. Wang, H. et al. The antibiotic resistome of swine manure is significantly metagenome-assemble genomes. T.Z.W., Z.Y.Z., and Q.Z. performed all the metage- altered by association with the Musca domestica larvae gut microbiome. ISME nomic analysis. N.H.X. and T.L. were responsible for the Machine learning model J. 11, 100–111 (2017). construction and related data analysis. Z.Y.Z. and Q.Z. performed the visualization of all 45. Lewnard, J. A., Lo, N. C., Arinaminpathy, N., Frost, I. & Laxminarayan, R. data and the artistic design of all figures. H.F.Q., T.L., and J.P. acquired funding for this Childhood vaccines and antibiotic use in low- and middle-income countries. project. Nature 581, 94–99 (2020). 46. Argimón, S. et al. Integrating whole-genome sequencing within the National Antimicrobial Resistance Surveillance Program in the Philippines. Nat. Competing interests Commun. 11, 2719 (2020). The authors declare no competing interests. 10 NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-022-29283-8 ARTICLE Additional information Open Access This article is licensed under a Creative Commons Supplementary information The online version contains supplementary material Attribution 4.0 International License, which permits use, sharing, available at https://doi.org/10.1038/s41467-022-29283-8. adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Correspondence and requests for materials should be addressed to Haifeng Qian. Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless Peer review information Nature Communications thanks the anonymous reviewer(s) for indicated otherwise in a credit line to the material. If material is not included in the their contribution to the peer review of this work. Peer reviewer reports are available. article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from Reprints and permission information is available at http://www.nature.com/reprints the copyright holder. To view a copy of this license, visit http://creativecommons.org/ licenses/by/4.0/. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. © The Author(s) 2022 NATURE COMMUNICATIONS | (2022)13:1553 | https://doi.org/10.1038/s41467-022-29283-8 | www.nature.com/naturecommunications 11