Whole-genome and pangenome insights into Mycobacterium colombiense clinical isolates from human infections

Bacterial isolates

M. colombiense clinical isolates were obtained from the collection at the Clinical Microbiology Laboratory at Srinagarind Hospital, Khon Kaen University, Khon Kaen, Thailand. These isolates were collected from patients diagnosed with MAC infections caused by NTM. Diagnosis was made in accordance with the guidelines of the American Thoracic Society and the Infectious Diseases Society of America (Griffith et al., 2007). All MAC isolates were initially identified using the INNO-LiPA MYCOBACTERIA v2 line-probe assay (Fujirebio, Ghent, Belgium) and subsequently sub-cultured on Löwenstein–Jensen medium for further analysis.

The study protocol was approved by the Institutional Review Board of the Khon Kaen University Ethics Committee in Human Research (no. HE671519; approval date: September 24, 2024). All experiments were conducted in accordance with the Declaration of Helsinki and ICH Good Clinical Practice guidelines. Informed consent was waived for the use of medical data as patient information was anonymized and de-identified prior to analysis.

Drug susceptibility testing

Minimum inhibitory concentrations (MICs) were determined using the Sensititre™ Myco SLOMYCOI assay (TREK Diagnostic Systems, West Sussex, UK) following the manufacturer’s protocol. The assay included serial two-fold dilutions of 13 lyophilized antimicrobial agents: clarithromycin (0.06–64 μg/mL), amikacin (1–64 μg/mL), rifabutin (0.25–8 μg/mL), rifampicin (0.12–8 μg/mL), ethambutol (0.5–16 μg/mL), trimethoprim/sulfamethoxazole (0.12/2.38–8/152 μg/mL), ciprofloxacin (0.12–16 μg/mL), moxifloxacin (0.12–8 μg/mL), ethionamide (0.3–20 μg/mL), isoniazid (0.25–8 μg/mL), doxycycline (0.12–16 μg/mL), linezolid (1–64 μg/mL), and streptomycin (0.5–64 μg/mL). MIC interpretations followed the Clinical and Laboratory Standards Institute (CLSI) guidelines (CLSI, 2018).

Genomic DNA extraction and whole-genome sequencing

Genomic DNA was extracted using the cetyl-trimethyl-ammonium bromide-sodium chloride method (Larsen et al., 2007), followed by RNase treatment and ethanol precipitation for purification. DNA concentration and quality were assessed using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Carlsbad, CA, USA). High-quality DNA was sent to NovogeneAIT (Hong Kong) for whole-genome sequencing. Paired-end libraries were prepared and sequenced on the Illumina HiSeq 2500 platform with 150-bp read lengths. Sequence data are deposited in the Sequence Read Archive under BioProject accession no. PRJNA1271787.

Species identification and genome assembly

