The International Commission for the Protection of the Danube River organized the JDS (Peter et al., 2017) campaign (Kittinger et al., 2016b). Sixty-eight water samples were collected from August 12 to September 26 of 2013. The sampling locations studied in this study were 14 sites along the river: JDS02, JDS03, JDS08, JDS10, JDS22, JDS28, JDS36, JDS38, JDS49, JDS57, JDS59, JDS63, JDS067, JDS68, and are indicated in Supplementary Table S1. Water samples (0.5 ml) were plated on Endo Agar, Lysine Deoxycholate Agar and Chromocult Coliform Agar and were incubated at 37^oC for 18–24 h as described (Kittinger et al., 2016b).

All Pseudomonas strains (611) isolated in the Joint Danube Survey 3 (JDS3) campaign identified by MALDI-TOF MS (Kittinger et al., 2016b) and stored at the collection of the Institute of Hygiene, Microbiology and Environmental Medicine, Medical University of Graz isolated at the JDS3 campaign were sent to the Microbiology laboratory of the University of the Balearic Islands (UIB). They were checked for viability and purity after culturing on LB medium (Conda) at 30°C for 24–48 h. To standardize procedures at the UIB laboratory the strains were cultured routinely at 30°C.

An UPGMA dendrogram using the MALDI-TOF MS protein profiles of the 611 strains was constructed with the raw Biomérieux data obtained with the VITEK MS identification system in order to group them in putative species. To obtain the spectra, a colony from blood agar plate (incubated at 37°C) was directly spotted on the MALDI plate, and then overlaid with 1 μl of matrix solution and air-dried. The loaded plate was then placed in the instrument according to the manufacturer’s instructions (Kittinger et al., 2016a). The list of peaks obtained from each mass spectrum was converted to a matrix of 0 and 1, according to the absence or presence of a specific peak. Peaks that differed in less than 1 unit of m/Z between them were considered the same protein. From a list of 6,870 peaks, only 747 different major proteins were considered. The matrix was used to calculate the Pearson distance between all pairs of data and to construct an UPGMA dendrogram. Hundred and ninety strains were selected as representatives of the 53 clearly defined clusters of strains that were at least 69% similar in the dendrogram.

Strains were cultured overnight in 4 ml LB (Conda) at 30°C. The DNA extraction, PCR amplification, primers used, purification of the amplified products and DNA sequencing conditions, as well as the sequence analysis procedures, have been previously described (Mulet et al., 2012). PCR conditions for the 16S rRNA gene are described in Mulet et al. (2008). Identifications were performed at the EzTaxon database.

The partial sequences of the rpoD gene (708 nucleotides) of 190 strains representatives of the 53 clusters were aligned and compared with the 220 Pseudomonas species type strains rpoD nucleotide sequences of an indoor database at the laboratory of Microbiology of the University of the Balearic Islands. The database was updated with all Pseudomonas species type strains described until January 2020. The neighbor joining (NJ) (Saitou and Nei, 1987), maximum likelihood (ML) (Felsenstein, 1981) and maximum parsimony (MP) algorithms were used for the comparisons with the integrated tool MEGA5 (Nei and Kumar, 2000; Tamura et al., 2011). A neighbor joining tree (Saitou and Nei, 1987) based on the rpoD nucleotide sequences of the 190 representative strains isolated from the Danube River was constructed. Detailed analysis of the specific Pseudomonas groups of species, like P. putida or P. stutzeri, included also additional strains to identify more precisely the strains or to propose them as putative new species. Each different sequence was considered a single phylotype, and several phylotypes were grouped in a single phylospecies when they affiliated in the same phylogenetic branch and the similarity in the rpoD gene sequence was equal or greater than 95% (Mulet et al., 2009; Scotta et al., 2013). Phylospecies was abbreviated PS.

Sequence information obtained was used to calculate phylospecies coverage with the formula C=1-(nN) where C is the coverage index, n the number of singletons (phylospecies with only 1 representative) and N the total number of strains. Rarefaction curves were constructed with the PAST software package version 3.20 (Hammer et al., 2001).

