Description of a Novel Fish Pathogen, Plesiomonas shigelloides subsp. Oncorhynchi, Isolated From Rainbow Trout (Oncorhynchus mykiss): First Genome Analysis and Comparative Genomics

1. Introduction

Plesiomonas shigelloides, formerly known as Pseudomonas shigelloides, has a debatable taxonomic history and is currently classified in the family Enterobacteriaceae as the only member of the genus Plesiomonas [1]. This bacterium naturally inhabits freshwater and marine environments as well as animals dwelling those ecosystems including fish, amphibia and insects [2]. P. shigelloides occasionally causes systemic infections, such as acute diarrhea and extraintestinal diseases, mainly in immunocompromised individuals. The natural habitat of P. shigelloides is associated with its transmission to human body since this agent enters the body either through contact with contaminated water or seafood. Asymptomatic carriage of the agent is common among freshwater and marine species such as oysters, shrimp, trout, sea bass or sea bream.

A large number of freshwater plesiomonads strains have been commonly reported from rainbow trout [3,4], carp [5], and tilapia [6,7,8,9]. Studies of the intestinal tracts of several freshwater fish suggest that the genus Plesiomonas is one of the most common species composing the bacterial microbiota of these vertebrates (58, 59). Plesiomonas has not yet been associated with post-natural disaster infections (45). However, indirect evidence suggests that plesiomonads could be involved in illnesses after major natural aquatic disasters [9]. High mortality rate in trout due to P. shigelloides alone [10] and along with other bacteria like Flavobacterium spp. and Aeromonas hydrophila were also reported [11]. This pathogen was also identified as one of the main pathogens in cultured sturgeons [12]; it was also isolated from fishes with clinical symptoms during mass mortality of Ctenopharyngodo nidellus [5] and Orechromis nilotica [13] providing more evidence that P. shigelloides might be highly pathogenic for these cultured fishes. Besides cultured species, 100% mortality rate was recorded in cichlid ornamental fish due to infection caused by P. shigelloides [14].

While the agent is assumed to be a natural member of microbiota or responsible for disease, especially found along with other fish pathogens, the information about the genomic characteristics of P. shigelloides is very limited. Almost all present data strongly suggest that the genus and species designation P. shigelloides is composed of a collection of homogeneous bacteria at the phenotypic and molecular levels [9]. A population study of a diverse collection of 77 P. shigelloides strains from different geographic as well as environmental settings indicated a monophyletic clade nested within the Enterobacteriaceae [15,16]. Complete genome information has noted P. shigelloides strains isolated from humans or terrestrial animals but not in fish [9,17,18]. The lack of information causes a gap in the relation of the agent isolated from fish to humans in zoonotic potency, transmission routes and virulence. The Plesiomonas genus strains, has commonly been reported from aquaculture and fish by 16S rRNA or 23S rRNA gene sequencing or housekeeping genes such as hugA [19,20]. Successful bacterial identification reports show P. shigelloides prevalence in fish, but the paucity of genomic information limits its exact role of it by means of virulence characteristics and antimicrobial resistance gene transfer.

To the best of our knowledge, there was only one report on the genome analysis and deep molecular characterization of P. shigelloides isolates recovered from rainbow trout [21]. We aimed to characterize P. shigelloides strains isolated from fish by genome analysis and annotate the first genome of strain P. shigelloides V-78 as a novel subspecies. In addition, virulence and antimicrobial resistance genes detected in the genome of strains will help understand the pathogenicity of this agent for fish, notably found in other fish pathogens.

2. Material and Method

2.1. Bacteria Isolation and Biochemical Characterization

Isolates were collected from rainbow trout (Oncorhynchus mykiss) farms using conventional flow-through systems with non-disinfected surface water and spring water in the Central Anatolia Aegean region of Turkey. Besides rainbow trout isolates, another group of strains were also collected from ornamental freshwater fish which were transferred to the fish disease laboratory between 2018-2021. The bacterial isolates were collected during a fish health surveillance program from internal organs of diseased fish (showing darkening in color, exophthalmia and ulceration on skin and fin) of different weights, such as 0,5-1 g and 200 g rainbow trout, and from species of freshwater ornamental fish. Two strains, designed as 63 and V-78 isolated from the samples collected from Aegean region in 2018 and Central Anatolia in 2013, respectively, were selected for comprehensive genome analyses.

