ARTICLE https://doi.org/10.1038/s41467-019-08853-3 OPEN Global monitoring of antimicrobial resistance based on metagenomics analyses of urban sewage Rene S. Hendriksen1, Patrick Munk 1, Patrick Njage1, Bram van Bunnik 2, Luke McNally3, Oksana Lukjancenko1, Timo Röder1, David Nieuwenhuijse4, Susanne Karlsmose Pedersen1, Jette Kjeldgaard1, Rolf S. Kaas1, Philip Thomas Lanken Conradsen Clausen1, Josef Korbinian Vogt1, Pimlapas Leekitcharoenphon1, Milou G.M. van de Schans5, Tina Zuidema5, Ana Maria de Roda Husman6, Simon Rasmussen 7, Bent Petersen7, The Global Sewage Surveillance project consortium#, Clara Amid8, Guy Cochrane8, 1234567890():,; Thomas Sicheritz-Ponten9, Heike Schmitt6, Jorge Raul Matheu Alvarez10, Awa Aidara-Kane10, Sünje J. Pamp1, Ole Lund7, Tine Hald1, Mark Woolhouse2, Marion P. Koopmans4, Håkan Vigre1, Thomas Nordahl Petersen1 & Frank M. Aarestrup 1 Antimicrobial resistance (AMR) is a serious threat to global public health, but obtaining representative data on AMR for healthy human populations is difficult. Here, we use meta- genomic analysis of untreated sewage to characterize the bacterial resistome from 79 sites in 60 countries. We find systematic differences in abundance and diversity of AMR genes between Europe/North-America/Oceania and Africa/Asia/South-America. Antimicrobial use data and bacterial taxonomy only explains a minor part of the AMR variation that we observe. We find no evidence for cross-selection between antimicrobial classes, or for effect of air travel between sites. However, AMR gene abundance strongly correlates with socio- economic, health and environmental factors, which we use to predict AMR gene abundances in all countries in the world. Our findings suggest that global AMR gene diversity and abundance vary by region, and that improving sanitation and health could potentially limit the global burden of AMR. We propose metagenomic analysis of sewage as an ethically acceptable and economically feasible approach for continuous global surveillance and pre- diction of AMR. 1 NationalFood Institute, Technical University of Denmark, Kgs. Lyngby 2800, Denmark. 2 Usher Institute, University of Edinburgh, Edinburgh EH8 9AG, UK. 3 Centre for Synthetic and Systems Biology, School of Biological Sciences, University of Edinburgh, Edinburgh EH9 3JD, UK. 4 Viroscience, Erasmus Medical Center, Rotterdam 3015, The Netherlands. 5 RIKILT Wageningen University and Research, Wageningen 6708, The Netherlands. 6 National Institute for Public Health and the Environment (RIVM), Bilthoven 3721, The Netherlands. 7 Department of Bio and Health Informatics, Technical University of Denmark, Kgs. Lyngby 2800, Denmark. 8 European Molecular Biology Laboratory, European Bioinformatics Institute, Hinxton CB10 1SD, UK. 9 Centre of Excellence for Omics-Driven Computational Biodiscovery, AIMST University, Kedah 08100, Malaysia. 10 World Health Organization, Geneva 1202, Switzerland. #A full list of consortium members appears at the end of the paper. Correspondence and requests for materials should be addressed to F.M.A. (email: fmaa@food.dtu. dk) NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 A ntimicrobial resistance (AMR) is a cross-cutting and Results increasing threat to global health1,2, and it threatens to Global distribution of AMR genes. Domestic sewage was col- undermine decades of progress in the treatment of lected from 79 sample locations, covering 7 geographical regions infectious diseases. AMR is a complex problem with multiple and from 74 cities in 60 countries (Fig. 1a, Supplementary Data 1). interconnected drivers, which may include changing dynamics in Each sample was sequenced using Illumina HiSeq and the travel, trade, climate change, and populations. Reliable informa- resulting data (>1.4 Tb) processed using MGmapper17. The tion that accurately describes and characterizes the global average number of reads per sample was 120 million reads (range: occurrence and transmission of AMR is essential to address this 8 million –398 million). An average of 0.03% of the reads were challenge and to support national and global priority setting, assigned to AMR genes, while on average 29%, 1%, 0.4%, and public health actions, and treatment decisions. 0.2% were assigned to bacterial, protozoa, plants, and human Current surveillance of AMR is often focusing on a few genomic material, respectively (Supplementary Data 2). Sixty- pathogens only and mainly based on passive reporting of phe- eight percent of the reads could not be assigned to any reference notypic laboratory results for specific pathogens isolated from sequence, and other metagenomic studies have also found larger human clinical infections1,3–5. This procedure leads to significant number of un-mapped reads (42%–48%)14,18. Rarefaction of the time delays, often incomparable data, and a narrow pathogen reads mapping to bacterial genomes showed a tendency toward spectrum not capturing all relevant AMR genes, where the major saturation in the sequence data (Supplementary Fig. 1). part might be present in the commensal bacterial flora of healthy Analyses of duplicate samples from separate days from eight individuals. However, obtaining fecal samples from healthy sites showed a high degree of within-site reproducibility humans is logistically difficult. (Supplementary Fig. 2). Comparison of the data from the same From a surveillance point of view, urban sewage is attractive countries showed much less (permutation test, p < 0.0001) because it provides sampling material from a large and mostly variance across sites within countries than across sites between healthy population, which otherwise would not be feasible to different countries (Supplementary Fig. 3), which suggests that a monitor. Globally, a rapidly increased proportion of the human single sample taken from one large city is representative of the population live in urban areas6 and an increasing proportion overall occurrence of AMR in a country. is connected to a sewer system7,8. In addition, analyzing sewage A total of 1546 genera were detected across all samples (range: samples does not require informed consent, thus limiting 942–1367 genera per sample), but a limited number of bacterial ethical concerns and has limited practical and logistical genera dominated (Supplementary Data 3, Supplementary Fig. 4). barriers for sampling. Most microbiological studies on Several of the dominant bacterial genera were typical fecal, such sewage have focused on the risk of discharge of insufficiently as Faecalibacterium, Bacteroides, Escherichia, Streptococcus, and treated sewage or problems related to heavy rainfall overflow, Bifidobacterium. However, other highly abundant bacterial but some evaluations on the surveillance of pathogens have genera, such as Acidovorax and Acinetobacter, are most likely also been performed9,10. Additionally, sewage has proven environmental bacteria. Thus, the bacterial composition of useful for surveillance in the global polio eradication sewage is complex and does not only reflect human feces but program11,12. also the changes occurring in the sewer. A comparison with Metagenomic techniques, using short-read next-generation publicly available metagenomic data, although generated using sequencing data, benefit from the ability to quantify thousands different DNA-purification methods, suggested that our urban of especially transmissible resistance genes in a single sample. sewage samples resemble more the human fecal microbiome than Moreover, it can provide additional information about the the animal fecal microbiome from chickens, pigs, or mice presence of bacterial species, pathogens, and virulence genes (Supplementary Fig. 5). and the data can be reanalyzed if novel genes of interest are The total AMR gene abundances varied across sites and identified. It should, however, also be acknowledged that short- continents (Fig. 1b, Supplementary Data 4). The highest AMR read metagenomics might provide limited information gene levels were observed in African countries (average: 2034.3 regarding the host of the genes or the genetic environment. fragments per kilo base per million fragments (FPKM)), although Metagenomics has been found to be superior to conventional Brazil had the highest abundance of all (4616.9 FPKM). At the methods for AMR surveillance in pig herds13 and has also been lower end of the spectrum were Oceania (New Zealand and utilized for the surveillance of global AMR gene dissemination Australia) (average: 529.5 FPKM). To the best of our knowledge, through international flights14. Additionally, an extensively comparable data on the global occurrence of AMR genes of shared resistome was observed across urban sewage samples predominantly healthy people do not exist. In agreement with our within China15, as well as between individuals and environ- findings, a previous study on the resistome from toilet waste from mental samples in Lima, Peru16. Interestingly, Pehrsson et al.16 long-distance flights14 suggested that the AMR levels in South showed that even though changes in the bacterial composition Asia were higher than in Europe. Data on AMR in bacteria were observed between feces and sewage, this was not the case isolated from clinical infections in humans, collected by the for AMR genes. World Health Organization (WHO)1, suggest a high prevalence Here we use metagenomic analysis of untreated sewage to of AMR in many developing countries, even though several of the characterize the bacterial resistome from 79 sites in 60 coun- national studies reported by the WHO give contradictory results1. tries. We find systematic differences in abundance and diversity A total of 1625 different AMR genes belonging to 408 gene of AMR genes between Europe/North-America/Oceania and groups were identified, including several that have emerged Africa/Asia/South-America. Antimicrobial use data and bac- recently, such as CTX-M, NDM, mcr, and optrA (Supplementary terial taxonomy only explain a minor part of the AMR variation Data 4). Several different AMR genes might encode resistance to the that we observe. However, AMR gene abundance strongly same antimicrobial agent. Thus the relative abundance of AMR correlates with socio-economic, health, and environmental genes was aggregated to the corresponding antimicrobial class level factors, which we use to predict AMR gene abundances in all for each sample to explore major trends across countries (Fig. 1c). countries in the world. Our findings suggest that global AMR AMR genes encoding resistance toward macrolides, tetracyclines, gene diversity and abundance vary by region and that aminoglycosides, beta-lactams, and sulfonamides were the most improving sanitation and health could potentially limit the abundant. Most samples from Europe and North America had a global burden of AMR. high relative proportion of macrolide resistance genes, while Asian 2 NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 ARTICLE a b 4000 Total AMR per sample [FPKM] 3000 2000 Africa 1000 Asia Europe Middle East 0 North America a ca e a a ia st i ic op ic As an Ea ri er er Oceania Af r ce Eu Am Am e O dl South America id th th M or u So N c Africa Asia Europe ME North America Oc. S. Amer. AMR class 1.00 AmGlyc/Quin Aminoglycoside Beta−Lactam 0.75 Colistin Relative abundance Fosfomycin Mac/Oxa/Phen Macrolide 0.50 Nitroimidazole Oxa/Phen Phenicol Quinolone 0.25 Rifampicin Sulphonamide Tetracycline 0.00 Trimethoprim Vancomycin BWA.19 CIV.13 ETH.24 GHA.4 KEN.72 NGA.50 SEN.8 TGO.44 TZA.15 ZAF.39 ZMB.49 ZMB.49b CHN.64 IND.11 KAZ.6 KHM.21 LKA.40 MYS.54 NPL.33 PAK.7 SGP.52 VNM.48 ALB.17 AUT.70 BGR.66 CHE.67 CZE.23 DEU.27 DNK.71_RA DNK.71_RD DNK.71_RL ESP.75 FIN.25 GEO.59 HRV.68 HUN.61 IRL.69 ISL.28 ITA.30 LUX.32 LVA.31 MDA.65 MKD.62 MLT.63 NLD.43 NOR.34 POL.36 SRB.37 SVK.9 SVN.38 SWE.41 SWE.41a XK.60 IRN.12 ISR.29 TUR.46 CAN.22 CAN.22a CAN.22b CAN.22c USA.74 USA.74a USA.74b USA.74c USA.74d USA.74e USA.74f USA.74g USA.74h USA.74i AUS.18 AUS.18a NZL.56 BRA.53 BRA.53a COL.2 ECU.14 ECU.14a PER.35 d Africa Asia Europe ME North America Oc. S. Amer. AMR gene 1.00 aadA_clust1 blaOXA_clust21 blaOXA_clust8 0.75 erm(B)_clust Relative abundance erm(F)_clust mef(A)_10 mph(E) 0.50 msr(D) msr(E) strA strB 0.25 sul1_sul3_clust tet(39) tet(Q) 0.00 tet(W) Other BWA.19 CIV.13 ETH.24 GHA.4 KEN.72 NGA.50 SEN.8 TGO.44 TZA.15 ZAF.39 ZMB.49 ZMB.49b CHN.64 IND.11 KAZ.6 KHM.21 LKA.40 MYS.54 NPL.33 PAK.7 SGP.52 VNM.48 ALB.17 AUT.70 BGR.66 CHE.67 CZE.23 DEU.27 DNK.71_RA DNK.71_RD DNK.71_RL ESP.75 FIN.25 GEO.59 HRV.68 HUN.61 IRL.69 ISL.28 ITA.30 LUX.32 LVA.31 MDA.65 MKD.62 MLT.63 NLD.43 NOR.34 POL.36 SRB.37 SVK.9 SVN.38 SWE.41 SWE.41a XK.60 IRN.12 ISR.29 TUR.46 CAN.22 CAN.22a CAN.22b CAN.22c USA.74 USA.74a USA.74b USA.74c USA.74d USA.74e USA.74f USA.74g USA.74h USA.74i AUS.18 AUS.18a NZL.56 BRA.53 BRA.53a COL.2 ECU.14 ECU.14a PER.35 Sample Fig. 1 Global sewage sampling sites and overview of antimicrobial resistance (AMR) abundance and composition. a Map of the sampling sites. b Boxplots of the total AMR fragments per kilo base per million fragments per sample, stratified by region. Each sample is represented by a dot with horizontal jitter for visibility. The horizontal box lines represent the first quartile, the median, and the third quartile. Whiskers denote the range of points within the first quartile − 1.5× the interquartile range and the third quartile + 1.5× the interquartile range. c Relative AMR abundance per antimicrobial class (AmGlyc aminoglycoside, Mac macrolide, Oxa oxazolidinone, Phen phenicol, Quin quinolone). d Relative abundance of the 15 most common AMR genes (mef(A) _10: mef(A)_10_AF376746) and African samples had a large proportion of genes providing Global diversity and clustering of AMR genes. We analyzed the resistance to sulfonamides and phenicols. Fifteen AMR genes AMR abundances on both the gene and antimicrobial class contributed >50% of the total AMR abundance (Fig. 1d). This levels using both principal coordinate analyses (PCoAs) and proportion was especially prominent for Europe, North-America, heat maps (Fig. 2, Supplementary Fig. 6). With regard to the and Oceania. None of the dominant AMR genes are known to be sample resistome dissimilarities, there was a clear geographical restricted to specific bacterial genera19–21. separation along the first principal coordinate of samples from NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 a 0.2 PCo 2: 15.66 % variation Africa 0.1 Asia Europe Middle East 0.0 North America Oceania South America −0.1 −0.2 −0.1 0.0 0.1 0.2 PCo 1: 21.79 % variation b Region 7 Macrolide Vancomycin 6 Nitroimidazole Rifampicin 5 AmGlyc/Quin Trimethoprim 4 Phenicol Sulphonamide 3 Tetracycline Aminoglycoside 2 Beta−Lactam Mac/Oxa/Phen 1 Oxa/Phen Colistin Fosfomycin 0 Quinolone BWA.19 ISR.29 SEN.8 VNM.48 SGP.52 MYS.54 ZAF.39 NGA.50 ZMB.49 ZMB.49b KEN.72 PAK.7 IND.11 TZA.15 NPL.33 CIV.13 KHM.21 GHA.4 TGO.44 CHN.64 MLT.63 BRA.53a ETH.24 COL.2 LKA.40 BRA.53 IRN.12 ECU.14a PER.35 HUN.61 DNK.71_RL SVN.38 USA.74 SWE.41a USA.74a USA.74d CAN.22b AUS.18a USA.74f USA.74h USA.74e USA.74b USA.74c USA.74g GEO.59 MDA.65 AUS.18 ISL.28 NZL.56 AUT.70 IRL.69 NOR.34 ITA.30 TUR.46 KAZ.6 ESP.75 USA.74i CZE.23 DNK.71_RD SRB.37 MKD.62 ALB.17 BGR.66 LUX.32 XK.60 ECU.14 NLD.43 DNK.71_RA CHE.67 HRV.68 CAN.22 CAN.22a CAN.22c POL.36 SVK.9 FIN.25 SWE.41 DEU.27 LVA.31 Region: Africa Asia Europe Middle East North America Oceania South America Fig. 2 Resistome clustering in sewage samples across regions. a Principal coordinate analysis (PCoA) performed on the resistome Bray–Curtis dissimilarity matrix. The amount of variation explained by coordinates 1 and 2 is included in the axis labels. b Antimicrobial resistance class-level heat map. Relative abundances of genes (fragments per kilo base per million fragments (FPKM)) were summed to drug classes (AmGlyc aminoglycoside, Mac macrolide, Oxa oxazolidinone, Phen phenicol, Quin quinolone). Colors represent log (ln) transformed relative abundances (FPKM). Complete-linkage clustering of Pearson correlation coefficients was used to hierarchically cluster both samples and drug classes Europe/North-America/Oceania and samples from Africa/Asia/ Analyses of the alpha diversity of the bacterial taxonomic South-America. Regional groupings explained 27% of the dis- composition showed less separation according to geographical similarity between sample resistomes (adonis2, p < 0.001; regions compared with the AMR genes (Supplementary Figs. 9– Fig. 2a). The separation among the groups was mainly driven 11), and regions explained less of the variation in the resistome by higher levels of resistance to tetracycline, aminoglycosides, (adonis2, 27%) compared with the bacterial composition beta-lactams, sulfonamides, and trimethoprim in the Africa/ (adonis2, 17%). To test the degree to which bacterial genus- Asia/South America cluster, whereas macrolide resistance was level composition of the microbiota is associated with the more evenly distributed among all samples (Fig. 2b). A clearer resistomes, Procrustes analyses were performed (Supplementary clustering was observed based on regions, compared with Fig. 12). We compared ordinations of the bacterial taxonomic clustering based on diet, income, or the Human Development composition with the resistome and found that they correlated Index (HDI) (Supplementary Figs. 6–7). A stronger regional significantly (protest, p < 0.001). separation was observed on the AMR class level, compared with the gene level (Fig. 2b, Supplementary Fig. 6). On the class level, a very clear separation of all samples in two groups was AMR genes and drug use association. Several studies have observed, where only a single European sample (Malta) did not shown that antimicrobial use (AMU) selects for AMR22,23 and cluster with all other samples from Europe/North America/ that reducing AMU often results in reduced AMR24,25, except Oceania (Fig. 2b). Furthermore, only one sample from Asia when there is interference between genetically linked AMR genes (Kazakhstan), one from the Middle East (Turkey), and one conferring resistance to different antimicrobial classes26. It has from South America (Galapagos Islands) clustered with the also been suggested that AMU explains only some of the Europe/North-America/Oceania group. variation22,23 and that other factors such as diet, cultural tradi- Several alpha diversity indices for each sample resistome were tions and occupation also have an influence22,27. calculated (Supplementary Fig. 8). The resistome median In this study, the association between AMU and the occurrence evenness was higher in African and Asian countries compared of AMR in the sewage samples was estimated using a generalized with resistomes from other geographical locations. Thus not only linear mixed-effects model, with the counts of genes in the do these regions have a higher prevalence of AMR genes, but also different antimicrobial classes as an outcome (Poisson) adjusting a more equal distribution of the different AMR genes. for sequencing depth and gene length. As AMU data, we included 4 NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 ARTICLE 2015 data from Europe [www.ecdc.dk] and IQVIA, formerly variables, we were able to explain up to 89% of the observed Quintiles IMS Holdings, Inc. (see Methods). In the regression variation across the samples (Supplementary Fig. 14, Supplemen- model, we accounted for the potential effects of cross-resistance tary Table 2). Most of the variables associated with AMR levels by fitting fixed effects of both usage of the antimicrobial class that were related to sanitation and general health (Fig. 3, Supplemen- a resistance gene primarily confers resistance to (direct selection tary Fig. 15, Supplementary Table 3). Subsequently, the identified for resistance) and the total AMU (indirect selection via cross- variables were used to predict the occurrence of AMR in 259 resistance). While our model showed a significant increase in the countries and territories. The three countries predicted to have abundance of AMR genes belonging to a specific antimicrobial the lowest level of AMR were The Netherlands, New Zealand, and class with increasing usage of that antimicrobial class, we found Sweden, whereas the highest predicted AMR levels were for no significant effect of total usage of all antimicrobials on Tanzania, Vietnam, and Nigeria (Fig. 4, Supplementary Data 5). abundance of the different classes (Supplementary Fig. 13A, The predicted global country-level resistance levels were multi- Supplementary Table 1). This suggests that, while AMU of a plied with the latest national population estimate and used to specific class is an important driver of AMR genes encoding create global maps of healthy human-associated AMR (Supple- resistance to that class, the effects of cross- and/or co-resistance mentary Fig. 16). appear to have a relatively minor contribution to AMR abundances. Furthermore, our model showed that the countries with a lower HDI (i.e., higher in rank) have lower abundances of Discussion AMR genes (glmm, p = 0.01) and that the number of passengers A reliable base of evidence that accurately describes and char- arriving in a country via flights has no effect on the abundance of acterizes the global burden and transmission of AMR is essential AMR genes (glmm, p = 0.62). to support national and global priority setting, public health A second (similar) model was developed to test the association actions, and treatment decisions. Compared to samples directly between the abundance of AMR genes on the antimicrobial class obtained from humans, a major advantage of sewage is that such level and the antimicrobial residue levels. Again we accounted for samples can be easily obtained and analyzed without ethical the potential effects of cross-resistance by fitting fixed effects of concerns. Thus, while current efforts to improve surveillance of both residues of the antimicrobial class that a resistance gene AMR in human clinical pathogens should be continued1,3–5, we primarily confers resistance to (direct selection for resistance) and do suggest that our study provides the foundation for a flexible, the total antimicrobial residue levels (indirect selection via cross- simple, affordable, and ethically acceptable global real-time sur- resistance). As before, this model showed a significant increase in veillance of AMR that could be immediately implemented glob- the abundance of AMR genes belonging to a specific antimicro- ally also in low- and middle-income countries. The study design bial class with increasing levels of drug residues of that can furthermore be used for other infectious disease agents. antimicrobial class; however, we found no significant effect of Our study represents, to the best of our knowledge, the first total residue levels (Supplementary Fig. 13B). As with the attempt to monitor and predict the occurrence of AMR in the previous model, this model also showed that the countries with global predominantly healthy human population. Even though a lower HDI (i.e., higher in rank) were less abundant in terms of the study suffers from the limitation that only a single sample was AMR genes (glmm, p < 0.001), and that the number of passengers analyzed from each site, our study suggests a strong systematic arriving in a country via flights had no effect on the abundance of separation of regions of the world according to AMR gene AMR genes (glmm, p = 0.746), indicating that the results are abundance with high-income countries in Europe/North-Amer- robust. ica/Oceania constituting one cluster and low-income countries in Of interest is that there was no correlation between AMU and Africa/Asia/South-America constituting another cluster. This antimicrobial residue in our data (lm, R2 = 0.00098, p = 0.55). separation was mainly driven by a relatively high abundance of a One possible explanation for this could be that AMR abundances limited number of AMR genes encoding macrolide resistance accumulates on a long timescale because of the average usage in a genes in Europe/North-America/Oceania, compared with a gen- country (i.e., (long term) AMU impacts the level of resistance); eral high abundance of several different AMR genes from dif- however, on a much shorter timescale fluctuations in AMU will ferent classes in the rest of the world. The countries standing out change the AMR composition (diversity) somewhat, which is as having the most divergent distribution of AMR genes were captured by the residue levels, which can be interpreted as a Vietnam, India, and Brazil, suggesting that these countries could (short) temporal snapshot of the actual usage. be hot spots for emergence of novel AMR mechanisms. In this study, we used metagenomics that benefit from the ability to quantify thousands of genes simultaneously and that the Prediction of AMR based on population-level health data. We data can be reanalyzed if novel genes of interest are identified. observed that AMU only explained a minor part of the occur- However, other methodologies might also be useful such as cul- rence of AMR across the world. In addition, AMU data are dif- ture and PCR-based methods that might have better sensitivity. ficult to obtain and likely subject to limitations due to the lack of Comparative studies evaluating the usefulness of various tech- an effective prescription system in many countries. Measuring nologies, including evaluation of sensitivity, specificity, and antimicrobial residue levels in sewage as a proxy for AMU is also number of targets detected, are warranted. associated with uncertainties and the rapid degradation in the In our study, we have focused on total AMR abundance. environment of the heavily used beta-lactamase-sensitive anti- However, the resistance to the different antimicrobial classes are microbials makes it difficult to reliably measure them. not equally important28 neither are all AMR genes29. Thus fur- Because the HDI was strongly associated with AMR, in the ther studies could also benefit from addressing specific resistance results of our model, we hypothesized that a number of other genes. factors could be either drivers or indicators for AMR. To An analysis of the same IQVIA AMU data from 2000 to investigate this further, we used 1503 variables from the World 2015, as we have used here, was recently published30, and even Bank’s Health, Nutrition and Population as well as Development though a significant increase in AMU was observed, especially indicator data sets collected between the years 2000 and 2016 for for countries in Africa, Asia, and South-America, the quantities 259 countries and territories to study the potential association are still below most countries in Europe and North-America. with the observed AMR gene abundances. Using country-specific We could not in our study find significant associations between NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 a Female child mortality rate b Water and sanitation investment c Total death reporting d Informal employment e Time to import 1000 1000 1000 1000 1000 750 750 750 750 750 500 500 500 500 500 250 250 250 250 250 0 0 0 0 0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 f External debt grace period g Educational attainment h Infection and malnutrition i Maternal death risk j Rural improved water source 1000 1000 1000 1000 1000 750 750 750 750 750 500 500 500 500 500 250 250 250 250 250 0 0 0 0 0 FPKM 0.90 0.95 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 k Number of surgical procedures l Female life expectancy m Number of physicians n Open defecation practice o Child diarrhea treatment 1000 1000 1000 1000 1000 750 750 750 750 750 500 500 500 500 500 250 250 250 250 250 0 0 0 0 0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 p Life expectancy q Surgical care expenditure risk r Skilled birth attendants 1000 1000 1000 750 750 750 500 500 500 250 250 250 0 0 0 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 Value Fig. 3 World Bank variables significantly associated with the observed antimicrobial resistance abundances. Detailed information concerning the variables in a–r are presented in the same order in Supplementary Table 3 Fig. 4 Global predictions of antimicrobial resistance (AMR) abundance in all countries and territories in the world. Map colored according to predicted abundance of AMR from light blue (low AMR abundance) to dark blue (high AMR abundance). Global resistance predictions for the 259 countries and territories are shown in Supplementary Data 5 obtainable AMU data and the concentrations of residues Historically, most strategies to reduce AMR have focused on measured, but we did observe the highest concentrations of reducing AMU, which relies upon AMR imparting a fitness cost antimicrobial agents in the African samples, which could sug- on the bacterial host, an effect that has been relatively weak in gest that measuring residues in sewage could provide alternative horizontally transferred AMR genes compared with chromoso- data for monitoring AMU compared to obtainable sales data. mal mutations31. Other factors, such as those related to 6 NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 ARTICLE transmission, including infection control, sanitation, access to For each batch of DNA extractions, a DNA extraction blank control was processed clean water, access to assured quality antimicrobials and diag- in parallel with the sewage samples to monitor background DNA. nostics, travel, and migration, have also been suggested to con- tribute significantly to AMR27. In this study, we found that, Whole-community sequencing. DNA was shipped on dry ice for library pre- irrespective of the diversity of AMR genes, the total AMR paration and sequencing to the Oklahoma Medical Research Foundation (OMRF). DNA from all samples was mechanically sheared to a targeted fragment size of 300 abundance was highly correlated with a limited number of World bp using ultrasonication (Covaris E220evolution). Library preparation was per- Bank variables, mainly concerning sanitation and health. formed with the NEXTflex PCR-free Library Preparation Kit (Bioo Scientific). The In contrast, human air travel had no significant influence on Bioo NEXTflex-96 adapter set (Bioo Scientific) was used, and in batches of roughly AMR abundance. Importantly, this suggests that the total AMR 60 samples, the libraries were multiplexed and sequenced on the HiSeq3000 platform (Illumina), using 2 × 150-bp paired-end sequencing per flow cell. The raw abundance is mainly influenced by local/national parameters, and sequencing data have been deposited at the European Nucleotide Archive under even though all AMR genes might rapidly disseminate and be accession number ERP015409. found in all corners of the world, local selection is required for them to reach appreciable frequencies. These findings Trimming and mapping of sequencing reads. The reads were trimmed, including suggest that improving sanitation, health, and perhaps education adaptor removal, using BBduk [BBMap—Bushnell B.—https://sourceforge.net/ as part of the Sustainable Development Goals [www.un.org/ projects/bbmap/] with a quality threshold at 20 and minimum length of 50 bp. sustainabledevelopment/sustainable-development-goals/] would Trimmed reads were used as input to the reference-based mapping and taxonomy- assignment tool MGmapper17, which is based on BWA-MEM33 version 0.7.12 and be effective strategies for limiting the global burden of AMR. SAMtools34 version 1.6. Reads were aligned against reference sequence databases for the best hit (Bestmode, i.e., a read can only map to 1 reference sequence). An acquired AMR gene database (ResFinder)35 was used to annotate properly paired Methods reads (MGmapper option fullmode) where each read pair had sequence coverage of Ethics. It was confirmed that this and similar studies using human sewage in at least 80% compared with the length of the trimmed reads. This was the only accordance with the Danish Act on scientific ethical treatment of health research filter that was applied to discard a read pair. The AMR genes were of bacterial (Journal no.: H-14013582) do not require preapproval from ethical committees. origin and could therefore align to both bacteria databases and the Resfinder database. To enable multiple database hits, AMR genes were mapped using the Fullmode option in MGmapper for the most optimal abundance calculation. Collection of urban sewage samples. The National Food Institute, Technical Genomic annotation was performed by identifying the best hit (MGmapper option University of Denmark (DTU Food) launched an open invitation 25 September bestmode) for a pair of reads when aligned against a range of reference sequence 2015 seeking potential collaborators to participate in the pilot study of the Global databases. Databases were primarily downloaded via NCBI genbank clade specific Sewage Surveillance Project (GSSP). The invitation was sent by electronic mail to assembly_summary.txt files unless another ftp site is provided below. The list of the following networks and individual research collaborators for further dis- databases used by MGmapper includes: Human (GRCh38.p3), bacteria (closed semination: WHO Global Foodborne Infections Network (GFN) [http://www.who. genomes), MetaHitAssembly (PRJEB674—PRJEB1046), HumanMicrobiome int/gfn/en/], WHO Advisory Group on Integrated Surveillance of Antimicrobial (genomes assemblies), bacteria_draft, plasmid, archaea, virus, fungi, protozoa, Resistance (AGISAR) [http://www.who.int/foodsafety/areas_work/antimicrobial- vertebrates_mammals, vertebrates_other, invertebrates, plant, and nt. For the resistance/agisar/en/], European Union Reference Laboratory for AMR network bacteria and bacteria_draft databases, sequences were selected from the assem- [www.eurl-ar.eu/], and to the European Food- and Waterborne Diseases and bly_summary.txt file, when annotated with the tags version_status=‘latest’ and Zoonoses Network (FWD-Net) [https://ecdc.europa.eu/en/about-us/partnerships- genome_rep=‘Full’. Furthermore, assembly_level= ‘Complete genome’ or ‘Chro- and-networks/disease-and-laboratory-networks/fwd-net]. mosome’ were required for entries in the bacteria database and refseq_categor- Participation was cost neutral and managed by online registration that also y=‘representative genome’ for entries in the bacteria_draft database. The plasmid included requests for proper packaging and sample containers and the need to sign database was constructed as the subset of bacteria and bacteria_draft sequences a material transfer agreement or similar approval/agreement to protect any having the word ‘plasmid’ in the fasta entry header line. intellectual property rights. All expenses related to the shipments were paid by The total bacteria read count for a sample was calculated as the sum of read DTU Food, and the shipments complied with the IATA regulation as per SP A197 counts from each of the bacteria-related databases (bacteria, bacteria_draft, because the content of the parcels was identified as UN3082 “Environmentally MetaHitAssembly, and HumanMicrobiome). The total fraction of unmapped reads hazardous substance, liquid, not otherwise specified” but did not exceed 5 L of for all sample sites were used to translate the percentage of unmapped reads into Z- sewage. scores, i.e., the number of standard deviations from the mean. Setting an absolute A protocol that instructed how to collect the urban sewage samples as well as Z-score threshold at 3 retained data from 79 sample sites. Data from Chad were epidemiological, demographic, and geographical information was provided to each excluded based on the Z-score threshold together with data from Gambia with participant in the study (Supplementary Data 8). Participants were requested in suspiciously low resistance read counts; i.e., lower than those observed from DNA conjunction with collecting the samples to submit via an online survey the captured extraction control samples. information related to epidemiological, demographic, and geographical information as well as a digital image of the sampling site, if possible the GPS Bacterial and AMR gene distribution. Inspection of the count tables and the coordinates of the sampling location, the temperature of the sample at the time of mapped reads distribution on the genomes revealed an overestimation of some sampling, pH of the sewage, and storage temperature of sample. From each genomes due to plasmids with high copy numbers. The issue only occurred in the location, on 2 consecutive days between 25 January and 5 February 2016, one included draft genomes, as the plasmid DNA were left out of the database with representative, non-processed, unfiltered urban sewage sample of 2 L was collected complete bacterial genomes. The distribution of mapped reads across contigs was from the respective main sewage pipeline(s) prior to the inlet of the wastewater expected to be evenly distributed, with some variation. The plasmids were revealed treatment plant or from the main outlet to rivers or similar according to the by large hit counts to specific contigs compared with the associated contigs in a protocol provided. Flow proportion sampling over 24 h was preferred. draft genome. For each draft genome, the hits to each contig was normalized with Alternatively, three crude point samples were collected in a short time interval, i.e., respect to contig size. The median of the contig hits was found and the third at least 5 min between each individual sample, to ensure as much randomness as quartile and the interquartile range was calculated. If a hit count was above the possible. The collected urban sewage was not treated, neither with additives nor third quartile plus 1.5 times the interquartile range, then the hit count was inter- DNA stabilizers and was recommended to be stored at −80 °C for at least 48 h preted as an overestimation and adjusted by replacing the hit count with the prior to shipment to avoid the use of dry ice. The sample-specific data are provided median. in Supplementary Data 1 and an explanation of the content of this file provided as Relative abundance of AMR genes and bacterial genera were calculated as Supplementary Note 1. FPKM. This was done to account for both sample-wise sequencing depth differences and size-dependent probability of observing a reference. For bacterial genome assemblies, FPKM was calculated based on the adjusted sum of fragments Sample handling and DNA extraction. At DTU Food, a photograph of each assigned to a genome assembly, whether or not the genome was closed. For AMR sample container was taken upon arrival combined with a short remark describing genes, FPKM was calculated on each individual ResFinder reference sequence. the condition of the sample, e.g., solid frozen, thawed, coloration, etc. The samples FPKMs were subsequently summed up across categories to bacterial genera (NCBI of the first of the consecutive days were thawed for 12 h at approximately 20 °C taxid), drug class level (NCBI tax ID), and highly homologous AMR gene groups before processing. After thawing, 250 mL of each sample were pelleted in a cen- (CD-HIT-EST, 90% similarity)36. trifuge at 10,000 × g for 10 min. The pellet was stored at −20 °C or −80 °C before DNA extraction and metagenomics analysis, and the supernatant was stored at −80 °C for subsequent antimicrobial residue and virus analysis. DNA was extracted Within-site reproducibility. In order to test reproducibility of sewage samples, a from the sewage pellets according to an optimized protocol using the QIAamp Fast dendrogram of resistome composition from all sewage samples including samples DNA Stool Mini Kit including twice the input material and initial bead beating32. from eight sites that were double-sampled sewage was generated using Bray–Curtis NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 (BC) dissimilarities and hierarchical clustering (Supplementary Fig. 2). The replicated Procrustes analysis. The vegan package was used for comparing the resistome samples were taken by 2 days apart. The day 2 samples were taken twice. The day dissimilarities with the bacteriome dissimilarities. The protest function was used to 2 samples were kept in freezer for 2 years prior sequencing. Samples from day 1 and scale and rotate the principal coordinates of the bacterial PCoA onto the principal day 2 were sequenced using the same DNA extraction and sequencing protocols. The coordinates of the resistome PCoA and testing the strength of the association. replicated samples had higher resistome similarity to its own replicates and all eight sets of replicate samples were clustered with their own replicate. This result shows a strong reproducibility of sewage samples despite having different day of sampling and Graphics and statistics. All plots and statistical analyses were carried out in 2 years of storage. The repeats were not included in subsequent analyses. Microsoft R Open 3.3.2. Within-country representativeness. To assess whether samples from individual Correlation between AMU and external explanatory variables. A multilevel sites are representative of other sites in that country, we compared the BC dis- Poisson model was developed to investigate the sources of variance for the similarities for pairs of sites within the same country and in different countries for abundance of AMR genes and the relationship between AMU and abundance of both resistome and bacteriome compositions. We assessed the significance of these AMR genes found in the samples. The counts of the individual AMR genes in each differences using permutation tests. We permuted the country labels for each of the samples (see “Collection of urban sewage samples” and “Whole-community sample and reassessed the dissimilarities for pairs of sites within the same country sequencing“ sections) aggregated at the antimicrobial class level was used as the and in different countries with permuted labels. We repeated this procedure 106 dependent variable. times to build up a null distribution of the differences in dissimilarity within and An observation-level random effect was used to model the over dispersion among countries. We found that resistome dissimilarities were on average 34% inherent to count data41. Because several samples were sequenced more than once higher for pairs of sites in different countries than for pairs of sites within the same on separate sequencing runs, we were able to estimate and correct for the noise in country (permutation, p < 0.0001, Supplementary Fig. 3A and B), while bacteriome mean AMR gene abundance owing to the sampling process. Therefore, a dissimilarities were 46% higher for pairs of sites in different countries than for pairs categorical variable indicating which sample was used was included as a random of sites within the same country (p < 0.0001, Supplementary Fig. 3C and D). These effect (sample). A categorical variable identifying each location to allow for the results show that there is less variance across sites within countries than across sites estimation of variation in abundance of resistance genes between the different within different countries. Thus individual sites in this study are representative of sampling locations (location) was included. Furthermore, a variable identifying the other sites in that country. resistance class a resistance gene is a member of was included to estimate the variance due to differences in abundance between antimicrobial classes (class). For the fixed effects, we included a measure of AMU at the country level. Data from Sample composition comparison. Metagenomes with the following accession the ECDC database (Supplementary Data 1) and data from the IMS database numbers were downloaded from SRA (December 2017): ERR011089, ERR011114, (Supplementary Data 1) were used to calculate a new variable from the two data sets. To ERR011344, ERR1104480, ERR1104481, ERR1121453, ERR1121455, ERR1121556, this end, we used the data from the ECDC data set where available to predict missing ERR1135427, ERR1135431, ERR1135437, ERR1135693, ERR1278103, values from the IMS data set by using a linear regression model. Because the data from ERR1278104, ERR1278105, ERR1414230, ERR1414260, ERR1527239, the two data sets are highly correlated (r = 0.97, p < 0.01, Pearson’s product moment ERR1527247, ERR1558700, ERR1559789, ERR1560016, ERR1560024, correlation), we could infer the (approximate) corresponding ECDC value from the ERR1560100, ERR1655116, ERR1682090, ERR1682101, ERR1698980, ERR186217, IMS value for that particular country. The antimicrobial usage data were then log- ERR1950597, ERR1950599, ERR1950601, ERR469632, ERR469644, ERR469650, transformed. Because it has been argued that many antimicrobial classes show cross- ERR675524, ERR675555, ERR675557, ERR675560, SRR1182511, SRR1202091, resistance, where resistance to one antimicrobial class also confers resistance to another SRR1267595, SRR2175658, SRR2175725, SRR2175750, SRR2751194, SRR2891615, class, we accounted for the potential effects of cross-resistance by fitting effects of both SRR2891618, SRR605600, SRR605634, SRR873603, SRR873608, and SRR924749. usage of the antimicrobial class that a resistance gene primarily confers resistance to The quality of the data was assessed and if necessary adapters were removed and (direct selection for resistance) and total AMU (indirect selection via cross-resistance). trimmed with a quality threshold at 20 and minimum length of 50 bp. All pairwise We also included the total number of passengers arriving in a country in 2015 as a fixed Jaccard distances of the aforementioned metagenomes, all global sewage samples, effect in the model. Data on the number of passengers were extracted from the ICAO and the control samples were calculated with Mash37. The heat map was drawn in international flight database [https://portal.icao.int, downloaded April 2016] and log- R 3.4.4 (pheatmap). transformed. Lastly, we included the United Nations HDI as a fixed effect. HDI data from 2015 were extracted for the United Nations Development Programme website39 Sample dissimilarities. Matrices with relative abundances (FPKM) of AMR genes and log-transformed and scaled before including them in the final model. The same and bacterial genera were Hellinger-transformed using the decostand function in model set-up was used to investigate the association between drug residue levels and vegan. The BC dissimilarity matrices were then calculated using the vegdist func- AMR gene abundance. In this model, data on the drug residues in the samples was tion, also in vegan (Supplementary Data 6). included as a fixed effect instead of the AMU data. All other effects (random and fixed) were kept the same. Others models including temperature at the collection site at the day of sampling and Gross Domestic Product showed no significant association with Heat maps. The relative abundances (FPKM) of AMR genes were log-transformed AMR genes abundance (glmm’s, all p > 0.05). and visualized in a heat map (pheatmap), showing the 50 most abundant genes. The gene-wise dendrogram is based on Pearson product moment correlation coefficients (PPMC), while the sample dendrogram is based on the aforementioned Resistance prediction. World Bank Health, Nutrition and Population as well as BC dissimilarity matrices, not just for the visualized genes, but all genes. Both Development indicator data sets collected between the years 2000 and 2016 for 259 dendrograms are the result of complete linkage clustering. The bacterial genera countries and territories were downloaded from http://databank.worldbank.org/ heat map was produced in the same way, while the AMR class-level heat map data/home.aspx in October 2017 and used to formulate AMR predictive models. differs by using PPMC for hierarchical clustering of both AMR classes and samples. Imputations of missing data were conducted using the missForest R package, which Continent association was included as metadata for all heat maps. For the AMR is a random forest-based technique that is highly computationally efficient for genes, we included World Bank income levels [http://apps.who.int/gho/data/node. high-dimensional data consisting of both categorical and continuous predictors42. metadata.COUNTRY?lang=en], GEMS cluster (identifies countries with similar The final data set consisted of 1503 economic and health indicator variables dietary intake)38, the HDI39, and GBD 2015 super-regions40. (Supplementary Data 7). The most important variables predicting total resistance (FPKM) were selected from the World Bank data set and a recursive feature elimination method from the R library caret (Supplementary Table 3). The model Sample ordination. Dissimilarity matrices were subject to classical multi- was also run on all 1503 variables for comparison (Supplementary Table 4) and a dimensional scaling (PCoA) to obtain the first two principal coordinates as well as good correlation was observed. The machine learning algorithms Support Vector the variance explained by each, ignoring negative eigenvectors. This was done using Machine and random forest were compared for accuracy in predictions based on the cmdscale command in R. their R2 and Root Mean Square Error43,44. Random forest was the best choice of model (Supplementary Table 2). Random forest is suitable for data sets with many Testing of sample dissimilarities. The BC dissimilarity matrices used for PCoA features, especially where each of the features contributes little information45. The were also used for permutational multivariate analysis of variance (adonis2 func- prediction model was trained for the 60 countries where resistance data were tion in vegan). The geographical group assigned to each sample was used as a available from the current project followed by global predictions of resistance for all predictor for dissimilarity. 259 countries and territories (Supplementary Data 5). Overfitting remains a major hurdle when applying predictive models especially involving many predictors. Breiman45 proved that random forest is protected from Alpha diversity. Diversity and richness were calculated on rarified versions of the overfitting by the incorporation of out-of-bag (OOB) estimates and from the law of resistance and bacterial count matrices. For resistance genes, the two samples with large numbers. This followed from earlier proposals on the use of OOB estimates as a <1000 read pairs were excluded. Subsequently, count matrices were subsampled to key part of estimation of generalization error46,47. In the study of error estimates for the lowest samples’ depth, using the Vegan rarefy function. The Simpson diversity bagged classifiers, Breiman48,49 provided empirical evidence demonstrating same index (1-D), Pielou’s evenness, and the Chao1 richness estimates were calculated accuracy from using the OOB estimate as using a test set of the same size as the using the diversity function in the vegan package. training set and further indicated that the use of the OOB error estimate removes the 8 NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 ARTICLE need for a set aside test set. OOB is the mean prediction error on each training sample therefore requires information on elimination rates and a precise estimate of the xi, using only the trees that did not have xi in their bootstrap sample50. During number of persons connected to a sewer and the total water flow at this specific random forest training, approximately one third of the instances are left out from location (or the average water consumption per person as surrogate). Here each bootstrap training set. This means that the OOB estimates are computed from information on elimination rates and water consumption were not included. In only about one third as many classifiers as in the ongoing main combination. addition, we can only calculate DDD for human drugs used for systematic infec- However, the error rate decreases as the number of combinations increases, which tions. As seen from Supplementary Data 1, the amount of dihysdrostreptomycin is means the OOB estimates will tend to overestimate the current error rate. To get almost absent, the use of Dapson very low, and the contribution of enrofloxacin unbiased OOB estimates, random forests are run past the point where the test set and tylosin 0.1% and 0.3% of the residues of quinolones and macrolides, respec- error converges. This has an added advantage that, unlike cross-validation, where bias tively. For the sulfonamides, the excluded drugs constitute 11%. is present but its extent unknown, the OOB estimates are unbiased51. The concentrations of antimicrobial residues were transformed into DDDs Random forests also aggregate many decision trees to limit overfitting as well as of antimicrobials for systemic use using the WHO DDD database [www.whocc. error due to bias because of the law of large numbers. Random forests limit no/atc_ddd_index]. Some residues were excluded as their predominant usage overfitting without substantially increasing error due to bias due to their ability to was not systemic. This was the case for Dapson (main treatment indication: mitigate the problems of high variance and high bias45. Lepra treatment), sulfacetamide (main indication: acne), sulfadoxine (main However, a ten-fold cross-validation was included during our model building indication: malaria), sulfathiazole (main indication: topical application, with by randomly partitioning model input samples into ten sets of roughly equal size varying dosages), enrofloxacin (animal use—no WHO/ATCC DDD available), followed by estimation of accuracy based on held-out samples. This held-out tylosin (animal use—no WHO/ATCC DDD available), and sample was each time returned to the training set and the procedure was repeated dihydrostreptomycin (animal use—no WHO/ATCC DDD available). The DDD with the second subset held out and so forth. Cross-validation and checking for were summed over the respective antimicrobial classes. Calculation of the sum valid accuracy were also performed with an implication that accuracy scores reduce of DDD across classes was not deemed meaningful, because excretion rates if there is overfitting and only “valid accuracy” is finally used. differ largely between antimicrobial classes and also within antimicrobial classes. Creation of global figures and maps. QGIS 2.18.11 using the cartogram package was used to create colored and distorted maps. Data availability The raw sequencing data have been deposited at the European Nucleotide Archive under Analysis of tetracyclines, sulfonamides, macrolides, and quinolone. The ana- accession number ERP015409. Data on the number of passengers were obtained under lysis of tetracyclines, sulfonamides, macrolides, and quinolones in the sewage license from the ICAO international flight database and is available from the authors samples were performed as in Beredsen et al.52, with adaptations as summarized in with restrictions. Data on antimicrobial use were obtained from ECDC and IMShealth the following text. For pretreatment, two 10-mL portions of each sample of sewage and the data used in our analyses are provided in supplementary Table 1. IMShealth data supernatant were weighed into separate 50-mL tubes, after which internal stan- were obtained under license and has restricted use. All source data underlying the dards were added. To one of these portions, antimicrobials were added at a level of multivariate analyses, the machine learning and Figs. 1–4 and supplementary Figs. 1–16 25 ng/L for the sulfonamides and 100 ng/L for the tetracyclines, quinolones, and are included in supplementary tables or supplementary data. macrolides. Four mL of EDTA-McIlvain buffer (0.1 M, pH 4.0) were added, after which the samples were shaken for 5 min head-over-head. The residue was taken up in 100 μL MeOH, after which 400 μL of water was added. For liquid chroma- Received: 16 November 2018 Accepted: 31 January 2019 tography tandem mass spectrometry (LC-MS/MS), the following gradient was applied: 0–0.5 min 1% B; 0.5–2.5 min linear increase to 25% B; 2.5–5.4 min linear increase to 70% B; 5.4–5.5 min linear increase to 100% B, with a final hold of 1.0 min. The injection volume was 5 μL. Detection was carried out by an AB Sciex (Ramingham, MA, USA) Q-Trap 5500 or a Q-Trap 6500 mass spectrometer in positive electrospray ionization (ESI). The parameters used for the Q-Trap 5500 References and the Q-Trap 6500 were: capillary voltage, −4.0 kV; declustering potential, 10 V; 1. World Health Organization. Antimicrobial Resistance: Global Report on source temperature, 450 °C; GAS 1 and 2, 50 (arbitrary units). Surveillance. [https://apps.who.int/iris/bitstream/10665/112642/1/ 9789241564748_eng.pdf] (WHO Press, World Health Organization, Geneva, Analysis of aminoglycosides. The analysis of aminoglycosides in the sewage 2014). samples was performed as in Bello53, with adaptations as summarized in the fol- 2. Aarestrup, F. M. The livestock reservoir for antimicrobial resistance: a lowing. For sample pretreatment, two 10-mL portions of each sample of sewage personal view on changing patterns of risks, effects of interventions and the were weighed into separate 50-mL tubes, after which internal standards were way forward. Philos. Trans. R. Soc. Lond. B Biol. Sci. 370, 20140085 (2015). added. To one of these portions, aminoglycosides were added at a level of 50 μg/L. 3. Zignol, M. et al. Twenty years of global surveillance of antituberculosis-drug Twenty mL of extraction liquid (10 mM KH2PO4 with 0.4 mM EDTA and 2% resistance. N. Engl. J. Med. 375, 1081–1089 (2016). TCA) were added, and samples were mixed by means of a vortex and shaken head- 4. Weston, E. J., Wi, T. & Papp, J. Strengthening global surveillance for over-head for 10 min. The extract was then brought to pH 7.6–7.9 and centrifuged antimicrobial drug–resistant Neisseria gonorrhoeae through the Enhanced (15 min, 3600 × g). The complete extract was transferred to a conditioned CBX Gonococcal Antimicrobial Surveillance Program. Emerg. Infect. Dis. 23, cartridge, followed by washing with 4 mL of water and drying. The aminoglyco- https://doi.org/10.3201/eid2313.170443 (2017). sides were eluted with 3 mL of acetic acid (10% in MeOH). The eluate was dried at 5. Hay, S. I. et al. Measuring and mapping the global burden of antimicrobial 60 °C, evaporated under N2 and taken up in 400 μL of HFBA (0.065%). For LC- resistance. BMC Med. 16, 78 (2018). MS/MS, the following gradient was applied: 0–0.5 min, 0% B; 0.5–5 min, linear 6. United Nations, Department of Economic and Social Affairs, Population increase to 45% B; 5–8 min, linear increase to 60% B; 8–10 min, linear increase to Division (2014). World Urbanization Prospects: The 2014 Revision, Highlights 100% B. The injection volume was 40 μL. Detection was carried out by a Waters (ST/ESA/SER.A/352) (WHO Press, World Health Organization, Geneva). (Milford, MA, USA) Quattro Ultima mass spectrometer in positive ESI mode. The 7. UNICEF/WHO. Progress on Sanitation and Drinking Water – 2015 Update parameters used were: capillary voltage, 2.7 kV; desolvation temperature, 500 °C; and MDG Assessment. United Nations. [https://esa.un.org/unpd/wup/ source temperature, 120 °C; cone gas, 150 L/h; and desolvation gas 550 L/h. publications/files/wup2014-highlights.pdf] (accessed April 2018). 8. Baum, R., Luh, J. & Bartram, J. Sanitation: a global estimate of sewerage Analysis of β-lactams. The analysis of β-lactams in the sewage samples was connections without treatment and the resulting impact on MDG progress. performed as in Beredsen et al.54, with adaptations as summarized in the following. Environ. Sci. Technol. 47, 1994–2000 (2013). For sample pretreatment, two 10-mL portions of each sample of sewage were 9. Fernández, M. D. et al. Environmental surveillance of norovirus in Argentina weighed into separate 50-mL tubes, after which internal standards were added. To revealed distinct viral diversity patterns, seasonality and spatio-temporal one of these portions, β-lactams were added at a level of 50 μg/L for the penicillins diffusion processes. Sci. Total Environ. 437, 262–269 (2012). and 500 μg/L for the cephalosporins and carbapenems. Detection was carried out 10. Madico, G. et al. Active surveillance for Vibrio cholerae O1 and vibriophages by a Waters model Xevo TQS or an AB Sciex (Ramingham, MA, USA) Q-Trap in sewage water as a potential tool to predict cholera outbreaks. J. Clin. 6500 mass spectrometer in positive ESI mode. The parameters used for the QTrap Microbiol. 34, 2968–2972 (1996). 6500 were: capillary voltage, 2.0 kV; cone voltage, 25 V; source offset, 20 V; source 11. Asghar, H. et al. Environmental surveillance for polioviruses in the temperature, 150 °C; desolvation temperature, 550 °C; cone gas flow, 150 L/h; and Global Polio Eradication Initiative. J. Infect. Dis. 210 Suppl 1, S294–S303 desolvation gas, 600 L/h. (2014). 12. Hovi, T. et al. Role of environmental poliovirus surveillance in global polio Calculation of defined daily dosages (DDDs) based on residues. After con- eradication and beyond. Epidemiol. Infect. 140, 1–13 (2012). sumption, antimicrobials undergo (1) metabolization in the body, (2) are elimi- 13. Munk, P. et al. A sampling and metagenomic sequencing-based methodology nated from the body with urine and/or feces, and might (3) further degrade in the for monitoring antimicrobial resistance in swine herds. J. Antimicrob. sewer. A proper calculation of DDD/person/day from concentrations in sewage Chemother. 72, 385–392 (2017). NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 14. Petersen, T. N. et al. Meta-genomic analysis of toilet waste from long distance 44. Zhou, Z.-H. Ensemble Methods: Foundations and Algorithms (Chapman and flights; a step towards global surveillance of infectious diseases and Hall/CRC, New York, 2012). antimicrobial resistance. Sci. Rep. 5, 11444 (2015). 45. Breiman, L. Random forests. Mach. Learn. 45, 5–32 (2001). 15. Su, J. Q. et al. Metagenomics of urban sewage identifies an extensively shared 46. Tibshirani, R. Bias, variance and prediction error for classification rules. antibiotic resistome in China. Microbiome 5, 84 (2017). Monogr. Soc. Res. Child Dev. 79, 1–14 (1996). 16. Pehrsson, E. C. et al. Interconnected microbiomes and resistomes in low- 47. Wolpert, D. H. & Macready, W. G. Efficient method to estimate Bagging’s income human habitats. Nature 533, 212–216 (2016). generalization error. Mach. Learn. 35, 41–55 (1999). 17. Petersen, T. N. et al. MGmapper: reference based mapping and taxonomy 48. Breiman, L. Bagging predictors. Mach. Learn. 24, 123–140 (1996). annotation of metagenomics sequence reads. PLoS ONE 12, e0176469 (2017). 49. Breiman, L. Bias, variance and arcing classifiers. Ann. Stat. 26, 801–849 18. Afshinnekoo, E. et al. Geospatial resolution of human and bacterial diversity (1998). with city-scale metagenomics. Cell Syst. 1, 72–87 (2015). 50. Clark, T. E. Can out-of-sample forecast comparisons help prevent overfitting? 19. Roberts, M. C. Update on acquired tetracycline resistance genes. FEMS J. Forecast. 23, 115–139 (2004). Microbiol. Lett. 245, 195–203 (2005). 51. Kursa, M. B. Robustness of random forest-based gene selection methods. BMC 20. Roberts, M. C. Update on macrolide-lincosamide-streptogramin, ketolide, and Bioinformatics 15, 8 (2014). oxazolidinone resistance genes. FEMS Microbiol. Lett. 282, 147–159 (2008). 52. Berendsen, J. A. et al. The analysis of animal faeces as a tool to monitor 21. Schwarz, S., Cloeckaert, A. & Roberts, M. C. in Antimicrobial Resistance in antibiotic usage. Talanta 132, 258–268 (2015). Bacteria of Animal Origin (ed. Aarestrup, F. M.) 73–98 (ASM Press, 53. Bello, G. T. J. et al. Study of the aminoglycoside subsistence phenotype of Washington, D.C., 2006). bacteria residing in the gut of humans and zoo animals. Front. Microbiol. 6, 1550 22. Davies, J. & Davies, D. Origins and evolution of antibiotic resistance. (2016). Microbiol. Mol. Biol. Rev. 74, 417–433 (2010). 54. Berendsen, J. A. et al. Comprehensive analysis of ß-lactam antibiotics 23. van de Sande-Bruinsma, N. et al. Antimicrobial drug use and resistance in including penicillins, cephalosporins, and carbapenems in poultry muscle Europe. Emerg. Infect. Dis. 14, 1722–1730 (2008). using liquid chromatography coupled to tandem mass spectrometry. Anal. 24. Aarestrup, F. M. et al. Effect of abolishment of the use of antimicrobial agents Bioanal. Chem. 405, 7859–7874 (2013). for growth promotion on occurrence of antimicrobial resistance in fecal enterococci from food animals in Denmark. Antimicrob. Agents Chemother. 45, 2054–2059 (2001). Acknowledgements 25. Dutil, L. et al. Ceftiofur resistance in Salmonella enterica serovar Heidelberg We are grateful to Hanne Nørgaard Nielsen, Christina Aaby Svendsen, Jacob Dyring from chicken meat and humans, Canada. Emerg. Infect. Dis. 16, 48–54 (2010). Jensen, Birthe S. Rosenqvist Lund, Kate Vina Vibefeldt, Inge Marianne Hansen, Gunhild 26. Aarestrup, F. M. Characterization of glycopeptide-resistant Enterococcus Larsen, Hanne Mordhorst, and Carsten Bidstrup for technical assistance. We also thank faecium (GRE) from broilers and pigs in Denmark: genetic evidence that Jeffrey Edward Skiby for internal review. This study has received funding from the persistence of GRE in pig herds is associated with coselection by resistance to European Union’s Horizon 2020 research and innovation programme under grant macrolides. J. Clin. Microbiol. 38, 2774–2777 (2000). agreement no. 643476 (COMPARE), the World Health Organization, The Villum 27. Holmes, A. H. et al. Understanding the mechanisms and drivers of Foundation (VKR023052), and The Novo Nordisk Foundation (NNF16OC0021856: antimicrobial resistance. Lancet 387, 176–187 (2016). Global Surveillance of Antimicrobial Resistance). 28. Collignon, P., Powers, J. H., Chiller, T. M., Aidara-Kane, A. & Aarestrup, F. M. World Health Organization ranking of antimicrobials according to their Author contributions importance in human medicine: A critical step for developing risk F.M.A. conceived the idea and ensured the funding. F.M.A., R.S.H., and M.P.K. elabo- management strategies for the use of antimicrobials in food production rated the idea and planned the study. R.S.H. conducted the sampling with assistance animals. Clin. Infect. Dis. 49, 132–141 (2009). from J.A., A.A.-K., S.K., and J.K. and in collaboration with The Global Sewage Surveil- 29. Martínez, J. L., Coque, T. M. & Baquero, F. What is a resistance gene? Ranking lance Project Consortium. S.J.P., M.S., T.Z., A.M.d.R.H., H.S., F.M.A., and P.M. generated risk in resistomes. Nat. Rev. Microbiol. 13, 116–123 (2015). the data and planned the analyses. P.M., P.N., O. Lukjancenko, T.R., B.v.B., L.M.K., D.N., 30. Klein, E. Y. et al. Global increase and geographic convergence in antibiotic R.S.K., P.T.L.C.C., J.K.V., P.L., C.A., G.C., O. Lund, T.H., M.W., H.V., T.S.-P., B.P., S.R., consumption between 2000 and 2015. Proc. Natl. Acad. Sci. USA 115, T.N.P., and F.M.A. performed the bioinformatics, statistical, and epidemiological ana- E3463–E3470 (2018). lysis. F.M.A. wrote the first draft. All authors have seen and accepted the final version. 31. Vogwill, T. & MacLean, R. C. The genetic basis of the fitness costs of antimicrobial resistance: a meta-analysis approach. Evol. Appl. 8, 284–295 (2015). 32. Knudsen, B. E. et al. Impact of sample type and DNA isolation procedure on Additional information genomic inference of microbiome composition. mSystems 1, e00095–16 (2016). Supplementary Information accompanies this paper at https://doi.org/10.1038/s41467- 33. Li, H. Aligning sequence reads, clone sequences and assembly contigs with 019-08853-3. BWA-MEM. ArXiv https://arxiv.org/abs/1303.3997 (2013). 34. Li, H. et al. The Sequence Alignment/Map format SAMtools. Bioinformatics Competing interests: The authors declare no competing interests. 25, 2078–2079 (2009). 35. Zankari, E. et al. Identification of acquired antimicrobial resistance genes. J. Reprints and permission information is available online at http://npg.nature.com/ Antimicrob. Chemother. 67, 2640–2644 (2012). reprintsandpermissions/ 36. Munk, P. et al. Abundance and diversity of the fecal resistome in slaughter pigs and broilers in nine European countries. Nat. Microbiol. 3, 898–908 (2018). Journal peer review information: Nature Communications thanks Jack Gilbert for his 37. Ondov, B. D. et al. Mash: fast genome and metagenome distance estimation contribution to the peer review of this work. Peer reviewer reports are available. using MinHash. Genome Biol. 17, 132 (2016). 38. Global Environment Monitoring System (GEMS). GEMS/Food consumption Publisher’s note: Springer Nature remains neutral with regard to jurisdictional claims in database (2012). (dataset: [https://extranet.who.int/sree/Reports? published maps and institutional affiliations. op=vs&path=/WHO_HQ_Reports/G7/PROD/EXT/ GEMS_cluster_diets_2012]; click on “Show Data” in top-right of world-map, then “Export” in top-right of the table that has now appeared) Open Access This article is licensed under a Creative Commons 39. United Nations Development Programme (UNDP). Human Development Attribution 4.0 International License, which permits use, sharing, Report 2016: Human Development for Everyone (2017). (dataset: [http://hdr. adaptation, distribution and reproduction in any medium or format, as long as you give undp.org/en/composite/trends]) appropriate credit to the original author(s) and the source, provide a link to the Creative 40. Institute for Health Metrics and Evaluation (IHME). Global Burden of Disease Commons license, and indicate if changes were made. The images or other third party Study 2015 (GBD 2015) Location Hierarchies. Seattle, United States: Institute for material in this article are included in the article’s Creative Commons license, unless Health Metrics and Evaluation (IHME) (2017). (dataset: [http://ghdx.healthdata. indicated otherwise in a credit line to the material. If material is not included in the org/record/global-burden-disease-study-2015-gbd-2015-location-hierarchies]) article’s Creative Commons license and your intended use is not permitted by statutory 41. Harrison, X. A. Using observation-level random effects to model regulation or exceeds the permitted use, you will need to obtain permission directly from overdispersion in count data in ecology and evolution. Peer J. 2, e616 (2014). 42. Stekhoven, D. J. & Bühlmann, P. Missforest-non-parametric missing value the copyright holder. To view a copy of this license, visit http://creativecommons.org/ imputation for mixed-type data. Bioinformatics 28, 112–118 (2012). licenses/by/4.0/. 43. Ren, Y., Zhang, L. & Suganthan, P. N. Ensemble classification and regression: recent developments, applications and future directions. IEEE Comput. Intell. © The Author(s) 2019 Mag. 11, 41–53 (2016). 10 NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 ARTICLE The Global Sewage Surveillance project consortium Artan Bego11, Catherine Rees12, Susan Cassar12, Kris Coventry12, Peter Collignon13, Franz Allerberger14, Teddie O. Rahube15, Guilherme Oliveira16, Ivan Ivanov17, Yith Vuthy18, Thet Sopheak18, Christopher K. Yost19, Changwen Ke20, Huanying Zheng20, Li Baisheng20, Xiaoyang Jiao21, Pilar Donado-Godoy22, Kalpy Julien Coulibaly23, Matijana Jergović24, Jasna Hrenovic25, Renáta Karpíšková26, Jose Eduardo Villacis27, Mengistu Legesse28, Tadesse Eguale28, Annamari Heikinheimo29, Lile Malania30, Andreas Nitsche31, Annika Brinkmann31, Courage Kosi Setsoafia Saba32, Bela Kocsis33, Norbert Solymosi34, Thorunn R. Thorsteinsdottir35, Abdulla Mohamed Hatha36, Masoud Alebouyeh37, Dearbhaile Morris38, Martin Cormican38, Louise O’Connor38, Jacob Moran-Gilad39, Patricia Alba40, Antonio Battisti40, Zeinegul Shakenova41, Ciira Kiiyukia42, Eric Ng’eno43, Lul Raka44, Jeļena Avsejenko45, Aivars Bērziņš45, Vadims Bartkevics45, Christian Penny46, Heraa Rajandas47, Sivachandran Parimannan9, Malcolm Vella Haber48, Pushkar Pal49, Gert-Jan Jeunen50, Neil Gemmell50, Kayode Fashae51, Rune Holmstad52, Rumina Hasan53, Sadia Shakoor53, Maria Luz Zamudio Rojas54, Dariusz Wasyl55, Golubinka Bosevska56, Mihail Kochubovski56, Cojocaru Radu57, Amy Gassama58, Vladimir Radosavljevic59, Stefan Wuertz60, Rogelio Zuniga-Montanez60, Moon Y.F. Tay61, Dagmar Gavačová62, Katarina Pastuchova62, Peter Truska62, Marija Trkov63, Kerneels Esterhuyse64, Karen Keddy65, Marta Cerdà-Cuéllar66, Sujatha Pathirage67, Leif Norrgren68, Stefan Örn68, D. G. Joakim Larsson69, Tanja Van der Heijden70, Happiness Houka Kumburu71, Bakary Sanneh72, Pawou Bidjada73, Berthe-Marie Njanpop-Lafourcade74, Somtinda Christelle Nikiema-Pessinaba75, Belkis Levent76, John Scott Meschke77, Nicola Koren Beck77, Chinh Dang Van78, Nguyen Do Phuc78, Doan Minh Nguyen Tran78, Geoffrey Kwenda79, Djim-adjim Tabo80, Astrid Louise Wester81,10 & Sara Cuadros-Orellana82 11 Institute of Public Health, Tirana 1000, Albania. 12Melbourne Water Corporation, Melbourne 3008 VIC, Australia. 13Canberra Hospital, Canberra 2605 ACT, Australia. 14Austrian Agency for Health and Food Safety (AGES), Vienna 1220, Austria. 15Botswana International University of Science and Technology, Private Bag 16, Papapye, Botswana. 16Vale Institute of Technology, Sustainable Development, Belém 66055-090, Brazil. 17National Center of Infectious and Parasitic Diseases, Sofia 1504, Bulgaria. 18Institut Pasteur du Cambodge, Phnom Penh 855, Cambodia. 19University of Regina, Regina S4S0A2, Canada. 20Guangdong Provincial Center for Disease Control and Prevention, No. 160, Guangzhou 511430, China. 21Shantou University Medical College, Shantou 515041, China. 22Corporacion Colombiana de Investigacion Agropecuaria (AGROSAVIA), Mosquera 250040, Colombia. 23Institut Pasteur de Côte d’Ivoire, Abidjan 01 BP 490 Abidjan 01, Côte d’Ivoire. 24Andrija Stampar Teaching Institute of Public Health, Zagreb 10000, Croatia. 25University of Zagreb, Zagreb 10000, Croatia. 26Veterinary Research Institute, Brno 621 00, Czech Republic. 27Instituto Nacional de Investigación en Salud Pública-INSPI (CRNRAM), Quito 170136, Ecuador. 28Addis Ababa University, Addis Ababa P.O. Box 1176, Ethiopia. 29University of Helsinki, Helsinki FI00014, Finland. 30National Center for Disease Control and Public Health, Tbilisi 177, Georgia. 31Robert Koch Institute, Berlin 13353, Germany. 32University for Development Studies, Tamale 233, Ghana. 33Semmelweis University, Institute of Medical Microbiology, Budapest 1089, Hungary. 34University of Veterinary Medicine, Budapest 1078, Hungary. 35Institute for Experimental Pathology, University of Iceland, Keldur 112, Iceland. 36Cochin University of Science and Technology, Cochin 682016, India. 37Foodborne and Waterborne Diseases Research center, Research Institute for Gastroenterology and Liver Diseases, Shahid Beheshti University of Medical Sciences, Tehran 1985717413, Iran. 38National University of Ireland Galway, Galway H91 TK33, Ireland. 39School of Public Health, Ben Gurion University of the Negev and Ministry of Health, Beer Sheva 8410501, Israel. 40Istituto Zooprofilattico Sperimentale del Lazio e della Toscana, Rome 178, Italy. 41National Center of Expertise, Taldykorgan 40000, Kazakhstan. 42Mount Kenya University, Thika P. O. Box 342-01000, Kenya. 43Kenya Medical Research Institute, Nairobi 254, Kenya. 44University of Prishtina “Hasan Prishtina” & National Institute of Public Health of Kosovo, Prishtina 10000, Kosovo. 45 Institute of Food Safety, Riga LV-1076, Latvia. 46Luxembourg Institute of Science and Technology, Belvaux L-4422, Luxembourg. 47Environmental Health Directorate, St. Venera SVR 9018, Malta. 48Agriculture and Forestry University, Kathmandu 44200, Nepal. 49University of Otago, Dunedin 9016, New Zealand. 50University of Ibadan, Ibadan 200284, Nigeria. 51VEAS, Slemmestad NO-3470, Norway. 52Aga Khan University, Karachi 74800, Pakistan. 53National Institute of Health, Lima 15072, Peru. 54National Veterinary Research Institute, Pulawy 24-100, Poland. 55Institute of Public Health of the Republic of Macedonia, Skopje 1000, Republic of Macedonia. 56State Medical and Pharmaceutical University, Chişinău MD- 2004, Republic of Moldova. 57Institut Pasteur de Dakar, Dakar BP 220, Senegal. 58Institute of Veterinary Medicine of Serbia, Belgrade 11000, Serbia. 59Nanyang Technological University, Singapore, Singapore Centre for Environmental Life Sciences Engineering (SCELSE), Singapore 637551, Singapore. 60Nanyang Technological University Food Technology Centre (NAFTEC), Nanyang Technological University (NTU), Singapore 637551, Singapore. 61Public Health Authority of the Slovak Republic, Bratislava 826 45, Slovakia. 62National Laboratory of Health, Environment and Food, Ljubjana 1000, Slovenia. 63Daspoort Waste Water Treatment Works, Pretoria 1, South Africa. 64University of the Witwatersrand, Johannesburg 1485, South Africa. 65IRTA, Centre de Recerca en Sanitat Animal (CReSA, IRTA-UAB), Bellaterra 8193, Spain. 66Medical Research Institute, Colombo 800, Sri Lanka. 67Swedish University of Agricultural Sciences, Uppsala 75007, Sweden. 68The Sahlgrenska Academy at the University of Gothenburg, Gothenburg SE-413 46, Sweden. 69Ara Region Bern Ag, Herrenschwanden 3037, Switzerland. 70Kilimanjaro Clinical Research Institute, Moshi 255, Tanzania. 71National Public Health Laboratories, Ministry of Health and Social Welfare, Kotu Layout, Kotu Layout, Kotu 1863, The Gambia. 72National Institute of Hygiene, Lome 228, Togo. 73Agence de Médecine Préventive, Dapaong, Togo. 74Division of Integrated Surveillance NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-019-08853-3 of Health Emergencies and Response, Lome, Togo. 75Public Health Institution of Turkey, Ankara 6100, Turkey. 76University of Washington, Seattle 98105-6099 WA, USA. 77Institute of Public Health in Ho Chi Minh City, Ho Chi Minh City 700000, Vietnam. 78University of Zambia, Lusaka 15101, Zambia. 79University of N’Djamena, N’Djamena 1117, Chad. 80Norwegian Institute of Public Health, Oslo 0456, Norway. 81Universidad Católica del Maule, Av. San Miguel 3605, Talca, Chile. 82Centro de Biotecnología de los Recursos Naturales, Facultad de Ciencias Agrárias y Forestales, Universidad Católica del Maule, Av. San Miguel 3605, Talca, Chile 12 NATURE COMMUNICATIONS | (2019)10:1124 | https://doi.org/10.1038/s41467-019-08853-3 | www.nature.com/naturecommunications