Polyphasic identification of 
Vibrio are Gram-negative rod-shaped bacteria often found in aquatic matrices and in close association with marine animals. This bacterial genus has been described as one of the most diverse marine bacterial groups owing to their ubiquity and metabolic and ecologic versatility . Out of more than 150 species described up to date, at least 12 are known for their pathogenic potential – . Apart from the well-known human infections, such as cholera, other forms of gastroenteritis, as well as tissue and wound infections, vibrios are also recurrent pathogens in aquaculture settings. In these environments, Vibrio infections are of great concern as they may trigger mass mortality events and financial losses, thereby posing a risk to food safety , . Therefore, strategies to rapidly and accurately identify these pathogens are of utmost importance, as they will help to investigate the occurrence and distribution of the various species, determining their epidemiology and better define their potential roles in disease and as pollution indicators. Additionally, the use of robust identification methods may aid in designing strategies for treatment and prevention of Vibrio diseases among humans and animals thus having impacts in public health and economy . Identification and differentiation of Vibrio species tend to be complex and sometimes controversial, given the high genetic relatedness between species, and different approaches have been employed for this purpose – . For instance, Matrix Assisted Laser Desorption/Ionization - Time of Flight Mass Spectrometry (MALDI-TOF MS) rises as a promising candidate for fast and relatively accurate identification of Vibrio species . Other identification strategies, such as sequence-based approaches, are considered valuable tools for more reliable identifications. The 16 S rRNA ( rrs ) sequence-based bacterial identification is currently one of the main strategies employed for this purpose. However, previous studies have documented the poor discriminatory power of the rrs gene, notably because rrs gene identities between different Vibrio species were shown to be higher than 97%, thus greatly hindering species differentiation, especially among those more closely-related , . Due to these limitations, the search for alternative phylogenetic markers that could allow for a more reliable identification is warranty. In this context, the pyrH gene, encoding the uridylate kinase, has been suggested for this purpose , , . It is a single-copy gene that belongs to the core genome, displays variation among species and tends to be stable due to a high proportion of synonymous mutations, favoring the conservation of the nucleotide sequence . In the present study, we evaluated three different approaches to identify 40 Vibrio isolates recovered from environmental sources comprising MALDI-TOF MS, and rrs and pyrH gene sequencing. A subset of eight isolates was selected for further characterization by whole genome sequencing. Through this comprehensive analysis, we sought to assess and compare the efficiency of these methods in accurately identifying and distinguishing Vibrio species. Our findings revealed discrepancies in the identification among the methodologies evaluated, suggesting that pyrH gene sequencing provided a superior discrimination between Vibrio species compared to the other methods. 

