Systematic Analysis of Chemokines Reveals CCL18 is a Prognostic Biomarker in Glioblastoma

Introduction

Glioblastoma (GBM) is the most common primary malignant tumor in the brain, with an incidence rate of approximately 3.47 per 100,000 people.¹ GBM, the World Health Organization (WHO) grade IV glioma, has the worst prognosis.² Though there are varieties of approved treatments for GBM, including surgical resection, radiotherapy, and chemotherapy, the clinical outcome remains dismal.³ For primary GBM patients, the median overall survival (OS) time was eight months with 95% confidence intervals, and the five‐year survival rate was only about 6.8%.¹ The possible reasons are the highly invasive nature of GBM cells, the chemo‐ and radio‐resistance, the high level of vascularization, complex cell composition, and decreased effusion of chemotherapeutic drugs due to the blood-brain barrier (BBB).^4,5 Many prognosis biomarkers for GBMs have been developed, such as mutations in isocitrate dehydrogenase (IDH) genes, MGMT gene promoter methylation status, chromosome 1p/19q deletion, and tumor protein P53 mutation.^1,2 In the WHO 2016 classification, secondary glioblastomas with an IDH mutation originate from lower grade precursor lesions.² Significantly, mutations in IDH1 are correlated with a better prognosis. However, these mutations are rarer in primary GBM than in most common adult gliomas.⁶ Moreover, antiangiogenic treatments like anti-VEGFR/VEGF drugs do not have significant improvement in survival. Monoclonal antibody bevacizumab might improve life quality, but not OS or progression-free survival (PFS) of GBM patients.⁷ Consequently, further efforts are urgently needed to understand the biological behavior of GBM and develop novel methods for diagnosis and treatments.

The tumor microenvironment (TME), one of the critical factors of GBM development and treatment, is the primary interaction location between tumor cells and the host immune system.⁸ Through interactions between chemokines and chemokine receptors, different immune cell subsets are recruited into the TME.⁹ The chemokines are the largest subfamily of cytokines. They can be further subdivided into four main classes depending on the first two cysteines (C) residues in their protein sequences: namely, the CC-chemokines, the CXC-chemokines, the C-chemokines, and the CX3C-chemokines.^8,10 Chemokines are also described functionally as inflammatory, homeostatic, or dual-function chemokines.¹¹ Their ligand-receptor relationships are promiscuous, with a single ligand-binding different receptor and vice versa.⁸ Recently, chemokine signaling pathways associated with gliomas have been reviewed.^12,13 Specific chemokines have distinct effects on tumor growth, metastasis, the low to high-grade gliomas transition, and therapeutic outcomes.^8,14,15 Especially, high-risk GBM patients exhibited higher expression level of CCL18, which correlated with high levels of Tim-3, indicating an association between immune cell infiltration and GBM prognosis.^16,17 However, a systematically contrastive analysis of the role of chemokines in human glioma prognosis is limited. Considering that the chemokine network is highly complicated, it seems unlikely that any single chemokine could be a tumor marker sufficiently effective for GBM diagnosis. Therefore, using a novel combination of multiple biomarkers was proposed as an alternative way to diagnose of patients with GBM.

The aim of this study is to investigate the prognostic value of chemokines in GBM precisely and provide a model for clinical diagnosis. Our analytic workflow follows the standard bioinformatics studies¹⁸ and is briefly outlined in Figure S1. Using the cancer genome atlas (TCGA) data, we explored the comprehensive roles of chemokines in GBM. For the first time, the study comprehensively provided the expression of the whole family of chemokine members, and compared their impacts on GBM. Our results indicated that seven chemokines, including CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13, could serve as prognostic factors for predicting the prognosis of GBM patients. Furthermore, the ELISA assay and the Kaplan-Meier plotter analysis for protein expression level in a cohort of glioma patients confirmed that CCL18 could be a potential independent biomarker in GBM clinical applications, including glioma diagnosis and drug development.

Materials and Methods

Data Processing and Expression Analysis

