Rifampicin has been a cornerstone drug in first-line tuberculosis treatment since the early 1970s, often used in conjunction with isoniazid, ethambutol, and pyrazinamide (Fox et al., 1999). Resistance mutations to rifampicin occur in the rpoB gene encoding the beta subunit of RNA polymerase III (Ovchinnikov et al., 1983; Jin and Gross, 1988; Telenti et al., 1993). These mutations are found in four distinct segments known as the rifampicin resistance determining region, which comprises the exit tunnel of the RNA-polymerase complex (Campbell et al., 2001; Molodtsov et al., 2017) where the drug binding occurs (Lin et al., 2017).

Many studies of clinically resistant isolates have revealed that the mutation most frequently encountered is S450L (codon 531 in Escherichia coli numbering) (Goldstein, 2014). This mutation is believed to have the least impact on bacterial fitness, and strains carrying it often exhibit extremely high rifampicin MIC. However, many other mutations are consistently found in clinical isolates (Walker et al., 2022). We decided to analyze the spectra of potential individual amino acid substitutions and reviewed publications starting from a comprehensive TBDreamDB database published in 2009 (Sandgren et al., 2009).

Before the era of whole genome sequencing, various methods were used to detect resistance-associated mutations, such as allele-specific or realtime PCR for already known changes in the genome (Pholwat et al., 2011) or Sanger sequencing for the identification of new ones (Ovchinnikov et al., 1983). We had to exclude data obtained by widely used strip hybridization technology by Hain Lifescience (Hillemann et al., 2007) due to partial information loss on specific mutations, detected by the absence of hybridization signal of the amplicon with wild-type probe (Miotto et al., 2018). Some studies using targeted Sanger sequencing lacked raw data and provided only amino acid substitutions, and thus were also excluded from codon substitution pattern analysis. The general data on the proportion of substitutions is shown in Figure 1. We also tried to exclude isolates with multiple mutations that occur in different places in the RpoB to avoid the effects of epistatic interactions, which would require a larger number of isolates to be properly evaluated. We analyzed mutations in both resistant and susceptible isolates because low-level resistance to rifampicin caused by substitutions in codons such as 435 was probably missed. In 2020, the critical rifampicin concentration was lowered from 1 mg/L to 0.5 mg/L for the Bactec MGIT 960 system, and therefore, isolates previously identified as susceptible could indeed have been resistant (Antimycobacterial Susceptibility Testing Group, 2022). Additionally, we were interested in analyzing all protein sequence variations that do not significantly affect the activity of RNA polymerase or the pathogen’s adaptation properties.

In total, our study included 11,730 isolates from 43 studies (Supplementary Table S1). More than half of the isolates had mutations in codon 450, followed by codons 445 and 435. Remarkably, a diverse range of amino acid substitution types was observed in all codons. Even at codon 450, where the S450L mutation dominated and was found in 67% of the cases, eight other amino acid substitutions were identified. Seven of these could not be achieved through single nucleotide substitutions, and the total proportion of 2N mutations for this codon was 2.5% (see Figure 1). The most diverse codon was 445, which codes for histidine in the wild-type isolate. Thirteen different amino acids were tolerated by RNA polymerase in this position, and 2N mutations were found in 7% of isolates with an altered 445 codon.

In general, it can be concluded that 2N substitutions provide amino acid diversity that cannot be achieved through 1N mutations (Belinky et al., 2019). However, it is worth noting that, in some cases, the same amino acid substitution could also result from a single nucleotide change, highlighting the inherent stochastic nature of evolution. The frequencies of alternative 2N are two orders of magnitude lower, as is evident from the analysis H445L and L452P mutations (see Figure 1).

