Comparative genomic analysis of established and emerging anti-tubercular drug targets reveals evolutionary conservation and antimicrobial resistance patterns in Mycobacterium tuberculosis isolates from Indonesia

Multidrug-resistant tuberculosis (MDR-TB) remains a major public health threat in Indonesia, which ranks second globally in tuberculosis burden. Whole-genome sequencing (WGS) studies of drug-resistant Mycobacterium tuberculosis in Indonesia have so far been analyzed largely in isolation, and no previous study has systematically compared the evolutionary conservation of classical, clinically validated resistance-associated genes with newly proposed metabolic drug targets within a single, geographically diverse Indonesian isolate collection. We compiled 444 publicly available Mycobacterium tuberculosis whole genomes from four Indonesian whole-genome sequencing studies spanning Java and Papua and analyzed them using TBProfiler for lineage assignment and genotypic drug-resistance prediction. Following exclusion of nine isolates for which an unambiguous single lineage could not be assigned, gene-level sequences for ten classical resistance-associated genes (rpoB, katG, inhA, fabG1, embB, pncA, gyrA, gyrB, rpsL, rrs) and ten emerging metabolic drug-target genes (guaA, guaB2, glnA1, gltB, dapB, folP1, bioA, panC, aroG, ilvB1) were reconstructed for 435 isolates and compared by multiple sequence alignment, variant calling, and functional annotation. (East-Asian) predominated (46.2% of isolates) and was disproportionately associated with drug-resistant genotypes (χ² = 49.07, df = 2, p < 0.001; Lineage 2 versus Lineage 4 odds ratio for any genotypic resistance 2.69, 95% CI 1.77–4.07). Overall, 28.2% of isolates were genotypically multidrug-resistant and 5.4% were pre-extensively drug-resistant, with no mutations associated with bedaquiline, linezolid, or clofazimine resistance detected. All twenty target genes were highly conserved (89.4–99.6% invariant sites), consistent with the clonal population structure of M. tuberculosis. However, the two panels differed substantially: classical resistance-associated loci were significantly less conserved than emerging targets (mean 96.47% versus 99.24%; p < 0.001), carried a significantly higher SNP density (20.70 versus 6.58 SNPs per kilobase; p = 0.002) and per-site Shannon entropy (p = 0.009), and accounted for 537 of the 670 variable sites detected across the panel despite comprising a similar total coding length. The non-synonymous-to-synonymous substitution ratio was also significantly higher among classical loci (median 1.67 versus 0.96; p = 0.033), driven by the prodrug-activating enzymes pncA (49.0) and katG (10.4) and by rpsL (6.0) and rpoB (2.08). Within the emerging panel, absolute variable-site counts were uniformly low (4–47 sites per gene) and did not support discrimination between individual candidate genes. Across a geographically diverse Indonesian Mycobacterium tuberculosis population, therapeutically unexploited metabolic drug targets are markedly more conserved and carry a substantially lower mutational burden than classical resistance-associated loci, whose elevated non-synonymous variation is concentrated in genes whose inactivation confers resistance. These findings support the emerging metabolic panel collectively as a low-variability target space for future anti-tubercular development, and establish evolutionary conservation combined with substitution-pattern analysis as a quantitative framework for target-space prioritization and genomic surveillance in high-burden tuberculosis settings.

Tuberculosis (TB) remains one of the leading causes of death from a single infectious agent worldwide, and drug-resistant disease continues to undermine global control efforts; it is estimated to cause approximately 13% of all antimicrobial-resistance-attributable deaths, driven by both de novo resistance acquisition and ongoing transmission of already-resistant strains (1). Multidrug-resistant tuberculosis (MDR-TB), which is resistant to at least rifampicin and isoniazid, requires prolonged, costly, second-line treatment and carries substantially worse outcomes than drug-susceptible disease. Indonesia bears a disproportionate share of this burden, ranking second worldwide in absolute TB case numbers and contributing substantially to the global MDR/RR-TB caseload (2). Despite this, whole-genome sequencing (WGS)-based genomic epidemiology and molecular surveillance of drug-resistant Mycobacterium tuberculosis (Mtb) has historically been limited in Indonesia, with published genomic studies concentrated in Java and involving comparatively small isolate collections (3, 4), and genomic characterization of drug-resistant Mtb from other Indonesian regions, such as Papua, remains even more limited (5). This gap constrains the ability of national and regional TB programmes to track circulating lineages, monitor resistance-mutation emergence, and identify candidate loci for future diagnostic or therapeutic development — precisely the kind of surveillance gap that genomic epidemiology approaches to emerging tropical pathogens are intended to close, particularly in resource-diverse settings where scalable, sequencing-based tools can substitute for slower phenotypic methods.

