The genomes used in this study and the main features of the corresponding strains are detailed in Table 1. The strains whose genome was sequenced are marked in Table 1. These strains were routinely grown in Tryptone Soy Broth or agar plus 5 g/l NaCl (TSB-1 or TSA-1, Pronadisa, Spain) at 28°C for 24 h. The strains were maintained both as lyophilized stocks and as frozen stocks at −80°C in marine broth (Difco) plus 20% (v/v) glycerol.

DNA was extracted using GenElute™ Bacterial Genomic DNA (Sigma, Spain) from bacteria grown with shaking at 28°C for 12 h. Samples with a DNA concentration of 10–15 ng/μl were used for sequencing with Illumina Genome Analyzer technology GAII (Illumina MiSeq) flow cell in the Genome Analysis Centre in Norwich (UK) and the SCSIE of the University of Valencia (Spain). To this end, unique index-tagged libraries for each sample (up to 96 strains) were created using TruSeq DNA Sample Preparation for subsequent cluster generation (Illumina cBot), and up to 12 separate libraries were sequenced in each of eight channels in Illumina Genome Analyser GAII cells with 100-base paired-end reads. The index-tag sequence information was used for downstream processing to assign reads to the individual samples (Harris et al., 2010).

The genomes of the Bt1 strain YJ016 (NC_005139 and NC_005140) and the plasmids pR99 (AM293858) and p4602-2 (AM293860) were used as templates for all the genomic analysis. The genome of the strain YJ016 was selected because it is one of the few genomes of V. vulnificus that has been accurately closed and annotated (Chen et al., 2003). The YJ016 genome contains 5,097 genes (3,387 in chromosome I and 1,710 in chromosome II), pR99 contains 71 genes and p4602-2 contains 67 genes.

Reads for each genome were done using SPAdes version 3.6.1 (Bankevich et al., 2012) with kmers from 21 to 127 and the careful option to reduce mismatches and short indels. Then, multiple sequence alignments were obtained by using Progressive Mauve software with default options (Darling et al., 2004). Locally Collinear Blocks (LCBs) with a size larger than 1 kb that were present in all the genomes were used to define the core genomes. The selected LCBs were concatenated to be used in subsequent analyses to build a core genome multiple alignment.

We selected seven Vibrio species that had at least one fully sequenced and annotated genome to define the core genome of the genus (CGG). Table S1 summarizes the characteristics of the closed Vibrio genomes selected for the study.

The identification of each gene was performed using local BLAST searches (Wang et al., 2003). Three independent searches, one per chromosome and plasmid, were performed, and a database was generated for each genome. The resulting sequences were mapped onto reference genes using Geneious 6.1.6 (Biomatters) software. As the sequences used for the generation of the databases were not closed genomes (draft genomes/contigs), only the genes present in all the genomes, with a minimum length of 80% of the total length of the homologous gene in the reference genome and with a DNA identity higher than 70%, were included in the CGS, the CGP or the CGG for further analysis.

In order to refine the preliminary alignment for the entire core, individual alignments for each gene were performed using the program MAFFT (Katoh and Standley, 2013). The elimination of unaligned ends and genes that did not match the above requirements was performed using an in-house Python script. Once all the genes were aligned, a concatenated sequence was generated for the CGS, the CGP and the CGG.

Functional annotation was performed using GPRO (Futami et al., 2011). Gene descriptions were obtained by blasting the predicted protein sequences against the NCBI non-redundant proteins (NR) and the Clusters of Orthologous Groups (COGs) databases (Tatusov, 2000) as reference subjects. Protein accessions obtained from the annotation based on NR were used to add additional Gene Ontology (GO) (Gene Ontology Consortium, 2008) and Enzyme Commission (EC) (Bairoch, 2000) designations. Information about metabolic pathways was retrieved via the web from the KEGG database (Kotera et al., 2012) based on EC numbers relation.

