Phylogenetic Analysis of Clinical Non-Tuberculous Mycobacteria Using 16S rRNA, hsp65 and rpoB Gene Markers: Implications for Species Identification and Clinical Management

1. Introduction

Nontuberculous mycobacteria species (NTM) are a group of environmental organisms that can be found in water, including drinking water, soil, animals, plants, food, vegetation, and human feces. These species have the potential to cause diverse pathologies, ranging from pulmonary, skin and soft tissue infections to deeper infections with or without systemic dissemination. Some NTM species are well-recognized pathogens, while others are only recently emerging, and their pathogenic potential remains unknown [1,2].

Currently, diseases caused by NTM have increased compared to M. tuberculosis. This may be related to the introduction of more sensitive laboratory techniques, the development of genomic techniques that allow faster diagnosis of mycobacterial diseases and a rapid differentiation of NTM [3]. More than 150 Mycobacterium species have been identified, with a fast pace increase in species identified due to the advancement of genetic techniques [4]. The progress of genetics and the development of mathematical models for analysis, in addition to the contribution of the global network of genetic sequences available in real time in public domain databases, have contributed to the advancement and more common use of phylogenetic analyses [5].

Phylogenetic and taxonomic studies of Mycobacterium species have been largely based on the analysis of the 16S RNA gene. Its high information content, conserved nature and presence in all organisms have made this gene a suitable candidate for phylogenetic analysis and identification. However, some limitations persist, mainly the genetic similarity (greater than 94%) that exists between species of the genus Mycobacterium [6]. For this reason, the combination of molecular marker analysis, such as rpoB and hsp65, together with the analysis of the 16S rRNA gene, has become an effective system for identifying mycobacterial species [7]. The 65 kDa heat shock protein gene (hsp65), highly conserved among mycobacterial species with hypervariable regions [8], and the rpoB gene, which encodes the β-subunit of RNA polymerase and contains regions of conserved sequences flanking highly variable regions [9], are used as alternative or complementary tools to identify mycobacteria [10].

Correct identification of NTM species is important for selecting the appropriate treatment for patients or epidemiological purposes [11]. Treatment depends on the specific NTM species and may range from the use of one or more antimicrobials, anti-tuberculous agents, to surgical debridement. Therefore, most NTM isolates should be identified at the species level and not only as species groups [12], both for treatment and to establish whether it is a colonizer of the clinical specimen or a true pathogen.

NTM species are traditionally classified by their growth rates into rapidly growing mycobacteria (RGM) and slowly growing mycobacteria (SGM) based on the time required for visible colony formation. Species that produce visible colonies on subculture plates within 3–7 days are classified as RGM, whereas those that take longer than 7 days are classified as SGM [13].

Because NTMs are a heterogeneous group of microorganisms, drug susceptibility greatly varies within species [14]. Currently, treatment of almost all RGM infections includes antibiotics such as macrolides (clarithromycin, azithromycin), sulfas (trimethoprim-sulfamethoxazole), aminoglycosides (amikacin, tobramycin), quinolones (moxifloxacin, ciprofloxacin), oxazolidinones (linezolid), cephamycins (cefoxitin), carbapenemics (imipenem), and tetracyclines (minocycline, doxycycline and tigecycline). While for infections caused by the SGM group, the regimen additionally includes ethambutol, isoniazid, streptomycin, and rifampicin (anti-tuberculous) [15].

Differences in antimicrobial susceptibility that determine treatment options among different mycobacteria mean that identification of NTMs to the species level is becoming increasingly important clinically. Mycobacterial susceptibility testing, therefore, is important for appropriate patient management and should be performed on clinically significant isolates of certain NTMs [12].

Although intraspecific and interspecific variations among mycobacteria from different geographic regions have been extensively studied using DNA sequences of 16S rRNA, rpoB and hsp65 genes, for ecological genetic studies, phylogenetic and evolutionary analyses, there is a lack of information on genetic variation among mycobacteria isolated from populations from various geographic regions, where mycobacterial infection is a major health problem [16]. Similarly, Ecuador has limited information on the genetic and epidemiological diversity of mycobacteria isolated from human samples [17] and in greater proportion in environmental NTM [18].

Therefore, in this study, we propose to develop a phylogenetic analysis using 16S rRNA, hsp65 and rpoB gene sequences for taxonomic categorization and accurate species identification, genotype-phenotype correlation and understanding of the prevalence and phylogenetic structure of NTM isolated from clinical samples from a private laboratory in the city of Quito, Ecuador.

