VV2018 was isolated from the blood culture specimen of a clinical patient with hepatitis B virus cirrhosis in 2018, in Nantong, Jiangsu Province, China (Wu et al., 2023). The clinical blood sample was used with the approval of the Ethics Committee of Affiliated Nantong Hospital 3 of Nantong University. The strain was identified using the bioMérieux VITEK 2 compact instrument (bioMérieux, Marcy-l’Étoile, France) and average nucleotide identity (ANI) analysis. The genomic DNA of VV2018 was extracted using a TIANamp Bacteria DNA Kit (Tiangen Biotech Company Ltd., Beijing, China), according to the manufacturer’s protocol.

The antibiotic resistance profiles were assessed through minimal inhibitory concentration (MIC) assays (Liu and Crosa, 2012). Briefly, a final suspension of 10 cfu/mL in broth supplemented with 2% NaCl and 1 mM CaCl₂-H₂O were distributed in triplicate throughout a 96-well microtiter plate. Escherichia coli ATCC 25922 was used as the susceptible-control reference bacterial strain for MIC assays. Cells were challenged with 0.25–1,024 g/mL antibiotics. MICs were determined by detection of cell pellet formation in the bottom of the wells of the 96-well plate by turbidometry at 600 nm using Multiskan GO (Thermo Fisher Scientific, USA). Drug susceptibility was determined according to the Clinical and Laboratory Standards Institute (CLSI) drug susceptibility test standard from 2018.

Whole genome sequencing (WGS) and assembly were conducted at Azenta Life Sciences (Suzhou, China). Sequences of VV2018 were obtained using PacBio Sequel platform (Pacific Biosciences, Menlo Park, CA, USA) and Illumina HiSeq X Ten platform (Illumina, San Diego, CA, USA). The PacBio reads were assembled by Hifiasm v0.13-r308 and Canu v2.2 (Koren et al., 2017; Cheng et al., 2021), and then the Illumina reads were mapped onto the assembled contigs to correct the primary assembly and control assembly quality using Pilon 1.22 and Quiver (Chin et al., 2013; Walker et al., 2014). The genome completeness and contamination of all V. vulnificus strains we used were further evaluated by checkM with default settings (Parks et al., 2015). Prokka 1.14.6 was used to predict potential CDSs (Seemann, 2014). The functional annotation of these CDSs was performed by DIAMOND (Buchfink et al., 2021) against the non-redundant protein sequence (NR) database of the National Center for Biotechnology Information (NCBI). Kyoto Encyclopedia of Genes and Genomes (KEGG), Cluster of Orthologous Groups (COG), UniProt/Swiss-Prot, Pfam, CAZymes, virulence factors of pathogenic bacteria (VFDB) and Antibiotic Resistance Genes Database (ARDB) were also used to annotate the functions of CDSs (Liu and Pop, 2009; Levasseur et al., 2013; Kanehisa et al., 2017; Galperin et al., 2021; Mistry et al., 2021; Liu et al., 2022; UniProt Consortium, 2023). The rRNA and tRNA sequences were annotated by RNAmmer (Lagesen et al., 2007) and tRNAscan-SE (Lowe and Chan, 2016), respectively. The mobile genetic elements (MGEs) were annotated using ISfinder (Siguier et al., 2006). Genomic islands (GIs), prophages, and CRISPR-Cas systems were identified using online tools IslandViewer 4, PHAge Search Tool (PHAST), and CRISPRCasFinder software, respectively (Arndt et al., 2016; Bertelli et al., 2017; Couvin et al., 2018). Multilocus sequence typing (MLST) was performed by analyzing the housekeeping genes on the MLST website.¹ The basic characteristics of the chromosomes were visualized by the CGView Comparison Tool (Petkau et al., 2010).

All the available V. vulnificus complete genome sequences (n = 25) were downloaded from the NCBI database with checkM values (Supplementary Table 1). The genome sequences were re-annotated with Prokka 1.14.6 and pan-genome analysis was conducted based on the output of Prokka using Roary with a BLASTP identity cutoff of 90% (Page et al., 2015). For genome similarity assessment, digital DNA-DNA hybridization (dDDH) values were computed using web tool GGDC 3.0 (formula 2, identities/HSP length) (Meier-Kolthoff et al., 2022). Whole genome ANI between pairwise V. vulnificus strains was calculated with Pyani software available at https://github.com/widdowquinn/pyani. The core genome of these strains was produced by Harvest software v1.1.2 (Treangen et al., 2014) using the V. vulnificus CMCP6 genome as a reference. Recombination events were removed from the core-genome alignment using Gubbins v2.2.0 (Croucher et al., 2015). Single nucleotide polymorphisms (SNPs) were then extracted from the recombination-free core genome alignment using the script available at https://github.com/sanger-pathogens/snp-sites. The multi-alignments were aligned with the ClustalW in MEGA 11.0 and analyzed using GeneDoc 2.7.0 (Nicholas and Nicholas, 1997; Tamura et al., 2021). The maximum likelihood (ML) phylogenetic tree of SNPs was constructed using RAxML in the GTRGAMMA model (1,000 bootstrap) (Stamatakis, 2014) and was visualized using Figtree v1.4.4.² The neighbor-joining phylogenetic trees of PglC and PglD performed by MEGA 11.0. CD-HIT used to cluster the retained sequences using the genome sequence of VV2018 as the reference with identity of 80% and coverage of 90% (Li and Godzik, 2006).

