Genomic Analysis of Three Cheese-Borne Pseudomonas lactis with Biofilm and Spoilage-Associated Behavior

Psychrotrophic pseudomonads cause spoilage of cold fresh cheeses and their shelf-life reduction. Three cheese-borne Pseudomonas sp., ITEM 17295, ITEM 17298, and ITEM 17299 strains, previously isolated from mozzarella cheese, revealed distinctive spoilage traits based on molecular determinants requiring further investigations. Genomic indexes (ANI, isDDH), MLST-based phylogeny of four housekeeping genes (16S rRNA, gyrB, rpoB and rpoD) and genome-based phylogeny reclassified them as Pseudomonas lactis. Each strain showed distinctive phenotypic traits at 15 and 30 °C: ITEM 17298 was the highest biofilm producer at both temperatures, whilst ITEM 17295 and ITEM 17299 showed the strongest proteolytic activity at 30 °C. A wider pattern of pigments was found for ITEM 17298, while ITEM 17295 colonies were not pigmented. Although the high genomic similarity, some relevant molecular differences supported this phenotypic diversity: ITEM 17295, producing low biofilm amount, missed the pel operon involved in EPS synthesis and the biofilm-related Toxin-Antitoxin systems (mqsR/mqsA, chpB/chpS); pvdS, required for the pyoverdine synthesis, was a truncated gene in ITEM 17295, harboring, instead, a second aprA involved in milk proteolysis. This work provided new insight into the food spoiler microbiota by identifying these mozzarella cheese spoilers as P. lactis; molecular targets to be exploited in the development of novel preservative strategies were also revealed.


