1. Introduction

Tuberculosis (TB) continues to be the most challenging, constantly present infectious disease worldwide, with 7.5 million newly reported cases and 1.3 million deaths per year(1). The resurgence of TB due to the SARS-CoV-2 pandemic(1) underscores the interconnected nature of global health and economic issues with TB incidence and control.

The causative agents of TB are members of the Mycobacterium tuberculosis (M. tuberculosis) complex. These obligate pathogen bacteria can incur a sustained threat to humanity thanks to their long term latency(2) and highly unresponsive nature to antibiotics(3, 4). Understanding the treatment evasion mechanisms and the outstanding stress tolerance of Mycobacteria are in the spotlight of TB research(5),(6). Drug tolerance arises when certain bacterial populations are temporarily able to survive antibiotic pressure in the absence of drug resistance-conferring mutations. Upon exposure to bactericidal drugs, tolerant mycobacteria are eliminated at a lower rate than the fully susceptible population(7). Several interconnected biological pathways are involved in the emergence and establishment of a drug-tolerant state(8, 9) including metabolic slowdown, metabolic shifting, cell wall thickening and transcriptional regulation guided adaptation(10). For example, several efflux pumps are upregulated under antibiotic stress(11, 12). In addition to temporary drug tolerance, the occurrence of genotypic resistance against the few useable antituberculotics is also recurrent(13). Interestingly, horizontal gene transfer, which is a major contributor to antibiotics resistance in other species does not appear to function in members of the M. tuberculosis complex(14, 15). Therefore, any resistant genotype can only emerge by de novo mutagenesis.

It is now commonly accepted that the M. tuberculosis population within individual TB patients can be more heterogeneous than was traditionally thought(16, 17). The coexistence of both drug-resistant and drug-sensitive strains in a single patient, or even several drug-resistant strains with discrete drug resistance-conferring mutations has been described in clinical isolates(18–20). Warren et al. found that the occurrence of mixed infections reached 19 % of the examined patients in South Africa by using a PCR-based strain classification method(21). Mixed infections can result from i) simultaneously or sequentially acquired infections by different strains or ii) genomic evolution of a strain under mutagenic pressure within the host (termed microevolution) and consequent coexistence of several populations. Accordingly, the emergence of genetically encoded resistance may either be due to microevolution or to the spreading of already existing variants from polyclonal infections under drug pressure. The difference between these two underlying mechanisms for the emergence of drug resistance is highly relevant to the treatment of TB. The investigation of stress induced mutagenesis in Mycobacteria are based on fluctuation assays(22),(23) besides several indirect evidence from descriptive studies(24–28). Some studies demonstrate the simultaneous presence of several subpopulations within the same host which they interpret as an indication of being prone to microevolution(19, 29). It is also possible that certain strains have intrinsically higher mutability. For example, the lineage 2 strains of the Beijing genotype were shown to exhibit a higher mutation rate(30). On the other hand, others found stable M. tuberculosis genomes with no or only few emerging genomic changes over prolonged periods of treatment(26). Genotyping has enabled researchers to describe cases of co-infection by ≥2 different strains (mixed infection) or coexistence of clonal variants of the same strain(29, 31, 32). Introducing whole genome sequencing into this field still leaves the distinction between mixed infections with multiple similar strains and strains arisen by microevolution elusive. Depending on the elapsed time between two sample collections, the stepwise acquisition of mutations might be missed and the observed diversity may reflect concurrently existing subclones rather than newly emerged mutations(27). In addition, a single sputum sample usually does not represent the whole genomic diversity of the infection(16, 31). Cell culturing can also lead to additional artefacts(33, 34). The lack of standardized reporting of genome sequencing analyses also limits our ability to draw conclusions on within-host microevolution(27). Therefore, although several factors such as drug pressure and disease severity have been suggested to drive within-host microevolution and diversity(35, 36) and it is now accepted that the M. tuberculosis population within individual patients can be heterogeneous, we could not find any unequivocal proof for explaining the mechanism of emergence of the observed genomic diversity which gives rise to drug resistance.

Therefore, to advance our knowledge on the effect of antibiotics on mycobacterial mutability, we conducted experiments under controlled laboratory conditions. We used Mycobacterium smegmatis (M. smegmatis) for our investigations. This non-pathogenic relative of the medically relevant Mycobacterium species shares most DNA metabolic pathways with the medically relevant strains. Davis and Forse compared the sequences of proteins involved in base excision repair and nucleotide excision repair pathways in E. coli and their homologs in M. smegmatis and M. tuberculosis and found that there is a high degree of conservation between the DNA repair enzymes in M. smegmatis and M. tuberculosis(37, 38). Bioinformatic analyses of completely sequenced mycobacterial genomes including M. tuberculosis(39), M. leprae(40), M. bovis(41, 42), M. avium, M. paratuberculosis, and M. smegmatis(43) also demonstrated through the comparison of genes participating in many of the DNA repair/recombination pathways that the basic strategy used to repair DNA lesions is conserved(44, 45). Durbach et.al, investigated mycobacterial SOS response and showed that the M. tuberculosis, M. smegmatis and M. leprae LexA proteins are functionally conserved at the level of DNA binding(46). In our earlier paper, we also compared the enzymes of thymidylate biosynthesis in M. tuberculosis and M. smegmatis and found high conservation (47). Considering M. smegmatis is non-pathogenic and fast-growing, it provides an attainable model to obtain information on genomic changes under drug pressure in M. tuberculosis.