In the specific context of drug-induced selective pressure, these substitutions decrease rifampicin binding to RNA polymerase without adversely affecting polymerase activity (Trinh et al., 2006). The wide range of substitution types observed can be attributed to the mechanism of drug action. Rifampicin binds to the DNA/RNA channel at a distance that exceeds 12 Å from the active site, effectively obstructing RNA elongation beyond 2–3 nucleotides (Campbell et al., 2001). It could be supposed that the variability would be lower for the drugs that bind directly to more conservative active sites of essential bacterial proteins.

Bedaquiline and linezolid were introduced into tuberculosis treatment 40 years after rifampicin, in 2014 and 2010, respectively. Since 2014, they have been used in the BPaL combined therapy scheme for drug-resistant TB around the world (Conradie et al., 2020). The safety and efficacy of this scheme led to the redefinition of tuberculosis with extended drug resistance, now defined as resistance to first-line drugs plus fluoroquinolones and bedaquiline or linezolid (Viney et al., 2021; Veziris et al., 2021). Although acquired drug resistance to these new drugs remains relatively low, there are a growing number of reports of clinical isolates with resistance (Timm et al., 2023).

The primary mutations selected during bedaquiline treatment emerge in the rv0678 gene, which encodes the repressor of the efflux operon (Andries et al., 2014). Loss-of-function mutations in this repressor have a minimal fitness cost and can be distributed throughout the open reading frame (Hartkoorn et al., 2014). Although most of them are frameshift mutations, amino acid substitutions in critical regions are also observed (Zimenkov et al., 2017). However, the huge number of possibilites to inactivate the gene by frameshifts limits the population frequency of complex and rare substitutions.

In contrast, substitutions in the bedaquiline target, the c-subunit of ATP synthase, which comprises the rotor part, are more conservative and rare, limited by the essentiality of the gene and, probably, by the fitness cost (Ghodousi et al., 2022). Mutations in the corrresponding atpE gene were initially identified in vitro and later confirmed in clinical strains of Mycobacterium intracellulare (Alexander et al., 2017) and M. tuberculosis (Zimenkov et al., 2017). To date, there are ten studies reporting 19 strains with substitutions in codons 25, 28, 61, 63, and 66 of atpE in clinical strains of patients who received bedaquiline in treatment schemes and developed decreased susceptibility or resistance (Zimenkov et al., 2017; Martinez et al., 2018; Peretokina et al., 2020; Mokrousov et al., 2021; Chesov et al., 2022; Ghodousi et al., 2022; Le Ray et al., 2022; Shang et al., 2023; Umpeleva et al., 2024; Peng et al., 2024). The WHO ‘Catalogue of mutations’ reported eleven isolates with non-synonymous mutations in atpE, but these sources of data overlap significantly (Walker et al., 2022). AtpE I66M and A63P are the most frequent substitutions found in six (Mokrousov et al., 2021; Chesov et al., 2022; Le Ray et al., 2022; Umpeleva et al., 2024) and four cases (Ghodousi et al., 2022; Chesov et al., 2022; Umpeleva et al., 2024; Peng et al., 2024), respectively. Additional substitutions in these codons, A63V and I66V were also described (Zimenkov et al., 2017; Martinez et al., 2018); however, all reported substitutions resulted from a single nucleotide change.

During the study of bedaquiline resistance in tuberculosis cases, we identified an isolate with two simultaneous mutations in atpE. This strain had two single nucleotide mutations at different positions, resulting in G25S and D28G amino acid substitutions (Zimenkov et al., 2017). The simultaneous presence of these two mutations in a single cell was confirmed by whole genome sequencing. We believe such occurrences are less frequent than single-nucleotide mutations, either it was a sequential or simultaneous mutational event.

However, the most intriguing case in the context of the hypothesis was the discovery of bedaquiline resistance caused by a rare dinucleotide mutation in codon 63 (GCA changed to TTA), leading to an A63L substitution (Ushtanit et al., 2021). No mutations in rv0678 or other proposed determinants were found in this isolate, yet the minimum inhibitory concentration (MIC) of bedaquiline was more than 2 mg/L (Bactec MGIT 960). A microevolution study of this case revealed a general pattern: the acquisition of the frameshift mutation t141tc in rv0678 initially occurred in a mixed state, followed by the extinction of the wild-type strain, and ultimately its elimination with simultaneous development of the AtpE substitution. A similar microevolutionary pathway has been proposed based on clinical studies and in vitro mutagenesis (Ismail et al., 2019).

