Reclassification of former V. toranzoniae R17 as V. kanaloae R17
According to genome sequence similarity and genomic indexes, the genomes of V. toranzoniae strains separated in two well defined clusters: on the one hand, the six strains isolated from clams and seawater in Europe, and on the other hand, the three strains isolated in Chile together with V. kanaloae (strains CCUG 56968T and 5S149)(Tables 2; Figs. 1, 2), a Vibrio species that was first isolated from diseased oyster (Ostrea edulis) larvae in France (Thompson et al., 2003a). Our results confirmed also that the three Chilean isolates are clones, with a dDDH value of 100%, and an OrthoANI value of 99.99–100% (Table 2, Fig. 1). Besides, OrthoANI and dDDH results showed that the chilean isolates are in fact V. kanaloae. OrthoANI and dDDH values between these isolates and V. toranzoniae strains were below the cut-off values proposed for the delineation of new species (< 96% and < 70%, respectively) (Konstantinidis and Tiedje, 2005; Goris et al., 2007). On the contrary, values for these genomic indexes compared to type strain V. kanaloae CCUG 56968T, were higher than 98.0% and 86%, respectively. Accordingly, the core-genome -based phylogenetic tree (Fig. 2) reinforced the existence of two separate monophyletic branches. Thus, based on these results, we proposed the assignation of Chilean isolates to V. kanaloae.
Revising the genetic sequences available at NCBI, we discovered that one of the two sequences deposited as the 16S rRNA gene of V. kanaloae type strain LMG 20539T was poorly named. Therefore, the 16S rRNA gene sequence with accession number AJ316193 (Thompson et al., 2001) coincided with V. kanaloae with 100% of similarity, followed by V. toranzoniae with 99.66%. Conversely, the other 16S rRNA gene sequence available, with accession number AM162657 (deposited by Le Chevalier, P. et al., unpublished), corresponded to V. atlanticus (99.93% of similarity), followed by V. tasmaniensis (99.86%), V. lentus (99.78%) and then V. toranzoniae (98.78%) and V. kanaloae (98.77%).
The wrong sequence AM162657 was deposited in 2005, when the second most similar species V. tasmaniensis had been already described (Thompson et al., 2003b). In addition, an identical sequence to AM162657 was submitted in 2011, with accession number NR_042468 and processed by NCBI staff, when both V. tasmaniensis and V. atlanticus 16S rRNA gene sequences were available. For the latter, 16S rRNA gene sequence was deposited in 2007, with accession number EF599163 (Beaz-Hidalgo et al., 2008). This last mislabelled sequence (AM162657) was uploaded by the National Center for Biotechnology Information for its NCBI RefSeq Targeted Loci Project, which includes curated RefSeq records and selected validated GenBank sequences for curated BLAST databases.
It has been highlighted previously that sequences wrongly deposited as type strains may lead to errors in further studies that depend on public databases. That was the case for the so-called Lelliottia nimipressuralis type strain SGAir0187 (Heinle et al., 2018), that was misclassified due to a false type strain and was not a strain of the species (Salvà-Serra et al., 2019). Also, Beaz-Hidalgo and coworkers (2015) detected at least 12 misidentified Aeromonas genomes among the 44 deposited at the NCBI, insisting these authors in the need of measures to prevent this kind of chaining errors.
In our case, the bad deposit of sequences led to the misassociation of the Chilean isolates with V. toranzoniae rather than V. kanaloae (Lasa et al., 2015). Considering all the results together and to avoid future problems, we have updated the taxonomic assignation of strain R17 and its deposited sequence (accession number GCA-001995825.2) to V. kanaloae.
Genomic Indices
The genome size of the V. toranzoniae strains studied ranged from 4.3 to 4.7 Mb, being 4.5 Mb the average size of the species (Table 3). This genome size is in accordance with the expected for a species of Vibrio genus (Thompson et al., 2009). A minimum of 3,826 and maximum of 5,184 coding sequences were predicted using RAST annotation server for the different strains. For RNAs amount, the number oscillated between 126 and 188.
Table 2
Values of DDH among V. toranzoniae strains.
| CECT 7225T | CMJ 9.4 | CMJ 9.11 | Cmf 13.9 | 96–373 | 96–376 | R17 | R18 | R19 | V. kanaloae CCUG 56968T |
CECT 7225T | 100.0 | | | | | | | | | |
CMJ 9.4 | 81.2 | 100.0 | | | | | | | | |
CMJ 9.11 | 77.9 | 78.5 | 100.0 | | | | | | | |
Cmf 13.9 | 80.7 | 79.6 | 77.1 | 100.0 | | | | | | |
96–373 | 82.3 | 80.9 | 77.5 | 80.30 | 100.0 | | | | | |
96–376 | 80.6 | 80.4 | 77.1 | 78.90 | 80.00 | 100.0 | | | | |
R17 | 58.7 | 58.9 | 60.6 | 59.5 | 58.6 | 58.6 | 100.0 | | | |
R18 | 58.5 | 58.9 | 60.6 | 59.4 | 58.5 | 58.5 | 100.0 | 100.0 | | |
R19 | 58.5 | 58.9 | 60.6 | 59.4 | 58.5 | 58.5 | 100.0 | 100.0 | 100.0 | |
V. kanaloae CCUG 56968T | 61.6 | 62.0 | 61.0 | 62.5 | 61.8 | 61.8 | 86.4 | 86.3 | 86.3 | 100.0 |
Table 3
Genome statistics for V. toranzoniae strains.
| CECT 7225T | CMJ 9.4 | CMJ 9.11 | Cmf 13.9 | 96–373 | 96–376 |
Genome size (Mb) | 4.61 | 4.76 | 4.56 | 4.71 | 4.64 | 4.64 |
G + C content | 44 | 43.9 | 43.8 | 43.9 | 43.9 | 43.9 |
Number of contigs | 2 | 299 | 311 | 192 | 132 | 135 |
Coding sequences | 4164 | 4236 | 4053 | 4158 | 4258 | 3826 |
RNA genes | 172 | 158 | 181 | 181 | 126 | 146 |
G + C content was practically the same among the strains, varying from 43.8 to 44 mol%, in the range for Vibrio species (Table 3). Values of OrthoANI among V. toranzoniae ranged from 94.73 to 100% (Fig. 1) and dDDH oscillated between 58.50 and 100% (Table 2).
Although several studies reported that genome size and G + C content show a correlation with the ecological strategies of marine bacteria (Giovannoni et al., 2014; Luo and Moran, 2015), our results did not show differences between the free-living bacteria (V. toranzoniae 96–373 and 96–376) and those associated to a host (V. toranzoniae CECT 7225T, CMJ 9.4, CMJ 9.11, Cmf 13.9).
Complete genome sequencing of type strain Vibrio toranzoniae CECT 7225T
V. toranzoniae was first described based on four isolates from cultured clams in Galicia (NW Spain), designating the strain CECT 7225T as the type strain of the species (Lasa et al., 2013, 2016). For this strain, the read depth obtained sequencing by PacBio technology was 92x, determining the genome size in 4,605,941 bp assembled in two contigs, consistently with the possession of two chromosomes by many species of the Vibrio genus, one larger and one smaller of approximately 3.2 and 1.4 Mb, respectively. G + C content was 44 mol% and no plasmids were identified.
The complementation between short Illumina and long PacBio reads did not significantly improved the genome assembly compared to PacBio-only sequencing, since the corrected genome size was 4,605,997 bp, only 56 bp longer than Pac-Bio-only sequenced genome.
With respect to strain 96–376, we were unable to close the genome, and the complementation between Illumina and PacBio reads yielded 6 contigs with a total genome size of 4,370,366 bp, that is, 31,016 bp less than PacBio assembly.
Core-genome analysis of V. toranzoniae
Core-genome analysis of the V. toranzoniae strains with GET_Homologues revealed a pan genome of 6,287 gene clusters, that is, shared by the six isolates included in the study. Of these 6,287 genes, 3,404 were shared by 5 isolates or more (soft core), 2,489 genes were only present in 2 or less taxa (cloud genome), 395 genes were remaining genes shared by several taxa (shell genome) and 2,953 genes were showed by all strains studied (core genome) (Fig. 3A). As seen in Fig. 3B, core genome moderately decreases when more genomes are included, whilst the pan genome acts in reverse. Phylogenomic tree based on pangenomic matrix of V. toranzoniae strains is represented in Fig. 4. Using Roary software, from the total pan genome of 6,132 genes, 2,766 belonged to the core genome (shared by 99–100% of taxa), whereas 3,366 formed the shell genome (shared by 15–95% of taxa).
The phylogenomic analysis of the core genome (Fig. 2) did not reveal a differentiation between strains according to the lifestyle either. However, when looking at the phylogenomic tree (Fig. 4), we observed a site-specific differentiation of the three strains isolated from clams in Camariñas (Galicia, Spain), thus sharing the same growing area. Since the pangenome comprises more genes, including those not shared by all strains, this could indicate a local episode of horizontal gene transfer. Consequently, geographical conditions appear to be more decisive than lifestyle or host in V. toranzoniae strains. Further studies are needed to confirm such hypothesis.
Genomic Features
Despite it was not observed a phylogenetic divergence between strains with different lifestyles, some notable differences in gene content were observed.
All strains except one, the environmental strain 96–376, showed the presence of the genes related to flagellar synthesis and regulation. Absence of motility in 96–376 strain was similarly observed in soft agar and checked by optical microscopy. The rest of strains exhibited motility in both soft agar and optical microscopy, being stained flagella observed in bacterial preparations at 100x optical microscopy (data not shown).
Likewise, the environmental strain 96–376, together with the other strain isolated from seawater 96–373, did not exhibit the genes for the rhamnose synthesis pathway, involved in the synthesis of the capsule. This biosynthetic pathway is common and highly preserved across both Gram-positive and Gram-negative bacteria, involving four distinct enzymes that transform glucose into dTDP-L-rhamnose. The initial enzyme in this pathway, glucose-1-phosphate thymidylyltransferase, is responsible for attaching a thymidylmonophosphate nucleotide to Glu-1-P. The resulting dTDP-glucose is further oxidated an d dehydrated by the enzyme dTDP-d-glucose 4,6-dehydratase. Subsequently, a third enzyme, dTDP-6-deoxy-d-xylo-4-hexulose 3,5-epimerase, facilitates a double epimerisation at the C3 and C5 positions. In the final step, the dTDP-6-deoxy-l-lyxo-4-hexulose reductase reduces the C4 keto group to produce the final product, dTDP-l-rhamnose. On the other hand, the type strain CECT 7225T lacked the reductase gene in this dTDP-rhamnose pathway, and the strain Cmf 13.9 was the only one hosting the thymidylyltransferase.
Presence of capsule was also assessed by growing in CRA plates. After 48 h of incubation, all the strains presented black colonies indicating the production of capsule although, according to the absence of rhamnose-synthesis pathway, strains 96–373 and 96–376 showed the lowest production, indicating also that rhamnose is important but not exclusive for capsule production. The presence of capsule in all the strains could be explained by the advantages that the extracellular polysaccharides confer for environmental survival, but also for host invasion, colonization, persistence and eventually pathogenesis (Bian et al., 2021). Contrary to what it was initially thought, capsule provide protection from physical and chemical stresses without detriment of a high transference of genetic materials between bacteria (Rendueles et al., 2018).
All the strains showed the presence of genes for the transport of iron and for the siderophore aerobactin, despite aerobactin synthase protein IucC was only present in strain 96–373. Nevertheless, only the strain Cmf 13.9 contained the kit of genes for the siderophore assembly. According to this, Cmf 13.9 was the only isolate capable to form orange halo around blue around the colonies in CAS plates, which is indicative of siderophore production.
Related to virulence factors (Table 4), all the strains hosted a vibriolysin and a haemolysin, putative for the case of CMJ 9.11. Besides, all the strains exhibited T1SS secreted agglutinin RTX toxin proteins, witj the exception of the type strain CECT 7225T. The strains also manifested the presence of the related Ca2+ binding proteins and the type I secretion system, and components necessary for the extracellular secretion, such as a TolC outer membrane protein, an ATP-binding cassette (ABC) and a LapC membrane fusion. Despite the presence of vibriolysins and haemolysins, which in other vibrios have been described as virulence factors (Yuan et al., 2020; Galvis et al., 2021), none of the strains of V. toranzoniae cause mortality for clam or turbot (data not shown). This led us to speculate that vibriolysins might not be expressed or that some of the regulation factors are absent. These observations suggest that the V. kanaloae strain R17 (reclassified in this work) isolated from moribund red conger eel in Chile could have been the responsible etiological agent, so that V. toranzoniae would remain only as a potential pathogen.
Table 4
Summary of genetic traits present in V. toranzonaie strains.
| CRISPR sequences | Incomplete prophage sequences | Secondary metabolites | Virulence factors | Phage defense elemente | |
CECT 7225T | 1 | 1 | PUFAs, ectoine, arylpolyene, bacteriocine, betalactone | Hemolysin, VirK, virulence-associated E family protein, Iron-regulated protein IrgB | RM (3), Druantia, Zorya, Cas, dGTPase, Viperin |
CMJ 9.4 | 3 | 2 | PUFAs, ectoine, arylpolyene, bacteriocine, betalactone | Hemolysin, Iron-regulated protein IrgB | RM (2), Cas (2), dGTPase, BstA |
CMJ 9.11 | 1 | 3 | PUFAs, ectoine, arylpolyene, bacteriocine, betalactone | Hemolysin, VirK | RM (2), dGTPase, BREX, DTR, Cas, Rst-sirtuin-like |
Cmf 13.9 | 2 | 1 | PUFAs, ectoine, arylpolyene, bacteriocine, betalactone, siderophore | Hemolysin, probable RTX,Iron-regulated protein IrgB, siderophore assembly kit | RM (2), dGTPase, Cas, Rst-sirtuin-like |
96–373 | 1 | 1 | PUFAs, ectoine, arylpolyene, bacteriocine, betalactone, siderophore | Hemolysin, Iron-regulated protein IrgB | RM (3), dGTPase, Septu, Hachiman |
96–376 | 1 | 0 | PUFAs, ectoine, arylpolyene, bacteriocine, betalactone | Hemolysin, Iron-regulated protein IrgB | RM (2), Nhi, Zorya, Kiwa, dGTPase, Rst-ATPase |
All strains harbour the antimicrobial peptides adeF, CRP, QnrS2, drfA6 |
Genomic differences between closely related strains are usually concentrated in strain-specific regions of the chromosomes known as genomic islands, that are generally acquired by HGT and that contain adaptive traits that can be linked to niche adaptation (Dobrindt 2004, Penn 2009). Using IslandViewer 4, Genomic Islands (GIs) were identified by SIGI-HMM and Island-Path-DIMOB methods, but not by IslandPick method (Table 5). For all the strains, the highest number of GIs was found by the SIGI-HMM method. The strain showing the highest GIs number was CMJ 9.4. Mobile elements, phage proteins, glycosyltransferases, lipid metabolism proteins and hypothetical proteins were the most found proteins within identified GIs. Iron acquisition system proteins, L-ectoine synthase or MSHA pilin proteins were also found.
Table 5
Number of identified GIs in V. toranzoniae strains.
| CECT 7225T | CMJ 9.4 | CMJ 9.11 | Cmf 13.9 | 96–373 | 96–376 |
| S | I | T | S | I | T | S | I | T | S | I | T | S | I | T | S | I | T | |
V. anguillarum | 29 | 8 | 37 | 33 | 12 | 45 | 25 | 13 | 38 | 26 | 12 | 38 | 29 | 10 | 39 | 26 | 12 | 38 | |
V. splendidus | 27 | 9 | 36 | 34 | 10 | 44 | 28 | 15 | 43 | 26 | 9 | 35 | 30 | 12 | 42 | 26 | 10 | 36 | |
V. vulnificus | 30 | 9 | 39 | 34 | 9 | 43 | 26 | 14 | 40 | 25 | 10 | 35 | 28 | 12 | 40 | 27 | 10 | 37 | |
S = SIGI-HMM method, I = IslandPath-DIMOB method, T = total. |
Horizontal Gene Transfer evidences
The evidenced high gene transfer was assessed by different indicators. For example, the abundance of secondary metabolites, which is indicative of genomic exchange since many of them are acquired by horizontal gene transfer (Khaldi et al., 2008). A total of six secondary metabolites were identified using AntiSMASH (Table 4). From them, five were distributed in all strains (polyunsaturated fatty-acid (PUFA) cluster, ectoine, bacteriocin, arylpolyene and betalactone). These secondary metabolites are related with the adaptation of the bacteria to marine environment (Jensen 1996, de Carvalho 2010), found in different marine bacterial genus. Thus, PUFAs are produced by different marine bacteria such as Vibrio, Photobacterium, Psychromonas or Shewanella, enabling the transportation of nutrients through the membrane and maintaining its fluidity in the deep-sea cold environment in which these genera inhabit (Moi et al., 2018). Aryl polyenes were described as natural bacterial products which protect bacteria from reactive oxygen species (Schöner et al., 2016). Also, bacteriocin and betalactone are compounds produced by bacteria which show inhibitory or killing activities against other cells (Manivasagan et al., 2014; Yang et al., 2014). Finally, ectoine is an organic compound whose accumulation within the cell allows bacteria to keep turgor pressure under high osmolarity, thus proportioning the cell resistance against saline stress (Gregory et al., 2019). Here, we found genes coding for ectoine synthesis in genomic islands, which are usually enriched in secondary metabolites genes, providing evidence that secondary metabolism is linked to functional adaptation (Penn 2009). Also, a siderophore cluster was only recognized in two strains, namely Cmf 13.9 and 96–373, consistent with what we observed in the genome browser.
Antiphage systems, whose variable possession in closely related strains, as it is our case, indicate high rate of horizontal gene transfer (Tesson et al., 2022). Those systems were identified for all the strains (Table 4), in a number from five to eight, as the average number for prokaryotic genomes which is 5 (Tesson et al., 2022). Among them, all the strains encode for RM and dGTPase, the most common antiphages systems together with Cas which, interestingly, is only present in the strains isolated from clams.
All the strains presented CRISPR sequences, in a variable number from 1 to 3 (Table 4). For the case of Cas cluster gene sequences, only the strains CMJ 9.4 and CMJ 9.11 presented 2 and 1, respectively. None intact prophage sequence was detected, although he majority of the strains showed 1 to 3 incomplete prophage sequences (Table 4). Only the environmental strain 96–376 did not present any prophage sequence, neither intact or incomplete or questionable.
Using CARD, four antibiotic resistance gene sequences were identified for all the strains (Table 4), coding for quinolone resistance protein QnrS2, two resistance-nodulation-cell division antibiotic efflux pumps (adeF and CRP) and a trimethoprim resistant dihydrofolate reductase dfrA6.