Identification and Validation of Genes Related to Macrophage Polarization and Cell Death Modes Under Mycobacterium tuberculosis Infection

Introduction

Tuberculosis (TB), caused by Mycobacterium tuberculosis (Mtb) infection, is a serious infectious disease. Until the coronavirus disease 2019 (COVID-19) pandemic, TB has been the leading cause of death from a single infectious agent worldwide. According to the Global TB Report 2022, World Health Organization (WHO) estimates that about a quarter of the world’s population has been infected with Mtb.¹ Between 2020 and 2021, the burden of drug-resistant TB increased by 3% to 450,000 cases.

During Mtb infection, macrophages (Mφ) are the main target and primary barrier, as well as the main site of immune evasion against innate immunity by Mtb. Mφ can be polarized into multiple states after infection, with M1 (typically activated Mφ) or M2 (alternatively activated Mφ) being predominant.² M1 Mφ exhibit Mtb-killing and pro-inflammatory effects through high expression of pro-inflammatory and antimicrobial molecules such as interferon-β (IFN-β), interleukin-1β (IL-1β), and tumor necrosis factor (TNF).³ In contrast, M2 Mφ control body inflammation and supports Mtb growth by expressing high levels of anti-inflammatory cytokines such as IL-4, IL-10, and chemokine (CCL22).^4,5 In the early stage of Mtb infection, the polarization state of Mφ is mainly M1-type, while in the late stage of infection, M2-associated markers gradually increase, and the Mφ polarization is dominated by M2-type.⁶ However, up to now, the regression of Mφ subtypes under Mtb infection is still not determined.

Mtb establishes an in vivo infection by parasitizing Mφ, spreads between cells and organs, and ultimately leads to host cell death.^7,8 Manipulation of the Mφ death mode in infected macrophages is a potential mechanism for host defense against Mtb.⁹ However, it is also one of the mechanisms by which Mtb evades host defense.¹⁰ Mtb-infected Mφ exhibit multiple modes of death, such as apoptosis, autophagy, ferroptosis, and pyroptosis.¹¹ Apoptosis is the first definitive mode of programmed cell death that reduces intracellular Mtb activity without causing any inflammatory responses. CASP3, as a key enzyme controlling apoptosis, marks the onset of apoptosis.¹² Autophagy promotes the maturation of Mφ phagosomes, a defense mechanism by which Mφ removes Mtb and inhibits Mtb growth and reproduction. LC3, as a structural protein of the autophagosome, is an indicator of the occurrence of autophagy.^13,14 Degradation of ferritin FTH1 increases intracellular iron levels, which triggers ferroptosis.¹⁵ The study by Eduardo et al identified ferroptosis as a major mechanism of Mtb-induced cell death in vitro and highly correlated with the course of Mtb infection.¹⁶ Pyroptosis is a lytic cell death mode that is dependent on CASP1, GSDMD, and NLRP3 inflammasomes.¹⁷ Mtb infection induces an inflammatory response in Mφ, which leads to pyroptosis.¹⁸ However, the correlation between the Mtb-infected Mφ death modes and Mφ polarization has not been illustrated, which greatly affects the understanding of immune regulatory mechanisms under Mtb infection and the immune escape mechanisms of Mtb. Moreover, exploration of biomarkers involved in Mtb-infected Mφ death modes and polarization would also benefit the diagnosis and treatment of TB.

This study screened genes that intersected Mφ polarization with different cell death modes under Mtb infection. Next, Hub genes were identified by constructing PPI networks and ROC curve analyses. Subsequently, we detected the marker genes’ expression of Mtb-infected M1/M2 Mφ death modes in vitro experiments.

Materials and Methods

Data Sources

Original gene expression profiles (GSE5099, GSE52819, GSE54992) were collected from GEO database. Apoptosis, autophagy, ferroptosis, and pyroptosis related genes were collected from the KEGG, FerrDb and GSEA databases. The workflow of the present study is shown in Figure 1.

