Whole genome sequences of Pseudomonas strains with species designation excluding those of P. aeruginosa were downloaded from the NCBI FTP server (ftp.ncbi.nih.gov) in February, 2016. All genomes were filtered for assemblies with contig number no greater than 300 to avoid the inclusion of overly fragmented genome into the analysis. In addition, ten Acinetobacter, four Cellvibrio species and two Azotobacter species whole genomes were also included as the outgroup for phylogenomic analysis.

Whole proteome was predicted using Prodigal V2.6.2 (Hyatt et al., 2010) and piped into PhyloPhlAn 0.99 (Segata et al., 2013) to construct a maximum likelihood tree using FastTree2 (–bionj –slownni –mlacc 2 –pseudo –spr 4 setting) (Price et al., 2010) based on the identification, alignment and concatenation of 400 universally conserved proteins (Edgar, 2004, 2010). Local branch support values were computed by FastTree2 using the Shimodaira-Hasegawa test. The tree was subsequently visualized and manually collapsed based on genospecies (Supplemental Table 1) using MEGA6 (Tamura et al., 2013).

Pairwise ANIm was calculated with Jspecies v1.2.1 using the standard MUMmer algorithm (Richter and Rosselló-Móra, 2009). To reduce computational calculation, we separated the ANIm calculation into 4 groups based on phylogenomic clustering e.g., Clades 1, 2, 3 and strains that are not part of the labeled clade (Figure 1). The matrix obtained from Jspecies for each major group was clustered and visualized using Rstudio 0.99.902, pheatmap package.

Nucleotide similarity search on housekeeping genes (rpoB, recA, rpoD, and/or gyrB) was performed with basic local alignment search tool (BLASTN) using an e-value cutoff of 1 × e⁻¹⁰ (Altschul et al., 1990).

A total of 373 Pseudomonas genome sequences were selected for phylogenomic tree reconstruction. The genome size among Pseudomonas strains is highly variable, ranging from 3,022,325 bp (P. caeni DSM 24390) to 8,608,769 bp (P. bauzanesis W13Z2). Similarly, the GC content of their genomes ranges from 56.51% (P. endophytica BSTT44) to 67.43% (P. syringae pv tomato T1). Based on current classification, the most commonly deposited non-P. aeruginosa pseudomonad species is P. syringae (N = 67), followed by P. fluorescens (N = 62), P. putida (N = 28), P. stuzeri (N = 17), P. amygdali (N = 15), P. psychrotolerans (N = 14), P. savastanoi (N = 13), and P. denitrificans (N = 11). Out of a total of 79 species outcomes, surprisingly 43 species have only one whole genome sequence representative in the database.

Maximum likelihood tree as inferred from 400 universal conserved proteins clustered the Pseudomonas strains into several major clades with maximal SH-like branch support (Figure 1 and Supplemental Figure 1). Support values of some inner clades are low suggesting that shallow relationships between strains within the same genospecies are not resolved, presumably due to the lack of divergence at the amino acid level among closely related strains. A wgMLST analysis restricted to strains from closely related Pseudomonas species will likely identify more conserved loci which can then be used to refine their evolutionary relationships which unfortunately is beyond the scope of this study. Such wgMLST analysis can be performed using existing bioinformatic tools such as BIGSdb, Roary and PGAdb-builder (Jolley and Maiden, 2010; Page et al., 2015; Liu et al., 2016). Strikingly, several whole genome sequences labeled as the same Pseudomonas species were located in different major clades or distantly apart taxa (Supplemental Figure 1). For example, while P. protegens CHA0^T clustered within Clade 2, P. protegens 231 PPRO formed a monophyletic group with other sequences labeled as P. denitrificans in a distant clade. Similarly, P. marginalis ICMP 11289 belongs to Clade 1 and clusters with other P. viridiflava whereas P. marginalis ICMP 9505 belongs to Clade 2 (Supplemental Figure 1). The lack of taxonomic congruence as reflected by the inconsistent phylogenetic clustering among Pseudomonas strains with the same species designation raises suspicion about the species validity of Pseudomonas genomes submitted to public database.

Using 95% ANIm as the cutoff point for species delineation, a total of 145 genomic clusters were formed and when possible, assigned to a valid Pseudomonas species (Figure 1 and Supplemental Table 1). Two-hundred and seven genomes were wrongly assigned at the species level as evidence by the lack of genomic clustering with their expected type strain genome and/or low nucleotide similarity (<95%) to the type strain Multilocus Sequence Typing (MLST) sequences.

