Introduction
Tuberculosis (TB) is an infectious disease caused by the bacillus Mycobacterium tuberculosis (M. tb), one of the most successful bacteria that infects humans (1). TB has existed for millennia and remains a major global health problem (2). In recent years, the emergence and spread of (multi)drug-resistant M. tb along with the limited protection afforded by the current BCG vaccine in a number of populations has made it difficult to effectively control the TB epidemic. TB is one of the top 10 causes of mortality and the leading cause of mortality from a single infectious agent, more than HIV/AIDS, and millions of individuals continue to become infected each year (3). The clinical infection of TB is complicated. The aim of the present study was to reveal the general host response mechanism to TB infection and examine its potential immune mechanism effects through the analysis of the M. tb H37Rv strain.
Macrophages, being professional phagocytes, ingest and degrade dead cells, debris and foreign material, and orchestrate inflammatory processes, which serve an essential role in the control of bacterial infections (4,5). When M. tb infects the host, it first needs to resist the immune response from macrophages (6,7). Paradoxically, macrophages may be the primary long-term environment for M. tb in the host, as M. tb is able to evade the principle antimicrobial mechanism of macrophages through a variety of bacterial immune subversive mechanisms, allowing the bacteria to utilize various intracellular cell resources (8). However, the mechanisms by which M. tb enters the host cell, circumvents host defenses and spreads to neighboring cell are not completely understood. Therefore, an improved understanding of host-pathogen interaction is essential if the global TB pandemic is to be controlled (9). The interaction between macrophages and M. tb alters cellular gene expression profiles, and studies of changes in expression of these genes may assist to understand the immunomodulatory effects of M. tb in macrophages.
With the beginning of the era of big data, the mining and application of high-throughput big data for genes has received more and more attention. Microarray technology is fundamental to generating large amounts of gene expression data. Microarray technology is widely used as it is able to simultaneously detect and express high-throughput expression levels of 10,000's of gene transcription sets (10). Microarray technology can be used to study gene expression more efficiently and in-depth. With the increasing availability of microarray technology, gene expression data associated with various physiological and disease conditions are provided on public platforms, including the Gene Expression Omnibus (GEO) database and The Cancer Genome Atlas database (11). The meta-analysis tools for the systematic and comprehensive quantitative analysis of the results of multiple independent studies investigating the same disease have been developed (12). A number of studies have used meta-analysis methods to analyze small sample experimental studies, improving the accuracy of statistical performance and effect estimation (13-19).
The present study used meta-analysis methods to identify overlapping DEGs in THP-1-derived macrophages infected by M. tb, analyze protein-protein interaction (PPI) networks, identify potential central genes, complete functional and pathway enrichment studies and verify pivotal gene differences expression through experiments. The potential effects of the immune mechanism of TB infection was estimated using the aforementioned method.
Materials and methods
Microarray data collection and quality control
The present study used the GEO database at the National Center for Biotechnology Information (http://www.ncbi.nlm.nih.gov/geo/) to retrieve expression profile datasets. The search terms used were: 'Mycobacterium tuberculosis', 'THP-1' and 'Homo sapiens'. Series Matrix Files of these datasets was downloaded respectively, and the probe IDs were transformed into official gene symbols. In order to obtain more credible results, the MetaQC package in R was applied to implement quality control analysis and directly used the principal component analysis chart to display statistical results. There are 6 quantitative quality control measures in this package, including: i) internal homogeneity of co-expression structure among studies; ii) external consistency of co-expression structure correlating with a pathway database; iii) accuracy of differentially expressed genes (DEG) detection and pathway identification; and iv) consistency of differential expression ranking of genes and pathways. Through quality control, only the data with high quality was included in the subsequent analysis.
DEGs identification
DEGs were screened out with MetaDE v.1.0.5 software package (20), which provides 12 major meta-analysis methods. In the present study, 3 methods, including Fisher, rth-ordered P-value (roP) and random effects model (REM) were selected to filtrate DEGs in high quality data. A false discovery rate (FDR) <0.05 was considered as the cut-off criterion. In addition, the intersection of the required datasets under the three statistical methods was extracted and considered as the most reliable DEGs. DEGs were regressed into high quality data to complete cluster analysis.
PPI network construction
The initial PPI network for the protein products of identified DEGs was constructed using the Search Tool for the Retrieval of Interacting Genes/Proteins Database (STRING v10; https://string-db.org/) (21), and then the network was visualized and analyzed with Cytoscape v.3.6.1 software (22). To extract the characteristic genes from the DEGs, the CentiScape and MCODE plugins in the Cytoscape software were used. Firstly, the Degree Centrality (DC) of CentiScape reflects hubness of the node in the PPI network by calculating the degree. Nodes with higher centrality were extracted. Secondly, constructing sub-networks through k-core analysis of the MCODE afforded the selection a stable structure from a high-centrality network, and narrowed the scope of the study.
Function and pathway enrichment analysis of characteristic genes
To assess the prospective functions of characteristic genes identified by the meta-analysis and PPI network analysis, Gene Ontology (GO) biological pathways and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways were identified using CluGO from Cytoscape software. P<0.05, subjected to Bonferroni adjustment, were defined as the cut-off criterion.
M. tb H37Rv and THP-1 cell line propagation and infection
The M. tb H37Rv strain (American Type Culture Collection; product no., 27294) was grown at 35°C for 15 days in a Middlebrook 7H9 medium (BD Biosciences), supplemented with 10% (v/v) OADC (BD Biosciences), with the addition of 0.2% glycerol (v/v) and 0.025% Tween-80 (v/v). The human acute monocytic leukemia THP-1 cell line (American Type Culture Collection; cat. no. TIB-202) was cultivated at 37°C in 5% CO2 in RPMI-1640 (Hyclone; GE Healthcare Life Sciences) medium supplemented with 10% FBS (Gibco; Thermo Fisher Scientific, Inc.). Prior to infection, THP-1 cells were treated with 200 nM phorbol-12-myristate-13-acetate (Sigma-Aldrich; Merck KGaA) for 24 h to allow for the differentiation into macrophages in 6-well plates. Infections were performed with the M. tb H37Rv strain at a multiplicity of infection of 5. The 7H9 broth was used as a blank control.
M. tb infection confocal microscopy
The success of the M. tb infection model was verified by observing the co-localization of M. tb with THP-1 cell lysosomes. Firstly, M. tb was labeled with fluorescent isothiocyanate (FITC) solution (Sigma-Aldrich; Merck KGaA) according to the manufacturer's protocol. M. tb was washed 3 times with PBS buffer, resuspended in carbonate buffer, and added to FITC (10 mg/ml) to a final concentration of 100 µg/ml. Following incubation for 2 h at room temperature, the supernatant was isolated by centrifugation at 1,000 × g for 2 min at 4°C, and inoculated into THP-1 cells previously induced to adhere to slides in a 6-well plate, and the whole procedure was performed in the dark. Then, according to the manufacturer's protocol, the Lyso Tracker Red99 (Sigma-Aldrich; Merck KGaA) solution was prepared with RPMI-1640 medium, mixed and stored at 4°C, and preheated at 37°C prior to use. After 24 h, the cell culture solution of infection was discarded, washed 3 times with PBS, and 1 ml Lyso Tracker Red99 solution was added to each well to label the lysosomal component of the macrophage, and incubated at room temperature for 30 min. Finally, the supernatant was directly discarded using a pipette and the slides were washed 3 times with PBS. Then, the slides were removed, fixed on a glass slide with 4% paraformaldehyde for 20 min at room temperature, observed under a confocal microscope (magnification, ×20) and visualized by FV10-ASW 2.1 Viewer software (Olympus Europa Holding GmbH) to observe the co-localization.
Reverse transcription quantitative polymerase chain reaction (RT-qPCR)
To validate the accuracy of DEGs microarray analysis results, the total RNA of cells from the infection model was extracted using Trizol® (Invitrogen; Thermo Fisher Scientific, Inc.) and 2 µg RNA was transcribed into cDNA by PrimeScript™ RT reagent kit with gDNA Eraser (Perfect Real Time) (Takara Biotechnology Co., Ltd.). Relative expression levels of mRNA were quantitatively normalized and analyzed by RT-qPCR using FastStart Universal SYBR Green Master (ROX; Roche Diagnostics) against the expression of β-actin, using the 2−ΔΔCq method (23). The primer sequences were derived from PrimerBank (https://pga.mgh.harvard.edu/prim-erbank/) (24) (Table SI). The thermocycling conditions for RT-qPCR are presented in Table SII.
Statistical analysis
Receiver operating characteristic (ROC) curves were used to analyze the differentially expressed genes [interferon (IFN)-stimulated gene 15 (ISG15), 2′-5-oligoad-enylate synthetase like (OASL), IFN regulatory factor 7 (IRF7) and DExD/H-box helicase 58 (DDX58)] between TB infected group and control group. Quantitative data were expressed as the mean ± standard deviation. GraphPad Prism 6 software (GraphPad Software, Inc.) was utilized to draw the ROC curves and calculated the sensitivity, specificity and area under the curve (AUC). The AUC of the ROC curve >0.7 was regarded as significant. Data from multiple groups were analyzed using analysis of variance followed by a Student-Newman-Keuls post-hoc test. P<0.05 was considered to indicate a statistically significant difference.
Results
Identification of DEGs with microarray meta-analysis
A total of 6 microarray datasets were downloaded from the GEO using the aforementioned search method (Table I). These datasets were analyzed by MetaQC package and visualized by principle component analysis (PCA; Table II; Fig. 1). The results of the comprehensive PCA indicate that the datasets GSE17477, GSE51029, GSE57028 and GSE89391 were rated as high quality data. The raw data of these 4 high quality datasets were analyzed using the R package MetaDE, for a total of 76 independent samples (30 H37Rv-infected samples and 46 controls). The shared genes from the 4 datasets were analyzed by 3 methods: Fisher; roP; and REM (Fig. 2A-C). During the analysis, the non-DEGs, and those exhibiting small variation and expression intensities across the majority of datasets were filtered out. Through the same three statistical methods (Fisher, roP and REM), 1,190, 426 and 927 differentially expressed genes were identified, respectively. The intersection of the 3 datasets yielded 306 common DEGs (Fig. 3). These 306 DEGs were arranged into 4 high-quality datasets for hierarchical cluster analysis (Fig. 4A-D). It was clear that these genes were well-differentiated between the M. tb infection and control groups. The top 20 shared DEGs identified in the meta-analysis are summarized in Table III, ranked by average log2 fold change.
Prioritization of DEGs by PPI network analysis
The PPI network of these DEGs was established using STRING database, and the node and edge relationships were imported into Cytoscape software for visualization, and then the DC evaluation was completed. The degree result is demonstrated by the size of the point and the thickness of the edge (Fig. 5A). The first 56 genes were selected as hub genes (Fig. 5B) as a focus for subsequent examination. Using the k-core decomposition function of MCODE with average MCODE score =28.774, nodes =32 and edges =446, a variety of sub-networks were obtained and a new network was created (Fig. 5C). This resulted in 32 remaining differentially expressed genes in the new network. These 32 differentially expressed genes were therefore selected for further research to narrow the scope of the present study. The association between the three PPI networks is demonstrated in Fig. S1. The results of DC and MCODE of DEGs analysis are summarized in Table IV.
Function and pathway enrichment analyses of hub genes
The GO and KEGG term enrichment analyses were performed using CluGO. There were a total of 35 GO terms within the 'Biological Process' (BP) group and 3 KEGG pathways with P<0.05 (including Bonferroni adjustment) identified within the analysis of the 32 characteristic genes.
Among the top significant GO BP terms identified were those associated with response to IFN-α and IFN-β, the IFN-γ-mediated signaling pathway and the type I IFN signaling pathway (Fig. 6A). Only 3 significant KEGG pathways, the retinoic acid-inducible gene-I (RIG-I)-like receptor signaling pathway, cytosolic DNA-sensing pathway and influenza A were identified (Fig. 6B). Detailed enrichment information for the GO and KEGG pathways is summarized in Tables V and VI. Based on the results of the function and pathway enrichment analyses, the result that M. tb infection was closely associated with virus-associated pathways, including the type I IFN signaling pathway, is consistent with the characteristics of M. tb as an intracellular pathogen. Among these pathways, it was hypothesized that genes which serve a significant role include bone marrow stromal cell antigen 2 (BST2), IRF7, IFN induced protein with tetratricopeptide repeats 1 (IFIT1), C-X-C motif chemokine ligand 10 (CXCL10), ISG15, DDX58 and OASL.
M. tb successfully enters THP-1-derived macrophage lysosomes
Co-localization of the lysosome specific-acido-tropic dye Lyso Tracker Red99 with live virulent M. tb in macrophages was examined by confocal microscopy (Fig. 7). A total of 2 sets of images were recorded: Firstly, the localization of the FITC-labeled Mycobacterium-positive compartment (green) and secondly, the Lyso Tracker Red99 lysosome-positive compartment (red). Following merging of the 2 image sets, the sections in which M. tb was co-localized with the lysosomes appeared yellow. Finally, a three-dimensional scan of THP-1 derived macrophages was obtained. These results demonstrated that THP-1 cells were successfully cultured and induced into macrophages, and that an M. tb infection model was established to verify the reliability of the previous bioinformatics analysis.
Analysis of relative expression levels of ISG15, OASL, IRF7 and DDX58 in the infection model
A total of 4 genes, ISG15, OASL, IRF7 and DDX58, were selected from the important genes identified, and their differential expression levels between the infected and control groups was verified by RT-qPCR. Using the 2−ΔΔCq method, the difference in expression levels of the 4 genes relative to β-actin was calculated, and it was identified that the expression of these 4 genes was significantly upregulated following M. tb infection. Then, GraphPad Prism 6 software was used to visually demonstrate the comparison between RT-qPCR results and microarray meta-analysis results, and it was identified that the experimental verification results were consistent with bioinformatics analysis data (Fig. 8). These genes were differentially expressed following M. tb infection of the host and are associated with the response mechanism of the host to M. tb, which requires additional study.
ROC curve analysis of the expression of ISG15, OASL, IRF7 and DDX58
The ROC curve analysis of ISG15, OASL, IRF7 and DDX58 identified that the AUC values were 0.9691, 0.9228, 0.8488 and 0.679, respectively (Fig. 9). Among the 4 genes, only the AUC of the ROC curve of DDX58 was <0.7, indicating that DDX58 did not have good specificity and sensitivity, while ISG15, OASL and IRF7 were statistically significant, and may be able to distinguish between infected and uninfected cells.
Discussion
In the present study, a meta-analysis was conducted using microarray-based gene expression studies that were previously published and are available for public reuse. Under the premise of reducing research bias and individual heterogeneity, a total of 306 DEGs with high credibility were identified, including 178 upregulated and 128 downregulated genes with FDR<0.05. Based on these DEGs, a PPI network was constructed, including 275 nodes and 2038 edges, and 32 characteristic genes were obtained through DC sorting and k-core analysis, including ISG15, OASL, IRF7, DDX58, IFN induced with helicase C domain 1, IFIT1, IFN induced protein with tetratricopeptide repeats 3, HECT and RLD domain containing E3 ubiquitin protein ligase 5, IFN regulatory factor 9, MX dynamin like GTPase 2, IFN induced protein with tetratricopeptide repeats 2, radical S-adenosyl methionine domain containing 2, ubiquitin conjugating enzyme E2 L6, IFN induced protein 35, XIAP associated factor 1, HECT and RLD domain containing E3 ubiquitin protein ligase family member 6, IFN induced protein 44 (IFI44), IFI44 like, DExD/H-box helicase 60, receptor transporter protein 4, IFN stimulated exonuclease gene 20, ubiquitin specific peptidase 18 and sterile α motif domain containing 9.
A THP-1-derived macrophage model of M. tb infection was constructed, and the co-localization association between M. tb and lysosome was confirmed by confocal microscopy to determine the success of infection. Among the 32 characteristic genes identified, the genes with the top 4 DC values in the PPI network and associations with type I IFN were selected as the verification targets. The expression of these 4 genes extracted by bioinformatics analysis was verified, and the RT-qPCR results were consistent with the microarray data analysis, which suggested that the results of the analyses were reliable.
These characteristic genes were primarily enriched in BP pathways including responses to IFN-α and IFN-β, the IFN-γ-mediated signaling pathway, the type I IFN signaling pathway, defense response to virus and regulation of viral genome replication. Analysis of the KEGG pathways revealed that the characteristic genes were involved in the RIG-I-like receptor signaling pathway, cytosolic DNA-sensing pathway, and influenza A pathways. These results are consistent with the intracellular parasite characteristics of M. tb. This species is known to possess a specific secretion system, early secretory antigenic target-6 secretion system 1 (ESX-1), which serves a key role in mycobacterial virulence by promoting the release of M. tb products into the cytoplasm of macrophages (25-28). We hypothesized that the double-stranded DNA of M. tb enters the host cell via the ESX-1 system, stimulating the response recognition of the cytosolic DNA-sensing pathway, which triggers downstream type I IFN-associated signaling pathways.
The type I IFN system, as the first important line of defense against viral infection, detects viruses by pattern-recognition receptor recognition of pathogen-associated molecular patterns, and then expresses >300 IFN-stimulated genes to mediate multiple antiviral mechanisms to exert antiviral and immunomodulatory effects (29). However, viruses have developed an impressive array of strategies to circumvent IFN responses, for example, by disrupting Toll-like receptors or RIG-I-like receptors to prevent initial virus detection, and the degradation of transcription factors involved in IFN expression, in order to prevent IFN production (30). The interaction between viruses and the type I IFN system requires additional study. The present study only focused on the role of the type I IFN system in M. tb infection, as it has been suggested previously that the type I IFN pathway may promote the survival and replication of M. tb in host cells, and that this is a common strategy of a number of intracellular pathogenic bacteria to promote virulence (28,31-38). Understanding the role of type I IFN signaling is critical to understanding the pathology and outcome of TB. Therefore, the present study focused on the most important first-line genes identified in the PPI network, which were associated with type I IFNs and validated by biological experiments and ROC curve analysis, ISG15, OASL and IRF7.
ISG15 is one of the most commonly induced proteins by type I IFN (39) and during M. tb infection (40). It is a ubiquitin-like protein that is conjugated to intracellular target proteins upon activation by IFN-α and IFN-β (41). It has been demonstrated to have several functions, including promoting T cell-dependent natural killer cell proliferation, augmenting lymphokine activated killer activity, serving as a chemotactic factor that enhances neutrophil chemotactic activity, directing ligated target proteins to intermediate filaments, cell-to-cell signaling and antiviral activity during viral infections (42-45).
It was suggested that type I IFN and ISG15 serve transient roles in promoting bacterial replication. However, as the disease progresses, ISGylation deviates from the overall effect of type I IFN and, ultimately, mice deficient in ISGylation are significantly more susceptible than IFNAR mice (46). In the study conducted by Keller et al (47), differences in gene expression in M. tb-infected susceptible and resistant mouse strains were investigated, and ISG15 expression was also significantly increased following infection with M. tb for 28 days. However, to the best of our knowledge, no studies have been conducted in patients with TB; ISG15 may lead to different disease outcomes during the development of TB infection in humans, and the actual role requires additional study.
The data from the study by Sambarey et al (48), involving whole blood samples from patients with TB, indicated that OASL was significantly upregulated in blood samples from patients with active TB compared with patients with latent TB, and the expression level in the blood samples of patients after 6 months of treatment was significantly downregulated. However, to the best of our knowledge, no in-depth study of the role of OASL in the infection process has been performed as of yet. OASL belongs to the type I IFN-induced antiviral protein OAS family consisting of OAS1, OAS2, OAS3 and OASL (49). OASL is associated with the OAS family through the N-terminal OAS-like domain but lacks synthetase activity (50). Although OASL lacks enzymatic activity, recently identified functions appear to be critical for viral innate and adaptive immune responses (51). In response to viral infections, OASL serves 2 distinct roles. In the case of RNA virus infection, OASL may enhance the ability of type I IFN signaling to exert antiviral effects (51). However, when challenged with DNA viruses, OASL exerts a pro-survival function by inhibiting type I IFN signaling (52,53). Nonetheless, its biological role in microbial infection is not clear. As M. tb is phagocytosed by macrophages, it is able to release contents including dsDNA into host cell solute, which is similar to DNA virus attack. A recent study identified that when the host was infected with M. leprae, OASL was able to inhibit antimicrobial peptides expression, impair the mechanism of bacterial kill activation and inhibit the autophagy process of host cells, thereby preventing microbial clearance (54). In view of these data and the similarity between the physiological state and infection process of M. tb and M. leprae, it is reasonable to conclude that OASL may serve a similar role in M. tb infection, relying on the type I IFN pathway to inhibit autophagy of macrophages, provide sites and raw materials for M. tb, and promote the progression of TB.
Finally, IRF7 encodes IFN regulatory factor 7, a member of the IFN regulatory transcription factor family. IRF7 was also demonstrated to be upregulated in a host RNA sequencing study in response to M. tb infection (55). Cheng and Schorey (56) indicated that M. tb may rely on cytoplasmic DNA and RNA sensing pathways to mediate IFN-βproduction by IRF7, thereby assisting the host in the identification of bacterial infection. In the study by Lin et al (57), involving whole blood gene expression profiles in patients with active and latent TB, a transcription factor gene regulatory network was constructed and IRF7 was revealed to be significantly upregulated and enriched in the immune and infection pathways, which was hypothesized to help patients resist TB infection (57).
In conclusion, the present study demonstrated that the type I interferon signaling pathway serves a significant role in the infection of macrophages by M. tb. In addition, the mechanism of a number IFN-stimulated genes in TB infection remains uncertain, but ISG15 and OASL are likely to become a potential target for TB treatment, creating new possibilities for TB treatment. However, the results from the present study were obtained via bioinformatics analysis, and require additional confirmation using animal model studies.