Isolates were biochemically characterized by Gram-staining, motility, oxidase and catalase reaction and glucose fermentation of O/F medium (Merck 103865). Hemolytic characteristic was tested on 5% sheep blood-added agar, and sensitivity to Vibriostatic agent were also tested (Oxoid DD0014 and DD0015). Growth in the presence of NaCl was tested in tryptic soy broth (TSB, Merck 105459) supplemented with 0-3% NaCl (w/v). Optimum growth temperature range was tested in the range of 4-50°C on TSB medium. In addition, the ability to grow on Thiosulfate Citrate Bile Sucrose (TCBS) agar (Merck, 103854) and to hydrolyse gelatin were determined by following published methods [22].

2.2. DNA Extraction, PCR and Sequence Analysis

The genomic DNA of the isolates was extracted using a QIAamp DNA mini kit (Qiagen, Hilden, Germany) according to the manufacturer’s instructions. The amount and purity of DNA in each sample were measured at 260 nm and 260/280 nm wavelengths with a spectrophotometer (Multiskan Go, Thermo).

The isolates were identified by MLSA analysis using 16S rRNA, atpA (ATP synthase α subunit), gyrB (DNA gyrase β subunit), pryH (Uridylate kinase) and recA (Recombinase A) gene sequences [22,23]. Polymerase chain reaction (PCR) was performed as described in previous works [22]. All PCR products were sequenced by Macrogen Korea (Republic of Korea). The DNA fragments used for the analysis are 1400 bp for 16S rRNA, 570 bp for gyrB, 501 bp for pyrH, 462 for recA and 489 for atpA.

2.3. Genome Sequencing, Assembly and Identification

