GDF15 expression in glioma is associated with malignant progression, immune microenvironment, and serves as a prognostic factor

Abstract Aims Growth differentiation factor 15 (GDF15) is involved in lots of crucial inflammatory and immune response. The clinical and immune features for GDF15 in glioma have not been specifically investigated so far. Methods Gene expression profiles obtained from public glioma datasets were used to explore the biological function of GDF15 and its impact on immune microenvironment. Interference with GDF15 in several glioma cell lines to verify its functions in vitro. Survival data were used for the survival analysis and establishment of a nomogram predictive model. Results GDF15 was up‐regulated in various malignant phenotypes of glioma. Function analysis and in vitro experiments revealed that GDF15 was associated with malignant progression and NF‐κB pathway. GDF15 was closely correlated to inflammatory response, infiltrating immune cells, and immune checkpoint molecules, especially in lower grade glioma (LGG). High expression level of GDF15 predicted poor survival in LGG, while the effect on glioblastoma (GBM) was not significant. A nomogram predictive model combining GDF15 and other prognostic factors was constructed and showed ideal predictive performance. Conclusions GDF15 could serve as an interesting prognostic biomarker for LGG. Regulating the expression of GDF15 may help solve the dilemma of immunotherapy in glioma.


| BACKG ROU N D
Glioma is the most common primary brain tumor of the central nervous system. 1 Despite receiving standard therapeutic regimens, the prognosis of glioma is still unsatisfactory. 2 Identification of key driven factors governing glioma progression could not only provide in-depth understanding of the disease but also may predict prognosis and serve as target for precise treatment. Meanwhile, the success of immunotherapy in solid tumors has brought new opportunities for glioma treatment. 3 However, the current relevant clinical trials in glioma have not shown encouraging results. 4,5 The unique and complex immune microenvironment of glioma is the crucial obstacle to immunotherapy. 6 Thus, it is also urgent to figure out the potential mechanism and factors that affect the immune microenvironment of glioma, which may help improve efficacy of immunotherapy.
Growth differentiation factor 15 (GDF15), a divergent member of the transforming growth factor-beta (TGFβ) superfamily, has received appreciable attention in the last two decades for its multiple key role in several diseases, including obesity, cachexia, and cardiovascular disease. [7][8][9] Nowadays, GDF15 is regarded as a marker of oxidative stress and an inflammation-induced central mediator of tissue tolerance. 8,10 Recently, several studies have demonstrated the close links between GDF15 and the development of cancers. [11][12][13] GDF15 performs pleiotropic functions in tumor progression, acting either as suppressor or as promoter in early and late stages of tumors, respectively. 14,15 It has been shown that GDF15 may influence glioma cell invasion, and elevated level of GDF15 in the cerebrospinal fluid is associated with worse outcome of GBM patients. 16,17 However, the clinical and immune features for GDF15 in glioma have not been specifically investigated so far.
In this study, we found that GDF15 correlated with malignant progression in glioma using clinical and transcriptome (RNA-seq) data from two large glioma datasets, including The Cancer Genome Atlas (TCGA) and Chinese Glioma Genome Atlas (CGGA) datasets.
The biological characteristics and functions of GDF15 were identified based on RNA-seq and verified in cell experiments. Besides, we showed that GDF15 closely related to immune response, immune infiltration cells, and immune checkpoint molecules, especially in LGG, highlighting its important role in immune microenvironment.
Survival analysis revealed GDF15 predicted poor survival for LGG and could be served as a novel prognostic biomarker.