Resistance to first- and second-line antitubercular drugs arises almost exclusively through chromosomal mutations, and the molecular basis of resistance to each major drug class is well characterized (6), enabling WGS-based tools, such as TBProfiler, to predict genotypic drug resistance and phylogenetic lineage directly from sequencing reads (7, 8). However, because these classical resistance-associated genes have been under sustained clinical selective pressure for decades, each is a site of continual resistance mutation acquisition, and the anti-TB drug discovery pipeline has disproportionately revisited this same small, increasingly compromised target set (9). Expanding the target space to genetically validated but therapeutically unexploited metabolic pathways, including de novo purine biosynthesis (guaA, guaB2), nitrogen assimilation and amino acid metabolism (glnA1, gltB), the diaminopimelate pathway supplying peptidoglycan precursors (dapB), folate biosynthesis (folP1) (10), biotin biosynthesis (bioA) (11), pantothenate/coenzyme A biosynthesis (panC), the shikimate pathway (aroG), and branched-chain amino acid biosynthesis (ilvB1), is therefore a priority, particularly as several of these pathways have recently been shown to be essential for Mtb survival during mammalian infection using genome-wide transposon mutagenesis (12).

Several of these pathways also carry independent clinical or virulence-related relevance that strengthens the case for their prioritization, beyond bare genetic essentiality. folP1 encodes the mycobacterial dihydropteroate synthase and is the physiological analogue of the classical sulfonamide target exploited in other bacterial pathogens; mycobacterial folate metabolism is also directly implicated in the mechanism of the second-line agent para-aminosalicylic acid (10), providing an existing chemical starting point. bioA -dependent biotin biosynthesis supplies the cofactor required by the acyl-CoA carboxylases of mycolic acid biosynthesis, linking this pathway directly to cell-envelope integrity (11). dapB catalyzes a committed step of the diaminopimelate pathway that supplies the meso-diaminopimelate cross-linking residue of the mycobacterial peptidoglycan, a pathway that is absent from humans and has a long history as an antibacterial target. glnA1 -encoded glutamine synthetase contributes to the poly-L-glutamate/glutamine layer of the mycobacterial cell wall and is required for efficient intracellular growth in human macrophages and for virulence in the guinea pig model of pulmonary TB, distinguishing it functionally from its human orthologue (13). ilvB1 -dependent branched-chain amino acid biosynthesis is dispensable in mammals, which obtain these amino acids from the diet; genetic inactivation of ilvB1 both confers auxotrophy and attenuates virulence in a murine infection model, underscoring its potential as a target with a wide therapeutic window (14). Amino acid biosynthesis more broadly has attracted renewed interest as an anti-tubercular target space, and inhibitors of several of these pathways are in early development (15). The de novo purine biosynthetic pathway is likewise attractive (IMP dehydrogenase) catalyze consecutive committed steps in guanine nucleotide synthesis, are required for Mtb growth in the absence of an exogenous purine source, and differ sufficiently from their human counterparts to offer a selectivity window (12). Together with the shikimate pathway (aroG), which is absent from humans, this panel spans a functionally diverse set of pathways with plausible routes to selective, low-toxicity inhibition.

A key, comparatively underexplored criterion for prioritizing such candidates is evolutionary conservation across genetically diverse clinical isolates: a target gene that tolerates little sequence variation across a large, phylogenetically diverse population is less likely to already harbor or rapidly acquire resistance-conferring polymorphisms once therapeutically exploited (9). While large-scale mutation catalogues exist for classical resistance genes, comparable population-level conservation and mutation frequency data for the emerging metabolic target panel are scarce. To our knowledge, no previous study has directly compared the evolutionary conservation and mutational burden of classical and emerging anti-tuberculosis drug targets within a single, well-characterized Indonesian Mtb isolate collection. We addressed this gap using publicly available Mtb whole genomes from Indonesia (Java and Papua), first characterizing lineage composition and genotypic drug-resistance profile, then comparing sequence conservation, SNP density, substitution patterns, and mutation frequency between ten classical and ten emerging drug-target genes, to build a quantitative, population-genomic evidence base for target prioritization and to provide a comparative genomic framework complementing existing surveillance of MDR-TB across these Indonesian settings.

Study population and genomic dataset

A total of 444 Mtb whole-genome sequencing datasets from Indonesian clinical isolates were identified across four publicly available NCBI BioProjects — PRJNA1433044 (Koesoemadinata et al. (16), Bandung, West Java), PRJNA1150781 (Djunaedy et al. (5), Jayapura, Papua), PRJNA633244 (Tania et al. (3), Jakarta metropolitan area, Banten, West Java, and Jogjakarta), and PRJEB24463 (Aung et al. (17), Bandung, West Java) — summarized in Supplementary Table S1. BioProjects were included if they comprised Mtb complex genomes with a recorded country of origin of Indonesia and publicly accessible paired-end short-read sequencing data; the resulting cohort therefore spans a range of genotypic drug-resistance categories, from pan-susceptible to pre-extensively drug-resistant, as characterized independently using TBProfiler in this study. No isolates from Bali were available in, or included from, any of the four contributing BioProjects.

