The Effect of Necroptosis-related Gene Signature on the Progression, Prognosis, and Immune Microenvironment of Low-grade Glioma

Despite the incorporation of various clinical and molecular criteria in the diagnosis and prognosis prediction of low-grade glioma, individual variation and risk stratication have not been completely explored. Necroptosis is considered closely related to different types of cancers, including low-grade gliomas. In this study, we obtained the necroptosis genes from the Kyoto Encyclopedia of Genes and Genomes website, extracted necroptosis genes from The Cancer Genome Atlas, and established a necroptosis-related gene signature (NECSig) through hazard analyses. Then we established a prognostic risk model consisting of four NECSig (BID, H2AFY2, MAPK9, and TNFRSF10B). Result Based on the model, the high-risk group is signicantly associated with poorer overall survival. The accuracy of this model is further supported by the receiver operating characteristic curve. Then, we constructed a prognostic nomogram combining NECSig and clinical features, which shows good predictive power for stratication of survival risk. We discovered variations in the kind of immune inltration, immune cells, and functions between the high-risk and low-risk groups using this risk model. We also showed that drug therapy is more sensitive in high-risk populations. The results revealed a prognostic indicator of NECSig, which may provide information for immunological research and treatment of low-grade gliomas.


Introduction
Diffuse glioma is the most common type in terms of malignant brain tumor [1]. Glioblastoma, or grade IV glioma, is the most aggressive kind, with a ve-year overall survival (OS) rate slightly around 5% [2].
Meanwhile, low-grade glioma (LGG) is generally a less aggressive tumor with a median survival rate of 5-10 years [3,4]. Total tumor excision, postoperative radiation, chemotherapy, and other comprehensive therapies are among the traditional therapeutic options. In recent years, many studies focused on integrating clinicopathological factors as prognostic biomarkers. However, a comprehensive knowledge of the pathomechanism and a more de ned prognostic model remain critical for improving the overall survival and quality of life in LGG patients.
Necroptosis, which depends mainly on the activity of receptor-interacting protein 1 (RIPK1) and receptorinteracting protein 3 (RIPK3) kinases, is a type of cell death that occurs in the absence of apoptosis [5,6].
The death signal induces the activation of RIPK3 kinase, phosphorylates the speci c executive protein mixed lineage kinase domain-like protein (MLKL) for cell necroptosis [6,7]. Phosphorylated MLKL may then be controlled by oligomerization and translocates to the cell membrane structure, resulting in the degradation of cell membranes and organelle membranes, cell death, and intracellular material leakage [8][9][10].
In recent years, the signaling pathway of necroptosis has been widely studied. It has some common upstream signal elements with apoptosis, among which tumor necrosis factor receptor 1 (TNFR1) is the most deeply studied [11,12]. TNFα binds to TNFR1 on the plasma membrane, causing the TNF receptorassociated death domain (TRADD) to signal RIPK1 and recruit RIPK3 to form a necrosome [13]. Caspase-8 is a key inhibitor of necroptosis as it may cleave and inactivate RIPK1 and RIPK3 [14]. If caspase-8 is inhibited, RIPK1 and RIPK3 interact with the RIP isotype domain, triggering necroptosis. This results in a cascade of automated and cross-phosphorylation of RIPK1 and RIPK3 [15][16][17]. Studies have demonstrated that some drugs are able to induce necroptosis in glioma cells through multiple molecular mechanics [18][19][20]. However, there is still no complete understanding of the mechanics and prognostic roles of necroptosis in glioma.
In this work, we built a prognostic model based on necroptosis-related mRNA features (including BID, H2AFY2, MAPK9, and TNFRSF10B), and comprehensively examined the relationship between the model and the prognostic features of LGG patients. The model has strong predictive performance, according to survival analysis and the receiver operating characteristic (ROC) curve. Multivariate Cox and strati ed analyses suggest that this necroptosis-related gene signature (NECSig) feature acts as an independent prognostic factor of OS. To estimate the survival rate of these individuals, we created a nomogram that combined NECSig features with clinical parameters. The immunological statuses of the high-risk and low-risk group, as well as their responses to immunotherapy were compared. This four-gene-NECSigbased risk pro le might be used to predict the prognosis and immunotherapy response of LGG patients.

