A total of 100 distinct M. marinum genome sequences (including contigs, scaffolds, chromosomes, and complete sequences) are currently available in the NCBI genome database. On September 12, 2024, 8 complete genomes of M. marinum were downloaded. To further refine the data, the CheckM program (version 1.2.1) (Parks et al., 2022) was employed for genome quality filtration. The criteria for genome completeness were set at a minimum of 90%, while the maximum allowed contamination was 5%. Recent studies have emphasized the importance of selecting a high-quality and appropriate dataset for pangenome research (Wu et al., 2021; Wu et al., 2022). This careful selection ensures the reliability of downstream analyses, leading to more robust insights into pangenome and potential SNPs within the species. Supplementary Table 1 presents genome assembly, accession numbers, and genomic information, including genome size, gene count, and protein-coding genes. Only whole genomes were analyzed to ensure the most accurate representation of strain virulence genes and potential identifying biomarkers.

To verify species relationships, we performed nucleotide-level comparisons of all potential genome combinations using the ANI approach. The Python script PyANI v0.2.722 was used to compute pairwise ANI values using two different methods, the MUMmer (Kurtz et al., 2004) and the BLAST + method (Camacho et al., 2009) employing the ANIm and ANIb options, respectively. All strains classified as M. marinum were chosen for further analysis, while those misidentified were discarded. This decision was based on setting a threshold value indicating that the ANI of the same species exceeds 94% (Goris et al., 2007; Richter and Rosselló-Móra, 2009).

To maintain consistency and standardization, the genomes of the M. marinum strains examined were re-annotated for their functional attributes using the Prokka suite version 1.14.6 (Podrzaj et al., 2022). The software Roary (Page et al., 2015) and Panaroo (Tonkin-Hill et al., 2020) with a default identity threshold of 95%, was employed to construct the core and pan-genome of M. marinum. Briefly, genes have been categorized into four distinct groups: core genes, present in more than 99% of genomes; softcore genes, found in most strains, specifically in over 95% but fewer than 99%; shell genes, identified in 15% to less than 95% of the strains; and cloud genes, which are observed in fewer than 15% of the total strains. In addition, the pangenome was analyzed using Heap’s Law (Rajput et al., 2023) and Power Law fit (Tettelin et al., 2005) in conjunction with the custom script¹ and the “ggcaller v1.3.0” in Python programing (Horsfield et al., 2023), respectively. This approach allowed for the calculation of constant variables and the application of the least squares method to fit an exponential regression decay model to both the core genome and singletons. Heap’s Law was used to determine the fixed parameters from the pangenome with the formula n = κN^γ, where n represents the number of pangenome genes and N is the number of genomes, least-squares fit is represented by the equation n = k × exp [-x/t] + tgθ, where n is the number of genes, and k, t, and tgθ are independent variables (Islam et al., 2025; Mwamburi et al., 2024). This algorithm is utilized to ascertain the number of genes that will compose the core genome upon stabilization (Soares et al., 2013) and to offer an approximate calculation of the genes contributed by each freshly sequenced genome. Additionally, the complete pangenome for the entire dataset of genomes was generated using anvi’o v8 (Eren et al., 2015), followed by the pangenomics procedure described in this guide: https://merenlab.org/2016/11/08/pangenomics-v2/2/. In brief, the following scripts were executed: anvi-gen-contigs-database, utilized to create a database, with Prodigal v2.6.3 (Hyatt et al., 2010) for identifying open reading frames in contigs, and anvi-run-ncbi-cogs, which provided gene annotations by utilizing NCBI Clusters of Orthologous Groups database (Tatusov et al., 2000). The genome database was created using anvi-gen-genomes-storage and anvil-pan-genome for visualization. Anvi’o employs the DIAMOND tool to compute the similarity between each amino acid sequence in every genome and all other amino acid sequences across all genomes in the dataset and then uses the MCL algorithm (van Dongen and Abreu-Goodger, 2012) to detect clusters in the results of amino acid sequence similarity.

Gene synteny analysis was performed using the Mauve program version snapshot_2015_02_13 (Edwards and Holt, 2013) and its progressive Mauve algorithm to detect potential gene rearrangement events. The software segmented the genomes into predefined fragments and used this information to perform numerous genome alignments. These alignments identified local collinear blocks (LCB) which were then visualized in a figure. This figure serves as a tool to identify and analyze rearrangements within the genomes. Subsequently, a customized Python script called msa2snp.py² was employed to detect all SNPs within the target genes. The results of SNP identification were manually curated to identify the most informative SNP, found exclusively in all M. marinum complete genomes.

