Enterobacter hormaechei subsp. hoffmannii subsp. nov., Enterobacter hormaechei subsp. xiangfangensis comb. nov., Enterobacter roggenkampii sp. nov., and Enterobacter muelleri is a later heterotypic synonym of Enterobacter asburiae based on computational analysis of sequenced Enterobacter genomes.

Introduction

The name Enterobacter hormaechei was created for a taxon at the rank of species that had previously been called Enteric Group 75. O’Hara et al¹. defined the type strain to be ATCC 49162^T from the 23 strains they studied. Twelve of the strains were shown to be closely related via DNA-DNA hybridization (DDH) and less closely related to other Enterobacter species. Numerous biochemical assays were performed on the 23 strains to characterize and differentiate the new species.

Hoffmann and Roggenkamp² investigated the genetic structure of the E. cloacae complex (the set of species included in this complex has varied over time) by a combination of sequencing of the three housekeeping genes hsp60, rpoB, and hemB; and PCR-restriction fragment length polymorphism (PCR-RFLP) analysis of ampC. They defined 12 genetic clusters (I-XII) based most exhaustively on the hsp60 sequencing. Three of the clusters (cluster III, 58 strains; cluster VI, 28 strains; and cluster VIII, 59 strains) accounted for 70% of the 206 strains studied. The authors noted that “Only 3% of our study strains clustered with the type strain of E. cloacae.” (cluster XI), “We found that 3% of our study strains clustered around the E. hormaechei type strain.” (cluster VII), and “Our clusters VI and VIII were closely related to E. hormaechei cluster VII. DDH studies are needed to verify whether these clusters form a common DNA relatedness group allowing emending and broadening of the species description of E. hormaechei.”.

Hoffmann et al³. followed up with a characterization of clusters VI, VII, and VIII asserting based on DDH that these clusters were subspecies of the same species. Since cluster VII contained the type strain for E. hormaechei Hoffmann et al. named cluster VII E. hormaechei subsp. hormaechei, cluster VI E. hormaechei subsp. oharae, and cluster VIII E. hormaechei subsp. steigerwaltii. Forty-eight strains were characterized using 129 biochemical tests showing that there were phenotypic differences between the subspecies. Unfortunately the authors did not decide to include the other predominant cluster (III) in their analysis, nor did they validly publish these subspecies names. This was rectified recently in Validation List no. 172⁴.

Gu et al⁵. defined E. xiangfangensis using a phylogenetic tree based upon concatenated partial rpoB, atpD, gyrB and infB gene sequences from a novel isolate and existing type strains where E. xiangfangensis grouped closest to E. hormaechei. Biochemical assays were performed and E. xiangfangensis strains were differentiable from the E. hormaechei type strain.

During analysis of the E. cloacae complex and E.(now Klebsiella⁶) aerogenes strains looking at antimicrobial resistance patterns⁷, many of the Hoffmann et al. clusters were rediscovered using whole genome comparisons such as SNP analysis and average nucleotide identity (ANI). The clusters were identifiable by the hsp60 sequences deposited by the Hoffmann group. The three subspecies of E. hormaechei defined by Hoffmann et al. fell within the expected ANI range for bacterial species, being greater than 95% ANI between subspecies and greater than 98% ANI within a subspecies. Unexpectedly Hoffmann cluster III also met the ANI criteria to be an E. hormaechei subspecies. Further, genomes named E. xiangfangensis in GenBank fell within the E. hormaechei subsp. steigerwaltii cluster rather than a separate cluster. Moreover, most of the genomes in these clusters were mistakenly identified as E. cloacae when they were submitted to GenBank. To resolve the naming inconsistencies of these genomes the type strains for E. hormaechei subsp. steigerwaltii, E. hormaechei subsp. oharae, E. xiangfangensis, Hoffmann cluster III, and Hoffmann cluster IV were sequenced.

Background