2. Materials and Methods

2.1. Samples/Sequences Selected for the Study

All 2018-2022 samples (57) of the Mycobacterium genus stored in the Microbiology service of Zurita & Zurita Laboratories, with respiratory and non-respiratory origin, were selected for PCR amplification and sequencing. Previously, NTMs were identified using acid-fast staining (Ziehl-Neelsen). Samples were seeded in Loewenstein-Jensen medium and Middlebrook 7H9 on the BACTEC 320 automated bacteriological system. For the purpose of this study, the previously isolated and identified strains were thawed and subcultured on chocolate agar medium for fast-growing strains and on Lowenstein-Jensen medium for slow-growing strains and incubated at 35°C with 5% CO2.

Definitive identification was performed by molecular methods using 16S rRNA, rpoB, and hsp65 gene sequencing. The sequences were entered into GenBank for preliminary identification; however, in certain cases where the species are closely related, the sequences could not be definitively identified. Therefore, all species were subjected to phylogenetic analysis.

Fisher's test was used to evaluate differences in the data for the different mycobacterial species analyzed. Statistical significance was set at p < 0.05. Analyses were performed in Past4.06b [19].

2.2. Antimicrobial Resistance Analysis

Positive samples were detected by the BACTEC 960 system (Becton Dickinson, USA). Antibiotic resistance testing was evaluated using the broth microdilution technique (MIC) with Sensititre system plates (Thermo Fisher Scientific, USA) for slow-growing NTM (SLOMYCO) and fast-growing NTM (RAPIDMYCO), according to the manufacturer's recommendations [11].

The antibiotics tested for RGM were trimethoprim/sulfamethoxazole, amikacin, clarithromycin, moxifloxacin, ciprofloxacin, linezolid, doxycycline, tigecycline, tobramycin, imipenem, minocycline, amoxicillin/clavulanic acid, cefepime, ceftriaxone, and for SGM were: trimethoprim/sulfamethoxazole, amikacin, clarithromycin, moxifloxacin, ciprofloxacin, linezolid, doxycycline, ethambutol, streptomycin, ethionamide, isoniazid, rifabutin, and rifampicin. The interpretation was performed according to CLSI guidelines [20]. Many mycobacteria do not have established breakpoints, therefore, in addition to the CLSI, other guidelines were used [21]. The antibacterial resistance index (ARI) [22] was calculated to assess resistance in the different mycobacterial species identified. This index relates pathogen susceptibility and antibiotic use to represent the overall level of resistance and efficacy of antibiotic therapy [23].

2.3. Genomic Extraction, Amplification and Sequencing

Genomic DNA was extracted using the QIAamp DNA Minikit kit (QIAgen, Hilden, Germany) directly from chocolate agar medium or Lowenstein-Jensen's medium once NTMs formed visible colonies. Procedures and reagents were used according to the manufacturer's protocol. Eluted DNA was stored at -20°C until PCR analysis. Amplification of 16S rRNA, rpoB and hsp65 genes was performed. The primers for amplification and the expected amplicon sizes are listed in Table 1.

The following cycling conditions were used: initial denaturation of 4 min at 95°C; 35 cycles of 30s at 95°C, 1 min at 60°C (hsp65), 55°C (rpoB), 57°C (16S rRNA), 1 min at 72°C; and a final extension of 10 min at 72°C. PCR products were visualized on 1.5% agarose gel electrophoresis under ultraviolet light.

The PCR products of the three genes were sequenced by the Sanger method using the same amplification primers, with the exception of the 16S rRNA gene, where primers 518F (5'-CCAGCAGCAGCCGCGGGGTAATACG-3') and 800R (5'-TACCAGGGTATCTAATCC-3') were used to cover the entire region [26].

2.4. Sequence Alignment and BLAST

The sequences produced by the forward and reverse primers were aligned using Sequencher 5.4.6 (Gene Codes Corporation, Ann Arbor, MI, USA) to obtain a consensus sequence [27], labeled with a respective code. Consensus sequences were aligned against the NCBI GenBank database using the BLASTn local alignment tool to obtain preliminary identity assessments, including percent nucleotide similarity and coverage values. Mycobacterial sequences for each marker (16S rRNA, rpoB, and hsp65) were retrieved from NCBI. A total of 167 sequences were included: 53 for 16S rRNA, 57 for rpoB, and 57 for hsp65. Sequences from related reference groups and outgroups, such as Tsukamurella species, Gordonia, and Nocardia [28], were also incorporated in the analysis. The Mycobacterium species and reference-group sequences used are listed in Supplementary Table 1, along with their respective localities, sample origins, and GenBank accession numbers.

