Tumor-associated macrophages based signaling pathway analysis and hub genes identification in glioma.

ABSTRACT: Tumor-associated macrophages (TAMs) play a crucial role in the immune response to many malignancies, but the signaling pathways by which the glioma microenvironment cross-talk with TAMs are poorly understood. The aim of this study was to uncover the potential signaling pathways of the regulation of TAMs and identify candidate targets for therapeutic intervention of glioma through bioinformatics analysis.Chinese Glioma Genome Atlas (CGGA) and The Cancer Genome Atlas (TCGA) datasets were used to download RNA-Seq data and microarray data of human glioma specimen. Differentially expressed genes (DEGs) between CD68-high samples and CD68-low samples were sorted. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis of the DEGs was conducted. Protein-protein interaction (PPI) network were formed to identify the hub genes.The prognostic value of TAMs in glioma patients was confirmed. A total of 477 specific DEGs were sorted. The signaling pathway was identified in pathway enrichment and the DEGs showed prominent representations of immune response networks in glioma. The hub genes including C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8 were identified form the PPI network and they were all correlated positively with the expression of CD68 and showed the excellent prognostic value in glioma patients.TAMs can be used as a good prognostic indicator in glioma patients. By analyzing comprehensive bioinformatics data, we uncovered the underlying signaling pathway of the DEGs between glioma patients with high and low expression level of CD68. Furthermore, the 6 hub genes identified were closely associated with TAMs in glioma microenvironment and need further investigation.

Introduction