The four contributing BioProjects also differed in their original study design and inclusion criteria: the source studies for PRJNA633244 and PRJEB24463 specifically recruited patients with phenotypically or genotypically drug-resistant tuberculosis, whereas the source studies for PRJNA1150781 and PRJNA1433044 sampled isolates on the basis of lineage characterization and active tuberculosis-case-contact finding, respectively, without restricting inclusion to drug-resistant disease. Consequently, the compiled cohort analyzed here is not uniformly drug-resistant by original study design, and the resistance categories reported in Table 1 and section 3.1 reflect genotypic classification performed independently in the present study using TBProfiler on the raw sequencing reads, rather than any phenotypic or genotypic resistance designation used for inclusion in the original source studies; the two classifications need not agree in every case, since phenotypic and genotypic resistance calls can diverge when the causal resistance variant is not captured in TBProfiler’s mutation catalogue, or vice versa.

Paired-end FASTQ files were retrieved from the NCBI Sequence Read Archive (SRA) using the NCBI SRA Toolkit v3.2.1 (prefetch followed by fasterq-dump) (18). No separate read quality-control or trimming step was applied prior to analysis; read trimming was performed internally by TBProfiler as described in section 2.2. Per-isolate BioSample and SRA run accession numbers are available from each BioProject via the NCBI SRA Run Selector and are summarized in Supplementary Table S1. Because all data were previously deposited, de-identified, and publicly accessible sequences, no new specimens were collected, and a formal institutional ethical review was not required.

Isolate exclusions. Nine of the 444 isolates were excluded from all sequence-level analyses because TBProfiler could not assign an unambiguous single main lineage: eight returned no lineage call and one returned a four-way mixed call spanning two species (M. orygis; lineages 2, 4 and 6). These nine isolates also carried markedly elevated numbers of non-resistance-associated variants relative to the cohort, consistent with mixed or contaminated sequencing libraries; two carried between 80,000 and 210,000 whole-genome variant calls, approximately 60–160 times the cohort median. Retaining them would have introduced substantial artefactual sequence diversity into the gene alignments. The nine excluded isolates, with their lineage calls, resistance classifications and variant counts, are listed in Supplementary Table S5. All sequence-level results below therefore describe 435 isolates. The cohort description in section 3.1 and Table 1, which derives directly from TBProfiler output rather than from the alignments, reports all 444 isolates as profiled; the nine isolates with mixed or undetermined lineage calls are shown as a separate category in Table 1 and were likewise excluded from the lineage–resistance association analysis.

Lineage assignment and resistance-genotype prediction

Reads for all 444 isolates were analyzed with TBProfiler v4.0.3 (7, 8), using the bundled tbdb mutation library (bioconda distribution, dated 23 November 2021), run with default parameters (tb-profiler profile -1 <R1> -2 <R2> -t 20 -p <sample> --txt). Within this workflow, reads were trimmed with Trimmomatic v0.40 (19) using TBProfiler’s default settings, aligned to the H37Rv reference genome (GenBank/RefSeq NC_000962.3) (20) with BWA-MEM v0.7.18 (21), and variants called with freebayes v1.3.8, the default caller in this TBProfiler version; no --mapper, --caller or --no_trim options were specified. Detected variants were compared against TBProfiler’s curated resistance-mutation library to assign predicted genotypic drug-resistance profiles and phylogenetic lineage and sub-lineage classifications. No masking of PE/PPE gene family regions was applied; none of the twenty target genes analyzed in this study lies within a PE/PPE region.

Because TBProfiler restricts variant calling to drug-resistance candidate regions unless the --call_whole_genome option is specified, genome-wide variant calls required for the emerging target panel were generated separately from the TBProfiler alignments using SAMtools v1.10 and BCFtools v1.10.2 (22): bcftools mpileup -f H37Rv.fasta <sample>.bam | bcftools call -mv -Oz, followed by bcftools consensus to generate a per-isolate whole-genome consensus sequence. No minimum-depth mask was applied at the consensus step; positions with insufficient read coverage therefore inherit the H37Rv reference base rather than being reported as missing (see Limitations, section 4.3). bcftools consensus applies IUPAC ambiguity codes at heterozygous genotype calls by default, and such codes were retained in the consensus sequences; their treatment in downstream analyses is specified in Section 2.3.

Target gene panel and comparative genomic analysis