Raw reads were quality checked using FastQC version 0.11.9 (Andrews, 2010). Adapter sequences and low-quality reads were removed using Trimmomatic version 0.39 (Bolger, Lohse & Usadel, 2014). High-quality paired-end reads were de novo assembled using Unicycler version 0.4.8 (Wick et al., 2017). Assemblies were scaffolded using RagTag version 2.1.0 (Alonge et al., 2022) and polished with Pilon version 1.24 (Walker et al., 2014) to improve base accuracy. Assembly quality was evaluated using QUAST version 5.3.0 (Gurevich et al., 2013). Genome completeness and contamination were evaluated using CheckM version 1.2.3 (Parks et al., 2015), applying thresholds of >97% completeness and <3% contamination. Final species identification and genome validation were performed using NTM-Profiler version 0.2.0 (Phelan, 2021), the Type Strains Genome Sever (TYGS) (Meier-Kolthoff et al., 2022), Pathogenwatch (https://pathogen.watch/), and GTDB-Tk (Chaumeil et al., 2022). Pairwise average nucleotide identity (ANI) was calculated using FastANI (Jain et al., 2018) for accurate taxonomic classification. ANI-based clustering and heatmap visualization were generated using ANIclustermap version 1.4.0 (Shimoyama, 2022).

Genome characterization

Gene prediction and functional annotation were performed using the RAST toolkit (RASTtk) through the Bacterial and Viral Bioinformatics Resource Center (BV-BRC) with a comprehensive genome analysis pipeline (Olson et al., 2023). Antibiotic resistance genes (ARGs) and virulence factors were identified using ABRicate version 0.7 (Seemann, 2020) with default parameters (a minimum sequence identity threshold of 80% and a coverage of 80%). The Comprehensive Antibiotic Resistance Database (CARD) was used to detect resistance genes and the Virulence Factor Database (VFDB) was used to identify virulence-associated genes. ARGs were identified using Resistance Gene Identifier (RGI) version 6.0.5 restricted to “Perfect and Strict hits only” for high-confidence predictions against the most recent version of the CARD (Alcock et al., 2023). The RGI predicts resistomes by utilizing curated homology and single nucleotide polymorphism (SNP) models to detect both acquired resistance genes and chromosomal resistance mechanisms. CRISPR loci and CRISPR-associated (Cas) proteins were detected using CRISPRCasFinder (Couvin et al., 2018) with high-sensitivity parameters to identify both confirmed and putative CRISPR arrays, as well as Cas subtypes.

Detection of genetic mutations associated with drug resistance

Genetic mutations associated with drug resistance were analyzed in key target genes, including 16S rRNA (rrs), 23S rRNA (rrl), gyrA, gyrB, large ribosomal protein L3 (rplC), and large ribosomal protein L4 (rplD). The rrs gene is the resistance target for amikacin, whereas rrl is the target for clarithromycin. The gyrA and gyrB genes serve as targets for moxifloxacin, and rplC and rplD correspond to linezolid. Gene sequences were extracted from the assembled M. colombiense genomes using BLASTn followed by BEDtools getfasta, and subsequently aligned with the reference genome M. colombiense CECT 3050 (accession no. NZ_AFVW02000005.1) using the MUSCLE algorithm in MEGA (Kumar et al., 2024) to identify mutations associated with resistance.

SNP-based phylogenetic analysis

SNP-based phylogenetic analysis was conducted to infer evolutionary relationships among M. colombiense isolates. SNP calling was performed using Snippy version 4.6.0 (Seemann, 2015) which aligned reads to the M. colombiense reference genome CECT 3035 (NZ_CP020821.1). Core genome SNPs were extracted using Snippy-core, generating a core alignment of 5,581,643 bp, of which 147,025 positions were identified as polymorphic. Recombination was filtered using Gubbins version 3.1.6 (Croucher et al., 2015), resulting a final non-recombinant alignment of 111,537 SNP positions used for downstream phylogenetic reconstruction. All tools were run with default parameters unless otherwise specified. Maximum likelihood trees were constructed using IQ-TREE version 2.0.7 (Minh et al., 2020) with 1,000 bootstrap replicates and the optimal substitution model selected using ModelFinder Plus. The phylogenetic tree was visualized using iTOL version 7 (Letunic & Bork, 2021). Pairwise SNP distance matrices were analyzed using snp-dists version 0.8.2 (Seemann, Klötzl & Page, 2021).

Public genome retrieval and species confirmation

The 18 publicly available genome assemblies were retrieved from the National Center for Biotechnology Information (https://www.ncbi.nlm.nih.gov/datasets/genome/) (accessed February 2025). Detailed information is provided in Table 1. All genomes were quality checked using CheckM, applying thresholds of >97% completeness and <3% contamination. Species confirmation was performed using NTM-Profiler, TYGS, Pathogenwatch, GTDB-Tk, and FastANI. The publicly available genomes identified as M. colombiense were combined with the 12 clinical isolates from this study for pangenome analysis.

Pangenome analysis

All M. colombiense genomes were re-annotated using Prokka version 1.14.6 (Seemann, 2014) to generate GFF3 files. Pangenome analysis was conducted using Roary version 3.12.0 (Page et al., 2015) with default parameters, including a 95% identity threshold for BLASTp and a core gene definition of ≥99% presence across genomes. Partial coding sequences were excluded and translated protein sequences were used for downstream analysis. Core genome alignments were refined using MAFFT (Katoh & Standley, 2013) for improved alignment accuracy. The gene presence–absence matrix and summary statistics were generated. Genes were categorized as core (present in all genomes), soft-core (present in most but not all genomes), accessory (present in several genomes), or unique (strain-specific). Visualization of the pangenome analysis was performed using the Python script [roary_plots.py] and the R program with the ggplot2 package. Functional annotation of core, accessory, and unique genes was performed using the Clusters of Orthologous Groups (COGs) system via the Pan-Explorer platform (Dereeper, Summo & Meyer, 2022), allowing for comparative analysis of gene functions related to metabolic pathways and functional roles.

Statistical analysis

Statistical analysis was performed using IBM SPSS Statistics for Windows, version 28.0 (IBM Corp, Armonk, NY, USA). Descriptive statistics were used to summarize the study population characteristics. Normality of the data distribution was evaluated using the Shapiro–Wilk test. Continuous variables with a normal distribution are presented as the mean ± standard deviation, followed by the range (minimum–maximum), while non-normally distributed data are presented as the median. Categorical variables are reported as frequencies and percentages.

Characteristics of the studied isolates

A total of 12 M. colombiense clinical isolates were identified using bioinformatics tools, including NTM-Profiler, TYGS, Pathogenwatch, GTDB-Tk, and FastANI (Table S1). Most isolates (75%) were obtained from male patients with a mean age of 44.25 ± 15.33 years (range: 27–79 years). Isolates originated from various clinical specimens, including sputum, tracheal suction, skin, pus, synovial fluid, and others (Tables 1, S2).

Drug susceptibility testing was performed on 10 of the 12 isolates using the Sensititre™ Myco SLOMYCOI assay (Table 2), with two isolates excluded due to inadequate growth. MIC breakpoints for amikacin, clarithromycin, linezolid, and moxifloxacin were interpreted according to CLSI guidelines. Results revealed universal resistance to moxifloxacin (100%) and high resistance to linezolid (90%). Moderate resistance was observed for clarithromycin and amikacin (30% each). The MIC50 and MIC90 values demonstrated that resistance to moxifloxacin was consistent at ≥8 µg/mL. Linezolid had elevated MICs at 32 µg/mL (MIC50) and ≥64 µg/mL (MIC90). Susceptibility to clarithromycin was variable (MIC50 of 4 µg/mL and MIC90 of 64 µg/mL); amikacin demonstrated intermediate efficacy (MIC50 of 16 µg/mL and MIC90 of 64 µg/mL).

Gene prediction and functional annotation

The average genome size of the 12 isolates was 5.79 Mbp (Table 3). The number of predicted coding sequences ranged from 5,350 to 5,665, with an average of 5,471. All isolates had a GC content of approximately 67.84%. Each genome contained 3 rRNA genes, while tRNA gene numbers ranged from 49 to 51.

Genome annotation using the RASTtk on the BV-BRC platform identified consistent numbers of annotated antibiotic resistance and virulence genes across isolates. The Victors database reported four resistance genes per genome, the CARD database detected 11–13 genes, and the PathoSystems Resource Integration Center (PATRIC) database identified the highest number (43–46 genes). Similarly, virulence gene counts varied by database: Victors reported 98–108, VFDB identified 13–22, and PATRIC detected 149–172.

Functional categorization revealed that most genes were involved in metabolism, followed by energy production, stress response, defense mechanisms, virulence, and protein processing. These distributions highlight the organism’s adaptability and survival strategies (Fig. 1).

Antibiotic resistance, virulence factors, and CRISPR-Cas profiling

Analysis using the PATRIC and BV-BRC pipelines with the CARD database via ABRicate revealed the presence of multiple resistance-associated genes, including those related to antibiotic target modification (alr, ddl, dxr, EF-G, EF-Tu, embA, embB, embC), activation (katG), protection (mfpA), and target replacement (fabG, htdX). Efflux pump genes (mmpL5, Rv1634, mtrA, efpA), cell wall modification genes (gdpD, pgsA), and various regulatory genes (embR, ethR, lpqB, marR, mtrA, mtrB, oxyR) were broadly distributed. The Rv2994 gene was absent in some isolates (Fig. 2, Table S3). Targeted interrogation of canonical resistance loci using RGI/CARD detected no resistance-associated mutations.

Virulence gene analysis using VFDB showed conserved presence of genes associated with toxin production, ESX secretion systems (esxH, esxM, esxN, ecc), adhesion (hbhA, fbp), immune evasion (phoP, icl), and nutrient acquisition (mbtN, absent in two isolates). Stress- and persistence-related genes (relA, ideR) were universally present (Fig. 3, Table S3).

CRISPR-Cas system analysis revealed variability in the number of CRISPR loci among isolates, ranging from six to nine per genome.

Mutations in genes associated with phenotypic drug resistance

For amikacin, resistance is typically associated with mutations in rrs. No rrs mutations were detected in any M. colombiense isolates. For clarithromycin, mutations in rrl confer macrolide resistance. In this study, rrl mutations at positions A2270G/C (A2059G/C based on E. coli numbering) were identified exclusively in clarithromycin-resistant isolates.

For linezolid, resistance-associated mutations were assessed in rplC and rplD. Because all included isolates were linezolid-resistant, comparisons between susceptible and resistant groups could not be performed. Six synonymous mutations were identified in rplC, and eight mutations were detected in rplD, including one nonsynonymous substitution (Val218Ala).

For moxifloxacin, mutations within the quinolone resistance–determining regions (QRDRs) of gyrA and gyrB were analyzed. Several mutations were identified; however, the absence of moxifloxacin-susceptible isolates prevented definitive correlation between genotype and phenotype. Notably, three GyrA substitutions (Ala218Val, Val453Ile, Gly833Ser) and two GyrB substitutions (His250Arg, Thr522Ala) were observed (Tables 4, S4 and S5).

SNP-based phylogenetic tree

Phylogenetic reconstruction based on SNP data revealed three main clades among the 12 M. colombiense clinical isolates (Fig. 4). Clade divisions were made according to pairwise SNP distances (Table S6). Clade I comprised one isolate (MCO121383869) which exhibited median SNP distances of 43,009 and 52,145 from Clades II and III, respectively. Clade II comprised two isolates (MCO102081859 and MCO062281265) which had a median of 7,598 SNPs from each other and median distances of 43,009 and 62,434 SNPs from Clades I and III, respectively. Clade III comprised nine isolates with median SNP distances of 52,145 and 62,434 from Clades I and II, respectively. Within Clade III, two subgroups were evident: Subgroup 1 included MCO021680268, MCO101582585, MCO061681546, and MCO051681265, with median pairwise distances of 47 SNPs (range: 31–57). Subgroup 2 included MCO031580717, MCO061481978, MCO061481970, MCO051681160, and MCO051681085, with median pairwise distances of 8 SNPs (range: 3–16).

Public genome retrieval and species confirmation

Species confirmation was performed on 30 M. colombiense genomes, comprising the 12 clinical isolates from the current study and the 18 publicly available genomes (Table S2). The genome analysis results from NTM-Profiler, Pathogenwatch, and GTDB-Tk were classified as M. colombiense. However, TYGS analysis using digital DNA–DNA hybridization (dDDH) values (species threshold ≥70%) indicated that only 16 genomes (12 from this study and 4 from public databases) exceeded 80% similarity (dDDH-d₄), consistent with M. colombiense classification. FastANI analysis (species threshold ≥95%) revealed that the same 16 genomes had ANI values above 97% (Fig. S1). To ensure accurate representation, only genomes meeting both criteria (dDDH ≥70% and ANI ≥95%) were included in the pangenome analysis.

Pangenome analysis

Pangenome analysis of the 16 confirmed M. colombiense genomes identified 7,771 gene clusters comprising 4,468 core genes (57.50%), 1,834 accessory genes (23.60%), and 1,469 unique or cloud genes (18.90%) (Fig. 5A). The high proportion of core genes confirmed a closed pangenome structure, indicating evolutionary conservation and low genetic variability among isolates. Gene cluster distribution (Fig. 5B) followed a typical pangenome pattern, with many genes shared across all genomes (core) and a smaller subset found only in a few genomes (accessory or strain-specific). The plateauing of the conserved gene curve (Fig. 5C) supports the closed pangenome structure, reflecting limited genetic acquisition. Conversely, the rapid decline in new gene discovery (Fig. 5D) indicates genomic stability across isolates, consistent with a conserved species. The gene presence–absence matrix and heatmap (Fig. 5E) suggests extensive gene content similarity among isolates, dominated by core genes. Only a small proportion of clusters varied between isolates, reflecting the low genomic diversity within the species.

Functional classification using the COGs system revealed gene assignments across 22 categories (Fig. 6). Categories R (general function prediction only) and S (function unknown) were the most abundant across core, accessory, and unique gene groups, indicating a large proportion of uncharacterized functional genes. Among classified functions, Category E (amino acid transport and metabolism) was most prevalent in the core genome (645 genes), followed by Category I (lipid transport and metabolism, 643 genes) and Category Q (secondary metabolite biosynthesis, 605 genes). Accessory and unique genes were relatively enriched in Category Q, although their absolute numbers were significantly lower than those in the core genome. These findings suggest that while unique genes may contribute to niche-specific adaptations, most essential metabolic and functional pathways are conserved within the core genome.