SNPs were identified in the genomes from the coding regions as described previously (Harris et al., 2010) with appropriate SNP cutoffs to minimize the number of false-positive/negative calls. SNPs were filtered to remove those at sites with a SNP quality score lower than 40 and/or with read coverage below 25X in this region. SNPs at sites with heterogeneous mappings were also filtered out if the SNP was present in fewer than 85% of reads for that position.

We identified human virulence-related genes in the CGS from previously published data to define the human virulence-related core genome (VCGS) (Chakrabarti et al., 1999; Chen et al., 2004; Horstman et al., 2004; Gulig et al., 2005; Bogard and Oliver, 2007; Alice et al., 2008; Oh et al., 2008; Brown and Gulig, 2009; Liu et al., 2009; Lee et al., 2010; Chen and Chung, 2011; Kim et al., 2011).

All the phylogenetic trees for CGG, CGS, CGP, and VCGS were reconstructed using the maximum-likelihood (ML) method with PhyML software (Guindon et al., 2009). The best evolutionary model for each dataset was determined with jModelTest (Posada, 2008) and considering the Akaike information criterion (AIC) (Akaike, 1974). The selected models were the generalized time-reversible model (GTR) with Gamma (+G) distribution and invariant positions (+I) (Tavaré, 1986) for SNP analysis. The pVvBt2 phylogeny was evaluated using the Hasegawa-Kishino-Yano model (Hasegawa et al., 1985). Support for the nodes derived in these reconstructions was evaluated by bootstrapping using 1,000 replicates (Felsenstein, 1985).

The congruence among the different phylogenetic reconstructions was evaluated using Shimodaira–Hasegawa (SH) (Shimodaira and Hasegawa, 1999) and expected-likelihood weight (ELW) tests as implemented in TreePuzzle version 5.2 (Schmidt et al., 2002; Strimmer and Rambaut, 2002).

We obtained almost complete genome sequences of 38 V. vulnificus strains (Table 1). The CGS was inferred from these genomes and additional genomes taken from the database (80 in total). The analysis included strains of the three Bts, isolated from a variety of hosts and habitats all over the world (Table 1). The main characteristics of the analyzed genomes and the CGS as well as the number of SNPs identified per chromosome are detailed in Tables 2, 3, respectively. The average gene identities in the CGS were 91 and 89% for chromosome I and II, respectively, which supports previous observations pertaining to the greater variability of chromosome II in Vibrio spp. (Kirkup et al., 2010).

The genes of the CGS and the associated metabolic pathways are shown in Tables S2, S3, respectively. The genes present in all the strains that did not match the criteria used to define the CGS (spanning less than 80% of the length and showing less than 70% identity with respect to the homologous gene in the reference YJ016 genome) are shown in Table S4. The CGS includes practically all of the genes for glycolysis, TCA cycle and pentose phosphate pathway, aerobic and anaerobic respiration, nitrate respiration, as well as for biosynthesis of metabolic intermediates, cofactors, nucleotides, amino acids and cell building blocks (Tables S2, S3). The CGS also includes genes involved in survival in different environments (water, animal intestine, chitinous surfaces), such as genes for resistance to different stressors (e.g., cold, heat, toxic oxygen forms, antimicrobials, and tellurite), genes for chitin degradation (i.e., various chitinases) and genes for surface colonization [e.g., pilus MSHA (mannose sensitive hemagglutination) and flagellum biogenesis]. It is generally assumed that genes on the Vibrio chromosome II have specific environmental functions related to habitat adaptation (Xu et al., 2003). Regarding the metabolic pathways, all of the common genes for general metabolic pathways such as butanoate metabolism, glycerolipid metabolism, peptidoglycan biosynthesis and pyrimidine metabolism are located on chromosome I while the genes for biosynthesis of siderophores, which could be related to habitat adaptation, are in the chromosome II but together with other metabolic genes not clearly related to habitat adaptation (i.e., fatty acid biosynthesis/elongation, glycerophospolipid metabolism, glyoxylate/dicarboxylate metabolism, nicotinate/nicotinamide metabolism, porphyrin metabolism, and valine/leucine/isoleucine metabolism) (Table S3). Regarding to survival genes, genes for the flagellum and MSHA pilus are located in chromosome I while those for the Flp/Tad pilus are located in chromosomes I and II (Table S2). In conclusion, there is no clear relationship between chromosome localization and habitat adaptation for the CGS-genes.