Interestingly, some other mycobacterial species are intrinsically resistant to bedaquiline. It was suggested that the presence of methionine in position 63 of AtpE, compared to alanine in M. tuberculosis, is responsible for higher MICs in species such as Mycobacterium flavescens and Mycobacterium xenopi (Aguilar-Ayala et al., 2017; Petrella et al., 2006). Alanine 63 interacts directly with bedaquiline (Guo et al., 2021) and is substituted with proline or valine in resistant clinical isolates and in vitro selected strains (Peretokina et al., 2020; Andries et al., 2005). Both A63P and A63V are the result of 1N substitutions in codon GCA. Other residues, that could be obtained by 1N include E, G, S, and T. The wild-type A, also as P, V and G residues have similar biochemical properties being small, non-polar and hydrophobic. The absence of A63G substitution in clinical strains, resistant to bedaquiline could be explained by the steric properties of glycine, which breaks the alpha helix structure. Dinucleotide substitutions add R, D, Q, I, L, K to the diversity of available amino acids in this position. R, D, Q are hydrophilic, which are unlikely to present of the ATP synthase rotor part exposed to the membrane. Further, K is a positively charged residue, which could have a consequence during the rotation and contact with the stator. Thus, L and I remains the most probable candidates for the 2N substitution of A63 that preserve the function of the enzymatic complex. The substitution of A63 codon for M could be obtained by the change of all three nucleotides in codon (GCA>ATG).

A second example of a rare 2N mutation was observed during a clinical trial in a case of ineffective treatment, resulting in resistance to linezolid. Linezolid inhibits translation by interacting with the ribosome, and resistance mutations are found in distant regions of the 23S rRNA, forming the peptidyl-transferase center. However, most resistant isolates of M. tuberculosis carry the substitution C154R in the ribosomal protein RplC, which is typically caused by a single nucleotide change t460c. We identified a resistant isolate with mixed chromatogram data, and Sanger sequencing was unable to distinguish between the simultaneous presence of the nonstandard C154R (CGG codon) and C154S (AGT codon) or two different types of C154R replacements (Zimenkov et al., 2017). Further whole genome sequencing confirmed the latter, where one subpopulation had the usual t460c mutation, while the other had the TGT codon substituted with AGG by the 2N mutation (Figure 2). The simultaneous presence of two different subpopulations was confirmed by the analysis of subsequent isolates; Even after 15 months, a mixed chromatogram was identified in the clinical isolate obtained from the patient.

Both these cases illustrate the rapid emergence of resistance under the pressure of drug treatment. Similarly to resistance to rifampicin, the 2N mutation in atpE expands the diversity of the primary protein structure, while in rplC, such a mutation does not lead to novel amino acid substitution compared to the commonly found single nucleotide mutation t460c. Although the number of mutated isolates obtained during and after ineffective treatment is not sufficient to provide strong statistical power for correlation analysis (Nimmo et al., 2022), these cases of 2N mutations undeniably provide evidence of ongoing selection in these loci, and directly point to resistance mechanisms. Consequently, it could be proposed that mutations in these codons are clinically relevant and reflect the development of resistance upon the antibacterial treatment.

Although the described examples of bacteria adaptation to drug action unequivocally demonstrate the elevated level of 2N mutations in genetic determinants of resistance, it was interesting to perform a genome-wide analysis of dinucleotide mutation frequencies for the identification of all genes under selection pressure. We assumed that additional resistance determinants and host-pathogen interaction could be found by this approach. For the analysis, we used a data set consisting of 9,941 of 12,288 isolates that were sequenced and analyzed by the CRYpTIC Consortium (The CRyPTIC Consortium, 2022b; The CRyPTIC Consortium, 2022a). The remaining 2,347 had an ambiguous description of amino acid substitutions and were omitted.