The complete genome of VV2018 comprised two chromosomes, Chr I and Chr II. The genome completeness of VV2018 was 100%, and the contamination was 0.05% accessed by checkM. Chr I consisted of 3,264,146 bp with a GC content of 46.60% containing 2,874 predicted CDSs, 106 tRNA genes and 31 rRNA genes. Chr II consisted of 1,816,653 bp with a GC content of 47.19% containing 1,515 predicted CDSs, 13 tRNA genes and 3 rRNA genes (Table 1 and Figure 1). MLST revealed that VV2018 belonged to a novel ST620 and was very close to ST387, with eight loci, isolated from humans in China.

The distribution of VV2018 CDSs into COG functional categories is shown in Supplementary Figure 1. Except for genes with unknown functions (8.57%), most genes were related to signal transduction mechanisms, amino acid transport and metabolism, transcription, and carbohydrate transport and metabolism. The annotation of genes of VV2018 in KEGG pathway analysis showed that the most genes were involved in metabolism, including carbohydrate metabolism, amino acid metabolism, metabolism of cofactors and vitamins, and energy metabolism (Supplementary Figure 2).

A total 151 putative VFs were predicted among Chr I (106, 70.20%) and Chr II (45, 29.80%). These genes were mainly associated with motility (polar flagellar proteins), immune evasion (capsular polysaccharide and iron uptake), secretion system (type II secretion system proteins), adherence (type IV pilus, lipooligosaccharide and OmpU) and toxin (RTX toxin) (Table 2). Iron uptake from host cells plays a key role in the survival of V. vulnificus (Dittmann et al., 2019). RTX toxin (rtxABCD) and OmpU have been shown to play important roles in the infection and pathogenesis, respectively, of V. vulnificus (Goo et al., 2006; Liu and Crosa, 2012). The vvhA and tlh were another two toxin genes, encoding cytolysin-hemolysin and thermolabile hemolysin, which induced acute cell death and were important in the pathogenesis and dissemination of these bacteria (Wang et al., 2015; Song et al., 2016). In this case, the patient’s temperature was 40.2°C after infection, which was accompanied by chills, unbearable low back pain, and forced position. Thus, the serious infection of this case might have a strong relationship with the mixing effect of multiple VFs (Wu et al., 2023).

Six antibiotic resistance-related genes were identified in the genome of VV2018, including dfrA3 (encoding a dihydrofolate reductase), qnrVC1 (encoding a pentapeptide repeat protein), catB9 (encoding a type B-5 chloramphenicol O-acetyltransferase), tet (Bertelli et al., 2017) (encoding a tetracycline efflux pump), crp (encoding a cAMP-receptor protein) and varG (showing resistance to penicillin, carbapenems and cephalosporins in vitro). CRP is a global regulator that not only regulates the expression of the multidrug efflux pump but also impacts the expression of multiple VFs (Nishino et al., 2008; Zhan et al., 2008). The antibiotic resistance pattern of VV2018 is shown in Table 3. This isolate was susceptible to most tested antibiotics, including tetracycline and chloramphenicol, with the exception of cefazolin. The fact that VV2018 showed susceptibility to tetracycline may be the result of acetylation-mediated down-regulation of tetA gene (Pang et al., 2020). The study also reported that all of Vibrio cholerae strains harboring catB9 gene were susceptible to chloramphenicol (Lepuschitz et al., 2019). Further work needs to be done to study catB9 in Vibrio. V. vulnificus has been reported to show complete resistance against cefazolin (Pan et al., 2013). The varG gene has been shown to have beta-lactamase activity against penicillin, carbapenems, and cephalosporins in vitro (Lin et al., 2017), which might account for the cefazolin resistance of VV2018.

