Vibrio metschnikovii as an emergent pathogen: analyses of phylogeny and O-antigen and identification of possible virulence characteristics

ABSTRACT

Vibrio metschnikovii is an emergent pathogen that causes human infections which may be fatal. However, the phylogenetic characteristics and pathogenicity determinants of V. metschnikovii are poorly understood. Here, the whole-genome features of 103 V. metschnikovii strains isolated from different sources are described. On phylogenetic analysis V. metschnikovii populations could be divided into two major lineages, defined as lineage 1 (L1) and 2 (L2), of which L1 was more likely to be associated with human activity. Meanwhile, we defined 29 V. metschnikovii O-genotypes (VMOg, named VMOg1–VMOg29) by analysis of the O-antigen biosynthesis gene clusters (O-AGCs). Most VMOgs (VMOg1 to VMOg28) were assembled by the Wzx/Wzy pathway, while only VMOg29 used the ABC transporter pathway. Based on the sequence variation of the wzx and wzt genes, an in silico O-genotyping system for V. metschnikovii was developed. Furthermore, nineteen virulence-associated factors involving 161 genes were identified within the V. metschnikovii genomes, including genes encoding motility, adherence, toxins, and secretion systems. In particular, V. metschnikovii was found to promote a high level of cytotoxicity through the synergistic action of the lateral flagella and T6SS. The lateral flagellar-associated flhA gene played an important role in the adhesion and colonization of V. metschnikovii during the early stages of infection. Overall, this study provides an enhanced understanding of the genomic evolution, O-AGCs diversity, and potential pathogenic features of V. metschnikovii.

KEYWORDS:

• Phylogeny
• O-antigen biosynthesis gene clusters
• pathogenicity
• genomics
• Vibrio metschnikovii

Introduction

The genus Vibrio consists of a group of bacterial species widely distributed in aquatic habitats. Vibrio species can cause diseases in fish, shellfish, mammals, and humans [1]. Of these, the species Vibrio cholerae, the etiological agent of cholera, is perhaps the best known, with illness associated primarily with strains in O-groups O1 and O139 that carry the genes for cholera toxin (CTX) [2]. “Non-cholera” Vibrio spp., including non-CTX producing V. cholerae within O-groups other than O1 and O139, V. parahaemolyticus, V. vulnificus, V. mimicus, V. fluvialis, and V. furnissii have been associated with mild to moderate gastroenteritis in humans; these species, together with V. alginolyticus, have also been linked with wound infections and, in susceptible hosts, bacteremia and sepsis [3–6]. While markers for virulence have been identified for some of these species/strains (i.e. the presence of genes for CTX in V. cholerae and for the thermostable direct hemolysin [TDH] in V. parahaemolyticus), the factors responsible for virulence in others are poorly understood, making it difficult to clearly identify strains that have an increased likelihood of causing human illness.

V. metschnikovii is an emerging environmental pathogen, capable of colonizing shellfish and causing infections in fish. Humans can become infected by ingesting the organism in fish or shellfish, or by exposure to water containing the organism, with human isolates reported from stool, blood, bile, and wounds [7–11]. In patients who are immunocompromised or have underlying liver disease, the microorganism is capable of causing bacteremia, sepsis, and death, despite what would appear to be appropriate antimicrobial therapy [8, 10, 12].

In the post-genome era, genome-based phylogeny has become a powerful tool for providing insights into the evolution and spread of Vibrio species [13, 14]. However, unlike other pathogens in the genus of Vibrio, the global genomic epidemiology, transmission dynamics, and pathogenicity determinants of V. metschnikovii are poorly understood. In a previous study, whole genome sequence data were obtained from strains of V. metschnikovii isolated from bird and water samples in China [15]. V. metschnikovii has also been isolated from a temperate-climate river basin of central Italy [16] and from surface water in southern and eastern Slovakia [17].

O-antigen analysis is a key element in the serotyping of Gram-negative bacteria. A systematic and large-scale analysis of V. metschnikovii O-antigen biosynthesis gene clusters (O-AGCs) is needed to better understand the diversity of their O-genotypes. In terms of pathogenicity, V. metschnikovii has been found to be pathogenic to freshwater-cultured hybrid sturgeon and pigs, horses, and cattle [7]. Some V. metschnikovii isolates from herring in the Norwegian Sea [18] were positive for cytolysin A (hlyA) and thermolabile hemolysin (tlh) genes, raising questions about the potential role of these and other virulence-associated genes in pathogenicity.