Two main approaches were used for the analysis of mutation frequencies: allele counting and homoplasy counting (Chen and Shapiro, 2015). Although the former uses frequencies observed in the population and does not require intensive calculations of phylogeny and ancestor reconstruction of mutation events, it has significant disadvantages in its application to M. tuberculosis. The clonal nature of the pathogen led to complete linkage disequilibrium and population stratification, which could result in false associations obtained by uncorrected allele counting methods (The CRyPTIC Consortium, 2022b).

The homoplasy count derived from the number of independent events in different branches (Farhat et al., 2013) could be a signal for positive selection (Shapiro et al., 2009). Following this method, we constructed the phylogenetic tree using mutations from all genes, excluding highly repetitive PE/PPE, PE-PGRS genes and insertion elements (Supplementary Figure S1). The reliability of the phylogenetic tree was verified by analyzing the distribution and frequency of selected lineage-specific SNPs (Coll et al., 2014).

The maximal parsimony approach was applied to the set of all variations at a specific codon. All mutations were mapped on the tree and traced by joining neighbors with this mutation and mapping it to the common ancestor. It was supposed that reverse mutations could have occurred if it allowed the joining of branches with mutations split by smaller branches with the wild-type codon. Thus, the final number of all events was minimized for each codon across the whole tree.

As an example, the mutation E112K in the rv0078A gene (codon change GAG>AAG, nucleotide substitution c87468t), which could be used as a phylogenetic marker of lineage 2 (Coll et al., 2014), was present in 3,838 of 9,941 isolates. Indeed, they all belong to this lineage (n = 3,842 isolates) of the constructed M. tuberculosis phylogenetic tree; however, four isolates in L2 did not have this mutation. The phylogeny-adjusted frequency of independent events in this case for the entire population was equal to 23. The introduction of the maximal parsimony approach, based on the proposal that in these four isolates the reverse substitution occurred recently, allowed the decrease in homoplasy counts to two for direct and reverse mutations (four genomes with ‘wild-type’ Lys at 112 constitute one branch). Therefore, the direct mutation was traced to the common ancestor of lineage 2.

In such an approach, the consequent 1N mutations leading to final 2N substitution were treated as separate, thus underestimating the final diversity of codon substitutions and decreasing the numbers of 2N events. However, counting the second mutation as 2N will lead to an overestimation due to cases where the first mutation was specific to the branch, for example, the whole lineage 2. Thus, the analysis was limited to the calculation of simultaneous 2N substitutions, providing a more specific but less sensitive result in terms of codon diversity at particular position.

The resulting number of events for 275,531 codons in 3,680 genes analyzed were the following: 142,321 events of synonymous substitutions, 252,702 events of 1N substitutions and 1,645 events of 2N substitutions. The latter occurred in 862 codons in the genome, with a maximum value of 350 for the substitution of cyp138 P114F. Excluding the genes with single events of 2N substitutions results in 251 codons in 97 genes. The comparison of the number of codons with 2N substitutions and 2N events for each gene is shown in Figure 3. The two most important genes, that cause resistance, rpoB and katG, have elevated parameters and are clearly distinguishable from the entire set of genes.

Splitting of the counts of nonsynonymous mutations into 1N and 2N mutations allowed the introduction of two metrics, the ratio of dinucleotide to single-nucleotide nonsynonymous substitutions d2N/d1N, and ratio of dinucleotide to synonymous substitutions d2N/dS. In the commonly used GWAS approach, the comparison of dN/dS in resistant and susceptible sets of isolates is performed for identification of resistance determinants or genetic traits that co-evolve with resistance (The CRyPTIC Consortium, 2022b). We used another approach for identification of all genes under selection by the use of codon-based models of evolution and estimation of the difference between measured and expected frequencies of codon substitutions (Delport et al., 2009; Anisimova and Kosiol, 2009).