For clarity, we refer to the ten established resistance-associated genes throughout as “classical resistance-associated genes” rather than “drug targets,” since several act indirectly on drug activity (e.g., via prodrug activation for katG / pncA, or ribosomal binding for rpsL / rrs) rather than as the direct enzymatic target of the corresponding drug; the ten metabolic genes are referred to as “emerging targets,” consistent with their status as genetically validated but therapeutically unexploited prospective targets. H37Rv locus tags and the anti-tubercular drug or metabolic pathway associated with each of the twenty target genes are summarized in Supplementary Table S4.

Coding sequence coordinates for all twenty genes were extracted from the H37Rv RefSeq annotation (GCF_000195955.2_ASM19595v2), including strand, and used to construct a single coordinate set applied identically to the reference and to every isolate. For each isolate and gene, the corresponding H37Rv reference interval was extracted with SAMtools faidx and the isolate’s variants within that interval applied with bcftools consensus, so that all gene sequences remain anchored to H37Rv coordinates and are unaffected by insertion or deletion events elsewhere in the genome. Sequences for genes on the minus strand were reverse-complemented, together with the reference, after extraction. Each gene-specific dataset (435 isolates plus the H37Rv reference; 436 sequences per gene) was aligned using MAFFT v7.525 (23) with the --auto strategy, and single-nucleotide polymorphisms were identified from each alignment using SNP-sites v2.5.1 (24). Sixteen of the twenty alignments were exactly the length of the corresponding coding sequence; the remaining four (katG, pncA, rpoB, rrs) exceeded it by 1–9 columns, corresponding to short insertions present in a small number of isolates.

Alignment gaps were excluded from per-column allele counts, so that conserved/variable classification and Shannon entropy were computed over non-gap characters only. IUPAC ambiguity codes arising from heterozygous consensus calls were retained as distinct character states in the conservation and entropy analyses, so that a position at which the only non-reference character is an ambiguity code is classified as variable; conservation percentages should therefore be read as conservative (lower-bound) estimates. Conservation was calculated as the proportion of invariant columns among alignment columns at which the H37Rv reference carries a base, so that insertion columns do not contribute to either numerator or denominator and per-gene values are directly comparable. For the transition/transversion analysis, which is most sensitive to ambiguous calls, the ratio was computed only over biallelic positions at which every observed allele was an unambiguous A, C, G or T.

For each gene we calculated the proportion of conserved (invariant) reference-anchored sites, SNP density (unambiguous SNP sites per kilobase of coding sequence), transition/transversion (Ti/Tv) ratio, per-site Shannon entropy, and the proportion of variable sites at which the alternate allele reached ≥50% frequency across the cohort (referred to here as a high-frequency variant site rather than a “mutational hotspot,” since a high-frequency variant can also arise through clonal expansion of a lineage-associated polymorphism rather than through recurrent, independent mutation; distinguishing these possibilities would require examining variant frequencies stratified by lineage, which we identify as a direction for future work).

Variant coordinates were translated from alignment positions to H37Rv chromosomal coordinates by mapping through the gapped H37Rv row of each alignment, with alleles complemented for minus-strand genes, and verified base-by-base against the H37Rv reference sequence before annotation. All variants were then annotated with SnpEff v5.3a (25) using the Mycobacterium_tuberculosis_h37rv database to classify each as synonymous, non-synonymous (missense), or loss-of-function (stop-gained, stop-lost, start-lost). Only annotations assigned to the target gene itself were counted, and each consequence class was counted once per variable site rather than once per alternate allele. A non-synonymous-to-synonymous (NS:S) substitution-count ratio was calculated for each protein-coding gene as a proxy for the direction and intensity of selective pressure. We note that this NS:S variant-count ratio is a simplified metric and is not equivalent to a codon-normalized dN/dS estimate (26, 27): it does not correct for the number of possible synonymous versus non-synonymous sites at each codon, gene length, codon usage, or the underlying trinucleotide mutation spectrum. We therefore interpret an elevated NS:S ratio as evidence of relative enrichment of non-synonymous variation rather than as direct, quantitative evidence of positive selection.

Statistical analysis

Locus-level metrics were summarized separately for the classical and emerging gene panels (mean ± standard deviation, with medians reported where distributions are skewed) and compared using the two-tailed Mann–Whitney U test, following confirmation of non-normal distributions with the Shapiro–Wilk test; significance was set at p < 0.05. For each panel-level comparison we additionally report the rank-biserial correlation as a standardized effect size and the Hodges–Lehmann estimate of the median difference between panels with its 95% confidence interval (Moses method), so that the precision of each comparison can be assessed alongside the p-value.

Association between phylogenetic lineage and drug-resistance category (any genotypic resistance versus pan-susceptible) was assessed using a chi-square test of independence, with pairwise post-hoc comparisons between individual lineages performed using Fisher’s exact test, reporting the corresponding odds ratio and 95% confidence interval (Woolf’s method). The nine isolates returning mixed or undetermined main-lineage calls were excluded from this analysis, which was therefore performed on 435 isolates; Cramér’s V is reported as an effect size for the overall association.

