Abstract
Vibrio vulnificus is an emergent marine pathogen and is the cause of a deadly septicemia. However, the evolution mechanism of antibiotic-resistant genes (ARGs) is still unclear. Twenty-two high-quality complete genomes of V. vulnificus were obtained and grouped into 16 clinical isolates and 6 environmental isolates. Genomic annotations found 23 ARG orthologous genes, among which 14 ARGs were shared by V. vulnificus and other Vibrio members. Furthermore, those ARGs were located in their chromosomes, rather than in the plasmids. Phylogenomic reconstruction based on single-copy orthologous protein sequences and ARG protein sequences revealed that clinical and environmental V. vulnificus isolates were in a scattered distribution. The calculation of non-synonymous and synonymous substitutions indicated that most of ARGs evolved under purifying selection with the Ka/Ks ratios lower than one, while h-ns, rsmA, and soxR in several clinical isolates evolved under the positive selection with Ka/Ks ratios >1. Our result indicated that V. vulnificus antibiotic-resistant armory was not only confined to clinical isolates, but to environmental ones as well and clinical isolates inclined to accumulate beneficial non-synonymous substitutions that could be retained to improve competitiveness.
Introduction
Vibrio vulnificus is an opportunistic pathogen for severe human infection, 1 and is one of the leading causes of non-Cholera, Vibrio-associated deaths globally. 2 Apart from clinical isolates, V. vulnificus inhabits a range of different environmental sources, including estuarine water, sediment and seafood produce. 3 Moreover. V. vulnificus isolates have been classified into 3 biotypes based on biochemical traits. 3 Biotype 1 as a human pathogen is the most common worldwide 4 and biotype 2 is regarded as an eel pathogen, 5 while biotype 3 is recorded as sharing biochemical properties of biotype 1 and 2. 6 The advances of sequencing and bioinformatic technology facilitate researchers to obtain V. vulnificus genomes. Recent studies elucidate genetic and evolutionary mechanisms of infections and pathogenesis of V. vulnificus at the genomic level.2,7,8
Antibiotic resistance in the V. vulnificus is a challenge that is associated with high morbidity and mortality. 9 Particularly, the presence of antibiotic-resistant genes (ARGs) in environmental isolates can be a huge risk to the public health. 10 However, the evolution mechanism of ARGs is still unclear. Moreover, the complete genomes can provide a more detailed gene information than the draft genomes. 11 In this study, the complete V. vulnificus genomes were obtained to demonstrate their ARGs distribution, phylogeny and evolutionary driving force, which could broaden our understandings of antibiotic-resistance in the V. vulnificus.
Materials and Methods
Collection of V. vulnificus and its relatives genomes
The complete genomes of V. vulnificus were obtained from the NCBI GenBank database. 12 In addition, other Vibrio type strain genomes were acquired from the gcType database 13 and their accession numbers in the NCBI GenBank database were listed in Table 1. Moreover, Escherichia coli ATCC 11775T was used as an outgroup in the phylogenomic analysis. Detailed information for the complete genomes in this study was shown in Table 1.
Genomic attributes of Vibrio genomes in this study.
Genomic quality estimation and annotations
The quality of obtained genomes was estimated by CheckM v1.10.3 with the typical workflow, and the genome with the completeness of >90% and contamination of <5%. were regarded as a high-quality genome as recommended by Bowers et al 14 Open reading frames (ORFs), rRNA and tRNA genes were annotated by using Prokka v.1.14.6 15 with the command “–gram neg –kingdom Bacteria –gcode 11.” Non-redundant 16S rRNA genes were generated by CD-HIT program v4.8.1 16 with the sequence identity of 99%. The classification of biotypes was analyzed based on the 16S rRNA gene phylogeny as described by Hoihuan et al 4 .
The ARGs were annotated against the Comprehensive Antibiotic Resistance Database (CARD) by using BLASTP with an e-value ⩽ 10, 5 with an identity and query coverage thresholds of 50% and 50%.17,18 In addition, those potential ARGs were double-checked by using the Resistance Gene Identifier (RGI) v.5.2.1 17 with the command. The potential of horizontal gene transfer was carried out by using BLASTP against the NCBI RefSeq select proteins database based on the best hit taxon.
Comparative genomic analysis and phylogenomic reconstruction
Based on the annotation, protein sequences translated from ORFs were compared by using OrthoFinder v.2.5.4 19 with the default setting. Single-copy were chosen into the following phylogenomic reconstruction as described previously. 20 In brief, each were aligned by MAFFT v.7.490 21 with the parameter “–auto.” Aligned protein sequences were refined by using trimAL v.10.2rev59 22 with the parameter “-automated1” to remove spurious sequences or poorly aligned regions and then concatenated. The maximum-likelihood phylogenomic analysis were carried out by using IQ-TREE v.1.6.12 23 with the parameter “-st AA -m LG+F+R10 -bb 1000,” and the best amino acid substitution model was also inferred by the same software with the parameter “-st AA -m MFP.”
Evolutionary analysis of ARGs
Protein sequences of each group of ARGs were aligned as described in the phylogenomic analysis. And the best amino acid substitution model for those sequences were also estimated by using IQ-TREE v.1.6.12. 23
The aligned codon sequences were generated based on converting a multiple sequence alignment of proteins and the corresponding DNA by using PAL2NAL v.14 24 And then, the resulting codon alignments were subjected to the calculation of synonymous and non-synonymous substitution (Ka/Ks) rates by KaKs_Calculator package v.1.2 25 through model selection and model averaging. Furthermore, the codon usages of those ARGs were summarized by using the sequence manipulation suite. 26
Statistics and visualization
The regression analysis were performed by using the formulation lm implanted in R version. 27 The phylogenetic trees were visualized by using MEGA7 software 28 and edited by PowerPoint 2019 software (Microsoft Cooperation, Redmond, WA, USA). Other figures were generated by and PowerPoint 2019 software (Microsoft Cooperation, Redmond, WA, USA).
Results and Discussion
General genomic features and phylogenomic relationship of V. vulnificus
A collection of twenty-four V. vulnificus isolates and other 18 Vibrio type strains genomes were obtained in this study (Table 1). Genomic quality estimations indicated that all genomes were high-quality with the completeness of >90.0% and contamination of <5.0% (Table 1), except for V. vulnificus FORC_036 (completeness of 100.0% and contamination of 5.1%), V. vulnificus FORC_053 (completeness of 100.0% and contamination of 5.4%) and V. cholerae CECT 514T (completeness of 55.3% and contamination of 0.8%). Based on isolation sources, twenty-two V. vulnificus isolates could be divided into 2 categories including clinical isolates(n = 16) isolated from hospital and patients,1,29-32 and environmental isolates (n = 6) cultivated from estuarine, Anguilla anguilla, Konosirus punctatus, Mya arenaria, Oreochromis sp. etc. 33 The neighbor-joining phylogenetic analysis based on 16S rRNA gene sequences revealed that those strains were classified into 2 clades including biotype 1 clade (strains 2142-77, 93U204, CMCP6, FDAARGOS_119, FORC_009, FORC_016, FORC_017, FORC_077, MO6-24/O, Vv180806, VV2014DJH, and YJ016) and biotype 2 clade (strains 06-2410, 07-2444, 2015AW-0208, 2497-87, ATCC 27562, CECT 4999, Env1, FDAARGOS_663, FORC_037, and FORC_054), as shown in Supplemental Figure S1.
Genomic sizes of 22 V. vulnificus isolates varied from 4.95 to 5.36 Mbp, while their DNA G + C contents were in a narrow range with 46.5% to 46.9% (Table 1). Commonly, their genomes constituted of 2 chromosomes, while several isolates genomes contained one plasmid as an accessory genetic material that was both present in the clinical (n = 3) and environmental (n = 4) isolates. Genomic annotation revealed those V. vulnificus isolates encoded 4447 to 5739 genes in their genomes. Their gene counts were mostly positively correlated with their genomic size (r2 = 0.96, P = 5.1e−15) except for the isolate 2015AW-0208, which harbored the most genes with the genomic size of 5.12 Mbp (Figure 1).