The mRNA expression data of 42 chemokines (CCL1-5, CCL7-8, CCL11, CCL13-28, CXCL1-3, CXCL5-6, CXCL8-14, CXCL16-17, CX3CL1, XCL1-2) were downloaded from the TCGA database (http://cancergenome.nih.gov/abouttcga), and the respective normal tissue samples were downloaded from the Genotype-Tissue Expression (GTEx, http://gtexportal.org/home/datasets) database. TCGA provided 168 GBM samples containing the prognostic information, while GTEx provided 1157 normal brain tissue samples. For RNA-seq data, expression levels were transcripts per million (TPM) normalized. Expression data for all chemokines were Log2 transformed, and the Wilcoxon rank-sum test was conducted on these tumor types. P < 0.05 was considered to indicate differential expression between tumor and normal tissues. Data analysis was conducted using R software (Version 3.6.3), and the R package “ggplot2” was used to draw box plots.

Survival Analysis

The relationship between each chemokine expression and patients’ prognosis (OS: overall survival) interval in 23 cancers was visualized with forest plots. The hazard ratio (HR) and 95% confidence intervals were calculated via univariate survival analysis. The Kaplan-Meier survival analysis for OS, DSS (disease-specific survival), and PFI (progression-free interval) of seven chemokines (CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13) in GBM via R packages “survminer” and “survival” was conducted to compare the survival difference. We selected the clinical characters of the age, gender, IDH status, and risk score to construct a prognostic nomogram to help predict the probability of 1- and 2-, OS for GBM patients via R packages “rms” and “survival”. The discrimination of the nomogram was calculated by the concordance index (C-index), and the time ROC analysis was performed by R packages “timeROC” and “ggplot2” to compare the predictive accuracy of gene and risk score. In addition, the prognostic value of seven chemokines (CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13) in GBM was validated through the Chinese Glioma Genome Atlas (CGGA, http://www.cgga.org.cn).

Immune Response Prediction

The correlation of the prognosis-related chemokines with levels of immune cell infiltration (including activated dendritic cells, aDCs, B cells, CD8 T cells, cytotoxic cells, DC, eosinophils, immature dendritic cells, iDCs, macrophages, mast cells, neutrophils, NK CD56bright cells, NK CD56dim cells, NK cells, plasmacytoid dendritic cells, pDCs, T cells, T cells, T helper cells, Tcm, Tem, Tfh, Tgd, Th1 cells, Th17 cells, Th2 cells, and TReg) was detected in GBM using the R package “GSVA”. Immune cell types were obtained from the article published by Bindea et al.¹⁹ The ssGESA analysis was introduced to quantify the relative infiltration of 24 immune cell types in the GBM microenvironment as described before.²⁰ The correlation analysis between the prognosis-related chemokines with each immune cell infiltration was performed with the Spearman’s test, and the correlation score was normalized to unity distribution, for which zero is the minimal, one is the maximal score for positive correlation and minus one is the maximal score for negative correlation.

GO and KEGG Analysis

We performed Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analyses of the co-expression genes by R packages “clusterProfile”, “org.Hs.eg.db” and “ggplot”. The clusterProfiler package offers a gene classification method, namely groupGO, to classify genes based on their projection at a specific level of the GO terms, and provides functions, enrichGO and enrichKEGG, to calculate enrichment test for GO terms and KEGG pathways based on hypergeometric distribution.²¹ In this study, we used p-values < 0.05 and q-values < 0.2 as the cutoff criteria through R package clusterProfiler.²¹

Protein-Protein Interaction (PPI) Network

A PPI network was conducted based on DEGs using the STRING database (http://string-db.org)²² and visualized by R packages “igraph” and “ggraph”. The cut-off value was defined as an interaction score (median confidence) of 0.4.

Clinical Glioma Samples

We used 72 glioma samples (15 WHO grade II gliomas, 14 WHO grade III gliomas, and 43 WHO grade IV GBM) to validate the results found from the database materials. In order to keep the comparative analysis consistent and focus on the validation of prognosis biomarkers in GBM (G4), we grouped grade II and grade III tumors together as G2 & G3, comparing with G4 in both mRNA and protein levels. The study complied with the Declaration of Helsinki and was reviewed and approved by the Ethics Committee of the Wenzhou Medical University, China. Written informed consent was obtained from all patients, their next of kin, or another surrogate decision-maker, as appropriate. Clinical specimens were obtained from glioma patients who underwent surgery at the Department of Neurosurgery, First Affiliated Hospital of Wenzhou Medical University, between July 2015 and March 2020. All involved patients were 18 to 77 years old, had detailed clinical history and follow-up information, and had no prior radiotherapy to the brain and no intracranial abscess within six months before surgery.

ELISA Analysis

The enzyme-linked immunosorbent assay (ELISA) kits for CCL2, CCL18, and CXCL5 were purchased from ABclonal Biotechnology Co., Ltd. Total protein was extracted with a RIPA buffer (50 mM Tris-HCl, pH 7.4, 150 mM NaCl, 1mM EDTA, 1% Triton X-100 and 1% sodium deoxycholate) supplemented with protein inhibitors (Roche, Switzerland). The relative levels of chemokines (CCL2, CCL18, and CXCL5) in tissue lysates were firstly quantified by the ELISA kits according to manufacturer’s instructions, and then normalized to that of β-actin. We compared the expression level for each chemokine at different grades (G2 & G3 vs G4), IDH1 status (wild-type vs mutant), gender (female vs male), and Ki67 values (< 30% vs ≥30%).

Statistical Data

All the data of gene expression was normalized by log2 transformation. Comparison of normal tissue and cancer tissue were used the Wilcoxon rank-sum test. The Kruskal–Wallis test was adopted to analyze the associations between clinical phenotypes and expression levels of 42 chemokines in GBM. The correlation analysis between the two variables was used with the Spearman’s or Pearson’s test. In the survival analysis, the HRs and p values were calculated by the univariate Cox regression analysis or Log rank test. Kaplan-Meier curves were used to compare the survival of patients stratified according to different levels of each chemokine expression. p < 0.05 was set as the significance threshold for all statistical analyses.

Results

Increased Chemokines Expression Levels in GBM

We compared the transcriptional levels of all chemokines in GBM with those in normal samples by using the GTEx and TCGA databases. Patient demographics and clinicopathological characteristics for all patient cohorts used in our study are shown in Table S1. The mRNA expression levels of five members of the CCs subfamily, including CCL15, CCL19, CCL21, CCL24, and CCL28, were significantly downregulated in patients with GBM, while 16 of 24 CCs were overexpressed in GBM versus normal tissues (Figure 1A). As for the other subfamilies of chemokines, the mRNA expression levels in most members of CXCs, XCs, and CX3CL1 were significantly elevated in GBM, except that CXCL17 was reduced significantly (Figure 1B). The pathological stage analyses revealed the expression levels of CCL1-3, CCL5-14, CCL18-27, CXCL1-3, CXCL6-11, CXCL13-16, CX3CL1, and XCL1-2 regularly fluctuated between each of the two different grades (Figure 1C). Compared with grade 2 and grade 3 gliomas (G2, G3), higher expression of CCL2, CCL5-8, CCL13-14, CCL18, and CCL22-26 occurred in GBM (G4) (Figure 1C). Moreover, the results showed that 11 of 14 CXCLs (except CXCL5, CXCL12, and CXCL17), CX3CL1, and XCLs mRNA expressions in patients with GBM were higher than those in patients with G2 and G3 gliomas (Figure 1D). IDH gene status has been considered an effective factor in predicting the overall survival (OS) of glioma patients.⁶ Therefore, we compared each chemokine expression in different IDH statuses of GBM from the TCGA database. Higher expressions of CCL2, CCL26, CCL28, CXCL5, CXCL10-11, CXCL14, and XCL1 were found in IDH wild-type GBM (Figure 1E and F). However, only two chemokines, CCL19 and CCL21, had lower expression levels in IDH wild-type gliomas (Figure 1E). The above results showed that the expression of almost all chemokines upregulated significantly with the increase of tumor malignancy.

Figure 1 Higher chemokine expression correlates with glioma malignancy. (A and B) 41 chemokines expression profiles at mRNA level in the normal and GBM tissues from the GTEx and TCGA databases. (C and D) 41 chemokines expression profiles at the mRNA level by clinical characteristics in the TCGA database were shown based on the grade in gliomas. Grade II: G2; Grade III: G3; Grade IV: G4. (E and F) 41 chemokines expression levels were shown based on IDH mutant status in GBM. The pink rectangles marked the chemokines that had increased mRNA expression levels, while the blue rectangles marked the chemokines that had decreased mRNA expression levels. *p < 0.05, **p < 0.01, ***p < 0.001.

Identification of Seven Prognosis-Related Chemokines in GBM

Following the univariate Cox regression analysis, CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13 were found to have significant prognostic correlations with OS (overall survival) (Figure 2A). Next, the Kaplan-Meier Plotter tool was used to analyze the correlation between the mRNA levels of chemokines and the survival of patients with GBM. The results revealed that the increased mRNA levels of CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13 were significantly associated with lower disease-specific survival (DSS) (p < 0.05) of all of the patients with GBM (Figure 2B). Analysis of progress-free interval (PFI) data revealed associations between high expression levels of CCL2, CCL8, CCL18, and CXCL5 with the poor prognosis in patients with GBM (Figure 2C). Furthermore, the factors including age, gender, IDH status, the expression level of CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, CXCL13, and risk score, were enrolled to construct a prognostic prediction nomogram to predict the 1- and 2-year survival probability of patients with GBM (Figure S2A). The nomogram showed that IDH status in GBM contributed mostly to prognosis, followed by CCL18 expression level. Each level of every variable was assigned a score on the point scale, and a total score was obtained by adding the scores for each of the selected variables (Figure S2A). Next, to evaluate the predictive efficiencies of CCLs and CXCLs in the 1- and 2-year survival rates in GBM, we performed a receiver operating characteristic (ROC) curve analysis using the TCGA dataset. For CCL2, CCL8, CCL18, and CCL28 at 1-year stage, the area under the curve (AUC) was 0.65 (95% CI, 0.561 to 0.738), 0.594 (95% CI, 0.501 to 0.687), 0.631 (95% CI, 0.541 to 0.722), and 0.6 (95% CI, 0.508 to 0.691), respectively (Figure S2B). For CXCL1, CXCL5, and CXCL13 at 1-year stage, AUC was 0.651 (95% CI, 0.561 to 0.740), 0.643 (95% CI, 0.553 to 0.733), and 0.57 (95% CI, 0.473 to 0.667), respectively (Figure S2C). These results suggested the appreciable reliability of chemokines as biomarkers for GBM prognosis, and CCL18 exhibited the best accuracy among these seven chemokines.

Figure 2 Prognostic values of different chemokines in GBM. (A) Forest plot of associations of 35 chemokines expression and OS showed that seven chemokines (CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13, marked in pink) were independently associated with the OS of GBM patients. (B) The Kaplan-Meier analysis of the association between seven chemokines and DSS. (C) The Kaplan-Meier analysis of the association between seven chemokines and PFI. Abbreviations: OS, overall survival; DSS, disease specific survival; PFI, progress free interval.

Relations Between the Prognosis-Related Chemokines with Immune Infiltration Levels and Immune Checkpoints

Immune cells in the TME affect the overall survival of cancer patients.⁸ Therefore, we investigated the correlation between the prognosis-related chemokines (CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13) expression levels and immune infiltration in GBM. The results showed a positive correlation between the prognosis-related chemokines expressions with poorer prognosis and higher immune infiltration in GBM. As shown in Figure 3A, the strongest positive correlation existed between the CCL2 or CXCL5 expression and infiltrating levels in macrophages (correlation = 0.742 and 0.707, respectively; all p < 0.001). In addition, the mRNA expression levels of all seven chemokines were appreciably positively correlated with the infiltration levels in macrophages, neutrophils, iDC, DC, Th1 cells, cytotoxic cells, T cells, mast cells, and eosinophils in GBM (marked in green rectangles). In contrast, CCL28 and CXCL13 had weaker correlations with the infiltration levels of immune cells than other chemokines (Figure 3A). Because chemokine networks affect tumor immunity and tumorigenesis by regulating tumor microenvironment,²³ we next calculated the correlations of chemokines with each other in patients with GBM. The results indicated significant (correlation > 0.6 and p < 0.001) and strong positive correlations in the following prognosis-related chemokines: CCL2 and CCL8, CCL5 and CCL8, CCL7 and CCL18, CCL13 and CCL18, CXCL1 and CCL7, CXCL1 and CCL20, CXCL5 and CXCL1-3, CXCL5 and CCL2, CCL7, CCL13, or CCL20 (Figure 3B). Furthermore, we calculated the correlations between prognosis-related chemokines expression and immune checkpoint molecules using the spearman’s method. The results showed that both CCLs and CXCLs expression levels were significantly associated with most checkpoint molecules, including TNFRSF14, LAIR1, TNFSF4, CD244, ICOS, CTLA4, CD48, CD28, CD200R1, HAVCR2, and CD80 in GBM (all p < 0.05, Figure 3C). Therefore, these results further confirmed that the prognosis-related chemokines were significantly correlated with immune infiltrating cells in GBM, suggesting that chemokines might play vital roles in the GBM microenvironment.

Figure 3 Correlations between seven prognosis-related chemokines expression and immune cell pathway or immune checkpoint molecules in GBM from the TCGA dataset. (A) Correlations between the infiltration of immune cells and the expression of the indicated chemokines. (B) Correlations between different chemokines in GBM. (C) Correlation analysis of the indicated chemokines expression levels with 40 common immune checkpoint genes in GBM. Abbreviations: aDC, activated DC; iDC, immature DC; pDC, plasmacytoid DC; Tcm, T central memory; Tem, T effector memory; Tfh, T follicular helper; Tgd, T gamma delta.

Enrichment Analyses and the PPI Network Construction Among Correlated Genes

As CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13 could be potential prognostic biomarkers for GBM patients, we next identified the differentially expressed mRNAs (DEmRNAs) in GBM samples with these chemokines^high and chemokines^low expression groups as well as in GBM and adjacent normal tissues using the TCGA database with p < 0.05 and |log fold change [FC]| > 1 as the mRNA threshold. Volcano plots visually displaying the distribution of DEmRNAs were generated for seven prognosis-related chemokines, respectively (Figure 4A–C). The relative expression values of the top 15 differentially expressed genes (DEGs) between the two cohorts are shown in Figure 4B–D. Moreover, the Venn diagrams showed that 89 genes were differentially expressed in CCL2, CCL8, CCL18, and CC28 datasets (Figure 4E). 229 overlapped DEGs were identified that were differentially expressed both in CXCL1, CXCL5, and CXCL13 datasets (Figure 4F).

Figure 4 Volcano plots and heatmap plots of DEmRNAs between the expression of chemokinehigh and chemokinelow in GBM samples. (A–D) Differentially expressed genes for high expression of the prognosis-related chemokines vs low expression of the prognosis-related chemokines in GBM were shown in the volcano plots. Red represented upregulated genes, and blue indicated downregulated genes. Heatmaps of differentially expressed genes were shown in the lower panel. 15 significant DEGs were shown for each chemokine. (E) Common differentially expressed genes representation through Venn diagrams. 89 genes were found commonly from CCL2, CCL8, CCL18, and CCL28 datasets. (F) 229 genes were found commonly from CXCL1, CXCL5, and CXCL13 datasets. *p < 0.05, **p < 0.01, ***p < 0.001.

Next, the common DEGs were compared between CCLs and CXCLs datasets, which were visualized through a Venn diagram (Figure 5A). The number of common DEGs was 64, which accounted for 25.19% of a total of 254 differentially expressed genes. Also, the functional enrichment analysis (GO and KEGG) was performed for 64 common DEGs to explore the potential functions associated with the common regulatory network. It showed that the 64 DEmRNAs participating in the network were mainly enriched in the “cytokine-cytokine receptor interaction”, “leukocyte migration”, “viral protein interaction with cytokine and cytokine receptor”, “myeloid leukocyte migration”, “leukocyte chemotaxes”, and “collagen-containing extracellular matrix” (Figure 5B). Then, to determine whether the 64 genes were associated with GBM prognosis, we used the Kaplan-Meier analysis and a Log rank test to perform OS analysis of GBM patients. In total, 17 DEmRNAs (VDR, MLPH, IL6, C5orf46, RARRES1, MCEMP1, LRRC15, RFX8, ADTRP, AQP9, PF4V1, HAS1, IL7R, HTR3A, SAA2, MARCO, and CD300E) were found to be related to prognosis based on p < 0.05 (Figure S3). The subsequent functional enrichment analysis of the 17 genes revealed that they were strongly associated with “cytokine-cytokine receptor interaction”, “hematopoietic cell lineage”, “viral protein interaction with cytokine and cytokine receptor”, “FoxO signaling pathway”, and “JAK-STAT signaling pathway” (Figure 5C). In the protein-protein interaction (PPI) network, the hub genes were positively correlated with the prognosis-related chemokines (interaction score > 0.4) (Figure 5D). MARCO, VDR, LRRC15, SAA2, MCEMP1, PF4V1, IL7R, IL6, and AQP9 were highlighted in the module network as these nine genes were the prognosis-related common DEGs in GBM (Figures 5D and S3).

Figure 5 Identification of DEGs, functions, pathways, and protein-protein network (PPI) of the prognosis-related chemokines in GBM. (A) The Venn diagrams showed a total of 64 overlapped DEGs identified from the CCLs (CCL2, 8, 18, 28) and CXCLs (CXCL1, 5, 13) datasets. (B) GO and KEGG analysis showed multiple biological processes and several pathways of the overlapped DEGs. (C) Functional enrichment analysis showed five pathways of the prognosis-related overlapped DEGs. (D) Protein-protein interactions (PPIs) network identified common differentially expressed genes shared by the seven prognosis-related chemokines datasets. Nodes in red indicated the seven prognosis-related chemokines, and nodes in pink showed the prognosis-related common genes in GBM.

High CCL18 Protein Expression in GBM Tissues Led to Poor Prognosis in Glioma Patients

Among the seven prognosis-related chemokines discussed above, CCL2, CCL18, and CXCL5 exhibited stronger correlations with clinical outcomes and better prognostic accuracy for glioma patients. Therefore, we analyzed the protein expression levels and prognostic values of CCL2, CCL18, and CXCL5 in 72 glioma patient samples (Table S1). These glioma tumors included 15 Grade II (G2), 14 Grade III (G3), and 43 Grade IV (G4) samples. We found that CCL18 and CXCL5, but not CCL2, were more highly expressed in GBM (G4) than in lower grades (G2 and G3) using the ELISA method (Figure 6A–C). Next, we analyzed the protein expression levels of CCL2, CCL18, and CXCL5 with the factors including IDH1 status, gender, and the human tumor cell proliferation marker Ki67 for glioma. Higher CCL18 expression was closely overlapped with IDH1 mutation in the cohort (Figure 6A–C). Besides, the rate of Ki67 proliferation was 1–25% in 35 cases, whereas it was 30–80% in 30 cases with higher CXCL5 expression (Figure 6A–C). Further, the Kaplan-Meier curve and Log rank test analyses revealed that only increased CCL18 protein level was significantly associated with the overall survival (OS) of all patients with glioma (Figure 6D). We next calculated the correlation of three chemokines (CCL2, CCL18, and CXCL5) with each other by analyzing their protein expressions levels in our cohort, which indicated significant and positive correlations in each tested pair (Figure 6E).

Figure 6 The protein expression levels of CCL2, CCL18, and CXCL5 in 72 human glioma tissues. (A–C) Quantification of indicated protein expression levels in patients grouped by G2 & G3 gliomas and GBM (G4), IDH1 wild type and mutant, female and male, low Ki67 rate (< 30%) and high Ki67 rate (≥ 30%), respectively. The relative CCL2, CCL18, and CXCL5 protein levels of human glioma tissue lysates were firstly measured by ELISA and then normalized to that of β-actin. (D) Kaplan-Meier OS curves (High vs Low) for patients stratified by different protein levels of CCL2, CCL18, and CXCL5 in gliomas. (E) Correlation between CCL2, CCL8, and CXCL5 expression in our cohort. All data are representative of three independent experiments. *p < 0.05, **p < 0.01, ****p < 0.0001, ns, no significance.

Discussion

Chemokines participate in anticancer immune response and form a gradient that can chemoattract leukocytes to the site of damage or infection.^24–26 On the contrary, chemokines affect tumor progression through various mechanisms that directly affect cancer cell proliferation, or indirectly regulate angiogenesis and recruitment of immune cells that facilitate tumor growth and metastasis.^11,27 For example, CCL2 inhibits apoptosis of endothelial cells by directly binding receptors CXCR4 and CCR2 that expressed on tumor vessels.²⁸ CCL18 directly influences tumor cells by promoting invasion, metastasis, and EMT in breast cancer, pancreatic cancer, ovarian cancer, and prostate cancer.^8,29 Among heterogeneous primary tumors of the central nervous system (CNS), gliomas are the most frequent type, with glioblastoma multiforme (GBM) characterized with the worst prognosis.¹ Previous studies showed that some chemokines were significantly related to prognosis and might be potential biomarkers for GBM, including CXCL1, CXCL5, CXCL10, CXCL11, CXCL14, CCL2, CCL11, CCL13, CCL21, CCL22, CCL27^14,30–32 (Table S2). However, different numbers of GBM samples from the TCGA database were used to analyze the expression of chemokines, making it challenging to identify promising biomarkers for GBM. Thus, accurate survival prediction through a comprehensive analysis of chemokines is essential for GBM patients.

We here analyzed the high-throughput RNA-seq data from the TCGA database to identify 42 differential expressed chemokines genes. The results showed that the expression levels of 32 chemokines were higher than that in normal tissues, and 27 chemokines were positively correlated with tumor stages in patients with gliomas (Figure 1C and D). The increased expression levels of seven chemokines, CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13, were significantly correlated with poor OS and DSS in all of the patients with GBM, indicating that these chemokines were oncogenic events in GBM (Figure 2). Higher CCL2, CCL28, or CXCL5 expression was also found in IDH-1 wild-type GBM (Figure 1E and F). Moreover, higher expression levels of CCL8, CCL18, CCL28, CXCL1, CXCL5, or CXCL13 predicted shorter 1- and 2-year OS in GBM patients, and each chemokine was an independent predictor of prognosis in GBM (Figure S2A). Importantly, CCL18 might be the most accurate biomarker for GBM than other chemokines, and the Kaplan-Meier curve for overall survival in the Chinese Glioma Genome Atlas (CGGA) database also confirms its power for independently predicting prognosis in GBM (Figure S4).

The microenvironment (TME) plays an essential role in the progression and metastasis of GBM.³³ Generally, TME is enriched for pro-inflammatory cytokines, chemokines, and their complex cross-talk with their receptors influences the development of cancers.²³ For example, macrophages can be recruited into the tumor microenvironment by the CCL2-CCR2 signaling, associated with poor patient prognosis in breast cancer.^34,35 As a significant cytokine of TAMs, CCL18 immunopositive cells represent a subset of macrophages in breast cancer tissues.³⁶ Besides, CXCL5 binds to its receptors CXCR2 in TME to participate in the recruitment of immune cells, promoting tumor growth and metastasis.³⁷ In our study, immune cell infiltration analysis demonstrated that the prognosis-related chemokine members (CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, CXCL13) were significantly associated with different immune cells, especially with macrophages among 24 immune cells in GBM (Figure 3). Previous studies showed that high-grade murine gliomas harboring a mutant IDH1 allele exhibited reduced macrophage infiltration, which correlated with lower chemokine expression levels.³⁸ Also, CCL2 is produced by the tumor cells, which attract macrophages in experimental glioblastoma models.^38,39 Combined with our results, we proposed that the prognosis-related chemokines might be involved in the recruitment of immune cells or partly derived from macrophages, deteriorating tumor microenvironment, inhibiting prolonged survival of GBM patients.

The functional enrichment analysis showed that the 64 common gene sets associated with the seven prognosis-related chemokines in GBM were mainly enriched in cytokine-cytokine receptor interaction and leukocyte migration. Our Kaplan-Meier survival analysis revealed that 17 of 64 common gene expressions were significantly associated with a poor prognosis for OS, involved in the JAK-STAT signaling pathway and the FoxO signaling pathway (Figure 5). This result further confirmed the efficacy of CCL2, CCL8, CCL18, CCL28, CXCL1, CXCL5, and CXCL13 as prognostic indicators of GBM.

Many pan-cancer analysis studies revealed some hallmarks of different cancers, including the genetic and metabolic alterations.¹ However, these prognostic biomarkers might not be accurate without considering protein translation or post-translation profiles. We here detected the protein expression levels of CCL2, CCL18, and CXCL15 in 72 glioma tissues, and discovered a remarkably higher expression of CCL18 in the subgroup of G4 glioma compared to lower grade glioma cases (G2 and G3), accompanied by a poor prognosis of glioma patients (Figure 6). Also, the expression of CXCL5 was associated with several clinicopathological parameters in glioma patients (eg, tumor grade and high Ki67 rate), suggesting the critical role in promoting tumor growth (Figure 6).

In conclusion, we comprehensively analyzed the differential expression of four subfamilies of chemokines in GBM and evaluated their clinical and prognostic values. We demonstrated that CCL18 could be an accurate, independent prognostic factor for glioma patients, while other prognosis-related chemokines (CCL2, CCL8, CCL28, CXCL1, CXCL5, and CXCL13) might be an additional diagnostic biomarker. It can be speculated that these prognosis-related chemokines, especially CCL18, could become therapeutic targets for GBM treatment.