An alignment matrix was constructed for each gene or marker using ClustalW in MacVector software 18.5.8 with Gap Creation Penalty-GOP and Gap Extension Penalty-GEP parameters of 30.0 and 10.0, respectively, looking for the highest degree of positional homology.

2.5. Phylogenetic Analysis

The respective aligned matrices were used for the phylogenetic hypothesis reconstruction as a second level of identification, using PAUP 4.0a (build 169) software for Maximum Parsimony (MP) optimization [29] and MEGAX for Maximum Likelihood (ML) were conducted using MEGA version X [30,31] of the 16S rRNA, rpoB and hsp65 genes.

The nucleotide substitution model of each matrix for ML analysis was calculated by Model Test under MEGAX, calculated by Akaike's criterion, and Bayesian Information Criterion (BIC). The ML tree was constructed with an initial MP tree and the respective model for each matrix-gene-marker.

MP construction was performed with 1,000 replicates of random addition of taxa and characters and branch swapping with TBR (Tree Bisection and Reconnection). Subsequently, homologous character reweighting was performed using the recalculated consistency index (RCI) and values. The branch support (in percentages) was obtained using 1,000 bootstrap pseudoreplicates of the matrix as a posteriori statistical support of clades or branches [32]. The sequence of analyses is based on previous work [28,33,34]. Using the same methodological approach, after performing individual analyses for each gene, a concatenated analysis of all three genes was conducted to confirm the final topology and assess the contribution of the information of each gene to the combined phylogeny.

As a third level of verification of molecular identification and for purposes of classification into complexes, groups, and verification of phylogenetic species (phylogenetic species concept), each monophyletic clade was evaluated by inter- and intra-clade divergence through the construction of distance matrices using the uncorrected p-distance model, as well as between sequences obtained from NCBI and those obtained in the laboratory.

2.6. Ethical Considerations

The mycobacteria samples come from the Zurita & Zurita Laboratories Biomedical Research Unit strain, which approved the use of mycobacteria from the strain, and has no information related to patient data.

3. Results

Sensitivity profile data, the nature of the clinical specimen, growth type, and year of collection were collected for all isolates. Most mycobacteria were isolated from respiratory specimens (31/57; 54.4%), followed by skin and soft tissues (25/57; 43.8%) and an effusion (1/57; 1.8%).

Of 57 samples, three (3/57; 5.26%) were identified as M. tuberculosis and 54 (94.74%) as NTM by phenotypic methods. In 24.56% (14/57) of the cases, NTMs were classified as SGM, and 70.17% (40/57) as RGM.

3.2. Antimicrobial Resistance Analysis

3.2.1. Analysis of Antimicrobial Resistance in RGMs

All isolates were subjected to susceptibility testing. Amikacin was the most active drug against RGM. No isolates were resistant to amikacin. A single isolate was resistant to tigecycline. For clarithromycin, 7.5% (3/40) of RGM were resistant. The resistance rate to linezolid was 42.5%. Resistance rates to amoxicillin/clavulanic acid, cefepime, and ceftriaxone showed 100% resistance in all isolates in the RGM group. Resistance rates to minocycline (32/40; 80%), imipenem (25/40; 62.5%), tobramycin (27/40; 67.5%), doxycycline (32/40; 80%), ciprofloxacin (33/40; 82.5%), moxifloxacin (20/40; 50%), and trimethoprim/sulfamethoxazole (25/40; 62.5%) were high for RGM (Table 2).

3.2.2. Analysis of Antimicrobial Resistance in SGMs

Resistance to clarithromycin was not observed in any isolate. Resistance rates to ciprofloxacin (6/14; 42.8%), linezolid (6/14; 42.8%), trimethoprim/sulfamethoxazole (5/14; 35.7%), moxifloxacin (5/14; 35.7%), ethionamide (4/14; 28.6%), rifabutin (3/14; 21.4%), amikacin (2/14; 14.3%), rifampicin (2/14; 14.3%), ethambutol (2/14; 14.3%), were less than 50% with the exception of streptomycin (7/14; 50%), isoniazid (7/14; 50%) and doxycycline (11/14; 78.6%) (Table 3).