Genomic size and gene counts of Vibrio vulnificus. Red and blue indicated clinical and environmental isolates, respectively. Circle and triangle represented absence and presence of the plasmid, respectively.
Clinical and environmental V. vulnificus isolates were in a scattered distribution as indicated by the maximum-likelihood phylogenomic tree, and they were also clustered into 2 clades classified as biotype 1 and 2 clade (Figure 2), which was identical with the phylogenetic reconstruction based on 16S rRNA gene sequences. Moreover, the absence/presence of plasmid(s) did not affect their phylogenetic relationship either (Figure 2). Our phylogenomic reconstruction is similar with those of others about V. vulnificus.2,8,34 The recent phylogenomic reconstructions of V. vulnificus indicated that its isolates were divided into 4 or 5 clades, among which nearly 90% of isolates were clustered into 2 clades.2,34 The biotype 1 clade appeared to contain a significantly higher proportion of clinical isolates, while biotype 2 clade contained both of clinical and environmental isolates, that was also reported by López-Pérez et al 8 Despite the genomic divergence among clusters, a distinct pattern linking strain phylogeny, source of isolation, and virulent capabilities was not identified. This mixed distribution of clinical and environmental V. vulnificus isolates suggested that the genomic difference between 2 groups was subtle, which was also detected in other pathogens, such as Escherichia coli, 35 Legionella pneumophila, 36 Pseudomonas aeruginosa, 37 and Pseudomonas putida. 38 Moreover, this distribution indicated V. vulnificus antibiotic-resistant armory was not only confined to clinical isolates, but to environmental ones as well.