BLAST searches were performed for all protein-coding genes in the National Centre for Biotechnology Information (NCBI) database, which was accessed at https://www.ncbi.nlm.nih.gov/on 20.09.2024. The entire genome, including repetitive elements, was annotated using eggNOG-mapper v2 (Cantalapiedra et al., 2021) and InterProScan (Quevillon et al., 2005) in a Linux environment. The sequence alignment process entailed mapping each sequence with either the hidden Markov model (HMM) or DIAMOND to match it with the eggNOG database. The ideal matching sequence of the target sequence is classified according to its taxonomy and further categorized and annotated utilizing gene ontology (GO) (Gene Ontology Consortium, 2015) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways (Kanehisa et al., 2017) through the application of the clusterProfiler v.3.19 package in R (Xu et al., 2024). Furthermore, BLAST analysis was performed to examine the 80 genomes for the presence of genes linked to virulence and pathogenicity. We employed BLASTN via ABRicate v1.0.1³ (de Man and Limbago, 2016) and the Virulence Factor Database (VFDB) (Liu et al., 2022) to examine the occurrence of virulence factors and antibiotic resistance gene in all strains of M. marinum. The analysis was performed using a 90% identity threshold and a 30% coverage threshold.

To investigate the RGP in eight complete genomes of M. marinum, PPanGGOLiN software⁴ (Gautreau et al., 2020) was utilized with “panrgp” command. Genomic plasticity regions, which include genomic islands, prophages, and other variable elements, were identified using PPanGGOLiN’s comprehensive pangenome analysis capabilities. Initially, the genomes were prepared and formatted according to requirements based on the sample input dataset of the software. PPanGGOLiN was then employed to perform a detailed analysis, identifying regions of genomic variability and plasticity across the M. marinum strains. The tool provided insights into the distribution and composition of these plastic regions, allowing to map and characterize their presence and variability among the genomes. This analysis revealed the dynamic regions contributing to the genomic diversity and adaptability of M. marinum.

Horizontal gene transfer is a vital factor in the evolution and diversification of numerous microbial species. The ensuing dynamics of gene acquisition and loss can significantly impact the emergence of antibiotic resistance and influence the development of vaccines and therapeutic interventions (Tonkin-Hill et al., 2023). Panstripe⁵ improves the post-processing of bacterial pangenome analyses (Tonkin-Hill et al., 2023). Panstripe accepts a phylogeny and a gene presence/absence matrix in Rtab format, which can be generated by tools like Panaroo and Roary. Utilizing a Generalized Linear Model (GLM) framework allows for straightforward comparison of slope terms across different datasets. The “compare_pangenomes” function analyzes the interaction term between the datasets alongside the core, tip, and depth terms. A significant p-value for the tip term indicates that the two pangenomes differ in the rates of gene presence and absence observed at the tips of the phylogeny. This variation is often attributed to differences in annotation error rates between datasets or to fluctuations in the gain and loss of highly mobile genetic elements that may not persist long enough to be detected across multiple genomes. A significant p-value for the core term suggests differing rates of gene gain and loss between the two datasets. Although the depth term is less critical, it reflects discrepancies in the ability to detect older gene exchange events between the two pangenomes. The dispersion parameter reveals whether there is a notable difference in the dispersion of the two pangenomes, implying that the relationship between gene exchange rates and the size of each event varies between the two datasets. The p-value is derived using a Likelihood Ratio Test, and a subset of inferred GLM parameters can be analyzed. The resulting p-values and bootstrap confidence intervals can help assess the significance of each term in the model to gene gain and loss (Tonkin-Hill et al., 2023).