Differences were observed between the two types of mycobacteria in terms of the ARI (0.54 vs. 0.25). RGM showed a higher rate of resistance compared to SGM. The 3 subspecies of M. abscessus showed the highest rate of resistance among RGM in terms of the ARI (0.65). Within the SGM group, strains belonging to M. avium exhibited a higher resistance rate, as indicated by the ARI (0.41). After analysis of the ARI, M. abscessus-chelone complex showed a higher value than the other complexes, with a value of 0.61.

3.3. Molecular Identification by Phylogenetic Analysis

Phylogenetic trees of the three 16S rRNA, rpoB and hsp65 genes were constructed, first from the genes separately and then using the concatenated data from a single alignment data set. All trees were rooted using the sequences of Nocardia, Gordonia and Tsukamurella species as outgroups, forming a monophyletic and basal group concerning Mycobacterium species, which were internally placed in a monophyletic group. The topologies and most of the interrelationships between Mycobacterium species were very similar and consistent in all the trees constructed, as well as the MP and ML results.

3.3.1. 16S rRNA Gene Phylogenetic Analysis

For the construction of the phylogenetic tree based on the 16S rRNA gene, an alignment matrix of 1400 base pairs or positional characters by 102 sequences or taxa was used. The resulting tree is shown in Figure 1. It was possible to observe the grouping of the Mycobacterium members into four main clades. Two of these clades are composed of slow-growing species, while the other two clades are composed of fast-growing species. Within the clade of fast-growing species, the basal subclade corresponds to the group we refer to as "Fortuitum", which encompasses most of the other fast-growing species and the inner or derived subclade as "Abscessus-Chelonae". Relationships among species in the "Fortuitum" subclade were difficult to discern, unlike members of the "Abscessus-Chelonae" subclade, which formed a monophyletic lineage with a bootstrap value of 100% within the other fast-growing Mycobacterium species and comprises the deepest-branching lineage among Mycobacterium species. For the case of M. phlei, M. thermoresistible and M. monacense species, which compose the fast-growing species and belong to the subclade "Fortuitum", the formation of an earlier branching lineage within the genus Mycobacterium was observed.

Of the two clades of slow-growing mycobacteria, the subclade designated as "Tuberculosis-Avium" contains most of the clinically important mycobacterial species, such as M. tuberculosis and M. avium; and the "Terrae" subclade, which forms a sister and derivative clade with the "Tuberculosis-Avium" subclade.

The percentage of divergence of the 16S rRNA sequences within each subclade was low; the species of the "Terrae" subclade had a value of 1.46%, of the "Tuberculosis-Avium" subclade a value of 1.23%, of the "Abscessus-Chelonae" subclade a value of 0.93%, and of the "Fortuitum" subclade a value of 1.15%. For the case of complexes, the divergence percentages were for M. tuberculosis complex 0.35%, M. avium complex 0.079%, M. intracellulare complex 0.20%, M. abscessus-chelonae complex 0.20% and M. fortuitum complex 0.1%.

The divergence between M. tuberculosis and M. avium complexes was 1.80%, between M. tuberculosis and M. intracellulare was 2%, between M. avium and M. intracellulare was 0.6% and the divergence between M. abscessus-chelonae and M. fortuitum complexes was 3.97 %.

Of the 54 samples analyzed using the 16S rRNA gene, only three isolates were identified to species level: sample L572 as M. marinum, L535 as M. chelonae and L39 as M. fortuitum. The other isolates were identified at the complex or group level, the most predominant being M. abscessus-chelonae complex

3.3.2. rpoB Gene Analysis

For the rpoB gene, an alignment matrix comprising 720 nucleotide positions across 102 sequences (taxa) was used. Phylogenetic analysis of this region yields a different topology than that obtained from the 16S rRNA gene. The rpoB gene tree based on the mycobacterial is shown in Figure 2.

In the rpoB tree topology, Mycobacterium is shown as a monophyletic clade with respect to the reference groups. The formation of two supported and differentiated monophyletic groups was observed. Group A is formed by the species belonging to the M. tuberculosis, M. avium-intracellulare and M. fortuitum complexes and group B is formed by the species of the M. abscessus-chelonae complex, which shows an earlier branching within the genus Mycobacterium. This species location on the tree is related to antibiotic susceptibility; thus the species belonging to M. abscessus-chelonae complex with a higher ARI value form a separate clade with respect to the other mycobacterial species (Figure 2).