All statistical analyses were performed in Python 3 using SciPy and pandas; analysis scripts are available from the corresponding author on request. Raw sequencing reads for all 444 isolates are publicly available through the NCBI Sequence Read Archive under BioProjects PRJNA1433044, PRJNA1150781, PRJNA633244, and PRJEB24463 (Supplementary Table S1). The H37Rv reference genome is available under the GenBank/RefSeq accession number NC_000962.3.

Cohort composition, lineage distribution, and drug-resistance profile

TBProfiler-based genotypic drug susceptibility prediction and lineage assignment were completed for all 444 Mtb genomes (Table 1). Overall, 229 isolates (51.6%) were genotypically sensitive at all tested loci, while 125 (28.2%) met the criteria for MDR, 34 (7.7%) for RR-TB (rifampicin-resistant, isoniazid-susceptible), 24 (5.4%) for pre-extensively drug-resistant TB (Pre-XDR), 14 (3.2%) for isoniazid mono-resistance (HR-TB), and 18 (4.1%) for other mono-resistance patterns. All resistance categories reported here reflect the genotypic classification generated independently in this study; this is not necessarily identical to any phenotypic or genotypic resistance designation used in the original source studies for patient or isolate selection.

By individual drug, the predicted resistance frequency was highest for rifampicin (41.2%) and isoniazid (36.7%), followed by ethambutol (21.2%), pyrazinamide (16.2%), streptomycin (13.5%), and fluoroquinolones (6.5%); no isolate carried a resistance-conferring mutation for bedaquiline, linezolid, clofazimine, or delamanid.

Lineage 2 (East-Asian) predominated, accounting for 205 isolates (46.2%), followed by Lineage 4 (Euro-American, 174 isolates, 39.2%) and Lineage 1 (Indo-Oceanic, 56 isolates, 12.6%). Nine isolates (2.0%) returned mixed or undetermined lineage calls and were excluded from all subsequent analyses (section 2.1; Supplementary Table S5). The full sub-lineage-level breakdown of the cohort is shown in Supplementary Figure S1. Lineage 2 was disproportionately associated with resistance at the isolate level: only 35.6% of Lineage 2 isolates were pan-susceptible, compared with 59.8% of Lineage 4 and 83.9% of Lineage 1 isolates. Within Lineage 2, 40.0% (82/205) of isolates were genotypically MDR and 8.3% (17/205) were Pre-XDR, compared with 21.8% (38/174) and 3.4% (6/174) of Lineage 4 isolates and 5.4% (3/56) and 1.8% (1/56) of Lineage 1 isolates. Expressed with the resistance category rather than the lineage as the denominator, Lineage 2 accounted for 65.6% (82/125) of all MDR isolates and 70.8% (17/24) of all Pre-XDR isolates, despite comprising 46.2% of the cohort.

A chi-square test of independence confirmed a statistically significant association between lineage and resistance category (any genotypic resistance versus pan-susceptible; χ² = 49.07, df = 2, p < 0.001; Cramér’s V = 0.34; n = 435). Pairwise Fisher’s exact tests showed that Lineage 2 isolates were substantially more likely than Lineage 1 isolates (odds ratio 9.44, 95% CI 4.38–20.36, p < 0.001) and than Lineage 4 isolates (odds ratio 2.69, 95% CI 1.77–4.07, p < 0.001) to carry any genotypic resistance; Lineage 4 isolates were in turn more likely than Lineage 1 isolates to do so (odds ratio 3.51, 95% CI 1.62–7.63, p = 0.001).

The canonical resistance mutations rpoB Ser450Leu (56.8% of rifampicin-resistant isolates) and katG Ser315Thr (65.6% of isoniazid-resistant isolates) predominated, and 54 isolates (12.2% of the cohort) carried both simultaneously, the single most common composite MDR-defining genotype. The full isolate-by-drug co-resistance architecture is shown in Supplementary Figure S2.

Classical resistance loci are significantly less conserved than emerging targets

Across all twenty target genes and 435 isolates, sequence conservation was uniformly high, ranging from 89.39% (rrs) to 99.62% (bioA), with a total of 670 variable sites detected across 33,871 bp of aligned coding sequence (Table 2; full per-gene statistics in Supplementary Table S2). This is consistent with the highly clonal population structure of M. tuberculosis, in which genome-wide diversity is low relative to most bacterial pathogens.