Tools for bacterial species assignment have changed over time^8,9. Initially, morphology as viewed through a microscope and later aided by staining such as Gram staining¹⁰ to distinguish cell wall differences was used. Biochemical assays and other methods to determine phenotype followed. Use of the genome started with DNA-DNA hybridization (DDH) where a 70% threshold for species followed later by a 79% threshold for subspecies were proposed. Widespread use of marker genes in particular the 16S rRNA gene made assays easier. A threshold of less than 97% identity for the 16S rRNA gene was used to determine a new species but values above 97% could not guarantee that isolates were the same species. The sequence of other less conserved marker genes such as hsp60 has also been used to differentiate species. More recently multiple marker genes are sequenced and a combined alignment is used. With the advent of inexpensive genome sequencing, computing average nucleotide identity (ANI), which correlates very closely with DDH, has largely supplanted other methods. Studies have shown that an ANI threshold between 94-96.5% correlates well with existing species definitions and 97-98% for subspecies^11–19. DDH has been shown to not only correlate with ANI but also with how many of the genes or what fraction of the genomes are shared in common so some ANI based tools take this measurement into account as well^17–19. Most definitions of new species involve sequencing the genome and taking ANI and shared gene content into account in some fashion but many species definitions predate genome sequencing and some type strains have not been sequenced. There is no generally accepted method for reconciling older species definitions with genome comparisons but usually ANI and shared gene content form a basis for the analysis.

As Hoffmann^2,3 and others^20–26 discovered the predominant species in clinical Enterobacter isolates is E. hormaechei. Unfortunately many articles, clinicians, and GenBank submissions misname these strains as E. cloacae perhaps as a short hand for the E. cloacae complex and possibly due to the E. hormaechei subspecies not being validly published until recently. Another issue was the lack of sequenced type strains or named species/subspecies for some clades. The definition of what species/subspecies make up the E. cloacae complex has been in flux^2,27,28 and even what species are in the genus Enterobacter^29–31.

The E. cloacae complex was shown to have 18 clades (A-R)⁷, 12 of which corresponded to 11 of the 12 clusters defined previously by Hoffmann². Hoffmann cluster X is E. nimipressuralis which has been reclassified as Lelliottia nimipressuralis²⁹. Table 1 incorporates more recently sequenced genomes and published papers adding four clades (S-V) and incorporating the latest literature. For example, clade R (Hoffmann cluster IX) was recently defined to be E. bugandensis³¹.

Results

All RefSeq genomes labelled as being in the genus Enterobacter were downloaded from NCBI RefSeq resulting in 1,249 genomes. A fast approximate ANI tool, called MASH³², was used to generate a pairwise ANI based distance matrix and average linkage hierarchical clustering was used to generate the tree shown in Figure 1. 1,216 genomes were assigned to 22 clades (A-V Table 1) in the E. cloacae complex (Supplemental Table 1) while 30 genomes were deemed to be outliers and not in the Enterobacter genus (best MASH matches in Supplemental Table 2) as well as 2 E. lignolyticus genomes and 1 E. timonensis genome deemed to be outside of the E. cloacae complex. Two species of Enterobacter: E. siamensis and E. tabaci do not have sequenced genomes and their type strains’ 16S rRNA sequences while having full length matches at 98% and 99% respectively to some E. cloacae complex genomes did not have definitive matches to any particular clade. The type strains for E. asburiae and E. muelleri fall within the same clade (J – Hoffmann cluster I). All 78 genomes in this clade are above the 95% ANI species cut-off (Table 2) but using a 98% ANI subspecies cut-off produces 8 subclades of sizes 1, 1, 2, 2, 2 (E. muelleri), 3 (E. asburiae), 24, and 43. Thus E. muelleri³³ is a later heterotypic synonym of E. asburiae³⁴ which should take precedence. Whether the 8 subclades of E. asburiae should be treated as subspecies is beyond the scope of this paper.

Figure 1. Average nucleotide identity (ANI) based tree for 1249 NCBI RefSeq Enterobacter labelled genomes.

Five clades (A-E) are above the 95% ANI cut-off to be considered the same species (Table 2). Almost all within-clade pairwise ANIs are greater than between-clade ANIs (Table 2) and all genomes within a clade had the highest pairwise ANI to the type strain for that clade, supporting that these are distinct subspecies. Based on hsp60 sequences, clade A containing the E. xiangfangensis type strain is Hoffmann cluster VI; clade B containing the E. hormaechei subsp. steigerwaltii type strain is Hoffmann cluster VIII; clade C containing the E. hormaechei subsp. oharae type strain is also Hoffman cluster VI; clade D containing the Hoffmann cluster III type strain (proposed name E. hormaechei subsp. hoffmannii subsp. nov.) is Hoffmann cluster III; and clade E containing the E. hormaechei subsp. hormaechei type strain is Hoffmann cluster VII.