The rpoB gene could not differentiate among subspecies or variants within the different mycobacterial complexes. An exception was the M. abscessus complex, in which M. chelonae was clearly distinguished from the M. abscessus complex. However, differentiation among subspecies within the M. abscessus complex was not achieved, with M. bolletii and M. massiliense clustering together as M. bolletii/massiliense. Only M. abscessus subsp. abscessus formed a well-supported monophyletic group with the related isolates (bootstrap = 100%). The percentage of sequence divergence among species within each complex was 1.54% for the M. tuberculosis complex, 0.33% for the M. avium complex, 0.80% for the M. intracellulare complex, 2.76% for the M. abscessus–chelonae group, and 3.20% for the M. fortuitum complex.

The divergence between M. tuberculosis and M. avium complexes was 14.89%, between M. tuberculosis and M. intracellulare was 14.61%, between M. avium and M. intracellulare was 5.6%, and the divergence between M. abscessus-chelonae and M. fortuitum complexes was 12.28%.

Of the 57 clinical isolates, 39 isolates were identified to species level: 28 isolates were identified as M. abscessus, of which 22 were identified to subspecies level as M. abscessus subsp. abscessus; seven isolates were identified as M. chelonae; three as M. fortuitum; and one as M. marinum, the rest were only identified to complex or group level.

3.3.3. hsp65 Gene Analysis

Finally, the analysis was performed with the hsp65 gene, using an alignment matrix of 400 base pairs or positional characters for 105 sequences or taxa. It is possible to visualize the formation of two clearly defined groups: those belonging to the M. abscessus-chelonae complex, M. tuberculosis, M. avium, M. intracellulare, and the second group formed by the M. fortuitum complex that presents a more basal branching (Figure 3). After the statistical analysis of the demographic characteristics, it was identified that the clade differentiation could be related to the type of sample, obtaining that the species belonging to the clade A are more frequently isolated from respiratory samples, while those of the clade B are mainly isolated from non-respiratory samples (p=0.0157).

The percentage of divergence of species within complexes for M. tuberculosis had a value of 1.68%, M. avium a value of 0.11%, M. intracellulare a value of 1.10%, M. abscessus-chelonae a value of 2.98% and M. fortuitum a value of 2.10%. The divergence between M. tuberculosis and M. avium complexes was 7.02%, between M. tuberculosis and M. intracellulare was 7.46%, between M. avium and M. intracellulare was 4.16% and the divergence between M. abscessus-chelonae and M. fortuitum complexes was 6.90%.

Thirty-five (61.4%) isolates were species of the M. abscessus complex, with M. abscessus subsp. abscessus being the most predominant of the subspecies with 22 isolates (62.8%), followed by M. abscessus subsp. massiliense with 5 isolates (14.3%) and M. abscessus subsp. bolletii with 1 isolate (2.8%). Seven (20%) isolates were from the M. chelonae group. Three isolates were identified in the M. tuberculosis complex, but failed to differentiate their individual species. In the case of M. intracellulare complex, only two isolates were differentiated as M. intracellulare, the other 3 isolates were identified only at the complex level. In the case of M. avium complex, only the subspecies M. avium subsp. hominissuis was differentiated among the other subspecies, forming a branch with 7 isolates with a bootstrap of 86%. One isolate failed to differentiate beyond the M. avium complex level.

For M. fortuitum complex, it was possible to differentiate between species of this complex, obtaining 2 isolates belonging to the species M. conceptionense and 3 isolates belonging to M. fortuitum. In addition, 1 isolate was identified as M. marinum.

3.3.4. Analysis of Concatenated Genes

Another phylogenetic tree was constructed based on concatenated sequences of the 3 genes analyzed, achieving a better recovery of the species level, including those belonging to the clade "Fortuitum", which was not achieved with the analysis of the 3 individual genes (Figure 4). It was not possible to differentiate between the subspecies M. bolleti and M. massillense. The topology was similar to that of the phylogenetic tree obtained with the rpoB gene.

Sequence divergence within each complex was 1.09% for M. tuberculosis, 0.15% for M. avium, 0.62% for M. intracellulare, 1.29% for M. abscessus-chelonae and 2.26% for M. fortuitum. The divergence between M. tuberculosis and M. avium complexes was 5.72%, between M. tuberculosis and M. intracellulare was 5.87%, and between M. avium and M. intracellulare was 2.37%. Finally, the divergence between M. abscessus-chelonae and M. fortuitum complexes was 6.4%.

In the concatenated analysis, M. chelonae was identified with seven isolates (bootstrap support 78%), M. fortuitum with three isolates (bootstrap support 58%), M. conceptionense with two isolates (bootstrap support 96%), M. abscessus subsp. abscessus (bootstrap support 70%), M. marinum with one isolate (bootstrap support 100%), and M. intracellulare with three isolates (bootstrap support 91%).

