Tuberculosis (TB) is a global health threat caused by the highly successful human pathogen Mycobacterium tuberculosis (M. tuberculosis). According to a report by the World Health Organization (WHO), an estimated 10.6 million new TB cases occurred worldwide in 2022, resulting in over 1.3 million deaths (World Health Organization, 2023). Despite the substantial global burden of TB, our knowledge regarding the factors influencing its transmission remains limited. Therefore, it is imperative to delve deeper into the mechanisms underlying the spread of M. tuberculosis.

The toxin-antitoxin (TA) system plays a critical biological role in M. tuberculosis. Composed of toxins and antitoxins, this system forms a small genetic unit that is widely present in prokaryotes (Schuster and Bertram, 2013; Dai et al., 2022). TA systems have been shown to assist cells in stress adaptation, antibiotic resistance, biofilm formation, persisters, and disease development. Toxins are typically translated into proteins, while antitoxins can be either proteins or RNA (Ogura and Hiraga, 1983; Aizenman et al., 1996; Magnuson, 2007; Fineran et al., 2009; Wang and Wood, 2011; Lobato-Márquez et al., 2016). Based on the nature of antitoxins and the mechanisms which inhibit toxin activity, TA modules can be classified into six distinct types (Page and Peti, 2016). Among these types, type II TA systems are well-characterized, where antitoxins directly interact with toxins to neutralize their effects. Bioinformatics and phylogenetic analyses have revealed the presence of numerous TA systems encoded in the M. tuberculosis genome. The retention of these TA systems in members of the M. tuberculosis complex suggests their crucial role in regulating metabolic pathways essential for bacterial pathogenesis. Type II TA systems predominate in M. tuberculosis. The abundance of TA loci in the M. tuberculosis genome raises important questions about their functional diversity (Ramage et al., 2009; Tandon et al., 2019). Previous studies have extensively investigated the various functions of TA systems in M. tuberculosis and their potential impact on pathogenic mechanisms (Schippers et al., 2005; Guo et al., 2016). These systems are believed to play a key role in M. tuberculosis ‘s response to stressors such as nutrient starvation and antibiotic treatment, promoting its survival and drug resistance (Kim et al., 2018). Additionally, TA systems are associated with the formation of persistent cells, a subpopulation exhibiting drug tolerance that plays a crucial role in establishing chronic infections in M. tuberculosis (Merfa et al., 2016). While the importance of toxin-antitoxin systems in M. tuberculosis has been acknowledged, our understanding of their specific mechanisms and functions within this bacterium remains limited. Therefore, comprehensive research is required to explore the roles of TA systems and gain deeper insights into the complex biology of M. tuberculosis.

Driven by the need to better understand the mechanisms underlying M. tuberculosis transmission, we conducted an extensive study investigating the impact of mutations in TA system genes on its spread. Our research aims to elucidate how genetic variations within this system can influence M. tuberculosis strain transmission dynamics. Utilizing whole-genome sequencing (WGS), we analyzed the genetic variations present in M. tuberculosis isolates at a high-resolution level. This enabled us to identify specific mutations within the TA system genes that may be associated with M. tuberculosis transmission. Advanced statistical and bioinformatics techniques, including random forest, gradient boosting decision tree, and generalized linear mixed models, were employed for comprehensive analyses to identify key genetic variants linked to transmission dynamics. We acknowledge challenges posed by confounding factors and population dynamics in our analysis. Future research should incorporate social networks and regional interactions for a more comprehensive understanding. Limitations of our study include a focus on gene analysis, potentially overlooking other important genetic influences such as drug resistance mutations or virulence determinants. Therefore, more comprehensive studies are needed to address these limitations adequately. Our study has yielded significant results, identifying multiple single nucleotide polymorphisms (SNPs) within the toxin-antitoxin system genes that positively correlate with clustering, suggesting their potential role in M. tuberculosis transmission. Furthermore, some of these SNPs were found to be associated with transmission clades across different geographical regions, indicating their potential global impact on the spread of M. tuberculosis. These findings provide valuable insights into the transmission dynamics of this pathogen and contribute to a more thorough understanding of M. tuberculosis transmission.