Introduction
The genus Pseudomonas comprises species naturally widespread in the environment, including human pathogens, such as the clinically relevant P. aeruginosa [1,2]. Pseudomonas species displayed a large array of ecological functions, i.e., plant-growing promoter, bioremediation, biodegradation, and biosorption, conferring them high ecological and biotechnological importance [3][4][5]. Its taxonomy was lastly updated in 2018 by Peix et al. [6] who listed further 70 species described in the last few years in addition to the 128 previously described the same authors in 2009 [7]. However, genus taxonomy is currently discussed; indeed, according to the List of Prokaryotic Names with Standing in Nomenclature (https://lpsn.dsmz.de/genus/pseudomonas), the genus Pseudomonas (family Pseudomonadaceae, class Gammaproteobacteria) comprises 278 child taxa with a validly published name, including synonyms, among which 114 were represented by at least one type strain (WFCC Global Catalogue of Microorganisms; http://gcm.wfcc.info/overview/).

Phylogenetic Analysis
Sequences of the four housekeeping genes (16S rRNA, gyrB, rpoB and rpoD) were extracted from each genome by BLASTn search and compared to type strains sequences retrieved from Gomila et al. [10] and Mulet et al. [11]. The alignment of the concatenated dataset was performed using MUSCLE implemented in MEGA version 10.1.7 [39]. Phylogenies was inferred by using the Maximum-Likelihood (ML) method and Tamura-Nei model [40]. The tree was graphically generated by iTOL version 5.5 [41].
Genetic divergence was calculated by the ANI/AAI calculator [42,43] which estimates the average nucleotide/amino acid identity (ANI/AAI) using both best hits (one-way ANI) and reciprocal best hits (two-way ANI) between genomic datasets. The Genome-to-Genome Distance Calculator (GGDC) [44,45] web service was used to report digital DDH for the accurate delineation of prokaryotic subspecies and to calculate differences in G+C genomic content (available at ggdc.dsmz.de). Formula 2 alone was used for analysis, providing an estimation of DDH independent of genome lengths, as recommended by the authors of GGDC for use with any incomplete genomes [44,46].
Genome-based phylogeny was constructed by using The Phylogenetic Tree Building Service implemented in Patric platform (www.patric.org) considering the Cellvibrio japonicus Ueda 107 as outgroup according to Gomila et al. [10] and Maximum Likelihood method RAXML with progressive refinement [47]. For each strain, genomic sequences deposited in GenBank were used for the analysis.

Phenotypic Analysis under Two Temperatures of Incubation
Due to the high versatility and adaptation mechanisms of Pseudomonas spp., phenotypic traits were investigated both at the optimum growth temperature (30 • C) and at lower temperature (15 • C) in order to mime cold storage conditions that could favor their persistence and spoilage traits in dairy industry.

Static Biofilm Formation and Motility Assays
Biofilm formation was assayed in 96-well microtiter plates and quantified as described by O'Toole [48]. Briefly, overnight cultures of Pseudomonas spp. ITEM 17295, ITEM 17298, and ITEM 17299 were diluted 1:100 into fresh LB (100 µL; 8 replicates for each time point) and incubated at 15 and 30 • C for 192 h. At 24,48 h,72 h,and 192 h, planktonic cells were carefully removed, and wells were washed twice with distilled water; biofilm cells adhering to the bottom and side of each well were stained with crystal violet (CV; 0.1%, w/v). After a second washing step, biofilm-associated crystal violet (CV) was solubilized with 30% acetic acid (v/v) and its optical density was measured at 570 nm. Swarming and swimming motility were performed in Petri dishes (polystyrene, diameter of 50 mm) containing 10 mL of LB solidified with 0.5 and 0.3% (w/v) of agar, respectively. Swim and swarm plates were inoculated with 2.5 µL of bacterial broth culture representing approximately 1 × 10 8 CFU/mL (corresponding to 0.3 OD λ = 600 nm ; [17]) as reported by Quintieri et al. [22]. All plates were incubated at 15 • C and 30 • C. The diameters of the swarming and swimming motility zones were measured at 24, 48, and 72 h of incubation. By contrast, twitching motility was evaluated on M8 minimal medium [49] supplemented with 1 mM MgSO 4 , 0.2% glucose (w/v), 0.5% casamino acids (w/v) and 1% agar (w/v). Bacterial cells were inoculated at the bottom of the agar-dish interface. The plates were incubated at 15 • C and 30 • C. At selected times (24, 48 and 72 h), the agar layer was carefully removed and the plates were stained with 0.1% of CV (w/v). After the washing step, the biofilm was solubilized and quantified as described above.

Exopolysaccharides (EPS) Production by Congo Red Binding Assay
Pseudomonas spp. cultures were spotted (5 µL) onto Petri dishes containing Congo red agar (10 g/L tryptone, 40 µg/mL Congo red (Sigma), 15 µg/mL Coomassie brilliant and 1.5% agar) and incubated at 15 • C for 72 h. The interaction of Congo red with EPS produced the appearance of red colonies (biofilm producers). EPS production was also detected in the swimming assay performed in LB supplemented with 40 µg/mL Congo Red and 15 µg/mL Coomassie brilliant, as described above.

Proteolytic and Lipolytic Activity
Fresh cultures of Pseudomonas spp. ITEM 17295, 17298, and 17299, routinely grown on PCA medium (Biolife) for 16 h, were patched in triplicate by toothpick onto Petri dishes with Milk agar (90% mL UHT skimmed milk added to 10% ml of sterile 1.6% agar) and Rhodamine Agar [50] using olive oil (3.3% v/v) and 2 mL of filter-sterilized (0.22 µm) Rhodamine B solution (0.1% w/v) for screening proteolytic and lipolytic activity, respectively. The agar Petri dishes were incubated at 15 and 30 • C for 48 h. The diameters of clear or orange fluorescent halos (revealed under UV light at 360 nm) surrounding strain colonies on Milk Agar and Rhodamine Agar, respectively, were registered with a caliper.

Pigment Production
Cell suspensions of each Pseudomonas spp. strain were spotted (2.5 µL; 1 × 10 8 CFU/mL) in triplicate onto King A and King B agar Petri dishes (Sigma Aldrich) and Potato Dextrose Agar (PDA; [17]) and incubated at 15 and 30 • C for 1 week; pigmented colonies were also examined under Wood's lamp.

Statistical Analyses
After assessing homogeneity variance by Levene's test (p < 0.05), a three-or two-ANOVA was carried out using R software v. 3.6.3 with Rcmdr and RcmdrPlugin.aRnova packages to examine the effect time and temperature on biofilm biomass, swarming, swimming, and twitching produced by each assayed strain. Partial eta squared values (η 2 ) were computed to analyze the effect size of the main factors and their interactions. Multiple comparisons among individual means were performed by HSD Tukey test (p < 0.05). When required, one-way ANOVA was performed followed by post-hoc HSD Tukey test (p < 0.05). In the case of unequal variances, the non-parametric Kruskal-Wallis H test was applied with Wilcoxon rank sum test and p-value adjustment Bonferroni's method using the RcmdrPlugin.EZR package.

General Features of Pseudomonas sp. Genomes
Pseudomonas spp. strains genomes were sequenced using a whole genome shotgun approach. Genomes were assembled using Spades v5.0 software for a total of 89, 117, and 137 contigs (>500 bp) for ITEM 17295, ITEM 17298, and ITEM 17299, respectively. The average GC content is of 59.8% for all strain and the assembly length of 6.3 Mb (Table 1). The quality of the assemblies is excellent with 93.7% of genome completeness and N50 of 351 kb, 147 kb and 139 kb for ITEM 17295, 17298, and 17299, respectively). The whole genome shotgun projects have been deposited at DDBJ/ENA/GenBank under the accessions WJRV00000000 (for ITEM 17295), WJRU00000000 (for ITEM 17299) and NPKB00000000 (for ITEM 17298; [31]). The versions described in this paper are WJRV01000000 (for ITEM 17295), WJRU01000000 (for ITEM 17299), and NPKB01000000 (for ITEM 17298). Since 16S rRNA is not adequate for discriminating Pseudomonas species as also previously reported [11], we performed an MLST analysis using partial sequences of four core housekeeping genes (16S rRNA, gyrB, rpoB and rpoD) of type strains of Pseudomonas [10,11] as well as close related P. fluorescens species. Cluster analysis and phylogenetic tree ( Figure 1) showed that the tested ITEM strains can be reasonably placed within the subgroups formed by the P. lactis and P. paralactis species recently described by von Neubeck et al. [12].
Amending the previous description of Pseudomonas spp. ITEM 17298 [22,24], this isolate should be correctly classified as Pseudomonas lactis. Hereinafter, these bacteria will be referred to as P. lactis.
In Supplementary Tables S1 and S2, ANI and DDH matrices of Pseudomonas spp. Genome are reported. According to ANI analysis, the closest relatives for ITEM 17295 is P. fluorescens Ps_22, with 98.57% of ANI, while ITEM 17298 and ITEM 17299 share 100% of nucleotide identity and are highly similar to P. fluorescens 11293 and Ps_22 with 98.55 and 98.54% of ANI, respectively. The highest value of DDH for ITEM 17298 and ITEM 17299 is the reciprocal 99.8%, while the highest value of DDH for ITEM 17295 is towards P. fluorescens Ps_20 (89.5%). Phylogeny based on partial sequences of four housekeeping genes (16S rRNA, gyrB, rpoB and rpoD ) was inferred by using the Maximum-Likelihood (ML) method and Tamura-Nei model [40]. Initial tree(s) for the heuristic search were obtained automatically by applying the neighborjoining and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood approach, and then selecting the topology with a superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Circles are representative of bootstrap values.
Amending the previous description of Pseudomonas spp. ITEM 17298 [22,24], this isolate should be correctly classified as Pseudomonas lactis. Hereinafter, these bacteria will be referred to as P. lactis.
In Supplementary Tables S1 and S2, ANI and DDH matrices of Pseudomonas spp. Genome are reported. According to ANI analysis, the closest relatives for ITEM 17295 is P. fluorescens Ps_22, with 98.57% of ANI, while ITEM 17298 and ITEM 17299 share 100% of nucleotide identity and are highly similar to P. fluorescens 11293 and Ps_22 with 98.55 and 98.54% of ANI, respectively. The highest Figure 1. Phylogeny based on partial sequences of four housekeeping genes (16S rRNA, gyrB, rpoB and rpoD) was inferred by using the Maximum-Likelihood (ML) method and Tamura-Nei model [40]. Initial tree(s) for the heuristic search were obtained automatically by applying the neighbor-joining and BioNJ algorithms to a matrix of pairwise distances estimated using the Maximum Composite Likelihood approach, and then selecting the topology with a superior log likelihood value. The tree is drawn to scale, with branch lengths measured in the number of substitutions per site. Circles are representative of bootstrap values.
Both phylogenetic trees (Figures 1 and 2) obtained by MLST analysis and Randomized Axelerated Maximum Likelihood (RAxML) analysis confirmed the clustering of Pseudomonas strains isolated from milk or dairy environment, validating the genetic similarity of the strains and the proposed novel species described by von Neubeck et al. [12].
value of DDH for ITEM 17298 and ITEM 17299 is the reciprocal 99.8%, while the highest value of DDH for ITEM 17295 is towards P. fluorescens Ps_20 (89.5%).
Both phylogenetic trees (Figures 1 and 2) obtained by MLST analysis and Randomized Axelerated Maximum Likelihood (RAxML) analysis confirmed the clustering of Pseudomonas strains isolated from milk or dairy environment, validating the genetic similarity of the strains and the proposed novel species described by von Neubeck et al. [12]. According to Chun et al. [8], who defined genome-based criteria for the taxonomy of prokaryotes, we can assume that Pseudomonas sp. strain ITEM 17295, strain ITEM 17298, strain ITEM 17299, strain Ps_22, strain Ps_40, strain Ps_20, strain Ps_77, strain 11293, and strain 74954 should all be classified as P. lactis. In addition, although von Neubeck et al. [12] described P. lactis and P. paralactis as different species according phenotypic characteristics, they did not perform any According to Chun et al. [8], who defined genome-based criteria for the taxonomy of prokaryotes, we can assume that Pseudomonas sp. strain ITEM 17295, strain ITEM 17298, strain ITEM 17299, strain Ps_22, strain Ps_40, strain Ps_20, strain Ps_77, strain 11293, and strain 74954 should all be classified as P. lactis. In addition, although von Neubeck et al. [12] described P. lactis and P. paralactis as different species according phenotypic characteristics, they did not perform any phylogenetic tree based on core genome alignment. The analysis we performed (Figures 1 and 2) clearly showed that all the strains cited above, with, additionally, P. lactis type strain WS 4992 and P. paralactis type strain WS 4672, are all included in the same paraphyletic clade, so they should be classified as the same species. The availability of a wider genomic dataset of Pseudomonas strains from dairy environment, with similar or different phenotypes, would help in the definition and resolution of this classification.
Predicted genes were assigned to the clusters of orthologous groups (COG) classification (Supplementary Figure S1). Despite the group of general function prediction, which was the largest, the highest count for both strains was related to amino acid transport and metabolism and to inorganic ion transport and metabolism followed by signal transduction mechanisms. Figure 3 shows the pairwise heatmap which visualizes the overlapping cluster numbers for nine Pseudomonas isolates. ITEM 17298 and ITEM 17299 shared the highest number of orthologue cluster (6055), followed by Pseudomonas spp. Ps_40 and P. lactis WS 4992 (5863) and Pseudomonas spp. Ps_77 and Ps_22 (5861). The species form 7004 clusters, 2882 orthologous clusters (at least contains two species), and 4122 single-copy gene clusters. All nine isolates shared 4171 clusters with a total of 37,689 proteins counted.
phylogenetic tree based on core genome alignment. The analysis we performed (Figures 1 and 2) clearly showed that all the strains cited above, with, additionally, P. lactis type strain WS 4992 and P. paralactis type strain WS 4672, are all included in the same paraphyletic clade, so they should be classified as the same species. The availability of a wider genomic dataset of Pseudomonas strains from dairy environment, with similar or different phenotypes, would help in the definition and resolution of this classification.
Predicted genes were assigned to the clusters of orthologous groups (COG) classification (Supplementary Figure S1). Despite the group of general function prediction, which was the largest, the highest count for both strains was related to amino acid transport and metabolism and to inorganic ion transport and metabolism followed by signal transduction mechanisms. Figure 3 shows the pairwise heatmap which visualizes the overlapping cluster numbers for nine Pseudomonas isolates. ITEM 17298 and ITEM 17299 shared the highest number of orthologue cluster (6,055), followed by Pseudomonas spp. Ps_40 and P. lactis WS 4992 (5,863) and Pseudomonas spp. Ps_77 and Ps_22 (5,861). The species form 7,004 clusters, 2,882 orthologous clusters (at least contains two species), and 4,122 single-copy gene clusters. All nine isolates shared 4,171 clusters with a total of 37,689 proteins counted. Belonging to the clusters unique to ITEM 17298, OrthoVenn2 analysis indicated two Phagerelated minor tail proteins, and, not shared with the other ITEM strains or with P. lactis or P. paralactis, several hypothetical proteins with no GO (Gene Ontology) associated, one protein involved in response to arsenic-containing substance, phage related proteins, proteins for the viral genome integration into host DNA and conjugation proteins. The same occurs for ITEM 17295, which has a wide repertoire of viral tail assembly and integration proteins and, unique to all other isolates, two Epoxide hydrolase-like proteins (GIB65_13865 and GIB65_16720).