To explore the gene content differences of the E. cloacae complex and the E. hormaechei subspecies in particular, the pan-genome of the 1216 E. cloacae complex genomes was determined using PanOCT³⁵. The pan-genome generates orthologous gene clusters that delineate which genes are in common between the clades and which genes differentiate the clades (Supplemental Table 3 and Supplemental Table 4). There were 2966 genes in “common to all” of the clades (present in 90% of the genomes of each clade). The number of genes “specific to” a clade (present in 90% of the genomes of that clade and in less than 10% of genomes from any other clade) varied from 0 (L) to 465 (V). The number of genes “missing from” a clade (present in less than 10% of the genomes of that clade and present in at least 90% of the genomes of all other clades) varied from 0 (A,C,H,K,O) to 40 (U). While ANI is the primary determinant of drawing distinctions between species and subspecies, gene content plays a role in generating phenotypic differences which might rationalize segregating a clade from species into subspecies. The clades which represent named species and subspecies show no qualitative difference in gene content from clades with no named species (Supplemental Table 4). In particular, clade D which is the proposed E. hormaechei subsp. hoffmannii has more genes specific to it than 3 of the 4 recognized subspecies. The gene content numbers need to be looked at carefully since they depend on the number of genomes in a clade (T has 187 clade specific genes but this is based on a single genome which means it is really strain specific genes rather than species specific), the distance from other clades (V the most distant clade has 465 specific genes and also has only 3 genomes), and sampling bias such as if most genomes in a clade are from a clonal outbreak. ANI appears to have less of these subjective issues to deal with.

Biochemical and other properties of the E. hormaechei subspp. clades have been previously published^3,5 except for clade D. With the availability of whole genome sequences and pan-genome analysis tools some of the observed phenotypic traits can be assigned to genetic features, such as the presence or absence of protein coding genes for known metabolic pathways. E. hormaechei subsp. hormaechei was previously distinguished from E. hormaechei subsp. oharae and E. hormaechei subsp. steigerwaltii by growth on dulcitol (a.k.a. galactitol) as the sole carbon source³. This phenotype can be explained by the presence of a gat operon^7,36 exclusively within the hormaechei subsp.. Also, in the same genomic location, between the D-galactarate dehydratase gene and the 16S rRNA methyltransferase gene, the steigerwaltii, oharae, xiangfangensis, and hoffmannii subspp. have a related, but different operon, encoding for N-acetyl galactosamine metabolism (a.k.a., the aga locus)^7,37. Similarly, steigerwaltii isolates can be distinguished from hormaechei, oharae, xiangfangensis, and hoffmannii by their ability to grow on adonitol (a.k.a. ribitol) and D(+)-arabitol; both 5 carbon sugar alcohols known as penitols. The rbt and dal operons known from Klebsiella aerogenes, which metabolize ribitol and D(+)-arabitol respectively^7,38, account for this difference and are found almost exclusively in steigerwaltii. E. hormaechei subsp. hoffmannii has 25 clade specific genes 10 of which (clusters 28856-28865 Supplemental Table 3) occur as a unit between core clusters (16694-5) and another 6 (15153-15156, 27141-2) occur between core clusters (17653-4). These clusters have no or vague annotation but are intriguing targets to provide functional phenotypic differences.

Methods

MASH³² is a very fast tool for determining approximate pairwise ANI values given sequenced genomes. A PERL script was used to invoke the following command to generate a set of MASH (version 2.0) sketches of k-mer size 16 for the 1249 downloaded Enterobacter genomes:

mash sketch -k 16 -o Enter.Sketch.file [List of the Genomes]

The resulting sketches file was then used to compare all the genomes against each other with an additional PERL script which calls MASH (version 2.0) with the command:

Mash dist Enter.Sketch.file [List of the Genomes]

