In total, 518 isolates (from 14 countries, 1964 to 2018) were used in this study, including 325 newly sequenced isolates (29 clinical and 296 environmental) and 193 publicly available isolates (58 clinical, 124 environmental, and 11 unspecified). Detailed sampling information for the V. vulnificus isolates used in this study is listed in Supplementary Table S1 . The genomes of the publicly available isolates were downloaded from the NCBI database (https://www.ncbi.nlm.nih.gov/genome/browse#!/prokaryotes/189/).

All isolates were taken from the freezer (-80°C), streaked on Columbia blood agar (CBA) plates, and grown at 37°C for 24 hours. A single colony of a strain was transferred to 5 mL Luria-Bertani (LB) medium containing 2% NaCl and cultured at 37°C to the exponential phase with shaking.

DNA was extracted from samples using the QIAGEN UltraClean® Microbial DNA Isolation kit, Catalog no. 12224-50, as per the manufacturer’s instructions.

Sequencing was performed on the Illumina MiSeq platform with 150 bp paired-end read length. Raw reads with low quality were trimmed with the FASTQ Quality Filter (FASTX-Toolkit) (Pearson et al., 1997). De novo assembly of the filtered reads was performed using Shovill version 1.0.4 (Bankevich et al., 2012) with standard parameters. Sequenced isolates had an average genome size of 5.03 Mb and GC content of 46.71%.

SNPs identified from aligning the assemblies to the reference genome YJ016 with the NUCmer module in the software MUMmer 3.0 (Kurtz et al., 2004) were used to describe genetic relationships between isolates. SNPs located in repeated regions with low sequence quality (quality score <20 or covered by <10 reads) were excluded to identify SNPs in the core genome (regions presented in all isolates). After filtering, core-genome SNPs were then used in the maximum likelihood tree (MLTree) construction by RaxML (Stamatakis, 2006), and the MLTree was visualized with iTOL (Letunic and Bork, 2016).

All genomic sequences were annotated using Prokka (Seemann, 2014). We used GFF3 files generated by Prokka passed to Roary (Page et al., 2015) to create a pangenome and output gene presence/absence for each isolate. To obtain additional annotation, we used the pan-gene protein sequences of Roary to BLAST (BLASTP) against the COG and KEGG databases.

The 518 V. vulnificus isolates, including 87 clinical and 420 environmental isolates, were analyzed using the software Pyseer (Lees et al., 2018) and DBGWAS (Jaillard et al., 2018), using the isolate source (clinical or environmental) of V. vulnificus as a phenotype.

In Pyseer, the effects of SNPs, insertion and deletion of accessory genes, or k-mers on phenotype was evaluated, and the corresponding P-value calculated. SNP-based GWAS captures variants in the bacterial core genome. Insertion and deletion of accessory genes-based GWAS captures variants in the bacterial accessory genes insertions and deletions. Finally, k-mers are DNA sub-sequences of length k (typically 3-100 base pairs), and k-mers-based GWAS can reflect diverse genetic events, including SNPs, insertion and deletion of accessory genes, and cover the noncoding regions, including those related to transcriptional and translational regulation, overcoming limitations of analysis at the level of SNPs and non-core gene insertions/deletions (Jaillard et al., 2018).

DBGWAS is a k-mer-based GWAS approach that generates interpretable genetic variants linked to diverse phenotypes. Utilizing compacted De Bruijn graphs (cDBG), this method groups cDBG nodes pinpointed by the association model into subgraphs defined by their neighborhood in the initial cDBG. DBGWAS does not require prior annotation or reference genomes. Importantly, it is also computationally efficient.

Co-selection analysis was performed separately on the sequence alignment of 518 V. vulnificus isolates using the GWES tool SpydrPick (Pensar et al., 2019). SpydrPick hinges on mutual information (MI), which is a general measure of the dependence between two variants. Pairwise analysis of variants by SpydrPick can be performed using core-genome SNPs and pan-genome-wide analyses.

We performed whole genome variation detection of 518 V. vulnificus isolates from 14 countries between 1964 and 2018. Isolates included 193 isolates from the NCBI database – 58 clinical, 124 environmental, and 11 unspecified; and 325 isolates from newly collected samples – 29 clinical and 296 environmental ( Supplementary Table S1 ). From the analysis, 373,704 core-genome SNPs were identified, including 37,386 intergenic SNPs, 211,731 synonymous SNPs, 121,162 nonsynonymous SNPs, and 3,425 nonsense SNPs. The median pairwise SNP distance between all 518 isolates was 24,403, implying a remarkably high genetic diversity among V. vulnificus isolates.

A maximum likelihood (ML) tree of the 518 isolates based on core-genome SNPs with V. vulnificus YJ016 as the reference sequence divided the species into six well-defined lineages ( Figure 1 ). The isolates were mostly in lineages L1 and L2 (93.1%). L1 primarily contained isolates from China (86.8%), India, Mexico, and Bangladesh (100%), with environmental isolates accounting for a large proportion (85.9%). L2 primarily contained isolates from America (50.7%), Denmark, Australia, and France (100%), with environmental isolates accounting for a large proportion (75.4%). L3 only contained biotype 3 isolates from Israel. The number of L3, L4, and L5 isolates was small, but these lineages had a high proportion of clinical isolates (75%). L6 only contained three isolates, of which two were environmental isolates and one was unknown. The overlap in the geographic distributions of the six evolutionary lineages, for both environmental and clinical isolates was considerable, indicating lineages are not type- or geographic location-specific.

To explore novel genetic variants of V. vulnificus, we analyzed GWAS results of the 518 V. vulnificus isolates. We used Pyseer to identify variants in the whole genome sequences of V. vulnificus associated with clinical and environmental phenotypes based on SNPs, insertion and deletion of accessory genes, and k-mers, respectively. To reduce the false positive rate, we used a strict threshold to judge the analysis results in combination with prior knowledge (Lees et al., 2018), i.e., only genetic variation with a statistical test P-value below 1.74 x 10^-9 was included [^-log₁₀(P-value) >8.76; Figure 2 ]. We discovered 567 genetic variants related to the clinical phenotypes of V. vulnificus, located on 567 coding genes at the k-mers level. Twenty-eight genetic variants were identified at both the SNPs level and k-mers level, and eleven genetic variants were identified at both the insertion and deletion of accessory genes level and k-mers level. These variants involved genes that encode structural proteins, transport proteins, metabolic proteins, and signal proteins; thus, they can directly or indirectly affect the pathogenicity of V. vulnificus through their roles in energy production and conversion, cell cycle control, translation, ribosomal structure and biogenesis, replication, recombination, and repair, cell wall/membrane/envelope biogenesis, cell motility, intracellular transport, nucleotide/amino acid/carbohydrate/coenzyme/lipid transport and metabolism, inorganic ion transport and metabolism, post-translational modification, protein turnover, chaperone functions, and signal transduction mechanisms. Moreover, the purH (encoding bifunctional purine biosynthesis protein PurH) (Kim et al., 2003), and pldA (encoding phospholipase A1) (Koo et al., 2007) genes identified in this analysis have previously been experimentally verified as associated with the pathogenicity of V. vulnificus, thus confirming the effectiveness of our method.

We also used DBGWAS to test the association between k-mers and clinical vs. environmental phenotypes. Overall, we found 100 nodes related to the clinical phenotypes of V. vulnificus, located on 20 coding genes [P-value = 2.10 x 10^-7, -log₁₀(P-value) >6.68; Figure 3 ]. These variants also involved genes that encode structural proteins, transport proteins, metabolic proteins, and signal proteins. The gene products in this group are involved in energy production and conversion, cell membrane biogenesis, nucleotide/amino acid transport and metabolism, protein turnover, and signal transduction mechanisms. In addition, purH (Kim et al., 2003), and pldA (Koo et al., 2007) genes identified by Pyseer in this analysis were also identified by DBGWAS, which have previously been experimentally verified as associated with the pathogenicity of V. vulnificus, thus confirming the effectiveness of our method.

Based on Pyseer and DBGWAS, the clinical virulence or pathogenicity of V. vulnificus isolates is associated with thirteen genes, of which eleven (pglA, gmr, mbtB, tycC, yiaV, dsbD, hsdR, YJ016_00159, ramA, wbpA, and codA) were newly discovered in this study ( Supplementary Table S2 ).

Current GWAS methods are not accurate enough to identify causal variants, due to extensive genome-wide linkage disequilibrium and they ignore the effect of mutations on protein function (Chen and Shapiro, 2015). Therefore, we performed GWES to identify mutations throughout the genome that are likely to have co-evolved, i.e., potential epistatic interactions between genes.

We used SpydrPick on the whole genome sequences of 518 V. vulnificus isolates to assess potential epistatic interactions between genes. Most strong associations occurred on the chromosome between sites within 3 kb. We thus eliminated all sets of associations that spanned less than 3 kb, including those between core/accessory genome elements, to remove associations caused solely by physical linkage, which mask the co-adaptive signals we seek. After screening (extreme outlier threshold ≥0.28), 28,508 co-adaptation groups were identified that could be further subdivided into 166 co-adaptation networks containing 464 core-genome SNPs (53 intergenic SNPs, 221 synonymous SNPs, 177 nonsynonymous SNPs, and 13 nonsense SNPs) on 34 core genes and 750 accessory genes ( Supplementary Table S3 ; Supplementary Figure 1 ). The genes included in each co-adaptation group, as well as the corresponding annotation information, are available in Supplementary Table S4 .

The network of epistatic interactions was examined to identify potential V. vulnificus virulence genes. The high-frequency genes in the largest co-adaptation network (N1) were flgK, flgE, and flgL ( Supplementary Table S4 ). As a result, we considered flgK, flgL, and flgE genes to be relevant prospective targets of co-selection, i.e., potential V. vulnificus virulence genes. The genes flgK (Kim et al., 2008), flgL (Kim et al., 2008), and flgE (Lee et al., 2004) identified in our co-adaptation network have previously been shown (through deletion mutants) to affect the lethality of V. vulnificus in mice, lose the ability to form flagella, lose mobility, and exhibit serious defects in cell adhesion and biofilm formation, thus confirming the effectiveness of our method. Our GWES results combined with GWAS results showed that a total of 6 genes, namely purH, gmr, yiaV, dsbD, ramA, and wbpA, were associated with the pathogenicity of V. vulnificus clinical isolates ( Figure 4 ). The purH gene encodes PurH, which is involved in purine nucleotide biosynthesis and catalyzes the last two steps in the de novo biosynthesis of IMP (the first nucleotide in the de novo purine biosynthesis pathway). The purH gene has previously been shown (by deletion mutants) to alter the lethality and cytotoxic activity of V. vulnificus in mice and is a virulence gene of V. vulnificus (Kim et al., 2003). Five genes, namely purH, gmr, yiaV, dsbD, ramA, and wbpA, were newly discovered in this study. The gene gmr (encodes Cyclic di-GMP phosphodiesterase Gmr) regulates the enzyme for synthesis of cyclic di-GMP. Cyclic di-GMP has emerged as one of the most common and essential bacterial second messengers, with the ability to influence biofilm formation, motility, virulence, the cell cycle, differentiation, and other processes, hence influencing V. vulnificus pathogenicity (Römling and Amikam, 2006; Römling et al., 2013; Jenal et al., 2017). The gene yiaV (encodes inner membrane protein yiaV precursor) (Hu and Coates, 2005; Shimada et al., 2022) and dsbD (encodes disulfide interchange protein dsbD precursor) affect antimicrobial drug tolerance of Escherichia coli ( Missiakas et al., 1995). The gene ramA (encodes R-stereoselective amidase) affects antimicrobial drug tolerance of Klebsiella pneumoniae (Veleba and Schneiders, 2012). The gene wbpA encodes UDP-N-acetyl-D-glucosamine 6-dehydrogenase, which is involved in the synthesis of LPS. It has been previously demonstrated that the wbpA deletion mutant can affect the synthesis of LPS in Pseudomonas aeruginosa. LPS can evade host defenses, resist phagocytosis and serum-mediated killing and is also a crucial virulence factor for V. vulnificus (Burrows et al., 2000; Jones and Oliver, 2009).

The largest co-adaptation network (N1) included most of the interacting SNPs, which challenged interpretation due to its abundance of pairwise interactions. We thus used hierarchical clustering based on N1 variants to classify variants into two ecological groups (EGs), EG1 and EG2 ( Figure 5A ). EG1 contained L1 and L2 variants, while EG2 included variants within L1, L2, L3, L4, L5, and L6. EG1 isolates (1/59) had a lower percentage of clinical samples than EG2 isolates (86/459), indicating a lower virulence potential of this ecological group in humans. Despite the substantial number of coadaptation-related differences between EG1 and EG2, both ecological groups had no obstacle to gene exchange between groups in most regions of the genome ( Figure 5B ).

Given the complex structure of the co-adaptation network N1, we also did a cluster analysis of the coadaptation loci, which separated the coadapted complex loci into two tiers. Tier 1 loci variants were stable in EG1 but highly variable in EG2, containing genes csgD, glmM, glmU, mrdB, flgK, and flgL. Tier 2 loci variants broadly distinguished the two EGs, containing most of the significant coadapted loci (genes glmU, mrdB, flgK, flgL, flgE, spovD, alr, clpB, and glxR) identified by SpydrPick.

We detected 4 co-adaptation networks (N54, N78, N79, and N80) involving SNPs that were interactions of the type “incompatibility”, meaning that when one or more genes exist, at least other genes are absent. We detected a further 161 co-adaptation networks (N2-N53, N55-N77, and N81-N166) involving accessory genes that were “genome island-like”.

To assess the potential role of coadapted gene complexes on ecologically important phenotypic traits, we conducted a preliminary investigation of isolates from the two EGs (n = 5, 2 isolates for EG1, 3 isolates for EG2) to measure relevant phenotypes: growth rate under temperature challenges (4°C, 37°C, and 45°C), osmotic challenges (2, 4, 6, and 8% NaCl), and pH challenges (pH 4, 5, 6, 7, 8 and 9); motility assay, biofilm formation ability, survival under anaerobic conditions, biochemical parameters, and antimicrobial susceptibility. Neither EG1 nor EG2 isolates grew at 4°C or 45°C. EG2 isolates had a faster growth rate than EG1 isolates at 37°C. All V. vulnificus grew best at 2% NaCl concentration, but EG2 isolates were more tolerant to high salt environments than EG1 isolates ( Figure 6 ). The growth states of the two ecological groups were similar at acidic (pH 4-5) and alkaline (pH 9) conditions. EG2 isolates grew faster than EG1 isolates at pH 6-8. In addition, acidic conditions (pH 4-5) are not suitable for the growth and reproduction of V. vulnificus. Neutral conditions (pH 6-7) favor the growth and reproduction of V. vulnificus. Alkaline conditions (pH 8-9) not only impede the growth and reproduction of V. vulnificus but also inhibit its activity ( Figure 7 ). The growth variation of two ecological groups under temperature challenges, osmotic challenges and pH challenges occurred mainly during the stabilization and senescence phases of bacterial growth, with little change during the logarithmic growth period. In motility assay, both EG1 and EG2 isolates had strong swimming and swarming abilities, but there was no difference ( Figure 8 ). The biofilm formation ability of EG2 isolates was stronger than that of EG1 isolates ( Figure 8 ). As expected, all isolates failed to grow under strict anaerobic conditions ( Supplementary Figure 2 ). Biochemical assays detected no difference among isolates ( Supplementary Table S5 ). Minimal inhibitory concentration (MIC) values were employed to assess antimicrobial susceptibility. Isolates were sensitive to most of the antibiotics tested, but insensitive to cefazolin ( Supplementary Table S6 ). Cefazolin, as a first-generation cephalosporin, has limited effectiveness against Gram-negative bacteria, making them more susceptible to developing resistance. This is attributed to the relatively weak stability of cefazolin to β-lactamase produced by Gram-negative bacteria. Previous studies have reported the resistance of V. vulnificus to cefazolin (Wang et al., 2009; Pan et al., 2013).