The maximum-likelihood phylogenetic tree based on single-copy orthologous protein sequences. Red and blue indicate clinical and environmental isolates, respectively. Escherichia coli ATCC 11775T was used as an outgroup.
Distribution of ARGs in the V. vulnificus
Twenty-three ARG orthologous proteins classified into 5 resistance mechanisms including antibiotic efflux, antibiotic inactivation, antibiotic target alteration, antibiotic target protection and antibiotic target replacement, were annotated in the genomes of V. vulnificus (Table 2). Besides, 2 ARGs annotated as ATP-dependent lipid A-core flippase MsbA and tetracycline efflux Na+/H+ antiporter family transporter Tet35 were compared into 2 different orthologous genes. Those ARGs were located in the V. vulnificus chromosomes, rather than in the plasmid, which was not common in other pathogens.39,40 Furthermore, core and accessory ARGs showed highest sequence identities with V. vulnificus or other Vibrio species (Supplemental Table S1), indicating that the low horizontal gene transfer frequencies of those ARGs. 41
Detailed information of antibiotic-resistant genes annotated in the V. vulnificus genomes.
Among those ARGs, 61% (14/23) of them shared by V. vulnificus and other Vibrio type strains were chloramphenicol acetyltransferase CatB9, translational regulators CRP and RsmA, dihydrofolate reductase Dfr-A3, fosfomycin resistance phosphotransferase FosC2, histone-like nucleoid structuring protein H-NS, ABC-type macrolide antibiotic exporter MacB, mobile colistin resistance phosphoethanolamine transferase MCR-9, ATP-dependent lipid A-core flippase MsbA, polymyxin resistance phosphoethanolamine transferase PmrE, quinolone resistance protein QnrVC1, redox-sensitive transcriptional activator SoxR, tetracycline efflux Na+/H+ antiporter family transporter Tet35 and AcrAB-TolC multidrug efflux pump YajC (Figure 3). Furthermore, the metallo-beta-lactamase VarG were exclusively in all of V. vulnificus isolates, other than the isolate VV2014DJH. Other exclusive genes encoding the response regulator ArlR, multidrug efflux RND transporter permease subunit MexK, ATP-dependent lipid A-core flippase MsbA, H+-coupled multidrug efflux pump PmpM, rifampin-resistant beta-subunit of RNA polymerase RpoB2, tetracycline efflux Na+/H+ antiporter family transporter Tet35, tetracycline-resistant ribosomal protection proteins TetT and TetWNW were mostly in the isoalte 2015AW-0208, which showed genomic differences compared with other V. vulnificus isolates.