The 611 isolates were classified in 53 similarity clusters by their MALDI-TOF MS protein profiles as depicted in Figure 1, Table 1, and Supplementary Figure S1. Isolates in each cluster were at least 69% similar in their protein profiles and each cluster was numbered, from 1 to 53. Six groupings comprised a high number of isolates (from 54 to 137), while the others included only 1 to 20 isolates. Isolates identifications with the VITEK MS system were coincident in most cases with the main clusters: groups 1–3 (7 strains) to P. stutzeri; groups 4–26 (401 isolates) to P. putida; groups 27–33 (161 isolates) to P. fluorescens; groups 34–39 (13 isolates) to P. fluorescens or P. viridiflava; groups 41–43 (16 isolates) to P. mendocina or P. oleovorans; less abundant groups 44–53 (12 isolates) included isolates identified in 5 different species. Some of the more abundant groupings, based on the rpoD gene sequence analyses could be identified in more than one phylospecies and were subdivided later in subgroups labeled a, b, c, d, and e.

Hundred ninety strains were selected as representatives of each MALDI-TOF MS cluster for their rpoD gene nucleotide sequence analysis to assign each grouping to a Pseudomonas species, as depicted in Figure 1 and Supplementary Figure S1 and Supplementary Table S2. For each of the 53 groupings at least two representative strains were selected when possible and are indicated in the Supplementary Figure S1 and Supplementary Table S2. The partial rpoD gene sequence of 190 strains was determined (Supplementary Table S2) and, as described in Materials and Methods, sequences were aligned together with an in house rpoD gene sequence database which contained 220 sequences of Pseudomonas type strains. Data obtained are shown in Table 1 and Supplementary Figure S1. Phylogenetic relationships are represented in Figures 2, 3 and Supplementary Figures S2, S3. Strains were identified at the species level when the similarity values with a species type strain was higher than 95%, the established species cutoff.

Following the phylogenetic classification proposed by Mulet et al. (2010) and Gomila et al. (2015), 108 strains (57%) affiliated to 34 phylospecies that were assigned to Pseudomonas known species. They were mainly distributed in the P. fluorescens lineage (Figure 2), in the P. fluorescens group (G), in six different phylogenetic subgroups of species (SG): P. fluorescens SG (four species), P. fragi SG (four species), P. koreensis SG (two species), P. protegens (two species), P. corrugata SG (two species) and P. jessenii (one species); 55 strains affiliated with the P. putida phylogenetic group of species (Figure 3) and were distributed in 12 known species; 19 strains affiliated to known species in the P. aeruginosa lineage, in three different groups: six in the P. stutzeri G (in genomovars 1 and 3), eleven in the P. oleovorans G (in three species) and two in the P. aeruginosa G (P. resinovorans and P. citronellolis) (Figure 2).

The remaining 82 strains (43% of the isolates sequenced) could not be assigned to any known Pseudomonas species, because their rpoD sequence similarities were below the cutoff value 95%, and were considered representatives of 32 putative new phylospecies. The putative new species were designated “New PS” followed by the corresponding number. When the whole genome sequence of the reference strain was known, it was designated “New GS” (GS for genomospecies) and the number of the reference strain was indicated in brackets. Strains of the putative new species were located in the P. fluorescens lineage, in the P. fluorescens G (P. koreensis SG, P. corrugata SG, and P. protegens SG), in the P. syringae G, and in the P. putida G. Strains in the P. aeruginosa lineage were distributed in the P. aeruginosa and P. oleovorans groups (Figures 1–3 and Table 1). Several groups or subgroups contained known and putative new species: P. protegens SG, P. koreensis G, P. corrugata SG, P. aeruginosa and P. putida G, P. oleovorans G.

Combining the identification based on the rpoD sequence and the groupings of strains in the MALDI-TOF MS analysis, we were able to classify 400 strains: 262 of them were assigned to known species, 65 to phylospecies for which a representative strain is well known, and 73 strains were assigned to new phylospecies not yet described. For instance, 31 of the 405 strains identified as P. putida with the VITEK MS system clustered in group 5a of the MALDI-TOF dendrogram. The rpoD gene of nine representative strains of this group were sequenced and showed 95–99.5 similarity with the corresponding sequence of the P. soli type strain, a member of the P. putida phylogenetic group; therefore, the 31 strains were identified as P. soli (Supplementary Figure S1). However, we have to consider that three other strains in the MALDI-TOF groups 5b, 7, and 51 were also identified by the rpoD sequence as P. soli. To verify the protein profiles, the 12 P. soli strains were reanalyzed in a Bruker autoflex III mass spectrometer and the Biotyper database as described in Peña et al. (2016). All strains clustered together with the P. soli type strain and were separated from the other species type strains in the corresponding dendrogram (results not shown).