Among the three genes analyzed, the hsp65 typing method showed better agreement for species- and subspecies-level identification, whereas the concatenated analysis provided superior overall phylogenetic topology.

This study identified 35 (64.8%) isolates belonging to the Mycobacterium abscessus complex, making it the most prevalent group. These included Mycobacterium abscessus subsp. abscessus (n = 29), Mycobacterium abscessus subsp. massiliense (n = 5), and Mycobacterium abscessus subsp. bolletii (n = 1). In addition, seven Mycobacterium chelonae isolates (12.9%) were identified as a separate species. The remaining isolates were identified as Mycobacterium avium subsp. hominissuis (n = 7, 12.9%), Mycobacterium fortuitum (n = 3, 5.5%), Mycobacterium intracellulare (n = 2, 3.7%), Mycobacterium conceptionense (n = 2, 3.7%), and Mycobacterium marinum (n = 1, 1.8%). Four isolates were identified only at the complex level: one (1.8%) Mycobacterium avium complex and three (5.5%) Mycobacterium intracellulare complex.

4. Discussion

The incidence of NTM infections in clinical specimens has increased in several regions of the world, so the immediate diagnosis and clear and accurate identification of NTM species causing the disease has become a necessity [35]. In Ecuador, the incidence of atypical mycobacteria in clinical specimens is unknown and the only available data are case reports of some species, including M. scrofulaceum [36], M. terrae [37], M. abscessus subsp. abscessus [38], M. timonense [39] and M. monacense [40]. Therefore, the present study proposes a way to accurately identify mycobacteria and determine the diversity of species isolated from clinical samples.

Molecular techniques have enabled accurate, rapid, and reliable identification of mycobacterial species, with whole-genome sequencing regarded as the gold standard for NTM identification [41]. However, these approaches are expensive, require specialized expertise, and are not readily available in routine laboratories in resource-limited contexts [42]. As an alternative, sequencing of highly conserved yet hypervariable regions, such as the 16S rRNA (~1500 bp), rpoB (~3500 bp), and hsp65 (~1600 bp) genes, has become more widespread in diagnostic workflows [43]. Many studies have focused on partial regions within these genes that balance sufficient conservation with intraspecies variability to achieve accurate identification while reducing costs. In our study, we analyzed the full-length 16S rRNA gene (1500 bp), a 400 bp fragment of hsp65 [24], and a 720 bp fragment of rpoB [10].

Phylogenetic analysis of 16S rRNA, rpoB and hsp65 generally showed topologically similar trees. In the case of the 16S rRNA tree, clades comprising most of the slow-growing mycobacteria were clearly differentiated from the fast-growing species. This differentiation is evident due to the insertion of 12 to 14 nucleotides in the 16S rRNA sequences in SGM, which are often used as a marker to differentiate between RGM and SGM [44]. In the analysis performed in this study, a 12-nucleotide insertion was observed in SGM mycobacteria belonging to the "Tuberculosis" subclade, and a 14-nucleotide insertion belonging to mycobacteria of the "Terrae" subclade (Figure 1).

Regarding the phylogenetic analysis of the rpoB gene, the M. abscessus-chelonae complex separated from the M. fortuitum complex, forming a separate clade. This differentiation may be related to resistance characteristics, as species in the M. abscessus–chelonae complex exhibit higher ARI values than other mycobacteria. Being members of the M. abscessus-chelonae complex, belonging to group B, the bacteria with higher resistance rates compared to mycobacteria grouped with group A [45,46] (Figure 2). Statistical analyses of antimicrobial resistance profiles showed that members of the M. abscessus-chelonae complex presented the highest ARI (0.61) with respect to the other species. This phenotype has generally been attributed to mycobacterial envelope impermeability, with the cell wall of M. abscessus and M. chelonae being 10 to 20 times less permeable than other mycobacteria, due to the higher lipid content, which coincides with M. abscessus-chelonae complex being the most resistant to antibiotics [47].

The rpoB gene contains additional useful information related to rifampicin susceptibility of a particular species (or strain) of mycobacteria [48]. Rifampicin resistance is largely associated with mutations in an 81 bp fragment of the rpoB gene, in addition to other mechanisms such as reduced cell wall permeability or an enhanced efflux pump [15].