Distribution of antibiotic-resistant genes in the V. vulnificus genomes.
Antibiotic resistance determinations revealed that V. vulnificus could resist various antibiotics including fluoroquinolones, β-lactam combination agent, lipopeptide, macrolide, nitrofurans, penicillin and phenicols,63-66 which were consistent with those ARGs found in their genomes (Table 2). Furthermore, the genomes of environmental isolates 93U204, ATCC 27562T, CECT 4999, Env1, FORC_037, and FORC_054 encoded 13 or 14 ARG orthologous genes, which were similar with the ARGs profile of clinical isolates (Figure 3). Therefore, environmental V. vulnificus isolates may serve as reservoirs for transmission of their antibiotic resistance.
Evolution of ARGs in the V. vulnificus
Except for FosC2 only annotated in the V. vulnificus VV2014DJH and V. campbellii ATCC 25920T, other 13 shared ARG orthologous proteins and VarG were processed into phylogenetic analysis. V. vulnificus isolates were clustered into an independent clade, which was separated from other Vibrio strains, of phylogenetic trees based on most ARG orthologous proteins including CatB9, CRP, Dfr-A3, H-NS, MacB, MCR-9, MsbA, QnrVC1, SoxR, Tet35, and YajC (Supplemental Figure S1). However, V. vulnificus isolates were in a scattered distribution in the phylogenetic trees of PmrE and RsmA, that were classified into antibiotic target alteration or antibiotic efflux resistance mechanisms. Clinical and environmental isolates in those scattered distribution phylogenetic trees were still mixed (Supplemental Figure S1), demonstrating that environmental ones had a similar evolutionary history with clinical ones which were also observed in other pathogens.67,68
The calculation of non-synonymous and synonymous substitutions indicated that most of ARGs evolved under purifying selection with the Ka/Ks ratios lower than one 69 (Figure 4). While several pairwise gene sequences of ARGs had the Ka/Ks ratios higher than 1 indicating that those genes evolved under positive selection 70 (Figure 4), especially for those antibiotic efflux ARG proteins H-NS, RsmA, and SoxR among which the group of Ka/Ks > 1 accounted for 6.5% to 10.0%. And those high Ka/Ks ratios were mostly confined to isolates 2015AW-0208, FDAARGOS_119, FDAARGOS_663, VV2014DJH, and YJ016, which were all isolated from clinical sources. Compared with environmental sources, clinical settings, where most antibiotics are prescribed, are hypothesized to serve as a major hotspot, 71 that forced several beneficial non-synonymous substitutions could be retained to improve competitiveness.

Ka/Ks ratios of each orthologous ARG in the V. vulnificus.
Conclusions
With regard to uptake of antibiotic resistance factors, marine environments with highly variable ecological niches provide an unrivaled gene pool with a diversity that considerably exceeds that of the human and marine animal microbiota. Indeed, the most remarkable feature of marine microbiome is its enormous diversity, providing numerous genes that potentially could be acquired and used by pathogens to counteract the effect of antibiotics. 72 Moreover, metals and antibiotic pollutions co-selecting for antibiotic-resistant strains via cross-resistance or co-resistance should also be taken into consideration seriously to retard the rapid evolutionary expansion and spread of antibiotic resistance factors. 73 Clinical and environmental V. vulnificus isolates were in a scatter distribution based on the genomic size and constitutes as well as the phylogenomic relationship. Genomic annotation and comparative genomic analysis also indicated that this mixed distribution of ARGs in clinical and environmental isolates. Unexpectedly, those ARGs were located in their chromosomes, rather than in the plasmids of them, suggesting that those genes were conserved in the V. vulnificus. The calculation of non-synonymous and synonymous substitutions indicated that most of ARGs evolved under purifying selection with the Ka/Ks ratios lower than one, while h-ns, rsmA and soxR in several clinical isolates evolved under the positive selection with Ka/Ks ratios >1. Therefore, V. vulnificus antibiotic-resistant armory was not only confined to clinical isolates, but to environmental ones as well and clinical isolates inclined to accumulate beneficial non-synonymous substitutions that could be retained to improve their competitiveness.
Supplemental Material
sj-docx-1-evb-10.1177_11769343221134400 – Supplemental material for Distribution, Phylogeny and Evolution of Clinical and Environmental Vibrio vulnificus Antibiotic-Resistant Genes
Supplemental material, sj-docx-1-evb-10.1177_11769343221134400 for Distribution, Phylogeny and Evolution of Clinical and Environmental Vibrio vulnificus Antibiotic-Resistant Genes by Nan Geng, Guojin Sun, Wen-Jia Liu, Bin-Cheng Gao, Cong Sun, Cundong Xu, Ertian Hua and Lin Xu in Evolutionary Bioinformatics
Footnotes
Funding:
The author(s) disclosed receipt of the following financial support for the research, authorship, and/or publication of this article: This study was funded by the Natural Science Foundation of Zhejiang Province (LQ19C010006), the Key Technology Research and Development Program of Zhejiang (2021C03019, 2021C03196), the Zhejiang Basic Public Welfare Research Project (LGF19E090001), and the National Science and Technology Fundamental Resources Investigation Program of China (2019FY100700, 2021FY100900).
Declaration of Conflicting Interests:
The author(s) declared no potential conflicts of interest with respect to the research, authorship, and/or publication of this article.
Supplemental Material
Supplemental material for this article is available online.
References
Supplementary Material
Please find the following supplemental material available below.
For Open Access articles published under a Creative Commons License, all supplemental material carries the same license as the article it is associated with.
For non-Open Access articles published, all supplemental material carries a non-exclusive license, and permission requests for re-use of supplemental material or any part of supplemental material shall be sent directly to the copyright owner as specified in the copyright notice associated with the article.
