DNA extraction was performed using Wizard genomic DNA purification kit (Promega, USA) following the manufacturer’s instructions. The Nanopore sequencing library was prepared using the Rapid Barcoding Kit (SQK-RBK004) and sequenced using MinION device (Oxford Nanopore Technologies, Oxford, UK). Raw reads were basecalled using Guppy 1.1. The Illumina DNA library was prepared using Nextera XT DNA Library Preparation Kit (Illumina) and sequenced with the Illumina MiSeq platform. Then, the complete genome sequences of Halioticoli clade type strains were assembled by means of the hybrid assembly approach using both Nanopore and Illumina reads by Unicycler 0.4.7 (Tanaka et al., 2018; Jiang et al., 2022). Finally, the genome sequences were annotated using DDBJ Fast Annotation and Submission Tool (DFAST) (Tanizawa et al., 2018) and deposited in the DDBJ/GenBank/ENA under BioProject PRJDB11924 with accession numbers as Table 2 .

MLSA was performed according to the previous description (Sawabe et al., 2013). Briefly, the entire nucleotide sequences of the eight housekeeping genes (ftsZ, gapA, gyrB, mreB, pyrH, recA, rpoA, and topA) were obtained after genome annotation. The sequences were aligned using MUSCLE (Edgar, 2004). Split decomposition analysis using the concatenated sequence was performed using SplitsTree 4.14.8 with a neighbor net drawing and a Jukes-Cantor correction. Phylogenetic analysis using the same concatenated sequence was constructed using Maximum Likelihood (ML), Neighbor-Joining (NJ), and Minimum-Evolution (ME) methods with 500 bootstraps by MEGA-X v10.1.8 (Kumar et al., 2018).

Genomic comparisons were performed based on chromosomes. The ten genomes from the Halioticoli clade were compared with the genome of V. ishigakensis JCM 19231^T by BLASTn and visualized using BRIG v.0.95 (Alikhan et al., 2011). Synteny between the same genomes was analyzed using the Artemis Comparison Tool (ACT) v.18.1.0 (Carver et al., 2005).

Pan-genome analysis was performed using Halioticoli clade genomes by the Anvi’o program ver. 7 (Eren et al., 2015). Firstly, each genome sequence file was converted to an anvi’o contigs database (anvi-gen-contigs-database) using Prodigal (Hyatt et al., 2010), these contigs databases were decorated with hits from HMM models (anvi-run-hmms). An anvi’o genome storage was generated (anvi-gen-genomes-storage) using prepared contigs databases, and then, the pan-genome was analyzed (anvi-pan-genome) using NCBI’s blastp for amino acid sequence similarity search and the MCL algorithm (Van Dongen and Abreu-Goodger, 2012) for cluster identification in amino acid sequence similarity search results. In addition, Average Nucleotide Identity (ANI) values were calculated using the PyANI with ANIb method (anvi-compute-genome-similarity) (Pritchard et al., 2016). Finally, it was visualized and decorated (anvi-display-pan). Core genes were filtered (anvi-get-sequences-for-gene-clusters) and extracted in fasta files (anvi-get-sequences-for-gene-clusters) for further analysis.

Gene annotation was performed using Clusters of Orthologous Groups 2020 (COG20) (Galperin et al., 2021) for function estimation (anvi-run-ncbi-cogs), and Kyoto Encyclopedia of Genes and Genomes (KEGG) (Aramaki et al., 2020) for metabolism estimation (anvi-estimate-metabolism) (Muto et al., 2013). In addition, the enrichment scores of function/metabolism in different groups were identified using an R script developed by Amy Willis (anvi-compute-functional-enrichment) (Shaiber et al., 2020).

Carbohydrate-Active enZYmes (CAZy) were annotated using the dbCAN2 meta server (HMMdb v9) with HMMER, DIAMOND, and eCAMI tools (Zhang et al., 2018), and domains supported by more than two tools were used in this study. Genomic islands (GEIs) predictions were calculated by IslandViewer4 with IslandPick, IslandPath-DIMOB, and SIGI-HMM methods (Bertelli et al., 2017) using the GenBank files after DFAST annotation, predictions supported by at least one method were used in this study. Results were visualized using ggplot2.

Genomes of all species consisted of two chromosomes and four of them (V. comitans LMG 23416^T, V. halioticoli IAM 14596^T, V. inusitatus LMG 23434^T, and V. rarus LMG 23674^T) had one plasmid ( Table 2 ). The genome sizes of Chromosome 1 (Chr. 1) ranged from 2,528,163 to 3,070,129 bp, and those of Chromosome 2 (Chr. 2) ranged from 989,994 to 1,816,305 bp. These genomes showed 42.9-46.2% GC content, identified 3,101-4,318 CDS, 25-31 rRNA, and 95-107 tRNA. V. ishigakensis JCM 19231^T had the biggest genome with the highest GC content, the largest numbers of CDS, rRNA, and tRNA. In contrast, V. gallicus LMG 21878^T had the smallest genome with the lowest GC content and the smallest numbers of CDS and tRNA. Compared to the genome size and GC content among 189 other Vibrionaceae species, Halioticoli clade species showed a relatively narrow GC content range but variable genome size ( Figure 1 ).

