Comparative genomic analysis of Flavobacterium species causing columnaris disease of freshwater fish in Thailand: insights into virulence and resistance mechanisms

Abstract

Background

Columnaris disease, a prevalent disease among farmed and wild freshwater fish, is caused by the Flavobacterium columnare group, which includes four distinct species: F. columnare, F. oreochromis, F. covae, and F. davisii. Among these, F. oreochromis, F. covae, and F. davisii are particularly prevalent in farmed freshwater fish in Thailand. In this study, a comparative genomic analysis of 22 isolates was conducted to elucidate virulence factors, antibiotic resistance genes (ARGs), genomic islands (GIs), phages, insertion elements (ISs), and clustered regularly interspaced short palindromic repeats (CRISPRs).

Results

A total of 212 putative virulence genes were predicted across three species with F. oreochromis exhibiting the highest number of unique virulence genes, followed by F. davisii, and F. covae. Moreover, 195 genes were predicted as ARGs, with F. oreochromis and F. covae showing an abundance of unique genes associated with resistance to quinolone, fluoroquinolone, and tetracycline antibiotics. Antimicrobial susceptibility testing, assessed with epidemiological cut-off values (ECVs), revealed decreased susceptibility to quinolones, fluoroquinolones and tetracycline in several isolates of F. oreochromis and F. covae. F. oreochromis and F. covae exhibited notable decreased susceptibility to quinolones, with mutations observed in the quinolone resistance-determining regions (QRDRs) of gyrA, including Ser83Phe, Ser83Val, Ser83Ala, and Asp87Tyr, the latter representing a novel mutation among isolates from Thailand. As a result, these findings suggest that gyrA is major target for quinolone in F. oreochromis, F. covae, and F. davisii, while gyrB, parC, parE might be less important to the decreased phenotypic susceptibility to this class of antimicrobials. Moreover, a tetracycline resistance gene (tetA_2) was found in only one F. covae isolate, which exhibited decreased phenotypic susceptibility to this drug, marking the first report of decreased susceptibility in this species.

Conclusions

This study provides insights into the genetic and pathogenic diversity of Flavobacterium species, aiding in the development of strategies to manage columnaris disease in farmed freshwater fish in Thailand.

Background

The aquaculture industry is experiencing the most rapid growth among all sectors of global animal food production. In 2019, Thailand ranked among the top 10 aquaculture producers, with an output of one million tons [1, 2]. The intensification of production to meet the demand for fish products has resulted in epizootics of various bacterial and viral pathogens [1, 3, 4]. Columnaris disease is a serious bacterial infection caused by Flavobacterium columnare group, which are Gram-negative, slender, motile, and rod-shaped bacteria that effect a wide range of freshwater fish species, leading to symptoms such as skin lesions, gill inflammation, and high mortality rates. This disease poses a major threat to aquaculture industries worldwide, resulting in significant economic losses and necessitating effective management and prevention strategies [4,5,6,7,8,9].

Flavobacterium columnare group has been extensively studied regarding its colony morphology, biochemical and physiological characteristic. Numerous studies have revealed phenotypic homogeneity, a large degree of genetic diversity, and differing levels of virulence exhibited by Flavobacterium columnare group across various hosts, using multiple molecular techniques and genomic analysis [5, 6, 10,11,12,13]. Based on phenotypic and genomic analyses, Flavobacterium columnare group causing the columnaris disease were recently categorized into four distinct species, F. columnare, F. covae, F. davisii, and F. oreochromis [13]. These four species have been collectively referred to as ‘columnaris-causing bacteria (CCB)’ [9]. They are commonly isolated from a wide range of freshwater fish species, including red tilapia (Oreochromis sp.), Nile tilapia (O. niloticus), Asian sea bass (Lates calcarifer), channel catfish (Ictalurus punctatus), rainbow trout (Oncorhynchus mykiss), and striped catfish (Pangasianodon hypophthalmus) [6, 13,14,15,16,17,18,19,20,21].

In addition to providing a comprehensive view of the virulome and resistome, comparative genomic analysis is a highly effective tool providing valuable insights into the genetics of virulence, antibiotic resistance, evolution, and potential strategies for management and control. The objective of this study was to identify different genomic features, such as potential virulence factors, antibiotic resistance genes, genomic islands, phages, insertion elements, and clustered regularly interspaced short palindromic repeats (CRISPRs), that contribute to the pathogenicity of F. oreochromis, F. covae and F. davisii. Research findings from this study may be used to develop effective control strategies for the management of columnaris disease in intensive freshwater fish farming in Thailand.

Methods

Bacterial isolates

A total of twelve Thai CCB isolates were obtained from the collection maintained at the center of excellence in fish infectious diseases (CE FID), Faculty of Veterinary Science of Chulalongkorn University (Bangkok, Thailand). Samples were isolated from three species of economically important diseased freshwater fish, including red tilapia, Asian seabass, and snakeskin gourami (Trichopodus pectoralis), from farms located across various regions of Thailand between 2019 and 2023, as detailed in Table 1. Fish showing clinical symptoms of Flavobacterium columnare group infection, such as gill necrosis, fin erosion, or skin ulceration, were collected as described in previous publications [5, 8] (Fig. 1). Bacteria isolates were primarily collected from necrotic gills and cultured on Anacker and Ordal (AO) agar and incubated at 28 °C for 48 h [5, 8, 22]. The CCB isolates were phenotypically identified by observing colony morphology and performing conventional biochemical tests before being preserved in an AO-based medium supplemented with 10% glycerol and 20% fetal bovine serum at -80 °C at the CE FID.

Table 1 Isolation of F. oreochromis, F. covae, and F. davisii from diseased farmed freshwater fish in Thailand

Fig. 1

Lesions caused by columnaris disease in red tilapia (Oreochromis sp.), showing necrosis in the gill filaments (red arrow)

Molecular identification and phylogenetic analysis