For comprehensive genome-based analyses, the whole-genome sequence of V-78 was obtained as described previously [17]. The genome of strain V-78 was sequenced with the Oxford Nanopore MinION system (Oxford Nanopore Technologies, Oxford, UK). The high-quality reads of the genome were assembled into contigs by de novo assembly using the Canu v. 1.7.1 for strain V-78. The draft genome sequence data was submitted to the GenBank. The annotation of the contigs longer than 1000 bp was achieved by the NCBI Prokaryotic Genome Annotation Pipeline (PGAP) [24]. In addition, the genome of strain V-78 was annotated on the Rapid Annotations Using Subsystems Technology (RAST) server (https://rast.nmpdr.org/) [25] using RASTtk pipeline [26]. The type strain genomes server pipeline (TYGS, https://tygs.dsmz. de/) was employed to perform phylogenomic analysis and determine the digital DNA-DNA hybridization (dDDH) values between the genome of strain V-78 and those of the type strains deposited in the DSMZ database.

2.4. Comparative Genome Analysis

The publicly available genomes of P. shigelloides in the NCBI GenBank database were retrieved for comparative genome analysis. From 37 genomes, a total of 19 P. shigelloides genomes whose isolation source is known were selected for further analyses. The Codon Tree method selects single-copy Bacterial and Viral Bioinformatic resource center (BV-BRC) programs and analyzes aligned proteins and coding DNA from single-copy genes using the program RAxML (https://www.bv-brc.org/app/PhylogeneticTree [27]. For comprehensive genomic characterization, secondary metabolite biosynthetic gene clusters and antibiotic resistance models encoded by the genome of strain V-78 were determined by using antiSMASH version 6.1.1 (https://antismash.secondarymetabolites.org/#!/start) [28] and Antibiotic Resistant Target Seeker version 2 (ARTS) (https://arts.ziemertlab.com/) [29], respectively. To determine the virulence genes of the strains, VFDB (virulence factor of pathogenic bacteria) online database [http://www.mgc.ac.cn/cgi-bin/VFs/v5/main.cgi] was used. The GenBank files of P. shigelloides genomes were subjected to BLASTn against the VFDB with default settings. In addition, the genomic islands, the clusters of genes of probable horizontal transfer origin, were predicted by IslandCompare software available at https://islandcompare.ca/ [30]. The compilation of taxonomic, functional and ecological features of strain V-78 was acquired by Protologger web server available at http://protologger.de/ [31].

2.5. Antimicrobial Susceptibility and Resistance Gene Characteristics

The isolates were tested for their susceptibility to the antimicrobials by an agar disk diffusion method as was essentially described in the Clinical and Laboratory Standards Institute (CLSI) guideline M42-A [32], using Mueller-Hinton Agar (MHA, BBL-Becton Dickinson). The antibacterial susceptibility patterns of isolates were performed using disks containing the antibacterial agents: Trimethoprim/Sulfamethoxazole (SXT, 25 µg), Amoxicillin (AML, 25 µg), Doxycycline (DO, 30 µg), Tetracycline (TE, 30 µg), Oxytetracycline (OT, 30 µg), Oxolinic Acid (OA, 2 µg), Enrofloxacin (ENR, 5 µg), Ciprofloxacin (CIP, 5 µg), Flumequine (UB, 30 µg), Cefalexin (CN, 10 µg), Ampicillin (AMP, 10 µg), Erythromycin (E, 15 µg), Lincomycin (MY, 15 µg) and Florfenicol (FFC, 30 µg). All disks were obtained from Oxoid Ltd. (Basingstoke, Hampshire, England). The bacterial isolates were suspended in sterile 0.85% saline to obtain a 0.5 McFarland (bioMerieux S.A.) turbidity and inoculated on the test media. The plates were incubated for 24–28 h at 28°C. The control strain E. coli ATCC 25922 was incubated at 28°C for 24-28h [32].

3. Results

3.1. Biochemical Characterization

The strains were Gram-negative, short bacilli, motile, oxidase and catalase positive, glucose fermentative, non-hemolytic on sheep blood and sensitive to Vibriostatic agent (O/129 150µg). The tolerance tests for growth in the presence of NaCl and in a range of temperatures showed that the strains could grow in the presence of up to 1.5% NaCl and at the temperature range of 4-45°C. Unlike Vibrios, the isolates V-63 and V-78 did not grow on TCBS medium and could not hydrolyze gelatin.

3.2. PCR and Sequence Analysis

The strains were identified as Plesiomonas spp. by the 16S rRNA gene sequence analysis. For a better identification the following genes were sequenced: gyrB, pryH, recA and atpA, suggested in PubMLST database for Vibrio sp. These genes were analyzed by Blast in the NCBI database and the closest match were P. shigelloides. The sequences were accessed to the PubMLST database and sequence types were attained 133, 103, 127, 93, for V-63 and 136, 104, 128, 95, for V-78, by gyrB, pryH, recA and atpA, respectively. The isolates were identified, ST: 177 and 178 for V-63 and V-78, respectively, in the PubMLST database due to the similarity of the Vibrioneaceae genus. The strains are closely related to P. shigelloides reference strain, sharing at least 99% similarity for the mentioned genes.

3.3. Phylogenomics

Genome-based identification of V-78 strain was done with the TYGS web server and identified as P. shigelloides (Figure 1). The dDDH value for strain V-78 compared with its most closely related type strain P. shigelloides NCTC 10360^T was calculated as 74.2%.

At the time of writing the manuscript (February 2023), a total of 37 Plesiomonas shigelloides genomes were available in the GenBank database. However, six of these genomes were assembled from metagenomics data, 2 genomes lacked information for the isolate source, and three genomes did not have information for the country of origin. Therefore, it was determined that eight genomes were obtained from the isolates from various fish species (Table 1). The codon tree was built by including 19 genomes. The codon tree reveals phylogenomic heterogeneity within the genus Plesiomonas inferred from the alignment of coding DNA from single-copy genes and amino acid sequences The strain V-78 were located in a distinct branch of the tree (Figure 2).

The genome of strain V-78 harbors 13 gene sequences coding for 16S rRNA whose pairwise identity level with that of P. shigelloides NCTC 10360^T ranged from 99.38% to 99.93%. These rRNA gene sequences differed also by length ranging from 1534 to 1541 and intra-strain similarity levels of 99.28-100%.

The genomes of closely related species were retrieved form the Genome Taxonomy Database (GTDB) by the Protologger to perform pairwise analyses for the average nucleotide identity (ANI) and percentage of conserved proteins (POCP). The ANI values were calculated by FastANI (v1.2) algorithm [33] while the POCP values were acquired by custom Python scripts available at the GitHub repository by the Protologger (Table S1). The genome of strain V-78 shared 97.1% ANI value with the type strain of P. shigelloides while that value with the other closely related taxa were found to be lower than 80% implying that the genus exhibit considerable distance to the other genera (Figure S1). The highest POCP value of 68.1% was calculated between the genome of strain V-78 and that of the type strain of P. shigelloides as expected. The other POCP values were below the threshold of 50% for genus demarcation [34] proving the distinctness of the genus Plesiomonas (Table 2).

3.4. Functional and Ecological Analyses

Functional and ecological features of strain V-78 was predicted by the Protologger software. Of 5454 coding sequences identified, 235 genes encode for transporter proteins, while 50 genes are responsible for secretion. The number of unique enzymes encoded by the genome of V-78 was determined as 920. Arbutin, glucose, salicin, sucrose and trehalose were predicted to be utilized by the strain as carbon sources. In addition, strain V-78 could produce acetate from acetyl-CoA (EC:2.3.1.8, 2.7.2.1) and propionate from propanoyl-CoA (EC:2.3.1.8, 2.7.2.1) as predicted by the Protologger. The flagellar proteins identified in the genome are FlhA, FlhB, FlgB, FlgC, FlgD, FlgE, FlgF, FlgG, FlgJ, FlgK, FlgL, FliD, FliF, FliG, FliM, MotA and MotB. The strain could produce L-cysteine and acetate by utilizing sulfide and L-serine (EC:2.3.1.30, 2.5.1.47). The genome also has pathways for L-glutamate production from ammonia via L-glutamine (EC:6.3.1.2, 1.4.1.-) and biosynthetic genes for biotin, folate and riboflavin synthesis. A total of 157 carbohydrate-active enzymes (CAZymes), mainly in glycoside hydrolase (GH) families, were identified within the genome of strain V-78. Glycoside transferase (GT), carbohydrate esterase (CE), polysaccharide lyase (PL) and carbohydrate-binding module (CBM) families were also predicted in the genome.

For an ecological analysis, both the 16S rRNA gene and whole genome sequences were submitted to the Protologger [31]. The genome of strain V-78 was compared with a database of 49,094 metagenome-assembled genomes (MAGs) collected from various habitats, as described by Hitch et al. [31]. Briefly, the comparison was performed using MASH [35], and the results were filtered by using a distance threshold of <0.05. In addition, the 16S rRNA gene was employed to perform a comparative analysis involving operational taxonomic units (OTUs) generated from 19,000 amplicon datasets from 19 different habitats obtained from the IMNGS database [36]. Only one MAG (MAG-ID: HeQ_2017__SZAXPI029476-50__bin_9.fa) collected from an adult human stool [37] clustered with the genome of strain V-78. Habitat distribution and preference analysis based on comparing 16S rRNA gene amplicons revealed that strain V-78 was detected in all 19 environments in the IMNGS database with a total of 38,163,501 OTUs (Figure 3). The highest detection ratio was observed within the rhizosphere amplicons, while the mean relative abundance in human gut amplicons was the highest compared to the other environments (Table 3).

3.5. Secondary Metabolites

Secondary metabolite biosynthesis gene clusters encoded by the genome of strain V-78 were identified by antiSMASH server. The strain has two gene clusters coding for a betalactone and a thiopeptide. The betalactone gene cluster shared 2% similarity with the gene cluster of gausemycin A/B, a lipoglycopeptide antibiotic, while the thiopeptide gene cluster showed 33% similarity to menaquinone biosynthetic gene cluster. In addition, the RAST annotation revealed that the genome of strain V-78 encodes genes for anthranilate phosphoribosyltransferase, phosphoribosylanthranilate isomerase, tryptophan synthase alpha and beta chains, which are responsible for synthesis of auxin-type plant hormones.

3.6. Antimicrobial Susceptibility and Resistance Gene Characteristics

Strains were fully resistant to Lincomycin, while full susceptibility were determined against Sulfamethoxazole+trimetprim and Ciprofloxacin. The quality control strains, A. salmonicida subps. salmonicida and E. coli fell the suggested values in CLSI (Table 4). However, cut-off values are not presented in EUCAST [38] for Amoxicillin, Doxycycline, Tetracycline, Oxytetracycline, Oxolinic Acid, Enrofloxacin, Flumequine, Cefalexin, Ampicillin, Erythromycin, Lincomycin, and Florfenicol so the strains could not be evaluated to be susceptible or resistant (Table 4).

3.7. Antimicrobial Resistance Genes and Virulence

A total of 24 antimicrobial resistance models and targets were determined by the ARTS server across 19 genomes. The resistance models detected include ATP-binding cassette (ABC) antibiotic efflux pump, CARB-PSE beta-lactamases (class a), chloramphenicol acetyltransferase (CAT), ABC transporter, membrane fusion proteins of the MexEF-OprN and MexGHI-OpmD multidrug efflux complexes, major facilitator superfamily (MFS) antibiotic efflux pump and resistance-nodulation-cell division (RND) antibiotic efflux pump. Of 24 resistance models and targets, class A beta-lactamase (RF0053) and tetracycline resistance MFS efflux pump were encoded only in the genome of strain V-78. The former is responsible for the resistance through hydrolysis of beta-lactam antibiotics such as penicillins, cephalosporins and carbapenems, while the latter enables the bacteria to pump out tetracycline or tetracycline derivatives selectively (Table S2).

Identifying the genomic islands encoding mainly virulence factors, antimicrobial resistance genes and other adaptation elements is a useful strategy for adopting a population-based approach to genomic epidemiology and characterization [30]. The IslandCompare tool was used to identify and compare the genomic islands residing in the genomes of P. shigelloides strains. All strains differed in terms of size, number, location and identity of the genomic islands encoded in their genomes (Figure 4).

When the virulence gene of the P. shigelloides strains was examined, a total of 83 virulence genes were detected in the various P. shigelloides genomes obtained from the GenBank database and our isolates. Into the 83 virulence genes, only 20 virulence genes were commonly found in all strains encoding motility, immune modulation, T2SS (Type II secretion system), adherence, immunogenic lipoprotein A, stress survival, biofilm formation, stress protein, early escape from the phagosome of macrophages, secreted toxins and enzymes into the extracellular fluid. The virulence genes according to the each analyzed isolate are presented in supplementary Table S3. Different from other P. shigelloides strains, tviB gene encoding to prevent antibody-mediated opsonization, increase resistance to host peroxide and resistance to complement activation by the alternate pathway and complement-mediated lysis was detected only the strain V-78 and MS-17-188 (catfish). The lgtF and flgG genes encoding for immune modulation and motility, respectively, were detected only in the genomes of strain V-78 and strain CAPA003 isolated from Oreochromis niloticus. Functioning to a tripartite multidrug efflux pump essential for resistance to β-lactams (penicillin G and nafcillin), macrolides (erythromycin) and host-derived compounds (peptide LL-37) and progesterone, and essential for growth of gonococci in the lower genital tract of experimentally infected female mice, the mtrD gene was found in V-78, Zfcc051 (isolated from zebra fish) and CAPA003. In addition to fbpC gene, encodes a periplasmic-binding protein-dependent iron transport system necessary for the utilization of iron bound to transferrin or iron chelates was found only in the genome of strain V-78.

4. Discussion

Plesiomonas shigelloides is a unique Gram-negative, polar-flagellated pathogenic bacterium naturally found in freshwater ecosystems, such as rivers, lakes, and surface waters, and marine estuaries in tropical and temperate climates [39]. Unlike other aquatic pathogens, the isolation of P. shigelloides might be enhanced by an incubation temperature of 42–44°C since most aeromonads or other psychrophilic microorganisms could not tolerate this temperature range. The temperature tolerance of bacteria, especially in high conditions (42–44°C), causes high mortality in freshwater ornamental fish species living in climate habitats; besides cold-water fishes such as rainbow trout [40]. In our study, the isolation of P. shigelloides from rainbow trout and freshwater ornamental fish species supports a broad host range and temperature-tolerant characteristics. Although P. shigelloides is grouped in Enterobacteriaceae, the agent is oxidase and catalase positive. It also shares high similarity with the genera Vibrio and Pseudomonas, especially regarding sensitivity for Vibriostatic agent, which we detected in our study. Besides susceptibility to Vibriostatic agent, we identified our isolates by gyrB, pryH, recA and atpA, which were specifically characterized for Vibrio species. Previously, the species were placed chronologically into the genera Pseudomonas, Aeromonas, and Vibrio [41]. Habs and Schubert [42] initially established the genus Plesiomonas in the family Vibrionaceae based on a number of its unique features and its similarity to the genus Aeromonas (“plesio,” neighbor, “monas,” Aeromonas). The genus Plesiomonas has several common characteristics with the members of the family Vibrionaceae, such as being polarly flagellated, facultatively anaerobic, and cytochrome oxidase positive. The similarity between our study and the taxonomical history of P. shigelloides enable the misidentification of this species as Vibrio sp. As commonly reported that P. shigelloides is the only member of Plesiomonas genus, the strains have a genetic homology and close relationship between each other. Unlike the close phylogenetic relationship among P. shigelloides strains, Gu et al. [43] reported that the strains had a significant variability based on isolation source revealed by Random Amplified Polymorphic DNA (RAPD) analysis. In our study, the strains obtained from the GenBank database clustered in a phylogenetic group; however, our strain was separated from all genomes uploaded in the GenBank database. Although the genome of strain V-78 encodes 5454 genes, higher than the genomes of the other strains analysed, it has the lowest number of the core genes among the others, which implies its flexibility for horizontal gene transfer events. The genomic island comparison also approved the versatility of the P. shigelloides genomes since each genome has various genomic islands scattered across the genome.

The members of P. shigelloides naturally occur in freshwater and marine habitats. However, the comprehensive analysis of taxonomic, functional and ecological features revealed that strain V-78 mostly matched with the 16S amplicons obtained from rhizosphere metagenomes followed by freshwater and insect gut metagenomes. Both CAZymes and secondary metabolites, such as plant hormones produced by P. shigelloides strains, may contribute to their adaptation to various environments, from plant rhizosphere to insect gut.

From a taxonomical point of view, strain V-78 is certainly a member of the genus Plesiomonas. However, the overall genome relatedness indices, i.e., dDDH, ANI and POCP, revealed that strain V-78 might be differentiated from the type strain and other strains of P. shigelloides. The dDDH, ANI and POCP values between strain V-78 and the type strain of P. shigelloides were calculated as 74.2%, 97.1% and 68.1%, respectively. These relatively low values imply that strain V-78 must be considered a subspecies of P. shigelloides.

In humans, most enteric Plesiomonas infections are self-limiting and do not need a specific antimicrobial treatment. Some severe dysenteric and chronic intestinal infections, however, benefit from antimicrobial therapy. In these cases, fluoroquinolones and trimethoprim are generally considered to be the best oral agents because they have been shown to shorten the clinical course of Plesiomonas diarrhoea [44,45]. For fish, because P. shigelloides is an agent that seems to be an opportunistic pathogen or host of microbiota, veterinarians do not need an antimicrobial treatment when the mortality rate is low. However, empiric treatments have been done in high mortality cases due to lacking identification or ignoring this species. Unfortunately, CLSI or EUCAST do not publish a breakpoint for antimicrobials commonly used in fish health rather than ciprofloxacin and Trimethoprim/Sulfamethoxazole. We have tested our strains on 14 different antimicrobial agents belonging to nine groups; unfortunately, we only determined our isolates as susceptible or resistant to only two agents. The isolates commonly tested and reported for antimicrobials used in human health as susceptible or not, but few reports suggested antimicrobials for fish isolates [46]. Due to the ubiquitous characteristics of P. shigelloides, we do not suggest an empirical treatment or antimicrobial usage for fish strains, because the strains could easily gain antimicrobial resistance, as revealed by the genomic island analysis. We also suggest that epidemiological cut-off values (ECVs) should be established to determine the antimicrobial susceptibility of P. shigelloides strains. While there is a significant gap in the antimicrobial susceptibility values of P. shigelloides, we detected 19 to 34 antimicrobial resistance genes in the genome of P. shigelloides strains we employed. Unlike other strains obtained from the GenBank database, the genome of strain V-78 harbours a relatively low number of antimicrobial resistance genes.

As an important opportunistic pathogen, the pathogenesis of P. shigelloides can be correlated to multiple virulence genes, which encode secreted proteins and toxins. A variety of virulence factors have been associated with the pathogenesis of this microorganism, including β-hemolysins [47,48,49], enterotoxins [50], cholera-like toxins [51] and possible endotoxins [52]. Among them, ast or act encode enterotoxins in many bacteria, which damage intestinal epithelial cells, lyse red blood cells, and cause gastroenteritis, water-like diarrhea and death [53,54]. In our studies, the analyzed P. shigelloides strains harbored 20 common virulence genes including cheD, flhA, fliI, fliN, flip, gmhA/lpcA, gspD, gspE, gspG, htpB, IlpA, katB, kdsA, lpxC, lpxD, luxS, mgtB, clpC, rfaD, xcpR, which may be related to the pathogenicity in humans, animals and fish. Although various virulence gene determinants were determined in P. shigelloides strains that originated from different sources, this is the first comprehensive report presenting virulence genes in detail. We also present all virulence genes, totalling 83, found in genomes of all employed strains. Our strain, V-78, has a difference for antimicrobial resistance genes from other employed strains, but all strains have a similar count of virulence genes between 39 to 43. The findings for virulence genes suggested that P. shigelloides has an important threat to all living organisms by means of similar counts but different virulence genes.