Bacterial strains and cultivation The 40 bacterial strains included in the present study were isolated from environmental sources in Guanabara Bay, Rio de Janeiro, Brazil as previously described , . Among them, 35 isolates were recovered from water and five, from floating plastic. All strains were grown on Luria Bertani (LB) agar (Kasvi, Madrid, Spain) and incubated at 25 ºC ± 2 ºC for the subsequent analyses. MALDI-TOF MS The strains were initially identified by Matrix-Assisted Laser Desorption/Ionization Time-of-Flight Mass Spectrometry (MALDI-TOF MS). Strains were grown overnight on LB agar as described above. Bacterial biomass was directly disposed onto an MSP 96 target polished steel plate (Bruker Daltonics, Bremen, Germany) and overlaid with 1 µL of formic acid (70%). Upon drying, 1 µL of α-cyano-4-hydroxycinnamic acid was spotted onto the samples. Each isolate was analyzed twice and Escherichia coli DH5α was used as a control in each experiment. Spectra were acquired using the MicroflexLT system (Bruker Daltonics) and the resulting mass spectra were compared to the references in the database using the MALDI Biotyper 7.0 (Bruker, 2019) program. In some instances, retrieving data on particular bacterial strains or their duplicates from the instrument’s database was hindered by either inaccessibility or suboptimal quality. The further analyses were conducted with the main 52 spectra comprising 20 strains with duplicates and 12 without duplicates using BioNumerics (version 7.1, Applied Maths, Kortrijk, Belgium). Raw data were preprocessed and normalized with default parameters, including baseline subtraction using Running Rolling disc and peak detection with a signal-to-noise ratio of 10. Following preprocessing, spectra underwent peak matching with a tolerance of ± 0.002 m/z. Similarity coefficients were calculated with the Pearson coefficient, and the Neighbor joining tree was used to perform hierarchical cluster analysis based on mass spectrometry information, with a permutation resampling of 500, as previously described by Pinto et al., 2017 , with minor adaptations. rrs and pyrH gene sequencing Amplification and sequencing of the rrs (16 S rRNA) and pyrH (uridylate kinase) genes were performed. For amplification of these genes, bacterial DNA was obtained using the thermal lysis method in the presence of Chelex 100 resin (Bio-Rad Laboratories, California, USA) . The primers 27 F 5′-AGAGTTTGATCMTGGCTCAG-3′ and 1492R 5′-CGGTTACCTTGTTACGACTT-3′ (V1-V9 hypervariable regions) ; pyrH-F 5′-GATCGTATGGCTCAAGAAG‐3′ and pyrH-R 5′-TAGGCATTTTGTGGTCACG-3′ were used for the amplification of the rrs and pyrH genes, respectively. The amplicons were purified and sequenced by the Sanger method (BPI, São Paulo, Brazil). All complete rrs and partial pyrH gene sequences were deposited in the GenBank database of NCBI. Accession numbers of all sequences are available in Table . Phylogenetic analyses Sequences were reviewed for quality and edited in the Bioedit software (v.7.2; https://thalljiscience.github.io/ ) . Subsequently, the resulting rrs gene sequences were aligned and taxonomically classified using the BLASTn algorithm from the EzBioCloud and NCBI rrs gene databases . A threshold of 98.7% nucleotide sequence identity was considered for species definition . Taxonomic classification based on the pyrH gene sequences was done exclusively by alignment with the nucleotide sequence database of GenBank of NCBI using the BLASTn tool, by the SINA tool of the SILVA rRNA gene database v1.2.12 , and EzBioCloud 16 S database . For phylogenetic analyses, reference sequences of Vibrio strains were initially selected for each gene. The best fitting model for each gene was determined by running the aligned sequences in ‘find the best DNA/protein model’ tool provided in MEGA v. 11 , followed by phylogenetic reconstruction using the maximum likelihood method based on the Hasegawa-Kishino-Yano nucleotide substitution model (HKY + G + I) for the rrs gene and the General Time Reversible (GTR + G + I) model for the pyrH gene. E. coli U5/41 strain (accession number: NR_024570.1) and E. coli CFT073 strain (accession number: NZ_CP051263.1) were adopted as outgroups. Whole genome sequencing Eight Vibrio strains, three isolated from plastic and five from water samples, were submitted to whole genome sequencing. Total DNA was extracted with the Wizard Genomic DNA Purification Kit A1120, according to the manufacturer’s instructions (Promega, Wisconsin, USA). Library preparation was performed with Nextera DNA Flex Library Prep (Illumina) kit and DNA concentration and quantification were assessed using a QuBit 4 fluorometer with the DNA High Sensitivity test. Genomic DNA was sequenced with the Illumina NextSeq 550. The quality of the read libraries was assessed using the FastQC v0.11.9 tool . Pre-processing was performed with the Trim Galore v0.6.7, wherein remaining adapter sequences, N bases from the read ends, and bases with quality lower than 20 (parameters: --nextera --trim_n --length 50 -q 20) were removed. Readsshorter than 50 bases were also filtered . Genome assembly was conducted with SPAdes v. 3.15.4 with the parameters --isolate --cov-cutoff auto -k 21,33,55 . Gene prediction was performed with Prokka v1.14.15 and the following parameters were used for genome assembly annotation: --compliant --genus [Genus] --usegenus --rfam --evalue 0.0001. The genome completeness was calculated with the CheckM v1.1.11 tool, with the following parameter: lineage_wf -t 8 -x fasta . Genome sequence data were uploaded to the Type (Strain) Genome Server (TYGS), available under https://tygs.dsmz.de , for a whole genome-based taxonomic analysis . Information on nomenclature, synonymy and associated taxonomic literature was provided by TYGS’s sister database, the List of Prokaryotic names with Standing in Nomenclature (LPSN, available at https://lpsn.dsmz.de ) . For the phylogenomic inference, all pairwise comparisons among the set of genomes were conducted using GBDP and accurate intergenomic distances inferred under the algorithm ‘trimming’ and distance formula d5 . 100 distance replicates were calculated each. Digital DDH values and confidence intervals were calculated using the recommended settings of the GGDC 4.0 , . The resulting intergenomic distances were used to infer a balanced minimum evolution tree with branch support via FASTME 2.1.6.1 including SPR postprocessing , . Branch support was inferred from 100 pseudo-bootstrap replicates each. The trees were rooted at the midpoint and visualized with PhyD3 . The genomes have been submitted to the NCBI GenBank database under the BioProjects PRJNA1080184 and PRJNA987140. 

The 40 bacterial strains included in the present study were isolated from environmental sources in Guanabara Bay, Rio de Janeiro, Brazil as previously described , . Among them, 35 isolates were recovered from water and five, from floating plastic. All strains were grown on Luria Bertani (LB) agar (Kasvi, Madrid, Spain) and incubated at 25 ºC ± 2 ºC for the subsequent analyses. 

The strains were initially identified by Matrix-Assisted Laser Desorption/Ionization Time-of-Flight Mass Spectrometry (MALDI-TOF MS). Strains were grown overnight on LB agar as described above. Bacterial biomass was directly disposed onto an MSP 96 target polished steel plate (Bruker Daltonics, Bremen, Germany) and overlaid with 1 µL of formic acid (70%). Upon drying, 1 µL of α-cyano-4-hydroxycinnamic acid was spotted onto the samples. Each isolate was analyzed twice and Escherichia coli DH5α was used as a control in each experiment. Spectra were acquired using the MicroflexLT system (Bruker Daltonics) and the resulting mass spectra were compared to the references in the database using the MALDI Biotyper 7.0 (Bruker, 2019) program. In some instances, retrieving data on particular bacterial strains or their duplicates from the instrument’s database was hindered by either inaccessibility or suboptimal quality. The further analyses were conducted with the main 52 spectra comprising 20 strains with duplicates and 12 without duplicates using BioNumerics (version 7.1, Applied Maths, Kortrijk, Belgium). Raw data were preprocessed and normalized with default parameters, including baseline subtraction using Running Rolling disc and peak detection with a signal-to-noise ratio of 10. Following preprocessing, spectra underwent peak matching with a tolerance of ± 0.002 m/z. Similarity coefficients were calculated with the Pearson coefficient, and the Neighbor joining tree was used to perform hierarchical cluster analysis based on mass spectrometry information, with a permutation resampling of 500, as previously described by Pinto et al., 2017 , with minor adaptations. 

Amplification and sequencing of the rrs (16 S rRNA) and pyrH (uridylate kinase) genes were performed. For amplification of these genes, bacterial DNA was obtained using the thermal lysis method in the presence of Chelex 100 resin (Bio-Rad Laboratories, California, USA) . The primers 27 F 5′-AGAGTTTGATCMTGGCTCAG-3′ and 1492R 5′-CGGTTACCTTGTTACGACTT-3′ (V1-V9 hypervariable regions) ; pyrH-F 5′-GATCGTATGGCTCAAGAAG‐3′ and pyrH-R 5′-TAGGCATTTTGTGGTCACG-3′ were used for the amplification of the rrs and pyrH genes, respectively. The amplicons were purified and sequenced by the Sanger method (BPI, São Paulo, Brazil). All complete rrs and partial pyrH gene sequences were deposited in the GenBank database of NCBI. Accession numbers of all sequences are available in Table . 

Sequences were reviewed for quality and edited in the Bioedit software (v.7.2; https://thalljiscience.github.io/ ) . Subsequently, the resulting rrs gene sequences were aligned and taxonomically classified using the BLASTn algorithm from the EzBioCloud and NCBI rrs gene databases . A threshold of 98.7% nucleotide sequence identity was considered for species definition . Taxonomic classification based on the pyrH gene sequences was done exclusively by alignment with the nucleotide sequence database of GenBank of NCBI using the BLASTn tool, by the SINA tool of the SILVA rRNA gene database v1.2.12 , and EzBioCloud 16 S database . For phylogenetic analyses, reference sequences of Vibrio strains were initially selected for each gene. The best fitting model for each gene was determined by running the aligned sequences in ‘find the best DNA/protein model’ tool provided in MEGA v. 11 , followed by phylogenetic reconstruction using the maximum likelihood method based on the Hasegawa-Kishino-Yano nucleotide substitution model (HKY + G + I) for the rrs gene and the General Time Reversible (GTR + G + I) model for the pyrH gene. E. coli U5/41 strain (accession number: NR_024570.1) and E. coli CFT073 strain (accession number: NZ_CP051263.1) were adopted as outgroups. 

Eight Vibrio strains, three isolated from plastic and five from water samples, were submitted to whole genome sequencing. Total DNA was extracted with the Wizard Genomic DNA Purification Kit A1120, according to the manufacturer’s instructions (Promega, Wisconsin, USA). Library preparation was performed with Nextera DNA Flex Library Prep (Illumina) kit and DNA concentration and quantification were assessed using a QuBit 4 fluorometer with the DNA High Sensitivity test. Genomic DNA was sequenced with the Illumina NextSeq 550. The quality of the read libraries was assessed using the FastQC v0.11.9 tool . Pre-processing was performed with the Trim Galore v0.6.7, wherein remaining adapter sequences, N bases from the read ends, and bases with quality lower than 20 (parameters: --nextera --trim_n --length 50 -q 20) were removed. Readsshorter than 50 bases were also filtered . Genome assembly was conducted with SPAdes v. 3.15.4 with the parameters --isolate --cov-cutoff auto -k 21,33,55 . Gene prediction was performed with Prokka v1.14.15 and the following parameters were used for genome assembly annotation: --compliant --genus [Genus] --usegenus --rfam --evalue 0.0001. The genome completeness was calculated with the CheckM v1.1.11 tool, with the following parameter: lineage_wf -t 8 -x fasta . Genome sequence data were uploaded to the Type (Strain) Genome Server (TYGS), available under https://tygs.dsmz.de , for a whole genome-based taxonomic analysis . Information on nomenclature, synonymy and associated taxonomic literature was provided by TYGS’s sister database, the List of Prokaryotic names with Standing in Nomenclature (LPSN, available at https://lpsn.dsmz.de ) . For the phylogenomic inference, all pairwise comparisons among the set of genomes were conducted using GBDP and accurate intergenomic distances inferred under the algorithm ‘trimming’ and distance formula d5 . 100 distance replicates were calculated each. Digital DDH values and confidence intervals were calculated using the recommended settings of the GGDC 4.0 , . The resulting intergenomic distances were used to infer a balanced minimum evolution tree with branch support via FASTME 2.1.6.1 including SPR postprocessing , . Branch support was inferred from 100 pseudo-bootstrap replicates each. The trees were rooted at the midpoint and visualized with PhyD3 . The genomes have been submitted to the NCBI GenBank database under the BioProjects PRJNA1080184 and PRJNA987140. 

Identification by MALDI-TOF MS All 40 strains were identified as belonging to the Vibrio genus using MALDI-TOF MS, revealing seven different species (Fig. ). The lowest identification score obtained was 1.708 for strain Vibrio parahaemolyticus B34F40, suggesting probable genus-level identification. Our results suggest that there was a tendency to affiliate environmental strains with pathogenic species, such as Vibrio alginolyticus and Vibrio harveyi . Therefore, discrimination between environmental strains was typically challenging, leading to low identification scores. Additionally, 60% of the strains ( n = 24) were identified with confidence scores indicative of genus-level identification only. This was likely due to the limitation of the database (Biotyper v. 3.1), which, at the time of the study, harbored 122 reference spectra for the Vibrionaceae family and 111 for the Vibrio genus, among which nine corresponded to the V. parahaemolyticus species, seven to V. harveyi , and seven to V. anguillarum , for instance. Despite being a relatively fast and reliable technique, identification by MALDI-TOF MS is highly dependent on the content of the database used for spectral analysis . Here, we highlight that 50% of the tested strains were identified as V. harveyi . The Harveyi clade currently comprises 12 species: V. harveyi , V. alginolyticus , V. parahaemolyticus , Vibrio diabolicus , Vibrio rotiferianus , Vibrio owensii , Vibrio campbellii , Vibrio jasicid a, Vibrio azureus , Vibrio sagamiensis , Vibrio natriegens , and Vibrio mytili , , . Difficulties in the differentiation of closely related species belonging to the same clade may explain the low species diversity obtained by MALDI-TOF MS. The use of specific databases enhances bacterial identification. This has been previously reported by Håkonsholm and colleagues (2020) , who were able to adequately identify 99% of the strains tested by using an in-house database containing spectra from Vibrio spp. known to be associated with marine fish. Application of such approach is especially useful in the context of animal mortality events in aquaculture settings, where Vibrio species are among the main bacterial pathogens associated with fish disease. As an aggravating factor, Vibrio infections often lead to significant production and financial losses . In these cases, species such as V. harveyi or V. alginolyticus are closely associated with infections in shrimps, often causing great losses in Penaeus spp. and Litopenaeus vannamei aquaculture . Additionally, other works have also demonstrated the higher performance of in-house databases in identifying potentially pathogenic Vibrio spp – . Thus, fast and accurate identification is important to develop and implement prophylactic and risk-management strategies, as well as to monitor the occurrence and spread of infections , . Further, cluster analysis based on MALDI-TOF MS spectra revealed that V. harveyi strains showed a tendency to cluster together (Fig. ). Herein, two main clusters, defined as V. harveyi groups, were observed, together encompassing almost all of the V. harveyi isolates (16/18) in the present study. On the other hand, the less represented species ( Vibrio pomeroyi and V. anguillarum ) showed a more dispersed distribution across the dendrogram, indicating that clustering might become more refined as more strains are included in the analysis. Overall, there was a tendency that strains belonging to the same species grouped together, with few exceptions. Namely, strains B34F40, identified as V. parahaemolyticus , and B1S46 and B34S11, both identified as V. alginolyticus , clustered together with the V. harveyi group. This is not entirely unexpected owing to the fact that V. parahaemolyticus and V. alginolyticus belong to the Harveyi clade according to a recent update of the molecular phylogenetics of the Vibrionaceae family . Also, one strain of V. alginolyticus (B1S27) was found in close proximity to V. pomeroyi (B34S44) rather than with other strains identified as V. alginolyticus . In some cases, these discrepancies in the distribution of Vibrio species across the dendrogram were also observed among duplicates of the same strain (Fig. ). For instance, the strain V. parahaemolyticus B7F71 clustered with other strains of V. harveyi or V. alginolyticus , depending on the duplicate analyzed. These findings suggest a considerable variation in the MALDI-TOF MS spectra obtained from technical replicates of identical isolates. The misidentification of species through MALDI-TOF MS analysis cannot be overlooked, particularly given that variance in identification outcomes depending on the methodology employed. Identification by rrs and pyrH gene sequencing and phylogenetic analyses Considering that the identities between rrs gene sequences among Vibrio species are very high and may lead to low discriminatory power among strains of closely related species, we investigated an additional gene, pyrH gene, as an alternative target for identification of environmental isolates. With the exception of two strains (B7F82 and B34F91) that failed to amplify the pyrH gene, B7F82 and B34F91, all were classified as belonging to the Vibrio genus ( n = 38; 95%). Strain B7F82 was indeed not identified as Vibrio spp. by rrs gene sequencing, contrary to B34F91. Additionally, a strain (B34F73) identified as belonging to the Vibrio genus by both MALDI-TOF MS and pyrH gene sequencing was identified as Providencia sp. by rrs gene sequencing. The rrs amplicon sequences were analyzed using three distinct tools, BLASTn NCBI, EzBioCloud, and SILVA SINA, as the use of different databases minimizes the risks of misidentification, especially considering non-clinical isolates . However, we found that identification could diverge when different databases were used. Using SILVA SINA alignment service of the ARB-Silva database and the BLASTn NCBI, a considerable number of strains were identified as “ Vibrio sp.” ( n = 30; 75% and n = 13; 32.5%, respectively), indicating that species-level identification was not possible with the identity threshold of 98.7%. On the other hand, through the analysis of the pyrH amplicon sequences, only two strains (5%) were identified as “ Vibrio sp.“, while all the others were identified at the species level. This is in agreement with previous studies that have demonstrated the limited capacity of the rrs gene in distinguishing Vibrio species , , and the higher discriminatory power of the pyrH gene. For most strains, identification obtained by pyrH gene sequencing diverged from the results obtained by both rrs gene sequencing and MALDI-TOF MS analysis. Notably, only three strains, OL1V-2, OL2S-1, and OL3V-2, all isolated from plastic and identified as Vibrio fluvialis demonstrated consistent identification across nearly all the methodologies used. The only exception arose with the SILVA SINA database, where these strains were labeled as Vibrio sp. On the other hand, certain species were detected only by specific methodologies, such as Vibrio pomeroyi , detected only using MALDI-TOF MS, and Vibrio sagamiensis or Vibrio owensii , detected only by rrs or pyrH gene sequencing. The identification obtained for each strain through each of the mentioned methodologies is described in Table . It is also worth noting that, through rrs gene sequencing, six strains analyzed by BLASTn NCBI and 14 analyzed by EzBioCloud presented identities lower than 98.7%, which is considered the cutoff value for species delimitation . To better understand how rrs and pyrH genes can assist in distinguishing between Vibrio species, phylogenetic analyses were conducted incorporating sequences from reference strains. The phylogenetic tree constructed based on rrs gene sequences is represented in Fig. . Except for the reference sequences of V. neocaledonicus and V. fluvialis , the reference sequences as a whole did not form clusters with the sequences derived from the isolates included in this study. Some test sequences, such as B1S51 and B1S46, both identified as V. alginolyticus (NCBI) or V. sagamiensis (EzBioCloud), were placed close to each other but distant from their corresponding reference sequences. In contrast, in the phylogenetic tree based on pyrH gene sequences (Fig. ), strains B1S51 and B1S46, as well as others identified as V. alginolyticus , formed a cluster with the reference sequence ( V. alginolyticus NBRC 15630/ATCC 17749). Overall, a better differentiation between species based on pyrH gene sequences was obtained. For example, strains with identical identifications formed clusters in the pyrH gene tree and were distant in the rrs gene tree (e.g., B7F71 and B34S58). This is in agreement with the results obtained by Amalina and colleagues (2019) , who have described the pyrH gene as a robust phylomarker due to its high discriminatory power. Additionally, the species V. harveyi , often identified by the MALDI-TOF MS, did not cluster with other strains also identified as V. harveyi by MALDI-TOF MS. Indeed, only in the phylogenetic tree corresponding to the pyrH gene, V. harveyi ATCC 14,126 was placed close to the sole strain identified as V. harveyi by all three approaches evaluated, thus suggesting that MALDI-TOF MS incorrectly identifies many strains as V. harveyi . Whole genome sequencing Among the eight strains submitted to whole genome sequencing (B1S33, B1S50, B1S-51, B34S44, B34S48, OL1V1, OL2S2, and OL3V3) five were identified as V. alginolyticus , two as V. fluvialis , and one as V. parahaemolyticus . Overall, results derived from the genomic sequences were in agreement with the data obtained through sequencing of the pyrH amplicons. The genomes of the two V. fluvialis strains grouped together with V. fluvialis NBRC 103,150, a synonym for V. fluvialis ATCC 33,809, reinforcing the results found using the pyrH gene sequencing. As to V. parahaemolyticus , both the genome and pyrH gene sequences grouped together with V. parahaemolyticus ATCC 17,802/NBRC 12,711. The same was observed for V. alginolyticus strains, which were grouped together with the reference strain V. alginolyticus NRBC 15,630 (Fig. ). 

All 40 strains were identified as belonging to the Vibrio genus using MALDI-TOF MS, revealing seven different species (Fig. ). The lowest identification score obtained was 1.708 for strain Vibrio parahaemolyticus B34F40, suggesting probable genus-level identification. Our results suggest that there was a tendency to affiliate environmental strains with pathogenic species, such as Vibrio alginolyticus and Vibrio harveyi . Therefore, discrimination between environmental strains was typically challenging, leading to low identification scores. Additionally, 60% of the strains ( n = 24) were identified with confidence scores indicative of genus-level identification only. This was likely due to the limitation of the database (Biotyper v. 3.1), which, at the time of the study, harbored 122 reference spectra for the Vibrionaceae family and 111 for the Vibrio genus, among which nine corresponded to the V. parahaemolyticus species, seven to V. harveyi , and seven to V. anguillarum , for instance. Despite being a relatively fast and reliable technique, identification by MALDI-TOF MS is highly dependent on the content of the database used for spectral analysis . Here, we highlight that 50% of the tested strains were identified as V. harveyi . The Harveyi clade currently comprises 12 species: V. harveyi , V. alginolyticus , V. parahaemolyticus , Vibrio diabolicus , Vibrio rotiferianus , Vibrio owensii , Vibrio campbellii , Vibrio jasicid a, Vibrio azureus , Vibrio sagamiensis , Vibrio natriegens , and Vibrio mytili , , . Difficulties in the differentiation of closely related species belonging to the same clade may explain the low species diversity obtained by MALDI-TOF MS. The use of specific databases enhances bacterial identification. This has been previously reported by Håkonsholm and colleagues (2020) , who were able to adequately identify 99% of the strains tested by using an in-house database containing spectra from Vibrio spp. known to be associated with marine fish. Application of such approach is especially useful in the context of animal mortality events in aquaculture settings, where Vibrio species are among the main bacterial pathogens associated with fish disease. As an aggravating factor, Vibrio infections often lead to significant production and financial losses . In these cases, species such as V. harveyi or V. alginolyticus are closely associated with infections in shrimps, often causing great losses in Penaeus spp. and Litopenaeus vannamei aquaculture . Additionally, other works have also demonstrated the higher performance of in-house databases in identifying potentially pathogenic Vibrio spp – . Thus, fast and accurate identification is important to develop and implement prophylactic and risk-management strategies, as well as to monitor the occurrence and spread of infections , . Further, cluster analysis based on MALDI-TOF MS spectra revealed that V. harveyi strains showed a tendency to cluster together (Fig. ). Herein, two main clusters, defined as V. harveyi groups, were observed, together encompassing almost all of the V. harveyi isolates (16/18) in the present study. On the other hand, the less represented species ( Vibrio pomeroyi and V. anguillarum ) showed a more dispersed distribution across the dendrogram, indicating that clustering might become more refined as more strains are included in the analysis. Overall, there was a tendency that strains belonging to the same species grouped together, with few exceptions. Namely, strains B34F40, identified as V. parahaemolyticus , and B1S46 and B34S11, both identified as V. alginolyticus , clustered together with the V. harveyi group. This is not entirely unexpected owing to the fact that V. parahaemolyticus and V. alginolyticus belong to the Harveyi clade according to a recent update of the molecular phylogenetics of the Vibrionaceae family . Also, one strain of V. alginolyticus (B1S27) was found in close proximity to V. pomeroyi (B34S44) rather than with other strains identified as V. alginolyticus . In some cases, these discrepancies in the distribution of Vibrio species across the dendrogram were also observed among duplicates of the same strain (Fig. ). For instance, the strain V. parahaemolyticus B7F71 clustered with other strains of V. harveyi or V. alginolyticus , depending on the duplicate analyzed. These findings suggest a considerable variation in the MALDI-TOF MS spectra obtained from technical replicates of identical isolates. The misidentification of species through MALDI-TOF MS analysis cannot be overlooked, particularly given that variance in identification outcomes depending on the methodology employed. 

Considering that the identities between rrs gene sequences among Vibrio species are very high and may lead to low discriminatory power among strains of closely related species, we investigated an additional gene, pyrH gene, as an alternative target for identification of environmental isolates. With the exception of two strains (B7F82 and B34F91) that failed to amplify the pyrH gene, B7F82 and B34F91, all were classified as belonging to the Vibrio genus ( n = 38; 95%). Strain B7F82 was indeed not identified as Vibrio spp. by rrs gene sequencing, contrary to B34F91. Additionally, a strain (B34F73) identified as belonging to the Vibrio genus by both MALDI-TOF MS and pyrH gene sequencing was identified as Providencia sp. by rrs gene sequencing. The rrs amplicon sequences were analyzed using three distinct tools, BLASTn NCBI, EzBioCloud, and SILVA SINA, as the use of different databases minimizes the risks of misidentification, especially considering non-clinical isolates . However, we found that identification could diverge when different databases were used. Using SILVA SINA alignment service of the ARB-Silva database and the BLASTn NCBI, a considerable number of strains were identified as “ Vibrio sp.” ( n = 30; 75% and n = 13; 32.5%, respectively), indicating that species-level identification was not possible with the identity threshold of 98.7%. On the other hand, through the analysis of the pyrH amplicon sequences, only two strains (5%) were identified as “ Vibrio sp.“, while all the others were identified at the species level. This is in agreement with previous studies that have demonstrated the limited capacity of the rrs gene in distinguishing Vibrio species , , and the higher discriminatory power of the pyrH gene. For most strains, identification obtained by pyrH gene sequencing diverged from the results obtained by both rrs gene sequencing and MALDI-TOF MS analysis. Notably, only three strains, OL1V-2, OL2S-1, and OL3V-2, all isolated from plastic and identified as Vibrio fluvialis demonstrated consistent identification across nearly all the methodologies used. The only exception arose with the SILVA SINA database, where these strains were labeled as Vibrio sp. On the other hand, certain species were detected only by specific methodologies, such as Vibrio pomeroyi , detected only using MALDI-TOF MS, and Vibrio sagamiensis or Vibrio owensii , detected only by rrs or pyrH gene sequencing. The identification obtained for each strain through each of the mentioned methodologies is described in Table . It is also worth noting that, through rrs gene sequencing, six strains analyzed by BLASTn NCBI and 14 analyzed by EzBioCloud presented identities lower than 98.7%, which is considered the cutoff value for species delimitation . To better understand how rrs and pyrH genes can assist in distinguishing between Vibrio species, phylogenetic analyses were conducted incorporating sequences from reference strains. The phylogenetic tree constructed based on rrs gene sequences is represented in Fig. . Except for the reference sequences of V. neocaledonicus and V. fluvialis , the reference sequences as a whole did not form clusters with the sequences derived from the isolates included in this study. Some test sequences, such as B1S51 and B1S46, both identified as V. alginolyticus (NCBI) or V. sagamiensis (EzBioCloud), were placed close to each other but distant from their corresponding reference sequences. In contrast, in the phylogenetic tree based on pyrH gene sequences (Fig. ), strains B1S51 and B1S46, as well as others identified as V. alginolyticus , formed a cluster with the reference sequence ( V. alginolyticus NBRC 15630/ATCC 17749). Overall, a better differentiation between species based on pyrH gene sequences was obtained. For example, strains with identical identifications formed clusters in the pyrH gene tree and were distant in the rrs gene tree (e.g., B7F71 and B34S58). This is in agreement with the results obtained by Amalina and colleagues (2019) , who have described the pyrH gene as a robust phylomarker due to its high discriminatory power. Additionally, the species V. harveyi , often identified by the MALDI-TOF MS, did not cluster with other strains also identified as V. harveyi by MALDI-TOF MS. Indeed, only in the phylogenetic tree corresponding to the pyrH gene, V. harveyi ATCC 14,126 was placed close to the sole strain identified as V. harveyi by all three approaches evaluated, thus suggesting that MALDI-TOF MS incorrectly identifies many strains as V. harveyi . 

Among the eight strains submitted to whole genome sequencing (B1S33, B1S50, B1S-51, B34S44, B34S48, OL1V1, OL2S2, and OL3V3) five were identified as V. alginolyticus , two as V. fluvialis , and one as V. parahaemolyticus . Overall, results derived from the genomic sequences were in agreement with the data obtained through sequencing of the pyrH amplicons. The genomes of the two V. fluvialis strains grouped together with V. fluvialis NBRC 103,150, a synonym for V. fluvialis ATCC 33,809, reinforcing the results found using the pyrH gene sequencing. As to V. parahaemolyticus , both the genome and pyrH gene sequences grouped together with V. parahaemolyticus ATCC 17,802/NBRC 12,711. The same was observed for V. alginolyticus strains, which were grouped together with the reference strain V. alginolyticus NRBC 15,630 (Fig. ). 

The precise identification of strains is a crucial step to get a comprehensive understanding of bacterial dynamics, particularly in the context of disease outbreaks. However, several aspects are usually taken into account for selection of methods, such as costs, time required for results, yielded information, ease of implementation, and the need of specialized expertise, e.g. bioinformatics . Despite the limited number of samples analyzed in our study, we reinforce that they encompass relevant pathogenic species and few previous studies have delved into the field of their identification. Our results indicate that pyrH gene sequencing is a promising alternative procedure for reliable identification of Vibrio strains, as corroborated by whole genome sequencing, which is currently the gold standard for pathogen source tracking. Nonetheless, despite its high accuracy, the widespread adoption of whole genome sequencing is hindered by its substantial costs and the need for expertise in bioinformatics, thus limiting its accessibility across various contexts. Although MALDI-TOF MS offers fast identification, it is highly dependent on the spectra database, and, in our study, there was a tendency to affiliate environmental strains with potentially pathogenic species. Therefore, we suggest that identification of Vibrio spp. would greatly benefit from standardization efforts and curation of databases, not only for MALDI-TOF, but also regarding pyrH and rrs gene databases. 

Below is the link to the electronic supplementary material. Supplementary Material 1 

