ARTICLE https://doi.org/10.1038/s41467-020-15123-0 OPEN The phylogenetic landscape and nosocomial spread of the multidrug-resistant opportunist Stenotrophomonas maltophilia Matthias I. Gröschel 1,2, Conor J. Meehan 3, Ivan Barilar 1, Margo Diricks 4, Aitor Gonzaga5, Matthias Steglich5, Oscar Conchillo-Solé 6,7, Isabell-Christin Scherer8, Uwe Mamat9, Christian F. Luz 10, Katrien De Bruyne4, Christian Utpatel1, Daniel Yero 6,7, Isidre Gibert 6,7, Xavier Daura 6,11, 1234567890():,; Stefanie Kampmeier12, Nurdyana Abdul Rahman13, Michael Kresken14,15, Tjip S. van der Werf 2, Ifey Alio16, Wolfgang R. Streit 16, Kai Zhou17,18, Thomas Schwartz 19, John W.A. Rossen 10, Maha R. Farhat 20,21, Ulrich E. Schaible 9,22,23, Ulrich Nübel 5,23,24,25, Jan Rupp 8,22,28, Joerg Steinmann26,27,28, Stefan Niemann1,22,23,28 ✉ & Thomas A. Kohl 1,22,28 Recent studies portend a rising global spread and adaptation of human- or healthcare- associated pathogens. Here, we analyse an international collection of the emerging, multi- drug-resistant, opportunistic pathogen Stenotrophomonas maltophilia from 22 countries to infer population structure and clonality at a global level. We show that the S. maltophilia complex is divided into 23 monophyletic lineages, most of which harbour strains of all degrees of human virulence. Lineage Sm6 comprises the highest rate of human-associated strains, linked to key virulence and resistance genes. Transmission analysis identifies potential outbreak events of genetically closely related strains isolated within days or weeks in the same hospitals. 1 Molecular and Experimental Mycobacteriology, Research Center Borstel, Borstel, Germany. 2 Department of Pulmonary Diseases & Tuberculosis, University Medical Center Groningen, University of Groningen, Groningen, The Netherlands. 3 School of Chemistry and Bioscience, University of Bradford, Bradford, United Kingdom. 4 bioMérieux, Applied Maths NV, Keistraat 120, 9830 St-Martens-Latem, Belgium. 5 Leibniz Institute DSMZ - German Collection of Microorganisms and Cell Cultures, Braunschweig, Germany. 6 Institute of Biotechnology and Biomedicine, Universitat Autònoma de Barcelona, Barcelona, Spain. 7 Department of Genetics and Microbiology, Universitat Autònoma de Barcelona, Barcelona, Spain. 8 Department of Infectious Diseases and Microbiology, University Hospital Schleswig-Holstein, Lübeck, Germany. 9 Cellular Microbiology, Research Center Borstel, Borstel, Germany. 10 Department of Medical Microbiology and Infection Prevention, University Medical Center Groningen, University of Groningen, Groningen, The Netherlands. 11 Catalan Institution for Research and Advanced Studies, Barcelona, Spain. 12 Institute of Hygiene, University Hospital Münster, Münster, Germany. 13 Department of Microbiology, Singapore General Hospital, Singapore, Singapore. 14 Antiinfectives Intelligence GmbH, Rheinbach, Germany. 15 Rheinische Fachhochschule Köln gGmbH, Cologne, Germany. 16 Department of Microbiology and Biotechnology, University of Hamburg, Hamburg, Germany. 17 Shenzhen Institute of Respiratory Diseases, the First Affiliated Hospital (Shenzhen People’s Hospital), Southern University of Science and Technology, Shenzhen, China. 18 Second Clinical Medical College, Jinan University, Shenzhen, China. 19 Karlsruhe Institute of Technology, Institute of Functional Interfaces, Eggenstein- Leopoldshafen, Germany. 20 Department of Biomedical Informatics, Harvard Medical School, Boston, MA, USA. 21 Division of Pulmonary and Critical Care, Massachusetts General Hospital, Boston, MA, USA. 22 German Center for Infection Research (DZIF), partner site Hamburg - Lübeck - Borstel - Riems, Cologne, Germany. 23 Leibniz Research Alliance INFECTIONS’21, Cologne, Germany. 24 Germany Center for Infection Research (DZIF), partner site Hannover - Braunschweig, Cologne, Germany. 25 Braunschweig Integrated Center of Systems Biology (BRICS), Technical University, Braunschweig, Germany. 26 Institute of Medical Microbiology, University Medical Center Essen, Essen, Germany. 27 Medical Microbiology and Infection Prevention, Institute of Clinical Hygiene, Paracelsus Medical Private University, Klinikum Nürnberg, Nuremberg, Germany. 28These authors jointly supervised this work: Jan Rupp, Joerg Steinmann, Stefan Niemann, Thomas A. Kohl. ✉email: sniemann@fz-borstel.de NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications 1 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 R ecently, local transmission and global spread of hospital- countries and four continents, which were either collected and acquired pathogens such as Mycobacterium abscessus and sequenced in this study or had sequence data available in public Mycobacterium chimaera were revealed by whole-genome repositories (Supplementary Data 3). All genome assemblies of sequencing (WGS), thereby challenging the prevailing concepts of the study collection passing quality thresholds (Supplementary disease acquisition and transmission of these pathogens in the Fig. 1, Supplementary Data 4) were analysed with the newly hospital setting1–3. Global genome-based collections are missing created wgMLST scheme. Upon duplicate removal, filtering for for other emerging pathogens such as Stenotrophomonas mal- sequence quality and removal of strains with fewer than 2000 tophilia, listed by the World Health Organization as one of the allele calls in the wgMLST scheme, our study collection leading drug-resistant nosocomial pathogens worldwide4. S. mal- comprised 1305 assembled genomes of majority clinical origin tophilia is ubiquitously found in natural ecosystems, and is of (87%) of which 234 were from public repositories and 1071 newly importance in environmental remediation and industry5,6. S. sequenced strains. Most strains came from Germany (932 strains), maltophilia is an important cause of hospital-acquired drug- the United States (92 strains), Australia (56 strains), Switzerland resistant infections with a significant attributable mortality rate in (49 strains) and Spain (42 strains) (Fig. 1d; Supplementary immunocompromised patients of up to 37.5%7. Patients under Data 3). WgMLST analysis resulted in an average of 4174 (range immunosuppressive treatment and those with malignancy or pre- 3024–4536) loci recovered per strain (Supplementary Data 5). existing inflammatory lung diseases such as cystic fibrosis are at Across the 1305 strains, most loci, 13,002 of 17,603, were assigned particular risk of S. maltophilia infection8. Although almost any fewer than 50 different alleles (Supplementary Fig. 2). Calculation organ can be affected, mere colonisation needs to be discriminated of the sample pan genome yielded 17,479 loci, with 2844 loci from infections that mainly manifest as respiratory tract infec- (16.3%) present in 95% and 1275 loci (7.3%) present in 99% of tions, bacteraemia or catheter-related bloodstream infections5. strains (Supplementary Fig. 3A). The pan genome at scale is not Yet, the bacterium is also commonly isolated from wounds and, at structured, likely due to extensive horizontal gene transfers25 lower frequency, in implant-associated infections9,10. Further- (Supplementary Fig. 3B). The genome sizes ranged from 4.04 Mb more, community acquired infections have also been described11. to 5.2 Mb. Treatment options are limited by resistance to a number of antimicrobial classes such as most β-lactam antibiotics, cepha- S. maltophilia complex comprises 23 monophyletic lineages. To losporins, aminoglycosides and macrolides through the intrinsic investigate the global diversity of the S. maltophilia complex, a resistome, genetic material acquired by horizontal transfer, as well maximum likelihood phylogeny was inferred from a concatenated as non-heritable adaptive mechanisms12,13. sequence alignment of the 1275 core loci present in 99% of the To date, no large-scale genome-based studies on the population 1305 S. maltophilia strains of our study collection (Fig. 1a). structure and clonality of S. maltophilia in relation to human Hierarchical Bayesian analysis of population structure (BAPS), disease have been conducted. Previous work indicated the pre- derived from the core single-nucleotide polymorphism (SNP) sence of at least 13 lineages or species-like lineages in the S. results, clustered the 1305 genomes into 23 monophyletic lineages maltophilia complex, defined as S. maltophilia strains with 16S named Sgn1–Sgn4 and Sm1–Sm18, comprising 17 previously rRNA gene sequence similarities >99.0%, with nine of these suggested and six hitherto unknown lineages (Sm13–Sm18). For potentially human-associated14–18. These S. maltophilia complex consistency, we used and amended the naming convention of lineages are further divided into four more distantly related lineages from previous reports14,16. In concordance with these lineages (Sgn1-4) and several S. maltophilia sensu lato and sensu studies14,16, we found a clear separation of the more distantly stricto lineages14,19. The S. maltophilia strain K279a, isolated from related lineages Sgn1–Sgn4 and a branch formed by lineages a patient with bloodstream infection, serves as an indicator strain Sm1–Sm18 (previously termed S. maltophilia sensu lato), with the of the lineage S. maltophilia sensu stricto19. largest lineage Sm6 (also known as S. maltophilia sensu stricto) To understand the global population structure of the S. mal- containing most strains (n = 413), including the strain K279a and tophilia complex and the potential for global and local spread of the species type strain ATCC 13637. Contrary to previous analyses, strains, in particular of human-associated lineages, we performed Sgn4 is the lineage most distantly related to the rest of the strains14. a large-scale genome-based phylogenetic and cluster analysis of a The division into the 23 lineages is also clearly supported by an global collection of newly sequenced S. maltophilia strains toge- average nucleotide identity (ANI) analysis (Fig. 1b, c). ANI com- ther with publicly available whole-genome data. parisons of strains belonging to the same lineage was above 95%, and comparisons of strains between lineages were below 95%. To evaluate structural genomic variation across the various Results lineages, we compiled a set of 20 completely closed genomes Strain collection and gene-by-gene analysis. To allow for stan- covering the 15 major phylogenetic lineages of both environ- dardised WGS-based genotyping and gene-by-gene analysis of our mental and human-invasive or human-non-invasive isolation data set, we first created an S. maltophilia complex whole-genome source. These genomes were either procured from the NCBI (n = multilocus sequence typing (wgMLST) scheme. This approach, 8) or newly sequenced on the PacBio platform (n = 12) implemented as core genome MLST, has been widely used in tracing (Supplementary Table 2). Interestingly, no plasmids were outbreaks and transmission events for a variety of bacterial spe- detected in any of the genomes. A genome-wide alignment of cies20–22. The use of a wgMLST scheme allows to analyse sequenced the 20 genomes demonstrated considerable variation in both strains by their core and accessory genome23. Using 171 publicly structure and size between strains of different lineages and even available assembled genomes of the S. maltophilia complex that strains of the same lineage (Supplementary Fig. 4). Several phage- represent its currently known diversity (Supplementary Data 1), we related, integrative and conjugative mobile elements were constructed a wgMLST scheme consisting of 17,603 loci (Supple- observed across the genomes. mentary Data 2). To ensure compatibility with traditional MLST/ gyrB typing methods, the wgMLST scheme includes the partial sequences of the seven genes used in traditional MLST as well as the Delineation of the S. maltophilia complex within its genus. We gyrB gene24 (Supplementary Table 1 and Supplementary Data 2). calculated a phylogenetic tree based on an alignment of 23 pre- To investigate the global phylogeographic distribution of S. dicted amino acid sequences from reference genes15 to visualise maltophilia, we gathered WGS data of 2389 strains from 22 the clade formed by strains of the S. maltophilia complex within 2 NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 ARTICLE a Tree scale: 0.01 Sm3 Sm2 Sm4b Sm4a Sgn3 Sm12 Sm14 Sm11 Sgn1 Sm13 Sgn4 Sm1 Sgn2 Sm10 Sm8 Sm7 S. maltophilia sensu lato Sm15 Sm16 S. maltophilia sensu stricto Sm9 Sm5 Sm17 Sm6 Sm18 b 100% c 0.7 0.6 87% 0.5 Density 0.4 0.3 0.2 0.1 0.0 88 90 92 94 96 98 100 Average nucleotide identity (%) d Lineages Sgn1 Sm8 Sgn2 Sm9 Sgn3 Sm10 Sgn4 Sm11 Sm1 Sm12 Sm2 Sm13 Sm3 Sm14 Sm4a Sm15 Sm4b Sm16 Sm5 Sm17 Sm6 Sm18 Sm7 Ungrouped Isolates count 1–10 10–50 50–100 100–1000 1000 or more None Fig. 1 The global population structure of the S. maltophilia complex is composed of 23 monophyletic, globally distributed lineages. a Unrooted maximum likelihood phylogenetic tree of 1305 S. maltophilia strains displaying the known population diversity of the S. maltophilia complex. The tree was built using RAxML on the sequences of 1275 concatenated core genome genes. Groups as defined by hierarchical Bayesian clustering are marked with shaded colours, and group numbers are indicated at the tree leafs of each corresponding group. orange shading = S. maltophilia sensu stricto; green shading = S. maltophilia sensu lato; 100% support values for the main branches are indicated with red circles. b Pairwise average nucleotide identity comparison calculated for 1305 S. maltophilia strains shown on a heatmap with blue indicating high and red indicating low nucleotide identity. c Histogram of pairwise average nucleotide identity (ANI) values, illustrating that strains of the same lineage are highly similar at the nucleotide level with ANI values above 95% (depicted in blue). Inter-lineage comparisons (in red colour) reveal low genetic identity between strains. The currently accepted species delimitation threshold at 95% is shown as a grey vertical line. d Geographic origin of the 1305 S. maltophilia strains comprising the study collection indicated on a global map. The green/ yellow colour code indicates the number of strains obtained per country. The distribution of phylogenetic lineages per continent is displayed as colour-coded pie charts. The map was created using the tmap package in R69. Source data are provided as Source Data Files. the genus Stenotrophomonas (Supplementary Fig. 5). When using only 982 loci were detected. Interestingly, the 16S rRNA gene the wgMLST scheme at the genus level, we recovered between 380 sequence of JCM9942 matched to that of S. acidaminiphila, and loci in S. dokdonensis and a maximum of 1677 in S. rhizophila, the JCM9942 16S rRNA sequence is only 97.3% identical and has with S. terrae, S. panacihumi, S. humi, S. chelatiphaga, S. dae- an average nucleotide identity (ANI) of 82.9% with that of S. jeonensis, S. ginsengisoli, S. koreensis and S. acidaminiphilia spe- maltophilia ATCC 13637. We note that this strain has been cies receiving allele calls between these two values. For the strain reclassified as S. pictorum per 21/12/2019. In contrast, the num- S. maltophilia JCM9942 (Genbank accession GCA_001431585.1), ber of recovered loci matches those of S. maltophilia strains for NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications 3 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 a b c Isolation source Sm18 Sm18 Sm17 Sm17 Anthropogenic Sm16 * Sm18 Sm16 Environmental Sm15 Sm15 Sm17 Sm16 Sm14 Sm14 Human–invasive Sm15 Phylogenetic lineage Sm13 Sm13 Sm14 Human–non–invasive Phylogenetic lineage * Phylogenetic lineage Sm13 Sm12 Sm12 Sm12 Sm11 Sm11 Percent * * Sm11 Human–respiratory Sm10 Sm10 Sm10 Sm9 10.0% Sm9 Sm9 * indicates a one-sided Sm8 Sm7 20.0% Sm8 * Sm8 p value of < 0.05 Sm6 Sm7 Sm7 Sm5 Sm4b 30.0% Sm6 * Sm6 Sm4a Sm5 Sm5 Sm3 Sm2 Sm4b * Sm4b Sm1 Sm4a * Sm4a Sgn4 Sm3 Sm3 Sgn3 Sgn2 Sm2 * Sm2 Sgn1 Sm1 Sm1 Asia Australia Europe North America Sgn4 Sgn4 Sgn3 * Sgn3 Sgn2 Sgn2 * Sgn1 * Sgn1 0 100 200 300 0 50 100 150 200 Number of strains (n) Number of strains (n) Fig. 2 Global distribution of lineages, their composition by isolation source and contribution to phylogenetic lineage total number of strains. a Bubble plot illustrating the proportion of lineages per continent. b Barplot showing the number of strains per lineage coloured by isolation source for the entire strain collection (see colour legend in c). c Barplot for the prospectively sampled representative collection of human-associated S. maltophilia complex strains per lineage coloured by human-invasive, human-non-invasive or human-respiratory. One-sided p-values for all within-lineage comparison of isolation sources can be found in Supplementary Table 3. n = 1179 S. maltophilia isolates where isolation source was known. ‘*' indicates a p-value of < 0.05 using one-sided Fisher’s exact test (for n < 5) or test of equal and given proportions corrected for multiple testing using the Benjamini–Hochberg procedure. P. geniculata (4060 loci), S. lactitubi (3805 loci), S. pavanii (3623 respectively). Sgn3 contained only few human-non-invasive loci), and S. indicatrix (3678 loci). Here, 16S rRNA sequence strains (8%, p = 0.02). Strains of lineages Sm6 (52%, p < 0.001), comparison to S. maltophilia ATCC 13637 would support the Sm2 (65%, p = 0.04) and Sm13 (74%, p = 0.03) were linked to the inclusion of the first three of the aforementioned species with 16S human-respiratory isolation source. Strains of Sgn3 (11%, p < rRNA sequence identity of >99.1% into the S. maltophilia 0.001) and Sm11 (10%, p < 0.001) were less likely to be isolated complex18,26. from the human-respiratory tract (Fig. 2b). The majority of strains sequenced within this study were prospectively collected through a hospital consortium across Lineages are globally distributed and differ in human asso- Germany, Austria and Switzerland (n = 741) (Fig. 2c). Restricting ciation. We next analysed the global distribution of strains of the to this collection of human-associated strains from hospitalised lineages defined above and found that eight (Sm2, Sm3, Sm4a, patients, the most common lineage was Sm6 (33%) and lineages Sm6, Sm7, Sm9, Sm10 and Sm12) are represented on all con- Sgn1-3 and Sm11, found to be environmentally associated in tinents sampled within this study, with strains of lineage Sm6 public data, were either not present or represented a minor accounting for the largest number of strains globally and the proportion (<0.1%) of the collection. largest proportion on each sampled continent (Figs. 1d and 2a). We attempted genome-wide association (GWAS) to investigate To further investigate whether the lineages correlate with isola- the genetic correlates of human niche specificity using an elastic tion source, particularly with regard to human host adaptation, net whole-genome model that has recently been shown to we classified the isolation source of the S. maltophilia strains into outperform univariate approaches in controlling for population five categories. Strains were considered environmental (n = 117) structure27. Using this approach, there was still considerable if found in natural environments, e.g. in the rhizosphere, and confounding due to residual population structure, likely related to anthropogenic if swabbed in human surroundings such as patient the strong niche specificity of some of the lineages described room sink or sewage (n = 52). Human-invasive (n = 133) was above (Supplementary Fig. 8). used for isolates from blood, urine, drainage fluids, biopsies or in cerebrospinal fluid, human-non-invasive (n = 353) refers to colonising isolates from swabs of the skin, perineum, nose, oro- Resistome and virulence characteristics of S. maltophilia. We pharynx, wounds as well as intravascular catheters, and human- next screened our collection to detect potential resistance genes, respiratory (n = 524) includes strains from the lower respiratory e.g. chromosomally encoded antibiotic resistance genes, tract below the glottis and sputum collected from cystic fibrosis including efflux pumps19,25,28. We could identify members of patients. For 126 strains, no information on their isolation source the five major families of efflux transporters with high frequency was available, and thus, these were not included in this analysis. in our strain collection (Fig. 3a)19,29. Aminoglycoside-modifying The more distantly related lineages Sgn1 (100%), Sgn2 (100%), enzymes were encoded in 6.1% of strains (aminoglycoside- Sgn3 (76%) and also Sm11 (38%) contained significantly more acetyltransferases) and 66% of strains (aminoglycoside-phos- environmental strains compared with all other isolation sources photransferases), respectively, with five strains also harbouring in this lineage (p < 0.001, one-sided test of equal or given aminoglycoside-nucleotidyltransferases. We observed that these proportions or Fisher’s exact test for n < 5, corrected for multiple enzyme families were unequally distributed among lineages, testing using the Benjamini–Hochberg procedure), whereas which preferentially contained either of the two major types. strains of lineages Sm4a and Sm6 (3% and 5%, p < 0.001) were Strains of lineage Sm4a had the lowest proportion (4%, p < minority environmental (Fig. 2b, c; Supplementary Table 3). 0.001) of aminoglycoside-phosphotransferases. Taken together, Anthropogenic strains were found at higher proportions in 69% of the strains of our collection featured aminoglycoside- lineages Sm11 and Sm12 (24% and 19%, p < 0.001). Strains of modifying enzymes. Other enzymes implicated in aminoglyco- lineage Sm4b were likely to be classified as human-invasive (33%, side resistance are the proteases ClpA and HtpX that p = 0.02), and strains of lineage Sm4a and Sm8 were more likely were present in 99.3% and 99.8% of the strains investigated, to be human-non-invasive (42%, p < 0.001 and 60%, p = 0.03, respectively30. The S. maltophilia K279a genome encodes two 4 NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 ARTICLE a β-lactamases b aph Aminp- Group Efflux pumps glycosides Sulfonamides Virulence 0.6 E) rM M ) su E ( S) aa (M R) no (S R AT gE SM F te U2 e e p em (M em CG s m -t y k a r1 pi R bl 1 su 2 rA P D Tree scale: 0.1 tA pX sm B tA A aL aL o Dim2 (14%) l1 l2 r lU h c m RN t c lp 0.4 ap su ca ht bl st pilU Sgn4 Sgn2 Sgn1 Sgn3 Sm1 Sm13 0.2 Sm12 blaL2 Sm14 aac katA origin_detailed sul1 blaL1 0.0 region qacE smoR Sm3 0.00 0.25 0.50 0.75 Sm4b Dim1 (28.8%) Sm2 c 3 Sm4a 2 Lineage Dim2 (14%) Sgn1 Sm17 Sgn2 Sm18 Sm11 Sgn3 Sm2 Sgn4 Sm3 Sm1 Sm4a Sm10 1 Sm10 Sm4b aph_1 Sm11 Sm5 Sm8 Sm12 Sm6 Sm13 Sm7 katA_0 Sm14 Sm8 Sm7 blaL2_1 blaL1_1 Sm15 Sm9 smoR_0 Sm16 Ungrouped smoR_1 katA_1 blaL1_0 0 Sm15 Sm16 blaL2_0 Sm9 aph_0 –1 Sm5 Sm17 –1 0 1 2 3 Sm18 Dim1 (28.8%) d aph_1 0.5 katA_0 blaL2_1 Sm6 blaL1_1 smoR_0 Dim2 (14.2%) Isolation_source smoR_1 0.0 blaL1_0 katA_1 Anthropogenic Environmental Human–invasive Human–non–invasive Human–respiratory blaL2_0 –0.5 aph_0 –1.0 –0.5 0.0 0.5 1.0 1.5 2.0 2.5 Dim1 (28.5%) Fig. 3 Resistance and virulence gene analysis. a Midpoint rooted maximum likelihood phylogenetic tree based on 1275 core gene sequences of 1305 S. maltophilia complex strains. The coloured shading of the lineages represents the groups found by Bayesian clustering, with lineage names given. Hundred percent branch support is indicated by red dots. The pattern of gene presence (blue coloured line) or absence (white) is displayed in columns next to the tree, showing, from left to right, selected efflux pump genes: resistance-nodulation-cell-division (RND)-type efflux pumps, smeU2 as part of the five-gene RND efflux pump operon smeU1-V-W-U2-X, tetACG, emrA of the major facilitator superfamily (MFS), emrE and sugE of the small-multidrug-resistance (SMR) efflux pump family, norM of the MATE family; the aminoglycoside acetyltransferase aac and phosphotransferase aph, clpA, htpX, the β-lactamases blaL1 and blaL2, the sul1 and sul2 genes encoding dihydropteroate synthases, catB, and the virulence genes smoR, pilU, stmPr1 and katA. b Variable correlation plot of a multiple correspondence analysis (MCA) visualising nine resistance and virulence genes as active variables in red and three supplementary variables region, origin and groupin blue. c Factor individual biplot map of phylogenetic lineages, indicated by their 99% confidence intervals (ellipses) across the first two MCA dimensions. The five highest contributing active variables are shown in red with 0 denoting absence and 1 presence of this variable. d Factor individual biplot map of the isolation source as indicated in the coloured legend. Source data are provided as Source Data file. β-lactamases, the metallo-β-lactamase blaL1 and the inducible harbouring the type B chloramphenicol-O-acetyltransferase Ambler class A β-lactamase blaL231. While blaL1 was found in CatB (0.6%). The sulfonamide-resistance-conferring sul1 was seen 83.2% of our strains, blaL2 was detected in only 63.2%. Inter- in 19 strains (1.4%), and sul2 was found in only 5 strains (0.4%), estingly, strains of some lineages lacked the blaL2 gene, i.e. Sgn4, mostly occurring in human-associated or anthropogenic strains. Sm1, Sm12, Sm13 and Sm16. Sm4a was the only lineage where This hints towards a low number of trimethoprim/sulfamethox- no blaL1 was found. Only one isolate encoded the oxacillin azole-resistant strains in our collection, which is the recommended hydrolysing class D β-lactamase OXA. We noted a few strains first-line agent for the treatment of S. maltophilia infection12. NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications 5 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 We investigated the presence of virulence genes in our (Fig. 4a). These thresholds were chosen based on the distribution collection. SmoR is involved in quorum sensing and swarming of allelic mismatches (Fig. 4b). We found 765 (63%) strains to motility of S. maltophilia, and was observed in 89.3% of our group into 82 clusters (median cluster size 6, IQR 6–11.7) within strains32. While SmoR was present in all Sm6 strains (proportion 100 alleles difference. A total of 270 (21%) strains were grouped of 100%, p < 0.001, test of equal or given proportions, corrected into 62 clusters within 10 alleles difference (median cluster size 4, for multiple testing using the Benjamini–Hochberg procedure), IQR 3–4.7). The maximum number of strains per cluster were 45 this gene was less prevalent in strains of lineages Sgn1 (7%, p < and 12 for the d100 and d10 clusters, respectively. Interestingly, 0.001), and absent in strains of Sgn2, Sgn3 and Sgn4 (all p < strains within d100 clonal complexes originated from different 0.001). PilU, a nucleotide-binding protein that contributes to countries or cities (Fig. 5a). Type IV pilus function, was found in 9% of strains and mainly in Some strains of lineages, notably those with primarily lineages Sm9 (proportion of 81%, p < 0.001,) and Sm11 (propor- environmental strains, did not cluster at d10 level at all (Sgn1- tion of 50%, p < 0.001)33. StmPr1 is a major extracellular protease 4, Sm1, Sm15) (Supplementary Table 4). The d10 clustering rate of S. maltophilia and is present in 99.8% of strains34. KatA is a ranged from 18% for strains of lineage Sm4a to 48% for Sm13, catalase mediating increased levels of persistence to hydrogen while 21% of lineage Sm6 strains were in d10 clusters. When day peroxide-based disinfectants and was found in 86.6% of strains35. and location of isolation were known and included in further While KatA is present at high proportions in strains of most investigations of the d10 clusters, we detected a total of 49 strains, lineages (i.e. Sm6 with 99% or Sm4a with 99%), Sm3 harbours grouped into 13 clusters (of at least two isolates), which were this gene in only 49% (p < 0.001) of its strains, whereas it is absent isolated from the same respective hospital in the same year. Of in Sgn1, Sgn2, Sgn3 and Sgn4 (all p < 0.001). Taken together, S. these, three d10 clusters consisted of strains isolated from the maltophilia strains harbour a number of resistance-conferring as respiratory tract of different patients treated in the same hospital well as virulence genes, some of which are unequally distributed within an 8-week time span or less (Table 1, Fig. 5b; over the lineages. Supplementary Fig. 7). We used multiple correspondence analysis (MCA) to investi- gate the correlation of the resistance and virulence profiles of the strains with geographic origin, isolation source and phylogenetic Discussion lineage. A total of nine genes, derived from virulence The findings of this study demonstrate that strains of the human databases36,37, that were either present or absent in at least ten opportunistic pathogen S. maltophilia can be subdivided into 23 isolates were selected to serve as active variables for the MCA monophyletic lineages, with two of these comprising exclusively (aac, aph, blaL1, blaL2, katA, pilU, qacE, smoR, sul1). As expected environmental strains. The remaining lineages contain strains with a complex data set, the total variance explained by the MCA from mixed environmental and human sources. Among these model was relatively low with the first four dimensions explaining strains, certain lineages such as Sm6 are most frequently found to 65.2% of variance in the model (Supplementary Fig. 6A). be human-invasive, human-non-invasive or human-respiratory Nevertheless, from examining the first two dimensions of the strains, pointing towards a potential adaptation to human MCA (accounting for 28.8% and 14% of variance), we noted that infection and enhanced virulence. This is supported by their the genes smoR, katA, blaL1, blaL2 and aac correlate with the association with antibiotic resistance genes, resulting in the first dimension of the MCA, while genes aph and pilU are multidrug resistance observed among human-associated lineages. corresponding to the second dimension (Fig. 3b; Supplementary Our data provide evidence for the global prevalence of particular Fig. 6B). When introducing geographic origin, isolation source circulating lineages with hospital-linked clusters collected within and phylogenetic lineage as supplementary variables to the model, short time intervals suggesting transmission. The latter empha- we observed a strong correlation of phylogenetic lineages with sises the need to instate or re-enforce hygiene and infection both dimensions, whereas little to no correlation was observed for control practices to minimise in-hospital spread of these isolation source and geographic origin (Fig. 3b). This indicates pathogens. that virulence and resistance profiles of the nine genes are largely In line with previous reports, our large genome-based study lineage-specific, with little impact of geographic origin or revealed that the S. maltophilia complex is extraordinarily diverse isolation source. However, we found a clear separation of the at the nucleotide level, representing a challenge for population- environmental strains from the rest of the collection when wide analyses and molecular epidemiology14–17. To address this, analysing the impact of human versus environmental habitat on we first developed a new genome-wide gene-by-gene typing the observed variance (Fig. 3d). A more detailed analysis of the scheme, consisting of 17,603 gene targets, or loci. This whole- observed lineage-specific variation, based on the explained genome MLST typing scheme provides a versatile tool for variance from the MCA analysis, reveals that the more distantly genome-based analysis of S. maltophilia complex strains and a related lineages Sgn1-4, Sm1 and Sm13 are characterised by the unified nomenclature to facilitate further research on the complex lack of smoR, katA and the presence of the aminoglycoside- with an integrative genotyping tool and sequence data analysis acetyltransferases aac. The human-associated lineages Sm2, Sm6 approach. Including the loci of the 7-gene classical MLST typing and Sm7 are associated with the presence of blaL2, sul1, blaL1, scheme as well as the gyrB gene enables backward compatibility smoR and katA (Fig. 3c). and comparison of allele numbers with sequence types obtained through the classical MLST scheme24. Applying the wgMLST approach to our extensive and geographically diverse collection of Possible local spread derived from genetic diversity analysis. S. maltophilia strains allowed us to infer a comprehensive phy- The identification of widely spread clonal complexes or potential logenetic population structure of the S. maltophilia complex, outbreak events of S. maltophilia complex strains would have including the discovery of six previously unknown lineages in significant implications for preventive measures and infection addition to those described previously14,17. control of S. maltophilia in clinical settings. We assessed our Altogether, we found 23 distinct phylogenetic lineages of the S. strain collection for circulating variants and clustered strains maltophilia complex, which are well supported by hierarchical using the 1275 core genome MLST loci, that were also used for Bayesian clustering analysis of the core genome and intra- and phylogenetic inference, and thresholds of 100 (d100 clusters) and inter-lineage average nucleotide identity. This genetic hetero- 10 mismatched alleles (d10 clusters) for single-linkage clustering geneity observed between the detected lineages is sufficient to 6 NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 ARTICLE a Tree scale: 0.01 Isolation source Environmental Anthropogenic Human-associated Unknown Sm2 b Detailed isolation source Sm4 Sm 4a Environmental 3 Anthropogenic Sm Human-invasive 14 Human-non-invasive Sm Sm Human-respiratory 11 Unknown 12 Sm City Sm Innsbruck (AUT) 10 13 1 Sm Antwerp (BEL) Sm Berlin (GER) 3 Sm8 Sgn Essen (GER) Frankfurt a.M. (GER) Sgn1 Groningen (NL) Homburg/Saar (GER) Kiel (GER) Sgn4 Sgn2 Sm7 Cologne (GER) Lübeck (GER) Madrid (ESP) b 0.015 Mainz (GER) 0.2 Density Marburg (GER) Sm16 15 0.1 Sm Munich (GER) 0.010 0.0 Münster (GER) Density 0 1000 2000 Nr. different alleles Paris (FR) Regensburg (GER) Seattle (USA) 0.005 Singapour (SIN) 9 Sm Year 0.000 2009 0 50 100 150 200 2010 5 2011 Nr. different alleles Sm 17 Sm 2012 2013 18 Sm Sm 2014 6 2018 Clustering based on mismatched alleles 100 alleles (d100) 10 alleles (d10) Fig. 4 Spatiotemporal cluster analysis of 1305 S. maltophilia complex strains. a The coloured ranges across the outer nodes and branches indicate the 23 lineages. The black dots indicate the location of the genome data sets used for wgMLST scheme generation. The rings, from inside towards outside denote (i) the isolation source of the strains classified as either environmental, anthropogenic, human or unknown; (ii) the detailed isolation source of strains similar to the first ring with the human strains subclassified into human-invasive, human-non-invasive and human-respiratory; (iii) the city of isolation; (iv) the year of isolation (where available), with light colours representing earlier years and darker brown colours more recent isolation dates. The outer rings in black-to-grey indicate the single-linkage-derived clusters based on the number of allelic differences between any two strains for 100 (d100 clusters) and 10 (d10 clusters) allelic mismatches. Red dots on the nodes indicate support values of 100%. b Distribution of the number of wgMLST allelic differences between pairs of strains among the 1305 S. maltophilia strains. The main figure shows the frequencies of up to 200 allelic differences, while the inset displays frequencies of all allelic mismatches. Source data are provided as Source Data Files. consider them as clearly separate lineages of the S. maltophilia potentially switch back and forth between surviving in the envir- complex, in line with previous results from classical typing onment and within a human host. These results do not support the methods and phylogenetic studies14–17,38. In fact, the average notion that the S. maltophilia sensu stricto strains of lineage Sm6 nucleotide identity between lineages was below the threshold represent the primary human pathogens14. We therefore propose generally considered to define a species, warranting further studies to continue using the term S. maltophilia complex and the on and possible revisions of the taxonomic assignments and respective lineage classification for all strains that are identified as nomenclature for this group. In parallel with these reports, human S. maltophilia by routine microbiological diagnostic procedures in adaptation is observed to vary, with strains from lineages Sgn1, hospitals and omit the use of sensu stricto or lato. Sgn2, Sgn3 and Sm11 mostly isolated from the environment and Beyond associating some of the lineages with either environ- strains from the other lineages mostly derived from human or ment or human, we were not able to identify the specific genetic human-associated sources. Apart from the purely environmental mechanisms that underlie this association due to the extent of lineages Sgn1 and Sgn2, our results indicate that strains from all stratification by population structure. S. maltophilia is believed to other lineages are able to colonise humans and cause infection, be a much less virulent pathogen relative to other nosocomials including lineage Sgn4 outside the “sensu lato” group, and such as P. aeruginosa or S. aureus39. The establishment of human NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications 7 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 a Sm6 b PEG-267 Cluster 42 (hospital A) Cluster 45 (hospital B) PEG-328, PEG-329 PEG-258 5 3 1 PEG-263, PEG-268 1 2 Country City PEG-331 Germany Innsbruck (AUT) Australia Antwerp (BEL) PEG-257 PEG-266 Austria Berlin (GER) Belgium Essen (GER) China Frankfurt a.M. (GER) France Groningen (NL) Cluster 47 (hospital C) Cluster 52 (hospital D) India Homburg/Saar (GER) PEG-351 943974Y, 944632W Kiel (GER) Italy Cologne (GER) Mexico Lübeck (GER) Singapore Madrid (ESP) Slovenia Mainz (GER) Spain Marburg (GER) Switzerland Munich (GER) The Netherlands Münster (GER) United Kingdom Paris (FR) USA Regensburg (GER) 1 Seattle (USA) 944796D Singapour (SIN) 5 PEG-349, PEG-350 2 ou r ry ity C uste nt C 1 cl 00 d1 PEG-353 945570W Fig. 5 Analysis of d100 clusters in lineage Sm6 and closely related d10 clusters across the study collection. a The d100 clusters in the largest human- associated lineage Sm6 consist of strains from various countries and, for strains from the same country, of various cities. The coloured bars represent, from left to right, the d100 clonal complexes, the country of isolation, and the city of isolation. b High-resolution analysis of four selected d10 allele clusters for which detailed metadata, i.e. day, source and ward of isolation, was available are shown as minimum spanning trees based on the 100% core genome MLST loci of the respective cluster. The number of loci used were 3734 for cluster 42 (hospital A), 4190 for cluster 45 (hospital B), 3637 for cluster 47 (hospital C) and 3714 for cluster 52 (hospital D). The number of mismatched alleles are shown in small numbers on the connecting branches. Node colours indicate isolation source, light blue = respiratory sample, dark blue = sputum, grey = wound swap, green = endoscope. Source data are provided as Source Data Files. Table 1 Site and date of isolation for the strains comprising the four d10 clusters isolated from the same geographic location within at most an 8-week time span. Strain Lineage Cluster Isolation date Isolation place Clinical source PEG-257 Sm2 42 October 4th, 2013 Hospital A Respiratory tract PEG-258 October 7th, 2013 Wound swap PEG-263 October 9th, 2013 Respiratory tract PEG-266 October 14th, 2013 Respiratory tract PEG-267 October 14th, 2013 Respiratory tract PEG-268 October 14th, 2013 Sputum PEG-328 Sm18 45 October 21st, 2013 Hospital B Respiratory tract PEG-329 October 21st, 2013 PEG-331 December 9th, 2013 PEG-351 Sm13 47 January 11th, 2014 Hospital C Sputum PEG-349 January 21st, 2014 Respiratory tract PEG-350 January 24th, 2014 Respiratory tract PEG-353 January 27th, 2014 Respiratory tract 943974Y Sm12 52 January 28th, 2014 Hospital D Endoscope 944632W February 17th, 2014 944796D February 21st, 2014 945570W March 24th, 2014 infection or colonisation with S. maltophilia is likely also strongly host immune status or other predisposing factors will enable driven by the host factors such as the immune status while the research in this area in the future. role of pathogen genetic background or specific virulence Our results further illustrate that strains of nearly all 23 mechanisms is still to be determined40. Collecting data on the lineages are present in sampled countries and continents, 8 NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 ARTICLE suggesting a long evolutionary trajectory of S. maltophilia. The prevention and control teams deal with S. maltophilia colonisa- finding that the more distantly placed lineages (Sgn1-4) as well as tion or infection. the other species of the genus Stenotrophomonas comprise pri- We are aware that our study is limited by our collection fra- marily environmental strains lends to speculations that this tra- mework. Molecular surveillance of S. maltophilia is currently not jectory took place from an exclusively environmental lifestyle routinely performed and no robust data on prevalence, sequence towards human colonisation and infection. This could be due to types or resistance profiles exist. The geographic restriction of our the emergence of individual strains adapted to survive in both prospective sampling is biased towards the acquisition of clinical niches or to multiple, independent events of pathoadaptation of and human-pathogenic S. maltophilia strains from a multi- environmental strains to human colonisation, as has been national consortium that mainly comprised German, Austrian observed for Legionella pneumophila41. A more recent study and Swiss hospitals. The inclusion of all available sequence data in expanded these findings on the entire Legionella genus, illus- public repositories compensates this restriction partially, how- trating that the capacity to infect eukaryotic cells can be acquired ever, for these strains information on isolation source and date independently many times42. The evolution within the S. mal- was incomplete or missing. More prospective, geographically tophilia complex might have been aided by the apparent genomic diverse sampling from different habitats is warranted to corro- plasticity as seen from quite distinct genome lengths and struc- borate our findings, especially concerning the apparent adapta- tural variation, even within individual lineages. In addition, tion to the human host. Ultimately, it will be highly interesting to multiple pathoadaption events along with extensive horizontal correlate genotype to patient outcomes to identify genomic gene transfer events could constitute one of the causes for the groups that might be associated with a higher virulence. relatively large and non-structured accessory genome we detec- Taken together, our data show that strains from several diverse ted25. A striking observation achieved by long-read PacBio S. maltophilia complex lineages are associated with the hospital sequencing was the absence of plasmids in the completed gen- setting and human-associated infections, with lineage Sm6 strains omes that hence did not play a role in gene exchange and resis- potentially best adapted to colonise or infect humans. Strains of tance development in the selected S. maltophilia strains. this lineage are isolated worldwide, are found in potential human- It is well established that S. maltophilia is equipped with to-human transmission clusters and are predicted to be highly an armamentarium of antimicrobial resistance-conferring resistant to antibiotics and disinfectants. Accordingly, strict mechanisms5,19. In our strain collection, we found several compliance to infection prevention measures is important to families of antibiotic efflux pumps ubiquitously present among prevent and control nosocomial transmissions especially of S. strains of all 23 lineages, as well as other genes implicated in maltophilia lineage Sm6 strains, including the need to ensure that aminoglycoside or fluoroquinolone resistance. In some cases, the commonly used disinfectants are effective against S. mal- resistance-related genes were only present in some lineages, such tophilia complex strains expressing KatA. Future anti-infective as the β-lactamase gene blaL2 or the aminoglycoside acetyl- and treatment strategies may be based on our finding of a very low phosphotransferases genes aac and aph. Interestingly, those prevalence of trimethoprim-sulfomethoxazole-resistance genes in lineages harbouring mostly environmental strains tended to our collection, suggesting that this antibiotic drug remains the harbour less resistance and virulence genes than lineages that drug of choice for the treatment of S. maltophilia complex comprised at a majority human-associated strains. For instance, infections. the four lineages most distantly placed from the remaining S. maltophilia complex, Sgn1–Sgn4, were associated with the lack of Methods key virulence and resistance factors. In contrast, the human- Bacterial strains and DNA isolation. All Stenotrophomonas maltophilia complex associated lineage Sm6 was linked to the presence of a β- strains sequenced in this study were routinely collected in the participating hos- lactamase (BlaL2) and KatA, involved in resistance to disin- pitals and identified as S. maltophilia using MALDI-TOF MS. The strains were grown at 37 °C or 30 °C in either lysogeny broth (LB) or Brain Heart Infusion fectants, pointing towards adaptation to healthcare settings and media. RNA-free genomic DNA was isolated from 1-ml overnight cultures using survival on and in patients. While other human-associated the DNeasy Blood & Tissue Kit according to the manufacturer’s instructions lineages also harboured resistance and virulence genes at high (Qiagen, Hilden, Germany). To ensure correct identification as S. maltophilia, the proportions, this finding might explain why strains of lineage 16S rRNA sequence of S. maltophilia ATCC 13637 was blasted against all strains. The large majority (1278 strains, 98%) of our data set had 16S rRNA similarity Sm6 were dominant in our investigation, both in our total study values ≥ 99% (rounded to one decimal). Twenty-seven strains, mostly from the collection as well as in the subset of prospectively collected strains more distant clades Sgn1-4, had 16S rRNA blast results between 98.8% and 98.9%. as the majority of strains were isolated from human-associated Where no 16S rRNA sequence was found (one study using metagenome assembled sources. This notion is also supported by our finding that we did genomes43 as well as accession numbers GCA_000455625.1 and GCA_000455685.1) we left the isolates in our collection if the allele calls were not detect any d100 clusters, or circulating variants, in the pri- above the allele threshold of 2000 (Supplementary Fig. 1H). marily environmental-associated lineages. Yet, in light of the low number of strains belonging to these lineages in our data set as Whole-genome data collection and sequencing. We retrieved available S. mal- well as the lack of systematic sampling for environmental isolates, tophilia sequence read data sets and assembled genomes from NCBI nucleotide these results should be interpreted with caution. databases as of April 2018, excluding next-generation sequencing (NGS) data from Importantly, our study indicates the presence of potential non-Illumina platforms and data sets from studies that exclusively described transmission clusters in human-associated strains, suggesting mutants. For studies investigating serial strains from the same patient, we chose only representative strains, i.e. one sample per patient was chosen from Esposito potential direct or indirect human-to-human transmission17. et al.44 and one strain of the main lineages found by Chung et al.45. In case of Indeed, we identified a remarkable number of closely related studies providing both NGS data and assembled genomes, we included the NGS strains (270) that congregated in 62 clusters as indicated by a data in our analysis. maximum of ten mismatched alleles in the pairwise comparison. In addition, we sequenced the genomes of 1071 clinical and environmental strains. NGS libraries were constructed from genomic DNA using a modified While no d10 clusters were found in the more distantly placed Illumina Nextera protocol46 and the Illumina NextSeq 500 platform with 2 × lineages Sgn1-4, all other lineages comprised of such clusters with 151 bp runs (Illumina, San Diego, CA, USA). NGS data were assembled de novo similar clustering rates. A common source of infection is sup- using SPAdes (v3.7.1) included into the BioNumerics software (v7.6.3, Applied ported in those cases where detailed epidemiological information Maths NV). We excluded assemblies with an average coverage depth < 30 × (Supplementary Fig. 1A), deviating genome lengths (< 4 Mb and > 6 Mb) concerning hospital and day of isolation was available. Further (Supplementary Fig. 1B), number of contigs > 500 (Supplementary Fig. 1C), >2000 studies looking into potential transmission events are warranted non-ACTG bases (Supplementary Fig. 1D), an average quality < 30 as this would have major consequences on how infection (Supplementary Fig. 1E) and GC content (< 63% or > 68%) (Supplementary NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications 9 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 Fig. 1F). Fifty-five data sets where assembly completely failed were excluded from assigned for this particular locus. If for a specific locus the allele call is only further analysis. For the phylogenetic analysis, we further excluded strains available for one algorithm, this allele call will be included. If multiple allele possessing <2000 genes of the whole-genome MLST scheme constructed in this sequences were found for a consensus locus, only the lowest allele number is study (Supplementary Fig. 1H). The resulting data set contained 1305 samples (234 retained. Genes of which the sequence was not yet in the allele database were only from public databases) with a mean coverage depth of 130 × (SD = 58; median 122, assigned an allele number in case the sequence had valid start/stop codons, had no IQR 92–152), consisted of, on average, 74 contigs (mean, SD = 44; median 67, IQR ambiguous bases or internal stop codons, had at least 80% homology towards one 47–93) and encompassed a mean length of 4.7 million base pairs (SD = 0.19; of the reference allele sequences and had no more than 999 gaps in the pairwise median 4.76, IQR 4.64–4.87) (Supplementary Data 4). All assemblies were assessed sequence alignment towards the closest allele sequence from the same locus. The for completeness (range 81.03–100, mean 99.7, SD = 1.3) and contamination loci of the scheme were annotated using the blast2go tool51 relying on NCBI blast (range 0–10.8, mean 0.38, SD = 0.59) using CheckM47. version 2.4.0 + 52 and InterProscan 5 online53 (Supplementary Data 2), and the Next-generation sequencing data generated in the study are available from November 2018 GO54,55 and NCBI nr databases were used. All loci of the wgMLST public repositories under the study accession number “PRJEB32355” (accession scheme in FASTA format can be accessed using this link (https://figshare.com/ numbers for all data sets used are provided in Supplementary Data 3). articles/Smaltophilia_wgMLST_all-alleles_fasta_gz/10005047). Generation of full genomes by PacBio sequencing. We used PacBio long-read Whole-genome MLST scheme validation. To validate the scheme, publicly sequencing on an RSII instrument (Pacific Biosciences, Menlo Park, CA, USA) to available sequence read sets from different publications44,56 were analysed with the generate fully closed reference genome sequences of S. maltophilia complex strains wgMLST scheme in BioNumerics (v7.6.3). In addition, wgMLST analysis was sm454, sm-RA9, Sm53, ICU331, SKK55, U5, PEG-141, PEG-42, PEG-173, PEG-68, performed three times on the same sequence read set of two samples56. These PEG-305 and PEG-390, which together with available full genomes, represent the technical replicates had the same number of consensus allele calls and the allele majority of the diversity of our collection. SMRTbellTM template library was prepared numbers were identical. The allelic profiles of three between-run replicates according to the Procedure & Checklist—20 kb Template Preparation using the (sequencing data obtained from different fresh cultures) and three within-run BluePippinTM Size-Selection System (Pacific Biosciences, Menlo Park, CA, USA). replicates (sequencing data obtained from different libraries made using the same Briefly, for preparation of 15-kb libraries, 8 μg of genomic DNA from S. maltophilia DNA extract of one fresh culture) of S. maltophilia strain ATCC 1363756 were strains was sheared using g-tubesTM (Covaris, Woburn, MA, USA) according to the identical, except for one locus (STENO00008). The difference in allele calling for manufacturer’s instructions. DNA was end-repaired and ligated overnight to hairpin this locus, a gene coding for a ferric siderophore transport system/periplasmic adapters applying components from the DNA/Polymerase Binding Kit P6 (Pacific transport protein tonB, is likely due to sequencing and assembly difficulties of this Biosciences, Menlo Park, CA, USA). BluePippinTM Size-Selection to 7000 kb was GC-rich gene. Replicating the core genome SNP tree from Esposito and collea- performed as instructed (Sage Science, Beverly, MA, USA). Conditions for annealing gues44 based on wgMLST results yielded a highly similar tree topology clustering of sequencing primers and binding of polymerase to purified SMRTbellTM template samples from each patient with few exceptions. were assessed with the Calculator in RS Remote (Pacific Biosciences, Menlo Park, CA, USA). SMRT sequencing was carried out on the PacBio RSII (Pacific Biosciences, Phylogenetic analysis. We characterised the core loci present in 99% of the data Menlo Park, CA, USA) taking one 240-minutes movie for each SMRT cell. In total, set based on loci presence, i.e. that genes received a valid allele call, amounting to one SMRT cell for each of the strains was run. For each of the 12 genomes, 1275 loci. For phylogenetic analyses, a concatenated alignment of the 1275 core 59,220–106,322 PacBio reads with mean read lengths of 7678–13,952 base pairs (bp) genes from all strains was created, and an initial tree was built using RAxML-NG were assembled using the RS_HGAP_Assembly.3 protocol implemented in SMRT with a GTR + Gamma model,using the site-repeat optimisation, and 100 bootstrap Portal version 2.3.048. Subsequently, Illumina reads were mapped onto the assembled replicates57. This alignment and the tree were then fed to ClonalFrameML to detect sequence contigs using BWA (version 0.7.12)49 to improve the sequence quality to any regions of recombination58. These regions were then masked using maskrc-svg 99.9999% consensus accuracy. The assembled reads were subsequently disassembled (https://github.com/kwongj/maskrc-svg), and this masked alignment was then used for removal of low-quality bases. The contigs were then analysed for their synteny to to build a recombination-free phylogeny using the same approach as above in detect overlaps between its start of the anterior and the end of the posterior part to RAxML-NG. iTOL was employed for annotating the tree59. The core gene align- circularise the contigs. Finally, the dnaA open-reading frame was identified and shifted ment length was 1,070,730 variants, amounting to 1,397,302,650 characters for the to the start of the sequence. To evaluate structural variation, genomes were aligned entire data set. Across all isolates, 593,506,119 positions (42% of all variants) were using blastn. PlasmidFinder50 was used to screen the completed genomes for plasmids. masked for recombination. For phylogenetic and BAPS analysis, all invariant sites Genome sequences are available under bioproject number; the accession numbers can were removed to obtain the final alignment length of 296,491 variants. The be found in Supplementary Table 2. assemblies were annotated with prokka60, and the pan genome was calculated and visualised using roary61. Construction of a whole-genome MLST scheme. A whole-genome multilocus The genus wide tree showing a comparative phylogenetic analysis of the sequence typing (wgMLST) scheme was created by Applied Maths NV (bioMér- lineages with Stenotrophomonas species data (Supplementary Fig. 5) was calculated ieux) using 171 publically available S. maltophilia genome data sets. First, an initial using IQtree62 based on an alignment of the concatenated predicted protein set of loci was determined using the coding sequences (CDS) of the 171 genomes sequence of 23 genes15 (dnaG, rplA, rplB, rplC, rplD, rplE, rplF, rplK, rplL, rplM, (Supplementary Data S1). Within this set, loci that overlapped >75% or that yielded rplN, rplP, rplS, rpmA, rpoB, rpsB, rpsC, rpsE, rpsJ, rpsK, rpsM, rpsS, tsf) that were BLAST hits at the same position within one genome were omitted or merged until extracted from the assemblies using blastn. only mutually exclusive loci were retained while preserving maximal genome We detected phylogenetic lineages within the tree using a hierarchical Bayesian coverage. Mutually exclusive loci are defined as loci for which the reference alleles Analysis of Population Structure (hierBAPS) model as implemented in R (typically one or two unique DNA sequences per loci) only yield blast hits at a (rHierBAPs) with a maximum depth of 2 and maximum population number of threshold of 80% similarity to their own genomic location and not to reference 10063. FastANI64 was employed to calculate the pairwise average nucleotide alleles of another locus, such as paralogs or repetitive regions. In addition, loci that identity (ANI) as a similarity matrix between all the strains with the option ‘many- had a high ratio of invalid allele calls (e.g. because of the absence of a valid start/ to-many’. The similarity matrix was imported into R and used together with the stop codon [ATG, CTG, TTG, GTG], the presence of an internal stop codon [TAG, group assignment obtained from hierBAPS to compare the ANI values in strains TAA, TGA] or non-ACTG bases) and loci for which alleles were found containing within and between groups. ANI values were plotted as a heatmap of all strains as large tandem repeat areas were removed. Lastly, multi-copy loci, i.e. repeated loci well as a composite histogram of identity between and within groups. for which multiple allele calls were retrieved, were eliminated to achieve 90% of the genome data sets used for scheme validation had <10 repeated loci. The resulting Resistome and virulence analysis. Resistome and virulome were characterised scheme contained 17,603 loci (including the seven loci from the previously pub- with abricate version 0.8.736 screened against the NCBI Bacterial Antimicrobial lished MLST scheme24, see Supplementary Table 1) (Supplementary Figs. 2 and 3) Resistance Reference Gene Database (NCBI BARRGD, PRJNA313047) and the and can be accessed through a plugin in the BioNumericsTM Software (Applied Virulence Factors of Pathogenic Bacteria Database (VFDB)37. All genes below 80% Maths NV, bioMérieux). On average, 4174 loci (range 3024–4536) were identified coverage breadth were excluded. In addition, literature was reviewed to identify per genome of our study collection. additional genes associated with antibiotic resistance, and virulence in S. mal- To determine the allele number(s) corresponding to a unique allele sequence for tophilia and the corresponding loci were extracted from the wgMLST scheme each locus present in the genome of a strain, two different algorithms were (emrA, emrE, sugE, norM, clpA, clpP, stmPr1, htpX, tetACG, smoR, smeU2, sul1 and employed: the assembly-free (AF) allele calling uses a k-mer approach (k-mers size sul2). Multiple correspondence analysis (MCA) was performed on nine genes (aac, of 35 with minimum coverage of 3) starting from the raw sequence reads while the aph, blaL1, blaL2, katA, pilU, qacE, smoR, sul1) that were present or absent in at assembly-based (AB) allele calling performs a blastn search against assembled least 10 isolates as these were most likely to explain data set variance. The analysis genomes with the reference alleles of each loci as query sequences. The word size was conducted using the factoextra and FactoMineR R packages65. for the gapped blast search was set at 11, and only hits with a minimum homology of 80% were retained. After each round of allele identification, all the available data from the two algorithms (AF and AB) were combined into a single set of allele Statistical analysis and data management. All statistical analyses and data assignments, called consensus calls. If both algorithms returned one or multiple management were performed in R version 3.4.366 using mainly packages included allele calls for a given loci, the consensus is defined as the allele(s) that both in the tidyverse67, reshape268 and rcompanion (https://cran.r-project.org/web/ analyses have in common. If there is no overlap, there will be no allele number packages/rcompanion/index.html). The map was created with the tmap package69. 10 NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 ARTICLE WgMLST data were analysed in BioNumerics v7.6.3 using the WGS and MLST 12. Chang, Y.-T., Lin, C.-Y., Chen, Y.-H. & Hsueh, P.-R. Update on infections plugins. For correlation testing between variables, Spearman’s rank test was caused by Stenotrophomonas maltophilia with particular attention to employed. The test of equal or given proportions or Fisher’s exact test for sample resistance mechanisms and therapeutic options. Front. Microbiol. 6, 893 sizes smaller than 5 was used to test for proportions. Association analysis was (2015). performed at the unitig level using the linear mixed model and whole-genome 13. Hu, L.-F. et al. Stenotrophomonas maltophilia resistance to trimethoprim/ elastic net models as implemented in pyseer70. Unitigs were counted with unitig- sulfamethoxazole mediated by acquisition of sul and dfrA genes in a plasmid- counter71, and a genetic similarity matrix was calculated from the core phylogeny. mediated class 1 integron. Int. J. Antimicrob. Agents 37, 230–234 (2011). All isolates with unknown or anthropogenic isolation source were excluded from 14. Vinuesa, P., Ochoa-Sánchez, L. E. & Contreras-Moreira, B. the analysis. Any unitig with a minor allele frequency of <1% was not tested. GET_PHYLOMARKERS, a software package to select optimal orthologous Control for population structure was verified with QQ plots using the packae clusters for phylogenomics and inferring pan-genome phylogenies, used for a qqman72 (Supplementary Fig. 8). We additionally ran a whole-genome elastic net critical geno-taxonomic revision of the genus Stenotrophomonas. Front. model where population structure was ensured by including sequence reweighting Microbiol. 9, 771 (2018). and supplying the BAPS groups as clonal complex designations. 15. Patil, P. P., Kumar, S., Midha, S., Gautam, V. & Patil, P. B. Taxonogenomics reveal multiple novel genomospecies associated with clinical isolates of Reporting summary. Further information on research design is available in Stenotrophomonas maltophilia. Microb. Genome 4, (2018). the Nature Research Reporting Summary linked to this article. 16. Ochoa-Sánchez, L. E. & Vinuesa, P. Evolutionary genetic analysis uncovers multiple species with distinct habitat preferences and antibiotic resistance phenotypes in the Stenotrophomonas maltophilia complex. Front. Microbiol. 8, Data availability 1548 (2017). The sequencing data generated and analysed during the current study are available in the 17. Steinmann, J. et al. Analysis of phylogenetic variation of Stenotrophomonas European Nucleotide Archive repository under the accession numbers PRJEB32355, maltophilia reveals human-specific branches. Front. Microbiol. 9, 806 (2018). PRJEB32585 and PRJNA543082 (Supplementary Table 2). The Bacterial Antimicrobial 18. Svensson-Stadler, L. A., Mihaylova, S. A. & Moore, E. R. B. Stenotrophomonas Resistance Reference Gene Database is available under accession PRJNA313047. All data interspecies differentiation and identification by gyrB sequence analysis. FEMS generated or analysed during this study are included in this published article and its Microbiol. Lett. 327, 15–24 (2012). supplementary information. The source data underlying Figs. 1a–c, 3, 4a, Supplementary 19. Crossman, L. C. et al. The complete genome, comparative and functional Figs. 3b, 6 are provided as a Source Data File, along with supplemental Supplementary analysis of Stenotrophomonas maltophilia reveals an organism heavily shielded Data 3, 5. Supplementary Data 3 contains the sample metadata used to create Figs. 1d, 2, by drug resistance determinants. Genome Biol. 9, R74 (2008). 5, and the tree annotation in Fig. 4a. Supplementary Data 5 comprises the allelic typing 20. Maiden, M. C. J. et al. MLST revisited: the gene-by-gene approach to bacterial information for all samples and was used to create Figs. 4b, 5, and S7. The Supplementary genomics. Nat. Rev. Microbiol. 11, 728–736 (2013). Figs. 1 and 2 are based on data provided in Supplementary Data 4 and Supplementary 21. Ruppitsch, W. et al. Defining and evaluating a core genome multilocus Data 2, respectively. The wgMLST scheme is available from https://figshare.com/articles/ sequence typing scheme for whole-genome sequence-based typing of Listeria Smaltophilia_wgMLST_all-alleles_fasta_gz/10005047. monocytogenes. J. Clin. Microbiol. 53, 2869–2876 (2015). 22. de Been, M. et al. Core genome multilocus sequence typing scheme for high- Code availability resolution typing of Enterococcus faecium. J. Clin. Microbiol. 53, 3788–3797 Custom python scripts used in the phylogenetic analyses are available from https:// (2015). github.com/conmeehan/pathophy. Custom R scripts used to calculate the statistical tests 23. Moura, A. et al. Whole genome-based population biology and epidemiological and perform the MCA analysis in Fig. 3b–d and Supplementary Fig. 5 are available from surveillance of Listeria monocytogenes. Nat. Microbiol. 2, 16185 (2016). https://github.com/ngs-fzb/SMaltophilia_phylo and can be replicated using the source 24. Kaiser, S., Biehler, K. & Jonas, D. A Stenotrophomonas maltophilia multilocus data files provided. sequence typing scheme for inferring population structure. J. Bacteriol. 191, 2934–2943 (2009). 25. Youenou, B. et al. Comparative genomics of environmental and clinical Received: 13 August 2019; Accepted: 15 February 2020; Stenotrophomonas maltophilia strains with different antibiotic resistance profiles. Genome Biol. Evol. 7, 2484–2505 (2015). 26. Weber, M., Schünemann, W., Fuß, J., Kämpfer, P. & Lipski, A. Stenotrophomonas lactitubi sp. nov. and Stenotrophomonas indicatrix sp. nov., isolated from surfaces with food contact. Int. J. Syst. Evol. Microbiol. 68, 1830–1838 (2018). References 27. Saber, M. M. & Shapiro, B. J. Benchmarking bacterial genome-wide 1. Bryant, J. M. et al. Emergence and spread of a human-transmissible association study methods using simulated genomes and phenotypes. Microb. multidrug-resistant nontuberculous mycobacterium. Science 354, 751–757 Genom. https://doi.org/10.1099/mgen.0.000337 (2020). (2016). 28. Brooke, J. S. New strategies against Stenotrophomonas maltophilia: a serious 2. van Ingen, J. et al. Global outbreak of severe Mycobacterium chimaera disease worldwide intrinsically drug-resistant opportunistic pathogen. Expert Rev. after cardiac surgery: a molecular epidemiological study. Lancet Infect. Dis. 17, Anti. Infect. Ther. 12, 1–4 (2014). 1033–1041 (2017). 29. Blanco, P. et al. Bacterial multidrug efflux pumps: much more than antibiotic 3. Coll, F. et al. Longitudinal genomic surveillance of MRSA in the UK reveals resistance determinants. Microorganisms 4, E14 (2016). transmission patterns in hospitals and the community. Sci. Transl. Med. 9, pii: 30. Huang, H.-H. et al. ClpA and HtpX proteases are involved in intrinsic eaak9745 (2017). aminoglycoside resistance of Stenotrophomonas maltophilia and are potential 4. World Health Organization. Public health importance of antimicrobial aminoglycoside adjuvant targets. Antimicrob. Agents Chemother. 62, resistance. http://www.who.int/drugresistance/AMR_Importance/en/ (2013). e00554–18 (2018). Accessed February 26th, 2020. 31. Sánchez, M. B. Antibiotic resistance in the opportunistic pathogen 5. Brooke, J. S. Stenotrophomonas maltophilia: an emerging global opportunistic Stenotrophomonas maltophilia. Front. Microbiol. 6, 658 (2015). pathogen. Clin. Microbiol. Rev. 25, 2–41 (2012). 32. Martínez, P. et al. Stenotrophomonas maltophilia responds to exogenous AHL 6. Mukherjee, P. & Roy, P. Genomic potential of Stenotrophomonas maltophilia signals through the LuxR solo SmoR (Smlt1839). Front. Cell. Infect. Microbiol. in bioremediation with an assessment of its multifaceted role in our 5, 41 (2015). environment. Front. Microbiol. 7, 967 (2016). 33. Comolli, J. C. et al. Pseudomonas aeruginosa gene products PilT and PilU are 7. Falagas, M. E. et al. Attributable mortality of Stenotrophomonas maltophilia required for cytotoxicity in vitro and virulence in a mouse model of acute infections: a systematic review of the literature. Future Microbiol. 4, pneumonia. Infect. Immun. 67, 3625–3630 (1999). 1103–1109 (2009). 34. Windhorst, S. et al. The major extracellular protease of the nosocomial 8. Zemanick, E. T. et al. Airway microbiota across age and disease spectrum in pathogen Stenotrophomonas maltophilia: characterization of the protein and cystic fibrosis. Eur. Respir. J. 50, 1700832 (2017). pii. molecular cloning of the gene. J. Biol. Chem. 277, 11042–11049 (2002). 9. Denton, M. & Kerr, K. G. Microbiological and clinical aspects of infection 35. Charoenlap, N. et al. Inactivation of ahpC renders Stenotrophomonas associated with Stenotrophomonas maltophilia. Clin. Microbiol. Rev. 11, maltophilia resistant to the disinfectant hydrogen peroxide. Antonie Van. 57–80 (1998). Leeuwenhoek https://doi.org/10.1007/s10482-018-1203-9 (2018). 10. Henning, S. & Teare, L. Stenotrophomonas maltophilia infection: an unusual 36. Seemann, T. abricate: Mass Screening of Contigs for Antimicrobial and complication of total elbow arthroplasty. JMM Case Rep. 1, e000661 (2014). Virulence Genes. (Github). https://github.com/tseemann/abricate. 11. Falagas, M. E., Kastoris, A. C., Vouloumanou, E. K. & Dimopoulos, G. 37. Chen, L., Zheng, D., Liu, B., Yang, J. & Jin, Q. VFDB 2016: hierarchical and Community-acquired Stenotrophomonas maltophilia infections: a systematic refined dataset for big data analysis—10 years on. Nucleic Acids Res. 44, review. Eur. J. Clin. Microbiol. Infect. Dis. 28, 719–730 (2009). D694–D697 (2016). NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications 11 ARTICLE NATURE COMMUNICATIONS | https://doi.org/10.1038/s41467-020-15123-0 38. Patil, P. P., Midha, S., Kumar, S. & Patil, P. B. Genome sequence of type 70. Lees, J. A., Galardini, M., Bentley, S. D., Weiser, J. N. & Corander, J. pyseer: a strains of genus Stenotrophomonas. Front. Microbiol. 7, 309 (2016). comprehensive tool for microbial pangenome-wide association studies. 39. Trifonova, A. & Strateva, T. Stenotrophomonas maltophilia—a low-grade Bioinformatics 34, 4310–4312 (2018). pathogen with numerous virulence factors. Infect. Dis. 51, 168–178 (2019). 71. Jaillard, M. et al. A fast and agnostic method for bacterial genome-wide 40. Kim, E. J. et al. Risk factors for mortality in patients with Stenotrophomonas association studies: bridging the gap between k-mers and genetic events. PLoS maltophilia bacteremia and clinical impact of quinolone-resistant strains. Genet. 14, e1007758 (2018). BMC Infect. Dis. 19, 754 (2019). 72. Turner, S. D. qqman: an R package for visualizing GWAS results using Q-Q 41. David, S. et al. Multiple major disease-associated clones of Legionella and manhattan plots. JOSS 3, 731 (2018). pneumophila have emerged recently and independently. Genome Res. 26, 1555–1564 (2016). 42. Gomez-Valero, L. et al. More than 18,000 effectors in the Legionella genus Acknowledgements genome provide multiple, independent combinations for replication in human We thank V. Mohr, F. Boysen and T. Ubben for technical assistance during next- cells. Proc. Natl Acad. Sci. USA https://doi.org/10.1073/pnas.1808016116 (2019). generation sequencings. Parts of the work have been funded by grants from German 43. Parks, D. H. et al. Recovery of nearly 8,000 metagenome-assembled genomes Center for Infection Research, Federal Ministry of Education and Research, Germany, substantially expands the tree of life. Nat. Microbiol. 2, 1533–1542 (2017). from Deutsche Forschungsgemeinschaft (DFG, German Research Foundation) under 44. Esposito, A. et al. Evolution of Stenotrophomonas maltophilia in cystic fibrosis Germany’s Excellence Strategy–EXC 22167-390884018, the Cluster Precision Medi- lung over chronic infection: a genomic and phenotypic population study. cine in Chronic Inflammation, and grants from the Leibniz Science Campus Evo- Front. Microbiol. 8, 1590 (2017). LUNG; I.A. and W.R.S. received funds from the University of Hamburg; C.F.L. and 45. Chung, H. et al. Global and local selection acting on the pathogen J.W.A.R. were funded through the European Commission Horizon 2020 Framework Stenotrophomonas maltophilia in the human lung. Nat. Commun. 8, 14078 Marie Skłodowska-Curie Actions (Grant agreement number: 713660 - PRONKJE- (2017). WAIL - H2020-MSCA-COFUND-2015); C.J.M. was funded by the Department of 46. Baym, M. et al. Inexpensive multiplexed library preparation for megabase- Economy, Science and Innovation of the Flemish Government and the Belgian Sci- sized genomes. PLoS ONE 10, e0128036 (2015). ence Policy (Belspo); K.Z. was funded by the National Natural Science Foundation of 47. Parks, D. H., Imelfort, M., Skennerton, C. T., Hugenholtz, P. & Tyson, G. W. China (grant number 81702045), and 13th Five-Year National Major Science and CheckM: assessing the quality of microbial genomes recovered from isolates, Technology Projects of China (grant number 2018ZX10712001); O.C.S., D.Y., I.G. single cells, and metagenomes. Genome Res. 25, 1043–1055 (2015). and X.D. were funded by the Spanish Ministry for Science, Innovation and Uni- 48. Chin, C.-S. et al. Nonhybrid, finished microbial genome assemblies from long- versities (grant reference BIO2015-66674-R); N.A.R. was funded by the SingHealth read SMRT sequencing data. Nat. Methods 10, 563–569 (2013). DUKE-NUS Pathology Academic Clinical Programme Clinical Innovation Support 49. Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows- Grant (09/FY2017/P1/06-A20); U.N. was funded by the EU Horizon 2020 pro- Wheeler transform. Bioinformatics 25, 1754–1760 (2009). gramme, grant agreement number 643476. 50. Carattoli, A. et al. In silico detection and typing of plasmids using PlasmidFinder and plasmid multilocus sequence typing. Antimicrob. Agents Author contributions Chemother. 58, 3895–3903 (2014). M.I.G., J.R., J.S., S.N. and T.A.K. conceived the study. M.I.G., I.B., M.G., O.C.S., J.S. and 51. Götz, S. et al. High-throughput functional annotation and data mining with T.A.K. curated the data. M.I.G., C.J.M., I.B., M.D., A.G., M.S., I.C.S., C.U., C.F.L., D.Y., the Blast2GO suite. Nucleic Acids Res. 36, 3420–3435 (2008). M.R.F., O.C.S. and T.A.K. performed the formal analysis. J.R., J.S. and S.N. acquired 52. Camacho, C. et al. BLAST+: architecture and applications. BMC Bioinforma. funding. S.K., N.A.R., M.K., T.S.v.d.W., W.R.S., K.D.B., K.Z., T.S., I.A., W.S., K.Z., U.E.S., 10, 421 (2009). J.W.A.R., U.N., D.Y., I.G., X.D., J.R., U.M., M.R.F. and J.S. provided resources for this 53. Jones, P. et al. InterProScan 5: genome-scale protein function classification. study. M.I.G., I.B., A.G. and C.F.L. visualised the results. M.I.G., J.S., S.N. and T.A.K. Bioinformatics 30, 1236–1240 (2014). wrote the initial draft. All authors critically reviewed and modified the paper. 54. Ashburner, M. et al. Gene ontology: tool for the unification of biology. The Gene Ontology Consortium. Nat. Genet. 25, 25–29 (2000). 55. The Gene Ontology Consortium. The Gene Ontology Resource: 20 years and Competing interests still GOing strong. Nucleic Acids Res. 47, D330–D338 (2019). The authors declare no competing interests. 56. Kozyreva, V. K. et al. Validation and implementation of clinical laboratory improvements Act-Compliant Whole-Genome Sequencing in the public health microbiology laboratory. J. Clin. Microbiol. 55, 2502–2520 (2017). Additional information Supplementary information is available for this paper at https://doi.org/10.1038/s41467- 57. Stamatakis, A. RAxML-VI-HPC: maximum likelihood-based phylogenetic 020-15123-0. analyses with thousands of taxa and mixed models. Bioinformatics 22, 2688–2690 (2006). Correspondence and requests for materials should be addressed to S.N. 58. Didelot, X. & Wilson, D. J. ClonalFrameML: efficient inference of recombination in whole bacterial genomes. PLoS Comput. Biol. 11, e1004041 (2015). Peer review information Nature Communications thanks Kelly Wyres and the other, 59. Letunic, I. & Bork, P. Interactive tree of life (iTOL) v3: an online tool for the anonymous, reviewer(s) for their contribution to the peer review of this work. Peer display and annotation of phylogenetic and other trees. Nucleic Acids Res. 44, reviewer reports are available. W242–W245 (2016). 60. Seemann, T. Prokka: rapid prokaryotic genome annotation. Bioinformatics 30, Reprints and permission information is available at http://www.nature.com/reprints 2068–2069 (2014). 61. Page, A. J. et al. Roary: rapid large-scale prokaryote pan genome analysis. Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in Bioinformatics 31, 3691–3693 (2015). published maps and institutional affiliations. 62. Nguyen, L.-T., Schmidt, H. A., von Haeseler, A. & Minh, B. Q. IQ-TREE: a fast and effective stochastic algorithm for estimating maximum-likelihood phylogenies. Mol. Biol. Evol. 32, 268–274 (2015). 63. Cheng, L., Connor, T. R., Sirén, J., Aanensen, D. M. & Corander, J. Open Access This article is licensed under a Creative Commons Hierarchical and spatially explicit clustering of DNA sequences with BAPS Attribution 4.0 International License, which permits use, sharing, software. Mol. Biol. Evol. 30, 1224–1228 (2013). adaptation, distribution and reproduction in any medium or format, as long as you give 64. Jain, C., Rodriguez-R, L. M., Phillippy, A. M., Konstantinidis, K. T. & Aluru, S. appropriate credit to the original author(s) and the source, provide a link to the Creative High throughput ANI analysis of 90K prokaryotic genomes reveals clear Commons license, and indicate if changes were made. The images or other third party species boundaries. Nat. Commun. 9, 5114 (2018). material in this article are included in the article’s Creative Commons license, unless 65. Lê, S., Josse, J. & Husson, F. FactoMineR: an R package for multivariate indicated otherwise in a credit line to the material. If material is not included in the analysis. J. Stat. Softw., Artic. 25, 1–18 (2008). article’s Creative Commons license and your intended use is not permitted by statutory 66. R Core Team. R: A Language and Environment for Statistical Computing (R regulation or exceeds the permitted use, you will need to obtain permission directly from Foundation for Statistical Computing, 2013). the copyright holder. To view a copy of this license, visit http://creativecommons.org/ 67. Wickham, H. et al. Welcome to the tidyverse. JOSS 4, 1686 (2019). licenses/by/4.0/. 68. Wickham, H. Reshaping data with the reshape package. J. Stat. Softw. 21, 1–20 (2007). 69. Tennekes, M. tmap: thematic maps in R. J. Stat. Softw., Artic. 84, 1–39 (2018). © The Author(s) 2020 12 NATURE COMMUNICATIONS | (2020)11:2044 | https://doi.org/10.1038/s41467-020-15123-0 | www.nature.com/naturecommunications