RGM, with common representatives including M. abscessus [49] and M. fortuitum [50], exhibit intrinsic resistance to rifampicin. In contrast, SGM such as M. tuberculosis often develop acquired resistance mediated by mutations in the rpoB gene or by overexpression of efflux pumps [51]. Rifampicin resistance in M. abscessus is attributed to its distinctive cell wall, characterized by higher lipid content and the presence of ADP ribosyltransferases, rather than the rpoB mutations that predominate in many other mycobacteria [52]. Although M. fortuitum is classified as an RGM, rifampicin resistance in this species has also been linked to rpoB gene mutations similar to those observed in M. tuberculosis [53], which makes M. fortuitum a useful surrogate for M. tuberculosis in susceptibility testing studies. These differences in resistance phenotypes may reflect the clear phylogenetic separation observed between these species.

The divergence in the phylogenetic analysis of rpoB could also be related to the G+C% content in mycobacterial species. It is known that the rpoB gene allows an efficient estimation of bacterial G+C% content [54], therefore the species found in group A, corresponding to the M. tuberculosis, M. avium, M. intracellulare and M. fortuitum complexes, which present higher G+C% contents, with values of 65.6%, 69.1%, 68% [55] 2021) and 66% [56], respectively, formed a monophyletic group differentiated from members of group B species corresponding to the M. abscessus-chelonae complex that present the lowest G+C% values, with a value of 64.2% [57].

Finally, the hsp65 gene tree showed a clustering of the M. abscessus-chelonae, M. tuberculosis, M. intracellulare and M. avium complex species into a single clade separate from the M. fortuitum complex members. This structured topology between clades could be related to sample type or origin, which upon statistical analysis showed a statistically significant difference of clades vs. sample type (p= 0.0157), suggesting that there is a higher prevalence of species isolated from respiratory samples in clade A in contrast to clade B. Species of the M. abscessus-chelonae complex have been isolated more frequently from respiratory infections, a characteristic shared with the slow-growing mycobacteria, as opposed to M. fortuitum, which produces a wide variety of extrapulmonary diseases, 60% of cases being cutaneous infections [45]. It has been established that three species make up the vast majority of pulmonary infections by atypical mycobacteria; these are: M. avium, M. kansasii and M. abscessus [58], similar to the results obtained in our analysis.

The hsp65 gene encodes a virulence-associated protein that plays an important role in mycobacterial latency [59]. Accordingly, divergence observed in phylogenetic analyses of hsp65 may reflect the pathogenic potential within the genus Mycobacterium, dividing more pathogenic species into group A and less pathogenic species into group B. Slow-growing species, most notably M. tuberculosis and M. leprae, are generally the most pathogenic, whereas most rapidly growing species exhibit lower pathogenicity. An exception is M. abscessus, which, owing to its extensive repertoire of virulence factors, phenotypically resembles slow-growing pathogens [60]. For example, a comparative study of pulmonary disease caused by M. abscessus and M. avium highlighted similarities in pathogenicity [45]. Moreover, at the cellular level, M. abscessus mimics features of pathogenic slow-growing mycobacteria, such as the ability to survive and replicate within macrophages. In contrast, M. fortuitum cannot multiply within macrophages and is rapidly cleared from infected cells [61].

Phylogenetic analyses yielded molecular identification of nine species using 16S rRNA, rpoB and hsp65 genes as molecular markers, with M. abscessus subsp. abscessus being the most prevalent species, followed by M. chelonae and M. avium subsp hominissuis. Not all molecular markers were able to reach species level identification, similar to the result reported by Ledesma et al., [28]. The hsp65 gene yielded optimal results in terms of species and subspecies identification.

Although the 16S rRNA gene is commonly used for bacterial identification and has a large amount of information in the NCBI database [62], for NTMs, the sequences show little interspecific variation, so it may not distinguish closely related species [28], with sequencing of the 16S rRNA region alone being insufficient [63]. In our study, only three isolates were identified to species level using the 16S rRNA gene: M. marinum, M. chelonae, and M. fortuitum.

In this study, the percentage of divergence between the M. abscessus-chelonae, M. tuberculosis and M. bovis, M. fortuitum complexes was lower using the 16S rRNA gene compared to rpoB and hsp65, which corroborates the difficulty in distinguishing the species using 16S rRNA sequences alone. Members of these complexes could be better distinguished by analysis of the hsp65 and rpoB genes.