Figure 1 presents the distribution of gene ontology terms. Most of the genes in the CGS encode proteins associated with cell membranes that are related with regulation of transcription, transport or metabolic/oxido-reduction process and present catalytic/hydrolase, transferase or transcription-factor activity (Figure 1). Again, many functions associated to CGS genes in chromosome I were not found in chromosome II such as kinase activity, metal ion transport and phospho-relay response regulator activity, etc. (Figure 1). In addition, no CGS genes related to proteolysis, phospho-relay signal transduction, carbohydrate transport, phosphorylation, ATP catabolic process, chemotaxis, biosynthesis, DNA recombination and repair, and phosphoenolpyruvate-dependent sugar phosphotransferase system were identified in chromosome II (Figure 1).

To determine the CGP, the plasmid was reconstructed from the sequenced genomes of the Bt2 strains indicated in Table 1. According to the phylogenetic tree (Figure 3), we found six variants of the plasmid, two previously described (types II and IV) (Lee et al., 2008; Roig and Amaro, 2009) and four new ones (Table 4). The list of genes in the CGP is shown in Table S5. The CGP includes three virulence genes, two host-specific, vep07 and ftbp (fish transferrin binding protein), involved in resistance to and growth in eel serum, respectively (unpublished results, Pajuelo et al., 2015), and one host-nonspecific, rtxA1. This is a mosaic gene related with resistance to phagocytosis by murine and eel phagocytes (Lee et al., 2013a; Satchell, 2015) that presents at least seven different forms in V. vulnificus (Roig et al., 2011; Satchell, 2015).

To root the V. vulnificus phylogenetic tree within the genus, we first reconstructed a phylogenetic tree from the closed genomes of 15 strains belonging to seven Vibrio species (V. parahaemolyticus, V. anguillarum, V. algynolyticus, V. campbellii, V. furnissii, V. cholerae, and V. splendidus-clade) (main genome characteristics shown in Table S1) together with 27 selected V. vulnificus genomes (Table 1). Maximum-likelihood (ML) trees for chromosome I, II and I+II were reconstructed based on the common genes of these Vibrio spp. (Figure S1). The ML trees showed V. vulnificus as a very compact and independent group. Next, we inferred the phylogeny of the species from the 80 V. vulnificus genomes indicated in Table 1 by using ML reconstruction obtained from the SNPs of coding regions.

The phylogenetic trees for V. vulnificus are shown in Figure 2 and Figure S2. The rooted V. vulnificus phylogenetic reconstructions for chromosomes I and II were compared by congruence tests (Table 5). The results clearly indicate that chromosomes I and II were not congruent with one another, which strongly suggests that both chromosomes have suffered inter and/or intra chromosomal rearrangement since the emergence of the V. vulnificus ancestor. In fact, several strains changed their relationships in the tree for each chromosome (Figure 2). Among them, we highlight the well-known human clinical isolates C7184, YJ016, and CMCP6 (Figure 2). Despite the global incongruence between the two chromosomal phylogenies, all the phylogenetic trees divided the species into five well-defined, highly supported lineages, two of them (L4 and 5) including only a few strains. No strain changed lineage between the two trees (Figure 2).

