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¹. The resurgence of TB due to the SARS-CoV-2 pandemic¹ 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² 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⁷. 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¹⁰. 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¹³. 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²¹. 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 exhibited a higher mutation rate³⁰. On the other hand, others found stable M. tuberculosis genomes with no or only a few emerging genomic changes over prolonged periods of treatment²⁶. Genotyping has enabled researchers to describe cases of co-infection by ≥2 different strains (mixed infection) or the 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²⁷. 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²⁷. 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³⁹, M. leprae⁴⁰, M. bovis^41,42, M. avium, M. paratuberculosis, and M. smegmatis⁴³ 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⁴⁶. In our earlier paper, we also compared the enzymes of thymidylate biosynthesis in M. tuberculosis and M. smegmatis and found high conservation ⁴⁷. 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⁴⁸. 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 to stress conditions. Notably, even with prolonged antibiotic exposure exceeding 230 bacterial generations, we observed no significant increase in the mutation rate, suggesting the absence of de novo adaptive mutagenesis.

2. Materials and methods

2.1. Bacterial strains, media and growth conditions

M. smegmatis mc²155 ⁴⁹ strains were grown in Lemco broth (5 g/l Lab-Lemco, 5 g/l NaCl, 10 g/l Bacto peptone, 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 peptone, 18.75 g/l Bacto agar).

2.2 Optimization of stress treatment conditions in liquid cultures and agar plates

The applied concentrations of drugs were optimized using serial dilutions of the compounds. In the case of liquid cultures, we monitored growth on a logarithmic scale by measuring the number of colony-forming units (CFU) or the optical density (OD) at 600 nm (Figure S1). The PZA treatments were done in acidic broth (pH=5.5 set using HCl). For agar plates, we determined the CFU of untreated mid-exponential phase (OD = 0.4-0.5) liquid cultures on both non-selective and drug-containing agar plates (Figure S2). We also monitored cell morphology in response to drug treatment. For further experiments, sublethal concentrations of drugs were chosen to obtain an 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 and/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 in 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 Td is the generation time, t is the time interval between measurements, and Nt and N0 are the final and initial population sizes, 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/)⁵⁰. The cell length 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 MA 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⁵¹. Treatments were performed for 60 days. We calculated 6.3±0.35 hours of generation time on the plate in this experimental setup. Thus, each line passed through ∼230 cell divisions. Some 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. Assessment of drug tolerance following MA experiments

The development of tolerance to the applied treatment was assessed by measuring the minimal inhibitory concentration (MIC) of both the mock-treated and stressed MA strains. Three randomly chosen strains from both the mock-treated and stress-treated groups were resuscitated on plates containing the same stress conditions as those used in the MA experiment. Liquid cultures were inoculated and diluted to an OD(600) of 0.001 in sterile, round-bottom 96-well plates (Sarstedt). The wells contained the specific drug in serial dilution for both the stressed strains and control samples. Cells were grown at 37°C without agitation. Plates were scanned and analyzed, and MIC values were determined based on the last well in which cell growth was observed.

2.7. DNA extraction

A single colony was inoculated into 10 ml liquid culture from all lineages, was grown until OD600 = 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.8. 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.9. 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 to analyse 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⁵². We mapped our paired end reads to M. smegmatis mc² 155 reference genome (GenBank accession number: NC_008596.1) by Bowtie2⁵³. PCR duplicates were removed with the use of Samblaster⁵⁴. We converted SAM files to BAM files, and sorted them with SAMtools⁵⁵. 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⁵⁵. 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 emerges 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, curated under the umbrella project at the ENA along with the present dataset.

2.10. RNA isolation and cDNA synthesis

For RNA extraction, cell pellets were resuspended in 2 ml RNA protect bacteria reagent (Qiagen; 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 Qiagen 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 Qiagen 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.11. 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⁵⁹. Using GeNorm^60,61 analysis, SigA and Ffh proved to be stably expressed in our experimental system (Figure S3).

2.12. 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 S3). 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 S4). The applied primers and their measured efficiencies are compiled in Table S3. The qPCR data were analyzed using the Bio-Rad CFX Maestro software and numerically shown in Table S4. 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.13. dNTP extraction

dNTP extraction and measurement were performed according to Szabo et al.⁶². 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.14. 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 S5) 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 MgCl2. 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/)⁶². 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 S6.

2.15. Tolerance assay

We used a modified version of fluctuation assays⁶³ for the estimation of the rate of emergence of tolerant cells upon preincubation with a sublethal dose of CIP (0.3 μg/ml). An initial 100 ml culture was grown to OD = 0.4-0.5 (3 biological replicates) , 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. 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 Bacto agar plates 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