In this benchmarking analysis, we compared the annotation performance of EggNOG-mapper and InterProScan for M. marinum, focusing on metrics such as GO term accuracy, average GO terms per protein, and proteome coverage (Huerta-Cepas et al., 2017). Annotation files generated by each tool were parsed using the pandas library in Python v3.8, extracting relevant data such as GO terms and Pfam domains. The data were structured into a standardized format, capturing sets of GO terms, functional categories, and proteome coverage for both tools. To evaluate the performance of each tool, key metrics were calculated. These included the proportions of true positive (TP) and false positive (FP) annotations, which were used to assess the accuracy of GO term predictions. Additionally, the average number of GO terms per protein was computed to reflect the specificity of the annotations. Finally, proteome coverage was evaluated based on the proportion of proteins that were annotated with only TP terms, a mix of TP and FP terms, or no TP terms. This analysis provides a comparative framework for assessing annotation quality, highlighting the strengths and limitations of each tool’s functional annotation capabilities for genomic studies.

The genomic features are represented in Supplementary Table 1, which lists the assembly accessions, organism infraspecific names, assembly release dates, number of scaffolds, and isolated species for eight complete genomes of M. marinum. Genetic variability between the fish isolates, such as E11, CCUG20998, and 1218R, and the human isolates, including H01, MMA1, and Mycobacterium, is minimal according to the ANI analysis (Figure 1). The heatmap shows that the ANI percentages are consistently high, nearing or above 98%, indicating close genetic relatedness between strains from fish and human sources. This suggests that despite the different host origins, these M. marinum strains share significant genomic similarity.

Pangenome analysis of M. marinum based on the eight complete genome sequences reveals a large proportion of core genes and significant genomic diversity. The gene presence-absence matrix (Supplementary Figure 1) showed a clear clustering of strains, with the core genome represented by 4634 gene clusters (Supplementary Figure 2). This suggests that a substantial fraction of genes is conserved across all strains, indicating evolutionary stability in essential functions. Additionally, the shell genome comprises 1390 gene clusters, reflecting moderate variability among strains. No soft-core genes were detected, while 1758 gene clusters belong to the cloud genome, representing strain-specific genes. These cloud genes highlight the presence of unique, possibly adaptive, traits linked to the strain’s environmental or host-specific pressures.

Figure 2 displays a circular representation of the core genome and pan-genome of M. marinum. The outermost circle represents the core genome, which is shared by all strains. The inner circles represent the pan-genome, which consists of genes that are present in at least one strain. The pangenome analysis reveals that the core genome of M. marinum is relatively small, suggesting that the species has a high degree of genetic diversity. The presence of many singleton genes in each strain suggests that M. marinum is rapidly evolving and adapting to different environments.

Figure 3 illustrates the open nature of the pangenome using different fitting models. A power-law fit to depict the cumulative number of genes discovered as more genomes are sampled (Figure 3A). The steady slope in the fitted curve indicates continuous gene discovery, confirming that the pangenome remains open, as suggested by the equation y = 0.0917x + 0.0012y. Heap’s law (Figure 3B) shows the number of unique genes discovered with increasing genome samples and further reinforces the open pangenome model, as new unique genes continue to emerge with additional genome sampling.

The Mauve software was utilized to create a visual representation of gene synteny across the genomic sequences of various M. marinum strains, allowing for the identification of potential genetic rearrangements based on sequence similarity (Supplementary Figure 3). The presence of gene blocks facilitated the detection of synteny, enabling the accurate identification of these regions. A comparative analysis of the genomes revealed several fragmented areas, characterized by multiple inversions and deletions. Notably, despite these disruptions, a high degree of uniformity was observed, with large, conserved segments exhibiting inversions and deletions when compared across the genomes.

The Gene Ontology (GO) enrichment (Figure 4) and KEGG pathway analysis (Figure 5) of the M. marinum pangenome annotation, conducted through eggNOG and InterProScan, provide a comprehensive overview of the functional capabilities of the genes identified. The GO enrichment analysis reveals that biological processes like positive regulation of RNA biosynthetic processes, DNA-templated transcription, carbohydrate derivative biosynthesis, mitochondrial function, and cell adhesion are highly significant. These processes underscore key roles in cellular organization, metabolic activity, biosynthesis, and cell signaling. The negative regulation of nucleobase-containing compound metabolism and RNA metabolic processes also suggests intricate regulatory mechanisms at play in gene expression and metabolic control.

The KEGG pathway analysis highlights critical pathways including the PI3K-Akt and MAPK signaling pathways, both essential for cell survival, growth, and immune responses. Other enriched pathways like neuroactive ligand-receptor interaction, cytokine-cytokine receptor interaction, and calcium signaling pathways point to significant communication and signaling mechanisms that might play a role in host-pathogen interactions.

