Pathogenomics of Virulence Traits of Plesiomonas shigelloides That Were Deemed Inconclusive by Traditional Experimental Approaches

One of the major challenges of modern medicine includes the failure of conventional protocols to characterize the pathogenicity of emerging pathogens. This is particularly apparent in the case of Plesiomonas shigelloides. Although a number of infections have been linked to this microorganism, experimental evidence of its virulence factors (VFs), obtained by traditional approaches, is somewhat inconclusive. Hence, it remains unclear whether P. shigelloides is a true or opportunistic one. In the current study, four publicly available whole-genome sequences of P. shigelloides (GN7, NCTC10360, 302-73, and LS1) were profiled using bioinformatics platforms to determine the putative candidate VFs to characterize the bacterial pathogenicity. Overall, 134 unique open reading frames (ORFs) were identified that were homologous or orthologous to virulence genes identified in other pathogens. Of these, 52.24% (70/134) were jointly shared by the strains. The numbers of strain-specific virulence traits were 4 in LS1; 7 in NCTC10360; 10 in 302-73; and 15 in GN7. The pathogenicity islands (PAIs) common to all the strains accounted for 24.07% ORFs. The numbers of PAIs exclusive to each strain were 8 in 302-73; 11 in NCTC10360; 14 in GN7; and 18 in LS1. A PAI encoding Vibrio cholerae ToxR-activated gene d protein was specific to 302-73, GN7, and NCTC10360 strains. Out of 33 antibiotic multi-resistance genes identified, 16 (48.48%) genes were intrinsic to all strains. Further, 17 (22.08%) of 77 antibiotic resistance islands were found in all the strains. Out of 23 identified distinct insertion sequences, 13 were only harbored by strain LS1. The number of intact prophages identified in the strains was 1 in GN7; 2 in 302-73; and 2 in NCTC10360. Further, 1 CRISPR element was identified in LS1; 2 in NCTC10360; and 8 in 302-73. Fifteen (78.95%) of 19 secretion systems and secretion effector variants were identified in all the strains. In conclusion, certain P. shigelloides strains might possess VFs associated with gastroenteritis and extraintestinal infections. However, the role of host factors in the onset of infections should not be undermined.

Furthermore, suitable models for studying the pathogenesis of Plesiomonas have not yet been established and traditional models proved to be inappropriate. The development of polymerase chain reaction and other rapid low-cost protocols, already used for the analysis of other enteropathogens, for the detection of Plesiomonas VFs is also urgently required. Presently, no specific primer-and molecular-based techniques for the detection of the major, if not all, predicted Plesiomonas VFs exist. Further, strain-specific and rapid systems for the delineation of pathogenic and non-pathogenic Plesiomonas strains are lacking.
Since the emergence of new pathogens (including Plesiomonas) is linked to complex interrelated factors, of which mutation and horizontal gene transfer are the focal driving elements (Che et al., 2014), it is not implausible that several such events occurred in Plesiomonas. In most cases, gene acquisition or loss could impact genome evolution (Che et al., 2014), converting a non-pathogenic strain into a pathogenic one (McDaniel and Kaper, 1997). Although these phenomena may have contributed to the emergence of pathogenicity in Plesiomonas, no studies focusing on this issue have been reported. The study of PAIs, resistance islands, metabolic islands, mobility genes, IS elements, phage-related genes, tRNA genes, and direct repeats will be important for the characterization of the virulence and genome evolution of Plesiomonas.
Consequently, rapid specific culture-independent tools for pathogenicity characterization should be developed for Plesiomonas. A comparative in silico profiling of the available genomes may provide important clues prior to wet-lab studies for the assessment of virulence marker candidates in Plesiomonas and answer other questions about its emerging pathogenicity capabilities. This would not only reveal the as yet unknown virulence traits associated with specific strains, but also inform the development of rapid and inexpensive tools for strain screening and diagnostic characterization of plesiomonads. In the current study, we therefore profiled the publicly available complete genome sequence assemblies of Plesiomonas in silico to identify the associated type III, IV, VI, and VII secretion systems (T3SS, T4SS, T6SS, and T7SS, respectively); integrative and conjugation elements (ICE elements); prophages; VFs; ARs; ARIs; type III, IV, VI, and VII secretion effectors (T3SE, T4SE, T6SE, and T7SE, respectively); class I integrons, IS elements; and CRISPR.

Whole-Genome Sequence Assembly Selection
The four publicly available sequence assemblies of Plesiomonas were downloaded from NCBI website via ftp.ncbi.nlm.nih.gov/ genomes/ASSEMBLY_BACTERIA (Table 1; last accessed [20 December 2017]).