which generated data that could be extracted into an all versus all ANI comparison (Supplemental Table 5). We used the GGRaSP R package (version 1.0) which generated an ultramateric tree by using the R hclust function with average linkage from the distance matrix calculated by subtracting 100 from the MASH ANI results. The result was translated into Newick format with the APE³⁹ R package (Supplemental File 1) rendered with metadata annotated using the Interactive Tree of Life⁴⁰ into Figure 1.

Based on the tree 30 genomes were deemed to be outliers and probably not in the Enterobacter genus as well as 2 E. lignolyticus genomes and 1 E. timonensis genome deemed to be outside of the E. cloacae complex. These 30 genomes were compared to all genome sequenced bacterial type strains from NCBI RefSeq (Supplemental Table 2) using MASH which confirmed that these genomes were likely misnamed as Enterobacter. The decision to leave E. lignolyticus and E. timonensis outside of the E. cloacae complex was based on two reasons: historically neither has been included in the complex, and there is a quantitative difference in the mean ANI values between genomes of these two species and genomes included in the 22 clades within the complex (last two rows of Table 2). The highest mean ANI for E. lignolyticus and E. timonensis to genomes included in the 22 clades within the complex is 86.2% for E. timonensis to clade S; whereas, the lowest mean ANI within the complex is 86.5% between clades P and U.

From the all versus all MASH ANI comparison GGRaSP was used to generate average linkage clusters and the medoids of those clusters at both the 95% (species) and 98% (subspecies) levels. If type strains existed at the subspecies level those clusters were used (E. hormaechei and E. cloacae) otherwise species level clusters were used resulting in 22 clades (A-V). If a type strain genome sequence existed for a clade it was selected otherwise the medoid was selected as a proxy. The one exception for this was clade J where two different type strains existed: E. asburiae and E. muelleri where both were retained for the typing. These 23 representative genomes were used to “type” all 1216 Enterobacter cloacae complex genomes (Supplemental Table 1). For typing the best MASH ANI match was used and resolved to either the species or subspecies level. As expected the typing was in complete agreement with the clades in the MASH ANI tree (Figure 1). The MASH sketches for these 22 clade representatives (after removing the redundant E. muelleri) can be used as a fast categorization tool for novel Enterobacter cloacae complex genomes.

GGRaSP was similarly used to select the 250 most diverse genomes including the outliers from the 1249 downloaded genomes while eliminating very closely related genomes. PanOCT^35,41 run at the nucleotide level was used to generate the orthologous clusters for a pan-genome. The primary use of this was to validate the approximate MASH ANI values. PanOCT determines pairwise ANI values by looking at every orthologous cluster shared by a pair of genomes. The percent identity of each match is weighted by the length of the match, summed over all relevant clusters, and divided by the sum of match lengths which is consistent with previous calculations of ANI. Supplemental Figure 1 shows that the MASH ANI estimate is very strongly correlated (98.9) with the PanOCT ANI measurement. For PanOCT ANI values greater than 94% the estimate is very tight (mean error 0.34±0.22) versus less than 94% (1.15±0.70). The clades and tree at the clade level remained the same using PanOCT ANI values.

For the PanOCT run with 1,216 genomes to determine gene content similarities, PanOCT was run as part of the JCVI pan-genome pipeline in hierarchical fashion with the following batches of genomes run by PanOCT at level 1: (combined 3 E. mori, 3 E. soli, 8 E. cancerogenus, 8 E. cloacae complex clade K, 13 E. cloacae complex clade L, 11 E. cloacae complex clade N, 4 E. cloacae complex clade O, 4 E. cloacae complex clade P, 5 E. cloacae complex clade S, 1 E. cloacae complex clade T); (combined 45 E. cloacae subsp. cloacae, 9 E. cloacae subsp. dissolvens); (randomly split into 4 groups 169 E. hormaechei subsp. hoffmannii); (7 E. hormaechei subsp. hormaechei); (68 E. hormaechei subsp. oharae); (randomly split into 8 groups 325 E. hormaechei subsp. steigerwaltii); (randomly split into 6 groups 255 E. hormaechei subsp. xiangfangensis); (78 E. asburiae); (30 E. bugandensis); (71 E. kobei); (29 E. ludwigii); and (70 E. roggenkampii). The level 1 clusters were then combined using PanOCT at level 2 and the final output generated using the PanOCT (version 3.27) command line:

panoct.pl -R matchtable.txt -f genomes.list -g combined.att_file -P combined.fasta -b final_panoct_run -c 0,95

The diverse 250 genome PanOCT run and the level 1 PanOCT batch runs used the PanOCT (version 3.27) command line:

panoct.pl -b results -t combined.blast -f genomes.list -g combined.att -P combined.fasta -S yes -L 1 -M Y -H Y -V Y -N Y -F 1.33 -G y -c 0,50,95,100 -T

The hierarchical PanOCT run of 1,216 genomes produced a matrix of orthologous gene clusters (Supplemental Table 3) where the rows are clusters and the columns are genomes with the cells containing the RefSeq IDs for the gene from the corresponding genome. This matrix was used to determine genes common to all, specific to, and missing from clades A-V. Individual PanOCT runs were also done for clade J, D, and M. Clade J to insure that PanOCT ANI values confirmed MASH ANI values that E. asburiae and E. muelleri are the same species which they did and these ANI values were used to determine the 8 subclades at 98% ANI using hierarchical clustering (hclust in R) average linkage. Clade D to confirm the MASH ANI values for E. hormaechei subsp. hoffmannii which they did. Clade M was done likewise to confirm E. roggenkampii which they did.

Discussion

The Background section reviews how the tools for defining a species have evolved. In a recent review of the genus Mycobacterium, the authors proposed that any newly defined bacterial species must have a genome sequence and an ANI comparison carried out against existing sequenced type strains to justify a novel species assignment⁴². ANI analysis should not be relied on in isolation for defining a species since historical or clinical phenotypic distinctions may be important for example in distinguishing between E. coli and Shigela which by ANI are the same species. However, genome sequencing appears to be outstripping the taxonomic definition of species within some genera. For the 22 clades of the E. cloacae complex identified here 9 do not have named type strains (7 if the two proposed here are adopted). For important pathogens where clinical practice may rely on proper classification the ability to name these clades/species and provide resources for identifying them could be pivotal. Unfortunately, the current established journal for validly publishing bacterial species’ names insists on phenotypic characterization and deposition of the type strain before naming is valid. This prevents computational based methods from moving quickly. Paradoxically almost all species identifying diagnostic tests are genotype not phenotype based so genotype is good enough for diagnosis but not species definition. Further, delineating what is acceptable to define as a new species is also genotype not phenotype based whether via DDH, marker genes, or more recently ANI. Worse there are no published standards for what defines the minimal set of phenotypic biochemical assays that must be performed. As the Mycobacterium review authors state: “The easy and affordable availability of reliable whole-genome sequences raises doubts about the real added value of investigating phenotypic traits when a new species is described. Actually, different taxonomists use their own panels of tests, often not standardized, to produce results of no use for colleagues and absolutely incomprehensible to the community of mycobacteriologists who have dismissed such approach since the ‘90s. For the genus Mycobacterium the major phenotypic traits that cannot be disregarded should include growth rate and pigmentation of colonies, while the classical investigation of biochemical activities is clearly obsolete.”. If there were accepted standards for minimal phenotypic characterization then culture collection repositories could choose to provide the characterization as fee for service or even for free for type strains as an incentive for deposition. With the rapid growth in synthetic genomics capabilities one could argue that the deposition of a high quality complete genome might suffice rather than a culture. We propose allowing “placeholder” species or subspecies names such as “E. cloacae complex clade S” in order to enable the most robust use of computational and genomic resources for clinical diagnosis while awaiting the designation and deposition of a type strain with a valid name possibly with some minimal phenotypic characterization.

Conclusion

Computational analysis supports the reassignment of E. xiangfangensis to E. hormaechei subsp. xiangfangensis. We propose to name clade D / Hoffman cluster III as E. hormaechei subsp. hoffmannii in honor of Harald Hoffmann’s work elucidating the phylogenetic structure of the E. cloacae complex² in particular the subspecies of E. hormaechei³. We propose to name clade M / Hoffmann cluster IV Enterobacter roggenkampii after Andreas Roggenkamp for his work on elucidating the phylogenetic structure of the E. cloacae complex². The analysis also shows that E. muelleri³³ is a later heterotypic synonym of E. asburiae³⁴ which should take precedence.

Data availability

All data underlying the results are available as part of the article and no additional source data are required