The virulence factor analysis of all the complete genomes identified several key virulence genes, highlighting important duplications: esxB, fbpC, fbpB, esxG, esxH, eccA3, eccC3, esxM, esxA, ideR, relA, hbhA, phoP, and mbtH (Table 1). Notably, the duplication of fbpB and fbpC, both associated with mycolyl transferase activity, suggests functional redundancy in adherence mechanisms. Additionally, the presence of multiple components of the Type VII secretion system, including esxB, esxG, esxH, eccA3, and eccC3, underscores the significance of this system in the virulence of M. marinum.

The RGP analysis of the pangenome of M. marinum revealed significant genomic diversity among the different strains. A total of 34 RGPs were identified in the study, with some RGPs present in multiple strains, indicating conservation among different strains. Notably, antimicrobial resistance genes, such as those conferring resistance to dibekacin, netilmicin, tobramycin, and Fosfomycin, were widespread among the strains. Additionally, antibiotic-resistant murA transferase, aminosalicylate resistant, and virulence factor genes were found to be co-localized in the same RGP, suggesting potential co-regulation. The analysis also identified several modules (or genes) that were present in multiple strains, including modules 26, 32, 28, 30, and 44. Strain-specific RGPs were also observed, with strain CCUG20998 having a unique RGP (RGP_0) and strain MMA1 having a unique RGP (RGP_14). Overall, the RGP analysis highlights the complex genomic landscape of M. marinum and the potential for co-regulation of antimicrobial resistance and virulence genes.

The analysis of the open pangenome model reveals critical insights into the temporal dynamics of gene gain and loss events. The significant negative estimate of the “core” term (–0.705, p = 1.02e-4), coupled with the narrow bootstrap confidence interval (–1.15, –0.255), indicates strong evidence for an open pangenome (Figure 6). This suggests that gene gain and loss events are continuously accumulating over time, characteristic of an open system where new genes are incorporated as the genome evolves. The accumulation is not temporally limited, signifying ongoing horizontal gene transfer or gene acquisition from various sources, possibly through mobile genetic elements, which are typically associated with open pangenomes. Furthermore, the non-significant estimates for “depth” (estimate = –0.0195, p = 4.62e-1) and “tip” (estimate = –0.153, p = 4.31e-1) align with the expected pattern for an open pangenome. These results suggest that there are no substantial differences in detecting older gene exchange events or gene presence at the tips, reinforcing the idea that gene gain and loss occur throughout the evolutionary history of the genomes, not confined to specific time points. The non-significance of the dispersion term also implies that the rate and size of gene exchange events remain consistent, rather than varying dramatically between different branches of the phylogeny. The plot comparing cumulative gene gain and loss events against core branch distance shows distinct patterns for fast and slow-evolving pangenomes. The rapid accumulation of genes in the fast-evolving pangenome supports the interpretation of an open system, where more gene exchange events occur over shorter evolutionary timescales, especially along terminal branches, as indicated by the denser clustering (Figure 6). This continuous acquisition of genes aligns with the expected behavior of an open pangenome, where there is no saturation point for gene gain, and new genes continue to enter the genome pool over time.

The benchmarking analysis presented in Figure 7 provides a comparative evaluation of EggNOG-mapper and InterProScan for annotating M. marinum genomes. Figure 7A showing a higher proportion of true positive annotations for both tools, with EggNOG-mapper achieving 80% true positives compared to 70% for InterProScan. False positive rates are comparatively lower, indicating that both tools provide a reliable level of annotation accuracy, though EggNOG-mapper has an advantage. The average GO terms per protein, reveals that EggNOG-mapper assigns a higher average of curated (true positive) terms per protein, reflecting greater annotation specificity and comprehensiveness (Figure 7B). InterProScan, while also providing curated terms, shows a slightly lower average, particularly in false positive and non-curated terms, which suggests a more conservative annotation approach. Figure 7C demonstrating that both techniques offer substantial coverage exclusively with true positive (TP) annotations. EggNOG-mapper covers 60% of the proteome with TP terms only, whereas InterProScan covers 55%, indicating robust proteome-level annotation by both tools, though EggNOG-mapper marginally outperforms in exclusive TP coverage. Together, these results indicate that EggNOG-mapper provides a more extensive annotation profile across metrics, making it slightly more effective in annotation accuracy, specificity, and proteome coverage.