The most abundant known phylospecies identified were: 88 strains of P. alloputida, 34 strains of P. soli, 20 strains of P. monteilii, 21 strains of P. fulva and 9 strains of P. capeferrum in the P. putida G; 14 strains of P. moraviensis in the P. koreensis SG and 8 strains of P. bubulae in the P. fragi SG, both in the P. fluorescens G. P stutzeri was the eighth more abundant species with 7 strains, 6 of them assigned to genomovar 1, and 1 to genomovar 3 (Supplementary Table S1). From the 34 known phylospecies, 18 are represented by only 1 or 2 strains; the same situation is shown in the 32 new putative phylospecies because 23 of them are also represented by 1 or 2 strains.

The most abundant potential new phylospecies corresponded to species described in the P. putida G. They have to be considered potential new species, because their rpoD sequence similarity was below the species cutoff stablished with any type strain. One of these new phylospecies (50 strains) was represented by a well-studied strain assigned to the species P. putida: strain W619. This result suggested that these strains have to be reclassified in a species different to P. putida.

One hundred and ninety-one strains were assigned to phylogenetic groups or subgroups, but not at the phylospecies level, because several phylospecies were detected in the same MALDI-TOF grouping. Only three strains were identified in genera distinct to Pseudomonas and were located distantly in the dendrogram. One of them (branch 52) was identified as Arthrobacter koreensis and two (group 53) as Brevundimonas terrae by sequencing the 16S rDNA gene.

Rarefaction curves and phylospecies coverage indices were calculated for each sampling point and for the total set of strains identified to determine the extent of the species richness detected (Tables 2, 3). In total, 400 strains were identified at the phylospecies level (66% of the strains in the study) and were distributed in 65 phylospecies. The strains assigned to phylospecies ranged between 37 and 100% depending on the sampling point. The rarefaction curves did not reach saturation, indicating that the total diversity in the samples was not detected. Supplementary Figure S4 shows, as example, the results for five sampling sites (Supplementary Figure S4A) and for all the sampling sites (Supplementary Figure S4B). However, the species richness assessed with the coverage index per sampling point reached values between 80 and 94% in those points with a high number of isolates. Points with less than 32 strains showed coverages lower than 69%. The sampling sites with the highest richness were JDS28, JDS36, and JDS63, with 22, 18 and 15 phylospecies, respectively (Table 2). Seventy-eight strains of the 118 isolates from sampling site JDS28 were assigned to 22 phylospecies (Table 2). They showed a wide diversity of species, distributed in 6 different species of the P. putida G (17 P. alloputida, 11 P. soli, 3 P. fulva, 3 P., monteilii, 2 P. inefficax, and 1 P. asiatica), 3 species in P. fluorescens G (6 P. moraviensis, 1 P. koreensis, and 1 P. saponiphila), one species in the P. aeruginosa G (P. citronellolis), one species in the P. oleovorans G (1 P. toyotomiensis/P. chengduensis) and in 11 new phylospecies. The remaining 40 strains were assigned to four different phylogenetic groups.

The rarefaction curve was also calculated for the whole course of the Danube River by considering the 400 strains assigned to phylospecies. Results indicate that 95% of the phylospecies diversity in the set of 611 strains was detected, but that a high Pseudomonas species richness present in the Danube River, at least 50% remains to be explored.