Within this narrow overall range, however, the two panels separated clearly and consistently. Classical resistance-associated loci were significantly less conserved than emerging metabolic targets (mean 96.47% ± 3.12% versus 99.24% ± 0.26%; median 97.97% versus 99.25%; Mann–Whitney U = 2.0, p < 0.001; Hodges–Lehmann estimate of the median difference −1.33 percentage points, 95% CI −2.16 to −0.86; rank-biserial r = −0.96; Figure 1). The classical panel also carried a significantly higher SNP density (mean 20.70 ± 19.86 versus 6.58 ± 2.75 SNPs per kilobase; median 15.22 versus 6.18; U = 92.0, p = 0.002; median difference +8.26, 95% CI + 3.48 to +14.60; r = +0.84) and a significantly higher mean per-site Shannon entropy (0.0020 ± 0.0018 versus 0.0006 ± 0.0005; U = 85.0, p = 0.009; r = +0.70).

In absolute terms, the ten classical loci accounted for 537 of the 670 variable sites detected (80.1%), despite comprising 17,611 bp of coding sequence against 16,260 bp for the emerging panel — a comparable total length. The four most variable genes in the entire panel were all classical loci or the ribosomal RNA gene: rrs (163 variable sites), rpoB (77), gyrA (66) and katG (60). The least variable were all emerging targets: dapB (4), bioA (5), panC (5) and aroG (7). Per-gene SNP densities across the full panel are shown in Figure 2.

The Ti/Tv ratio, computed over biallelic positions at which all observed alleles were unambiguous, was 1.66 pooled across all twenty genes (211 transitions, 127 transversions among 338 biallelic sites), consistent in direction with the transition bias reported from phylogenetically corrected analyses of the Mtb mutation spectrum (26). Per-gene ratios ranged from 0.50 to 8.00, but rest on small numbers of qualifying sites (3–44 per gene), and did not differ significantly between panels (classical 1.43 ± 0.47, n = 9; emerging 2.96 ± 2.42, n = 8; U = 22.5, p = 0.208). We therefore report Ti/Tv descriptively and do not interpret per-gene values as estimates of the underlying mutational spectrum.

Of the 670 variable sites, 568 (84.8%) had an alternate allele frequency below 1%, and only six sites across the entire panel reached ≥50% frequency (three in rpoB, one each in gyrA, katG and guaB2), underscoring that the variation detected is overwhelmingly rare and unevenly distributed (Supplementary Table S3).

Non-synonymous enrichment at classical loci is concentrated in prodrug-activating enzymes

SnpEff functional annotation of variants in the 19 protein-coding target genes (rrs, a non-coding 16S rRNA gene, was excluded from the coding-consequence analysis) classified 312 sites as non-synonymous (60.6%), 195 as synonymous (37.9%), and 8 as loss-of-function (1.6%; full per-gene breakdown in Supplementary Table S2).

The NS:S ratio was significantly higher among the nine protein-coding classical resistance-associated genes than among the ten emerging targets (median 1.67 versus 0.96; U = 71.5, p = 0.033; rank-biserial r = +0.59; Figure 3). The distribution within the classical panel was, however, strongly bimodal, and the effect is driven by a specific and mechanistically interpretable subset of genes. The two prodrug-activating enzymes showed by far the highest ratios in the entire panel — pncA (49 non-synonymous versus 1 synonymous variant; NS:S = 49.0) and katG (52 versus 5; NS:S = 10.4) — followed by rpsL (6.0), rpoB (2.08) and inhA (1.67). For pyrazinamide and isoniazid, resistance arises through loss of function of the activating enzyme, so amino-acid-altering and disruptive changes are directly selected; both genes additionally carried loss-of-function variants (three each), which together with rpoB and panC account for 7 of the 8 loss-of-function variants detected across the entire panel.

By contrast, the ten emerging targets formed a tight cluster around parity (range 0.25–4.00; median 0.96; interquartile range 0.50–1.00), with no gene showing the strong non-synonymous excess characteristic of the drug-activating loci. Three classical genes — gyrB (0.54), fabG1 (0.50) and embB (1.39) — fell within the emerging range, consistent with their status as direct enzymatic targets under strong purifying selection, at which resistance arises through a small number of specific substitutions that do not substantially shift a gene-wide ratio.

The emerging target panel shows uniformly low mutational diversity

All ten emerging metabolic targets were highly conserved, with conservation ranging narrowly from 98.81% (folP1) to 99.62% (bioA) and SNP density from 3.81 to 11.86 SNPs per kilobase (Table 2). Absolute variable-site counts were correspondingly low: four genes carried fewer than eight variable sites across 435 isolates (dapB 4, bioA 5, panC 5, aroG 7), and the most variable emerging gene, gltB, carried 47 variable sites across 4,584 bp — a density still below that of eight of the ten classical loci.