To explore the evolutionary history of Halioticoli clade species, the MLSA network and phylogenetic tree using concatenated eight house-keeping genes with V. cholerae ATCC 14035^T and E. coli K-12 MG1655 as outgroups were constructed ( Figures 2 , S1 ). Both methods indicated five evolutionary directions in the Halioticoli clade: sub-clade 1) V. breoganii, V. comitans, V. inusitatus, and V. superstes, sub-clade 2) V. ezurae, V. neonatus, and V. halioticoli, sub-clade 3) V. rarus, sub-clade 4) V. gallicus, and sub-clade 5) V. ishigakensis. In which, sub-clade 1 to 4 consist of symbiotic species, and sub-clade 5 consists of the only one planktonic species. In addition, genome size and GC content relationship also showed the clustering of sub-clades and the differentiation between symbiotic and planktonic species ( Figure 1 ), suggesting that genome expansions may occur during evolution.

Analyses of homologous gene conservation and gene order across two or more genomes of different species play a vital role in comparative genomics since they can provide further insights into evolutionary processes that contribute to diversity, chromosomal dynamics, and interspecies rearrangement rates (Bhutkar et al., 2006; Lee et al., 2018; de la Haba et al., 2019). The nucleotide identity comparison using BLASTn of both chromosomes (Chr. 1 and Chr. 2) for Halioticoli clade species was performed using V. ishigakensis JCM 19231^T as the reference genome ( Figure 3A ). The analysis revealed a higher nucleotide similarity in Chr. 1 but a lower similarity in Chr. 2, which indicates the genomes of Chr. 1 were highly conserved but those of Chr. 2 was relatively varied. Intra-sub-clade and inter-sub-clade genome rearrangement mappings of both chromosomes demonstrated that 1) similar gene arrangements amongst genomes of intra-subclade species in sub-clades 1 and 2, and 2) less similar of those arrangements among those of inter-subclade species ( Figures 3B , S2 ).

Relatively more GEIs and transposase/integrase were predicted in sub-clades 2 and 3, which were likely to be shared common ancestry ( Figures 2 , S1 , S3 ). The lowest number of GEIs and transposase/integrase were found in V. gallicus (sub-clade 4) and V. comitans (sub-clade 1), respectively. Interestingly, V. ishigakensis showed opposite results that much higher GEIs numbers but much lower transposase/integrase numbers.

Pan-genomics is capable of investigating the relationships between a given group of genomes by means of characterizing the core and accessory genes, providing a unique insight in the phylogeny and taxonomy analysis (Eren et al., 2015; Delmont and Eren, 2018). The ten genomes of Halioticoli clade species were used for pan-genome analysis using Anvi’o v7. In the Halioticoli clade pan-genome ( Figure 4 ), a total of 8,062 gene clusters (GCs) with 36,612 genes were defined, in which 2,130 GCs with 22,123 genes (60%) were recognized in the core-genome (1,973 GCs with 19,730 genes were recognized as the single-copy core-genome), and 2,456 GCs with 10,881 genes (30%) were recognized in the accessory-genome. The remaining genes were recognized as species-specific genes as in Figure 4 , among which, V. ishigakensis (sub-clade 5) possessed the most specific genes (922), 2-5 times more than other species. In addition, gene cluster analysis also showed that sub-clade 5 gained the highest number of GCs (4134), followed by sub-clade 1 (3,732 in average), sub-clade 2 (3,250 in average), sub-clade 3 (3193), and sub-clade 4 (2987). These results showed positive relationships with genome size.

Due to the difficulty to use a large number of single-copy genes (SCGs), a set of 125 better-SCGs was filtered using a custom setting (–min-geometric-homogeneity-index 1, –max-functional-homogeneity-index 0.9). The concatenated amino acid sequence of the 125 better-SCGs was used for constructing a more accurate phylogenetic tree, and the result showed that the topology was congruent with the one constructed by eight house-keeping genes ( Figures 2B , 4 ), the five sub-clades could be identified as well. Moreover, the five sub-clades could also be illustrated by clustering in the ANI matrix, and sub-clade species showed at least 85.5% intra ANI similarity.