The previous data obtained for multidrug resistant Pseudomonas species (MDR) identified the strains in the species P. putida and P. fluorescens (Kittinger et al., 2016b). The rpoD sequence analysis, however, indicated that these assignations are correct at the level of phylogenetic group, P. fluorescens G, and P. putida G as defined by Mulet et al. (2010), but they can be now identified more precisely at the phylospecies level. The P. fluorescens group contains at the moment more than sixty different species and the P. putida group is composed of at least by 20 species (Keshavarz-Tohid et al., 2019; Peña et al., 2019). Table 4 shows the assignation by rpoD of the 14 MDR strains resistant to three or more antibiotic classes. They were assigned to seven phylospecies. Seven MDR strains were assigned to six known species, but the predominant phylospecies was New PS-6 with five isolates, that clustered in the MALDI-TOF dendrogram in group 27a composed by 45 strains. The closest species in the rpoD sequence is P. protegens in the P. protegens SG. New PS-6 strains could harbor the combination of resistances to three (MDR3) or four (MDR4) antibiotic classes and were distributed in three different sampling points: JDS02 (downstream Regensburg), JDS22 (downstream Budapest) and JDS28 (upstream Drava). The five strains of the New PS-6 are resistant to the antibiotics meropenem and ceftazidime, three of them also to imipenem, two of them to ciprofloxacin and the remaining three strains also to one more antibiotic (piperacilin/tazobactam, cefepime or gentamicin). These result indicated a broad spectrum of different combination of antibiotics resistances in the same phylospecies. Non-MDR strains of these seven phylospecies have been also isolated from the same or different sampling points. It has to be pointed out that, to our knowledge, none of the species detected have been found previously in clinical specimens, and therefore have to be considered environmental.

The VITEK MS identification system clearly assigned the strains to the correct phylogenetic group of species, but failed at the species level in many strains. Forty-one strains were assigned by MALDI-TOF to a species different from the identification obtained by the combined MALDI-TOF/rpoD sequencing procedure. Additionally, other 42 strains have to be considered possible new species under a phylogenetic point of view. This significant discrepancy is probably due to the extent of the VITEK MS database, designed mainly for the identification of clinically relevant species and has to be updated with environmental species, not included in the commercial database. Other authors, even in the clinical laboratory (McElvania TeKippe and Burnham, 2014), have also addressed the difficulties of the MALDI-TOF MS system for the identification of species in the P. putida phylogenetic group. However, a fair good concordance was found in the assignation of strains in the 53 MALDI-TOF groupings and the corresponding phylogenetic group or subgroup. Three possible situations were detected:

At least 73 strains were classified in 25 phylospecies of not yet described bacterial species (labeled New-PS) because their rpoD gene sequences were less than 95% similar to species type strains in the genus: 14 in the P. putida, 6 in the P. fluorescens, 2 in the P. syringae, 2 in the P. oleovorans and 1 in the P. aeruginosa phylogenetic groups. These possible new species merit deeply taxonomic studies.

We analyzed in a previous publication the diversity of Pseudomonas in the Woluve River by analyzing the rpoD sequences from cultivated strains and by pyrosequencing rpoD amplicons from total DNA extracted from water samples (Sánchez et al., 2014). Two hundred and forty-six possible new phylospecies were detected in the amplicon analyses and were not found in the cultivated strains. Interestingly, two Danube River strains are 97% identical to two of the 246 putative novel phylospecies detected in the Woluve River by pyrosequencing an rpoD amplicon: strains JDS49PS004 (New PS-1) and JDS02PS014 (New PS-19), confirming in this way the value of the rpoD sequences to detect Pseudomonas diversity and predict new species by culture independent methods.

Correct and precise species identification is also a prerequisite for the assessment of antimicrobial susceptibilities of the isolates, because the intrinsic, or basal, resistance will vary significantly among species and even different strains. MALDI-TOF identified the 14 MDR strains as members of the species P. putida or P. fluorescens. However, the rpoD sequence analysis demonstrate that they belong to seven different Pseudomonas species, five known (P. protegens, P. soli, P. inefficax, P. asiatica, and P. bubulae) and two other possible new species (PS-6 and PS-16) all in the P. fluorescens or P. putida phylogenetic groups of species. For all of them, non-MDR strains were also detected, which suggests a process of resistance acquisition. These strains belong to species not described as pathogens and represent an environmental reservoir of antibiotic resistance genes that can be transferred to other species. Other authors have found a MDR strain of P. aeruginosa (opportunistic pathogen) in a water sample of the Woluve River (Belgium) (Pirnay et al., 2005) but also indicating that the majority of the Woluve strains were antibiotics susceptible. It is worth to mention, that we did not find any P. aeruginosa strain in the Danube River isolates, an abundant species in aquatic environments and considered indicator for water risk assessments (Mena and Gerba, 2009).