| Data collection
Gene expression and glioma survival data in TCGA were obtained from GlioVis (http://gliov is.bioin fo.cnio.es/), and CGGA were obtained from http://www.cgga.org.cn/. A total of 1362 samples of RNA-seq data were included. RNA-seq data from different datasets were uniformly log2 transformed. Samples without complete gene expression profiles were first excluded.

| Cell culture
Human astrocytes (HA) cell line was obtained from ScienCell Research Laboratories. Human glioma cell lines such as U87, T98G, LN229, and U251 were purchased from Chinese Academy of Sciences Cell Bank. Cell lines were cultured in DMEM with 10% fetal bovine serum and 1% penicillin/streptomycin. All cells were maintained in an incubator containing 5% CO2 at 37°C. Cells were passaged when reaching 90% confluence.

| Small interference RNA (siRNA) transfection and validation
According to the product instructions (Ruibo company), siRNA was transfected at an siRNA concentration of 100nM. Twenty-four hours after transfection, medium without siRNA was added to the cells. Western blot analysis was performed with rabbit anti-GDF15 antibody (1:1000, cat. no.8479; CST). Goat anti-rabbit IgG-HRP (1:5000, ab6721, Abcam) was used as secondary antibodies and GAPDH (1:1000, cat. no.5174; CST) was used for loading control. Total RNA from cells was extracted using TRIzol reagent (Invitrogen). The Reverse Transcription System kit (Takara Bio, Inc.) was used to perform RNA reverse transcription reactions. The real-time (RT) quantitative polymerase chain reaction (PCR) was performed using the SYBR Green real-time PCR kit (Takara). The forward and reverse primers sequences (5'-3') of GDF15 and β-actin were given in the following: GDF15-Forward, ACCTGCACCTGCGTATCTCT; GDF15-Reverse, CGGACGAAGATTCTGCCAG; β-actin-Forward, TGGC ACCCAGCACAATGAA; and β-actin-Reverse, CTAAGTCATAGTCCGCC TAGAAGCA. The relative quantitative data of mRNAs were normalized to β-actin and quantified using the 2 −ΔΔCt method.

| Western blot
Cells were lysed with protein lysis buffer containing a cocktail of protease and phosphatase inhibitors (Roche

| Cell scratch assays
U251 and LN229 cells were seeded in six-well plates (5 × 104 cells per well) and transfected with GDF15 siRNA or a NC siRNA. When the cell grew to 100% confluency, sterile 200 μl pipette tip was used to scratch and washed three times with PBS to remove the floating cells. The cells were cultured with serum-free medium and incubated in a 37°C, 5% CO 2 incubator. Cells were observed by microscopy and photographically recorded (0 h, 24 h, and 48 h cell scratches).

| Cell migration assays
Twenty-four-well transwell chambers (8 μm pore size, Corning) were used for the assay. After 24 h transfection with siRNA, 1 × 105 Cells in 100 μl serum free medium were seeded in upper chambers. 600 μl medium containing 10% FBS was added to the lower chamber. After 16 h, the migration cells were fixed using 4% paraformaldehyde, stained with 0.05% crystal violet, and counted from five random fields.

| Bioinformatics analysis
Wilcoxon test was used to investigate the expression pattern of apps.io/timer/) was employed to correct the effect of tumor purity on the expression of genes and performed correlation analysis. 23 The correlation between GDF15 and immune checkpoint molecules was evaluated by Pearson's correlation test and converted to chord diagram using circlize package (https://github.com/joker goo/circlize).

| Survival analysis and construction of predictive model
Patients with missing important clinical data were further eliminated here. A total of 588 patients from TCGA and 508 patients from CGGA were selected (Table S1). Survminer package (https://github. com/kassa mbara/ survm iner) was used for Kaplan-Meier survival analysis, and log-rank test was used to evaluate the difference between survival curves. Univariate and multivariate Cox regression analyses were used to identify the independent prognostic factors for glioma. The 5 years survival rate was selected as a representative to construct the nomogram using hdnom package (https://github. com/nanxs tats/hdnom/). A risk classification system was subsequently developed based on the nomogram. For data from cell experiments, test for normality was performed using Shapiro-Wilk test. Data followed a normal distribution were analyzed via Student's t-test and one-way ANOVA test. Otherwise, data were analyzed via Mann-Whitney test. The continuous variables were showed as mean ± SD. A two-sided p-value <0.05 was considered to be statistically significant.

| GDF15 was up-regulated in malignant phenotypes of glioma
We first compared the GDF15 expression levels across different grades in TCGA and found that GDF15 was positively correlated with tumor grades ( Figure 1A). The result was validated in CGGA ( Figure 1B). Since the status of IDH and chromosome 1p19q had important influence on glioma, the relationship between GDF15 and them was also analyzed. The results showed that patients with high GDF15 expression had higher proportion of IDH wild-type and 1p19q non-co-deletion, which were all signs of poor prognosis ( Figure 1C). Similar results could also be observed after distinguishing patients into LGG and GBM (Figure S1A-D). Subsequently, we found that classical and mesenchymal subtypes with poor prognosis had higher GDF15 expression compared with neural and proneural subtypes ( Figure 1D). Moreover, high GDF15 expression was significantly related to EGFR and chromosome 7 copy number amplification, as well as PTEN and chromosome 10 copy number reduction, which were both markers of poor prognosis for glioma ( Figure 1E).
Collectively, GDF15 was up-regulated in various malignant phenotypes of glioma, suggesting that it may affect glioma prognosis.

| GDF15 affected NF-κB related pathways and progression in glioma
To elucidate the biological functions of GDF15 in glioma, we performed GSEA in TCGA and CGGA (Table S2). The results showed  Figure S2G, H). We found that the migration ability of glioma cell lines was apparently suppressed after transfection with GDF15-siRNA. Overall, these findings demonstrated that GDF15 affected NF-κB related pathways in glioma, and played an important role in tumor progression.

| Specific relationship between GDF15 and inflammatory response
Our above finding suggested GDF15 was also closely correlated with inflammatory response. For the sake of studying the specific relationship between GDF15 and inflammatory response, we selected seven metagene clusters, which were previously estab-

| GDF15 predicted poor survival in lower grade glioma
To explore the prognostic value of GDF15 in glioma, Kaplan-Meier survival analysis were performed. We found that high expression level of GDF15 was significantly associated with worse overall survival (OS) of glioma both in TCGA (p < 0.0001) and CGGA (p < 0.0001) ( Figure 5A,D). The same association was observed while analyzing only the LGG with p < 0.0001 and p < 0.0001 for TCGA and CGGA, respectively ( Figure 5B,E). However, no significant correlation was found in OS of GBM in TCGA (p = 0.14) or CGGA Univariate and multivariate Cox regression analysis for LGG showed that GDF15, grade, and 1p19q status were significantly related to prognosis both in two datasets (Table 3). Since the p-value (p = 0.058) of the correlation between IDH status and prognosis in CGGA was close to the statistical threshold, combined with its important biological significance, it also needed to be considered as a prognostic factor. No significant association between GDF15 and OS was found in subsequent Cox analysis for GBM (Table S2). Taken together, these results revealed that GDF15 could be served as interesting prognostic biomarker for LGG.

| Construction and validation of a GDF15related prognostic nomogram
For the sake of taking advantage of the prognostic value of GDF15 in LGG, we constructed a GDF15-related nomogram and risk classification system for predicting survival. The cases from TCGA were used as primary cohort to construct the predictive model, and the cases from CGGA were chosen as validation cohort for verification of the model. Based on above results, the prognostic nomogram model was established ( Figure 6A). The model performance was assessed internally by time-dependent AUC (area under the receiver operating characteristic curve) ( Figure 6B). In the prediction of early survival rate (<3 years), the predictive ability of model was exceptional, and the AUC was all greater than 0.8. Over time, the AUC declined slightly but still remained at around 0.7. The calibration plot for the probability of survival at 5 years showed an admirable agreement between the prediction and observation in primary cohort ( Figure 6C). A risk classification system was subsequently developed based on the nomogram. After dividing the primary cohort into three risk groups, it can be seen from the Kaplan-Meier survival curves that the prognosis of low-risk group was significantly better than that of medium and high-risk groups (log-rank p < 0.001) ( Figure 6D). Similar to primary cohort, the predictive model could predict survival well in the validation cohort, and the AUC remains above 0.7 in all time periods ( Figure 6E). The calibration plot for validation cohort showed a valuable agreement between the prediction and observation in 5 years survival ( Figure 6F). Low-risk group of validation could also achieve favorable OS compared to medium and high-risk groups (log-rank p < 0.001) ( Figure 6G). Collectively, these findings further proved the prognostic value of GDF15, and combined with other prognostic factors could better predict the prognosis of LGG.

| DISCUSS ION
Traditional treatment for glioma did not achieve a breakthrough development in last decade. Although the optimization of traditional treatment showed certain advantages, further relevant

F I G U R E 4 Correlation of GDF15 expression with relatively abundance of immune cells in LGG (A-F) and GBM (G-L) of TCGA dataset.
Pearson's correlation test was used to calculate the correlation coefficient (r) and p-value, and generated regression lines fitted to each dot plot clinical researches are needed for verification. [24][25][26] As one of the most promising approaches in oncology field, immunotherapy have shown encouraging success for various types of cancer treatment. 3 However, the current prospects of immunotherapy for glioma are not optimistic. [4][5][6] The hindrances to make immune checkpoint inhibitors effectively therapeutic in glioma include

Immune Cells Gene markers
LGG GBM  GDF15 has been demonstrated its multiple important roles in several diseases comprising obesity, cachexia, and cardiovascular disease. [7][8][9] With the deepening of research, GDF15 has been shown to be closely related to the occurrence and development of some types of tumors. [11][12][13] For glioma, several studies have shown that GDF15 is associated with many cellular mechanisms, such as apoptosis, migration, and invasion. 16,29,30 Although GDF15 is closely related to the immune system and is considered to be the center of inflammation-induced tissue tolerance, 10   differed in the expression of oncogenic transcriptional programs. 28 In addition, the glioma gene expression data collected by the TCGA and CGGA datasets come from different regions, and the ethnic composition is quite different, which may further deepen the heterogeneity between GBM populations. However, there was no relevant research to analyze the difference of immune landscape in GBM between TCGA and CGGA. We would conduct deeper research in the future to verify our views.
In recent years, improvements in sequencing technology were facilitating increasingly deep studies of glioma. Researchers have reported that many genetic biomarkers, such as CAMKK2, 35 GINS2, 36 FOXD1-AS1, 37 CpGs methylation, 38 TRIB2, and MAP3K1, 39 have significant impact on the prognosis of glioma and could be used as prognosis predictors of glioma. These predictors were of great significance for in-depth understanding of the disease mechanism of glioma. In our study, Kaplan-Meier survival analysis and COX regression analysis revealed that high GDF15 expression level predicted poor overall and progression-free survival in LGG, highlighting its potential to be served as a novel prognostic biomarker.
Based on these findings, we established a GDF15-related nomogram and validation (E) cohorts. The calibration curve displayed the difference between the prediction of model of 5 years survival and the actual survival outcome in primary (C) and validation (F) cohorts. The risk classification system based on nomogram divided LGG patients into three risk groups and showed distinct survival curves in primary (D) and validation (G) cohorts predictive model and divided LGG into different risk groups. For high-risk group, more aggressive treatment may help prolong its unfavorable prognosis. On the contrary, the low-risk group may consider receiving more conservative treatment to improve the overall quality of life.

| CON CLUS IONS
In conclusion, this study confirmed the important role of GDF15 for the malignant progression and immune microenvironment of glioma, especially in LGG. Regulating the expression of GDF15 may help solve the dilemma of immunotherapy in glioma. GDF15 was significantly correlated to poor survival of LGG, suggesting its potential to be served as a novel prognostic biomarker.

CO N FLI C T O F I NTE R E S T
The authors declare no conflict of interest.

AUTH O R CO NTR I B UTI O N S
S Du, C Ren and L Guo conceived and designed the study. L Guo, N Tang, R Liu and Y Qin performed data collection and assembly. Y Chen, S Hu, and L Gao performed experiments in vitro. L Guo and Y Chen analyzed, interpreted and visualized data, and were major contributors in writing the manuscript. All authors read and approved the final manuscript.

DATA AVA I L A B I L I T Y S TAT E M E N T
TCGA dataset is available at GlioVis (http://gliov is.bioin fo.cnio. es/), and CGGA dataset is available at CGGA (http://www.cgga. org.cn/). Web tools and R packages used in this study are included in article.