All of the major species contain at least one member with incorrect species classification except for P. amygdali (Figure 2). Surprisingly, all 13 P. savastanoi genomes, majority of which originated from a single study (Mott et al., 2016) may have been misclassified as they were clustered into GS-1 containing P. amygdali ICMP3918^T (Figure 2 and Supplemental Table 1). ANIm matrix of 113 whole genome sequences within Clade 1 (Figure 1) containing a majority of P. syringae identified at least 14 different genomic clusters (Supplemental Figure 2). The presence of whole genome sequence for the P. syringae type strain KTCC 12500^T = DSM 10604^T in GS-4 suggests that only 23 out of 67 P. syringae genomes were correctly assigned (Supplemental Figure 2), indicating that more than 50% of the currently deposited P. syringae genomes should be taxonomically revised. At least 10 P. syringae genomes have to be assigned to other valid Pseudomonas species (Supplemental Table 1). For example, P. syringae CC1513 should be re-classified as P. coronafaciens given its high pairwise ANIm of 99.2% to P. coronafaciens LMG5060^T.

Prior to the genome availability of P. fluorescens DSM 50090^T, P. fluorescens strain WH6 has served as the reference genome for the taxonomic assignment of P. fluorescens (Duan et al., 2013; Feng et al., 2015). Unfortunately, ANIm calculation shows that strain WH6 belongs to GS-35 with a pairwise ANIm of only 87.94% to the recently available P. fluorescens DSM 50090^T genome (Supplemental Table 1 and Supplemental Figure 3). In another study, strain WH6 was also shown to be more similar to P. azotofomans LMG 21611^T than the representative species P. fluorescens DSM 50090^T using the MLSA method (Gomila et al., 2015). This finding along with our results here, suggest that the designation of P. fluorescens WH6 genome as the reference genome might have led to the subsequent mislabeling of newly-sequenced Pseudomonas genomes. As expected, out of a total of 62 whole genome sequences that were assigned as P. fluorescens, none of them belongs to the same genospecies as P. fluorescens DSM 50090^T. By re-evaluating the species designation based on ANIm clustering against other type strains of Pseudomonas, 22 of the 61 wrongly identified P. fluorescens strains were successfully re-assigned to the correct valid species name leaving 39 pseudomonad genospecies in Clade 2 (Figure 1) unassigned at the species level.

Only 1 strain of the total 28 P. putida strains belongs to the same genospecies with the P. putida type strain NBRC 14164^T (GS-102) which was similarly supported by their monophyletic clustering in the phylogenomic tree (Clade 3 in Supplemental Figure 1). Some of the mislabeled P. putida strains might be novel species given its independent formation of genospecies e.g., GS-103 in P. putida MTCC5279 (Supplemental Table 1 and Supplemental Figure 4).

With the current availability of defined genospecies within the genus Pseudomonas, we re-evaluated the taxonomic status of 5 Pseudomonas strains previously isolated from interior vine tissue of poison ivy by our group (Tran et al., 2015). Strain RIT-PI-g could be assigned to the species Pseudomonas libanensis forming GS-19 along with the P. libanensis DSM 17149^T corroborating our previous species assignment based on similarity to various housekeeping genes. Pseudomonas. sp. RIT-PI-q and P. sp. RIT-PI-r both belong to a single-member genospecies (GS-69 and GS-52, respectively), further underscoring the high genomic diversity of plant-associated Pseudomonas that remains to be explored and described with extensive taxon sampling effort. On the contrary, Pseudomonas sp. RIT-PI-o and P. sp. RIT-PI-a were placed in the same genomic cluster (GS-13) as the incorrectly labeled P. fluorescens MEP34 and Pseudomonas sp. leaf129 (Supplemental Table 1 and Supplemental Figure 5). Both P. fluorescens MEP34 and P.s sp. leaf129 were isolated from leaf of Arabidopsis thaliana, a popular flowering plant in North America that has been adopted as the flowering plant genetic model species (Bai et al., 2015; Perisin et al., 2015). Such similar plant origins, e.g., poision ivy vine and Arabidopsis, poses interesting links between bacterial genomic characteristics and plant-hosting capacity.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.