The codon substitution matrix for the Monte-Carlo modeling was derived from the whole dataset of phylogeny-adjusted events. Additional correction was introduced by proportional use of all variants of the particular codon in the whole set of isolates as the starting codon in modeling. Otherwise, lineage-specific codons with high prevalence in the genome set would not add to the calculated substitution frequencies. The difference between observed and estimated frequencies allows the identification of any genomic traits that are putatively under selection (Wilson and Consortium, 2020). The statistical approach for the identification of significant traits was analogous to that of the widely used ratio of non-synonymous to synonymous substitutions dN/dS (The CRyPTIC Consortium, 2022b). Therefore, three values for genes were analyzed and compared simultaneously presented as Manhattan plots (Figure 4), two for 2N substitution statistics and conventional dN/dS was used as reference. The significance threshold of -log(P value) was estimated by Bonferroni correction (Power et al., 2017).

The phylogeny adjusted counting of single-nucleotide mutations of different type and the comparison between genes resulted in distinct statistical signals for 68 genes, including the main resistant determinants rpoB, katG, rpsL, pncA, gid, embB, and ethA (Supplementary Table S2). Fourteen genes had a statistically significant change, either increase or decrease, in the dinucleotide events statistics, including the resistance determinants rpoB, katG, and pncA. Two genes, rv2024c and rv2082, with low p-values of d2N/d1N and d2N/dS, were not significant by dN/dS analysis (Table 1).

The statistically significant increase in 2N substitutions was observed for five genes only (Table 1). The rpoB gene had the most diverse mutation types and the largest number of isolates due to its role in resistance to rifampicin. Eighty-seven isolates with twelve 2N substitutions in six codons were found. Homoplasy counts of 2N mutations in the whole gene were 43, while that of 1N mutations was 2,216 (3%). Several types of dinucleotide substitutions were identified in codons 435, 445, and 450, similar to the pattern described in Section 2.1. An additional new 2N mutation, I542K, was also found.

In addition to rpoB, the main isoniazid resistance katG gene was found to be significant by dN/dS and d2N/dS ratios. The most prevalent S315T mutation found in resistant isolates is obtained by single nucleotide G>C substitution (codon AGC is changed to ACC). In this set of isolates, three other dinucleotide substitutions were identified that led to the same amino acid change: AGC codon was changed to ACA, ACG, or ACT.

The favored use of particular amino acids in drug-binding pockets with low fitness-cost could explain this observation for resistance genes rpoB and katG. Other genes associated with resistance to first- and second-line drugs in M. tuberculosis, such as rpsL, embB, gyrA, gid, and ethA, were significant only by dN/dS analysis.

Two more genes, that were not previously linked with resistance, had significant statistical signals for 2N substitutions, cyp138 and rv0988. Both amino acid substitutions, cyp138 P114F (CCT > TTT) and rv0988 L191A (CTG > GCG) cannot be obtained by single nucleotide substitution, confirming the importance of particular amino acid residues at these positions.

The predicted cytochrome P450 cyp138 is involved in pH tolerance mediated by MTS1338 sRNA regulation (Martini et al., 2023). Cytochromes P450s Cyp138 join metabolism, virulence, and susceptibility to antibiotics by altering cell wall composition (Lu et al., 2024). The gene rv0988 is a subunit of the ABC transporter (Rosas-Magallanes et al., 2006) and negatively regulated in the presence of toxic levels of CuCl2 (War et al., 2008).

