Comprehensive analysis of a TNF family based-signature in diffuse gliomas with regard to prognosis and immune significance

Several studies have shown that members of the tumor necrosis factor (TNF) family play an important role in cancer immunoregulation, and trials targeting these molecules are already underway. Our study aimed to integrate and analyze the expression patterns and clinical significance of TNF family-related genes in gliomas. A total of 1749 gliomas from 4 datasets were enrolled in our study, including the Cancer Genome Atlas (TCGA) dataset as the training cohort and the other three datasets (CGGA, GSE16011, and Rembrandt) as validation cohorts. Clinical information, RNA expression data, and genomic profile were collected for analysis. We screened the signature gene set by Cox proportional hazards modelling. We evaluated the prognostic value of the signature by Kaplan–Meier analysis and timeROC curve. Gene Ontology (GO) and Gene set enrichment analysis (GSEA) analysis were performed for functional annotation. CIBERSORT algorithm and inflammatory metagenes were used to reveal immune characteristics. In gliomas, the expression of most TNF family members was positively correlated. Univariate analysis showed that most TNF family members were related to the overall survival of patients. Then through the LASSO regression model, we developed a TNF family-based signature, which was related to clinical, molecular, and genetic characteristics of patients with glioma. Moreover, the signature was found to be an independent prognostic marker through survival curve analysis and Cox regression analysis. Furthermore, a nomogram prognostic model was constructed to predict individual survival rates at 1, 3 and 5 years. Functional annotation analysis revealed that the immune and inflammatory response pathways were enriched in the high-risk group. Immunological analysis showed the immunosuppressive status in the high-risk group. We developed a TNF family-based signature to predict the prognosis of patients with glioma. B-u8CpeWu2mbi9y1T26HjP Video abstract Video abstract

14.6 months [5,6]. In order to improve the prognosis of patients, it is urgent to further understand the pathogenesis of gliomas and develop novel treatment strategies, such as targeted therapy and immunotherapy. Immune checkpoint blockade therapy enhances the anti-tumor immune response by targeting the regulatory pathways in immune cells. In recent years, checkpoint blockade of the CTLA-4 or PD-1/PD-L1 axis has been proven to be an effective strategy for advanced melanoma, non-smallcell lung carcinoma, advanced renal cell carcinoma and Hodgkin lymphoma [7][8][9].
Tumor necrosis factor (TNF) family consists of 19 ligands of TNF superfamily (TNFSF) and 29 members of TNF receptors superfamily (TNFRSF), which affect many biological processes, including apoptosis, host defense, inflammation, and autoimmunity [10]. TNF is produced by many different immune and non-immune cell types, which plays an important role in the development and function of the immune system [11]. At present, in addition to blocking the immune checkpoints of the B7-CD28 family (such as PD1/PD-L1), anti-tumor immunity can also be augmented by biologics or genetic engineering techniques that modulate TNFSF/TNFRSF signaling [12,13]. Some cancer immunotherapy targets from the TNF family are very attractive and have entered the stage of clinical trials, such as 4-1BB, OX40, GITR, and so on [14][15][16].
In gliomas, some members of the TNF family have been widely studied for their role in regulating tumor genesis and growth. Yeung et al. revealed that TNF-α promoted the production of multiple inflammatory mediators by p38 MAPK signaling, thereby contributing to the expansion of GBM [17]. Ramaswamy et al. found that TNF-α enhanced the invasion ability of glioma cells through MEK-ERK signaling [18]. Shibahara et al. demonstrated that the OX40/OX40L (TNFRSF4/ TNFSF4) signaling pathway induced anti-tumor immunity in a mouse glioma model, and OX40 could also trigger regulatory T cells to cause immunosuppression under hypoxia [19]. The members of TNF family might activate or inhibit immune responses in the tumor microenvironment [13], so some members have been selected as potential targets for glioma immunotherapy. Woroniecka et al. revealed that 4-1BB (TNFRSF9) agonism reduced exhaustion of tumor-infiltrating lymphocytes and improved their function, thereby prolonging the survival of GBM in combination with anti-PD1 therapy [20]. Shoji et al. found that the local delivery of an anti-CD40 (TNFRSF5) agonistic antibody induced significant anti-tumor effects in mouse glioma models [21]. The agonist OX40 (TNFRSF4) immunotherapy combined with vaccination reversed T lymphocyte exhaustion and prolonged survival in the glioma mouse model [22]. The TNF family has shown great potential in targeted therapy. However, currently, the characteristic of TNF family-related gene set has not been systematically profiled in gliomas.
In this study, we systematically analyzed the expression patterns and clinical significance of TNF family-related genes in gliomas. We developed a TNF family-based risk signature that classified patients into the high-risk or low-risk group. There were significant differences in clinicopathological characteristics, genomic alteration, prognosis, and immune status between the two groups. Finally, we established an individualized nomogram model based on signature to predict the 1-year, 3-year, and 5-year overall survival rate of glioma patients.