Detection of Virulence Traits
Presence of PVFs within each genome assembly was predicted using the VRprofile 2.0 (Li et al., 2018). For gene clusters (including genes encoding T3SS, T4SS, T6SS, and T7SS), ICE, and prophage prediction, default parameter settings were used. BLASTP searches were done with a minimum e-value of 1 × 10 −4 ; HMMer search E-value of 1 × 10 −4 for ICE and prophage detection; the minimum number of significant hits for each prophage-like region set to 5; and a minimum length of direct repeats of an ICE-like region of 15 bp. Single genes related to VFs, ARs, T3SE, T4SE, T6SE, T7SE, class I integrons, IS elements, GI-like region, PAIs, and ARIs were predicted by setting the BLAST Ha-value to 0.64 (Li et al., 2018). The distribution of PVFs and related genes among the various strains was evaluated using a web-based Venn diagram program 1 .
In addition, a direct search of the WGS for presence of elements related to already suggested virulence mechanisms in P. shigelloides such as heme iron uptake systems or iron-siderophores (Daskaleros et al., 1991;Santos et al., 1999;Henderson et al., 2001;Gonzalez-Rodriguez et al., 2007;Oldham et al., 2008;Rodríguez-Rodríguez and Santos, 2018) was performed.

Detection of Prophage Sequences
To identify prophages within the Plesiomonas genomes, strain sequences were analyzed using the PHAge Search Tool Enhanced Release (Zhou et al., 2011;Arndt et al., 2016).