We had originally intended to prioritize individual candidate genes within this panel by combining above-median conservation with below-median, length-normalized SNP density. Applying this criterion to the present data, we do not consider gene-level discrimination to be statistically supportable: the range of variable-site counts across the emerging panel (4 to 47) is small in absolute terms, differences between the most and least conserved emerging genes amount to 0.8 percentage points of conservation, and the ranking is therefore sensitive to the addition or removal of individual variants. We report the per-gene values in Table 2 for completeness, but interpret the emerging panel as a collectively low-variability target space rather than nominating individual genes for prioritized validation on the basis of these data.

Principal findings

In this Indonesian Mycobacterium tuberculosis collection, compiled from four publicly available whole-genome sequencing studies spanning Java (Bandung, the Jakarta metropolitan area, Banten, and Jogjakarta) and Papua (Jayapura), Lineage 2 (East-Asian) dominated overall and carried a disproportionate resistance burden relative to Lineage 1 and Lineage 4, an association that was statistically significant and of moderate magnitude (Cramér’s V = 0.34). This distribution is consistent with the lineage patterns reported in the source studies themselves: the East-Asian lineage was likewise the dominant lineage among MDR isolates from Java, accounting for 63.3% of MDR strains in the source study of the Java isolates used here (3), whereas Lineage 1 predominated among the Papua isolates, where only a small minority carried resistance-associated variants (5), attributed by that study to inconsistent access to effective treatment rather than to reduced transmission (5).

The principal comparative finding is that classical resistance-associated loci and therapeutically unexploited metabolic targets differ substantially and consistently in their population-level sequence variation. Across five independent metrics — conservation, SNP density, Shannon entropy, absolute variable-site count, and the NS:S ratio — the classical panel carried more variation than the emerging panel, with effect sizes ranging from moderate to very large (rank-biserial r = 0.59 to 0.96). The classical loci accounted for four fifths of all variable sites detected despite representing a comparable total length of coding sequence.

Two features of the data indicate that this difference reflects drug-driven selection rather than a general property of the genes. First, the excess variation at classical loci is disproportionately non-synonymous, and the excess is concentrated precisely where pharmacology predicts: pncA and katG, whose products activate pyrazinamide and isoniazid respectively and whose inactivation confers resistance, show NS:S ratios of 49.0 and 10.4, orders of magnitude above every emerging target. Both also carry loss-of-function variants, which are otherwise almost absent from the panel. This is the expected signature of selection for enzyme disruption. A comparable enrichment of non-synonymous variation at loci under direct clinical drug pressure has been reported previously in M. tuberculosis (28). Second, classical genes that are themselves the direct enzymatic target of one — gyrB, fabG1, embB — fall within the emerging-panel range for NS:S, consistent with strong purifying selection on essential catalytic function, at which resistance arises through a small number of specific substitutions that do not shift a gene-wide ratio. The classical panel is therefore not uniformly variable; it is bimodal, and the split falls along a mechanistic line.

The absence of comparable variation at the emerging loci is pharmacologically plausible. None of the agents used in Indonesia’s drug-resistant TB regimens — the first-line agents rifampicin, isoniazid, ethambutol, and pyrazinamide; the injectable agents streptomycin, amikacin, and kanamycin; fluoroquinolones; or the repurposed and newer agents bedaquiline, linezolid, clofazimine, and delamanid — acts on the purine, nitrogen-assimilation, diaminopimelate, folate, biotin, pantothenate, shikimate, or branched-chain amino acid pathways represented by the emerging panel. This holds even in Papua, the most resource-limited setting contributing isolates to this cohort, where the source study reported markedly lower rates of resistance-associated variants than in the Java-derived isolates (5).

We note, however, that our NS:S metric is a simple variant-count ratio rather than a codon-normalized dN/dS estimate (26, 27), and this interpretation should be read as descriptive evidence of relative non-synonymous enrichment rather than a formal, quantitative test of positive selection.

Implications for drug-target prioritization and genomic surveillance

The practical implication of these results is at the level of the target space rather than the individual gene. The ten emerging metabolic targets examined here are, collectively, markedly more conserved and less mutationally diverse across a genetically heterogeneous Indonesian population than the loci that current anti-tubercular chemotherapy engages. For target discovery, this is the relevant comparison: a candidate drawn from this pathway set begins from a position of substantially lower standing sequence variation, and therefore lower pre-existing potential for target-site resistance, than the classical target set.

We deliberately do not nominate individual genes within the emerging panel. The absolute number of variable sites per emerging gene is small — four genes carry fewer than eight variable sites across 435 isolates — and differences of this magnitude cannot support a reliable ranking. This is worth stating explicitly, because conservation-based prioritization schemes are attractive precisely where the underlying counts are lowest, and the apparent precision of a percentage conservation figure can obscure how few observations lie behind it. In a pathogen as clonal as M. tuberculosis, distinguishing between genes that are 99.2% and 99.5% conserved requires either a much larger and more phylogenetically diverse cohort, or an approach that does not rest on standing variation alone, such as saturating mutagenesis or structure-guided assessment of the mutational tolerance of the ligand-binding site.