Identi cation of NECDEs and Functional Enrichment Analysis
Among 161 necroptosis-related genes extracted, a total of 12 genes were upregulated and 13 genes were downregulated when comparing the tumor tissues and normal tissues according to the following criteria: false discovery rate (FDR) < 0.05 and |log2(fold change) | > 1 ( Figure 1A, Table S2). The box plot depicts the expression levels of each differential gene in normal and malignant tissues ( Figure 1B). We then performed GO and KEGG enrichment analyses. As for biological process enrichment, the NECDEs are signi cantly involved in extrinsic apoptotic-signaling pathway, regulation of cytokine-mediated signaling pathway, regulation of response to cytokine stimulus, regulation of cysteine-type endopeptidase activity involved in apoptotic process, programmed necrotic cell death, necrotic cell death, activation of cysteinetype endopeptidase activity involved in apoptotic process, positive regulation of cysteine-type endopeptidase activity involved in apoptotic process, positive regulation of cysteine-type endopeptidase activity, and necroptotic process ( Figure 1C). The top 10 KEGG pathways from enrichment analysis of the NECDEs were necroptosis, in uenza A, lipid and atherosclerosis, tuberculosis, measles, Hepatitis B, Coronavirus disease -COVID-19, toll-like receptor signaling pathway, NOD-like receptor signaling pathway, and Kaposi sarcoma-associated herpesvirus infection ( Figure 1D). With the median risk score as the cutoff, there was a clear distinction between the high-risk and low-risk groups ( Figure 2B). Each sample survival status is shown via a scatterplot ( Figure 2C). The heat map shows how the NECSig gene is expressed in each sample ( Figure 2D). A notable difference in OS was detected between the high-risk and low-risk groups. High-risk patients have a much lower survival rate than low-risk individuals (p < 0.001; Figure 2E). Furthermore, the ndings show that the model has excellent sensitivity and speci city in LGG survival prediction using area under the curves (AUC). of 0.839 for one year, 0.815 for two years, 0.835 for three years, 0.688 for ve years, 0.713 for seven years, and 0.784 for nine years ( Figure 2F). Evaluation of NECSig as an independent prognostic factor for LGG We computed the hazard ratios (HRs) (95% CI) for the risk score to see if NECSig operates as an independent prognostic factor for LGG. In the univariate Cox regression analysis, the risk score was 1.286 (1.240−1.354) (p < 0.001; Figure 3A), with a score of 1.234 (1.176−1.294) in the multivariate Cox regression analysis (p < 0.001; Figure 3B), suggesting that NECSig is an independent prognostic predictor. Furthermore, the time-dependent ROC curve for the risk score was 0.828 ( Figure 3C), which was considerably greater than that for age, sex, radiation, medication, and tumor grade. This result indicates that the risk score based on the NECSig can predict survival in LGG patients.
(A) The univariate and (B) the multivariate Cox regression analyses of relations between clinical factors (including NECSig) and OS. (C) Time-dependent ROC curve analysis for predicting OS after three years based on risk score, age, gender, radiation, pharmaceutical therapies, and tumour grade. AUC: area under the curve.

Correlation between risk score and clinicopathological factors
We assessed the risk score and the phase between clinical pathological variables to better determine NECSig function in LGG development. There was a different link between risk score and several pathological variables, as indicated in the picture (p < 0.05). The risk score of patients older than 45 years was signi cantly higher than that of patients less than 45 years old (p = 8.618e−06; Figure 4A). The risk score of deceased patients was signi cantly higher than that of surviving patients in the high-risk and low-risk groups (p = 5.385e−07; Figure 4B). Similarly, Figure 4C shows that patients with grade III tumours had a considerably higher risk score than patients with grade II tumours (p = 9.512e−09). It was also found that the risk scores of patients receiving radiotherapy (p = 0.004; Figure 4D) and pharmaceutical treatments (p = 0.037; Figure 4E) were higher than those without radiotherapy and pharmaceutical treatments, which may be related to the stage of tumor. This result suggests that patients with severe clinical manifestations may have a higher risk score.
The necroptosis-related mRNA in the cohorts strati ed by (A) age, (B) survival state, (C) tumor grade, (D) radiation, and (E) pharmaceutical treatments.