The tolerance assay was repeated on 15 parallel plates for each biological replicate to obtain enough cells for genomic DNA extraction for WGS. Plates containing 0.3 μg/ml CIP were incubated for 96 hours at 37°C. Colonies were washed off the plates with 6 ml of Lemco broth. Genomic DNA was isolated and sent to WGS. CFU was also determined on Bacto Agar plates containing 0.5 μg/ml CIP and non-selective Bacto Agar plates (Figure S6).

2.16. 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 if not stated otherwise. F- and t-statistics were counted for the ΔCt values⁶⁴ of the qPCR results and for the concentrations normalized to the cell volume in case of the dNTP measurements. For the statistical analysis of the mutation rates, we used the t-test on the natural logarithm of the obtained mutation rate values.

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; ethambutol – EMB; rifampicin – RIF; pyrazinamide - PZA)³⁶. 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 measurable phenotypic effect while allowing to keep an adequate number of cells for the MA experiments on plate and for the downstream measurements in liquid culture (Figure S1-S2 and Table 1). In the first-line combination treatment, a 10-fold reduced concentration of each separately adjusted drug had to be applied in both liquid and agar media 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 in liquid cultures compared to the control (Table 1, Figure S1), consistent with their mechanism of action⁶⁷. To quantify the phenotypic effect of the applied drug treatments in liquid cultures, 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. The combination treatment did not affect the cell size (Figure 1 and Table 1).

Figure 1

We also assayed the clinically relevant drug pyrazinamide (PZA). However, M. smegmatis was reported to exhibit an intrinsic resistance to PZA⁶⁸. Indeed, PZA treatment alone, even at high concentrations in acidic conditions, did not affect cell viability in our experiments (Figure S1). Regardless of its inefficacy as a monotherapy, we included PZA in the combination treatment, as we could not rule out the possibility that PZA interacts with the other three drugs or that PZA elimination mechanisms are equally active in M. smegmatis under this regimen.

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 measured an average generation time of 6.3±0.35h on the plate within the timeframe of a single passage. Therefore, bacteria produced on average 230 generations during the 60-day treatment. 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 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.

A surprisingly few new mutations were detected after carefully cross-checking the sequencing data. 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). We calculated a 1 x 10^-10 mutation rate for our untreated M. smegmatis mc²155 strain. To our great surprise, the mutation rates of all treated lineages fell in one order of magnitude (4×10^-11 - 3×10^-10) except for the UV treatment used for positive control (Figure 2B).

Figure 2

We analyzed each mutation except for those obtained following the UV treatment and found no sign of adaptive changes (Table S2). The Excel file containing the positions of all obtained mutations, including those of the UV sample, is provided in the archive deposited for the article (DOI: 10.6084/m9.figshare.26585884).

We assessed the drug sensitivity of the MA strains by measuring the MIC of each drug on three randomly selected strains from both the mock-treated and stressed MA groups. Contrary to the mutation rate results obtained from genomic sequencing data, the MIC values for the MA strains were higher than those of the mock-treated strains (comparable data in line with⁶⁹, indicating phenotypic adaptation to the applied drugs (Fig 2C). However, for the EMB treatment, we observed no increase in MIC, despite repeating the experiment several times.

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

To investigate a possible reason for detecting so few newly generated mutations under antibiotic pressure, we studied whether the DNA repair pathways and other elements of the stress response potentially involved⁷⁰ were activated under drug pressure. The mycobacterial DNA repair system is highly redundant, many of its enzymes have overlapping functions^45,71,72. Although canonical mismatch repair proteins are thought to be missing, a recently described protein, NucS is encoded with similar function⁷³. 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 presented in Figure 3, grouped by functional relevance, with consistent heatmap coloring across all measurements. Figure S5 shows a clustered heatmap without prior functional grouping. Numerical data for expression level changes are provided in Table S4.

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 and S5). DNA polymerases DinB2 and DnaE2 involved in these DNA repair pathways are also strongly overexpressed (Figure 3 and S5). 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 and S5). 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 and S5). Interestingly, however, the DSB repair enzymes are exempt from this overall upregulation tendency (Figure 3 and S5). 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^74–77. Imbalanced or altered levels of dNTPs could cause an 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⁶². We used MMC treatment as a positive control as this is a generally used positive control for DNA damage^66,78. To calculate cellular concentrations, we used the cellular volumes determined from measured cell dimensions (Table S1). Interestingly, we found altered dNTP pools upon most treatments (Figure 4 and Table S6). 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 4F and H). RIF and MMC treatments promoted an increase in the dGTP and dCTP pools (Figure 4D-E). The INH treatment coincided with a decreased concentration of purine nucleotides (Figure 4B), while in EMB-treated cells we could measure very low levels of all dNTPs (Figure 4C and H). In the combination treatment, we could not measure significant differences (Figure 4A). The dGTP pool decreased in both absolute and relative terms across all treatments where dNTP pool changes were observed (Figure 4B-F and G, respectively). A smaller cell size coincides with a lower cellular dNTP concentration, while no clear correlation is observed between drug-induced cell length increase and dNTP pool expansion (Figure 4I).

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⁶³. 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 a 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. Figure 2B). Furthermore, following a 24-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

