Introduction
Mycobacterium tuberculosis is the causative agent of tuberculosis (TB), which remains a major public health threat. However, the treatment of
M. tuberculosis infections is a significant challenge. The World Health Organization (WHO) has estimated a treatment success rate of 85% in patients with drug-susceptible TB (DS-TB) who receive a 6-month treatment regimen consisting of 4 first-line TB drugs (i.e., isoniazid, rifampicin, pyrazinamide, and ethambutol) [
1]. In general, the WHO does not recommend extending the duration of treatment beyond 6 months, as it does not significantly increase the efficacy of the drugs [
2]. However, the treatment duration may need to be extended depending on the patient’s age, comorbidities, and treatment progress [
3].
Furthermore, the treatment is complicated by heterogeneity within
M. tuberculosis populations [
4,
5]. Genetic diversity results in the formation of subpopulations within
M. tuberculosis strains, and these subpopulations within a host can result from infection by multiple strains, reinfection with a new strain, or the microevolution of an infected strain within the host [
6–
8]. The heterogeneous
M. tuberculosis subpopulations exhibit different phenotypes—for instance, the formation of persister cells, asymmetric cell growth and division, and variations at the transcriptional and posttranslational levels [
9–
11]. During pathogenesis, phenotypically diverse subpopulations exhibit varying levels of drug susceptibility and resistance, which is the primary reason for the prolonged treatment duration [
9].
Persisters are a subpopulation without genetic resistance that temporarily halts their growth and sustains it in the presence of antibiotics and environmental stressors [
10–
15]. The presence of
M. tuberculosis persister cells within the host hinders the effective eradication ability of therapeutic drugs and affects drug efficacy [
16–
18]. Several studies have been conducted to elucidate the mechanisms of persister formation. Transcriptome analysis of
M. tuberculosis persisters revealed downregulation of metabolic and biosynthetic pathways and upregulation of toxin-antitoxin (TA) modules [
10,
19]. In addition, it has been shown that differential expression patterns of genes involved in stringent responses and catalase-peroxidase (
katG) affect drug tolerance [
20,
21]. Using whole-genome sequencing and transcriptome analysis of the high persister mutants and clinical isolates of
M. tuberculosis, a recent study showed that genes involved in lipid biosynthesis, carbon metabolism, TA systems, and transcriptional regulators are associated with persister formation [
10]. Specifically, genetic variation of the genes involved in outer membrane lipid phthiocerol dimycocerosate biosynthesis increases drug tolerance [
10,
22].
The persistence of the clinical isolates of
M. tuberculosis has been confirmed after isolation from sputum samples
in vitro. The isolation of
M. tuberculosis clinical isolates from sputum samples requires several stages of subculture. The culture process can influence the population structure, such as variant genotypes and clonal composition, resulting in the loss of genetic diversity and modification of clinical specimens [
23–
25]. Moreover, drug treatment influences the population structure of
M. tuberculosis [
26]. Therefore, it is important to select the time point for isolation to properly characterize infected
M. tuberculosis.
This study aimed to identify the molecular factors that contribute to prolonged DS-TB treatment using genomic and transcriptomic profiling. Clinical isolates of M. tuberculosis were collected from the first mycobacterial growth indicator tube (MGIT) culture medium inoculated directly into sputum samples to characterize the infected M. tuberculosis.
Materials and Methods
Patients, Samples, and Culture Conditions
Clinical isolates of
M. tuberculosis were obtained from patients’ sputum before the standard treatment regimen or 2 weeks after treatment.
M. tuberculosis clinical isolates were identified from sputum samples using both Löwenstein-Jensen medium and the MGIT culture system. Acid-fast bacilli (AFB) staining with Ziehl-Neelsen stain and MPT64 antigen detection using Bioline TB Ag MPT64 (Abbott) were used to confirm positive culture isolates. AFB-stained samples or MPT64 detection-negative samples were assessed by polymerase chain reaction (PCR) using a PowerChek
M. tuberculosis/nontuberculous mycobacteria (MTB/NTM) Real-time PCR Kit (Kogene Biotech). DS-TB and drug-resistant (DR)-TB were differentiated using GeneXpert MTB/RIF (Cepheid), GenoType MTBDR
plus (Hain Lifescience GmbH), and traditional minimum inhibitory concentration measurements. In addition, 13 and 9 clinical isolates were subjected to genome and transcriptome analyses, respectively. The patient details are presented in
Table 1. To maintain the population structure of the clinical isolates of
M. tuberculosis, preserved sputum samples were directly inoculated into the MGIT after identification. The MGIT culture medium was collected by centrifugation at 4,000 rpm for 3 minutes. The harvested cultures were used to isolate genomic DNA and total RNA.
Genome Sequencing
Genomic DNA was isolated using a Quick-DNA Fungal/Bacterial Miniprep Kit (Zymo Research) as per the manufacturer’s instructions. The TruSeq Nano DNA Kit (Illumina Inc.) was used to generate paired-end (2×150-bp) libraries. Whole-genome sequencing of the libraries was performed using an Illumina NovaSeq 6000 platform. Reads generated through sequencing were qualitatively validated and trimmed using Prinseq v0.20.4 [
27]. Sequencing reads were aligned to the
M. tuberculosis H37Rv reference genome (NC_000962.3) using CLC Genomics workbench 22.0 (Qiagen). This whole-genome sequencing project has been deposited in the National Center for Biotechnology Information (NCBI) Short Read Archive (SRA) and is accessible through the Bioproject accession number PRJNA887480.
RNA Sequencing
Total RNA was isolated using the RNeasy Plus Mini Kit (Qiagen) as per the manufacturer’s instructions. RNA sequencing libraries were created using the TruSeq rapid SBS kit v4 (Illumina), and sequencing of the libraries was conducted using the HiSeq 2500 (Illumina) sequencing protocol. Quality validation and trimming of the reads were performed using FastQC v0.11.7 and Trimmomatic v0.38 [
28], respectively. Paired-end reads (101 bp) were then mapped to the reference genome sequence of
M. tuberculosis H37Rv (NC_000962.3) using CLC Genomics workbench 22.0 (Qiagen). The RNA sequencing data have been deposited in the NCBI SRA and are accessible through the Bioproject accession number PRJNA1010683.
Ethics Approval
This study was approved by the Institutional Review Board of the Korea Disease Control and Prevention Agency (approval no: 2019-04-05-2C-A), and patients were recruited from a tuberculosis cohort from 2016 to 2018.
Results
Characteristics of DS-TB Patients
A total of 600 participants were screened from December 7, 2016, to October 8, 2019. Of these patients, 334 were excluded from the analysis because of NTM infection, TB without pulmonary involvement, death, transfer out, or withdrawal of consent (
Figure 1). Among the remaining 266 patients, 132 were diagnosed with DR-TB and were excluded from the analysis. Of the 134 patients with DS-TB, 48 (35.8%) completed the standard 6-month treatment and 86 (64.2%) received >6 months of treatment. This indicates that the number of patients who received >6 months of treatment was almost twice the number of patients who completed the standard 6-month treatment.
Patient characteristics are described in
Table 2. The mean age of patients treated with the 6-month regimen was 50.9±19.9 years, and there were 26 females (54.2%) and 22 males (45.8%). In patients who were treated for >6 months, the mean patient age was 64.0±19.5 years, with 37 females (43.0%) and 49 males (57.0%). In addition, 2 patients (4.2%) in the 6-month treatment group and 16 patients (18.6%) in the group treated for >6 months had a history of TB treatment. There were 24 patients (50.0%) with comorbidities in the 6-month treatment group and 55 patients (64.0%) in the group treated for >6 months. In the 6-month treatment group, 17 patients (35.4%) had pulmonary cavity positivity, 35 patients (72.9%) had chest X-ray positivity, and 2 patients (4.2%) had a positive AFB smear at 2 months. In the group treated for >6 months, 37 patients (43.0%) were positive for the pulmonary cavity, 72 patients (83.7%) were chest X-ray positive, and 6 patients (7.0%) had a positive AFB smear at 2 months.
Whole-Genome Sequencing Analysis of M. tuberculosis Isolated from Sputum Samples
We attempted to isolate
M. tuberculosis from patient sputum samples collected 2 weeks after treatment for which bacterial culture was possible.
M. tuberculosis could be isolated from the samples of 16 of the 134 DS-TB patients (data not shown). Analysis was conducted using samples from the 16 selected patients. Thirteen
M. tuberculosis isolates from 10 patients were subjected to whole-genome sequencing (
Table 1).
Table 3 shows the total number of non-synonymous single-nucleotide polymorphisms (SNPs) between the clinical isolates and the
M. tuberculosis H37Rv reference genome. The mean of the total number of non-synonymous SNPs was 676.5±264.9 and 875.6±28.3 in the 6-month and >6-month treatment groups, respectively (
Figure 2). Interestingly, 2 isolates in the 6-month treatment group and all isolates in the group treated >6 months were identified as belonging to the East-Asian lineage. These isolates had larger numbers of non-synonymous SNPs than the Euro-American lineage isolates (CB050 and IP004). The large standard deviation of the 6-month treatment group implied that clinical isolates identified as East-Asian lineage are more genetically diverse than the Euro-American lineage isolates in our data.
Non-synonymous SNPs were analyzed to identify candidate genes that influenced treatment duration. We found 1 non-synonymous SNP in
pe_pgrs9 (
rv0746) and
ppe34 (
rv1917c) in the group treated for >6 months (
Figure 3). The average frequencies of non-synonymous SNPs in
pe_pgrs9 and
ppe34 were 44.1%±4.4% and 40.0%±6.8%, respectively. This suggested that approximately 40% of the subpopulation had single non-synonymous SNP in
pe_pgrs9 and
ppe34, but only in the group treated for >6 months.
Transcriptome Analysis of M. tuberculosis Isolated from Sputum Samples
Nine clinical isolates from 8 patients were used for transcriptome analysis (
Table 1). The transcriptomes of the group treated for >6 months and the group treated for 6 months were compared using RNA deep sequencing. In the group treated for >6 months, 165 genes were differentially expressed (
p<0.05; fold change ≥2) relative to the group treated for 6 months, of which 26 genes were induced and 139 were repressed (
Figure 4A;
Table S1).
To identify the function of differentially expressed genes (DEGs), functional categories of the DEGs were annotated by Mycobrowser (
http://mycobrowser.epfl.ch/) and the database for annotation, visualization and integrated discovery (DAVID) (
https://david.ncifcrf.gov/). Genes involved in virulence, detoxification, adaptation, cell wall formation, and cell processes were induced (
Figure 4B). Gene ontology (GO) analysis revealed that ribonuclease activity and magnesium ion binding were enriched in the induced genes (
Figure 4C). Virulence, detoxification, and adaptation categories from Mycobrowser analysis, and the ribonuclease activity category from DAVID analysis included virulence-associated protein (Vap) family TA genes, such as
vapBC26,
vapB38, and
vapBC44. Repressed DEGs were associated with cell wall formation and cell processes, intermediary metabolism and respiration, and information pathways (
Figure 4B). In these categories, reductions in ribosomal protein genes and peptidoglycan biosynthesis involving
murE-X genes were observed. According to GO analysis, the reduced genes were enriched in the structural constituent of the ribosome, the acyl carrier protein phosphopantetheine attachment site binding involved in the fatty acid biosynthetic process, acyl binding, hydrolase activity, and hydrolyzing N-glycosyl compound categories (
Figure 4D). Based on the functional categorization results, the expression of TA module genes, rRNA and ribosomal protein genes, and proline-glutamate (PE)-proline-proline-glutamate (PPE) family genes are shown as a heatmap (
Figure 5). Among the TA module genes induced in the group treated for >6 months, Vap family genes were the most common, and various TA family genes were included in the reduced genes. The expression levels of rRNA genes (
rrs, rrl, and
rrf), 8 ribosomal protein genes (
rpmC, rpsQ, rpsN1, rplR, rpmD, rpmA, rplU, and
rpmJ), and 12 PE-PPE family protein genes (
pe_pgrs4, ppe29, ppe30, pe_pgrs32, ppe32, ppe38, ppe39, pe_pgrs42, pe29, lipY, ppe53, and
pe_pgrs 61) were reduced in the group treated for >6 months. This suggests that the TA module genes of several Vap families were upregulated in the group treated for >6 months, and peptidoglycan biosynthesis-involved rRNA and ribosomal protein genes were downregulated.
Discussion
The genetic diversity of
M. tuberculosis influences its physiological outcomes. For example, TB is caused by
M. tuberculosis, which includes 9 lineage classifications based on large sequence polymorphisms and SNPs. Genomic differences between lineages affect pathogenesis, host immune responses, and transmissibility [
29–
32]. In addition, the presence of
M. tuberculosis subpopulations within the host affects cell growth, metabolic status, responses to environmental conditions, and drug use [
9]. A subpopulation of drug-tolerant
M. tuberculosis called persisters exhibits variable drug susceptibility. Persisters can avoid drug pressure through growth arrest and regrowth after the removal of bactericidal drugs [
11,
33]. These properties lead to long-term antibiotic treatment to eradicate persister cells [
11]. Therefore, understanding the genetic diversity of
M. tuberculosis is important for effective treatment.
The involvement of TA modules and stringent response regulator A (RelA)-mediated stringent responses in persister formation has been well characterized in other bacteria [
34–
36]. Stringent response is a bacterial response to nutrient starvation and other stress conditions and is mediated by the accumulation of ppGpp, which is synthesized by RelA. Accumulated ppGpp stimulates protease Lon-dependent degradation of type Ⅱ TA module antitoxins, resulting in endonuclease toxin-mediated translational inhibition and persister cell formation [
35,
37]. Under these stress conditions, expression of toxin genes is upregulated, and activated toxins inhibit important cellular processes, leading to growth arrest [
38]. Several mechanisms have been proposed for the formation of mycobacterial persister cells. Like other bacteria, persister formation in
M. tuberculosis is associated with TA modules and RelA [
19,
20,
36,
39,
40]. In
M. tuberculosis, the expression of
higAB,
relFG,
vapBC3, and
vapBC31 TA module genes are upregulated in drug-tolerant subpopulations [
19]. In addition, cholesterol-induced activation of the RNase toxin VapC12 inhibits translation, resulting in growth regulation and increased persister populations [
41]. Three
relE toxin homologs are involved in the formation of drug-tolerant persisters
in vitro [
20]. In
M. tuberculosis, RelA is required for growth arrest and drug tolerance during nutrient deprivation [
40]. Furthermore, RelA downregulates ribosomal protein genes and mediates the long-term persistence of
M. tuberculosis in vivo [
36]. Recent studies have shown that the expression of
vapC30 and
mazF toxins and
relA is induced in single-cell dormant persisters [
39]. Our RNA sequencing results showed the upregulation of several Vap family TA module genes, especially
vapBC17,
vapBC26, and
vapBC44 in the group treated for >6 months (
Figures 4 and
5). Moreover, VapC26 has ribonuclease activity targeting 23S rRNA, and overexpression of the
vapC26 gene strongly inhibits cell growth in
M. tuberculosis [
42,
43]. The downregulated genes in our study included rRNA genes in the group treated for >6 months (
Figure 5). The expression of genes encoding ribosomal proteins is downregulated in persisters [
19,
36]. Considering these data, it can be assumed that
M. tuberculosis isolated from the group treated for >6 months showed a higher persistence rate in the subpopulation than in the group treated for 6 months.
Genomic profiling of the 2 groups revealed a single non-synonymous SNP in
pe_pgrs9 and
ppe34 only in the group treated for >6 months (
Figure 3). These 2 genes that encode PE-PPE family proteins account for approximately 10% of the
M. tuberculosis genome [
44]. The PE-PPE family is characterized by highly conserved N-terminal domains with PE and PPE motifs and is classified into PPE, PE, and PE Polymorphic GC-rich subfamilies (PGRS). Several studies have suggested that the PE-PPE family proteins function in mycobacterial pathogenesis, particularly in the modulation of host immune responses. For instance, the PE-PPE family of proteins binds to toll-like receptors (TLRs), which induces the maturation of dendritic cells (DCs) and macrophage activation, and modulate type 1 T helper (Th1) and type 2 T helper (Th2) responses [
45–
47]. Some PE-PPE family proteins induce immune evasion of intracellular
M. tuberculosis by eliciting the Th2 response and inhibiting autophagic killing mediated by interferon-γ secreted from Th1 cells [
48,
49].
Proteomic analysis was performed to identify important proteins during
M. tuberculosis infection in guinea pig lungs with PE_PGRS9 and PPE34; however, the functions of these 2 proteins have not been well studied [
50]. Although the function of PE_PGRS9 remains inconclusive, it has been demonstrated that PPE34 is involved in immune evasion. Surface-exposed PPE34 induces DC maturation in a TLR2-dependent manner, and PPE34-matured DCs secrete interleukin (IL)-10, IL-4, and IL-5 to skew the Th1/Th2 balance toward Th2-preferred cellular immunity [
45]. These PPE34-driven Th2 immune responses are induced by the Src homology 3-interacting domain of PPE34 (amino acids 348–824) [
45]. Our genome sequencing data showed that approximately 40% of the subpopulation in the group treated for >6 months had a substitution of isoleucine for methionine at position 449, located in the SH3-interacting domain of PPE34. This suggested that there is an association between a single non-synonymous SNP in the
ppe34 gene and the duration of treatment; however, a detailed study is required. In addition, RNA sequencing data revealed a significantly decreased expression of the
ppe38 gene in the group treated for >6 months (
Figure 5). In
M. tuberculosis, it has been demonstrated that PPE38 mediates the secretion of the PE-PGRS proteins and deletion of the
ppe38 gene increases virulence [
51]. In this respect, the decrease in
ppe38 gene expression in the group treated for >6 months appears to be responsible for the increased virulence of infected
M. tuberculosis, thus affecting treatment duration.
Our data showed that non-synonymous SNP in ppe34 and decreased ppe38 expression were associated with long-term treatment. Transcriptome analysis revealed upregulation of the ribonuclease toxin vapC26 gene and downregulation of rRNA and ribosomal protein-encoding genes in the prolonged treatment group. These characteristics of M. tuberculosis from the group treated for 6 months may affect treatment outcomes by modulating the host immune response and increasing the virulence and frequency of M. tuberculosis persister cells. However, our study demonstrated the molecular characteristics of drug-susceptible M. tuberculosis clinical isolates with extended treatment duration using a small number of samples. Further research is needed to determine whether the genetic variations and expression changes we identified influence the extended treatment.
Article information
Ethics Approval
We obtained approval from the Institutional Review Board of the Korea Disease Control and Prevention Agency (approval no: 2019-04-05-2C-A) and isolated M. tuberculosis from tuberculosis patients from 2016 to 2018. All clinical M. tuberculosis isolates were collected from the sputum of tuberculosis patients at Chungnam National University Hospital, Chungbuk National University Hospital, and Inje University Ilsan Paik Hospital in Korea.
Conflicts of Interest
The authors have no conflicts of interest to declare.
Funding
This research was supported by the National Institute of Health research project (project No. 2018-NI-003-01 and No. 2022-NI-012-01).
Availability of Data
The datasets generated during the current study has been deposited in NCBI and accessible through the Bioproject accession number PRJNA887480 and PRJNA1010683.
Authors’ Contributions
Conceptualization: EMK, JM, SL, SK; Data curation: EMK; Formal analysis: EMK; Funding acquisition: JAJ, SL, SK; Investigation: HK; Methodology: HK, EMK; Project administration: JAJ, SL, SK; Resources: JM, SL, SK; Software: EMK; Supervision: JM, SL, SK; Validation: EMK; Visualization: EMK; Writing–original draft: EMK; Writing–review & editing: authors. All authors read and approved the final manuscript.