Both cyp138(P114F) and rv0988(L191A) emerged in many unrelated strains of the M. tuberculosis lineage 2 (East-Asian) and were absent in all other lineages: lineage l (Indo-Oceanic), lineage 3 (East-African-Indian), lineage 4 (Euro-American). Cyp138 (P114F) is present in 3,154 of 3,527 isolates of sublineage L2.1. Such lineage-specific distribution differ from the distributions of resistance-associated variants, which are present in all human-adapted lineages–L1, L2, L3, L4, which is cause by the common therapy used in different regions of the world. Thus, it could be supposed that mutations in cyp138 and rv0988 are related to the pathogen adaptation to host pressures from immune system and antituberculosis drug use. Alternatively, these traits could be responsible for the evolutionary success of Lineage 2, which was widely discussed in the phylogeography studies on tuberculosis (Zhu et al., 2023).

Using the maximal parsimony approach, we suppose that P114F occurred only once in the common ancestor of 2.1, and then, in 373 isolates the reverse substitution F114P emerged (Figure 5). The phylogeny adjusted frequency of the reverse mutation is 349, resulting in 350 total events in the codon for this set of isolates. Alternatively, single-stage accumulation of P114F in the 3,154 genomes demands 1,130 mutational events.

The emergence of cyp138(P114F) was recently described as having evolved under the pressure of HIV-1 coinfection (Koch et al., 2017). However, taking into account the model of sequential forward and reverse mutations in Lineage 2 strains, the opposite situation could be supposed: the revertion F114P occured under the pressure of immunocompetent HIV-naïve hosts.

The same considerations are also relevant for rv0988(L191A): forward substitution in the common ancestor of a smaller sublineage within L2.1 and reversion to Leu in many unrelated strains, which is common for all other lineages (Figure 5). The phylogeny-adjusted frequency in the absence of reversion hypothesis was 742, compared to 241 events of reverse substitutions. Interestingly, the distribution of A191L revertion is not unified in L2 sublineages. For example, the Beijing BO/W148 genotype, which is a successive clone spreading in Russia (Mokrousov et al., 2012), has a low frequency of reversions, while some other branches are nearly completely reverted (Figure 5). Such a non-random purge of probably deleterious polymorpisms (Namouchi et al., 2012) demands further studies on correlation with geographical localization of clinical isolates and impact on virulence and pathogenicity.

Contrary to the previous, the substitutions N508T and C514L in rv2024c could be achieved also by simple 1N substitutions in corresponding codons, however, the codon usage in these cases is significantly lower compared to obtained by observed 2N substitutions. Thus, the observed AAT>ACG for N508T changes frequency from 5.3 to 35.1, while the codon Thr codon ACT has frequency of 3.7 (Naka et al., 2000). Similarly, for the C154L codon usage frequency increased from 8.1 to 50.4, while the alternative Leu CTA codon, which could be obtained by 1N substitution, has frequency of 4.8 (Naka et al., 2000). The exceptional use of 2N substitutions in unrelated strains also points to the preference for particular amino acid at these positions of the protein. Rv2024 has an unknown function but showed diversity in parallel and sequential isolates during the within-host microevolution analysis (Ley et al., 2019) and is believed to have a virulence-associated function, since it is was found to be induced in mouse lungs (Dubnau et al., 2005).

Two more classes of genes could be identified among those with a low number of 2N substitutions: under positive selection (dN >> dS) and under purifying selection (dS >> dN). A high number of non-synonymous substitutions were identified for pncA, cut, pstP, dipZ, and Rv2082 in the top-scored list of genes (Table 1). Mutations in pncA that lead to the loss of function of the gene are selected during pyrazinamide treatment and lead to resistance. Such mutations are frequent, distributed along the open reading frame, and include frameshifts. Therefore, the low number of 2N substitutions is not an unexpected phenomenon, also as for the rv0678c gene alterations leading to the resistance to bedaquiline. In addition to selection for loss of function, genes with high 1N and low 2N counts could encode proteins that are tolerant to substitutions, or the possible changes with low fitness cost are achieved by 1N codon substitutions.

Genes with a high number of synonymous substitutions have low counts of 1N and 2N substitutions, and the statistical significance of the d2N/dS ratio is the direct consequence of the significance of dN/dS ratio.