To confirm that the rapid increase in drug tolerance following short-term exposure to CIP is linked to non-genetic factors, we repeated the experiment using the 96-hour preincubation time for DNA isolation and WGS. After pretreatment, DNA was isolated from colonies on 5 parallel plates for each of the three biological replicates, followed by WGS (Figure 5A, Figure S6). In all measured samples, we detected a single mutation in a gene encoding an uncharacterized protein probably involved in lipid metabolism (MSMEG_6151; Table S2).

We also sequenced the genomes of colonies grown at the higher CIP concentration (0.5 µg/ml) and detected no mutations.

4. Discussion

One of the reasons 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 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^79,80,30,81. This low constitutive mutation rate by itself does not explain the biological diversity observed in clinical isolates²⁴. 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 possible molecular mechanism. We measured the appearance of drug-induced mutations in the genome in a mutation accumulation assay using WGS. 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 pretreatment, a tolerance rate on the order of 10^-7 /generation was observed for CIP, consistent with literature data from fluctuation assays^82,83). However, in the mutation accumulation assay, the number of mutations did not change significantly compared to the untreated control. The mutation rate increase was only significant in the case of the UV treatment serving as a positive control for the experiments (Figure 2B). 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,84. However, it should be noted that David’s study, which automatically classified bacteria growing in fluctuation assays as mutants without confirming genetic changes, also suggested that the term “acquired resistance” in tubercle bacilli has only practical meaning and lacks experimental foundation⁸². Our findings imply that the emergence of drug resistance in this study is solely attributed to phenotypic factors. Phenotypic changes upon antibiotic treatment have widely been investigated⁸⁵ including potential bistability⁸⁶ and/or the upregulation of efflux pumps^87,88. 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^85,89 and one study found a specific activation of the DNA repair system in response to CIP similar to ours⁶⁶. 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⁴⁵. 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⁹⁰. 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⁹⁰, 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⁹¹. High ROS levels are known to cause damage to nucleobases and the nucleotide pool is a major effector of oxidative stress-induced genotoxic damage⁹². 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 in 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⁹³ 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 4H-I) aligned well with the concurrent downregulated transcript levels (Figure 3). Resting states of bacteria have also been characterized by a decrease in cell size and dATP levels^95,96. These observations thus probably reflect the bacteriostatic effect of these drugs causing metabolic processes to enter a dormant state, accompanied by the downregulation of enzymes involved in dNTP synthesis. The combined treatment yielded the least significant alteration from the untreated control compared to all monotreatments (Figure 4). 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^97,98. 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⁹⁹ and mycobacterial genomes are characterized by a high GC content¹⁰⁰. 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 (but not significantly) 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. However, to confirm this hypothesis, these investigations should be conducted using genetically modified DNA repair mutant strains.

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¹. 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²⁷. Several studies suggest examples of microevolution^28,101 especially those involving the hypermutable Beijing Mtb lineage¹⁰². 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¹⁴, 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^103,104. The transmissibility of resistant variants has been confirmed through strain- specific PCR¹⁰⁵, and selective adaptation in a patient during chemotherapy has also been demonstrated ¹⁷. Furthermore, it has been shown that clonal complexity is reduced by culturing, leading to the underrecognition of polyclonal infections in culture-based diagnosis¹⁰⁶. The WHO estimates that a quarter of the world’s population is latently infected by M. tuberculosis, accumulating different TB strains throughout their lives¹⁰⁷. Consequently, patients may harbor high heterogeneity, facilitating the spread and fixation of a genetic variant with 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 may warrant the preservation of genomic integrity, while non-genetic factors convey quick adaptation to stress conditions.

6. Data availability

The data underlying this article are available in the Figshare repository (DOI: 10.6084/m9.figshare.26585884). Sequencing data are available at European Nucleotide Archive (ENA) with PRJEB71590 project number. Table S7 helps to identify the deposited files and provides information on the sequenced samples.

Acknowledgements

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

7. 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