Global population structure of the Serratia marcescens complex and identification of hospital-adapted lineages in the complex

Abstract

Serratia marcescens is an important nosocomial pathogen causing various opportunistic infections, such as urinary tract infections, bacteremia and sometimes even hospital outbreaks. The recent emergence and spread of multidrug-resistant (MDR) strains further pose serious threats to global public health. This bacterium is also ubiquitously found in natural environments, but the genomic differences between clinical and environmental isolates are not clear, including those between S. marcescens and its close relatives. In this study, we performed a large-scale genome analysis of S. marcescens and closely related species (referred to as the 'S. marcescens complex'), including more than 200 clinical and environmental strains newly sequenced here. Our analysis revealed their phylogenetic relationships and complex global population structure, comprising 14 clades, which were defined based on whole-genome average nucleotide identity. Clades 10, 11, 12 and 13 corresponded to S. nematodiphila, S. marcescens sensu stricto, S. ureilytica and S. surfactantfaciens, respectively. Several clades exhibited distinct genome sizes and GC contents and a negative correlation of these genomic parameters was observed in each clade, which was associated with the acquisition of mobile genetic elements (MGEs), but different types of MGEs, plasmids or prophages (and other integrative elements), were found to contribute to the generation of these genomic variations. Importantly, clades 1 and 2 mostly comprised clinical or hospital environment isolates and accumulated a wide range of antimicrobial resistance genes, including various extended-spectrum β-lactamase and carbapenemase genes, and fluoroquinolone target site mutations, leading to a high proportion of MDR strains. This finding suggests that clades 1 and 2 represent hospital-adapted lineages in the S. marcescens complex although their potential virulence is currently unknown. These data provide an important genomic basis for reconsidering the classification of this group of bacteria and reveal novel insights into their evolution, biology and differential importance in clinical settings.

Keywords: Serratia marcescens complex; antimicrobial resistance; genomics; hospital-adapted lineage; population structure.

Fig. 1.

Phylogenetic relationships between Serratia species. (a) An ML phylogenetic tree based on the core gene SNPs (296 627 sites in 801 genes) of the type strains of 17 Serratia species/subspecies and two completely sequenced S. marcescens strains (SM39 and Db11), with the Y. enterocolitica type strain as an outgroup. The core genes that are present in all strains were identified using Roary with a 60 % blastp identity cutoff. (b) ANI matrix among the nine strains of S. marcescens (Sma) and its close relatives. All-to-all ANI values were calculated using PYANI in ANIb mode, and cells are coloured according to ANI values. Strains other than Db11 and SM39 are type strains of each species/subspecies.

Fig. 2.

Regions and sources of the isolation of the Sma complex strains. (a) Distribution of the regions where strains were isolated. Regions were classified into seven geographic regions: the UK, Japan, North America, South America, Europe other than the UK (Europe w/o UK), Asia/Oceania other than Japan (Asia/Oceania w/o Japan), and Africa. (b) Distribution of isolation sources. Subcategories that contained <10 strains were grouped as ‘others’. The proportion of strains that were sequenced in this study is indicated in the most inner circles in both panels.

Fig. 3.

ANI-based clustering of the Sma complex strains. Complete-linkage unsupervised hierarchical clustering of the Sma complex strains was performed based on their all-to-all ANI analysis. The type strains of S. marcescens (Sma^T), S. nematodiphila (Sne^T), and two Sma strains (Db11 and SM39) are indicated by their names. Nodes whose minimum pairwise ANI scores were <97 % are indicated by coloured symbols according to minimum within-node ANI scores. Fourteen clades defined according to the minimum within-node ANI threshold of 97 % are indicated under the tree.

Fig. 4.

Phylogenetic relationships among the Sma complex strains. The phylogenetic relationship between the Sma complex strains analysed in this study (n=775) is shown. An ML tree was constructed based on the 182 520 SNP sites identified in 2488 core genes. The root was defined by a similar phylogenetic analysis using the type strain of S. ficalia as an outgroup. Blue diamonds indicate the strains whose complete genome sequences are available. The type strains of Sma and Sne and two Sma strains (Db11 and SM39) are indicated by their names and red dotted lines. For each strain, the isolation region, isolation source, clade, presence/absence of the pig gene cluster, total genome size and GC content are shown. Branches where the deletion of the pig gene cluster occurred are indicated by arrows.

Fig. 5.

Relationship between the core-gene-based phylogeny and the pangenome-based clustering of the Sma complex strains. The upper tree is the core-gene-based ML tree shown in Fig. 4, and the lower dendrogram was constructed based on the cluster analysis of the pangenome (n=34 342) identified by Roary. Ward’s clustering method with Euclidean distance was used for the cluster analysis. The same strains present in the two trees are connected by coloured lines according to their clades.

Fig. 6.

Interclade differences in genome size and GC content. Genome size ranges (a) and GC content ranges (b) of SNP₁₀ clusters belonging to each clade. For multimember SNP₁₀ clusters, one strain was randomly selected from each cluster. Open diamonds in the boxplot indicate the mean values of each clade. Clades that show significant differences (P<0.01) compared to the clades marked by diamonds are indicated by short vertical lines. Lines are coloured according to the clade indicated by diamonds. Only the combinations that show significant interclade differences are shown. The numbers of total strains and SNP₁₀ clusters in each clade are shown at the bottom of (d). (c, d) The numbers of plasmid replicons and integrase genes found in the SNP₁₀ clusters belonging to each clade. Diamonds indicate the averages in each clade. (e) Correlations between genome size and GC content (left) and between genome size and the numbers of plasmid replicons (middle) or integrase genes (right) in nine clades that contained >10 members. Linear regression lines and R-squared results are shown. In all panels, colours for each clade are the same as those in Figs. 3–5.

Fig. 7.

Proportions of the SNP₁₀ clusters isolated from clinical/hospital environmental sources in each clade and numbers of AMR genes in each SNP₁₀ cluster belonging to each clade. (a) The proportions of the SNP₁₀ clusters isolated from clinical/hospital environment samples in each clade. The dotted line indicates the proportion of clinical/hospital environment clusters in the entire strain set (73.8%). Note that the isolation sources of the strains belonging to the same SNP₁₀ clusters were the same for all multimember clusters. *; P<0.05, **; P<0.01. (b) The numbers of AMR genes in each SNP₁₀ cluster belonging to each clade. The numbers of AMR genes in multimember SNP₁₀ clusters were calculated as the average numbers in each SNP₁₀ cluster. The numbers of MDR strain-containing SNP₁₀ clusters, total SNP₁₀ clusters, MDR strains and total strains belonging to each clade are indicated at the bottom. In both panels, colours for each clade are the same as those in Figs. 3–6.

Fig. 8.

Distributions of AMR genes and FQR-conferring mutations in the Sma complex strains. The same ML tree and strain information on the clade and isolation source presented in Fig. 4 are shown here. Different mutations in the parC and parE genes are indicated by different colours. Carbapenemase and ESBL genes are marked by * and **, respectively.