Globally, there are about 100,000 people are diagnosed as having diffuse gliomas every year.[ Diffuse glioma contains glioblastoma (GBM) and low grade glioma (LGG) and they are related to substantial morbidity and mortality although it comprises <1% of all newly diagnosed cancers.[ GBM, the most common and lethal primary brain tumor in adults, accounts for 70% to 75% of all diffuse glioma diagnoses and has a median overall survival of around 15 months despite aggressive treatment.[ By several mechanisms such as the induction of immunosuppression, fast proliferation, the promotion of angiogenesis and its propensity to infiltrate vital brain structures, glioma can achieve rapid growth and dissemination within the brain.[ Among these mechanisms, escape from immune surveillance is gradually recognized as a landmark event in glioma biology.[ Similar to many other solid malignancies, the tumor microenvironment (TME) was supposed to play a critical role in the control of the immune response to glioma.[

Inside the tumor, tumor-associated macrophages (TAMs) are a key component of the local TME, and they can contribute greatly to tumor immune system evasion, suppress T-cell activity, and control cancer initiation, progression, metastasis in an array of malignancies.[ In GBM, TAMs constitute more than 30% of infiltrating cells and TAMs infiltration is closely related to glioma progression and the outcome of glioma patients.[ Even though TAMs-induced immunosuppression was found in the different phases of the anti-glioma immune response, the signaling pathways by which glioma microenvironment controls TAMs and glioma-specific immunity are not fully understood. In some malignancies which including glioma, the acquisition of a polarized phenotype resembling anti-inflammatory (M2) macrophage by TAMs has been associated with the suppression of tumor-specific immunity and tumor progression in tumor pathogenesis.[ While the acquisition of a polarized phenotype resembling pro-inflammatory (M1) macrophage by TAMs has been associated with the promotion of tumor-specific immunity and tumor regression in tumor pathogenesis.[ However, the M1/M2 transition of TAMs is a dynamic process and this model has its own drawbacks such as the intermediate phenotypes of TAMs between M1 and M2 do not fit this model although the intermediate states do exist.[ In view of the complicated biological function of TAMs in TME which may contribute to their pathogenic activities, the identification of the molecular mechanisms that control TAMs in glioma has crucial basic and clinical significance.

Here, using CD68, which is widely used as a cellular marker for macrophages,[ we confirmed the prognostic value of TAMs in glioma patients at the level of gene. Through Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis of the differentially expressed genes (DEGs) between glioma patients with the high and low expression levels of CD68, we found that these DEGs showed prominent representations of immune response networks in glioma. Then, the protein-protein interaction (PPI) network analysis elaborated the immune system process in which the DEGs are involved in. Finally, the hub genes including C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8 were identified from the PPI network. These hub genes were correlated positively with the expression of CD68 and showed their excellent prognostic value in glioma patients. These findings uncovered the potential signaling pathways of the regulation of TAMs and identified candidate targets for therapeutic intervention of glioma.

Materials and methods

Collection and processing of data sets

The healthy tissue samples were obtained from Genotype-Tissue Expression (GTEx) datasets. Two RNA-Seq data sets and clinical information containing 693 (batch_1) and 325 (batch_2) glioma samples were obtained from Chinese Glioma Genome Atlas (CGGA, http://cgga.org.cn/download.jsp), which were assigned as training sets (Train set 1 and Train set 2). After normalizing these data, CD68 high (top third) and low expression (last third) samples were selected for each data set according to the expression of CD68. Furthermore, a microarray data set including 301 glioma samples was downloaded from CGGA for a test set (Test set 1), and RNA-Seq data set including 702 gliomas was downloaded from The Cancer Genome Atlas (TCGA, http://cancergenome.nih.gov) for another test set (Test set 2). After preprocessed, each test set was divided into two groups of CD68-high and -low by the median expression of CD68. All these processes were conducted by R-studio software (https://rstudio.com/). The information of these data sets was shown in Table 1.

Table 1

Main information of the datasets used in the present study.

	Train set (1)	Train set (2)	Test set (1)	Test set (2)
Data type	RNA-Seq (CGGA)	RNA-Seq (CGGA)	Array (CGGA)	RNA-Seq (TCGA)
Cases (n)	693	325	301	702
CD68 high	231 (top)	109 (top)	150	351
CD68 low	231 (last)	109 (last)	151	351

Screening of differential expression genes (DEGs)

In order to clarify the differences of expression profiles between gliomas with high and low CD68 expression, Limma package was used to screen the DEGs by comparing CD68-high samples with CD68-low samples in Train set 1 or Train set 2, respectively. All genes with a P value < .01 and |log2fold change| > 2 were selected as DEGs for further analysis.

Function enrichment analysis of DEGs

Then, for understanding the biological processes and pathways in which DEGs involved, GO (Gene Ontology) and KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment analysis was performed on the DAVID database (https://david.ncifcrf.gov/home.jsp). The cut-off value was set as gene count >10 and P value < .05.

Construction and analysis of protein-protein interaction (PPI) network

STRING database (http://www.string-db.org/) was used to construct the PPI network of DEGs with high confidence (confidence score > 0.70). Significant clusters were obtained by MCODE plug-in with a score >10. Furthermore, the immune process enrichment of the PPI network was performed by Clue-Go plug-in with P < .01 and genes > 10%/ term. Additionally, hub genes were calculated by CytoHubba plug-in using two algorithms including Degree and Betweenness. All these data were visualized by Cytoscape software (http://www.cytoscape.org/) and R-studio software (https://rstudio.com/).

Statistical analysis

Two-tailed Student's t test was used for analyzing the difference between the two groups. Pearson Chi-square test was used to evaluate the correlation between genes. Overall survival (OS) was analyzed using the Kaplan-Meier method and Log-rank test. P < .05 was considered significantly different.

Results

The prognostic value of TAMs in glioma

Firstly, to investigate the prognostic significance of the infiltrating level of TAMs in glioma, we queried a TAMs associated gene, CD68, and compared its expression between glioma patients and healthy people. We found that there was a significant increase in CD68 expression in glioma including GBM and LGG, more so than normal tissues (P < .01) (Fig. 1A). When we compared the survival of glioma patients with the low or high expression level of CD68 using the TCGA dataset, we found significantly worse outcomes for patients that had higher expression of the CD68 gene (Fig. 1B). Similarly, we founded that the grades of glioma were increased with the increment of CD68 expression and the higher the expression of CD68 and worse outcome the glioma patients end up with using another three datasets (Fig. 1C–E). These results indicate that the role of TAMs in glioma is harmful and the presence of TAMs is unfavorable for survival.

Figure 1

The prognostic value of TAMs in glioma. (A) Expression analysis of CD68 gene between glioma patients (including GBM patients (n = 163) and LGG patients (n = 518)) and healthy people (n = 207) using TCGA and GTEx datasets. (B) Kaplan–Meier graphs of patients with low (n = 338, blue) and high (n = 338, red) CD68 expression using TCGA dataset. (C-E) Up: Expression analysis of CD68 gene among glioma patients of different grades using CGGA datasets including RNA-seq (batch 1) (c), RNA-seq (batch 2) (d) and microarray date (e). Down: Kaplan–Meier graphs of patients with low (blue) and high (red) CD68 expression using CGGA dataset of RNA-seq (batch 1) (c), RNA-seq (batch 2) (d) and microarray date (e). Log-rank Mantel–Cox test, P < 0.0001.

DEGs identification of gliomas with different expression of CD68

To study the difference of the transcriptional program between glioma with high expression of CD68 and low expression of CD68, we first analyzed their DEGs using two different Train sets. In Train set 1, there are 693 glioma patients, including 231 CD68-high cases and 231 CD68-low cases. In Train set 2, there are 325 glioma patients, including 109 CD68-high cases and 109 CD68-low cases (Table 1). In Train set 1, we screened 976 up-regulated genes (Fig. 2A) and 539 down-regulated genes, while in Train set 2 (Fig. 2B), we enriched 93 up-regulated genes and 323 down-regulated genes. Then, by Venn analysis, we identified 435 up-regulated DEGs and 42 down-regulated DEGs that Train set 1 and Train set 2 shared (Fig. 2C-D). Thus, we finally selected a total of 477 DEGs for further analysis.

Figure 2

Differential genes expression analysis. A-B Volcano plot of gene expression with P < .05 (fold change in the relative expression of genes as determined by log2 in Train set 1 (a) and Train set 2 (b), with up-regulated genes are shown in red and down-regulated genes in blue. C-D Venn diagram of differentially expressed up-regulated genes (c) and down-regulated genes (d) between Train set 1 and Train set 2, with 435 up-regulated DEGs are shown in brown and 42 down-regulated genes in dark blue.

GO and KEGG pathway enrichment analysis of DEGs

To further assess CD68-associated signaling pathway in glioma, we conducted GO and KEGG analysis of the DEGs. Specifically, in the biological process aspect of GO analysis, these DEGs showed prominent representations of immune response networks, signal transduction, innate immune response signaling, inflammatory response signaling and proteolysis (Fig. 3A). Of note, among all the selected signaling pathway, most of them are associated with immune response, suggesting the important role of TAMs in regulating glioma immune microenvironment (Fig. 3A). In the cellular component aspect, these DEGs showed prominent representations of plasma membrane, integral component of membrane, extracellular exosome, extracellular space and extracellular region (Fig. 3B). In the molecular function aspect, these DEGs showed prominent representations of serine-type endopeptidase activity, receptor binding, receptor activity, carbohydrate binding and cytokine activity (Fig. 3C). Specifically, in KEGG analysis, these DEGs showed prominent representations of staphylococcus aureus infection, tuberculosis, cytokine-cytokine receptor interaction, phagosome and osteoclast differentiation (Fig. 3D). Moreover, the enriched pathways showed a close relationship with each other (Fig. 3E). Together, these results suggested that the DEGs are mainly involved in the immune response in glioma.

Figure 3

GO and KEGG enrichment analysis of DEGs. (A-C) The biological process aspect (A), the cellular component aspect (b) and the molecular function aspect (c) of GO analysis. (D) KEGG enrichment analysis of DEGs. (E) The relationship between different enriched pathways from KEGG enrichment analysis. The circle represents the gene and the rectangle represents the pathway.

PPI network analysis of DEGs

To further determine the role of these DEGs in the molecular processes of glioma, the PPI network with 233 nodes and 1235 edges was constructed (Fig. 4A). Then, we recognize the most significant module of the network, cluster 1 (score = 17.46) and cluster 2 (score = 10.12). In cluster 1, we identified 49 genes (Fig. 4B) and in cluster 2, we identified 18 genes (Fig. 4C).

Figure 4

PPI network analysis of DEGs. (A) PPI network analysis of DEGs, with 233 nodes and 1235 edges. (B-C) The two most significant modules, including cluster 1 (score = 17.46) and cluster 2 (score = 10.12).

Immune system process enrichment analysis of genes in PPI network

In view of the above results of GO analysis of the DEGs which suggested that these genes are mainly involved in the immune response in glioma. We further analyzed the immune system process of the 233 genes form PPI network to elaborate how these DEGs were took part in the regulation of the glioma immune response. These enriched genes showed prominent representations of lymphocyte proliferation (including T cell proliferation, regulation of lymphocyte proliferation, regulation of T cell proliferation, positive regulation of T cell proliferation and positive regulation of lymphocyte proliferation), peptide antigen assembly with MHC class II protein (including peptide antigen assembly with MHC protein complex), macrophage activation (including microglial cell activation and regulation of macrophage activation), granulocyte activation (including neutrophil degranulation, neutrophil activation involved in immune response and neutrophil activation), neutrophil chemotaxis (including granulocyte chemotaxis, granulocyte migration and neutrophil migration), interferon-gamma-mediated signaling pathway (including response to interferon-gamma, cellular response to interferon-gamma and granulocyte migration), type II a hypersensitivity (including regulation of type II a hypersensitivity, type II hypersensitivity and regulation of type II hypersensitivity), antigen processing and presentation of exogenous peptide antigen, positive regulation of leukocyte migration and regulation of B cell proliferation (Fig. 5).

Figure 5

Significant immune immune system process was identified of all these 233 PPI genes.

Identification and validation of hub genes

Given the important role for TAMs in glioma pathogenesis, we determined to further detect the modulation mechanism of TAM activity in glioma. Through Venn analysis, we identified six Hub genes, namely C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8 (Fig. 6A). An analysis of CGGA data showed that there is higher expression level of these six hub genes in glioma patients with higher expression of CD68 (Fig. 6B). Moreover, a similar result was obtained through the analysis of TCGA data (Fig. 6C). These results suggested the potential role of these hub genes in TAMs associated signaling pathway in glioma.

Figure 6

Identification and characterization of hub genes. (A) Venn diagram of hub genes obtained by Degree and Betweenness algorithms, the hub genes are C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8. (B) Expression analysis of hub genes in Test 1(CD68-low, n = 151; CD68-high, n = 150; P < .05). (C) Expression analysis of hub genes in Test 2 (CD68-low, n = 351; CD68-high, n = 351; P < .05).

The prognostic value of hub genes in glioma patients

To further investigate whether there is positive correlation or negative correlation between these hub genes and cd68, we performed Pearson correlation analysis of the expression level of CD68 and six hub genes in glioma patients using two datasets. In Test set 1, there are 301 glioma patients, including 150 CD68-high cases and 151 CD68-low cases. In Test set 2, there are 702 glioma patients, including 351 CD68-high cases and 351 CD68-low cases (Table 1). As we can see, the expression level of all these six hub genes are correlated positively with the expression of CD68 with different correlation coefficient, namely C3 (R = 0.76), IL6 (R = 0.40), ITGB2 (R = 0.86), PTAFR (R = 0.80), TIMP1 (R = 0.60) and VAMP8 (R = 0.88) (Fig. 7A). Moreover, when we compared the survival of glioma patients with the low or high expression level of these hub genes, we found significantly worse outcomes for patients that had higher expression of these hub genes (Fig. 7B). Similarly, we validated these results using the TCGA dataset. The expression level of all these six hub genes are correlated positively with the expression of CD68 with different correlation coefficient, with C3 (R = 0.57), IL6 (R = 0.25), ITGB2 (R = 0.89), PTAFR (R = 0.83), TIMP1 (R = 0.44) and VAMP8 (R = 0.73) (Fig. 7C). And again, these hub genes showed excellent prognostic values in glioma patients (Fig. 7D), indicating their potential roles of therapeutic candidates in glioma.

Figure 7

The prognostic value of hub genes in glioma patients. (A) Pearson correlation analysis of the expression level of CD68 and 6 hub genes in glioma patients using the data of Test set 1 (n = 301), with C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8. (B) Kaplan–Meier graphs of patients with low (blue) and high (red) hub genes expression using the data of Test set 1. (C) Pearson correlation analysis of the expression level of CD68 and six hub genes in glioma patients using the data of Test set 2 (N = 702), with C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8. (D) Kaplan–Meier graphs of patients with low (blue) and high (red) hub genes expression using the data of Test set 2.

Discussion

Recent studies demonstrated that the interplay between TME components and cancer cells play critical roles in the maintenance of cancer characteristics, including treatment resistance, immune sequestration and angiogenesis.[ The TME can be molded by secreted factors such as different growth factors and cytokines and cancer cell-intrinsic signaling pathways.[ Besides tumor cells, the infiltrating immune cells and stromal cells within TME are dominant cell components that can influence the malignant progression of tumors.[ In the glioma microenvironment, the most common non-neoplastic cells are TAMs, which include macrophages of peripheral origin and brain-intrinsic microglia, that create a supportive stroma for glioma cell proliferation and invasion.[ Through specific TAMs-glioma iterative interactions, a particular glioma ecosystem is builded up, which offers new strategies for glioma therapeutic targeting.[ Therefore, we set out to determine the role of TAMs in the glioma microenvironment. First, we confirmed the prognostic value of TAMs in glioma patients at the level of gene and transcription using TCGA, CGGA and GTEx datasets. The higher the expression of CD68 and worse outcome the glioma patients end up with, which was consistent with previous studies, suggesting the pivotal role of TAMs in glioma microenvironment.[ Then, to further determine the underlying molecular difference between glioma with high expression of CD68 and low expression of CD68, we performed DEGs analysis and detected 976 up-regulated genes and 539 down-regulated genes in the data of Train set 1, and 93 up-regulated genes and 323 down-regulated genes in the data of Train set 2. Through Venn analysis, we identified 435 up-regulated DEGs and 42 down-regulated DEGs that train set 1 and train set 2 share together, and we finally selected total 477 DEGs for further investigation.

In support of an important role for TAMs in tumor pathogenesis, recent studies which were about the modulation of TAM activity in glioma have recognized some specific molecules and signaling pathways which has beneficial effects on GBM.[ Thus, we performed GO and KEGG analysis of the DEGs between glioma patients with high and low expression level of CD68 to uncover the mechanisms behind these DEGs. To identify the terms which contain more genes in different categories, we define the interested terms as the terms which constitute more genes. In the biological process aspect of GO analysis, we identified 22 interested terms. Among all the selected signaling pathways, most of them are associated with immune response, such as innate immune response, inflammatory response, complement activation, interferon-gamma-mediated signaling pathway and T cell proliferation, suggesting the crucial role of TAMs in regulating glioma immune microenvironment. In the cellular component aspect of GO analysis, we identified 12 interested terms, most of them are associated with membrane and extracellular functions, such as plasma membrane, integral component of membrane, extracellular exosome and extracellular space, indicating the function of these genes may be implemented in membrane structures and extracellular space. In the molecular function aspect of GO analysis, we identified 6 interested terms, namely serine-type endopeptidase activity, receptor binding, receptor activity, carbohydrate binding, cytokine activity and antigen binding. In the pathway category of KEGG analysis, we identified 18 interested terms, most of them are associated with pathogen infections, such as tuberculosis, staphylococcus aureus infection, cytokine-cytokine receptor interaction and phagosome, indicating the potential role of these genes in the body's resistance to microorganisms.

Then, to assess the immune system process of some key molecules, we performed PPI network analysis. Most of the enriched genes from PPI network were involved in the process of adaptive immune response, such as peptide antigen assembly with MHC class II protein, lymphocyte proliferation, antigen processing and presentation of exogenous peptide antigen and regulation of B cell proliferation. Through Venn analysis, we identified six hub genes, namely C3, IL6, ITGB2, PTAFR, TIMP1 and VAMP8. Previous studies have shown that along with complement C3, T cells can mediate photodynamic therapy-induced anti-glioma responses.[ IL6 can contribute to glioma progression through different signaling pathway such as NFAT1-regulated IL6 pathway, p-STAT3-MIR155–3p-CREBRF pathway, IL6/JAK/STAT3 Pathway.[ What's more, there were also other studies which shown that ITGB2, PTAFR, TIMP1 and VAMP8 are involved in glioma evolution.[ Even though, these studies were only the initial exploration of the function of these genes in glioma and the exact mechanisms behind are still remain in a state of incomplete understood. Combined with our results that the expression level of all these six hub genes are correlated positively with the expression of CD68 and higher expression of these hub genes signified worse outcomes for glioma patients, we may conclude that these hub genes may be potential candidate targets for therapeutic intervention of glioma.

Conclusion

In conclusion, we confirmed the prognostic value of TAMs in glioma patients. Through bioinformatics analysis, we uncovered the underlying signaling pathway of the DEGs between glioma patients with high and low expression level of CD68. Furthermore, the six hub genes identified were closely associated with TAMs in glioma microenvironment and need further investigation.

Author contributions

Conceptualization: Xiang Wang, Weihai Ning.

Data curation: Xiang Wang, Weihai Ning.

Formal analysis: Xiang Wang, Weihai Ning.

Funding acquisition: Hongwei Zhang, Chunjiang Yu.

Methodology: Chunjiang Yu.

Writing – original draft: Xiang Wang, Weihai Ning.

Writing – review & editing: Zhiqiang Qiu, Shenglun Li, Hongwei Zhang.