Large parts of the genome designated as genomic islands (GIs) and phages were transferred from one bacterium to another (Canchaya et al., 2003; Dobrindt et al., 2004). Twenty GIs and nine GIs were detected on Chr I and Chr II in VV2018, respectively (Supplementary Tables 2, 3). The length of GIs on Chr I ranged from 4 kbp to 163 kbp. In the GIs of Ch I, a total of 6 transposase genes were predicted, all of which were classified into the IS4, IS481 and IS5 families. Meanwhile, one integrases (intS) and two tyrosine recombinases (xerC and xerD) were encoded. Two genes were predicted to encode type I restriction enzyme proteins. The length of GIs on Chr II ranged from 4 kbp to 86 kbp, containing one tyrosine recombinase-encoding genes (xerC). Meanwhile, sulfate permease genes (cysTWA) which allowed the bacteria survive in selenite environment by decreasing the expression, were found on Chr II_GI2 (Tempel et al., 2022). Only one incomplete prophage sequence was predicted on Chr I with a length of 9.7 kbp encoding genes with unknown functions (Supplementary Table 4). Thus, further work needs to be done to investigate the functions of these genes.

One CRISPR locus was predicted without Cas genes on Chr I (Supplementary Table 5), the phenomenon that CRISPR locus without Cas genes was also found in other Vibrio strains, Listeria monocytogenes and Staphylococcus, indicating that it was unable to effectively exert adaptive immunity (Mandin et al., 2007; Zhang et al., 2019, 2021). There were four direct repeats with a length of 32 bp and three spacers. The sequences of spacers closely matched other V. vulnificus strains in the NCBI database.

The ANI and dDDH values of 26 V. vulnificus strains are summarized in Supplementary Table 6. The ANI value between VV2018 and other V. vulnificus strains was 97.14% (range 95.41 to 98.45%) (Figure 2A). The most similar strain compared to VV2018 was V. vulnificus FORC_017 (98.45% identity) isolated from human in South Korea. The heatmap showed that all 26 V. vulnificus strains were divided into two clusters, most strains including VV2018 in cluster 1 were isolated from human, and most strains in cluster 2 were isolated from seafood or unknown places. Meanwhile, the phylogenetic tree of all 26 V. vulnificus strains constructed on the basis of the core SNPs showed that VV2018 was grouped with two other V. vulnificus strains isolated from human (V. vulnificus FORC_009 and V. vulnificus FORC_016) (Figure 2B). The dDDH values among 26 V. vulnificus strains were more than 60%. An ANI cut-off of around 95% did not correspond to an absolute dDDH value (70% cut-off for dDDH). A previous study reported that the value of 70% dDDH could not be used as absolute boundary, but still a gap between 60 and 70% similarity seemed to embrace clear-cut clusters of organisms, given the large extent of diversity among prokaryotes (Richter and Rosselló-Móra, 2009).

VV2018 shared a total of 3,016 core genes and 314 soft core genes with other 25 V. vulnificus strains according to the pan-genome analysis. A total of 138 strain-specific genes, accounting for 3.14%, were identified in VV2018. The functions of the majority of VV2018 specific genes (78.26%) were unknown, the other specific genes were involved in functional categories of replication/recombination/repair (7.97%) and cell wall/membrane/envelope biogenesis (4.36%) (Supplementary Table 7).

A total of 180 specific VFs were identified in 26 V. vulnificus strains and 115 VFs were included in core genes. The heatmap based on the presence and absence of all VFs clearly showed that the distribution of virulence genes differed between V. vulnificus strains (Figure 3A), and the VFs in VV2018 were similar to those in V. vulnificus 07-2444. Meanwhile, the differences in VFs between strains were among adherence and the immune system. In addition, some V. vulnificus strains isolated from humans were closely clustered with those isolated from seafood, indicating that these strains may cause foodborne infection.

Due to the misuse of antibiotics, V. vulnificus in seafood and aquatic environments are exhibiting resistance to multiple antibiotics (Elmahdi et al., 2016). V. vulnificus resistance toward common antibiotics has reached alarming levels in many countries which has serious implications for the treatment methods for bacterial infections (Heng et al., 2017). The distribution of antibiotic resistance genes in the V. vulnificus strains was also investigated. All 26 V. vulnificus strains possessed the resistance-related genes dfrA3, varG, tet (Bertelli et al., 2017), qnrVC1 and crp (Figure 3B). The prevalence of dfrA3, varG, tet (Bertelli et al., 2017), qnrVC1, and crp genes in these strains suggested that these genes may increase the resistance of these strains to trimethoprim, penicillin, tetracycline, quinolone and oxacillin. The phenicol resistance gene catB9 were also in most genomes of V. vulnificus strains, except V. vulnificus CECT 4999 and V. vulnificus MO6-24/O. In addition, among these strains, V. vulnificus VV2014DJH carried the fosfomycin resistance gene fosC2. Moreover, V. vulnificus VV20-8B-2 isolated from seafood possessed the most antibiotic resistance genes than other V. vulnificus strains, indicating that its antibiotic resistance may be more extensive.