The species identity of the isolates was confirmed as CCB and assigned to specific species by amplifying and sequencing their 16S rRNA genes. A phylogenetic analysis was then performed, incorporating genes from CCB type and reference strains of F. oreochromis (Costa Rica 04–02-TN^T), F. covae (AL-02-36 ^T, 94–081, C#2, LV-359–01), F. davisii (90-106 ^T), and F. columnare (ATCC 23463 ^T, ATCC 49512). Genomic DNA extraction was carried out in accordance with the manufacturer’s instructions using the NucleoSpin® Tissue (Macherey–Nagel, Düren, Germany), and the isolated DNA samples were stored at -20 °C until needed. Amplification of the 16S rRNA gene was performed using a pair of universal primers, 20F (5′- AGA GTT TGA TCM TGG CTC AG -3′) and 1492R (5’- GGT TAC CTT GTT ACG ACT T-3’), with expected product size of 1,475 bp, and PCR conditions were as described by LaFrentz et al. [23]. The resulting DNA fragments were visualized through gel electrophoresis using 1% w/v agarose, stained with RedSafe™ (Intron, Gyeonggi-do, Korea), and observed under UV light. The PCR products were purified from agarose gel using NucleoSpin® Gel and PCR clean-up kit (Macherey–Nagel, Düren, Germany). Subsequently, the purified products were submitted for sequencing using the 20F and 1492R primers at First Base in Kuala lumpur, Malaysia.

Phylogenetic analysis based on the 16S rRNA gene was performed using the gene sequences from the 12 Thai CCB isolates from this study and an additional 13 Thai CCB isolates from previous studies [6, 8, 24]. The 16S rRNA genes from the type and reference strains of F. oreochromis, F. covae, F. davisii, and F. columnare were included, along with F. psychrophilum JIP02/86 as an outgroup (Table S1). The full-length 16S rRNA sequences were aligned using ClustalW tool in the Molecular Evolutionary Genetics Analysis (MEGA X) software [25]. The nucleotide substitution model with lowest Bayesian Information Criterion (BIC) scores was selected for 16S rRNA-based tree, which was then constructed using a maximum likelihood method based on the Kimura 2-parameter model (K2 + G + I) with 1000 bootstraps [26].

Antimicrobial susceptibility tests by minimum inhibitory concentration (MIC) assay

The broth microdilution method was used to determine the minimum inhibitory concentrations (MICs), following the protocol outlined in the CLSI VET04 guidelines [27, 28]. Briefly, three drugs were diluted two-fold into ten different concentrations of oxolinic acid (OA; Quinolone; 0.008–4 µg/mL), enrofloxacin (ENR; Fluoroquinolone; 0.002–1 mg/µL), and oxytetracycline (OT; Tetracycline; 0.03–16 µg/mL) using cation-adjusted Muller-Hinton broth (CAMHB) containing 3 g/L of Mueller–Hinton broth supplemented with cations (4 mg/L Ca2⁺ and 2 mg/L Mg2⁺). Fifty µL of drugs were dispensed into 96-well plates and 50 µL of each bacterial isolate (5 × 10⁵ CFU/mL) was added to each well. Afterward, the plates were incubated at 28 °C for 48 h. Positive and negative controls were established using only-CAMHB and bacterial suspension wells, respectively. Escherichia coli ATCC 25922 was used as a quality control strain for quality assurance [29]. The minimum inhibitory concentration (MIC) was evaluated using epidemiological cutoff values (ECVs) established in the previous investigation [30]. These ECVs were set at ≤ 0.25 µg/mL for OA, ≤ 0.03 µg/mL for ENR, and ≤ 0.25 µg/mL for OT [30]. Isolates with MIC values at or below the ECVs were classified as having no decreased susceptibility (wild-type, WT), while those with MIC values above the ECVs were categorized as having decreased susceptibility (non-wild-type, NWT). Nine isolates from three different species, including OA wild-type (MIC ≤ 0.25 µg/mL) and OA non-wild-type (MIC > 0.25 µg/mL) isolates with varying susceptibility to OA, were selected for next-generation sequencing.

Whole genome sequencing, assembly, and annotation

Nine isolates, selected based on their OA MIC values from the twelve tested, were combined with 13 isolates from previous studies, comprising F. oreochromis (1214, CC1808, CC1805, CC1803, CC1802, 1201, NK01), F. covae (SP1802, CF1, SP1809, 1362, 1372), and F. davisii (1215), bringing the total to 22 Thai isolates included in the comparative analysis [6, 8, 24]. CCB type strains of F. oreochromis (Costa Rica 04–02-TN^T), F. covae (AL-02-36 ^T), and F. davisii (90-106 ^T) were also included. In this study, the quality of nine genomic DNA samples was assessed using a QubitTM 4 Fluorometer (InvitrogenTM, Singapore) and submitted to a next-generation sequencing service (Vishuo Biomedical, Beijing, China). The NEBNext® UltraTM DNA Library Prep Kit for Illumina was used to create the sequencing library (New England Biolabs, Massachusetts, USA), and the samples were sequenced using an Illumina HiSeq platform using paired-end mode with a read length of 150 bp. Before analysis, all reads were filtered to remove low-quality bases (Q score < 30), short length (< 36 nt), and adapter sequences using Trimmomatic (ver. 0.32) [31]. The improvement of read quality was assessed using FastQC (ver. 0.11.8) [32]. Trimmed reads were then de novo assembled using SPAdes (ver. 3.15.5) [33], and assembly quality was evaluated using the Quast (ver. 5.2.0) [34]. Subsequently, Prokka (ver. 1.14.6) was utilized to annotate all genomic contigs and re-annotate previous genomes [35]. In total, nine complete genome sequences in this study (KCRT2301, CNRT2201, RCBRT1901, UBRT2201, SKNRT2101, KCRT2304, KCRT1906, PCBSB2203, KCRT2007) were deposited with GenBank through WGS submission under accession numbers GCA_045055415.1, GCA_045055375.1, GCA_045055575.1, GCA_045055515.1, GCA_045055495.1, GCA_045055455.1, GCA_045055395.1, GCA_045055595.1, and GCA_045055305.1, respectively.

Pan-genome analysis

An analysis was conducted to examine the distribution of genes constituting the core genome and accessory genome and to identify strain-specific genes from each genome. The pan-genome analysis of 22 isolates was performed using the BPGA pan-genome analysis pipeline [36]. The amino acid sequences obtained from Prokka annotation were utilized as input files. Protein homologs were grouped using USEARCH with an identity cut-off at 50%. The pan-genome and core genome were determined through iterative calculations, using exponential growth and decay models with each subsequent addition of a new genome. For the construction of core genome-based phylogenetic trees, the BPGA pipeline automated multiple sequence alignments using MUSCLE [36]. These alignments were then concatenated, and a maximum likelihood phylogenetic tree based on core genes was generated using FastTree (ver. 2.1.10) with 1000 bootstrap replicates [37]. Based on BLAST distance phylogenetics (GBDP) distances derived from genome sequences, Type Genome Server (TYGS) facilitated the construction of neighbor-joining phylogenetic trees using FastME (ver. 2.1.4) [38, 39].

Comparative genome analysis

The SEED subsystem categorization of 22 genomes was conducted using the RAST server [40, 41]. Virulence factors were predicted using a local Protein–Protein BLAST (ver. 2.2.28) against the Virulence Factor Database (VFDB) (http://www.mgc.ac.cn/VFs/download.htm, accessed in December 2023) [42, 43]. A cut-off E-value of 10^–4 was considered significant. Additionally, the presence of genes associated with secretion systems in F. oreochromis, F. covae, and F. davisii was investigated using MacSyFinder (ver. 2.0) and TXSScan [44].

The resistome, referring to ARGs, was analyzed through the Resistance Gene Identifier (RGI; ver. 6.0.3) web portal, applying selection criteria for perfect, strict, and loose hits against the Comprehensive Antibiotic Resistance Database (CARD ver. 3.2.8) (https://card.mcmaster.ca/, accessed in December 2023) [45]. Moreover, the amino acid sequences corresponding to quinolone resistance-determination regions (QRDRs) of gyrA, gyrB, parC, parE were extracted from these genomes. Subsequently, multiple sequence alignments were conducted using ClustalW within the MEGA X software. Amino acid residues were numbered according to the E. coli numbering system, and amino acid substitutions were interpreted by comparing them with QR isolates [25].

Genomic islands (GIs) were identified using the IslandViewer web tool (https://www.pathogenomics.sfu.ca/islandviewer/, accessed in December 2023) with four distinct methods: IslandPick, SIGI-HMM, IslandPath-DIMOB and Islander [46]. Plasmids were predicted using PlasmidFinder (ver 2.1) [47]. Prophages were detected using the default parameters of the PHASTER web server (https://phaster.ca/, accessed in December 2023) [48]. PHASTER evaluated DNA sequences by considering the count of coding sequences (CDS) and the presence or absence of the phage-associated genes. Insertion sequence (IS) families were examined using the ISFinder web tool (https://isfinder.biotoul.fr/, accessed in December 2023) [49]. CRISPR-Cas systems of the 22 isolates were identified utilizing CRISPRimmunity (http://www.microbiome-bigdata.com/CRISPRimmunity/, accessed on December 2023), which integrates three methods: CRISPRCasFinder, CRISPR Recognition Tool (CRT), and PILER-CR). The data regarding cas gene clusters of the isolates were extracted [50, 51].

Results

Phylogenetic analysis

A phylogenetic tree based on the 16S rRNA gene was constructed to explore the evolutionary relationship among 12 Thai CCB isolates from this study and an additional 13 Thai CCB isolates from previous studies of F. oreochromis, F. covae, and F. davisii. The result revealed four distinct clusters corresponding to the four species of CCB (Fig. 2). The Thai isolates were identified as F. oreochromis (n = 15) aligning with the type strain Costa Rica 04–02-TN^T, F. covae (n = 8) aligning with the type strain AL-02–36 ^T, and F. davisii (n = 2) aligning with the type strain 90–106 ^T. None of the Thai isolates were identified as F. columnare. The 16S rRNA sequences utilized in this study have been deposited in the GenBank database under accession numbers (Fig. 2, Table S1).

Fig. 2

Phylogenetic relationship of F. oreochromis, F. covae, F. columnare, and F. davisii generated from 16S rRNA sequences. Host species are represented by star and circular shapes, Flavobacterium species by triangles, and countries of isolates by rectangles. (Certain color references in this figure legend can be interpreted by the reader through the online version of this article)

Antimicrobial susceptibility tests by minimum inhibitory concentration (MIC) assay

The MIC values were ranged from 0.25 to > 4 μg/mL for OA, 0.015 to > 1 μg/mL for ENR, and 0.06 to > 16 μg/mL for OT (Table 2). Four F. oreochromis and two F. covae isolates exhibited decreased susceptibility to OA (MIC > 4 µg/mL), while one F. davisii isolate was no decreased susceptibility (MIC = 0.25 μg/mL). Five isolates of F. oreochromis were no decreased susceptibility to OA (MIC = 0.25 μg/mL). Additionally, five F. oreochromis isolates were no decreased susceptibility to ENR (MIC = 0.015 μg/mL), while four F. oreochromis isolates exhibited decreased susceptibility to ENR, with MICs ranging from 0.06 to > 1 μg/mL. Both F. covae isolates demonstrated decreased susceptibility to ENR, with a MIC greater than 1 μg/mL. Most isolates were no decreased susceptibility to OT, with MICs ranging from 0.06 to 0.125 μg/mL, while one F. covae isolate showed decreased susceptibility to OT with a MIC > 16 μg/mL. Subsequently, nine isolates representing three species were selected for next-generation sequencing based on their MIC values for OA, comprising OA non-wild-type isolates of F. oreochromis (KCRT2301, CNRT2201, RCBRT1901); OA non-wild-type isolates of F. covae (PCBSB2203); OA wild-type isolates of F. davisii (KCRT2007), and OA wild-type isolates of F. oreochromis (SKNRT2101, KCRT2304, UBRT2201, KCRT1906).

Table 2 The MIC values for quinolone, fluoroquinolone, and tetracycline antibiotics of F. oreochromis, F. covae, and F. davisii isolated from diseased farmed freshwater fish in Thailand

Genome features

A total of nine isolates, selected based on their OA MIC values from the twelve tested, were combined with 13 isolates from previous studies and the type strains of F. oreochromis (Costa Rica 04–02-TN^T), F. covae (AL-02–36 ^T), and F. davisii (90–106 ^T). These included seven F. oreochromis genomes (KCRT2301, CNRT2201, RCBRT1901, SKNRT2101, KCRT2304, UBRT2201, KCRT1906) and one genome of each F. covae (PCBSB2203) and F. davisii (KCRT2007) sequenced in the present study. Additionally, seven genomes of F. oreochromis (1214, CC1808, CC1805, CC1803, CC1802, 1201, NK01), five genomes of F. covae (SP1802, CF1, SP1809, 1362, 1372), and one genome of F. davsii (1215) from previous studies were included [8, 24]. The compiled sequences from these 22 Thai isolates exhibited notable variations in genome size. The genomes of F. oreochromis and F. davisii ranged from 3.3 to 3.5 Mbp, while those of F. covae ranged from 3.0 to 3.1 Mbp. The GC content of F. oreochromis and F. davisii genomes ranged from 29.5% to 30.0%, and 30.5% to 31.6%, respectively, whereas F. covae genomes averaged 30.5%. Moreover, F. oreochromis and F. davisii isolates exhibited a higher number of CDS compared to F. covae. F. oreochromis and F. davsii Thai genomes are generally consistent with the reference genomes in term of size, GC content, and CDS count. However, F. covae showed a smaller genome size and slightly different GC content when compared to the reference genome. Details regarding the genomic sequences of all Thai isolates are presented in Table 3.

Table 3 General genome features and assembly details of Thai F. oreochromis, F. covae, and F. davsii isolates and type strains used in genome comparative analysis

Core pan-genome analysis

The chart generated by the BPGA pipeline revealed the presence of 1,867 core genes (34.5%), 2,322 accessory genes (42.9%), and 1,221 unique genes (22.6%). The relatively low count of core genes observed across these isolates can be attributed to their genetically diverse nature. Moreover, F. davisii KCRT2007 possessed the highest number of strain-specific genes (346 genes) followed by F. davisii 1215 (257 genes), while F. oreochromis CC1802, CC1803, CC1805, CC1808 and F. covae SP1802, SP1809 had none. The distribution pattern of strain-specific genes among 22 genomes is illustrated in Fig. 3A. A classification using Clusters of Orthologous Group (COG) analysis revealed that core genes were primarily categorized into general function prediction only (R), amino acid transport and metabolism (E), and translation, ribosomal structure, and biogenesis (J). In addition, unique genes were predominantly found in categories such as general function prediction only (R), defense mechanisms (V), replication, recombination, and repair (L), transcription (K) as well as cell wall/membrane/envelop biogenesis (M) (Fig. 3B). In the context of pan-genome, the total gene of pan-genome would initially consist of 5,410 genes, with an expected continuous growth with the addition of new genomes. This growth follows a power law equation “f(x) = a × x^b” where the expansion rate (b) is calculated to be approximately 0.199627 (Fig. 3C). For core genome based phylogenetic tree, the isolates were classified into three groups, which corresponded to the three distinct species. Flavobacterium davisii and F. covae were the most closely related, whereas F. oreochromis showed the most distant relationship with the others (Fig. 3D).

Fig. 3

Pan-genome analysis from twenty-two isolates of F. oreochromis, F. covae, and F. davisii. A The number of core genes, accessory genes, and unique genes presented in each isolate. B COG functional categories of core genes, accessory genes, and unique genes among twenty-two isolates. C Core pan-genome curves of power-fit and exponential are indicated by brown lines, and purple lines. The corresponded equations are provided within the figure legend. The blue dots indicated total gene families; the pink dots indicated core gene families. D Phylogenetic analysis based on the core gene sequence of twenty-two isolates of F. oreochromis, F. covae, and F. davisii. (Certain color references in this figure legend can be interpreted by the reader through the online version of this article)

A whole genome phylogenetic analysis was conducted on 22 Thai isolates and seven reference isolates using the TYGS pipeline. Flavobacterium covae, F. columnare and F. davisii were closely clustered, while F. oreochromis diverged into a separate cluster (Fig. 4). These findings supported the core genome-based classification described above (Fig. 3D).

Fig. 4

Phylogenetic analysis of twenty-two isolates of F. oreochromis, F. covae, and F. davisii based on whole genome sequences. Reference isolates, including F. oreochromis (Costa Rica 04–02-TN^T), F. covae (AL-02–36 ^T, 94–081, C#2), F. davisii (90–106 ^T), and F. columnare (ATCC 49512) were used, with F. psychrophilum JIP02/86 as an outgroup

Comparative genomic analysis

RAST subsystem categorization

According to the SEED subsystem classification, all isolates showed a higher number of genes related to “amino acids and derivatives”, “cofactors, vitamins, prosthetic groups, pigments”, and “protein metabolism” (Fig. 5, Table S2). F. oreochromis exhibited the highest number, whereas F. covae and F. davisii showed the lowest number of genes in categories such as “Amino Acids and Derivatives”, “Sulfur Metabolism”, “Virulence, Disease and Defense”, “Nitrogen Metabolism”. F. davisii KCRT2007 had more “Cell Wall and Capsule”, and “Stress Response” genes than other isolates.

Fig. 5

SEED subsystem categorization of F. oreochromis, F. covae, and F. davisii genomes. Each colored bar indicated the number of genes assigned to specific categories

Virulence gene analysis

A total of 212 potential virulence genes were predicted across fourteen F. oreochromis, six F. covae, and two F. davisii isolates using the VFDB database based on E-values < 10^–4, with identities generally exceeding > 40%. The distribution of virulence genes varied across all isolates; however, a similar distribution within each CCB species was observed (Fig. 6). The 212 virulence genes were grouped into 13 categories based on their function, in which 6 categories (nutritional/metabolic factor, immune modulation, adherence, effector delivery system, motility, and exotoxins) collectively constituted 89.15% of total, while the remaining 7 categories, related to biofilm, stress survival, and other functions, accounted for the remaining 10.85% (Table S3). These categories show primary mechanisms of pathogenicity in these bacteria.

Fig. 6

Distribution of putative virulence factors among twenty-two Thai isolates of F. oreochromis, F. covae, and F. davisii

Furthermore, 88 out of 212 genes were shared by all isolates, while 37 genes were unique to F. oreochromis, six were unique to F. covae, and 23 were unique to F. davisii. Notably, capsule, lipopolysaccharide (LPS), and lipooligosaccharide (LOS) genes related to immune modulation were identified. Flavobacterium oreochromis showed a greater number of capsule genes and LPS genes compared to the other species. Several unique LPS genes were present in F. oreochromis (neuC1, rfbM), F. covae (wbaP/rfbP, kdtB, vexE), and F. davisii (kdsB). Unique LOS genes were identified in F. oreochromis (neuB1) and F. davisii (msbA, galU). Additionally, the exotoxin-related gene cylG was presented in all isolates, while hlyB was found only in F. davisii KCRT2007. Type IV pili were identified in all isolates as a major within the category of adherence. Type III, IV, VI secretion systems involved in the effector delivery system category, were also found in all isolates. Several unique genes annotated as flagella genes were present in F. covae (cheA, motB), F. davisii (pomB2, flrC), and F. oreochromis (fliA). icaA biofilm-related genes were present in all F. oreochromis isolates. The stress survival related genes, sodA and sodB, were identified in all isolates.

Secretion systems distribution

MacSyFinder and TXSScan were utilized to identify genes encoding protein components of various secretion systems. T1SS, T3SS, T6SS, and T9SS components were predicted in all isolates and the copy number of mandatory elements were varied among isolates. Flavobacterium oreochromis and F. davisii exhibited a high copy number of genes associated with the T1SS. Meanwhile, F. covae 1372 and F. oreochromis SKNRT2101 displayed a unique gene pattern in the T4SS (Fig. 7, Table S4).

Fig. 7

Secretion systems distribution among twenty-two Thai isolates of F. oreochromis, F. covae, and F. davisii

Analysis of ARGs and detection of mutations in QRDRs

A total of 195 well-characterized Antibiotic Resistance Genes spanning multiple drug classes were identified and classified (Fig. 8A). The predicted mechanisms underlying ARGs involved alterations in antibiotics target alteration, antibiotics target replacement, antibiotics protection, antibiotic inactivation, and antibiotics efflux pump (Table S5). An abundance of genes encoding resistance-nodulation-division (RND) superfamily were identified compared to other efflux families. Genes involved in resistance to a wide range of antibiotics, including fluoroquinolone, glycopeptide, peptide, tetracycline, cephalosporin, and multidrug resistance were observed in all isolates. Cluster analysis revealed that all isolates were grouped into three clusters, corresponding to the three CCB species. Each cluster correlated with specific drug groups, showing significant variation in the number of ARGs conferring drug resistance, indicating distinct patterns. Interestingly, the cluster of F. oreochromis displayed the highest number of unique genes associated with glycopeptide, fluoroquinolone, and cephalosporin resistance. Conversely, F. davisii exhibited unique genes related to streptogramin and macrolide resistance, while genes unique to F. covae were predominantly linked to fluoroquinolone and tetracycline resistance. Indeed, F. oreochromis and F. covae showed abundance of unique genes associated with resistance to fluoroquinolone, while F. covae showed a significant number of unique genes related to tetracycline resistance, correlating with elevated MIC to ENR, OA, and OT (Table 2). The QR-associated genes among all isolates consisted of 23 genes in F. oreochromis, 28 genes in F. covae, and 25 genes in F. davisii (Fig. 8B). The percent identity of these genes in each isolate was compared to QR-associated amino acid sequences curated in CARD database using protein blast, and detailed results are provided in Table S6.

Fig. 8

A Distribution of antimicrobial resistance genes predicted among twenty-two Thai isolates of F. oreochromis, F. covae, and F. davisii. B the percentage amino acid identity of quinolone resistance genes from F. oreochromis, F. covae, and F. davisii isolates compared to the Comprehensive Antibiotic Resistance Database. (Certain color references in this figure legend can be interpreted by the reader through the online version of this article)

Moreover, the QRDR sequences comparison of all tested isolates showed amino acid substitutions within gyrA at codon 83 in multiple isolates, including Ser83Val in CNRT2201, Ser83Ala in NK01, 1201, 1214, RCBRT1901, 1215, and Ser83Phe in KCRT2301, PCBSB2203, SP1802, and SP1809. Additionally, an amino acid substitution was observed at codon 87 in PCBSB2203 (Asp87Tyr). No significant amino acid substitutions were detected in gyrB, parC, or parE across all isolates (Table 4). Furthermore, the genome of F. covae PCBSB2203 contains two genes encoding tetracycline resistance (tetA_1,2), specifically the tetA_2 gene, which correlates with the observed decreased phenotypic susceptibility to oxytetracycline. In contrast, other isolates with only a single tetA_1 gene did not exhibit the same level of decreased susceptibility.

Table 4 Mutations detected in the quinolone resistance-determining regions of gyrA, gyrB, parC, parE from type strains and twenty-two Thai F. oreochromis, F. covae, F. davsii isolates

Genomic islands, prophage prediction, insertion elements and CRISPR system detection

To further investigate the genetic diversity of these strains, IslandViewer4 was utilized to identify GIs using multiple prediction methods (Fig. 9). F. oreochromis and F. covae isolates contained between 5 to 18 GIs and 10 to 15 GIs, respectively. Flavobacterium davisii 1215 harbored 12 GIs, while KCRT2007 contained 22 GIs. These predicted GIs across species were mainly comprised of genes encoding hypothetical proteins, tryptophan synthase, and transcriptional regulators. Notably, GIs of F. oreochromis isolates exhibited unique genes encoding Ribosome-binding factor, Ribose-phosphate pyrophosphokinase, ADP-ribosylarginine hydrolase, O-acetyl-ADP-ribose deacetylase, nicotinate-nucleotide adenylyltransferase, phosphoglycerate kinase, histidine decarboxylase, and other proteins. Similarly, GIs of F. covae contained unique genes encoding peptidoglycan O-acetyltransferase, lipoyl synthase, methylglyoxal synthase, biosynthetic arginine decarboxylase, glyceraldehyde-3-phosphate dehydrogenase. GIs of F. davisii contained unique genes encoding serine-tRNA ligase, L-asparaginase, adenine DNA glycosylase, and other proteins. Surprisingly, some resistance genes were found within the GIs of certain isolates. The tuf gene was identified within the GIs of SP1809, CF1, 1215, KCRT1906, PCBSB2203, 1362, 1372, SP1802. Similarly, the tetA, emrE, bcglB genes were identified in GIs of CF1. The gyrA and vat genes were present in SP1802 and SP1809. Additionally, several virulence genes were observed within the GIs of multiple genomes, including clpP in 1215, CF1, KCRT2007, SKNRT2101, 1214, and KCRT2301, lpxD in KCRT2007, msrA in 1214, clpB in KCRT1906, CF1, PCBSB2203, 1372, RCBRT1901, SP1802, 1362, and SP1809 (Fig. 9, Table S7). No plasmids were detected in any isolate.

Fig. 9

Analysis of genomic islands in twenty-two Thai isolates of F. oreochromis, F. covae, and F. davisii isolates. Representation of prediction methods: Red (integrated methods), Orange (SIGI-HMM), Blue (IslandPath-DIMOB). Identification of incomplete phage components using the PHASTER tool. (Certain color references in this figure legend can be interpreted by the reader through the online version of this article)

To ascertain the presence of phages, the PHASTER tool identified one to two prophage sequences in each genome. The prophages in the 22 isolates received scores between 10 to 50, indicating incomplete and questionable prophage regions. These incomplete prophage regions, ranging in length from 6.3 kb to 32.9 kb and one questionable prophage region measuring 21 kb, were distributed throughout different locations within each genome. Notably, F. davisii KCRT2007 exhibited the largest prophage region, which spans 32.9 kb and contains the attR gene. The attR gene encoding protein related to generate recombinant junctions with assistance of the integrase gene. Most isolates harbored incomplete prophages containing hypothetical proteins, phage-like proteins, and other proteins across all regions. Additionally, prophages in some isolates contained genes encoding coat protein, fiber protein, transposase protein, protease protein, tail shaft protein (Fig. 9; Table S8).

Insertion elements play a crucial role in understanding genome organization and evolution. A total of 17 distinct insertion sequence (IS) families were identified (Fig. 10, Table S9). Among these families, IS3 was the most prevalent. All isolates contained IS3 families, with F. davisii 1215 showing the highest number of IS3 families (n = 20). Annotation of insertion elements across all genome isolates revealed that F. oreochromis RCBRT1901 harbored the highest number of IS elements. Specifically, IS5 was the most abundant in RCBRT1901, with a count of 55. Interestingly, the number of IS elements in F. covae and F. davisii genomes exceeded those in F. oreochromis.

Fig. 10

Insertion elements in 22 Thai isolates of F. oreochromis, F. covae, F. davisii

Another notable aspect of these genomes was the presence of a CRISPR/Cas system. Using CRISPRminer2 web tool to predict the CRISPR arrays and cas genes. The CRISPR regions in all isolates contained multiple spacer elements, ranging from 1 to 55. Notably, F. davisii (1215 and KCRT2007) lacked cas genes flanking the CRISPR regions. On the other hand, F. oreochromis and F. covae carried various cas genes, including cas1, cas2, cas9, and cas13b. The CRISPR types associated cas protein of each isolate were identified and are presented in Table S10.

Discussion

Flavobacterium oreochromis, F. covae, and F. davisii have emerged as the primary pathogens responsible for columnaris disease in Thailand, leading to significant economic losses in the freshwater fish industry [3,4,5,6, 8, 24, 52]. In this research, a comprehensive comparative analysis was undertaken on the complete genomes of Thai F. oreochromis, F. covae, and F. davisii isolates.

The phylogenetic analysis of 25 Thai CCB isolates revealed that all Thai isolates were identified as F. oreochromis, F. covae, or F. davisii. Moreover, the core genome and whole genome phylogenetic trees were supported by trees constructed from previous phylogenetic analyses based on multilocus sequence analysis [24, 53]. Flavobacterium covae and F. davisii showed higher genomic similarity to each other than to F. oreochromis, potentially indicating a closer evolutionary linkage between these two species. F. oreochromis may have diverged earlier or undergone more genetic divergence than F. covae and F. davisii [54]. Otherwise, the presence of F. oreochromis and F. covae in red tilapia, Asian seabass, snakeskin gourami, and striped catfish underscores the potential for inter-species transmission. Previously, host associations were identified in CCB [53]. Flavobacterium columnare, previously classified as genetic group 1 (GG1), was more frequently recovered from infected cold water fish families such as Salmonidae and distributed in North and South America and Europe. F. covae (previously GG2) and F. davisii (previously GG3) were commonly reisolated from infected fish of the family Ictaluridae and reported in North America, Africa and Asia. F. oreochromis (previously GG4) was mostly recovered from tilapia (family Cichlidae), found in South America, Central America and Asia [53]. A recent study testing the virulence of CCB species in catfish revealed high mortality rates in fish infected with F. covae and F. oreochromis despite F. oreochromis not being known to occur in the United States [9]. The findings highlighted that F. oreochromis could potentially cause mortality rates comparable to those observed with F. covae. In contrast, lower mortality rates were recorded with F. columnare and F. davisii. Mortality rates increased when water temperature reached 27 °C during the experimental challenge, suggesting that relying solely on the number of virulence genes is insufficient for comprehensive conclusions. However, this is a preliminary step toward understanding Flavobacterium virulence. Understanding these variations can inform the development of strategies to control and manage Flavobacterium infections.

Furthermore, comprehensive comparisons regarding antibiotic resistance are limited among these species. MIC method was used to screen the difference in susceptibility patterns in F. oreochromis, F. covae, and F. davisii. According to the study by Chokmangmeepisarn et al. [8], F. oreochromis isolates were susceptible to both OA and ENR, while F. covae isolates were only susceptible to ENR. Since OA is a first-generation quinolone and ENR is a second-generation fluoroquinolone, our study suggests that F. oreochromis may be developing decreased susceptibility to fluoroquinolones. The F. oreochromis isolates in the previous study were collected in 2018, whereas those in our study, isolated between 2019 and 2023, showed decreased susceptibility to the second-generation fluoroquinolone, ENR. Moreover, this study identified a highly decreased susceptibility F. covae isolate that exhibited decreased susceptibility to all tested drugs, including OA, OT, and ENR. This represents the first report of decreased phenotypic susceptibility to tetracycline in this species in Thailand. Together, these findings underscore the increasing concern over antimicrobial resistance in these bacteria. Currently, OA, ENR, and OT are widely used in aquaculture, and consequently, there is a possibility of reduced efficacy of these antibiotics for treating columnaris disease [8, 52, 55]. A total 195 genes were identified as ARGs through protein homology blast analysis against CARD. The primary resistance mechanisms were identified included antibiotics target alteration, efflux pump activation, and antibiotics inactivation. RND associated with antibiotics efflux pump was predominantly discovered in three species, comparing to other resistance mechanisms. Efflux pumps are protein complexes spanning bacterial cell membranes, play a crucial role in reducing the accumulation of drugs within the bacteria cytoplasm, thereby contributing to multidrug resistance [56]. RND efflux was predominant in three species, leading to resistance against multidrug phenotypes as previous reported in F. columnare Pf1 [57]. Identified ARGs in this study indicated their widespread presence in the environment [8, 52]. Moreover, ARG clustering revealed that the three species indicated distinct drug resistance patterns. Surprisingly, the tetracycline resistance gene (tetA_1) was distributed across these species. However, only one isolate of F. covae PCBSB2203, possessed tetA_2 gene and exhibited decreased phenotypic susceptibility to OT. Tetracycline resistance has been observed in Flavobacterium spp. isolates from ornamental fish [17]. Class A tet determinants can confer high-level tetracycline resistance [58]. Tetracycline resistance tetA is encoded in the genome of HV carp isolate 04017018, which exhibited acquired phenotypic resistance to OT. However, this gene was absent in all other isolates that were susceptible to OT [59]. Further analyses are needed to investigate the phenotype and genetic mechanisms associated with this gene.

Nine isolates of F. oreochromis, F. covae, and F. davisii were selected for analysis of quinolones resistance and susceptibility, based on the amino acid substitution within the QRDRs. The findings indicated the highest presence of QR-associated genes in F. covae, demonstrating their influence in reducing susceptibility to quinolone. Several studies have suggested that four genes gyrA, gyrB, parC, and parE have a significant impact on the phenotype of quinolone resistance [8, 52, 60]. Among F. oreochromis and F. covae, mutations were found in the QRDR of gyrA (Ser83Phe, Ser83Val, Ser83Ala; Asp87Tyr), resulting decreased phenotypic susceptibility to quinolones and fluoroquinolones as documented in previous study [8, 52]. All OA non-wild-type isolates of F. covae exhibited the Ser83Phe amino acid substitution. Conversely, the Ser83Val, Ser83Phe, Ser83Ala substitutions were consistently found in isolates of F. oreochromis, while the Ser83Ala substitution was found specifically in F. davisii. These findings indicate distinct resistance mechanisms across different Flavobacterium species. Mutations in gyrA, specifically at codon 83, have been identified as the primary driver for the development of QR phenotypes in Gram-negative bacteria [52, 61]. Additionally, mutations at codon 87 may potentially contribute to a significant reduction in susceptibility to OA. Extensive studies are required to confirm this. Previous studies have shown that mutations at positions 83 and 87 in the gyrA gene are associated with high levels of antibiotic resistance (MIC > 16 µg/mL) in isolates from Vietnam [52, 62, 63].

The Gln458His mutation in gyrB was observed in all isolates of F. covae and F. davisii. Additionally, another mutation, Asn439Ala, was present in all F. davisii isolates, potentially indicating species-specific genes. In the previous study, F. oreochromis showed variable amino acids in the QRDR of parC and displayed both low and high levels for all tested quinolones. Thus, mutations in this position of parC gene were not directly linked to QR [8]. However, previous studies have indicated that mutations in both gyrA and parC genes confer a moderate level of resistance to fluoroquinolones [63]. In our investigation, isolates with mutations in both gyrA and parC exhibited decreased phenotypic susceptibility to OA and ENR. Therefore, while a single mutation in gyrA was adequate to induce a loss of fitness, the presence of double mutations in gyrA or a single mutation in gyrA coupled with a single mutation in parC significantly increased resistance in F. oreochromis and F. covae. Mutations in gyrB and parC were found to be less influential, as observed in studies of Enterobacteriaceae and F. psychrophilum [61, 64]. These results suggest that gyrA is the major target of the quinolone in F. oreochromis, F. covae, and F. davisii, while gyrB, parC, parE might not be directly related to QR phenotype. Variation in codons 83 and 87 of gyrA gene show a significant influence on QR in bacteria by changing the structure of the DNA gyrase protein, which is the target of quinolone antibiotics. Specially, the ser83 of the DNA gyrase protein is involved in forming an important interaction with the keto acid of the quinolone compound. This interaction is essential for quinolones to inhibit DNA gyrase activity [8, 52, 61]. In Gram-negative bacteria, mutations occur more often in DNA gyrase, whereas in Gram-positive bacteria, the primary target is topoisomerase IV, which consists of the ParE and ParC subunits [65, 66]. Consequently, bacteria carrying these mutations exhibit reduced susceptibility to quinolone antibiotics, thereby contributing to the development of quinolone resistance. Further studies are needed to determine whether specific protein-altering mutations contribute to decreased susceptibility.

Flavobacterium oreochromis isolates demonstrated a notable presence of genes associated with various metabolic pathways and defense mechanisms, which may contribute to their predominance in freshwater fish in Thailand. Specifically, these isolates showed a high number of genes involved in sulfur-containing compound metabolism, virulence, disease defense, and nitrogen metabolism compared to the other species. The active involvement of F. oreochromis in sulfur metabolism may be essential for cellular processes and environmental adaptation, similar to bacteria found in deep-sea hydrothermal vents, where bacteria utilize reduced sulfur compounds, such as hydrogen sulfide, as an energy source for chemosynthesis [67, 68]. This metabolic capability may provide F. oreochromis with a competitive advantage in diverse and potentially anaerobic environments, such as those found in freshwater systems. Seventy genes in F. columnare ATCC 49512 related to nitrogen metabolism, reflect its self-sufficient nutrient utilization, with the capability for denitrification enabling potential anaerobic growth in aquatic pond sediments [69, 70]. Additionally, the presence of genes associated with virulence, disease, and defense suggests that F. oreochromis has evolved mechanisms to enhance its pathogenicity and defense against environmental stressors or host immune responses [71]. Differences in virulence factors, including those related to fimbriae, flagella, toxins, various secretion systems, and iron ion uptake systems, are likely key for the heightened virulence observed in Aeromonas veronii [72]. Therefore, the genetic adaptations of F. oreochromis in sulfur and nitrogen metabolism, combined with enhanced virulence and defense mechanisms, likely contribute to its predominant isolation from freshwater fish in Thailand compared to other species [6, 8, 24]. Moreover, F. davisii KCRT2007 may possess an enhanced capacity for cell wall and capsule synthesis, and stress response mechanisms compared to other isolates. Genes involved in capsule synthesis in Streptococcus pneumoniae plays an important role in both colonization of nasopharyngeal cavity, where reduced capsule amounts facilitate survival, and in systemic infections like pneumonia and bacteremia, where high capsule levels are essential for evading host immune responses [73].

Additionally, the identification of pathogenic virulence factors is significant for understanding bacterial pathogenesis and identifying potential vaccine targets [74]. The application of comparative genomics in predicting virulence factors allows for the comparison of specific virulence-associated genes across species, thereby advancing our understanding of bacterial pathogenicity and host interactions [75, 76]. Virulence genes responsible for iron acquisition and utilization, vitamin and cofactor biosynthesis, peptide synthesis, and ABC transporter systems were classified into the nutritional/metabolic factor category, contributing to the bacterial growth, survival, and metabolic host–pathogen interactions [77]. In this study, genes related to iron acquisition systems were present in all species. Iron acquisition systems are crucial for F. columnare, relies on iron for virulence and metabolism [78, 79]. These systems include TonB-dependent and ferroxidase genes. Previous research showed significant variations in the expression of these genes among two F. columnare isolates with different levels of virulence. This upregulation correlated with increased virulence under iron-limiting conditions, indicating that iron acquisition systems play an important role in virulence differences among F. columnare isolates [79]. Genes for iron uptake components, including hmuY, outer membrane receptors (fhuA, fhuE, fecA), and an ABC transporter, have been identified. Deletion mutants showed growth defects under iron-limited conditions, though some retained wild-type-like virulence. However, mutants lacking in multiple iron uptake components exhibited decreased virulence [80]. This suggests that all species may possess a robust iron acquisition capability, contributing to its virulence and survival under iron-limited conditions. Understanding and targeting the different iron acquisition systems in each species could be a focus for further investigation.

Flavobacterium oreochromis exhibited a greater abundance of capsule genes associated with immune modulation compared to other species. The capsule is vital for F. columnare adherence to gill tissue and aids in identifying extracellular chemoattractant from catfish mucus [81, 82]. LPS and LOS genes were particularly abundant in all three species, serving as a major component of the outer membrane, known for significance in the host immune defenses and virulence of Gram-negative bacteria [83,84,85]. Analysis of LPS through immunoblotting demonstrated that the avirulent mutant lacked the presence of higher molecular bands within LPS profile. This result enabled to identify the attenuated mutant from other strains of F. columnare [86, 87]. The presence of different LPS genes in each species may contribute to variations in the structure of the outer membrane and their interactions with the environment. LOS plays a role in enabling Histophilus somni to evade host defense [88]. Variations in LOS genes serve as a significant mechanism for variation in LOS among various strains of Campylobacter jejuni [89]. These modifications of LOS structure led to the synthesis of new LOS structures that enhance the survival of pathogen [90].

Furthermore, biofilms play an important role as protective and adhesive structures, enabling bacteria to attach to surfaces, resist environmental stresses, evade host immune system, and enhance antibiotic resistance, thus enhancing survival and resilience in diverse environments [91]. The presence of icaA biofilm-related genes was detected in all F. oreochromis isolates, indicating their potential significance in biofilm formation similar to that observed in Staphylococcus aureus, especially under conditions involving NaCl exposure, as evidenced by a marked decrease in biofilm formation in an icaA-deletion mutant compared to the wild type of strain [92]. Moreover, sodA and sodB were identified in all isolates of three species. These genes play an important role in resisting oxidative damage, allowing the bacteria to survive within fish macrophages and evade further infection [93]. Surprisingly, motB was also present in all F. covae isolates, while motY was found in all F. oreochromis and F. davisii isolates. These genes are known components of the flagellar motor, essential for driving flagella rotation [94]. Similarly, two MotB-related proteins were identified in F. columnare strain Pf1, despite the absence of the flagella-encoding genes in the genome. This finding suggests that these proteins may play a role in movement or serve other functions [70]. Flavobacterium oreochromis and F. davisii exhibited a higher number of virulence genes in comparison to F. covae, showing a significant enrichment in genes associated with adherence, effector delivery systems, immune modulation, and motility.

Secretion systems play a crucial role for bacterial metabolism and pathogenicity, contributing significantly to processes such as host invasion, evasion of host defenses, tissue damage, and competition with other bacterial species [95, 96]. The T1SS is common among Gram-negative bacteria and aids in directly transporting substrates into the extracellular space through the periplasm, comprising three essential proteins that span the cell envelope [97]. The T3SS enables the delivery of effector molecules into host cells, allowing bacteria to manipulate various host cell responses, including signal transduction, cytokine production, and cytoskeleton rearrangement, thereby enhancing their pathogenicity [98, 99]. T6SS is a multifunctional system that may play roles in virulence, interbacterial competition, nutrient acquisition, and horizontal gene transfer [100]. The number of TssD (Hcp) and TssH (ClpV) encoded in three species varied. Strain NK01 exhibited higher number of Hcp and ClpV proteins compared to the other strains, potentially contributing to its pathogenicity in fish. Moreover, previous research found that the genome of a highly virulent salmonid isolate contained a unique gene encoding the TssI (VgrG) protein, which was absent in non-virulent isolate. VgrG is a key component of the T6SS, plays an important role in virulence [59]. The gld, por, and spr genes encoded proteins related to the T9SS, facilitating adhesion, the secretion of surface motility proteins, and the pathogenesis of F. johnsoniae and F. columnare [101, 102]. Alternatively, gldL and gldM is implicated in both motility and secretion, as the absence of gldL and gldM resulted in failed gliding and secretion processes [103, 104]. Moreover, deletion of porV displayed avirulent phenotypes of F. columnare [105]. Thus, T1SS, T3SS, T6SS, T9SS secretion systems presented in all of F. oreochromis, F. covae and F. davisii isolates are a major function in facilitating the transport of proteins across species, contributing to virulence by delivering toxins and effector proteins into host cells, enabling communication between bacterial cells, interacting with other organisms in the environment, and promoting adaptation to diverse ecological niches [71].

The pili genes in T4aP system are essential in signaling, motility, and adhesion [106, 107]. T4SS facilitates the transfer of mobile genetic elements among bacteria. T4SS is capable of releasing both proteins and DNA into host cells. To the date, among several types of secretion systems of different Gram-negative bacteria, only the T4SS allows the enlargement of DNA into host cytosol. This capability is significant for processes such as virulence factors transfer, bacterial conjugation, and the manipulation of host cellular functions [108,109,110]. However, the absence of the T4SS in all Flavobacterium genomes was reported in a previous study [69]. Interestingly, the presence of T4SS related to the high number of insertion elements detected in some of A. veronii [111]. However, isolates F. oreochromis SKNRT2101 and F. covae 1372 exhibited a high number of T4SS without increasing in insertion elements. Unique T4SS and T4aP genes in F. covae and F. oreochromis likely contribute to species-specific traits or adaptations. Variation in T4SS may affect pilus biogenesis and adherence to host cell receptors, which may influence the virulence of certain hosts [95]. These variations in secretion systems are believed to facilitate adaptation and survival across diverse habitats [112]. Research on the role of secretion systems, other than the T9SS, in Flavobacterium virulence is limited, underscoring the need for further studies to elucidate the functions of the proteins involved in these secretion systems.

Furthermore, the CRISPR-Cas system acts as a bacterial immune mechanism, enabling bacteria to identify and eliminate foreign genetic materials such as bacteriophages, and plasmids by storing short segments of non-hybridized DNA sequences in bacterial genome as “memory”, allowing Cas nucleases to precisely target and differentiate between matching sequences when regenerating [113,114,115]. Putative cas genes was detected in all isolates, except for the two F. davisii isolates (1215 and KCRT2007). This observation may explain the higher prevalence of phages and genomic islands in F. davisii isolates. These findings suggest that the two F. davisii isolates may lack the ability to utilize the CRISPR-Cas system for adaptive immunity against viruses or plasmids. Without the necessary CRISPR-Cas proteins to cleave foreign DNA, this species is more allowable to viral infections and horizontal gene transfer events [116]. Streptococcus pyogenes strains exhibited higher prophage counts with lacking CRISPR-Cas systems [117]. Moreover, encapsulated strains of S. pneumoniae were genetically modified by introducing CRISPR sequences targeting capsule genes into nonencapsulated strains, rendering them avirulent. These engineered strains exhibited almost no infection in mice. However, loss of CRISPR function in some engineered strains allowed for the acquisition of capsule genes, resulting in successful infections in mice. The phenomenon of CRISPR loss in bacterial pathogens may occur under intense selective pressure favoring virulence or antibiotic resistance [116, 118]. This could potentially explain the absence of CRISPR in F. davisii isolates, facilitating their survival. The presence or absence of CRISPR-Cas may correlate with environment variables. In addition, optimal growth temperature is correlated with the number of spacers [119, 120]. The elevated spacer count observed in F. oreochromis correlates with the wide temperature tolerance, ranging from 15 °C to 40 °C [13].

Conclusions

This comprehensive genomic analysis of F. oreochromis, F. covae, and F. davisii showed significant genetic variability that has shed light on the evolutionary paths and adaptive strategies. The phylogenetic analysis demonstrated that F. covae and F. davisii are closely related, while F. oreochromis diverged significantly, occupying a distinct branch. The unique genes of each species might enable them to adapt to specific environmental conditions and host interactions. A notable discovery was the absence of CRISPR systems in all F. davisii isolates, which might explain the high number of genomic islands and the presence of large phages observed in these isolates. A critical mutation at Ser83Phe, Ser83Val, Ser83Ala, and Asp87Tyr in the resistance-determining regions of gyrA played a significant role in antibiotic resistance. Moreover, an F. covae isolate exhibited decreased phenotypic susceptibility to OT, potentially due to the presence of the tetA_2 gene, which was absent in all OT wild-type isolates. These findings emphasize the importance of carefully considering the use of quinolone and fluoroquinolone antibiotics in treating infections caused by these bacteria, while tetracycline remains a viable option. However, further research is necessary to explore the pathogenic mechanisms of these species across different hosts, as well as the roles of virulence genes and mutations in antibiotic resistance. This study provides valuable insights into the genetic, evolutionary, and pathogenic diversity of these species, contributing to the development of effective strategies for managing columnaris disease in freshwater fish species, protecting aquatic ecosystems and mitigating associated human health risks in Thailand.