The Clusters of Orthologous Genes (COGs) database has been a popular tool for functional and comparative genomics of bacteria and archaea in recent decades with the newest update of COG20 (Galperin et al., 2021). COG20 function estimation for each genome of Halioticoli clade species showed that the same function structure was shared among the species, as well as in the core-genome ( Figure S4A ). However, the result in the species-specific genomes was diverse in some aspects ( Figure 5A ). Despite the poorly characterized and unknown functions, the planktonic species, V. ishigakensis-specific genome had the most diverse and abundant gene sets in each function category. Among them, V. ishigakensis showed more abundant in the metabolism functions, in particular in the “Carbohydrate transport and metabolism” function. In addition, the functions of “Transcription”, “Signal transduction mechanisms”, and “Cell wall membrane envelope biogenesis” were also abundant. It is likely that these functions were an adaptation for survival in planktonic environments. The symbiotic species-specific genomes showed poorer ability in metabolizing, with the V. superstes-specific one being the most powerful. To further investigate the metabolism functions in the Halioticoli clade, metabolism pathway modules were reconstructed using the KEGG database. As in the function structure, the Halioticoli clade species shared the common metabolism structure in complete genomes ( Figure S4B ), but different structures in species-specific genomes ( Figure 5B ). A wide range of genes responsible for diverse metabolism pathways was only detected in the V. ishigakensis-specific genome, in particular carbohydrate metabolism genes were highly detected. On the other hand, other species-specific genomes, demonstrated their advantages, such as a remarkable detection for “cofactor and vitamin metabolism” in the V. rarus-specific genome. In addition, the enrichment analyses of metabolism pathway modules indicated that the V. ishigakensis shared almost all enriched modules with other species with the exception of D-galactonate degradation (M00552), which is exclusively and most enriched in itself ( Figure 6 ).

Furthermore, the COG20 and KEGG annotations were also performed for gene clusters (GCs) of each Halioticoli clade species. Same COG20 functions were shared among the clade but with different abundance. V. ishigakensis, which has the biggest genome size, gained the highest number of GCs with a wide range of COG functions, in particular GCs classified into “Transcription (K)”, “Cell wall/membrane/envelope biogenesis (M)”, and “Carbohydrate transport and metabolism (G)”. In more detail, numbers of GCs were 1.2, 1.6, 1.7, and 1.7 folds in K, 1.1, 1.3, 1.3, and 1.4 folds in M, and 1.4, 1.9, 2.3, and 1.9 folds in G, compared with those numbers of sub-clades 1, 2, 3 and 4, respectively. On the contrary, the lowest number of GCs were observed in V. gallicus, numbers of GCs classified to “Cell cycle control, cell division, chromosome partitioning (D)”, “Posttranslational modification, protein turnover, chaperones (O)”, and “Defense mechanisms (V)” were reduced ( Figure 7A ). KEGG annotation of GCs which were lost or gained showed “Polyamine biosynthesis” (E6, in Figures 7B , S5 ) were gained in sub-clades 1, 3, and 5 but lost in sub-clades 2 and 4, and GCs encoded trans-2,3-dihydro-3-hydroxyanthranilate isomerase [EC:5.3.3.17] in “Biosynthesis of other bacterial compounds (I5)” were gained in sub-clades 1 and 5 but lost in other sub-clades. In more details, glycine/D-amino acid oxidase (deaminating) (dadA) (PDB:3AWI) and acyl-CoA reductase or other NAD-dependent aldehyde dehydrogenase (adhE) (PDB:1A4S), both of which can use putrescine to produce GABA (M00136), were gained in the sub-clades 1 and 3; and genes (rfbB, rfbC, rfbD, and rmlA1) involved in the biosynthesis of the dTDP-L-rhamnose (M00793) were lost in the sub-clade4 ( Figure S6 ).

Carbohydrate-Active enZYmes (CAZy) were predicted for each genome to describe the catalytic modules (enzymes) encoded in these genomes. Generally, each genome of the Halioticoli clade species contained 4 main enzymes classes: Carbohydrate Esterases (CEs), Glycoside Hydrolases (GHs), GlycosylTransferases (GTs), and Polysaccharide Lyases (PLs); and one associated module: Carbohydrate-Binding Modules (CBMs); only some of them contained few Auxiliary Activities (AAs) class enzymes ( Figure 8E ). Among them, Sub-clade 5 (V. ishigakensis) had the most CAZy (120), followed by sub-clade 1 (95-109), as a result of the abundance of GHs in these sub-clades ( Figure 8A ). Meanwhile, kinds of polysaccharide lyases (PL6, PL7, PL15, and P17) involved in alginate degradation were found enriched in the Halioticoli clade species, while V. ishigakensis had fewer numbers (15) than most symbiotic species (19-21) ( Figure 8C ). Compared to the planktonic sub-clade, symbiotic sub-clades showed higher ability on alginate degradation but with different advantages, which is a higher capacity of intracellular degradation (PL15 and PL17) in sub-clade 1, but a higher capacity of extracellular degradation (PL6 and PL7) in sub-clade 2 and 3.These finds indicate that V. ishigakensis obtained a powerful ability for degrading diverse glycosidic bonds but became weaker in degrading polysaccharides during the evolution from the gut environment to the planktonic environment.