Construction of a predictive nomogram
The clinically adaptable nomogram incorporating clinicopathological characteristics (age, gender, tumor grade) and the risk score was created in order to estimate the one-, two-, three-, ve-, and seven-year survival probabilities in patients with LGG ( Figure 5A). The calibration plots the nomogram for one, two, three, ve, and seven years suggest that the nomogram estimate of mortality was close to the actual mortality. and that our nomogram was suitable for clinical application ( Figure 5B Tumor microenvironment (TME) analysis Immune cells interact with tumour cells and with each other in the tumour microenvironment. Six algorithms, speci cally TIMER, CIBERSORT, quanTIseq, MCP-counter, xCell, and EPIC, were applied to estimate the abundance of immune cells in ltrating among high-risk and low-risk groups (p < 0.05; Figure  6A, Table S3). Comparative analysis of immune cells showed that myeloid dendritic cells (DCs), common lymphoid progenitors, macrophage M2s, macrophages, and neutrophils were the top ve different immune cells between the two risk groups. Given the importance of checkpoint-based immunotherapy, we then looked at the differences in immune checkpoint gene expression between the two groups. We found that most immune checkpoints have higher expressions in high-risk patients, such as PDCD-1 (PD-1), CTLA4, and CD274 (p < 0.05; Figure 6B).
Based on ssGSEA of TCGA-LGG data, a correlation assessment of immune cell subpopulations found that the scores of immune cells, including B cells, CD8+ T cells, DCs, induced DCs (iDCs), macrophages, neutrophils, pDCs, T-helper cells, Th2 cells, TILs, and Tregs were signi cantly higher in high-risk groups, while the score of NK cells was signi cantly lower (p < 0.05; Figure 7A). There were also substantial disparities in immune function scores including APC co-inhibition, APC co-stimulation, CCR, checkpoint, cytolytic activity, HLA, in ammation−promoting, MHC class I, parain ammation, T-cell co-inhibition, T-cell co-stimulation, type-I IFN response, and type-II IFN response between the low-risk and high-risk groups. All immune function scores were higher in the high-risk group (p < 0.05; Figure 7B). Six kinds of immune in ltrates in human malignancies have been discovered, ranging from tumour promoting to tumour suppressive. The six types are C1 (wound-healing), C2 (INF-γ dominant), C3 (in ammatory), C4 (lymphocyte-depleted), C5 (immunologically-quiet), and C6 (TGF-β dominant). Interestingly, the expression of immune cells and immune function is higher in high-risk patients (p < 0.05; Figure 7C). We also discovered that the immune in ltrates of the high-risk and low-risk groups in the model varied considerably, with immune in ltration type C5 being the most prominent (p < 0.05; Figure 7D). In summary, our ndings suggest that NECSig is, to some extent, connected to immune cell in ltration.

Differences of (A) immune cells and (B) immune function in high-risk and low-risk groups. (C)
Association of four NECSig expressions with different immune in ltrate subtypes tested with ANOVA in TCGA-LGG. (D) Heat map showing the distribution of LGG immune in ltrate subtypes between the lowrisk and high-risk groups. NS, not signi cant; *, p < 0.05; **, p < 0.01; ***, p < 0.001. C3: in ammatory, C4: lymphocyte-depleted, C5: immunologically-quiet, and C6: TGF-β-dominant. Drug treatment differences in the prognostic model Furthermore, we compared the sensitivity of numerous medicines in the high-risk and low-risk categories. The high-risk group had a distinctly lower IC 50 value of the drugs including Bortezomib (p < 2.22e−16), Cyclopamine (p = 4.6e−16), Dasatinib (p < 2.22e−16), Rapamycin (p < 2.22e−16), and Roscovitine (p = 7.4e−12) ( Figure 8A-E), indicating those patients may be more sensitive to these chemotherapy and targeted treatments. We used tumor immune dysfunction and exclusion (TIDE) to evaluate the potential clinical e cacy of immunotherapy between the two groups. Since a higher TIDE score represents a higher likelihood of antitumor immune escape, indicating less likely to bene t from immune checkpoint blockade (ICB) treatment. We calculated the TIDE score for each patient, which ranged from -0.5 to 1.23. The TIDE scores of patients in the high-risk group were higher than those in the low-risk group ( Figure 8F), indicating the low-risk group could be more sensitive to ICB therapy. Taken together, the NECSig may be a potential predictor of drug sensitivity.

Discussion
Extracellular stimuli are known to initiate necroptosis and cause in ammation and cell death. When cell apoptosis fails, necroptosis becomes a "failure prevention mechanism" [27]. Recent animal studies have shown that this mechanism can regulate the number of T-cells in peripheral tissues and is also necessary for clearing abnormal lymphocytes during T-cell development [28]. Studies have found that some diseases are also related to the abnormal regulation of necroptosis, such as breast cancer, ovarian cancer, glioma, lung cancer, and leukemia [29][30][31][32][33].
In this study, we constructed four necroptosis-related genes (BID, H2AFY2, MAPK9, and TNFRSF10B) to build a prognostic model. The Kaplan-Meier curve indicates that the high-risk group's survival rate is much lower than that of the low-risk group, and the ROC curve indicates that the model has strong predictive ability. There are also signi cant differences in clinical characteristics between two groups. In addition, our nomogram showed the potential to improve clinical decision-making and guide the formulation of treatment strategies. A detailed investigation found that the model we created is also related to the type of immune in ltration, high-risk immune checkpoints, and the expression of immune cells and immune functions. Therefore, this research identi es prospective indicators and therapeutic targets of necroptosis signalling pathways.
Four genes in the model are widely involved in tumor progression. BID plays a major role in the process of cell apoptosis and is related to the growth of speci c types of gastric cancers [34]. The BID gene is also one of our prognostic model markers [35]. Kun Zhang found that H2AFY2 is a poor prognostic gene in gliomas [36]. MAPK9 is upregulated in normal tissues. In NSCLC tissues, CircRNF20 aggravates the progression of non-small cell lung cancer by activating MAPK9 [37]. YIPF2 can enhance the recirculation of TNFRSF10B to the plasma membrane in non-small cell lung cancer cells, thus promote chemotherapeutic-mediated apoptosis [38]. The reports of these genes in LGG are not complete, and further analysis is needed to explore their effects.
Our study has proved the link between NECSig and the tumor immune microenvironment. It is important to note that the intricate interaction between tumor cells and TME both helps tumor formation and has a substantial in uence on immunotherapy, thus affecting OS [39]. The functional enrichment analysis was performed, showing that necroptosis-related genes are mainly involved in immune pathways. Signi cant differences between high-risk and low-risk groups were found in the immune in ltration analysis, including B-cell, T-cells, CD4, and NK cells, which con rmed the role of NECsig in regulating tumor immune in ltration. Tumor immune escape is a phenomena in which tumour cells can avoid detection while attacking the immune system through a number of methods, allowing the tumour cells to live and multiply in the body [40]. TIDE is a creative calculation method used to determine the factors of tumor immune escape [41]. Our results show that low-risk populations may have a more positive response to immunotherapy. In summary, these ndings may bring new insights into immunotherapy for the treatment of LGGs.

Conclusion
In short, we constructed the characteristics of NECSig (including BID, H2AFY2, MAPK9, and TNFRSF10B), and systematically evaluated the correlation between the model, prognosis, and clinicopathological characteristics of LGG patients. Survival analysis and the ROC curve indicate that the feature has good predictive performance. Multivariate Cox and strati ed analysis indicate that the four-NECSig feature is an independent prognostic factor of OS. In addition, we built a nomogram that combines the characteristics of NECSig and clinical factors to predict the survival rate of LGG patients. A variety of factors were compared between the high-risk group and the low-risk group, including immune status and response to immunotherapy. This risk pro le based on four NECSig may be expected to be utilized for clinical prediction of the prognosis and immunotherapy response of LGG patients.

Identi cation of Necroptosis-related Differentially Expressed Genes (NECDEs)
To identify NECDEs, the RNA sequencing (RNA-seq) data and the corresponding clinical features from The Cancer Genome Atlas (TCGA) database (https://portal.gdc.cancer.gov/) was obtained. After data screening and ltering, 529 LGG samples and ve normal samples from the TCGA-LGG cohort that had a follow-up time of 30 days or greater were selected. Necroptosis-related genes were collected from the Kyoto Encyclopedia of Genes and Genomes (KEGG) website (https://www.kegg.jp/). The data were download in June 2021. The "limma" package was used to nd NECDEs with p-values less than 0.05.

Function Enrichment Analysis of NECDEs
Among all the available mRNA expression data for LGG in TCGA, a total of 161 necroptosis-related genes were eventually identi ed (Table S1). The R packages ggplot2, enrich plot, and cluster pro ler packages were used to perform Gene Ontology (GO) enrichment analysis and KEGG pathway analysis.

Development and Validation of the NECSig
We used univariate Cox regression analysis to nd predictive genes. Lasso penalized regression, and multivariate Cox regression. The level of risk was calculated as . Based on their median risk score, patients were divided into high-risk and low-risk categories. To generate Kaplan-Meier survival curves, the R packages "survminer" and "survival" were used.

Immune Correlation Analysis
Based on NECSig, we utilized TIMER [21], CIBERSORT [22], quanTIseq [23], MCP-counter [24], xCell [25], and EPIC [26] algorithms are used to analyse the cell types of immune responses in heterogeneous samples to differentiate between high-risk and low-risk groups. The heatmap was used to visualize the differences in the abundance of tumor in ltrating immune cells under distribution algorithms. Immune checkpoints were collected from earlier research. Based on the NECSig, The algorithm was used to calculate the relative concentration of tumor-in ltrating immune cells between the two risk groups based on the NECSig. Finally, the "gsva" programme used ssGSEA to compare immune cells and pathways between the two groups. We analyzed the six types of immune in ltrations using the TCGA dataset and correlated them with the expression levels of four NECSig. Immunotherapeutic response was predicted by TIDE algorithms (http://tide.dfci.harvard.edu/)

The Prediction of Clinical Treatment Response
The IC50 value determines the sensitivity of certain anticancer medicines. The IC50 values of Bortezomib, Cyclopamine, Dasatinib, Rapamycin, and roscovitine were compared between groups in R with pRRophetic and ggplot2 using the Wilcoxon signed-rank test.
Statistical Analysis R was used for all statistical analysis (version 4.1.3). To explore the link between the prognostic model risk score and stemness, stromal score, immunological score, and drug sensitivity, the R programme "limma" was employed. To investigate differences across groups, the Kruskal-Wallis test was utilized. The nomogram was created with the R software "rms". Clinical data was analysed using the chi-squared test. The Student's t-test was performed to determine the statistical signi cance of the difference between the two groups, with p <0.05 indicating that the difference was statistically signi cant.

Declarations
Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and materials All data generated or analyzed during this study are included in this published article.       Immune in ltration and immune checkpoint analysis (A) Immune response heatmaps based on TIMER, CIBERSORT, quanTIseq, MCP-counter, xCell, and EPIC algorithms for high-risk and low-risk groups. (B) In the TCGA-LGG cohort, immune checkpoint expression differed between the low-risk and high-risk groups. *, p < 0.05; **, p < 0.01; ***, p < 0.001.  Comparison of responses to medical treatment between the high-risk and low-risk groups The high-risk scores were related to lower IC50 values for chemotherapeutics and target-directed treatments, such as (A) Bortezomib, (B) Cyclopamine, (C) Dasatinib, (D) Rapamycin, and (E) Roscovitine. (F) TIDE score in low-risk and high-risk groups in the cohort from TCGA-LGG. ***, p < 0.001.

Supplementary Files
This is a list of supplementary les associated with this preprint. Click to download.