Consistent with prior research findings, our study further emphasizes the diverse functionality of TA systems in M. tuberculosis. These redundant TA systems serve as a backup mechanism enabling cellular adaptation and survival under adverse conditions (Min et al., 2012). They play a critical role in M. tuberculosis’s stress response, including nutrient deprivation, by regulating essential cellular processes like DNA replication, protein translation, and cell division. Moreover, TA systems contribute to the formation of drug resistance and persistence in M. tuberculosis. However, it is important to acknowledge that certain studies have reported conflicting results regarding the specific contributions of TA systems to persistence formation and stress conditions (Yu et al., 2020; Sharma et al., 2021). These discrepancies may arise from variations in experimental setups or genetic differences among M. tuberculosis strains used in different investigations. Therefore, additional research is needed to precisely determine the roles of TA systems in persistence formation, stress responses, and their impact on M. tuberculosis pathogenesis. In our study, we focused on examining the relationship between gene mutations in toxin-antitoxin systems and the transmission dynamics of M. tuberculosis. The M. tuberculosis genome contains numerous toxin-antitoxin systems, including VapBC, MazEF, ParDE, and RelBE (Ramage et al., 2009; Tandon et al., 2019). To gain deeper insights into the significance of these toxin-antitoxin systems in M. tuberculosis transmission, we analyzed the prevalence and genetic variation of specific toxin-antitoxin system genes across various clusters and evolutionary branches. Our analysis detected multiple mutations in these genes, suggesting they could be involved in M. tuberculosis transmission.

In our study, we have found a strong association between SNPs in the VapB antitoxin-related genes and the transmission of M. tuberculosis. Specifically, we identified several significant SNPs that were linked to transmission, including vapB1 G34A, vapB24 A76C, vapB31 T112C, vapB14 A29G, and vapB15 (T6C, G237A). We observed that vapB24 A76C and vapB23 T2C were particularly associated with transmission, especially in lineage 2. Additionally, vapB2 T171C, vapB11 C54T, vapB14 A29G, vapB15 G237A, vapB17 G213C, and vapB20 A54C were significantly related to transmission, especially in lineage 4. Furthermore, we found that vapB43 G28C was associated with transmission in lineage 4, while vapB6 A29C, vapB31 T112C, and vapB37 G205C were correlated with cross-country and cross-regional transmission. We also found that vapB1 G34A, vapB11 C54T, vapB14 A29G, vapB15 T6C, and vapB46 G70A were related to clade size. The VapBC system is crucial for regulating the behavior and adaptation of M. tuberculosis under diverse environmental stresses. It comprises stable VapC toxins and labile VapB antitoxins, whose interplay is essential for bacterial growth, survival, and response to stress conditions (Robson et al., 2009; Winther and Gerdes, 2011). During periods of stress, antitoxin molecules are degraded, leading to the release of toxins, such as VapC, through their RNase activity (Min et al., 2012). Consequently, these toxins inhibit or slow down cellular metabolism, providing a survival advantage to the bacterium during adverse conditions. The delicate balance between VapB antitoxins and VapC toxins is crucial for maintaining bacterial homeostasis and ensuring appropriate responses to external stimuli (McKenzie et al., 2012). Overall, our study provides compelling evidence for the significant association between SNPs in VapB antitoxin-related genes and M. tuberculosis transmission. These findings shed light on the intricate role of the VapBC toxin-antitoxin system in regulating bacterial behavior and underscore the importance of genetic variations within this system in shaping transmission dynamics.

Our study has revealed the association between SNPs in other TA system genes and the transmission of M. tuberculosis. Specifically, we focused on the MazEF family, which consists of nine TA systems encoded in an operon (Ahmed et al., 2022). We found a close connection between the mazE6 G156A and mazF2 C85T gene polymorphisms and the transmission clusters, particularly within lineage 2. These variants exhibited significant correlations with the formation and expansion of transmission clusters. However, mazE6 G156A is a synonymous mutation (Arg52Arg), which does not directly alter the protein’s function but may still affect the protein through other mechanisms. For example, in certain situations, synonymous mutations can lead to changes in transcription regulatory elements, thereby influencing gene expression levels. However, further research is needed to confirm these effects. Similarly, the mazF2 C85T variant may alter the stability of the MazF and modulate the delicate balance between toxin and antitoxin interactions (Leplae et al., 2011). Furthermore, our study identified a strong correlation between the mazE2 G104A, mazE5 T344C, and mazF8 A97G gene polymorphisms and the transmission clusters, especially within lineage 4 and clade size. While it’s plausible that these genetic variations influence the MazEF system activity, stability, and domain structure, our ability to fully elucidate these mechanisms is currently limited. Therefore, it’s crucial to interpret these functional implications cautiously and consider other potential contributory factors to M. tuberculosis transmission. Furthermore, no SNPs in the MazEF system were found to be associated with cross-country and cross-regional transmission of M. tuberculosis in our study. Future investigations should aim to provide a more comprehensive understanding of these effects, confirm these hypotheses, and uncover the precise impact of these mutations on the dynamics of M. tuberculosis transmission.