In this study, we investigate the global evolution, population structure, and O-AGCs diversity of the species and explore its pathogenicity and transmission dynamics from aquatic niche to human host. Our findings provided data critical to understanding this underestimated pathogen.

Materials and methods

Bacterial strains, genome sequencing and annotation

A total of 103 V. metschnikovii genomes were analyzed in this study (Table S1). Among them, 30 strains were newly sequenced; this includes 29 isolated from various sources and geographic locations in China during years of 2000–2020, and V. metschnikovii type strain JCM 21189^T. Another 73 whole genome sequences were acquired from the GenBank (https://www.ncbi.nlm.nih.gov/genbank/) and SRA (https://www.ncbi.nlm.nih.gov/sra/) database of NCBI (accessed in March 2022). Public genomes with less than 200 contigs were selected.

Genomes were sequenced by using the Illumina HiSeq 2000 platform. The reads data were assembled into contigs using SOAP de novo v2.04 [19]. QUAST v5.0.1 [20] software was applied to evaluate the genome assembly. Open reading frame (ORF) prediction and annotation were done by Prodigal v2.6.3 [21] and Prokka v1.13.3 [22] software, with Rapid Annotation using Subsystem Technology (RAST [23], http://rast.nmpdr.org).

Phylogenetic and pangenome analysis of V. metschnikovii

Intergenomic distances were calculated using the GBDP algorithm [24]. A principal co-ordinates analysis (PCoA) was performed using the matrix of V. metschnikovii accessory-genome organized with the gene present or absent. Core genome SNPs (cgSNPs) were called using snippy v4.3.6 software (https://github.com/tseemann/snippy) against the complete genome sequence of V. metschnikovii strain 2012V_1020 (assembly accession: GCA_009665235) as the reference genome. Gubbins v3.2 software (https://github.com/nickjcroucher/gubbins) was used as a recombination-removal tool to mitigate the effect of horizontal sequence transfer on phylogenetic reconstructions [25]. The FastBAPS algorithm was used for rapidly identifying genomic lineages [26]. The maximum-likelihood (ML) tree was constructed based on the non-recombinant sites using IQ-TREE v1.6.12 [27] based on cgSNPs with 1000 bootstrap replicates. The 103 V. metschnikovii genomes were annotated by Prokka v1.13.3 [22] and then all the *.gff files were collected as a local database. The pan-genome matrix was defined using Roary pan-genome pipeline v1.12 (https://sanger-pathogens.github.io/Roary) with a nucleotide identity threshold of 90%.

Identification and annotation of O-AGCs of V. metschnikovii

The sequences of the O-AGCs regions between gmhD and rjg of each V. metschnikovii genome were extracted. The RAST server pipeline was used to predict ORFs and annotate the O-AGCs genes. Web-BLAST (https://blast.ncbi.nlm.nih.gov/Blast) was applied to infer putative functions. Protein domains were found using Pfam program (https://pfam.sanger.ac.uk) with an E value threshold of 0.01 [28]. TMHMM v.2.0 (http://www.cbs.dtu.dk/services/TMHMM-2.0/) program was used to predict transmembrane domains. Insert sequences (IS) searches were performed using the ISfinder database (https://www-is.biotoul.fr/). In silico O-genotyping system was performed by BLASTn search using the wzx gene (or wzt gene) as a reference gene. Only a perfect match with an identity cutoff of >99% was considered.

Distribution of functional virulence genes and antimicrobial resistance genes

The virulence-associated genes were annotated based on the core dataset of virulence factor database (VFDB [29], http://www.mgc.ac.cn/VFs/). All the protein sequences were searched against the database using BLASTp with a relatively relax cutoff (identity ≥ 60%, length coverage ≥ 30%) and an E-value of 1e-5. The antimicrobial resistance genes were compared by using the comprehensive antibiotic resistance database (CARD, http://arpcard.mcmaster.ca). The BLASTn cutoff values were chosen as follow: an E-value of 1e-5, the identity ≥ 90% and length coverage ≥ 60%. Differences in resistance genes between lineages were statistically compared using the Fisher’s Exact Test.

In vitro cellular infectious, motility and biofilm formation assays

Bacterial infection was performed using the human epidermoid carcinoma cell line Hep-2 (CCC0068; Peking Union Medical College cell resource center) [30]. Cell adhesion and invasion assay was performed in 24-well plates as described previously [31] with some modifications. The lactate dehydrogenase (LDH) released by cells was detected using Cytotox 96 kit (Promega, Madison, USA). The criteria for classifying the cytotoxicity ranks referred to the previous study [30]. All experiments were performed three times in duplicate.

Growth curves of V. metschnikovii strains were determined by measuring the OD600 values after inoculation with fresh medium at ratio of 1/2000 (v/v). Bacterial swimming motility test as well as biofilm formation assays were performed as previously described [32]. Detailed procedures for in vitro cellular infectious, motility and biofilm formation assays are described in the supplementary data.

Construction of flhA mutant and complementary mutant

A V. metschnikovii mutant strain was constructed from clinical strain ZJ20140528 as a wild-type (WT). The suicide plasmid pWM91 [33] was used to construct a flhA mutant (ΔflhA) by homologous recombination. E. coli SM10 λpir served as the donor strain to transfer the recombinant plasmid into the WT by means of conjugative transfer. Incubation culture on LB agar (with 10% sucrose but without NaCl at 22°C for 48 h) was used to select the mutant strain. To construct the complementary mutants (C-ΔflhA), the entire coding region of flhA gene was cloned into plasmid pSRKGm [34]. Then the complementary plasmids were transformed into the above mutant with 0.1% arabinose as an inducer. The primers for the mutation and complementary protocol were described in Tables S2 and S3.

Data availability

The genome sequences of the 30 V. metschnikovi strains sequenced in this study have been deposited at GenBank/DDBJ/ENA under the BioProject ID no. PRJNA900004.

Results

Phylogenetics and population structure of V. metschnikovii

The 103 V. metschnikovii strains were clustered into two distinct lineages (L1 and L2) on the phylogenetic tree (Figure 1(a)) using the FastBAPS algorithm based on non-recombinant cgSNPs. Similarly, the GBDP-based phylogenomic tree also divided V. metschnikovii populations into two major lineages (Figure S1). Most strains in the L1 were isolated from clinical and food samples, while strains in the L2 were isolated from migratory birds. Using the criteria for defining a clonal group (CG) [35], 17 CGs were identified in this study (Figure 1(a)). It is noteworthy that the four clinical isolates were nested in L1, which were sourced from three countries (China, the United States, and the United Kingdom), indicating that strains in L1 are more likely to be pathogenic to humans, although the number of clinical strains was limited.

Figure 1. Population structure and phylogenetic analyses of V. metschnikovi. (a) Phylogenetic ML tree based on cgSNPs with 1000 bootstrap replicates. The red and blue branches represented Lineage 1 (L1) and Lineage 2 (L2), respectively. The purple dots on the top of the branches represented the 17 identified clonal groups (CGs). The inner ring showed different isolation sources, and the outer ring showed different isolation regions. (b) Correlation analysis between genetic distances and geographical distances within L1 strains. (c) Correlation analysis between genetic distances and geographical distances within L2 strains. (d) Splits-tree of phylogenetic networks based on the single-copy genes. (e) Heatmap of paired-SNP variants among 103 V. metschnikovii genomes. The colours from blue to red indicated the increase of SNP number. (f) Comparison of the ratios of six mutagenesis types between L1 and L2. (g) Chord diagram between drug resistance genes (DRG) and V. metschnikovii strains. A total of 17 DRG-positive isolates were identified, including 13 in L1 and 4 in L2.

We then assessed whether L1 of V. metschnikovii was preferentially associated with human activity using all 98 strains with well-defined isolation sites (including 59 strains of L1 and 39 strains of L2). By calculating the genetic distances and the actual geographical distances among V. metschnikovii isolates, we found that strains in L1 (including strains isolated from the global human community, i.e. food sources and clinical samples) were more affected by the factor of human activity (Mantel statistic r: 0.1325, P = 0.18781, Figure 1(b)) than strains in L2 (Mantel statistic r: 0.3351, P = 0.010989, Figure 1(c)), which were mainly isolated from migratory birds [15].

We re-examined all orthologous sequences of single-copy genes, resulting in a splits-tree of phylogenetic networks (Figure 1(d)). The two lineages were clearly distinguished in the spatial networks diagram. The heatmap of paired-SNP variants among the 103 genomes also confirmed the phylogenetic differences between the two major lineages (Figure 1(e)). More transitions than transversions were found in both lineages, and strains in L2 had a higher quantity of mutations than those in L1 (P < 0.001, Figure 1(f)), especially in the following two types of mutations (shift to C > T/G > A and shift to A > G/T > C).

The distribution of antibiotic resistance genes of V. metschnikovii also strongly indicated a clear evolutionary difference between L1 and L2. Nine kinds of antimicrobial resistance genes belonged to six categories were identified (Figure 1(g)). Gene qnrVC6 associated with quinolone resistance accounted for the highest positive rate at 12.6% (13/103), followed by trimethoprim resistance gene (dfrA6) and ß-lactamase resistance gene (bla_CARB-4), accounting for 5.8% (6/103) and 2.9% (3/103), respectively. Regardless of the category or quantity of drug resistance genes, the positive strains in L1 were significantly more than that in L2 (P < 0.05, Figure 1(g)), suggesting again that L1 strains is more closely related to human activities, as compared with L2 strains.

O-antigen biosynthesis gene clusters (O-AGCs) of V. metschnikovii

Due to the draft genomes, only 49 contiguous O-AGCs (located in the same contig) were obtained in this study. The 49 O-AGCs were grouped into 29 clearly distinguishable O-antigen genotypes based on gene organization (Figure 2), defined as VMOg1-VMOg29 respectively, including two sub-groups within VMOg10 (VMOg10-a and VMOg10-b) and four sub-groups within VMOg29 (VMOg29-a, VMOg29-b, VMOg29-c, and VMOg29-d). These four VMO29 subgroups had one or two hypothetical protein gains or losses between the orf19 to orf22 and the orf29 to orf35. Of all VMOgs including six sub-groups, 39.4% of them belonged to L1 (13/33), and the rest were in L2 (n = 20, 60.6%) (Figure 2). The putative functions of the identified ORFs in these V. metschnikovii VMOgs were listed in Tables S4.1–S4.29.

Figure 2. Genetic structures of the 29 V. metschnikovii O-genotypes identified from the V. metschnikovii reference strains. The different colourful arrows indicated the function, location and direction of gene transcripts in the O-AGCs. For putative glycosyltransferases, in green, types of domains families (family 1, 2, 4, UT, untypable) were indicated instead of gene names.

O-AGCs exhibited a high diversity due to variations in their arrangement and organization, but they mostly contained three major classes of genes, namely nucleotide sugar biosynthesis genes, glycosyltransferase genes, and O-unit processing genes. For nucleotide sugar biosynthesis genes, as an example, rmlBADC genes which are involved in converting glucopyranose-1-phosphate to dTDP-Rha in the deoxythymidine diphosphate (dTDP) sugar pathway, were present in VMOg5, VMOg7, VMOg16 and VMOg29. While eight genotypes (VMOg2, VMOg3, VMOg4, VMOg11, VMOg15, VMOg17, VMOg23, and VMOg24) only contained rmlBA genes, more than half (17/29) contained no nucleotide sugar biosynthesis genes or only a limited number of these genes. For putative glycosyltransferase genes, each genotype contained two to six genes encoding glycosyltransferases, with genes in families 2 and 4 being widely distributed. For O-unit processing genes, wzx/wzy in the Wzx/Wzy pathway or wzm/wzt in the ABC transporter pathway, were identified in V. metschnikovii O-AGCs. Most O-AGCs (97%, 28/29) were assembled by the Wzx/Wzy pathway (Figure S2a, S2b), while only VMOg29-a/b/c/d used the ABC transporter pathway to assemble the O-unit (Figure S2c, S2d). We also developed an in silico O-genotyping system for V. metschnikovii according to the sequence variation of the wzx gene (or wzt gene). Table S5 showed the VMOg distribution of V. metschnikovii with known O-AGCs sequences.

Carriage and distribution of virulence-associated genes in V. metschnikovii

Based on the annotations and predictions of the VFDB database, 161 putative virulence-associated genes involving 19 virulence factors in 9 virulence factor classes were found within the 103 V. metschnikovii genomes, including genes associated with adherence, chemotaxis, motility, secretion system, and toxins (Figure 3(a), Table S6). Almost all V. metschnikovii isolates contained virulence factors like flagella, cytolysin (hlyA), thermolabile hemolysin (tlh), and mannose sensitive hemagglutinin (MSHA) type IV pili. The type VI secretion system (T6SS) also appears to have evolved in V. metschnikovii populations. More than half of the strains in L1 carry T6SS-related genes. Interestingly, nine strains in L1 which belonged to VMOg29 and formed into two distinct evolutionary branches contained four genes related to the lateral flagella. In contrast, these nine strains of VMOg29 all appear to lack the ugd gene coincidentally, which is related to core tetrasaccharide synthesis in the process of capsular polysaccharide biosynthesis (Figure 3(a)). The RTX clusters (rtxA, rtxB, rtxC, rtxD), an important virulence factor in V. cholerae encoding the RTX toxin, were only present in four strains in L2. Also, only three strains in L1 carried the T3SS gene cluster (Figure 3(a)). Putative virulence-associated genes in V. metschnikovii were found to have homologs in other 23 foreign species (or genera). Among these species, V. cholerae, Haemophilus influenzae, Aeromonas hydrophila, Francisella tularensis, and Legionella pneumophila were the top five species sharing homologous virulence genes (Figure S3).

Figure 3. Cytotoxic and genomic-based analyses of the pathogenicity factors in V. metschnikovii. (a) Distribution of virulence-associated genes, O-genotypes, and cytotoxicity ranks. The colours from blue to red represented the identity of the gene sequences. The detail of the different cytotoxicity ranks and O-genotypes are listed on the right. (b) Cytotoxicity value of 30 V. metschnikovii strains. Low cytotoxicity (below 0.3), medium cytotoxicity (between 0.3 and 0.5) and high cytotoxicity (above 0.5). Vibrio cholerae strain N16961 (O1 biovar El Tor) and MO45 (O139) served as cytotoxicity positive controls. (c) Principal component analysis (PCA) of virulence genes for the 30 V. metschnikovii strains. The colours from blue to red represented the values of cos2, and the sizes of the circles represented the contributions of each strain sample. Gene lfhA (accession: VF0474) was described in the VFDB database as a gene encoding lateral flagellar biosynthesis protein in Aeromonas genus. The nomenclature of the flhA gene in V. parahaemolyticus was used for the sake of uniform description [55]. (d) Relationships among cytotoxicity levels, O-genotypes, and virulence factors, according to the binary logistic regression model used in this study. Virulence factors associated with the lateral flagella and the T6SS are marked with red dashed boxes.

Cytotoxicity features among V. metschnikovii strains

We performed cytotoxicity assays on a selected set of 30 Chinese V. metschnikovii strains collected as part of this study. By cytotoxicity, the 30 V. metschnikovii strains were divided into three groups, i.e. high-, median- and low- cytotoxicity (Figure 3(b)). Most strains were classified as median- and low- cytotoxicity (MC and LC). Seven strains were identified as high toxicity (HC), including strain ZJ20140528 isolated from a patient with diarrhea.

We further investigated the relationships between virulence genes, population structure, cytotoxicity levels as well as O-genotypes. By the principal component analysis (PCA), four lateral flagella genes (i.e. lfgG, lfiQ, flhA, lfiN) were the variables that contributed the most to the observed variance in the PCA, followed by gene vgrG and hcp (related to T6SS, Figure 3(c)).

Using the cytotoxicity data of the 30 V. metschnikovii strains, a Binary Logistic Regression (HC vs. non-HC) was performed to predict correlation between virulence genes and cytotoxicity. Lateral flagella genes were found to be an independent risk factor for HC (Odds ratio = 11.244; 95% CI = 2.725–30.647), together with the prevalence of the complete T6SS gene clusters (including hcp and vgrG, Odds ratio = 15.833; 95% CI = 1.175–28.843). A relatively good accuracy (sensitivity 85.7%, specificity 95.6%) was achieved in identifying whether a given V. metschnikovii strain belonging to HC using this model. Sankey diagram was then used to depict the correspondence of virulence factors and O-genotypes with cytotoxicity ranks (Figure 3(d)). Only three sub-groups within VMOg29 (VMOg29-b, VMOg29-c, and VMOg29-d) showed an increasing susceptibility to HC group, while the rest of the O-genotypes were more likely to be related to non-HC groups.

Lateral flagellar biosynthesis gene (flhA) is essential for the pathogenicity of V. metschnikovii

The flhA gene was detected in 9.7% (9/103) of V. metschnikovii strains, which all belonged to VMOg29 groups (Table S7). Early studies on lateral flagella were investigated in Aeromonas [36], and the flhA gene cluster of V. metschnikovii was similar to that of Aeromonas hydrophila (Figure 4(a)). On analysis of a flhA mutant (ΔflhA) constructed as part of this study, together with and a complemented flhA mutant (C-ΔflhA), ΔflhA was found to be significantly impaired in swimming motility, while C-ΔflhA reverted to similar swimming motility as the WT (Figure 4(b)). Phenotypic validation of four representative strains (belonging to LC, MC, HC, and control) revealed that the deletion of the flhA gene affected swimming motility (Figure 4(c,d)) and also weakened the biofilm formation ability of the V. metschnikovii strains (Figure 4(e,f)). Next, we examined the ability of V. metschnikovii to adhere to Hep-2 cells and cause infection. ΔflhA exhibited little cell infections relative to the WT and the C-ΔflhA strains, for example, a reduction in cytotoxicity (Figure 4(g)), cell adhesion (Figure 4(h)) and invasion (Figure 4(j)).

Figure 4. Genetic and functional analyses of lateral flagellar related flhA gene. (a) Gene cluster organization of a 15.8 kb fragment surrounding the flhA gene in the V. metschnikovii chromosome. (b) Swimming motility patterns of the ZJ20140528 (WT), the flhA mutant (ΔflhA), the complemented flhA mutant (C-ΔflhA) and control (2020QD6) strains. (c) Swimming motility of four representative V. metschnikovii strains for 12 h on soft agar (0.35%) and (d) statistical analysis. Strain ZJ20140528 for high cytotoxicity, VMOg29-d and flhA positive; strain ZJ060925 for median cytotoxicity, VMOg29-c and flhA positive; strain 2020QD3 for low cytotoxicity, VMOg-undefined and flhA negative; strain 2020QD6 for control, low cytotoxicity and flhA negative. (e) Crystal violet staining for 48 h and (f) statistical analysis. (g) Cytotoxicity values of five V. metschnikovii strains. (h) Adhesion ability in Hep-2 cells. (i) Invasion ability in Hep-2 cells. (j) Growth curves. The four time points (T1 to T4) for subsequent RT-qPCR were marked with black arrows. (k) Transcription levels of flhA gene at T2. The 16S rRNA gene was used as an internal control. All primer sequences were shown in Table S8. (l) Transcription levels of flhA gene at T3. Those groups with significant differences were marked with * (P < 0.05), ** (P < 0.01), *** (P < 0.001), **** (P < 0.0001).

It should be emphasized that the positive phenotypes observed in the strain ZJ060925 (flhA+) were weaker than that in WT (Figure 4(c–i)). Four different time points (T1 to T4) were selected based on the growth curves (Figure 4(j)), and then reverse transcription quantitative PCR (RT-qPCR) was performed to estimate the transcript levels of flhA gene in different strains at the same growth state. During the mid and late stage of the exponential phase (T2 and T3), the transcription of flhA in the strain ZJ060925 was higher than that of the ΔflhA but remarkably less than that of the WT and C-ΔflhA strains (P < 0.05, Figure 4(k–l)). The higher the expression of flhA, the stronger the swimming and cell infection activity of V. metschnikovii strains. The transcription of the flhA gene at the four time points is shown in the Figure S4. These data indicate that the lateral flagella are essential for the initial infection stage of chemotaxis and adhesion to the cell in the cell infection model.

Discussion

V. metschnikovii is a rare human pathogen. Previous studies have shown that V. metschnikovii belongs to the Cincinnatiensis clade [37]. Our comparative genome analysis of V. metschnikovii identified two major lineages with a high genetic diversity of O-genotypes. Our finding also revealed a pathogenic link among virulence-related factors, O-genotypes, and cellular infection models. This study provides an important basis for further study on the evolution and potential pathogenicity of V. metschnikovii.

We defined two major lineages of V. metschnikovii. The strains we collected in China were isolated from different regions in different years. The Mantel test and phylogenetic analyses strongly suggested that the spread of strains from L1 was likely to be closely related to human activities, while strains from L2 tended to be locally endemic as clonal complex groups. A range of important evidence was found: (1) L1 strains were sourced from food and clinical settings, while the origin of L2 strains was not closely related to humans. The sources of food often involve multiple epidemiological factors of human activity, such as food supply, packaging, transportation, and retail trading; (2) By analyzing six different mutation patterns, we found that both C > T/G > A and A > G/T > C were more frequent in L2 strains relative to L1 strains. The cumulative effect of C > T/G > A and A > G/T > C mutations is a continuous process that shapes bacterial ecological-niche adaptation to host immune environments [38]; (3) L1 strains had more drug resistance genes compared to L2, both in terms of category and quantity of resistance genes. A growing trend over antibiotic resistance in human-pathogenic bacteria was observed under environmental selective pressures [39], and part of the antibiotic-resistance rate in L1 strains can be attributed to the potential impact of frequent antibiotic utilization or even drug abuse in human communities. Of course, the real resistance of strains to those antibiotics is related to the presence or absence of resistance genes in addition to a variety of other factors, such as the capsule. Ernst et al. reported the global evolution of pathogenicity in carbapenem-resistant K. pneumoniae, and they demonstrated that mutations disrupting capsule biosynthesis genes impair capsule production, which enhances drug-tolerant K. pneumoniae to persist to yield potentially untreatable, persistent infection [40].

The O-AGCs based Vibrio serotyping became the “gold standard” for the clinical detection and epidemiological surveillance of human Vibrio pathogens [41]. We defined 29 O-genotypes, which were named VMOg1–VMOg29, respectively. But one concern is that only 61% (63/103) of the strains were matched in the in silico O-genotyping, leaving the VMOgs of other 40 strains undefined. Repeats has always been the source of most mis-assemblies, although the vast majority of genes in the bacterial genome are single-copy. Assembly errors from short reads or variations resulting from different assembly methods deserve our attention [42, 43]. And a potential mechanism for recombination and deletion events could contribute to the antigenic diversity [44]. Our collection of V. metschnikovii has a considerable O-antigen diversity, suggesting that the genetic differences of these strains are also representative. On the other hand, nearly 40% of the strains still have O-genotypes to be identified in future work. We hypothesize that the real O-antigen diversity of V. metschnikovii that exists in nature is far greater than we currently know. It would be interesting to know much more about the O-AGCs diversity of V. metschnikovii.

The pathogenicity of V. metschnikovii to humans is often overlooked. In the present work, the pathogenicity of V. metschnikovii may be the outcome of multiple factors like some opportunistic pathogenic vibrios [45]. Virulence genes (e.g. hlyA, tlh and MSHA) were commonly present, which have been proved to cause infection in other vibrios [46, 47]. Some strains carried virulence genes associated with chemotaxis and motility, secretion systems, and toxins. The contribution of lateral flagella (the first principal component, P1) to the cell infections can explain more observed variance. The absence of flhA reduces the swimming activity and biofilm formation capacity of the V. metschnikovii. And their motility and biofilm-forming capacity in the early stages of infection [48] is associated with the transcript levels of flhA. The mixed flagellation system was initially described in V. metschnikovii in the early 1950s [49]. V. parahaemolyticus and some other members of Vibrios also contain elaborate lateral flagella that are used for rapid movement on surfaces [50]. A previous study has suggested that the transition from the motility phase to the biofilm formation phase is one of the crucial events in infection [51]. Inconsistent with previous studies [18, 52], virulence genes related to cytolysin (or hemolysin) were not found to be statistically significant in differentiating cytotoxicity ranks.

Our data also suggested a link among evolutionary lineage, O-genotypes, and virulence factors. To begin with, the evolutionary tree based on cgSNPs and heatmap (Figure 3(a)) revealed virulence heterogeneity between the two lineages. The lateral flagella, T6SS and T3SS were mainly found in the L1 strains. While several V. metschnikovii strains in L2 carried the RTX toxin, which can play a pathogenic role in cell adhesion and/or pathogenesis [53] in V. cholerae [54]. Next, V. metschnikovii strains carrying different virulence genes exhibited different ranks in cytotoxicity and pathogenic phenotypes. These observations suggest that pathogenicity in humans can be linked with the presence of specific virulence-associated genes; however, further studies, with additional human isolates, are needed to confirm this hypothesis. Thirdly, we found that the cytotoxic rank was related to the O-genotype of V. metschnikovii. The HC group of V. metschnikovii contained VMOg29 specific virulence-associated genes, for example, the lateral flagella, whereas it was absent in the non- VMOg29 (MC or LC O-genotypes) strains. In this study, we propose a potential pathogenic mechanism for V. metschnikovii in the early stages of infection: the lateral flagella may play a fundamental and necessary role in the chemotaxis and host adhesion; other virulence factors could also promote pathogenicity, acting synergistically.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This work was supported by the National Science and Technology Fundamental Resources Investigation Program of China (2021FY100904).