Biofilm and Motility
As reported in Figure 4, ITEM 17298 produced the highest amount of biofilm under all experimental conditions over the incubation time (Absorbance CV value > 4). In addition, ITEM 17299 Belonging to the clusters unique to ITEM 17298, OrthoVenn2 analysis indicated two Phage-related minor tail proteins, and, not shared with the other ITEM strains or with P. lactis or P. paralactis, several hypothetical proteins with no GO (Gene Ontology) associated, one protein involved in response to arsenic-containing substance, phage related proteins, proteins for the viral genome integration into host DNA and conjugation proteins. The same occurs for ITEM 17295, which has a wide repertoire of viral tail assembly and integration proteins and, unique to all other isolates, two Epoxide hydrolase-like proteins (GIB65_13865 and GIB65_16720).

Biofilm and Motility
As reported in Figure 4, ITEM 17298 produced the highest amount of biofilm under all experimental conditions over the incubation time (Absorbance CV value > 4). In addition, ITEM 17299 was a strong biofilm producer; however, its amount was significantly (p < 0.05) lower by 30 and 43%, on average, than ITEM 17298 both at 15 and 30 • C, respectively. In accordance with our previous data [22], both strains doubled biofilm biomass at the lowest temperature. By contrast, ITEM 17295 showed no detectable biofilm biomass at 30 • C, whilst absorbance values lesser than 0.3 were registered at 15 • C over the entire period of incubation. Statistical analyses indicated that biofilm production significantly depended on the strains and temperature, whilst the incubation time was irrelevant; in addition, the interaction between time and temperature showed no significant (p > 0.005) effect on biofilm formation by each tested strain. Therefore, the low temperature was more effective than that of the optimal growth to promote biofilm production only for ITEM 17298 and 17299.
As widely reported, biofilms are organized structures surrounded by extracellular matrix in which the cells are embedded and protected from the environmental stresses. Recently, biofilms and QS have been suggested as important factors in the deterioration of dairy products, opening new fields of investigation and challenges to counteract spoilage [21,24,51,52]. In particular, the production of hydrolytic enzymes (extracellular proteases, lipases, glycosidases, and phospholipases) have been demonstrated to be enhanced in QS regulated-biofilms by several bacterial spoilers, including P. fluorescens [51][52][53]; thus, innovative antibacterial strategies based on the application of QS inhibitors have been considered promising to prevent food spoilage.
Among factors affecting the biofilm phenotype, a positive correlation between low temperatures and biofilm production was previously demonstrated in foodborne P. fluorescens [54,55]. In accordance with these authors, our results highlighted that low incubation temperatures were more important than incubation time for biofilm formation, putatively favoring spread and persistence in the dairy sector. on average, than ITEM 17298 both at 15 and 30 °C, respectively. In accordance with our previous data [22], both strains doubled biofilm biomass at the lowest temperature. By contrast, ITEM 17295 showed no detectable biofilm biomass at 30 °C, whilst absorbance values lesser than 0.3 were registered at 15 °C over the entire period of incubation. Statistical analyses indicated that biofilm production significantly depended on the strains and temperature, whilst the incubation time was irrelevant; in addition, the interaction between time and temperature showed no significant (p > 0.005) effect on biofilm formation by each tested strain. Therefore, the low temperature was more effective than that of the optimal growth to promote biofilm production only for ITEM 17298 and 17299.
As widely reported, biofilms are organized structures surrounded by extracellular matrix in which the cells are embedded and protected from the environmental stresses. Recently, biofilms and QS have been suggested as important factors in the deterioration of dairy products, opening new fields of investigation and challenges to counteract spoilage [21,24,51,52]. In particular, the production of hydrolytic enzymes (extracellular proteases, lipases, glycosidases, and phospholipases) have been demonstrated to be enhanced in QS regulated-biofilms by several bacterial spoilers, including P. fluorescens [51][52][53]; thus, innovative antibacterial strategies based on the application of QS inhibitors have been considered promising to prevent food spoilage.
Among factors affecting the biofilm phenotype, a positive correlation between low temperatures and biofilm production was previously demonstrated in foodborne P. fluorescens [54,55]. In accordance with these authors, our results highlighted that low incubation temperatures were more important than incubation time for biofilm formation, putatively favoring spread and persistence in the dairy sector. In addition to biofilm formation, P. lactis motility (swimming, swarming and twitching) was also investigated. Indeed, cell surface characteristics (occurrence of chemotaxis systems or hydrophobicity) and motility are crucial in nutrient acquisition, cell-cell and cell-surface interactions, stress tolerance, sites colonization and the initial attachment during biofilm formation [56]. As expected, patterns of motility were different depending on agar percentage used (Supplementary Figures S2 and S3); low agar percentage (0.3%) favored swimming motility (Supplementary Figure In addition to biofilm formation, P. lactis motility (swimming, swarming and twitching) was also investigated. Indeed, cell surface characteristics (occurrence of chemotaxis systems or hydrophobicity) and motility are crucial in nutrient acquisition, cell-cell and cell-surface interactions, stress tolerance, sites colonization and the initial attachment during biofilm formation [56]. As expected, patterns of motility were different depending on agar percentage used ( Supplementary Figures S2 and S3); low agar percentage (0.3%) favored swimming motility (Supplementary Figure S2), a flagella-mediated movement in liquid media [57]; swimming motility was registered for each strain and mainly affected under higher temperature of incubation.
ANOVA analysis revealed that there was a statistically significant interaction between the effects of strain, time, and temperature on swarming (F(4, 36) = 3.719, p = 0.0012). In particular, swarming of all strains significantly increased as a time-dependent manner at the optimal growth temperature, whereas, at low temperatures, they swarmed later but to a lesser extent; ITEM 17295 showed smaller time-dependent increments at both incubation temperatures (Supplementary Figure S3). No appearance of elongated swarm edge was recorded over the time. Twitching motility, requiring type IV pili, was instead registered for ITEM 17298 at both temperature and significantly (χ 2 (5) = 52.89, p < 0.0001; Kruskal-Wallis rank sum test) increased depending on time of incubation reaching the highest value at 72 h; at the end of incubation, ITEM 17299 exhibited a twitching only at 15 • C but with an average value significantly lesser than those of ITEM 17299 (Supplementary Figure S3).
Based on these results, a genomic determinant related to the above reported phenotypic traits was identified. Table 2 listed most biofilm-related genes in all analyzed P. lactis strains; genes were catalogued with respect to their role in the consequential steps occurring during biofilm formation, including control mechanisms of motility apparatus. Even though results from biofilm and motility assay demonstrated a specific behavior for each strain, comparative analysis highlighted that all strains shared the majority of the biofilm-related genes. This is expected, considering the genetic similarity of the strains described in previous paragraphs. Differences occurred in a just few genes: among the most relevant, ITEM 17295 missed the pel operon, orthologue to the previously described in P. aeruginosa PA3058-PA3064. The operon codes for pellicle biofilm biosynthesis proteins, and it is present in ITEM 17298 (PROKKA_02829-PROKKA_02835) and ITEM 17299 (GIB64_23225-GIB64_23255) ( Table 2). In P. aeruginosa, pel operon is involved in the formation on the extracellular matrix that surround the cells in a mature biofilm [58]. In addition, it modulates colony morphology giving rise to wrinkled red colonies on Congo red agar; these traits were confirmed by pel deletion in mutant strains of P. aeruginosa appearing as orange and invariably flat and smooth colonies [59]. Differences in color were correlated to the production or the absence of a Congo-red-binding component of the biofilm matrix. Although both ITEM 17298 and ITEM 17299 genomes harbor pel operon, a wrinkled red morphology appeared only in ITEM 17298 when incubated at 15 • C (Supplementary Figure S4), suggesting that pel operon genes are differently regulated in these strains and putatively expressed only in ITEM 17298 in this latter condition. As expected, orange and flat colonies were instead obtained for ITEM 17295 (Supplementary Figure S4).
Among the genes shared by all three ITEM strains, we found the alginate and psl operons; alginate operon showed a genetic organization similar but not identical to that of P. aeruginosa, in which arn operon is not contiguous to the algD operon, as reported by Quintieri et al. [24].
Although their role in cell physiology needs to be investigated further, some TAs are involved in the switch from the planktonic to the biofilm lifestyle, as well as virulence and multidrug resistance in pathogenic pseudomonads [61][62][63]; in particular, MqsR/MqsA plays an important regulatory role in the persistence and biofilm formation by P. putida and E. coli [61]. Few differences were also retrieved in the content of biofilm related genes in the genomes of other P. lactis, such as P. lactis SS101 and P. lactis DSM 29,167 (Fanelli, unpublished results), highlighting certain genetic variability among these strains.
As reported in Supplementary Table S3, no differences were found among strains in the list of genes involved in flagella biosynthesis and twitching motility. This list was obtained by homology search towards P. aeruginosa, possessing a single polar flagellum that plays a variety of roles in the virulence in addition to its central role in swimming motility. In all strains, we retrieved quite the same genetic content and organization of P. aeruginosa, with genes encoding the flagellar basal body MS ring and motor switch complex, and genes coding for flagella assembly and the sigma factor regulating this process fliA and rpoN [64,65], as well as the fleQSR locus coding for the major flagella gene regulator in P. aeruginosa and the two-component sensor system fleSR involved in flagellin synthesis. By contrast, the flagellar glycosyltransferase FgtA, involved in the P. aeruginosa flagellin glycosylation that in turn exacerbates its virulence [66], was absent in all P. lactis. Twitching motility genes, such as pilT, pilG, pilH, pilI, and pilJ were also identified.

Microcolonies formation and maturation
QS-signals Cell-cell signals type VI secretion system T6SS GIB65_23600 GIB65_23520

Pigment Production
Pseudomonas spp. strains are able to release multifunctional pigments. Although in dairy products they are primarily implicated in spoilage [17], in the bacterial cells they exhibit specific activities including regulatory functions [67] and protective roles under stress response [68]. In addition, for this phenotypic trait, each strain showed a specific behavior (Supplementary Table S4); no pigments were released on the selective media by ITEM 17295, whilst pigmented colonies were found for the remaining strains (Supplementary Table S4). The list of identified genes referred to pyoverdines, pyomelanin, and indigo-derivatives synthesis is shown in Supplementary Table S5. As concern pyoverdines, genes are shared and contiguous, and similarly organized in all the ITEM strains, with the only exception of pvdS, coding for a sigma factor which is required for the expression of pyoverdine biosynthetic genes [69], in contrast to ITEM 17298 and 17299; indeed, this latter was predicted as a pseudogene in ITEM 17295. This could explain in ITEM 17295 the absence of pyoverdine on the selective medium. In previous studies, it has been revealed that pvdS encodes an alternative sigma factor which directs RNA polymerase to the promoters of pyoverdine biosynthesis genes; mutants of pvdS gene in P. aeruginosa, indeed, failed to make detectable pyoverdine when grown on King's B agar or in liquid medium [70].
Among the relevant differences we identified among our three P. lactis strains are the absence in ITEM 17295 of the trpC gene coding for the tryptophan biosynthesis protein (PROKKA_00895, GIB64_11525) and the tryptophan synthase alpha and beta chain (PROKKA_00904-00905; GIB64_11475-11480), which were demonstrated to be involved in the pathways leading to the formation of indigo-derivatives.
The locus comprising these genes is homologue to the c4_BAR (Contig 4 Blue Accessory Region) that Andreani et al. [71] identified as only present in Blue pigmenting pseudomonads strains. This region was described as containing 16 genes (16 Kb) including those coding for trp accessory genes: trpD, trpF, trpA, and trpC.
A deeper investigation into our pigmenting strains ITEM 17298 and ITEM 17299 confirmed the presence of this locus in both strains, a locus that is absent in ITEM 17295 that is not pigmenting; the genomic context of the flanking regions (data not shown) reveals the occurrence of genetic elements which suggest the presence of integrative conjugative plasmid. Indeed, we identified several integrating conjugative elements, as well as pilus formation and assembly proteins. Furthermore, we detected genes coding for one integrase as well as a locus comprising several t-RNA coding genes, which are a hotspot for genomic island integration.

Protease and Lipase
The tested strains exhibited hydrolytic activity against milk proteins showing a time-dependent manner increase of proteolytic activity as the following order: ITEM 17295 > ITEM 17299 > ITEM 17298 ( Figure 5); these results agreed with our previous results describing ITEM 17295 (called PS37) as strong proteolytic strain in milk samples [12]. The incubation at the temperature lower than that of optimal growth caused a significant decrease in halo diameters, much greater for ITEM 17298 ( Figure 5).
The genomic analysis identified several protease genes (aprA, prsDE, prtAB) in all tested strains (Supplementary Table S6). These genes are positioned in a QS-regulated operon previously associated with milk spoilage by P. fluorescens [72]. In particular, aprA, encoding a heat stable caseinolytic alkaline metalloprotease [73], is widespread in Pseudomonas, especially in the psychrotrophic strains [74] to the extent that it is recommended in the search of milk spoilage pseudomonads instead of rpoB [75]. In previous works, AprA isolated from several Pseudomonas spp. was described as 98% similar to the peptidase AprX secreted by one strain of P. fluorescens [72]. AprA and AprX generally exhibit activity in a large range of temperatures (0-55 • C) with optimal activity between 37 and 47 • C [76]. In addition to the improved growth at 30 • C, this could explain the diameter increase registered in all strains under this experimental condition ( Figure 5).  Although all the herein assayed strains shared the spoilage aprA-lipA operon, the genome of the most proteolytic strain, ITEM 17295, bears a second gene coding for a metallopeptidase AprA (family M10) GIB65_09350 (Supplementary Table S7). This genomic locus comprising aprA also included the gene coding for chemotaxis protein CheB, protein-glutamate O-methyltransferase CheR, several response regulators, transporters, and several fimbrial proteins, and is flanked by two recombinase/integrase. A similar locus was retrieved also in the other strains, leading to the assumption that this locus is derived from an integration event, but with relevant differences in the gene content, with aprA and other elements exclusive of ITEM 17295.
Supplementary Table S7 reports a list of the peptidases/transpeptidases identified in each P. lactis genome; beyond proteolytic activity, several enzymes were involved in biological processes, such as stress response and signal transduction. In accordance with the high similarity among strains, just a few differences were revealed.
The examination of P. lactis proteolytic halos under UV lamp revealed fluorescence around ITEM 17928 colonies, putatively due to pyoverdine release ( Figure 5, panel B, lower side) and most marked with regard to the growth at low temperature. These results agreed with the results reported above for pigment release. Indeed, only ITEM 17298 produced pyoverdine although in synthetic media the synthesis mainly occurred at the low temperature. As widely reported in literature, pyoverdines can be considered as signaling molecules mediating specific mechanisms as the spoilage induced by proteases and lipases activities [77][78][79]. Preservation techniques, indeed, modulating pyoverdines synthesis to counteract pseudomonads growth, spoilage, or pathogenesis are a current challenge [80].
Lipase activity was apparently registered after 72 h of incubation only for the strains ITEM 17295 and ITEM 17299 colonies displaying under UV lamp a brilliant glow pink-red color and putatively associated with intracellular release of fatty acids (Supplementary Figure S5). At a low temperature, only ITEM 17295 proved to grow on the Rhodamine agar with fluorescent smaller colonies (Supplementary Figure S5). However, evident blue fluorescent halos were registered around colonies of all tested strains grown at 30 • C putatively associated with the pyoverdine release. As reported above, two lipases were identified in the spoilage operon (Supplementary Table S6) in all assayed strains. The activity of lipase included in aprX-lipA operon are regulated by EnvZ-Ompr [60], herein identified and annotated as PROKKA_02877-PROKKA_02878 for ITEM 17298, GIB65_26250-GIB65_26255 for ITEM 17295, and GIB64_23460-GIB64_23465 for ITEM 17299. Although each strain exhibited different lipolytic activity, no differences were reported in the list of lipolytic enzymes by genomic analysis (Supplementary Table S8); this suggested that, in each P. lactis strain, different regulatory mechanisms could activate the expression of these genes.

Conclusions
In this study, we analyzed three cheese-borne spoilage Pseudomonas in order to identify the genomic determinants responsible for their phenotypic traits. Based on the results of validated genomic methods, these strains grouped in the paraphyletic clade of Pseudomonas lactis, and thus they were classified as so. Despite their high genetic similarity, phenotypic analysis showed distinctive traits in relation to biofilm formation, pigmentation, proteolytic, and lipolytic activities. Some of these characteristics (biofilm and pigmentation) are promoted at lower temperatures of incubation, suggesting their possible involvement in the spreading and persistence of these strains in the dairy environment. Several differences in the identified genetic elements were consistent with the phenotypic diversity among these strains. In this work, psychrotrophic pseudomonad strains spoiling mozzarella cheese were reclassified as P. lactis, previously associated with the spoilage of the refrigerated raw milk. The availability of two novel P. lactis genomic sequences and the analysis here performed will improve the knowledge about the genetic elements involved in their spread, persistence, and spoilage, thus paving the way towards novel preservative strategies.

Supplementary Materials:
The following are available online at http://www.mdpi.com/2076-2607/8/8/1208/s1. Figure S1. COGs functional classification of genes present in P. lactis ITEM 17295, ITEM 17298, and ITEM 17299 genomes. Figure S2. Swimming motility of P. lactis ITEM 17295, ITEM 17298 and ITEM 17299, in LB performed at 15 • C and 30 • C for 24, 48 and 72 h. Bars represent the mean diameters (±standard deviation, n = 3) of corresponding motility zones. Similar values (p > 0.05) are represented by the same subscript letters according HSD Tukey post-hoc test. Figure S3. Swarming and twitching motility of P. lactis ITEM 17295, ITEM 17298 and ITEM 17299, registered in LB and M8 medium, respectively for 24, 48, and 72 h at 15 • C and 30 • C. Bars represent the mean diameter (±standard deviation, n = 3) of corresponding motility zones. Twitching motility values were determined by measuring the absorbance of Crystal violet (CV) at 570 nm. Bars represent the average ± standard deviation (n = 3). Similar values (p > 0.05) of swarming and twitching are represented by the same subscript letters according HSD Tukey post-hoc test and Wilcoxon rank sum test and p-value adjustment Bonferroni's method, respectively. Figure     Acknowledgments: It is a pleasure to acknowledge Pasquale Del Vecchio for the assistance with journal publication policy and copyright.

Conflicts of Interest:
The authors declare no conflict of interest.