The ParDE toxin-antitoxin system in M. tuberculosis plays a crucial role in bacterial transmission dynamics. Our research has identified specific genetic variations in the parE and parD genes, such as parE1 C88T, parE2 C48G, parE1 C103T, parD2 A196G, parE1 C25G, and parD1 C134T, that are closely linked to transmission clusters, particularly within lineage 4 and lineage 2. These genetic variants impact cross-country and cross-regional transmissions, highlighting the significance of the ParDE system in the spread of M. tuberculosis. Variations in the parD gene, including those involving Rv2142A (parD2) and Rv1960c (parD1), can modify the activity and regulatory mechanisms of the ParD antitoxin. Similarly, variations in the parE gene, particularly those affecting Rv1959c (parE1), influence the function and stability of the ParE toxin, thus impacting its interaction with the ParE antitoxin (Xu et al., 2018). Understanding these genetic interactions is crucial for deciphering M. tuberculosis transmission dynamics and developing targeted interventions to effectively combat tuberculosis. Additionally, our research has identified a unique SNP, T226A, in the relB gene that is associated with transmission clusters in M. tuberculosis. This genetic variation further adds to the complexity of bacterial transmission dynamics, highlighting the intricate interplay between genetic factors and the spread of M. tuberculosis.

In terms of drug development and therapeutic interventions, our research findings could potentially have significant implications. The diverse functions of TA systems suggest potential targets for novel therapeutic strategies in M. tuberculosis. Understanding the relationship between genetic variations and functional consequences within these TA systems might help us discover new methods to disrupt or modulate their activity, thereby affecting the survival and transmission dynamics of the bacterium. Firstly, interventions targeting specific SNPs in TA systems such as VapBC, MazEF, ParDE, and RelBE could possibly directly alter the stability and activity of their toxins or antitoxins, thus impacting the growth, survival, and adaptability of M. tuberculosis (Robson et al., 2009; Leplae et al., 2011; Winther and Gerdes, 2011; McKenzie et al., 2012; Ahmed et al., 2022). The SNPs we discovered, including vapB24 A76C and vapB23 T2C, have the potential to serve as genetic markers for targeted drug design, allowing for more personalized treatment approaches. Additionally, mutations like parE1 C88T, parE2 C48G, and parE1 C103T show associations with cross-national and cross-regional transmissions of M. tuberculosis, which could aid in the development of more effective treatment plans to reduce global transmission. However, it’s important to note that while these genetic insights hold potential, they still require experimental validation to confirm their clinical significance and functional implications. Each mutation may lead to different functional impacts, and there might be other complexities involved, such as drug tolerance or adaptability of the bacterium under different environmental conditions. Therefore, further research is needed to delve deeper into the functional impacts of these genetic variations and precisely determine their roles in new drug development and treatment strategies. It is crucial to validate these findings through rigorous experimental studies and clinical trials before implementing them in clinical practice. Future research should aim to elucidate the specific mechanisms underlying these genetic variations and their contributions to drug response and transmission dynamics. By gaining a better understanding of the functional implications, we can more accurately tailor treatment strategies and contribute to the development of more targeted and effective interventions.

Our findings emphasize that both synonymous and non-synonymous mutations can influence the transmission of M. tuberculosis, suggesting that synonymous mutations in TA system genes are not universally neutral, in line with prior research by Shen et al. (2022). We believe that synonymous mutations may affect mRNA stability, splicing, or secondary structure formation. Changes in these regulatory elements can influence gene expression patterns and protein folding, thereby impacting bacterial adaptability and transmission capacity. Additionally, synonymous mutations may be part of a compensatory mechanism. While synonymous mutations themselves may not directly provide selective advantages, they may be associated with compensatory changes in other regions of the genome. These compensatory mutations could restore proper interactions between proteins, maintain enzyme activity, or optimize cellular functions affected by primary mutations, ultimately enhancing transmission capacity. Although the specific mechanisms and advantages of synonymous mutations in tuberculosis transmission are not yet fully understood, we cannot overlook their potential significance. Future research should consider the functional consequences of synonymous mutations and explore their interactions with other genetic factors, including non-synonymous mutations, drug resistance mutations, or virulence determinants. In our study, we combined local and global datasets to increase sample size for robust analysis of M. tuberculosis genetic variations. This approach helped identify shared and distinct variants across regions, enhancing our understanding of global pathogen diversity. Despite potential limitations such as variability from different protocols and sequencing technologies, stringent quality control measures, including SNP filtering within repetitive regions, were applied to minimize biases. Our novel findings contribute valuable insights into global M. tuberculosis genetic characteristics, advancing knowledge on tuberculosis pathogenesis and evolution. In future research, separate and comparative analyses of local and global data can be considered to highlight region-specific variations.