We systematically investigated the effects of currently used TB drugs on genome stability, on tolerance/ resistance acquisition, on the activation of the DNA repair system, and on the cellular dNTP pool. We focused particularly on drugs used in the standard treatment of drug-susceptible TB, comprising isoniazid (INH), rifampicin (RIF), ethambutol (EMB), and pyrazinamide (PZA), the so-called first-line antibiotics(48). We also used a second-line antibiotic, ciprofloxacin (CIP). We found that following exposure to these antibiotics, the activation of DNA repair pathways maintains genomic integrity, while non-genetic factors convey quick adaptation for the stress conditions. Notably, even with prolonged antibiotic exposure exceeding 360 bacterial generations, we observed no significant increase in the mutation rate, suggesting the absence of de novo adaptive mutagenesis.

2. Material and methods

2.1. Bacterial strains, media and growth conditions

M. smegmatis mc²155 (49) strains were grown in Lemco broth (5 g/l Lab-Lemco, 5 g/l NaCl, 10 g/l Bacto pepcone, 0.05 % Tween-80) or on solid Lemco plates (6.25 g/l Lab-Lemco, 6.25 g/l NaCl, 12.5 g/l Bacto pepcone, 18.75 g/l Bacto agar).

2.2 Optimization of stress treatments conditions

The optimal concentrations of drugs were searched for by using serial dilutions of the compounds on a logarithmic scale and by monitoring the number of colonies forming units (CFU) or the optical density (OD) at 600 nm on plates and in liquid cultures, respectively (Figure S1). We also monitored cell morphology in response to drug treatment. For further experiments, sublethal concentrations of drugs were used to obtain adequate quantity of research material (DNA, RNA, dNTP) for downstream analysis while the effect of the treatment was clearly indicated by a decrease of viability or change in cell size and morphology. The concentrations of applied drugs and stress conditions are compiled in Table 1.

Table 1

2.3. Stress treatment of liquid cultures

Cells were grown in 100 ml liquid culture until an OD (600) = 0.1±0.02 was reached, then the appropriate quantity of drug (Table 1) was added to half of the cultures. The other half of the same culture was used as a control. We conducted the treatments for 8 hours. The cultures were then centrifuged (20 min, 3220 g, 4°C) and the resulting pellets were used for downstream analysis. The total CFUs were determined for each culture. The generation time after the treatments was calculated using the formula:

where T_d is the generation time, t is the time interval between measurements, and N_t and N₀ are the final and initial population size, respectively.

2.4. Microscopic analysis of cell morphology upon treatments