Data collection
Our study collected 1749 glioma cases from four public datasets. As the training set, the Cancer Genome Atlas (TCGA, http:// cance rgeno me. nih. gov/) dataset contained RNA-seq data, somatic mutation, copy-number alterations (CNAs), clinical and pathological information of 702 glioma cases. The validation sets contained our Chinese Glioma Genome Atlas (CGGA) dataset, GSE16011 dataset, and Rembrandt dataset. In our CGGA dataset, we have collected RNA-seq data of 325 gliomas, which were generated by the Illumina HiSeq platform [23]. We also collected clinical and molecular information of CGGA patients. Our CGGA dataset was approved by the Beijing Tiantan Hospital Capital Medical University Institutional Review Board (IRB KY2013-017-01) [24]. The other two validation sets included 268 glioma cases from the GSE16011 microarray database (http:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE16 011) and 454 glioma cases from The Repository for Molecular Brain Neoplasia Data (Rembrandt, https:// www. ncbi. nlm. nih. gov/ geo/ query/ acc. cgi? acc= GSE10 8474).

Signature development
We included all well-defined TNF family genes, including 18 TNFSF genes and 29 TNFRSF genes. In the training set, univariate Cox proportional hazards regression analysis was used to screen TNF family genes related to overall survival. With selected genes, the Least Absolute Shrinkage and Selection Operator (LASSO) regression algorithm generated a Cox model with minimum average cross validation error based on tenfold cross validation [25][26][27]. The LASSO cox model included 8 genes and our signature risk score was developed with a linear combination of 8 gene expression level (expr) weighted by their LASSO regression coefficients: Risk Score = (expr gene1 × coefficient gene1 ) + (expr gene2 × coefficient gene2 ) + ⋯ + (expr gene8 × coefficient gene8 ).

DAVID functional annotation and Gene Set Enrichment Analysis (GSEA)
DAVID (https:// david. ncifc rf. gov/) is a comprehensive set of functional annotation tool for understanding the biological meaning behind gene sets. We first performed Pearson correlation analysis and screened out genes that were significantly positively correlated with signature (Pearson R > 0.6, p < 0.05). Then these genes were uploaded to DAVID for Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis. Gene set enrichment analysis (GSEA) is a widely used tool for assessing pathway enrichment with transcriptome data. In this study, we adopted FGSEA (Fast Gene Set Enrichment Analysis) method, which could estimate low GSEA P-values with a high accuracy in a short time [28]. |NES| > 1 and adjusted p-value < 0.05 were considered significant in GSEA.

Analysis of immune and inflammatory responses
The CIBERSORT algorithm allowed us to quantify the infiltrating immune cells in tumors with gene expression profiles [29]. We calculated 22 immune cell subtypes with the customized gene signature file "LM22". Differences in the proportion of immune cells between high-risk and low-risk groups were assessed by Student's t-test.

Statistical analysis
R language (v4.0.0, https:// www.r-proje ct. org/) was the main statistical analysis environment. Normalized gene expression values were log-transformed (based on 2). Univariate and multivariate Cox regression analysis were performed to evaluate the prognostic value. Differences in variables between the groups were assessed using Student's t-test or Chi-square test. The survival differences in Kaplan-Meier survival curves were evaluated by logrank test. Time-dependent ROC curve (timeROC) was used to predict one-, three-and five-year overall survival [32,33]. The nomogram model integrated signature and clinical indicators to predict 1-year, 3-year and 5-year survival probability with R package "rms" [34]. Other R packages involved in this study included pheatmap, ggl-pot2, pROC, ComplexHeatmap, Hmisc, circlize and corrgram. A two-sided test p value < 0.05 was considered statistically significant. P values were adjusted by Benjamini-Hochberg procedure (BH) in multiple hypothesis testing.

Identification of a TNF family based-signature in gliomas
After univariate Cox regression analysis, 39 genes were left for further analyzed in TCGA dataset. By LASSO regression algorithm, we screened eight genes as covariates to evaluate prognostic value (Fig. 1a). Our TNF family based-signature (risk score) was developed with a linear combination of the expression of eight genes weighted by their regression coefficients (Fig. 1b). Then, TCGA patients were divided into low-risk group and high-risk group according to the median risk score as the cutoff value. We found significant differences in clinical and pathological characteristics between the two groups ( Fig. 1c and Table 2). Patients in the high-risk group were older than those in the low-risk group (p < 0.001). And the incidence of GBMs (WHO grade 4) in the highrisk group was higher than that in the low-risk group (p < 0.001). Meanwhile, IDH-wildtype, 1p/19q non-codeletion, and MGMT promoter unmethylation were more common in the high-risk group (p < 0.001). Malignant molecular subtypes, including classical and mesenchymal subtypes, were significantly enriched in the highrisk group (p < 0.001). In the validation dataset (CGGA, GSE16011 and Rembrandt datasets), the same regression coefficients from TCGA (Fig. 1b) were used to calculate the risk score of each patient. In three validation datasets, we also observed differences in clinical and pathological  Log(λ) Partial Likelihood Deviance 39 39 38 38 36 31 28 22 19 11 9 8 7 5 4 3 2  features between high-risk and low-risk group (Fig. 1d, Additional file 1: Figure S2, Table 2 and Additional file 1: Table S1).

Correlation between TNF family based-signature and pathological features in gliomas
Since gliomas covered grade 2-4 and different molecular subtypes, we further studied the distribution characteristics of TNF family based-signature (Fig. 2a). As the WHO grade increased from 2 to 4, risk score increased significantly (p < 0.05). Meanwhile, risk scores were significantly higher in the IDH-wildtype group, 1p/19q non-codeletion group, and the MGMT promoter unmethylated group (p < 0.05). Among the three TCGA subtypes, we found that patients defined as mesenchymal subtype had the highest risk score (p < 0.05). In all three validation datasets ( Fig. 2b and Additional file 1: Fig. S3A and S3B), we also observed that the distribution of signature was consistent with the above results. Then, we utilized Receiver Operating Characteristic (ROC) curve to evaluate the predictive value of our signature for pathological indicators. Compared with age and gender, our signature showed superior predictive value in WHO grade (AUC = 0.934), IDH mutation status (AUC = 0.964), 1p/19q codeletion status (AUC = 0.850), MGMT promoter methylation status (AUC = 0.801) and mesenchymal subtype (AUC = 0.926) (Fig. 2c). Signature also showed high predictive power in all three validation datasets (Fig. 2d, Additional file 1: Fig. S3C and S3D).

Different patterns of genomic alterations between lowand high-risk gliomas
At the level of genomic alterations, somatic mutation and copy-number alterations (CNA) data from TCGA were included for further study. In Fig. 3, we found a significant enrichment of IDH1, ATRX, CIC, NOTCH1,  (Fig. 3).

Prognostic analysis of the TNF family based-signature
Next, we evaluated the clinical prognostic value of our signature using follow-up data. The Kaplan-Meier survival curve showed that the high-risk group had a significantly worse prognosis than the low-risk group in four datasets (Fig. 4, p < 0.001, log-rank test). When patients were divided into WHO grade 2, grade 3 and grade 4, we also observed significantly shorter survival in the high-risk group than in the low-risk group (p < 0.05).
In both univariate and multivariate Cox regression analyses, our signature risk score was significantly associated with overall survival in four datasets (    (Fig. 5a). In the three validation datasets, our signature also showed high time-dependent AUC ( Fig. 5b and Additional file 1: Fig. S4). Combining independent prognostic indicators (age and risk score) in the TCGA dataset, we then constructed a nomogram model to predict 1-year, 3-year and 5-year survival probability for glioma patients (Fig. 5c, C-index = 0.868). The calibration diagram showed satisfactory consistency between the nomogram model prediction and observations in survival (Fig. 5d).

Association between the TNF family based-signature and glioma immune and inflammatory response
To better understand the relationship between TNF family and glioma immunity, we included immune cells and immune checkpoints for analysis. By CIBERSORT algorithm, we calculated 22 immune cell components in glioma samples. We found that in the high-risk group, suppressive or resting immune cells were significantly increased, including macrophages M0, resting NK cells, and regulatory T cells (Fig. 6d,   lymphocytes-related inflammatory responses, but not B lymphocytes, in the high-risk group.

Discussion
As a new milestone in cancer treatment, modern immunotherapy has brought light to many tumors that previously had limited treatment options. Immunity checkpoint blockade with monoclonal antibodies targeting the B7-CD28 superfamily (CTLA-4, PD-1, and PD-L1) produces a durable anti-tumor immune response, and this strategy has been applied in many tumors and translated into clinical benefits [35]. However, some patients are resistant to immune checkpoint blockade therapy, and some patients will eventually relapse [36]. In lung cancer, the response to immune checkpoint blockade therapy is not universal, with more than half of patients showing primary resistance [37]. Among advanced melanoma that respond to CTLA-4 or PD-1 blockade therapy, approximately one-quarter to onethird of patients will relapse over time [38]. Also in gliomas, even though clinical trials of B7-CD28 superfamily inhibitors are now under active investigation, they have ended in failure so far and most patients have received little or no obvious benefit (CheckMate-143, Check-Mate-498, CheckMate-548) [39][40][41]. These results suggest that there may be other immunoregulatory signaling pathways in the tumor microenvironment.
In fact, there are many pathways to protect tumor cells from immune damage, which is the recognition of self by immune cells, leading to resistance to the first generation of immune checkpoint blockers [42]. Therefore, it is necessary to develop new immunostimulatory targets to overcome the primary and acquired resistance of immunotherapy. In addition to the B7-CD28 family, the TNF family also contains many immune checkpoints, such as OX40, 4-1BB, GITR, and so on. In the process of antigen recognition, OX40 (TNFRSF4) is induced to express on activated T cells. Agonists targeting OX40 can provide powerful co-stimulatory signals, thereby enhancing the expansion and proliferation of CD4 + and CD8 + T cells that recognize tumor antigens. Several agonistic antibodies targeting OX40 are currently undergoing cancer clinical trials [15] and combining OX40 antibody and other immune checkpoint antibodies is more effective than monotherapy [43]. 4-1BB (TNFRSF9) is another attractive cancer immunotherapeutic target in the TNF family, which is a co-stimulatory receptor expressed on activated T cells and NK cells. The preclinical results of tumor models indicate that the currently developed agonist antibodies targeting 4-1BB can clear tumors and maintain durable anti-tumor immunity [44,45]. Agonistic monoclonal antibodies targeting 4-1BB have shown good results in patients with lymphoma, and are undergoing combined therapy trials with other immunomodulators [14]. Due to the inflammatory hepatotoxicity of the first-generation 4-1BB antibody (urelumab), new 4-1BB agonists are under clinical development, seeking to maximize immune activation and avoid liver inflammation side effects [46]. In gliomas, immune checkpoints from the TNF family also serve as a co-stimulatory signal in regulating immunity. Woroniecka et al. found in mouse models that 4-1BB agonists eliminated the limitations of poor T cell activation and severe exhaustion, and combined anti-4-1BB and anti-PD1 treatments provided survival benefits [20]. Nusrat Jahan et al. found that the agonist anti-OX40 was effective for intracranial glioma and prolonged survival time in a mouse model of glioma [22]. In this study, we comprehensively analyzed the expression characteristics and prognostic significance of TNF family in gliomas. Pearson correlation analysis showed that the RNA expression of most TNF family genes was positively correlated. At the same time, univariate Cox regression analysis showed that the expression of 39 TNF family genes was significantly associated with overall survival. These results indicated that the TNF family members were closely related and had potential clinical value. Then using Lasso regression model, we developed    TNF family genes: TNFSF4, CD70, TNFSF14, TNFRSF19,  NGFR, TNFRSF11B, TNFRSF14, and TNFRSF12A. Based on the median risk score, we divided patients into highrisk groups and low-risk groups, and assessed the differences between the two groups. We found that patients with older age, WHO grade 4, IDH-wildtype, 1p/19q intact, MGMT promoter non-methylated, and mesenchymal subtype were more common in the high-risk group. Meanwhile, we found a significant correlation between risk score (RS) and clinical molecular features ( Figure S6, Spearman correlation, p < 0.05). At the level of genomic variation, we found that mutations in EGFR, NF1, PTEN, and RB1 were significantly enriched in the high-risk group. And in the high-risk group, we found more amplification regions such as EGFR, CDK4, PDGFRA, MDM2, and deletion regions such as CDKN2A, CDKN2B, MLLT3, PTEN. These differences suggest that TNF family-based signature may be associated with the malignant biological process, as well as poor prognosis in glioma patients. Next, survival curve analysis confirmed that patients in the high-risk group had a worse prognosis than low-risk group. And univariate and multivariate Cox regression analyses identified our signature as an independent prognostic indicator after adjustment of other clinicopathological factors. Then we used ROC curve to evaluate the survival predictive value of signature, which was superior to the traditional indicators (age and grade). Based on the superior predictive ability of signature, we combined with age to construct a nomogram survival prediction model. This model had good clinical application value in predicting the 1-year, 3-year and 5-year survival rates of individuals with gliomas.
In order to explore the potential biological mechanism of our signature, we performed DAVID functional annotation and GSEA enrichment analysis. We found differences in immune and inflammatory responses between the high-risk and low-risk groups, revealing the correlation between the TNF family and the immune microenvironment of gliomas. Then the immune cell infiltration analysis showed that the high-risk group had more suppressive or resting immune cells in, including macrophages M0, resting NK cells, and regulatory T cells. Meanwhile, signature was also positively correlated with the expression of immune checkpoints (PD -L1, PD1, LAG3, CTLA4, B7-H3, IDO1, CD80, TIM-3). These results all suggested the immunosuppressive status of the high-risk group. From these we could infer the close connection between the TNF family and the tumor immune microenvironment. TNF family members coordinately drove co-stimulation or co-inhibition of the immune response [47]. TNF family members usually exhibited the pro-inflammatory properties that were partly due to the activation of NF-kB signaling [48]. In addition, TNF members could activate immunosuppressive cells (regulatory T cells and myeloid-derived suppressor cells) through TNF receptor 2 (TNFR2), thus supporting immune escape and promoting tumor cell proliferation [49]. Although TNF members was initially found to mediate anti-tumor effects, recent studies have shown that they also promoted tumor progression. Lei et al. found that TNF-α treatment promoted the proliferation of glioma cells [50]. Wei et al. found that TNF-α secreted by macrophages could activate endothelial cells and promoted GBM angiogenesis [51]. In addition to the dual effects of TNF members, some members had synergistic effects. For example, CD27 (TNFRSF7), HVEM (TNFRSF14), 4-1BB (TNFRSF9) and OX40 (TNFRSF4) all had co-stimulatory effects on T cells, and the regulation of these co-stimulators might prolong T cell response and control the survival of T cells [52]. CD27 and HVEM expressed on resting T cells functioned early after initial activation of T cells, while OX40 and 4-1BB signals on T cells were delayed relative to initial activation and showed preferential effects on CD4 and CD8 T cells. The TNF family had diverse and complex interactions with the immune system [53], so the global evaluation of the TNF family was of great significance. Through bioinformatics analysis, we have systematically explored the TNF family-related genes and their potential relationship with immunity, but our study is preliminary and requires further in-depth analysis and biological experiment verification.

Conclusions
In conclusion, this is the first study in the expression profile and clinical prognostic significance of TNF family members in gliomas. We also identified a TNF familybased signature to stratify the risk of glioma patients. Our research contributes to the individualized prognostic management of glioma patients, and provides evidence for immunotherapy targeting TNF family members.