In our study, we investigated the impact of mutations in TA system genes on tuberculosis transmission. However, it is crucial to acknowledge that these correlations alone do not establish a causal relationship and should be interpreted with caution. Our modeling approach has limitations, notably in addressing potential confounding factors, such as population mobility, social networks, and inter-regional interactions. These elements may influence M. tuberculosis transmission but were not fully integrated into our models. We recognize that our primary focus on mutations within TA system genes may have led us to overlook other significant genetic influences, including SNPs related to drug resistance mutations or virulence determinants. While our findings contribute to the growing body of knowledge regarding the impact of toxin-antitoxin system gene mutations on tuberculosis transmission, further research is necessary to explore these intersections and understand their functional significance in detail. Limitations also arise from the sheer number of genes and computational resources required, which restricted our ability to analyze SNPs beyond the scope of our current investigation. Moreover, we lack a clear understanding of the cross-interactions and mutual regulation among the TA systems of M. tuberculosis, adding another layer of complexity to our study. Additionally, uncertainties inherent in the phylogenetic inference method used, such as homoplasy or recombination events, can present challenges when accurately determining evolutionary relationships. Therefore, future research should consider alternative methods to validate these findings and develop a more nuanced understanding of tuberculosis transmission. Further experimental validation is necessary to confirm the specific impact of TA system gene mutations. Future investigations should focus on refining our models to account for potential biases or shortcomings, and expanding research scope to explore the functional significance of these mutations and their direct influence on tuberculosis transmission.

We also discuss the limitations of using H37Rv as a single reference genome for analyzing M. tuberculosis WGS data, particularly regarding virulence and transmission. Recent studies suggest that relying solely on H37Rv may not fully capture the virulence characteristics of M. tuberculosis. H37Rv, commonly used as a reference genome in molecular epidemiology and drug resistance studies, does not represent the genetic diversity and variations present across all M. tuberculosis strains. Polymorphic loci involving genes associated with pathogenicity and host immune response, such as phospholipase C, membrane lipoproteins, adenylate cyclase gene family members, and PE/PPE gene family members, show significant differences between H37Rv and clinical isolates. Several gene families, including PE/PPE, exhibit higher substitution frequencies compared to the entire genome. Widespread genetic variability is observed at these polymorphic loci among M. tuberculosis clinical isolates (Fleischmann et al., 2002; O’Toole and Gautam, 2017). Phylogenetic and epidemiological analyses reveal independent occurrences of these polymorphisms, suggesting selective pressures driving these changes. Future research should incorporate genome sequences of additional reference strains, especially those directly obtained from clinical isolates, to comprehensively understand factors related to M. tuberculosis virulence and enable further investigations. For drug resistance inference, our analysis primarily utilized the TBProfiler platform. While incorporating additional tools/methods such as PhyResSE or bioinformatic SNP analysis could enhance robustness, resource constraints limited their implementation in this study. Thus, our results should be interpreted within the context of utilizing TBProfiler alongside the WHO-endorsed catalog. Future studies with expanded resources could consider alternative tools/methods for validation and complementation.

The results of this study suggest that mutations in toxin-antitoxin genes may increase the risk of M. tuberculosis transmission, underscoring the significance of conducting further research to explore the impact of these mutations on M. tuberculosis control and transmission. These findings offer new insights into the development of drug treatment strategies against tuberculosis.

The whole genome sequences have been submitted to the NCBI under the accession number PRJNA1002108.

YH: Conceptualization, Formal analysis, Software, Writing – original draft, Writing – review & editing. YiL: Conceptualization, Formal analysis, Methodology, Writing – review & editing. NT: Conceptualization, Formal analysis, Investigation, Project administration, Software, Validation, Writing – review & editing. XK: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Software, Validation, Writing – original draft. YamL: Conceptualization, Data curation, Formal analysis, Investigation, Methodology, Project administration, Software, Supervision, Validation, Writing – review & editing. YaoL: Data curation, Formal analysis, Methodology, Project administration, Supervision, Validation, Writing – review & editing. HL: Conceptualization, Data curation, Funding acquisition, Investigation, Resources, Visualization, Writing – original draft, Writing – review & editing. ZW: Conceptualization, Data curation, Formal analysis, Funding acquisition, Investigation, Project administration, Resources, Software, Supervision, Visualization, Writing – original draft.