Tuberculosis (TB) is one of the world’s most lethal illnesses. Every day, about 4000 people die as a result of TB. The enormous disruptions in health systems created by the COVID-19 pandemic are likely to escalate the dreadful toll by 1 million TB fatalities per year over the next 4 years . TB is a dangerous infectious illness that primarily affects the lungs and is caused by Mycobacterium TB (MTB). TB affects one-quarter of the world’s population, with 10.4 million new cases reported each year, despite the widespread use of a live attenuated vaccine and a variety of treatments. As an aerosol droplet, MTB penetrates the alveolar airways of exposed persons, where it is thought to establish its first contact with resident macrophages . The COVID pandemic, the emergence of multidrug-resistant (MDR) TB, and communities disproportionately restricted access to treatment all impact the worldwide burden of TB.
Both extensively drug-resistant TB (isoniazid, rifampin, and fluoroquinolone) and injectable drug (capreomycin or kanamycin) and complicated forms of MDR TB were shown to be ineffective in controlling threats of TB. Bedaquiline, a novel diarylquinoline, was added to the WHO-recommended all-oral regimen to replace injectable treatments for MDR-TB patients. The WHO consolidated guidelines on drug-resistant TB treatment . Bedaquiline has been reported to be a mycobacterial adenosine triphosphate (ATP) synthase inhibitor; it binds and disrupts the FO subunit interface, resulting in an inefficient proton cycle, which is lethal to Mycobacterium . When compared to comparable eukaryotic enzymes, the mycobacterial ATP synthase enzyme is highly selective (Selectivity Index >20, 000) . As a result, many people infected with TB remain misdiagnosed or are treated without test confirmation. Postmortem investigations demonstrated that TB is a common undetected cause of death, highlighting the need for better diagnostic tools. Microarray-based RNA expression analysis has evolved into an essential tool for studying disease biology. Many illnesses, such as cancer, infectious disorders, arthritis, HIV, and TB, are linked to distinct transcriptional patterns in blood or tissue .
It uses gene expression profiling to characterize complicated cellular responses, and new molecular pathways have been discovered under various circumstances. Because of the abundance of publicly available data, it is possible to pool gene expression datasets and boost sensitivity by increasing the number of data points. This method has already been utilized to discover gene signatures and pathways that are coregulated in various disorders. Similarly, a meta-analysis of gene expression has helped to discover new genes and pathways that are deregulated inactive TB . We conducted gene interaction network studies to investigate the role of differentially expressed genes (DEGs) at a molecular level in the present study, which utilized gene expression data during exposure to antituberculosis medications, bedaquiline. The gene expression omnibus GEO2R tool and R programming were used to perform differential expression analysis. The functional enrichment analysis was performed for the DEGs to find the hub genes and enhanced pathways, which will help researchers, better comprehend the MTB system-wide roles following bedaquiline exposure. The study’s goal was to evaluate and identify possible biomarkers that might aid in the diagnosis of TB from a worldwide TB pandemic era.
2.1. Retrieval of Datasets
The gene dataset was retrieved using GEO database using keywords: TB and bedaquiline and Homosapiens (https://www.ncbi.nlm.nih.gov/gds/?term=tb+and+bedaquiline) . The NCBI-Gene expression omnibus database (NCBI-GEO) is facilitated with the analyses by furnishing several microarray datasets. GSE43749 datasets were chosen for the study as per the inclusion criteria framed.
2.2. Datasets Processing
GEO2R is a web-based interactive tool to compare two different data sets (Normal and Patient group). It is effectively used to pick out DEGs from the given dataset. It will analyze by assigning the group types. For the dataset GSE43749, 251 DEGs were obtained. The results obtained consist of gene ID, gene symbol, gene title, P-value, adj. P-value, and log FC value of DEG (https://www.ncbi.nlm.nih.gov/geo/geo2r/).
2.3. Identification of DEGs
R programming was performed to obtain upregulated and downregulated genes for the given dataset. The cutoff value for screening the DEGs was considered by fixing the log fold change ranges between –0.5 and 0.5 and P ≤ 0.05. The result was obtained using an R programming code consisting of 38 upregulated genes and 37 downregulated genes.
2.4. Functional Annotation and Pathway Enrichment Analysis
The enrichment analysis was performed to explore the biological implications of DEGs. The functional enrichment of DEG in cellular functionality, biological processes. and molecular functions was identified using the PANTHER database. It is a bioinformatics tool to study the biological information for a large set of genes, and functional annotation can be done using an integrated gene ontology database . Pathway enrichment analysis was introduced from genome-scale trials to better understand its features, and it finds biological pathways that are enriched in a gene list provided as input. EnrichR is a comprehensive resource for curated gene sets and a search engine that accumulates biological information for subsequent biological discoveries of DEGs, which is used in this study .
2.5. Comprehensive Analysis of Protein-protein Interactions (PPI) Networks and Modules
After removing redundant genes, the collection of genes generated from the preceding stages was pooled together. PPIs are crucial for understanding how biological systems are regulated. PPI has been found in recent research to identify critical hub genes that cause illnesses and also act as therapeutic targets for the precise treatment of a drug over the corresponding disease. The hub genes mainly play a key role in the potential pathway analysis. The STRING database was used to perform PPI among uniquely discovered DEGs on the whole human genome . The DEGs were mapped to the STRING database with a confidence score of 0.6, and the results were visualized using Cytoscape. The molecular complex detection (MCODE) plugin was used to filter the key network modules with degree cutoff = 10, node score cutoff = 0.2, k core = 2, and max depth = 100 .
2.6. Construction of Target Genes miRNA – Regulatory Network
The potential microRNAs were evaluated for correlations with the candidate genes to uncover possible microRNA-gene correlations. The increased target-gene expression may be linked to increased microRNA expression and vice versa. To create an in-depth relationship, a link between the microRNA profile and gene regulatory networks of both mouse and human by collecting and integrating the documented regulatory interaction was found using the EnrichR tool. The Regnetwork tool was used to understand the details of the transcription factor associated with the human miRNA with the target score and the gene symbol . Finally, the Chea3 tool was used to get the clustering association of the transcription factors involved with the respective miRNA .
3. RESULTS AND DISCUSSION
3.1. Identification of DEGs
To find DEGs, the gene expression dataset related to TB was found by searching in the GEO database using the keywords: TB, bedaquiline, and Homosapiens. The dataset which satisfies our inclusion criteria was GSE43749; it was taken for the subsequent processing steps. The GEO datasets are compared using the GEO2R tool and the DEGs are identified from the whole gene sets [Figure 1]. According to the [log2FC] ±0.5 and P < 0.05, 38 are upregulated genes and 37 are downregulated genes that were observed out of 251 DEG’s by performing R programming.
|Figure 1: Volcano plot for control and patient samples of GSE43749.|
[Click here to view]
3.2. Functional Annotation
We used PANTHER to enrich the procured set of DEGs to get the potential GO categories by classifying them by molecular function, cellular function, biological process, and protein class to determine the importance of the detected DEGs. The identified DEGs were primarily associated with molecular functions divided into the binding, catalytic activity, and transporter activity. It also indicates that both upregulated and downregulated genes have a high amount of catalytic activity, as per Figure 2 and 3. A biological process component of GSE43749 dominantly belongs to the cellular process, metabolic process, localization, and biological response. It was predicted that both cellular and metabolic processes are involved in the essential life processes. The cell envelope of MTB is complex and it is mainly composed of peptidoglycans, mycolic acids, lipids, and carbohydrates. Thus, it is conceivable that the cellular processes of MTB genes belonging to this class could play a key role in mycobacterial intracellular functions and cellular anatomical entry, which can be predominantly used as potential drug targets .
|Figure 2: Functional annotation of upregulated gene set.|
[Click here to view]
|Figure 3: Functional annotation of downregulated gene set.|
[Click here to view]
The protein complex possesses a structural molecular activity in binding the target drug with the genes. The protein classes of GSE43749 include chemokine, cytokine, hydrolase, ribosomal protein, and RNA-binding protein that were enriched along with the signaling and cell adhesion molecule. The curated DEGs were enriched in the gene-specific transcriptional regulator, metabolite inter-conversion enzyme, nucleic acid metabolism protein, and transporters. The genes required for optimal growth of the MTB are identified to be involved in several central metabolic pathways.
3.3. Pathway Enrichment Analysis
The pathway analysis of DEGs was identified using the EnrichR tool. The result so served that Fc epsilon RI signaling pathway enriched for the dataset taken followed by other pathways in Table 1. Fc epsilon RI-mediated signaling pathways in mast cells are initiated by the interaction of antigen (Ag) with immunoglobulin (Ig)E bound to the extracellular domain of the alpha chain of Fc epsilon RI.
Table 1: Major pathways involved in differentially expressed genes.
|Serial number||Name||P||Adjusted P value||OR||Combined score|
|1||Fc epsilon RI signalling pathway||0.08479||0.1954||11.88||29.33|
|2||Fc gamma R-mediated phagocytosis||0.1119||0.1954||8.84||19.36|
|3||Th1 and Th2 cell differentiation||0.1130||0.1954||8.74||19.05|
|4||NF-kappa B signaling pathway||0.1165||0.1954||8.46||18.19|
|5||T-cell receptor signaling pathway||0.1234||0.1954||7.95||16.63|
|6||Th17 cell differentiation||0.1303||0.1954||7.50||15.28|
|7||Natural killer cell-mediated cytotoxicity||0.1572||0.2021||6.11||11.30|
|8||Rap1 signaling pathway||0.2361||0.2618||3.86||5.57|
|9||Ras signaling pathway||0.2618||0.2618||3.42||4.58|
OR: Odds ratio.
Fc epsilon RI-mediated signaling pathways are initiated by the interaction of Ag with antibody (IgE) bound to the extracellular domain of the alpha chain of Fc epsilon RI. The activation pathways are regulated positively and negatively by the interactions of many signaling molecules. Mast cells that are activated thus release granules that contain biogenic amines and proteoglycans. The activation of phospholipase A2 causes the release of membrane lipids, which develops lipid mediators such as leukotrienes (LTC4, LTD4, and LTE4) and prostaglandins (especially PDG2). There is also a secretion of cytokines, the most important of which are TNF-alpha, IL-4, and IL-5. These mediators and cytokines contribute to inflammatory responses .
3.4. Identification and Validation of Hub Genes
We used Cytoscape to execute a PPI network analysis based on the STRING interactome database to investigate the interactive interactions between the DEGs that are represented in Figure 4. The DEGs were formed as two major clusters during Cytoscape MCODE plug-in analysis. One cluster encodes the hub gene dnaA, dnaB, and dnaE1 [Figure 5a], whereas another cluster has atpE, atpF, and atpB [Figure 5b].
|Figure 4: Protein-protein interactions of differentially expressed genes|
[Click here to view]
|Figure 5: (a) Hub genes for upregulated dataset and (b) Hub genes for downregulated dataset.|
[Click here to view]
The dnaA protein interacts with repeated and non-palindromic DNA boxes located inside the oriC region to initiate bacterial chromosomal replication. The associations of dnaA protein with MTB DNA may provide new light on the protein’s function and aid in understanding the control of mycobacterial chromosomal replication start . In Mycobacteria, dnaB is the replicative helicase, which is essential in the replication and growth process. dnaB is also one of the most often seen intein-containing proteins in bacteria . dnaE1 and dnaE2, two DNA polymerases involved in replication in MTB, appear functionally identical. dnaE2 is required for damage-induced base-substitution mutagenesis in MTB. MTB becomes hypersensitive to DNA damage when dnaE2 activity is lost, induced mutagenesis is eliminated, virulence is reduced, and drug resistance is reduced in vivo.
Along with dnaE2 and recA, dnaE1 provides critical and high-fidelity replicative polymerase activity and is expressed in response to DNA damage . The genes atpE, atpF, and atpB belonged to ATP synthase. ATP synthase is a crucial enzyme in practically all living cells’ energy metabolism. In certain species, ATP synthase has distinct traits that might adapt to the habitats they encounter. Pathogenic bacteria may confront unique obstacles regarding ATP generation since they must deal with low oxygen tensions and nutritional scarcity .
3.5. miRNA Prediction
MicroRNAs are indigenous RNAs with roughly 22 bases that target mRNAs for cleavage or post-translational modifications and serve as an essential regulatory function in plants and animals. The miRNA identification is made using the Enrichr tool, and the results are shown to consist of human and mouse miRNA. The four human miRNA (hsa-miR-574-3p, hsa-miR-4787-5p, hsa-miR-601, and hsa-miR-1234) was preceded further for our study in Table 2.
Table 2: Details of human microRNA.
|miRNA||Target score||Target site||Inference|
|hsa-miR-574-3p||92||Retinoid and receptor alpha||The receptors function as transcription factors by to specific sequences in the promoters of target genes|
|hsa-miR-4787-5p||89||Cell migration inducing hyaluronidase2||It binds hyaluronic acid and catalyzes the depolymerization activity|
|has-miR-601||96||BCL2||MTB involves preventing apoptosis induced by up-regulation of Bcl2|
|hsa-miR-1234||87||Innate immunity activator||The development of a tuberculosis drug mainly targets the induction of adaptive immune responses|
MTB: Mycobacterium tuberculosis, miRNA: Micro RNA.
The hsa-miR-574-3p is found in the human chromosome 4 intron regions. A recent study has discovered that changes in miR-574-5p expression are linked to various disorders, notably myocardial infracted, colorectal, and lung cancer. TLR9 signaling may be enhanced by miR-574-5p in lung cancer, promoting tumor growth. MiR-574 expression was elevated in colorectal cancer and myocardial infracted leading to enhanced proliferation and invasion, whereas in this study, the expression of miR-574 was downregulated in response with TB . MicroRNAs and cardiac sarcoplasmic reticulum calcium ATPase-2 in human myocardial infarction: Expression and bioinformatics analysis. The study by Budak et al., 2018 determined that miRNA-4787 showed a similar expression pattern in both acute and chronic groups in interaction with CD8 cells, whereas another study by Pu et al., 2016 reveals that miR- 4787 was consistently expressed in tissues and plasma of patients with various stages of lung carcinoma [22,23].
The hsa-miR-4787-5p and hsa-miR-601 can be used to target gene-related cytokines. The miRNAs are involved in several physiological and pathological processes. It is a positive modulator of TLR signaling and is activated when murine macrophages are infected with Mycobacteria.
3.6. Comprehensive Analysis of miRNA with Transcription Factor
The Regnetwork was used to find possible transcription factors involved in miRNA (hsa-miR-574-3p, hsa-miR-4787-5p, hsa-miR-601, and hsa-miR-1234) that were then evaluated to determine TF coexpression networks using the chea3 tool.
The network analysis for TF with hsa-miR-574-3p was found to have 10 top TF, in which JUND, CAMTA2, and MEF2D were identified to have more interconnections with other TFs shown in Figure 6a. JUND belong to the JUN family of proteins. JUN and other TF can influence immunological activity in the body, and their gene function should be primarily extracellular and protein binding. The expression level JUN will be higher for people affected by SLE peripheral blood mononuclear cells infection .
|Figure 6: Clustering of transcription factors involved with miRNAs, (a) hsa-miR-574-3P, (b) hsa-miR-4787-5p, (c) hsa-miR-601, and (d) hsa-miR-1234-3p.|
[Click here to view]
The miRNA, hsa-miR-4787-5p, showed the higher expression of CEBP and JUND in Figure 6b. An integrated study by Delgado et al., 2019 revealed that TB patients tend to have disease severity with the presence of the CEBP gene in connection with IL6 . The network analysis of a set of TFs for the-miR-601 had shown the central clustering with Early Growth Response Protein-1 (EGR1) and GATA transcription factors in Figure 6c. EGR1 is involved in the control of cell physiology, which affects growth, division, and survivability. EGR1 is widely expressed in tissues and may be quickly activated by various environmental cues, including growth stimulants; shear stress, and reactive oxygen species. The study by Kumar et al., 2020 experimented that knocking out the EGR1 gene leads to prolonged survival of MTB in the host cell. Because MTB survival is strongly related to nitric oxide generation in murine macrophages, we examined nitrite production in infected macrophages . GATA protein family members operate as lineage-specific transcription factors for a range of hematopoietic cell type systems. GATA2 regulates the transcription of genes involved in hematopoietic and endocrine cell lineage formation and proliferation. Mutation in the GATA2 gene leads to the development syndrome and is prone to mycobacterial infections .
The transcription factor analysis of the-miR-1234-3p had shown central clustering of expression levels with NRF2 (NEF2) and GATA1 [Figure 6d]. NEF2, a crucial modulator of the antioxidant defensce system, was often found to protect against various lung illnesses. Human investigations have revealed that NEF2 plays an essential role in protecting against the oxidative damage caused by vigorous smokers. In healthy smokers, NEF2 was stimulated, and several antioxidant enzymes were boosted. The functionality of NEF2 in TB pathogenesis makes this transcription factor a unique entity and notable clinical gene for distinguishing TB patients from normal subjects .
We performed the analysis using gene expression of TB against antituberculosis medications, bedaquiline. The most significant contribution is the hub proteins that might be investigated as potential therapeutic targets and vaccination candidates. In this study, the hub genes identified are dnaA, dnaB, dnaE1, atpE, atpF, and atpB. During host infection, the increase of DNA gene expression directly contributes to mutation leads to DNA damage induced by immune-mediated reactive oxygen and nitrogen intermediates. Especially, ATP genes may play an essential role in the emergence of mutants that are better suited to surviving during infection and drug resistance evolution in this organism . The bacteria tend to reduce cellular ATP consumption while increasing the capacity of ATP-generating pathways, which adds to bacterial survival in the face of antibiotic stress. In Mycobacteria, ATP synthase is required for growth, and the drug Bedaquiline effectively inhibits the operation of this critical metabolic enzyme. As a result, inhibitors of respiratory ATP production may be beneficial for eradicating this difficult-to-kill MTB subpopulation in human host cells .
This study identified the DEGs for the dataset GSE43749 and their associated biological processes, interactions, and pathways among healthy controls and active TB groups. The four miRNAs (hsa-miR-574-3p, hsa-miR-4787-5p, hsa-miR-601, and hsa-miR-1234) were determined to have a higher expression when susceptible to TB infection and primarily improve immunological responses. JUND, EGR1, NEF2, and GATA1 transcription factors were explored as a potential target gene. The observed hub genes seemed to have the most significant impact on ATP generation, DNA replication, immunological activity, and oxidative damage. Taken as a whole, the scientific community’s efforts are expected to usher in a new age of diagnostic tests based on hub genes – miRNAs and TF can be offered as prospective biomarkers. These diagnostic techniques must not only meet quality standards for specificity and sensitivity, but they must also be biologically relevant in the pathophysiology of TB.
5. AUTHORS’ CONTRIBUTIONS
V. Anusuya and P. Madhan carried out the experiments and prepared the original manuscript. B. Nivetha helped in editing and drafting the manuscript. K Santhiya: Supervision, conceptualization, methodology, final review and approval.
There is no funding to report.
7. CONFLICTS OF INTEREST
The authors report no financial or any other conflicts of interest in this work.
8. ETHICAL APPROVALS
This study does not involve experiments on animals or human subjects.
9. DATA AVAILABILITY
All data generated and analyzed are included within this research article.
10. PUBLISHER’S NOTE
This journal remains neutral with regard to jurisdictional claims in published institutional affiliation.
1. WHO global lists of high burden countries for tuberculosis (TB), TB/HIV and multidrug/rifampicin-resistant TB (MDR/RR-TB), 2021–2025. Available from:https://cdn.who.int/media/docs/default-source/hq-tuberculosis/who_globalhbcliststb_2021-2025_backgrounddocument.pdf?sfvrsn=f6b854c2_9 [Last accessed on 2021 Dec 24].
2. Balashanmugam MV, Shivanandappa TB, Nagarethinam S, Vastrad B, Vastrad C. Analysis of differentially expressed genes in coronary artery disease by integrated microarray analysis. Biomolecules 2019;10:E35. [CrossRef]
3. World Health Organization. WHO Consolidated Guidelines on Drug-Resistant Tuberculosis Treatment. Geneva:World Health Organization;2019.
4. Wu LS, Lee SW, Huang KY, Lee TY, Hsu PW, Weng JT. Systematic expression profiling analysis identifies specific microRNA-gene interactions that may differentiate between active and latent tuberculosis infection. Biomed Res Int 2014;2014:895179. [CrossRef]
9. Mi H, Ebert D, Muruganujan A, Mills C, Albou LP, Mushayamaha T, et al. PANTHER version 16:A revised family classification, tree-based classification tool, enhancer regions and extensive API. Nucleic Acids Res 2021;49:D394-403. [CrossRef]
11. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11:Protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res 2019;47:D607-13. [CrossRef]
14. Keenan AB, Torre D, Lachmann A, Leong AK, Wojciechowicz ML, Utti V, et al. ChEA3:Transcription factor enrichment analysis by orthogonal omics integration. Nucleic Acids Res 2019;47:W212-24. [CrossRef]
15. Kalscheuer R, Palacios A, Anso I, Cifuente J, Anguita J, Jacobs WR Jr., et al. The Mycobacterium tuberculosis capsule:A cell structure with key implications in pathogenesis. Biochem J 2019;476:1995-2016. [CrossRef]
17. Zawilak A, Kois A, Konopa G, Smulczyk-Krawczyszyn A, Zakrzewska-Czerwi?ska J. Mycobacterium tuberculosis DnaA initiator protein:Purification and DNA-binding requirements. Biochem J 2004;382:247-52. [CrossRef]
19. McGuire AM, Weiner B, Park ST, Wapinski I, Raman S, Dolganov G, et al. Comparative analysis of Mycobacterium and related Actinomycetes yields insight into the evolution of Mycobacterium tuberculosis pathogenesis. BMC Genomics 2012;13:120. [CrossRef]
21. Lai Z, Lin P, Weng X, Su J, Chen Y, He Y, et al. MicroRNA-574-5p promotes cell growth of vascular smooth muscle cells in the progression of coronary artery disease. Biomed Pharmacother 2018;97:162-7. [CrossRef]
22. Pu Q, Huang Y, Lu Y, Peng Y, Zhang J, Feng G, et al. Tissue-specific and plasma microRNA profiles could be promising biomarkers of histological classification and TNM stage in non-small cell lung cancer. Thorac Cancer 2016;7:348-54. [CrossRef]
25. Delgobo M, Mendes DA, Kozlova E, Rocha EL, Rodrigues-Luiz GF, Mascarin L, et al. An evolutionary recent IFN/IL-6/CEBP axis is linked to monocyte expansion and tuberculosis severity in humans. Elife 2019;8:e47013. [CrossRef]
26. Kumar M, Majumder D, Mal S, Chakraborty S, Gupta P, Jana K, et al. Activating transcription factor 3 modulates the macrophage immune response to Mycobacterium tuberculosis infection via reciprocal regulation of inflammatory genes and lipid body formation. Cell Microbiol 2020;22:e13142. [CrossRef]
27. Mendes-de-Almeida DP, Andrade FG, Borges G, Dos Santos-Bueno FV, Vieira IF, da Rocha LK, et al. GATA2 mutation in long stand Mycobacterium kansasii infection, myelodysplasia and MonoMAC syndrome:A case-report. BMC Med Genet 2019;20:64. [CrossRef]
29. Chengalroyen MD, Mason MK, Borsellini A, Tassoni R, Abrahams GL, Lynch S, et al. DNA-Dependent Binding of Nargenicin to DnaE1 Inhibits Replication in Mycobacterium tuberculosis. ACS Infect Dis. 2022;612-25. [CrossRef]
30. Koul A, Vranckx L, Dhar N, Göhlmann HW, Özdemir E, Neefs JM, et al. Delayed bactericidal response of Mycobacterium tuberculosis to bedaquiline involves remodelling of bacterial metabolism. Nat Commun 2014;5:3369. [CrossRef]