The pairwise identity for the strains of all the lineages was 82.5% (Table 6). Lineage 1 (L1) (pairwise identity 93.1%) includes clinical (50% of isolates) and environmental (50% of isolates) Bt1 isolates from the USA, South Korea, Taiwan, Israel, and Spain. The clinical L1 isolates were recovered from human infections in the USA, South Korea, and Taiwan, with 80% of these derived from blood (related to primary septicemia, when the etiology is known), 10% from stools (related to gastroenteritis) and 10% from wound samples. The environmental L1 isolates were recovered from oysters (40%), water (30%), and fish (30%) in Spain, Israel, Taiwan, South Korea, and the USA. Interestingly, the strain yv158, although environmental, belongs to a previously described highly clonal group (clade A), which includes both environmental and clinical isolates (wound samples), all of vcg C-type (clinical type according to the polymorphism in virulence correlated gene; Rosche et al., 2005), supporting the potential virulence of this specific strain (Broza et al., 2012).

L2 (pairwise identity 87.4%) includes Bt1 and Bt2 strains from environmental sources (39.6% of the isolates; 21% from water, 21% from fish, 53.5% from seafood and 4.5% from sediment), diseased humans [29.2% of the isolates; 57% from human blood (none of them with known etiology) and 43% from wounds] and diseased fish (31.2%) origins. All of the zoonotic strains (Bt2-SerE), regardless of their source (environment, diseased animal or diseased human), country (France, Australia, Denmark and Spain) or year of isolation (from 1980 to 2004), cluster in a highly homogeneous group (pairwise identity 97.7%). The non-typeable Bt2 strain (960426-1 4C) formed a subgroup together with the SerE clonal-complex while SerI and A strains formed another subgroup both within L2 (Figure 2).

L3 includes all the Bt3 strains, which form a highly homogeneous group (pairwise identity 97.8%). These Bt3 strains were isolated in Israel from outbreaks of human vibriosis associated with the handling of farmed-tilapia (year of isolation, 1996–2003) and from aquaculture fishponds of tilapia (2002–2004). In consequence, all of them are associated to the aquaculture industry.

Interestingly, all human clinical cases of known etiology caused by L2 and L3 isolates in Europe and Israel were related to farmed-fish handling regardless the Bt of the isolate.

L4 (similarity 98.7%) is formed by two Spanish isolates from the Ebro-Delta (a nature park by the Mediterranean Sea) area, one from seawater and the other from a human clinical (leg wound) case, both of Bt1. Finally, L5 is formed by a unique isolate of Bt1 from Israel and clinical origin that is considered as representative of a highly virulent clone designated as Clade B (Raz et al., 2014; Efimov et al., 2015).

The phylogenetic tree for pVvBt2 based on the CGP shows relationships among isolates that did not match the previous phylogenetic relationships (Figures 2, 3). Further, the congruence analyses revealed that phylogenetic trees from the CGP and CGS were not congruent to each other, suggesting a different evolutionary history for the plasmid and the chromosomes (Table 5).

We searched for human virulence-related genes described in the literature for V. vulnificus, and we found that 75% of them were present in the CGS. This finding allowed us to define the human virulence-related core genome (VCGS) (Table S6). The VCGS includes genes for the flagellum, capsule, LPS- and cell-wall biosynthesis, motility, hemolysin, and proteases (such as VvhA and a collagenase), resistance to human serum (trkA; Chen et al., 2004), heme uptake, biosynthesis and uptake of vulnibactin, several genes for transcriptional regulators such as Fur, and finally genes for MARTX (Multifunctional Autoprocessive Repeat in Toxin) transport and modification systems, genes that are duplicated in pVvbt2 (Chakrabarti et al., 1999; Horstman et al., 2004; Yu-Chung et al., 2004; Gulig et al., 2005; Bogard and Oliver, 2007; Alice et al., 2008; Lee et al., 2008; Oh et al., 2008; Brown and Gulig, 2009; Liu et al., 2009; Chen and Chung, 2011; Kim et al., 2011; Oliver, 2015; Satchell, 2015). Remarkably, genes for several collagenases and chitinases were found in both chromosomes while the MARTX operon and vvhA were located on chromosome II.

The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.