For morphological studies, 200-200 µl stress-treated and control cells were retrieved before RNA or dNTP extraction and washed with PBS containing 0.1% Triton X-100. The cells were then fixed in 4% PFA dissolved in PBS for 30 min at 37°C. Cells were stained with 10 µg/ml DAPI for 30 min at 37°C, then streaked onto microscopy slides covered with 0.1% low melting agarose (Sigma). Imaging was done using phase-contrast and fluorescent modes on a Leica DM IL LED (Leica) microscope. The cell size and volume were quantified using the automated recognition of the BacStalk software (https://drescherlab.org/data/bacstalk/)(50). The cell volume distribution diagram was prepared using OriginPro 2018 (OriginLab Corporation, Northampton, MA, USA.). The sample size, calculated means and standard deviations are compiled in Table S1.

2.5. Mutation accumulation (MA) experiments

Sixteen independent M. smegmatis mc² 155 MA lines were initiated from a single colony for every treatment. The ancestor cell colony was generated by streaking a new single colony from plate to plate 5 times before the beginning of the treatments to ensure a single common ancestor. Lemco agar medium was used for the mutation-accumulation line transfers. The specific stress treatment conditions are summarized in Table 1. All MA lines were incubated at 37°C. Every 3 days, a single isolated colony from each MA line was transferred by streaking to a new plate, ensuring that each line regularly passed through a single-cell bottleneck(51). Treatments were performed for 60 days, therefore each line passed through ∼360-480 cell divisions. We calculated 3 hours generation time for mock treatments and 4 hours for stress treated lines. Some of the mock treatments were performed for 120 days to ensure a presumably sufficient number of mutational events without stress treatment. Following the MA procedure, a single colony was transferred from all strains to a new plate without stress-treatment and grew for another 3 days for expansion. Frozen stocks of all lineages were prepared in 20% glycerol at −80°C.

2.6. DNA extraction

A single colony at the end of the MA experiment was inoculated into10 ml liquid culture from all lineages, was grown until OD₆₀₀ = 0.8-1.0, and harvested. For genomic DNA purification, 5 or 6 grown cultures of individual lineages from the same treatment with identical estimated cell number (based on OD measurements) were pooled before isolation. For cell disruption, the cells were resuspended in 1 ml of 10 mM Tris, pH 7.5, and 0.1 mm glass beads were added to a final volume of 1.5 ml. The cells were disrupted using a cell disruptor (Scientific Industries SI-DD38 Digital Disruptor Genie Cell Disruptor) in a cold room (at 4 °C). After centrifugation for 10 min at 3220 g, and at room temperature, DNA was extracted from the supernatant by phenol:chloroform:IAA (25:24:1) extraction followed by isopropanol precipitation. The quality and quantity of the extracted DNA was evaluated using UV photometry in a Nanodrop-2000 instrument and by agarose gel electrophoresis.

2.7. DNA library preparation and whole genome sequencing

The DNA library preparation and whole genome sequencing (WGS) was done at Novogene Ltd., Beijing, China. Sequencing was executed on Illumina 1.9 instruments with 600-basepair (bp) fragments as 2 × 150 bp paired end sequencing. An average read depth of 267 was achieved across all samples.

2.8. WGS analysis and mutation identification

Three parallel pooled samples were sequenced for every treatment, each contained 5 or 6 individually treated MA lineages that add up to a subtotal of 15-18 individual lineages. FastQC was used for analysing the quality of the raw reads. In case if adapters and low-quality bases (Phred score < 20) were present in the samples, bases were trimmed with Trimmomatic(52). We mapped our paired end reads to M. smegmatis mc² 155 reference genome (GenBank accession number: NC_008596.1) by Bowtie2(53). PCR duplicates were removed with the use of Samblaster(54). We converted SAM files to BAM files, and sorted them with SAMtools(55).

Read groups were replaced by the Picard tool. Single nucleotide variations (SNVs), insertions and deletions were called from each alignment file using the HaplotypeCaller function of the Genome Analysis Toolkit(55). We analyzed the frequency of occurrence (% of all reads of a pooled sample) of each SNV, insertions and deletions (hits) with our in-house Python scripts and compared to the frequency of occurrence of the same hits in every other lineage. We considered mutations as spontaneously generated mutations only in case if no other lineages carried that variant in any depth and if hits reached at least 6% frequency of the reads at the corresponding position (theoretically, a spontaneously generated mutation in a pooled sample emerge with 20 % or 16.7 % frequency when 5 or 6 lineages are pooled, respectively, however, we allowed some variety when choosing 6 % as a lower limit and 39.9 % as an upper limit). Sequencing data are available at European Nucleotide Archive (ENA) with PRJEB71590 project number. Please note that we incorporated some of our additional sequencing data into the analysis, which is curated under the umbrella project at the ENA along with the present dataset.

2.9. RNA isolation and cDNA synthesis

For RNA extraction, cell pellets were resuspended in 2 ml RNA protect bacteria reagent (Quiagen; cat. no.:76506), incubated for 5 min at room temperature and centrifuged for 20 min at 3220 g and at 4 °C before storage at −80 °C. Total RNA extraction was performed with the Quiagen RNeasy Mini kit (cat. no.: 74524). To disrupt cells, 5 x 1 min of vortexing with glass beads in the manufacturer’s lysis buffer was performed followed by 1 min poses on ice. DNase digestion was performed on column with Quiagen DNase I (cat. no.: 79254), for 90 min at room temperature. For quantitative and qualitative RNA analysis, spectrometry by Nanodrop 2000 and non-denaturing 1 % agarose gel electrophoresis (50 min/100 V) were performed, respectively. cDNA synthesis was performed using the Applied Biosystems™ High-Capacity cDNA Reverse Transcription Kit with RNase Inhibitor (cat. no.: 4374967). 95-105 ng total RNA was used for each reaction.

2.10. Choosing the reference genes for the study

We tested SigA (MSMEG_2758), Ffh (MSMEG_2430), and ProC (MSMEG_0943) as possible reference gene candidates. SigA is a widely used reference gene in prokaryotes(56–58) Ffh and ProC genes are shown to be stably expressed in other pathogenes(59). Using GeNorm(60, 61) analysis, SigA and Ffh proved to be stably expressed in our experimental system (Figure S2).

2.11. Gene expression quantification

qPCR measurements were performed on a Bio-Rad CFX96 Touch™ Real-Time PCR Detection System. Primers were designed using IDT DNA oligo customizer (https://eu.idtdna.com/), and were produced by Sigma Aldrich (for sequences, see Table S2). The qPCR reaction mixtures contained 7-7 nmoles of forward and reverse primers, 0.25 µl of the cDNA, Bioline Mytaq PCR premix (cat. no.: 25046), and VWR EvaGreen (cat. no.: #31000) in a total reaction volume of 10 µl. The thermal profile was as follows: 95 °C/10 min, 50x (95 °C/10 s; 62 °C/10 s; 72 °C/10 s). Melting curves were registered between 55 °C and 95 °C with an increment of 0.5 °C (Figure S3). The applied primers and their measured efficiencies are compiled in Table S2. The qPCR data were analyzed using the Bio-Rad CFX Maestro software and numerically shown in Table S3. Non-reverse transcribed controls and no-template controls were used to account for any irrelevant DNA contamination. 3 technical, and 3 biological replicates were used for all measurement. All raw data can be found in the Supplementary ZIP file.

2.12. dNTP extraction

dNTP extraction and measurement were performed according to Szabo et al.(62). Briefly, the cell pellets were extracted in precooled 0.5 ml 60% methanol overnight at −20°C. After 5 minutes boiling at 95°C, the cell debris was removed by centrifugation (20 min, 13 400 g, 4°C). The methanolic supernatant containing the soluble dNTP fraction was vacuum-dried (Eppendorf) at 45°C. Extracted dNTPs were dissolved in 50 μl nuclease-free water and stored at −20°C until use.

2.13. Determination of the cellular dNTP pool size

Determination of the dNTP pool size in each extract was as follows: 10 pmol template oligo (Sigma), 10 pmol probe (IDT) and 10 pmol NDP1 primer (Sigma) (see sequences in Table S4) was present per 25 μl reaction. The concentration of each non-specific dNTP was kept at 100 μM. VWR® TEMPase Hot Start DNA Polymerase (VWR) was used at 0.9 unit / reaction in the presence of 2.5 mM MgCl₂. To record calibration curves, the reaction was supplied with 0–12 pmol specific dNTP. Fluorescence was recorded at every 13 seconds in a Bio-Rad CFX96 Touch™ Real-Time PCR Detection System or in a QuantStudio 1 qPCR instrument. The thermal profile was as follows: 95°C 15 min, (60°C 13 s) × 260 cycle for dATP measurement, and 95°C 15 min, (55°C 13 s) × 260 cycle for dTTP, dCTP and dGTP measurements. Results were analyzed using the nucleoTIDY software (http://nucleotidy.enzim.ttk.mta.hu/)(62). Results were given in molar concentrations for better comparison. To this end, cell volumes were calculated using the BacStalk software based on microscopic images for every treatment. Besides the graphical presentation of the result, numerical data can be found in Table S5.

2.14. Tolerance assay

We used a modified version of fluctuation assays(63) for the estimation of the rate of emergence of resistant cells upon preincubation with a sublethal dose of CIP. An initial 100 ml culture was grown to OD = 0.4-0.5, was centrifuged for 30 min at 800 g and at 4 °C, then resuspended in 5 ml Lemco. 100 μL from this stock solution was streaked and cultured on normal Bacto Agar plate, and Bacto Agar containing 0.3 μg/mL CIP at a standardized fashion. Parallel plates were incubated for 4, 24, 48, 72, and 96 hours at 37 °C. Colonies were washed off the plate with 6 mL Lemco broth by incubation for 30 mins on a rocking shaker. Then CFU was determined on selective (containing 0.5 μg/mL CIP), and non-selective Bacto agar plates. Tolerance rates were calculated using the following formula:

t: time of preincubation

: Number of colonies on non-selecting agar plates at reference time point

: Number of colonies on non-selecting agar plates at t hours

: Number of colonies on plates containing 0.5 μg/ml CIP at reference time point

: Number of colonies on plates containing 0.5 μg/ml CIP at t hours

− : Number of newly emerging resistant colonies

−: Number of growing colonies

: Generation time during pre-treatment

2.15. Statistics

We used an initial F-test to test the equality of variances of the tested groups. If the F-test hypothesis was accepted (p < 0.05), we used the two-way homoscedastic t-probe; if rejected, we used the two-way Welch’s t-probe to assess differences at a significance level p < 0.05 in each case. F- and t-statistics were counted for the ΔCt values(64) of the qPCR results and for the concentrations normalized to the cell volume in case of the dNTP measurements.

3. Results

3.1. Adapting stress conditions and assessing their impact on cell viability and morphology

For an efficient TB treatment, first-line antituberculotics are used in combination in the clinics (isoniazid – INH; ethambuthol – EMB; rifampicine – RIF; pyrazinamide - PZA)(36). To model this drug pressure in our study, we also combined the four first-line drugs in addition to applying them one by one. We added a second line antibiotic, ciprofloxacin (CIP). MitomycinC (MMC) and ultraviolet (UV) irradiation were used as positive controls for direct DNA damage(65, 66). We optimized the drug concentration for all applied treatments. First- and second line antituberculosis drugs were used in sublethal concentrations to convey a phenotypic effect while allowing to keep an adequate number of cells for the downstream measurements (Figure S1 and Table 1). In the first-line combination treatment, a 10-fold reduced concentration of each separately adjusted drug had to be applied to allow the survival of enough cells for the analyses (Table 1). The fact that a lower dose of antibiotics applied in combination resulted in higher CFU reduction indicates the synergistic effect of the first-line antibiotics on M. smegmatis growth inhibition (Table 1). After an eight-hour drug treatment, we determined the viable cell count by CFU measurements (Table 1). The bacteriostatic drugs INH and EMB caused moderate CFU decrease compared to the control (Table 1, Figure S1) consistent with their mechanism of action(67). To quantify the phenotypic effect of the applied drug treatments on the cells, we analyzed the cellular dimensions using microscopy (Figure 1 and Table 1). The observed morphological changes provided evidence of the treatments’ effectiveness (Figure 1 and Table 1). Specifically, following RIF, CIP, and MMC treatments, we observed cells elongating by more than twofold, whereas INH and EMB treatments led to a reduction in cell length (Figure 1, Table 1). The combination treatment did not affect the cell size.

Figure 1

We also assayed the clinically relevant drug pyrazinamide (PZA). However, M. smegmatis was reported to exhibit an intrinsic resistance to PZA(68). Indeed, PZA treatment on its own even at high concentration in acidic conditions had no effect on cell viability in our experiments (Figure S1). Therefore, we only used it in the combination treatment.

3.2. The genome of M. smegmatis remains stable even under antibiotic pressure

16 independent M. smegmatis MC² 155 lineages for each stress treatment condition and 56 lineages for the mock control were initiated and cultured from single colonies. The stress treated lines and some of the mock lines were maintained through 60 days on agar plates. The rest of the mock lines were maintained through 120 days on agar plates. Drug treated lineages were maintained for shorter times as more mutations were expected to arise under drug pressure. We estimated an average generation time of 4h for stress treated cells based on our experiments done in liquid media with CFU counting at the beginning and at the end of the treatments (for the calculation, see Materials and Methods). Therefore, bacteria produced ∼360-400 generations during the 60 days treatment. The untreated control (mock) cells were assessed to have an approximate generation time of 3 hours, resulting in around 480 or 960 cell divisions over 60 or 120 experimental days, respectively. Following the treatment on solid plates, we expanded each lineage in a liquid culture without drug pressure and isolated genomic DNA. All lineages were sent to WGS to reveal the mutational events induced by the drug treatments. We set conditions so to obtain at least 30-60x sequencing depth for all positions per independent lineage. The ancestor colony was also sent for sequencing to detect already existing variations compared to the reference genome. According to the WGS results, our M. smegmatis ancestor strain carried 151 various mutations compared to the M. smegmatis reference genome deposited in the GenBank. These mutation hits were also found in all treated and untreated lineages and were omitted from further calculations as these are specific variations of our laboratory strain. We also removed those mutation hits that were found in any other independent lineage at the same position in any depth.

After careful cross-checking of the sequencing data, a surprisingly low number of new mutations was detected. We found that a maximum of 1 mutation per lineage occurred during the 60-day drug treatments. Also, a maximum of 1 mutation per lineage was detected during the 60- or 120-day mock treatment (16 newly generated mutations for 56 lineages). Therefore, we calculated an approximately 5 x 10^-11 mutation rate for our untreated M. smegmatis mc²155 strain. To our great surprise, the mutation rate of all treated lineages was in the same range (2.5-20 x 10^-11) except for the UV treatment used for positive control (Figure 2 and Table 2).

Figure 2

Table 2

We analyzed each mutation except for those obtained following the UV treatment and found no sign of adaptive changes. The excel file containing the positions of the mutations in the genome and their identifiable consequences at the transcriptional or protein level is provided in the archive deposited with the article.

3.3. The DNA repair system shows a treatment-specific activation pattern

To reveal the reason for detecting so few newly generated mutations under antibiotic pressure, we investigated whether the DNA repair pathways and other elements of the stress response potentially involved(69) were activated under drug pressure. The mycobacterial DNA repair system is highly redundant, many of its enzymes have overlapping functions(45, 70, 71). Although canonical mismatch repair proteins are thought to be missing, a recently described protein, NucS is encoded with similar function(72). We investigated the expression pattern of DNA repair genes in all known DNA repair pathways in Mycobacteria including NucS using RT-qPCR, a method suitable to accurately show changes in transcript levels. The measured relative expression levels are shown in Figure 3 and Table S3.

Figure 3

Treatments with the two antibiotics affecting cell wall synthesis (INH and EMB) show similar patterns in the expression levels with an overall downregulation of DNA repair genes. On the contrary, CIP and MMC, drugs directly targeting DNA integrity induce a pattern marked by a moderate to strong overexpression of nucleotide excision and double-strand break (DSB) repair genes, respectively (Figure 3). DNA polymerases DinB2 and DnaE2 involved in these DNA repair pathways are also strongly overexpressed (Figure 3). RIF, the DNA-dependent RNA polymerase inhibitor does not seem to induce any change in the expression pattern of the investigated genes except for the Ahp peroxiredoxin (Figure 3). As a result of the 1^st line combination (COMBO) treatment, 14 out of 38 investigated genes are significantly (p <0.05) upregulated. More than 4-fold upregulation can be measured for 5 members of the base excision repair pathway. In addition, the MutT2 dNTP pool sanitization enzyme and the error-prone DNA polymerases are also strongly upregulated. (Figure 3). Interestingly, however, the DSB repair enzymes are exempt from this overall upregulation tendency (Figure 3). The strongest measured effect of all is the 17-fold expression increase of the KatG1 peroxidase (Figure 3). When the 1^st line antibiotics were used one by one, significant expression change could only be observed upon the INH treatment (4/38 genes) and in the opposite direction (downregulation).

3.4. All but the combination treatment alters the size and balance of dNTP pools

It was shown that dNTP pools are crucial for genome maintenance and proper DNA synthesis(73–76). Imbalanced or altered levels of dNTPs could cause increased rate of DNA lesions and therefore, may play a role in the development of drug resistance. Therefore, we measured cellular dNTP concentrations and ratios in function of the applied drug treatments using a fluorescent detection based method optimized in our lab(62). We used MMC treatment as a positive control as this is a generally used positive control for DNA damage(66, 77). To calculate cellular concentrations, we used the previously determined cellular volumes (Figure 1). Interestingly, we found altered dNTP pools upon most treatments (Figure 4 and Tables S4-5). The CIP treatment resulted in the most remarkable differences in particular for dATP and dTTP concentrations which increased ∼7-fold accompanied by a decrease in the dGTP concentration (Figure 4E). RIF, CIP and MMC treatments promoted an increase in the dCTP and dTTP pools (Figure 4D-F), while INH and EMB coincided with a decreased concentration of these nucleotides (Figure 4B-C). In the case of the combination treatment, the levels of these two nucleotides did not change (Figure 4A). The dGTP pool decreased in both absolute and relative values in all treatments (Figure 4A-G). Another major observation is that a larger overall dNTP pool size coincides with a larger cell size and vice versa (Figure 4H).

Figure 4

3.5. Stress induced drug tolerance is developed upon pretreatment with sublethal concentration of CIP

To compare the result of the mutation accumulation experiment to a phenotype-based drug resistance assay, we chose the fluctuation assay generally used in the literature(63). Mutation rates in these tests are calculated based on the difference in the number of CFU values between cultures grown in regular broth compared to those in selecting broths. These assays assume that the resistance exclusively occurs upon one mutation event. Since the genetic background of a drug-tolerant colony is not confirmed, this presumption potentially leads to a significant misinterpretation of the actual mutation rate. For clarity, we refer to the mutation rate estimations in our phenotype-based resistance assay as tolerance rate. For valid comparison with the results of our mutation accumulation assay, we installed similar experimental conditions. Specifically, culturing was done on agar plates, the applied drug concentrations were in the same range as used during the mutation accumulation process, then colonies were washed off and CFU counting plates were streaked from the resuspended bacteria (Figure 5A). We found that the estimated rate of emergence of the tolerance for CIP is three orders of magnitude higher than the mutation rate calculated based on WGS (10^-7 vs. 10^-10, cf. Table 2). Furthermore, following a 20-96 h exposure to a sublethal 0.3 μg/ml dose of CIP, a phenotypic tolerance appears in a significant portion of the cells to an otherwise lethal 0.5 μg/ml dose (Figure 5B). The tolerant cell population increased with the length of the preincubation time before reaching a maximum (Figure 5B).

Figure 5

4. Discussion

One of the reasons that TB is still a great medical challenge is the frequent incidence of resistant cases. The main goal of our research was to get a better understanding of the molecular mechanisms of drug resistance development in mycobacteria. We started from the hypothesis that the long-term exposure to first-line antitubercular drugs increases mutability.

4.1 Drug resistance in M. smegmatis does not arise from increased mutation rates under antibiotic pressure

Measured and estimated mycobacterial mutation rates in the earlier literature are in the order of 10^-10 /bp/generation(78),(79),(30),(80). This low constitutive mutation rate by itself does not explain the biological diversity observed in clinical isolates(24). This diversity might result from an elevated mutagenesis rate or the accumulation of different strains from the environment. We conducted a modelling study in M. smegmatis to investigate whether exposure to first-line antibiotics generates such biological diversity and if yes, by what molecular mechanism. We measured the appearance of drug-induced mutations in the genome in a mutation accumulation assay using whole genome sequencing. We also examined the rapid occurrence of phenotypic tolerance. The difference between the results of the phenotypic and the mutation accumulation studies was surprisingly large. Even without a pretreatment, a tolerance rate of 10^-7 /generation could be observed for CIP (consistently with literature data from fluctuation assays(81, 82)), while in the mutation accumulation assay, the number of mutations remained under 10^-10 /generation. In fact, the mutation rate increase was only significant in the case of the UV treatment serving as positive control for the experiments (Figure 2). Previous studies claiming mutation rate increase upon antibiotics treatment assessed mutation rates using fluctuation assays and no direct evidence of the change in the genetic material was shown(23, 83). Our findings imply that the emergence of drug resistance in this study is attributed to phenotypic factors. Phenotypic changes upon antibiotics treatment have widely been investigated(84) including potentially bistability(85) and/or the upregulation of efflux pumps(86, 87). It is noteworthy that spontaneous mutagenesis is easily induced through UV treatment. Considering that mycobacterial species spread through air droplets, it is conceivable that the exposure of these droplets to environmental UV radiation could potentially lead to the generation of new mutations.

4.2 The combination treatment with frontline drugs induces an overall upregulation in the DNA repair pathways aimed at eliminating misincorporations

The intracellular lifestyle of the TB pathogen implies that these bacteria must face various stress conditions and damaging agents including reactive oxygen and nitrogen species inside macrophages. Therefore, stress-induced transcriptional changes in Mycobacteria have been studied on genome-wide scales(84),(88) and one study found a specific activation of the DNA repair system in response to CIP similar to ours(66). Although M. smegmatis is not an intracellular pathogen, it shares the DNA repair pathways with M. tuberculosis and is often used to study how mycobacteria deal with DNA lesions(45). We focused our investigation on stress-induced transcriptional changes that may account for the protection of genomic integrity under the drug pressure of first-line antituberculotic drugs.

Redox potential change is a well-known and common phenotypic response to INH in Mycobacteria(89). The downregulation of KatG1 and Nei2 in response to our INH treatment (Figure 3) is in line with this and might indicate a reduced cellular redox potential. KatG1 is the enzyme that activates the prodrug INH(89), therefore, the downregulation of this enzyme decreases the active drug concentration and increases the tolerance of M. smegmatis against INH. In the case of the first-line combination treatment, however, KatG1 was highly upregulated, indicating high ROS levels in the cell(90). High ROS levels are known to cause damage to nucleobases and the nucleotide pool is a major effector of oxidative stress-induced genotoxic damage(91). In line with this, we observed upregulation in dNTP pool sanitation, base- and nucleotide-repair pathways which play crucial roles in preventing and repairing DNA damage caused by oxidative stress. The observed synergistic effect clearly results from the combination of first-line drugs, as we did not observe this effect when applying the drugs individually. The observed upregulation of the relevant DNA repair enzymes might account for the low mutation rate even under drug pressure. Notably, error-prone polymerases DinB2 and DnaE2 exhibited significant upregulation without inducing a mutator phenotype. This indicates that error-prone and error-free repair mechanisms are coactivated, predominantly resulting in error-free repairs.

4.3 dNTP pool alterations induced by frontline drugs neutralize each other in the combination treatment resulting normal DNA precursor pools

The building blocks of DNA constitute a critical component within the molecular aspects of mutability. It has been shown that increased or imbalanced dNTP pools induce mutagenesis in prokaryotes(92) and eukaryotes⁸⁶. To assess the impact of drug treatment on dNTP pools and its correlation with genome stability, we quantified the concentrations of dNTPs in cell extracts obtained from the drug-treated cells. When treating the cells with frontline drugs EMB and INH individually, the observed reductions in dNTP pool sizes and cell size (as illustrated in Figure 4E) aligned well with the concurrent downregulated transcript levels (Figure 3). Resting states of bacteria have also been characterized with a decrease in cell size and dATP levels(93, 94). These observations reflect the bacteriostatic effect of these drugs causing metabolic processes entering a dormant state, accompanied by the downregulation of enzymes involved in dNTP synthesis. This contrasts with the monotreatment with RIF, in which case dNTP concentrations, mainly dCTP and dTTP were elevated correlating with increased cell length. The combined treatment yielded the least significant alteration as if the opposite effects of the frontline drugs cancelled each other out. An elevation in the dNTP pools during cytostatic or cytotoxic treatment is more unexpected and suggests elevated DNA repair activity. This observation, particularly in the case of CIP treatment, aligns with the substantial increase in the expression of DNA repair synthesis genes, as depicted in Figure 3. Among all administered treatments, only the CIP treatment led to a notable dNTP imbalance and a substantial overall rise in dNTP pools, due to elevated levels of dTTP and dATP. This coincides with the largest changes in the expression of DNA repair genes, particularly those associated with the SOS response and homologous recombination (Figure 3). Interestingly, the dGTP level decreased with all drug treatments. This finding suggests that dGTP may play a role in a general stress response. It is noteworthy that not all dNTP imbalances are created equal. Specifically, an excess of dGTP has been identified as a significant contributor to mutations(95, 96). It must be noted that in these (and most) organisms dGTP is the least abundant among dNTPs. However, in mycobacteria, a unique scenario exists where dGTP is the most abundant dNTP species(97) and mycobacterial genomes are characterized by a high GC content(98). A reduction in dGTP levels in this context may contribute to minimizing DNA lesions by enhancing proofreading efficiency.

4.4 Our results do not support drug resistance acquisition through drug-induced microevolution

Our hypothesis that systematic antibiotics treatment induces mutation rate increase in M. smegmatis failed, as we did not observe any significant impact of antibiotics on mutability in laboratory conditions. Only in the case of CIP treatment, a second-line TB drug known for directly inducing DNA damage, could we detect a slightly elevated mutation rate. The treatment of M. smegmatis with the clinically used combination therapy drugs did not induce a mutator effect, quite the opposite. The observed activation of DNA repair processes likely mitigates mutation pressure, ensuring genome stability.

If there is no indication for a priori drug resistance, TB patients are treated with the combination therapy of first-line antituberculotics. In at least 17 % of the treatments, resistance to RIF or RIF+INH (called multidrug resistance) emerges(1). There are two models for the development of drug-resistant TB: acquired and transmitted drug resistance. The acquired drug resistance model suggests that resistance is developed within patients with active TB through microevolution(27). Several studies suggest examples for microevolution(28),(99) especially those involving the hypermutable Bejing Mtb lineage(100). However, it is crucial to note that distinguishing between acquired and transmitted resistance is not straightforward based solely on allele variants found in the sputum. In the transmitted resistance model, a patient accumulates a pool of Mycobacteria with different genotypes during latent infection. This population mix is essentially clonal, as M. tuberculosis strains possess a highly conserved core genome(14), but with several genetic allele variants having limited representation. The transition of the disease to an active phase, along with subsequent chemotherapy, leads to adaptive selection from the pre-existing pool of variants. The concept that certain TB cases involve mixed infections has been substantiated in clinical cases using phage typing and whole-genome sequencing(101, 102). Transmissibility of resistant variants has been confirmed through strain specific PCR(103), and selective adaptation in a patient during chemotherapy has also been demonstrated (17). Furthermore, it has been shown that clonal complexity is reduced by culturing, leading to the underrecognition of polyclonal infections in culture-based diagnosis(104). The WHO estimates that a quarter of the world’s population is latently infected by M. tuberculosis, accumulating different TB strains throughout their lives(105). Consequently, patients may harbor high heterogeneity, facilitating the spread and fixation of a genetic variant that holds some advantage in specific environmental conditions.

We acknowledge the limitations of using M. smegmatis as a model for the intracellular pathogen M. tuberculosis, which is associated with complex pathology. Nevertheless, given the conserved molecular mechanisms of genome maintenance in Mycobacteria, we can conclude that the mycobacterial genome is not prone to microevolution upon prolonged exposure to the antibiotics employed in our study and the clinics.

5. Conclusion

We showed that prolonged exposure to clinically used antibiotics did not lead to de novo adaptive mutagenesis in M. smegmatis under laboratory conditions. The activation of DNA repair pathways preserves genomic integrity, while non-genetic factors convey quick adaptation for stress conditions.

Supplementary Data are available at NAR online.

6. Data availability

The data underlying this article are available in the Figshare repository (using a private link shared with the Editor for now); DOI: 10.6084/m9.figshare.25028186 (reserved DOI until acceptation).

Acknowledgements

We would like to express our gratitude for the assistance provided by Ádám Póti in the analysis of the WGS data.

8. Funding

This work was supported by the National Office for Research and Technology, Hungary [grant number OTKA-K115993 and OTKA-K138318 to J.T., OTKA-PD128254 to R.H.]. Funding for open access charge: National Office for Research and Technology, Hungary