Figure 1 Workflow chart of the bioinformatics analysis.

Differential Expression of Death Related Genes (DEDRGs) Analysis

The data of GSE5099, GSE52819 and GSE54992 were analyzed using the “limma” package of R, and the differentially expressed genes (DEGs) were filtered with the thresholds of adj.P.Val<0.05, |logFC|>1. The volcano maps, venn maps, and Pheatmap were plotted using the “ggplot2”, “VennDiagram”, and “Pheatmap” packages of R, respectively.

GO and KEGG Pathway Enrichment Analysis of DEDRGs

Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses were performed using the “clusterProfiler” package with the standard setting of P<0.05. The GO analysis consisted of cellular component (CC), biological process (BP), and molecular function (MF).

Protein–Protein Interactions (PPI) Analysis of DEDRGs

The STRING database (http://string-db.org) was used to build the PPI network. The Hub genes with the highest degree were analyzed by Hubba plugin of Cytoscape (3.9.1).

Receiver Operating Characteristic (ROC) Curve Analysis

The ROC curve analysis was performed to calculate the area under the curve (AUC) and evaluate the diagnostic value of the Hub genes using gene expression data in GSE54992 with the “pROC” package. Genes with AUC >80% was selected as the candidate gene.

Genetic Correlation Analysis

The dataset GSE54992 was used to map the correlation network of Hub genes with the “corrr” package.

Immune Infiltration Analysis

The proportion of 22 immune cells in the samples of the GSE54992 dataset was obtained using the CIBERSORTx online website (https://cibersortx.stanford.edu/). Box plots and heat maps depicting the level of immune cell infiltration between active TB patients and healthy controls were plotted using the “ggplot2” package.

Macrophage Polarization

THP-1 cells were purchased from ATCC (USA, TIB-202) and cultured using RPM1 1640 and 10% fetal bovine serum (FBS). Logarithmic growth phase cells were taken and inoculated into six-well plates at a density of 1.0×10⁶ per well and treated with 100 ng/mL PMA (SIGMA) for 24 hours to activate to M0 Mφ. After that, M1 Mφ were generated by incubation with 150ng/mL LPS (SIGMA) and 30ng/mL IFN-γ (Sinobiological) for 48 hours. And M2 Mφ were generated by incubation with 30ng/mL IL-4 (Sinobiological) for 48 hours. This study was conducted in accordance with the Helsinki Declaration and was approved by the Institutional Review Ethics Committee of the Affiliated Hospital of Guangdong Medical University.

H37Ra Infects Mφ

H37Ra was cultured and proliferated in 7H9 broth containing OADC (1:10) and ampicillin (1:1000). H37Ra (MOI = 10) was inoculated into the above cultured Mφ and infected for 4 hours, followed by three PBS washes to remove extracellular bacteria for subsequent experiments.

Quantitative Real-Time Polymerase Chain Reaction (qRT-PCR)

Total RNA was extracted from Mφ using RNAiso Plus reagent (TakaRa). mRNA was reverse transcribed to cDNA using the PrimeScript RT kit (TakaRa). qPCR was performed in QuantStudio 6 Flex using SYBR premix Ex Taq II (TakaRa). Relative expression of the genes was calculated using GAPDH as an internal reference gene by 2^−ΔCt [2 ^{(−ΔCt sample)} / 2^{(−ΔCt of control)}] method. The qRT-PCR oligonucleotide primers are shown in Table 1.

Table 1 Primer Sequences Used for qRT-PCR

Western Blot

Mφ was lysed in RIPA lysis buffer (Beyotime). Protein concentrations were measured using the BCA protein assay kit (Solarbio). GAPDH (proteintech, 10494-1-AP), Tubulin (Affinity Biosciences, DF7967), IL-1β (proteintech, 16806-1-AP), IL-10 (abcam, ab133575), CASP3 (abcam, ab13847), CASP1 (abcam, ab179515), LC3B (abcam, ab225382), FTH1 (Cell Signaling Technology, 4393) antibodies were incubated at 4°C overnight. Further incubation with HRP-linked anti-rabbit IgG Antibody (Beyotime, A0208) was performed for 1 hour at room temperature. Protein bands were visualized by an enhanced chemiluminescence detection kit (Millipore) and quantified using ImageJ software.

Statistical Analysis

All experiments were repeated at least three times. The results were analyzed using GraphPad Prism 9.0 software. Multiple sample data comparison was analyzed using ordinary one-way ANOV. For two sample data, a paired t-test was used to determine P values. Statistical significance was defined as P<0.05. P values were assigned as *P<0.05; **P<0.01; ***P<0.001; and ****P<0.0001.

Results

Identification of DEDRGs Between Normal M0 Mφ and Polarized Mφ Infected with Mtb

We analyzed 3 cases of M0 Mφ, 3 cases of M1 Mφ, and 3 cases of M2 Mφ samples from GSE5099, and identified 2584 DEGs (1535 up-regulated genes and 1049 down-regulated genes) for M1 Mφ and 394 DEGs (327 up-regulated genes and 67 down-regulated genes) for M2 Mφ. In addition, we analyzed 3 cases of Mtb-infected Mφ and 3 cases of untreated Mφ samples from GSE52819, and identified 6417 DEGs (6168 up-regulated genes and 249 down-regulated genes) of Mtb-infected Mφ. These DEGs were presented in a volcano plot (Figure 2A and B).

Figure 2 Identification of DEDRGs. (A) Volcano map of dataset GSE5099 shows DEGs of polarized Mφ (left: M1 Mφ, right: M2 Mφ). The red dots represent significantly upregulated genes, while the blue dots represent significantly downregulated genes. (B) The volcanic map of dataset GSE52819 shows DEGs of Mtb-infected Mφ. (C) Venn plot shows DEGs of M1/M2 Mφ infected with Mtb. (D) Venn plot shows intersection genes of Mtb-infected M1/M2 Mφ and apoptosis, autophagy, ferroptosis, or pyroptosis related genes, which names differentially expressed death related genes (DEDRGs). (E) Heat map of differential expression of DEDRGs in GSE52819.

The up- and down-regulated genes of Mφ infected with Mtb intersected with the up- and down-regulated genes of M1/M2 Mφ, respectively. Obtaining 748 DEGs of M1 Mφ infected with Mtb(M1 DEGs, including 704 upregulated genes and 44 downregulated genes), as well as 107 DEGs of M2 Mφ infected with Mtb(M2 DEGs, including 106 upregulated genes and 1 downregulated gene) (Figure 2C). Seven hundred and one apoptosis genes, 574 autophagy genes, 502 ferroptosis genes, and 65 pyroptosis genes were obtained from the database collection. As a result, we obtain the differential expression of death related genes (DEDRGs) of M1/M2 Mφ infected with Mtb through the Venn intersection (Figure 2D, Table 2). The heatmap shows the expression pattern of the DEDRGs in GSE52819 (Figure 2E).

Table 2 Intersection Genes of Differentially Expressed Genes of M1 or M2 Macrophages Infected with Mtb and Death Related Genes (DRIGs)

GO and KEGG Enrichment Analysis of DEDRGs

We performed GO and KEGG enrichment analysis of DEDRGs from M1 Mφ and M2 Mφ, respectively, to determine the potential functions of DEDRGs in Mtb infection. GO enrichment analysis of M1 DEDRGs revealed that these genes were involved in the regulation of apoptotic signaling pathway, I-kappaB kinase/NF-kappaB signaling, and ubiquitin-like protein ligase binding (Figure 3A). GO enrichment analysis of M2 DEDRGs showed that these genes were involved in the regulation of apoptotic signaling pathway, growth factor receptor binding, and associated with autophagosome (Figure 3B). KEGG enrichment analysis of M1 DEDRGs showed that these genes were associated with signaling pathways such as Apoptosis and NOD like receiver signaling pathway (Figure 3C). KEGG enrichment analysis of M2 DEDRGs showed that these genes were associated with signaling pathways such as Apoptosis, Autophagy-animal, and NF-kappa B signaling pathway (Figure 3D).

Figure 3 Enrichment analysis of DEDRGs. (A) GO enrichment analysis of DEDRGs of M1 Mφ. (B) GO enrichment analysis of DEDRGs of M2 Mφ. (C) KEGG enrichment analysis of DEDRGs of M1 Mφ. (D) KEGG enrichment analysis of DEDRGs of M2 Mφ.

Identification of Hub Genes

To determine the interactions between DEDRGs, the intersecting genes of Mtb-infected M1 Mφ with apoptosis, autophagy, pyroptosis, and ferroptosis as well as the intersecting genes of Mtb-infected M2 Mφ with apoptosis were subjected to PPI analyses. These PPI networks were also analyzed using the CytoHubba plugin to obtain the Hub genes (M1-Hub, M2-Hub) associated with cell death and M1/M2 Mφ infected with Mtb (Figure 4A and B). Because the number of autophagy and ferroptosis related genes in M2 DEDRGs was too small, all of them were listed as Hub genes (Table 3). Box line plots were drawn to show the differential expression of these Hub genes between the Mtb-infected and uninfected groups in GSE52819 (Figure 4C and D).

Table 3 Hub Genes Were Screened by PPI Analysis

Figure 4 Screening for Hub genes. (A) The PPI network of DEDRGs in M1 Mφ. Use red to represent the Hub gene. (Upper left: apoptosis, upper right: autophagy, lower left: ferroptosis, lower right: pyroptosis). (B) The PPI network of DEDRGs of M2 Mφ and apoptosis related genes. (C) Box plot of M1 Hub gene expression in GSE52819. (D) Box plot of M2 Hub gene expression in GSE52819.

ROC Curve Analysis and Correlation Analysis

Peripheral blood expression data of 9 active TB patients and 6 healthy people from GSE54992 were analyzed. These data were utilized to perform ROC curve analysis of the above Hub genes to further validate the reliability of these Hub genes in the peripheral blood of patients (Figure 5A and B). The Hub genes with AUC >80% may have excellent specificity and sensitivity for the next experimental validation, including 8 genes in M1-Hub and 6 genes in M2-Hub (Table 4). Mapping the correlation network of Hub genes showed that in M1-Hub, NFKB1, IL6, and SOCS1 as the center of the network graph were closely correlated with each gene, and in M2-Hub, TNF, BIRC3, DEPTOR, and SOCS1 as the clusters of genes had a strong correlation (Figure 5C and D).

Table 4 Genes to Be Verified

Figure 5 ROC curve analysis and correlation analysis. (A) ROC curve screening for M1-Hub with AUC >80%. (B) ROC curve screening for M2-Hub with AUC >80%. (C) Correlation network diagram of M1-Hub. (D) Correlation network diagram of M2-Hub.

Immune Infiltration Analysis

Boxplots showed that the proportions of Macrophages M0 and Macrophages M1 were significantly lower in active TB samples compared with healthy controls (Figure 6A). Heatmaps showed the correlation of the three types of Mφ with M1-Hub and M2-Hub genes. It was evident that M1 Mφ had a strong correlation with M1-Hub (Figure 6B and C).

Figure 6 Immunoinfiltration analysis. (A) Box plot of differential infiltration of 22 immune cells in tuberculosis patients and healthy people. Using the wilcox.test function for Wilcoxon rank sum test, calculate the p-value between two sets of samples, *P<0.05 **P<0.01 ***P<0.001. (B) The correlation heatmap between M1-Hub and three types of Mφ. (C) The correlation heatmap between M2-Hub and three types of Mφ. Abbreviation: ns, nonsignificant.

Mφ Polarization

The Mφ polarization process was shown in Figure 7A. qRT-PCR was performed to detect the expression of CD80, IL-1β, TNFα, CD206, IL-10, CCL22 in the three groups (Figure 7B). Compared with M0 Mφ, the M1 Mφ group significantly highly expressed the M1 Mφ marker molecules (CD80, IL-1β, TNFα), while the M2 Mφ group significantly highly expressed the M2 Mφ marker molecules (CD206, IL-10, CCL22). Western Blot verified the expression of M0/M1/M2 Mφ of IL-1β, IL-10 (Figure 7C) and determined that the Mφ was successfully oriented to polarize.

Figure 7 Macrophage polarization. (A) The flowchart of inducing THP-1 polarization into M1/M2 Mφ. Created with BioRender.com. (B) PCR validation of M1/M2 Mφ marker genes. Compare data between multiple groups using one way ANOVA. (C) WB validation of M1/M2 Mφ marker genes. Compare data between multiple groups using one way ANOVA. The statistical chart shows the analysis results of three independent repeated experiments. **P<0.01 ***P<0.001 ****P<0.0001.

Experimental Proof of Mode of Death

Using M0 Mφ as a control, qRT-PCR and Western Blot were used to detect the expression of marker genes of the cell death modes of M1/M2 Mφ infected with Mtb (Figure 8A and B). It can be seen that the gene expression of CASP3, LC3B, FTH1, and CASP1 in M1 Mφ infected with Mtb was significantly higher than that of M0 Mφ (P<0.05), while there was no significant difference in M2 Mφ (P>0.05). It indicates that the occurrence of apoptosis, autophagy, ferroptosis, and pyroptosis in Mtb-infected M1 Mφ was significantly higher than that in M2 Mφ.

Figure 8 Death modes and Hub genes validation. (A) PCR detection of marker genes of four death modes in M1/M2 Mφ infected with Mtb. (B) WB validation of marker genes of four death modes in M1/M2 Mφ infected with Mtb. The statistical chart shows the analysis results of three independent repeated experiments. Compare data between multiple groups using one way ANOVA. (C) PCR validation of M1 Mφ Hub genes. (D) PCR validation of M2 Mφ Hub genes. The P value is calculated using paid T text. All experiments were conducted with three independent replicates. *P<0.05 **P<0.01 ***P<0.001 ****P<0.0001.

Validation Hub Genes by qRT-PCR Analysis

M1-Hub and M2-Hub were detected by qRT-PCR, respectively (Figure 8C and D). Compared with M0 Mφ, NFKB1, TNF, CFLAR, TBK1, IL6, RELA, AIM2 were significantly up-regulated in Mtb-infected M1 Mφ, as well as TNF, BIRC3 were significantly up-regulated and DEPTOR was significantly down-regulated in Mtb-infected M2 Mφ, which was in accordance with the results of bioinformatics analysis.

Discussion

As an infectious disease, tuberculosis seriously jeopardizes the life and health security of all human beings. Mφ, as the frontline of innate immune defense, can kill Mtb by controlling self-death.¹⁹ However, a large number of Mtb aggregated and internalized within the Mφ phagosome will induce M1 Mφ involuntary deaths and even lead to cell death cascade.^20,21 As a result, the Mφ in the organism is biased from M1 Mφ to M2 Mφ,⁶ promoting the immune escape of Mtb.

Our study obtained the intersecting genes (DEDRGs) of Mtb-infected M1/M2 Mφ with death-related genes from bioinformatics analysis. It was found that there was a difference in the DEDRGs between M1 and M2 Mφ, which indicated that there may be differences in their death modes. To test the role of these DEDRGs in the intrinsic immunity of the organism, we conducted KEGG and GO enrichment analysis and found that they were significantly correlated with the cell death modes. The activation of NF-kappaB signaling promotes BCG infection-induced pyroptosis of THP-1,²² and induces ferroptosis through the widespread effects of multiple factors.²³ In addition, NF-kappaB is involved in numerous protein kinase cascades that regulate the crosstalk between autophagy and apoptosis.²⁴ The NOD-like receptor signaling pathway and ubiquitin-like protein ligase binding have also been widely shown to be associated with autophagy.^25,26 To further determine the diagnostic value of these DEDRGs, we constructed PPI networks to access highly connected sub-networks in bioregulation and analyzed the AUC values by ROC curves to identify genes with high specificity and sensitivity, thus identifying the Hub genes. In addition, to explore the level of macrophage infiltration in overall immune cells after Mtb infection and the association of Hub genes with two major macrophage subtypes. We performed immune infiltration analyses that showed differences in infiltration of M1 macrophages between TB patients and healthy individuals, as well as a strong correlation between M1 macrophages and M1-hub.

Mφ subpopulations exercise differential control over Mtb infection due to different gene expression patterns, with M1 and M2 Mφ expressing pro-and anti-inflammatory cytokine and chemokines, respectively.²⁷ Our experiments determined that compared to M0 Mφ, the occurrence of apoptosis, autophagy, ferroptosis, and pyroptosis was significantly higher in Mtb-infected M1 Mφ than in M2 Mφ. Thus, inhibition of death in Mtb-infected Mφ may become a key component in controlling Mtb transmission. Exploring the biomarkers associated with Mφ polarization and different death modes under Mtb infection will greatly facilitate the development of this critical link.

Mtb infection promotes the production of the pro-inflammatory cytokine IL-1β, prompts activation of AIM2 inflammasomes to induce pyroptosis or autophagy,^28,29 and kills intracellular Mtb in a TBK1-dependent manner.³⁰ In addition, the inflammatory factor TNF activates the NF-kappaB signaling pathway centered on NFKB1 to induce apoptosis and inflammatory responses in Mφ.^31,32 In addition, NFKB1 activates autophagy by upregulating the transcription of Beclin 1 (a key autophagy gene).³³ Whereas the onset of ferroptosis is induced by activation of IL-6/ferritin signaling,³⁴ and negatively regulated by RELA.³⁵ Among the M2-hub genes, BIRC3 is a closely related member of the inhibitor of apoptosis proteins (IAP family) and plays an important role in apoptosis regulation;³⁶ and DEPTOR acts as an inhibitor of cell proliferation by enhancing autophagy.³⁷ M1 macrophages highly express NFKB1, TNF, AIM2, IL-6 and other Hub genes after Mtb infection, which promote or participate their own apoptosis, autophagy, pyroptosis and ferroptosis. On the contrary, although M2 macrophages highly express Hub genes such as BIRC3 and DEPTOR to lead to death, they have a much lower proportion of the deaths compared with M1 macrophages. As a result, the majority of macrophages in the body gradually changed from M1 type to M2 type. It was experimentally verified that NFKB1, TNF, CFLAR, TBK1, IL6, RELA, and AIM2 were dysregulated in M1 Mφ infected with Mtb. TNF, BIRC3, and DEPTOR were dysregulated in M2 Mφ infected with Mtb. These results showed that there are multiple different death modes in Mφ with different polarization states.

Conclusion

In this work, we validated the potential diagnostic efficacy of NFKB1, TNF, CFLAR, TBK1, IL6, RELA, AIM2, BIRC3, and DEPTOR in TB. We also investigated the correlation between Mφ subtypes infected with Mtb and different death modes and demonstrated that M1 Mφ is more susceptible to apoptosis, autophagy, ferroptosis, and pyroptosis upon Mtb infection, which partially explains the mechanism of the transition of Mtb-infected Mφ from M1-dominant to M2-dominant (Figure 9). This work is not only of great significance for the exploration of potential biomarkers of TB but also provide more deep insights into the innate immune responses against Mtb infection and the immune escape mechanisms of Mtb.

Figure 9 Visualization of the biological process of M1 macrophages undergoing apoptosis, autophagy, pyroptosis, ferroptosis, and gradually transforming into M2 macrophages. Created with BioRender.com.