The genetic context of the resistance-related genes of VV2018 was almost the same compared with other 25 V. vulnificus strains, except the varG gene. The gene varG might account for the cefazolin resistance of VV2018, and showed resistance to penicillin, carbapenems, and cephalosporins in vitro (Lin et al., 2017), however, the genetic context of varG was unknown in V. vulnificus strains. The varG gene was present in all 26 V. vulnificus strains, and the sequences containing varG were clustered into eight clusters with coverage of 90% and identity of 80% (Supplementary Table 8). The largest cluster was cluster 3 containing 12 sequences, and most of them (8/12, 75.0%) were isolated from humans including VV2018. Eight representative sequences were chosen for further analysis (Figure 4). The results of this gene neighborhood analysis of representative sequences revealed that a few transposase genes (IS30, IS5 and IS110) were upstream and downstream of varG. The genetic context of varG was conserved and the genes upstream and down of varG were ompV and nodD that were present in 96.15% of the sequences. These observations indicated that sequence rearrangement rarely occurs in the varG-encoding region, and genetic commonalities of ompV-varG-nodD across sources strongly suggested the structure of ompV-varG-nodD was conserved in V. vulnificus strains.

Two (LNNJENCE_00235 and LNNJENCE_00237) of six genes on Chr I_GI4 were glycosyl transferases, belonging to the GT4 family (Supplementary Table 9). LNNJENCE_00237 was compared with PglA in Neisseria including N. elongata subsp. glycolytica, with 37% identity. Meanwhile, four genes downstream of LNNJENCE_00237 also shared a high identity compared with Pgl-related proteins in Neisseria such as N. elongata subsp. glycolytica ATCC29315 (Table 4 and Figure 5). Broad-spectrum O-linked protein glycosylation (pgl) systems have been defined in Neisseria, such as N. gonorrhoeae, N. meningitidis, and N. elongata subspecies glycolytica (Ku et al., 2009; Vik et al., 2009; Naess et al., 2023). The pglA and pglBCD may be involved in pilin glycosylation (Power et al., 2003).

We selected six sequences of Neisseria for further analysis and the result of the multiple sequence alignment showed that the sequences of PglC and PglD were highly conserved (Figure 6A and Supplementary Figure 3A). The phylogenetic tree of PglC and PglD showed that the PglC of VV2018 was closest to the protein of N. weaveri LMG 5135, and PglD of VV2018 was closest to N. bacilliformis ATCC BAA-1200 and N. elongata subsp. glycolytica ATCC29315 (Figure 6B and Supplementary Figure 3B). Meanwhile, comparative genomic analysis revealed that the sequence of VV2018 containing pgl-like genes was similar to the sequences of V. vulnificus 07-2444, V. vulnificus YJ016 and V. vulnificus FORC_017 (Figure 7A). The result showed that the pglA gene was only found in V. vulnificus 07-2444 with coverage of 76% and identity of 85%. The other genes pglBa, pglBb, pglC, and pglD also had high homology with identity ranging from 90 to 100%. In addition, in the genome of V. parahaemolyticus, the pglB2 gene is associated with a pglC and pglD homolog, suggesting that a complete glycosylation system might also be present (Chamot-Rooke et al., 2007). In V. cholerae, O-glycosylation via PglL_Vc and possibly RbmD could represent a fine-tuned feedback mechanism controlling release of type II secretion system (T2SS) effectors by modulation of secretion efficacy (Vorkapic et al., 2019). Therefore, pgl-like genes (pglABCD) may play important roles in Vibrio species.

The total domain structures of PglA, PglBa, PglBb, PglC and PglD were also analyzed and compared with protein sequences in Neisseria. The domains of PglA, PglBa, PglC and PglD were similar, except for PglBb (Figure 7B). The PglD_N domain was found in the N-terminus of all PglBb sequences, besides the number and location of hexapeptides. There were three hexapeptides in PglBb of VV2018 and V. vulnificus 07-2444. The pglBb of V. vulnificus YJ016 and V. vulnificus FORC_017 had two hexapeptides in the second and third position. The pglBb of N. elongata subsp. glycolytica strain ATCC29315 had two hexapeptides in the first and third position. A number of different transferase protein families contain hexapeptide repeats, such as galactoside acetyltransferase-like proteins (Wang et al., 2002). It has been shown that most hexapeptide acyltransferases form catalytic trimers with three symmetrical active sites (Bergfeld et al., 2007). This is the first study to report that pglA and pglBCD genes were found in V. vulnificus strains. The pgl-like genes may affect the adherence of the strain, however, much more work is need to prove this hypothesis.