The rpoB gene, which encodes the b subunit of RNA polymerase, has emerged as a candidate gene for phylogenetic analyses and identification of bacteria, especially when studying closely related strains [54]. Several studies have indicated that the hypervariable region of the rpoB¸ gene located between positions 2300 and 3300 is best suited for phylogenetic identification and discrimination at the species and subspecies level [54]. In our study, 39 isolates were identified to species level, 7 isolates as M. chelonae, 3 as M. fortuitum, 1 as M. marinum and 28 were identified as M. abscessus, of which 22 were identified to subspecies level as M. abscessus subsp. abscessus.

In the case of the closely related subspecies M. abscessus, M. bolletii and M. masiliense, the rpoB region alone was not sufficient for good discrimination [64]. Finding divergence values of 0.0% between M. bolletii and M. masiliense subspecies and 2.03% between M. abscessus and M. bolletii, M. masiliense. Between-group divergence values <2% in rpoB sequences indicate that these groups belong to the same species [65]. This was more marked for M. massiliense and M. bolletii, than for M. abscessus. The rpoB sequences in the strains of the M. massiliense and M. bolletii groups did not diverge. Therefore, according to Adékambi et al. [65] criteria, M. massiliense and M. bolletii for this analysis do not constitute different species. The importance of correctly identifying the subspecies of M. abscessus lies in their antimicrobial profile. M. abscessus and M. bolletii are known to have inducible macrolide resistance due to the presence of the functional erm(41) gene, while M. massiliense is susceptible because this gene is truncated and non-functional [66]. Therefore, the correct identification of M. abscessus subspecies is of clinical relevance.

The hsp65 gene is a good option for the identification of RGM, because it shows more variability than the 16S rRNA and rpoB genes [67]. In this study, it was possible to identify the subspecies of M. abscessus by analyzing the hsp65 gene, obtaining a percentage of divergence of 0.64% between the subspecies M. bolletii and M. masiliense, 1.28% between M. abscessus and M. masiliense and 1.41% between M. abscessus and M. bolletii. This gene is part of the 60 kDa heat shock protein family and has been considered a useful phylogenetic marker in several bacterial genera, due to its highly conserved primary structures [68]. The presence of a single copy of the gene in the genome means that it is not easily transferred from one bacterium to another, unlike the rpoB gene [68], it is suitable for phylogenetic studies of species or closely related strains [70]. Analysis based on the hsp65 gene identified 53 isolates to species level, seven belonging to M. chelonae, 28 to M. abscessus (22 to M. abscessus; five to M. massiliense; one to M. bolletii), three to M. tuberculosis, two to M. intracellulare, seven to M. avium-hominissuis, one to M. marinum, two to M. conceptionense and three to M. fortuitum.

The partial paraphyletic location of some fast-growing species belonging to the subclade "Fortuitum" and the subclade "Terrae" present in the tree based on 16S rRNA, rpoB and hsp65 were resolved in the tree constructed with concatenated nucleotide sequences, therefore, the increasing availability of markers spanning the entire bacterial genome allows the construction of more robust phylogenetic trees, thus providing a reliable understanding of the relationships within the genus Mycobacterium [71]. The concatenated analysis yielded improved topology and enabled identification of most isolates at the species and subspecies levels, with the exception of two isolates that could only be resolved at the complex level, consistent with findings from other genome-wide studies. Thus, increased genomic information facilitates more accurate inference of evolutionary relationships, clearer demarcation of taxa, and enhanced understanding of interrelationships among these mycobacterial species [44].

5. Conclusions

Synapomorphic characters (shared positional homologies of an immediate common ancestor) for the identification and reliable demarcation of different monophyletic clades of mycobacteria varied among genes analyzed, with correlation observed between the hsp65 gene and pathogenicity phenotypes and sample type, rpoB with resistance profile, G-C% content, and 16S rRNA with growth rate.

The 16S rRNA gene was more conserved and showed less variability compared to the rpoB and hsp65 genes. The hsp65 and rpoB genes showed greater variability, which allowed an optimal identification approach at the species level and even at the subspecies level, as in the case of the hsp65 gene. These phylogenetic results emphasize that molecular markers other than the 16S rRNA gene should be applied to obtain an accurate phylogenetic and evolutionary scheme, and concatenated analyses of genes with greater phylogenetic information at specific levels.

Clinical samples showed phylogenetic homology with the groups of M. abscessus (Mycobacterium abscessus subsp. abscessus, Mycobacterium abscessus subsp. massiliense, Mycobacterium abscessus subsp. bolletii), M. avium (Mycobacterium avium subsp. hominissuis), M. fortuitum, M. intracellulare and M. marinum.