Determination of the CRISPR Elements
CRISPR elements in the Plesiomonas genomes were predicted using CRISPR Finder by Random Forest program (Wang and Liang, 2017

Detection of Plasmids
Plasmid elements in the Plesiomonas genomes were predicted using PlasmidFinder 1.3 (Carattoli et al., 2014) and a direct search for replicon information on NCBI website 2 .

Putative Virulence Factor
Out of 134 unique PVFs identified in all the strains, strain 302-73 harbored 98 PVFs; strain GN7 harbored 107 PVFs; strain LS1 harbored 97 PVFs; and strain NCTC10360 harbored 90 PVFs. Detail information on the ORFs, amino acid length of predicted encoded proteins, Ha-values, and functional description of the homologous or orthologous predicted genes are provided in the Supplementary Data Sheets S1-S4.
The distribution of IUSRVTs among the strains is presented in Supplementary Figure S1. Of these, 9 IUSRVTs were common to all the strains. The common IUSRVTs were multispecies ubiquinol-cytochrome c reductase iron-sulfur subunit, multispecies iron donor protein CyaY, multispecies iron-sulfur cluster insertion protein ErpA, multispecies succinate dehydrogenase/fumarate reductase iron-sulfur subunit, zinc/iron-chelating domain-containing protein, iron export ABC transporter permease subunit FetB, iron-sulfur cluster repair di-iron protein, ferrous iron transport protein B, multispecies iron-sulfur cluster assembly protein IscA. Two IUSRVTs namely multispecies zinc/iron-chelating domain-containing protein and iron transporter FeoA were only found in strains 302-73, LS1, and NCTC10360. Similarly, the iron ABC transporter permease was identified in strains 302-73, GN7, and NCTC10360. Further, strains 302-73, GN7 and LS1 shared 2 IUSRVTs, including multispecies succinate dehydrogenase iron-sulfur subunit and multispecies iron-sulfur cluster scaffold-like protein.
IUSRVTs that were specifically found only in strains 302-73 and LS1 were iron(III) ABC transporter ATP-binding protein.
The distribution of HUSRVTs among the P. shigelloides strains are summarized in Supplementary Figure S2. The various HUSRVTs common to all the strains ranged from heme ABC transporter ATP-binding protein, heme exporter protein CcmB, heme utilization protein HutZ, to heme
Antibiotic Resistance Island Figure 4 shows the distribution of ARIs among the strains. The various ARIs common to the strains were related to

IS Elements
The IS elements identified in the four P. shigelloides strains are presented in Table 6. The IS elements were shared by the strains or specific to a given strain, as follows: ISEhe3_PEP, and ISAbo1_PEP3 were identified in strains 302-73, GN7, and LS1; ISSm4_PEP2 was identified in strains NCTC10360 and LS1; IS1N_PEP was identified in strains 302-73 and LS1; ISIba1_PEP3 was identified in strains GN7 and LS1; ISEc16_PEP3, IS3F_PEP, and ISAs1_PEP were only identified in strain 302-73; ISErsp1_PEP and ISIba2_PEP were only identified in strain GN7; ISShes11_PEP5, ISEic1_PEP3, IS15DII_PEP, ISIba2_PEP3 ISSag9_PEP, ISStso1_PEP, ISStso1_PEP3, ISKpn25_PEP3, IS1381A_PEP3, IS1381A_PEP ISEic1_PEP ISSm4_PEP, and ISEhe3_PEP3 were only identified in strain LS1. The distribution of IS elements among the strains is depicted in Figure 6. LS1 harbored 56.52% (13) of all IS elements identified in the strains.

CRISPR Elements
No CRISPR element was identified in strain GN7. Various features of CRISPR elements harbored by the other three Plesiomonas strains are detailed in Table 8. Strain NCTC10360 harbored two CRISPR elements (172 and 185 bp). The 172-bp element harbored three repeated sequences and two spacers, with a repeat length of 25 bp. The 185-bp element harbored four repeated sequences and three spacers, with

Prophage Entities
The results of comparative phageomics of the four P. shigelloides strains are shown in Table 8. Strain NCTC10360 harbored two intact phages sized 21.9 kb and 53.8 kb, most similar to PHAGE_Entero_mEp235_NC_019708, PHAGE_Klebsi_phiKO2_NC_005857, PHAGE_Yersin_PY54_ NC_005069, and PHAGE_Entero_lambda_NC_001416. Strain GN7 harbored two intact and one incomplete phage sized 36.2, 30.3, and 29.5 kb respectively. The GN7 phages were most similar to the following phage species: PHAGE_Flavob_1H_NC_031911,

Plasmid Elements
No plasmid element was identified in the four strains. Replicon information of the WGS revealed chromosomes.

DISCUSSION
One major problem of the conventional elucidation of pathogenesis mechanisms of Plesiomonas is the lack of a suitable experimental animal model. In the current study, comparative pathogenomics approach was used to detect virulence traits of the emerging pathogen that, until now, had been unresolved by traditional approaches. As an example, using an adult rabbit model, the ileal loop system, Sanyal et al. (1980) had established the diarrheagenic (enterotoxigenic) potential of Plesiomonas but did not model invasive pathogenesis. Further, pathogenomics  studies revealed the presence of several genes associated with adherence and attachment to the villi of enterocytes (including thin aggregative fimbriae genes, such csgE, csgF, csgG, gmhA, pilT, and pilU) (Collinson et al., 1996;Sukupolvi et al., 1997;Römling et al., 1998;Comolli et al., 1999a;Bacon et al., 2001), characteristics necessary for diarrheagenic activity of pathogens.
Other core pathogenic traits that remained unidentified or were only weakly supported by conventional experimentation include the cytotoxic and hemolytic potential of Plesiomonas. According to cell culture and animal model studies of Plesiomonas cytotoxicity, some strains produce extracellular cytotoxins with activities against Vero cells (Olsvik et al., 1990). Similarly, vacuolation and cytotoxigenic activities of plesiomonad strains recovered from different environments, including human/animals and aquatic bodies, in Vero cell monolayers were described by Ekman (2003) and Falcón et al. (2003), respectively. Although there is some evidence (as mentioned), other studies yield contradicting data. Vitovec et al. (2001) were unable to demonstrate either cytotoxic or hemolytic activity of Plesiomonas isolates during both, monoinfection and C. parvum co-infection in experimental models. Likewise, Pitarangsi et al. (1982) and Holmberg and Farmer (1984) concluded that the tested Plesiomonas strains are not cytotoxic in the animal model used. Further, the plesiomonad strains studied by Maluping et al. (2004) exhibited no cytotoxicity or hemolytic potential in any of the experimental models explored. Nonetheless, in the current study, we identified some plesiomonad ORFs homologous or orthologous to genes, such as ssb, known to mediate toxin secretion similar to E. coli hemolysin secretion (Wong et al., 1998;Marcus et al., 2000), and xcpA/pilD, known to facilitate the internalization and cytotoxicity in P. aeruginosa (Hahn, 1997;Comolli et al., 1999b). These suggest the cytotoxicity potential of Plesiomonas.
The LOS gene from Haemophilus that showed homology with a Plesiomonas ORF is linked to endotoxin and immunogenic functions, and its phosphoryl choline impacts the invasion via interaction with PAF receptor; and stimulates inflammatory signals, evasion of antigen-specific host defenses, and colonization of diverse microenvironments within the host (Rao et al., 1999;Erwin et al., 2000;Swords et al., 2000). A homologous ORF might provide Plesiomonas with a related pathogenic ability. Other genes related to toxin production and enteropathogenesis with homology to Plesimonas ORFs include waaC and waaF. These genes function in antiphagocytosis and resistance to serum killing (Rocchetta et al., 1999;Lyczak et al., 2000;Schroeder et al., 2002). Further, Plesiomonas possesses ORFs that are homologous to Vibrio PAI genes (tagD, nanE, and nanA). These Vibrio VFs are associated with cholera toxins and a putative protease-modulating VF in toxigenic Vibrio species (Karaolis et al., 1998;Jermyn and Boyd, 2002;Faruque and Mekalanos, 2003;Schmidt and Hensel, 2004). However, Herrington et al. (1987), using early traditional approaches, were unable to establish Plesiomonas enteropathogenesis in an animal model.
The presence of ARIs and clinically relevant antibiotic (multi)resistance genes might limit the treatment and management options for the infections caused by Plesiomonas. This could have potentially serious consequences since Plesiomonas infections are not routinely diagnosed in the medical setting (Chen et al., 2013). Similar to many cases of antimicrobial resistance shown by various microorganisms, infections involving resistant Plesiomonas might be very difficult to manage, leading to extended hospitalization, high morbidity, and, ultimately, elevated mortality rates (Kain and Kelly, 1989b). Multi-resistance of zoonotic Plesiomonas could pose increased threat to human and veterinary health. This might further raise the associated possibility of dissemination and contact between human and animals, especially pets and livestock. Plesiomonas are by nature, β-lactamase producers (resistant to β-lactam drugs) (Meyers et al., 1976;Brenden et al., 1988;Fisher et al., 1988;Avison et al., 2000;Stock and Wiedemann, 2001a). Plesiomonad tetracycline resistance has been reported by Kain and Kelly (1989a); González et al. (1999), Wong et al. (2000), and Suthienkul et al. (2001). Co-resistance of plesiomonad strains to antibiotics, including ampicillin, amoxicillin/clavulanic acid, trimethoprim/sulfamethoxazole, chloramphenicol, streptomycin, tetracycline, erythromycin, sulphametoxazole, and neomycin has also been documented (Olsvik et al., 1990;González et al., 1999;Wong et al., 2000;Jun et al., 2011). Chloramphenicol-, co-trimoxazole-, erythromycin-, and rifampin-resistant Plesiomonas strains are frequently reported (Wong et al., 2000;Stock and Wiedemann, 2001b;González-Rey et al., 2004). Tobramycinand gentamicin-resistant strains have also been documented (Clark et al., 1990). Tetracycline-resistant strains of Plesiomonas veterinary isolates are not uncommon. González-Rey et al. (2004) documented a multi-drug resistant Plesiomonas and tetracycline-resistant strain recovered from a Cuban dog and an environmental sample, respectively. Widespread tetracyclineresistant strains have also been noted among plesiomonad isolates from catfish ponds (DePaola et al., 1993). Further, aminoglycoside-resistant and penicillin-resistant Plesiomonas isolates were reported (Clark et al., 1990;Bravo-Fariñas et al., 1998;Avison et al., 2000;Stock and Wiedemann, 2001b). Yeh and Tsai (1991) have documented vancomycin-resistant strains of Plesiomonas. Furthermore, Jun et al. (2011) reported Plesiomonas resistance to amikacin, cefotaxime, cefepime, gentamicin, and ciprofloxacin. Plesiomonas inactivation of aminoglycoside antibiotics has been attributed to the production of specific enzymes, such as aminoglycoside-modifying enzymes (Shaw et al., 1993). Homology of Plesiomonas ORFs to heavymetal resistance genes closely related to arsenic resistance in Acinetobacter spp. suggests the bacterium's increasing ability to evade potential antiseptic treatment. Previously, Chao et al. (2016) reported cadmium-resistance in Plesiomonas, with a maximum tolerance concentration as high as 150 mg/L. Although Chao et al. (2016) concluded that this trait could be applicable in wastewater treatment, employed as an indicator of heavy metal pollution for eco-monitoring and environmental impact assessment, heavy metal resistance in Plesiomonas also suggests possible adverse implications for infectious cases.
The presence of Plesiomonas ORFs homologous to genes encoding different secretory and effector systems (T6SS, T3SS, T4SE, and T6SE) found in other pathogens might be indicative of a similar virulence function in the microorganism. For example, a Yersinia and Salmonella effector molecule that depolymerizes actin plays a similar role in tissue culture cells (Rosqvist et al., 1995). Secretory systems directly transport effector molecules from the cytoplasmic compartment to the cell surface to interact with and modify mammalian host cell proteins (Michiels et al., 1990). Other functional consequences of secretory systems in enteric and many other pathogens, such as Yersinia spp., Shigella spp., and Salmonella spp., include abolition of the immune cell function and engulfment ability, and macrophage protein modification (Rosqvist et al., 1990;Bliska et al., 1991;Cornelis, 1992); and enhanced pathogen entry into non-phagocytic cells (Galán, 1996;Ménard et al., 1996). Similar virulence traits are probably linked to the extraintestinal infections caused by Plesiomonas. Many genes harbored on or related to TSS elements, such as gspD-J, play roles in protecting pathogens against human complement activity and serum (Aiello et al., 2010;Johnson et al., 2016;Cianciotto and White, 2017;Waack et al., 2017). TSS elements are also involved in the mediation of cascades of secreted virulence effectors in gram-negative pathogens (Aiello et al., 2010;Johnson et al., 2016;Cianciotto and White, 2017).
Almost all Plesiomonas strains analyzed in the current study possessed PAIs homologous to those found in most enteropathogens. One of the PAIs, for instance, contained cadavarine decarboxylase gene (cadA). Probable alteration of the Plesiomonas cadA gene could be responsible for the virulence of some strains. Deletion of the cadA gene in Shigella reportedly promotes virulence to evade host cell protection against Shigella enterotoxin that relies on the inhibition of cadaverine synthesis (Maurelli et al., 1998;Brosch et al., 2001). Virulence plasticity is attributed principally to insertions and deletions in a pathogen genome (Brosch et al., 2001). Plesiomonas harbors IS elements that share homology with those found in pathogenic bacteria. The identified Plesiomonas IS elements could have selected for the pathogenic potential of this microorganism in many instances. IS elements are known mutagenic systems of the bacterial genome that modify genes by deletion, disruption, or up-regulation of neighboring genes (Mahillon and Chandler, 1998;Gyles and Boerlin, 2014). The identified IS elements might have impacted the evolutionary advantage of virulence in Plesiomonas, with both diarrheagenic and extraintestinal pathogenesis of this microorganism increasingly emerging.
The prophages detected in Plesiomonas strain genomes contribute to their genetic identity and possible differences of virulence traits. Prophage entities are known to impact bacterial pathogen fitness, from increasing pathogenicity potential, to chromosomal rearrangement and prophage-encoded VFs (Canchaya et al., 2004;Bondy-Denomy and Davidson, 2014).

CONCLUSION
As demonstrated in the current study, comparative genomics of Plesiomonas strains revealed the presence of probable virulence traits that could not be unambiguously determined using traditional approaches. The remarkable homology or orthology of Plesiomonas ORFs with VFs and antibiotic resistance genes from species commonly identified in the healthcare and food industries are non-accidental and further attest to the pathogenic capacity of plesiomonads. We hypothesize that differential expression and regulation of virulence traits in the presence of different host factors or under different conditions might be responsible for the lack of reproducibility (i.e., dissimilar results) of data on Plesiomonas pathogenicity obtained during early investigations that relied on conventional experimental models. That probably stems from the facts that the conditions that warrant the (over)expression of a Plesiomonas virulence gene in the host might not be replicable under experimental conditions, compromising the efficiency of conventional protocols for the assessment of the bacterial virulence potential. Expression of such virulence genes or validation of our hypothesis might be of interest to future transcriptomic investigations. Further, the current study provides clues on an alternative exploratory approach that could be adapted for strain typing and delineation of pathogenic and non-pathogenic variants of the bacterium. In-depth exploration of numerical and structural variation of CRISPR systems in Plesiomonas in combination with other specific high-throughput techniques could provide low-cost but sufficient discriminatory power invaluable for strain typing and virulence profiling. Thus, the quest for developing rapid and inexpensive tools for strain and diagnostic characterization of Plesiomonas becomes imperative. We conclude that P. shigelloides possesses some VFs associated with gastroenteritis and extraintestinal infections; however, the role of host factors in the onset of infections cannot be undermined.

DATA AVAILABILITY STATEMENT
All datasets generated for this study are included in the manuscript and the Supplementary Files.

AUTHOR CONTRIBUTIONS
TE conceptualized idea for the study, responsible for data collection, analysis, interpretation, and wrote the first draft of the manuscript. AO participated in result interpretation, proofread the first draft, and supervised the study. All authors read and approved the final manuscript.
DATA SHEET S5 | Structural and compositional details of CRISPRs obtained in P. shigelloides strain.