More broadly, this analysis illustrates the value of applying population-scale WGS data to genomic epidemiology and molecular surveillance questions in tropical, high-burden settings that have historically been underrepresented in global pathogen surveillance efforts. Extending this comparative framework to additional Indonesian provinces and isolate collections, and linking it prospectively to phenotypic drug susceptibility and treatment outcome data, would strengthen its value for both drug target discovery and routine programmatic MDR-TB surveillance, and would offer a scalable template for genomic and molecular surveillance of other emerging pathogens in comparable resource settings — an extension that this study does not directly test.

Limitations

This study had several limitations. Conservation and mutation frequency metrics were derived exclusively from short-read WGS data of previously published clinical isolates from a limited number of Indonesian regions (predominantly Java, with a smaller contribution from Papua) and may not be generalizable to other Indonesian islands or lineage backgrounds. Because this cohort was assembled from four independently designed public studies rather than a single prospective collection, differences in sampling frame, patient inclusion criteria, and metadata completeness across the contributing BioProjects are an additional source of heterogeneity that this secondary analysis cannot fully control for.

The absolute number of variable sites underlying the gene-level metrics is low. Four emerging genes carry fewer than eight variable sites across the entire cohort, and per-gene Ti/Tv ratios rest on between 3 and 44 qualifying biallelic positions. Panel-level comparisons, which aggregate across ten genes, are robust to this; gene-level comparisons within a panel are not, and we have accordingly refrained from ranking individual emerging targets. With only ten genes in each panel, the panel-level tests also have limited power to detect differences smaller than those observed, and the reported confidence intervals should be consulted alongside the p-values.

The NS:S substitution-count ratio used to compare selective histories is a simplified metric that does not correct for the number of possible synonymous/non-synonymous sites, codon usage, gene length, or the mutation spectrum, and is therefore not equivalent to a codon-normalized dN/dS estimate (26, 27). Ratios computed on small denominators are additionally unstable: the pncA ratio of 49.0 rests on a single synonymous variant.

Three limitations arise from the consensus-sequence approach. First, no minimum-depth mask was applied when generating consensus genomes, so positions with insufficient read coverage inherit the H37Rv reference base rather than being reported as missing; low-coverage regions will therefore appear conserved, and conservation percentages should be regarded as upper-bound estimates in this respect. Second, IUPAC ambiguity codes arising from heterozygous consensus calls were retained as character states in the conservation and entropy analyses, which acts in the opposite direction, and were excluded from the Ti/Tv analysis. Third, nine isolates were excluded for failing lineage assignment and carrying variant loads incompatible with a single-strain infection; these isolates accounted for the large majority of apparent sequence diversity in a preliminary analysis, and their exclusion, while necessary, means the reported diversity applies to unambiguously single-strain infections only. Mixed infections are epidemiologically relevant in high-burden settings, and a dedicated analysis of within-host diversity would require a different methodological approach.

The TBProfiler mutation library used here dates from November 2021 and predates the current WHO catalogue of M. tuberculosis mutations associated with drug resistance; resistance classifications therefore reflect the catalogue available at the time of analysis, and re-profiling against the current WHO catalogue may reclassify a small number of isolates. Finally, phenotypic drug susceptibility testing was not available for this cohort to cross-validate the genotypic resistance predictions, and high-frequency variant sites were defined by cross-sectional allele frequency alone, which does not distinguish recurrent independent mutation from clonal expansion of a lineage-associated polymorphism.

Conclusion

This study directly compares evolutionary conservation between classical and emerging anti-tubercular drug targets within a single, geographically diverse Indonesian Mycobacterium tuberculosis collection spanning Java and Papua. All twenty genes examined were highly conserved, consistent with the clonal population structure of the species, but classical resistance-associated loci carried four times the variable sites of therapeutically unexploited metabolic targets and were significantly less conserved on every metric examined. Their excess variation was disproportionately non-synonymous and concentrated in the prodrug-activating enzymes pncA and katG, where loss of enzyme function is the resistance mechanism — a selection signature entirely absent from the emerging panel. These results support the emerging metabolic pathways collectively as a low-variability target space for future anti-tubercular development, while showing that population-level conservation data from a clonal pathogen cannot, on their own, discriminate reliably between individual candidate genes within that space. Beyond target prioritization, this framework illustrates the value of whole-genome sequencing for genomic surveillance of tuberculosis in high-burden, resource-diverse settings. Because the present analysis examined only clinical Mtb isolates, its extension to a broader, One Health-oriented surveillance of other emerging pathogens and environmental or animal reservoirs is presented here as a possible future application rather than a demonstrated one.