Immune Characteristics of eQTL and Gene Risk Model and the Inhibitory Effect of DCTD and RRAS on Ferroptosis in Glioblastoma

Abstract

Background:

Glioblastoma (GBM) is an extremely aggressive brain tumor, marked by restricted therapeutic possibilities and a generally unfavorable prognosis. GBM’s complexity and heterogeneity necessitate comprehensive genetic and immunological profiling to enhance therapeutic strategies.

Methods:

The study integrated The Cancer Genome Atlas (TCGA) and Integrative Epidemiology Unit Open Genome-Wide Association Studies (IEU OpenGWAS) data to identify genetic factors influencing GBM using expression quantitative trait loci (eQTL) and genome-wide association studies (GWAS). Mendelian randomization (MR) analysis revealed 250 GBM-associated genes. A GBM risk prediction model was built using Least Absolute Shrinkage and Selection Operator (LASSO) and Cox regression. The research examined immune infiltration, drug response, and mutation profiles to characterize GBM molecular features. Functional enrichment and in vitro experiments validated key findings.

Results:

The analysis uncovered significant genetic associations with GBM, emphasizing key genes such as follistatin-like 1 (FSTL1), FXYD domain-containing ion transport regulator 5 (FXYD5), Ras-related protein (RRAS), and ring finger protein 216 pseudogene 1 (RNF216P1). The risk model effectively categorized patients into low-risk and high-risk groups, showing significantly worse survival outcomes in the high-risk group. Immune profiling revealed differential infiltration of cancer-associated fibroblasts (CAFs), macrophages, and T cells, which correlated with the expression levels of the genes that were identified. Patients at high risk showed increased sensitivity to chemotherapeutic drugs such as dasatinib and lapatinib, while those at low risk were more responsive to elesclomol and lisitinib. Notably, key genes such as DCMP Deaminase (DCTD) and RRAS were found to regulate ferroptosis, underscoring their potential as therapeutic targets for GBM treatment.

Conclusion:

This study deepens the understanding of GBM by pinpointing critical genetic markers and elucidating their influence on the tumor immune microenvironment (TME) as well as treatment response. The risk model developed in this study holds promise for enhancing prognostic accuracy and facilitating the personalization of GBM therapy.

Graphical abstract

Keywords

GBM / TME / DCTD / RRAS / ferroptosis / prognostic

1. Introduction

Glioblastoma (GBM) ranks among the most malignant and deadly types of primary brain tumors [1, 2, 3, 4]. Although treatments like surgical removal, chemotherapy, and radiation therapy are currently available, their therapeutic efficacy remains limited [5]. The high level of heterogeneity and the complex molecular mechanisms associated with GBM present substantial challenges to research and treatment [6]. In recent years, driven by advancements in bioinformatics approaches and sequencing technologies, researchers have initiated investigations into the GBM utilizing expression quantitative trait loci (eQTL) analysis and genome-wide association study (GWAS) [7, 8, 9]. These studies offer novel insights into the identification of key genes associated with GBM and elucidate their functions in tumor initiation, progression, and therapeutic intervention.

In recent years, the understanding of tumor immune dynamics has continued to advance [10]. GBM is distinguished by density of effective tumor-infiltrating CD8+ T lymphocytes within the tumor immune microenvironment (TME), alongside a concomitant rise in immunosuppressive cell populations [11, 12, 13, 14, 15]. These findings underscore the “cold tumor” characteristics of GBM and highlight that immune activation constitutes a promising strategy for GBM treatment. Specifically, targeting distinct immune cell populations to counteract the immunosuppressive TME is increasingly acknowledged as a critical approach to enhancing GBM therapeutic outcomes [16, 17].

This study integrates data from The Cancer Genome Atlas (TCGA) and Integrative Epidemiology Unit Open Genome-Wide Association Studies (IEU OpenGWAS) project to explore the genetic factors influencing GBM, focusing on the identification of eQTLs and genetic variants associated with risk of GBM. Through Mendelian randomization (MR) analysis, we identify 250 genes significantly associated with GBM, many of which have potential implications for treatment. This work deepens our understanding of GBM biology and lays the groundwork for developing targeted genetic and immune-based therapies. Through this integrative approach, we aim to improve prognostic accuracy and identify novel therapeutic targets, particularly those involved in immune modulation and ferroptosis regulation. Our findings could significantly impact the development of personalized treatment regimens for GBM, contributing to more effective clinical management.

2. Materials and Methods

2.1 Data Set

The dataset comprises a total of 176 GBM samples, categorized into 5 non-cancerous adjacent tissues and 171 tumor tissues, along with their respective clinical details from TCGA. Additionally, we obtained 59 GBM samples with matched clinical data from the Gene Expression Omnibus (GEO) repository (accession number GSE4421, https://www.ncbi.nlm.nih.gov/geo/), which were also subjected to normalization procedures. eQTL data encompassing 19,942 genes (refer to Supplementary Table 1), where significant single nucleotide polymorphisms (SNPs) were identified based on a significance threshold of p 5 \times 10^-6. Furthermore, GWAS data for GBM from the FinnGen R10 release (https://www.finngen.fi/fi) were utilized, consisting of 253 GBM cases and 314,193 control subjects.

2.2 Mendelian Randomization Analysis

To guarantee the reliability of the MR analysis, potential instrumental variables (IVs) were chosen using a significance criterion of p 5 \times 10^-6. This threshold was chosen due to the relatively small sample size of eQTL data (less than 50,000 samples per eQTL), which warranted a more stringent filtering criterion of p < 5 \times 10^-6 to identify strong genetic instruments. To reduce the impact of linkage disequilibrium (LD) and ensure the independence of IVs, using parameters of R² = 0.001 and kb = 10,000. Additionally, palindromic SNPs, which could cause ambiguity in allele orientation, were excluded to avoid confounding results. Heterogeneous eQTLs were excluded. Pleiotropy in eQTLs was tested using the “mr_pleiotropy_test” function, and pleiotropic eQTLs were removed. Finally, SNP outliers were detected using the “run_mr_presso” function, and biased SNPs were excluded.

In the MR analysis, the core assumptions were carefully considered: the relevance of the IVs to the exposure, their independence from confounding factors, and the exclusion of direct effects on the outcome. The Wald ratio method was employed to obtain unconfounded MR estimates for each SNP. To further assess the robustness of the findings and account for potential pleiotropy, meta-analysis of the individual estimates was conducted using multiple methods. These methods allow for the assessment of causality while minimizing bias from pleiotropic effects and ensuring the reliability of the MR estimates.

2.3 Functional Enrichment Analysis

Functional enrichment of differentially expressed gene modules between TCGA-derived high- and low-risk groups was performed using Gene Ontology (GO, https://geneontology.org/) and Kyoto Encyclopedia of Genes and Genomes (KEGG, https://www.kegg.jp/kegg/) analyses via the clusterProfiler (v4.10.0, https://www.bioconductor.org/packages/release/bioc/html/clusterProfiler.html) package in R. GO terms were identified with the enrichGO function (q 0.05), and KEGG pathways with the enrichKEGG function.

2.4 Risk Model Construction

Expression data of intersecting genes were extracted from TCGA and combined with clinical data. The coxph function to select candidate genes significantly affecting GBM patient prognosis. LASSO regression analysis was then performed using the “glmnet” and “survminer” packages (v4.1-8, v0.4.9) in R language (v4.3.2) with the glmnet and cv.glmnet functions. Multivariate Cox regression analysis determined the final prognostic genes (p 0.05), and riskScore were calculated using the coef function (riskScore = \mathrm{\Sigma } (coefficient_i \times expression of signature gene_i)) (Supplementary Table 2). GBM samples were generated with the survival, regplot, and rms packages (v6.8-2, v1.1, v6.8-2) in combination with the coxph function. Differential expression of candidate genes was assessed.

2.5 Single Sample GSEA Analyses

All configured with the same parameter settings from a unified configuration file. The ssGSEA scores for each pathway were calculated using the gsva function and subsequently normalized via the normalize function.

2.6 Tumor Immune Dysfunction and Exclusion (TIDE) Prediction

GBM samples from TCGA with higher scores indicating stronger immune escape capability and poorer immunotherapy efficacy. Since the TIDE database (http://tide.dfci.harvard.edu/) does not include a specific category for GBM, the tumor type was designated as “Other” and the prior immunotherapy type was specified as “None”. Prior to uploading to the TIDE database, the TCGA GBM expression data were log2-transformed. The chi-square test was utilized to assess differences in immunotherapy response rates between these two groups.

2.7 Drug Sensitivity Analysis

The “pRRophetic” R package (v0.5) was employed to forecast sensitivity to chemotherapy drugs, and a correlation analysis was carried out between the drug sensitivity data and the riskScore, applying a filter with a significance threshold of p 0.001.

2.8 Reagents, Antibodies, siRNA, and shRNA

RSL3, ferrostatin-1, deferoxamine, JKE-1674, and atezolizumab (HY-100218C, HY-100579, HY-B1625R, HY-138153, HY-P9904) were purchased from MedChemExpress (Shanghai, China). The siRNAs and shRNA plasmids for DCTD and RRAS were constructed by Shanghai Genechem (China). Antibodies against RRAS (27457-1-AP), SLC7A11 (26864-1-AP), DCTD (68357-1-Ig), ACSL4 (22401-1-AP), GAPDH (10494-1-AP), and GPX4 (30388-1-AP) were obtained from Proteintech (Wuhan, China). All the primary antibodies were diluted at a ratio of 1:1000. Specifically, the sequence for both si-RRAS1 and sh-RRAS1 is 5^′-ACACGAAGATCTGCAGTGTGGAT-3^′, while the sequence for si-RRAS2 and sh-RRAS2 is 5^′-AGACGAAGATCTGCAGTGTGGAC-3^′. Additionally, the sequences for si-DCTD1 and sh-DCTD1 are 5^′-GGAUCUUCUGCUUCCAAUA-3^′, and those for si-DCTD2 and sh-DCTD2 are 5^′-CCAAGACAGAUCCUGUUAU-3^′.

2.9 Cell Culture and Transfection

Human glioblastoma cell lines U-87MG, A172, and TG-905 (ATCC, Shanghai, China) were cultured in a CO₂ incubator (TCI-80, Panasonic Healthcare Holdings Co., Ltd, Ora, Japan) at 37 °C and monitored every 8 or 24 h as required. Cells were subcultured at 80–95% confluence. All disposable materials were UV-sterilized in a laminar flow cabinet (Jiesen, Xuzhou, Jiangsu, China), and procedures were performed with sterile gloves. After trypsin digestion, cells were suspended in complete growth medium, centrifuged, and resuspended in fresh medium for continued culture. During cell culture, we regularly add mycoplasma scavengers to the medium to avoid mycoplasma contamination.

Cells were cultured in 6-well plates to reach 70–80% confluence. The medium was replaced with serum-free opti-MEM before transfection. Transfection was performed using Lipofectamine 3000 (L3000015, Thermo Fisher Scientific, Shanghai, China) according to the manufacturer’s instructions, with either plasmid DNA or siRNA. This complex was then added dropwise to the cells in the 6-well plates. For siRNA transfection, a similar procedure was followed, but with siRNA molecules instead of plasmids. The cells were exposed to the transfection complexes for a 6-hour period at 37 °C. Following this incubation, the existing medium was substituted with new complete medium supplemented with 10% fetal bovine serum (FBS), aiming to promote cell recovery. Post-transfection, cells were incubated for 48 hours under standard culture conditions before being used for further experimentation, such as RNA/protein extraction, phenotypic analysis, or functional assays. All cell lines were validated by short tandem repeat (STR) profiling and tested negative for mycoplasma.

2.10 Western Blotting

Protein extracts were prepared by lysing cells in cold lysis buffer with protease and phosphatase inhibitors, followed by centrifugation to obtain the soluble protein fraction. Protein concentrations were determined using a BCA assay (E-BC-K318-M, Elabscience Biotechnology Inc., Wuhan, Hubei, China), and samples were mixed with sodium dodecyl sulfate (SDS) buffer and denatured by heating. Proteins were separated by sodium dodecyl sulfate - polyacrylamide gel electrophoresis (SDS-PAGE) using a 10% resolving gel.

2.11 MTT, Colony Formation, Wound Healing, and Invasion Assays

2.11.1 MTT Assay

Cells (5000–10,000 per well) were seeded in 96‑well plates. On day 1, the medium was replaced with 100 µL 3-(4,5-Dimethylthiazol-2-yl)-2,5-diphenyltetrazolium bromide (MTT) solution (5 mg/mL in PBS) and incubated for 2–4 h at 37 °C. The solution was aspirated, and 100 µL dimethyl sulfoxide (DMSO) was added to dissolve formazan crystals. Absorbance was measured at 490 nm using a microplate reader. This procedure was repeated daily for 4 days to monitor proliferation. Optical density (OD) values were plotted with GraphPad Prism to determine cell viability.

2.11.2 Colony Formation Assay

Cells (500/well) were seeded in 6‑well plates and cultured for 1–2 weeks, with medium changes every 3–4 days. Colonies (>50 cells) were fixed with methanol, stained with crystal violet, and counted under a light microscope.

2.11.3 Wound Healing Assay

A scratch was created in confluent monolayers using a sterile pipette tip. Cells were rinsed with PBS to remove debris and cultured in serum‑free medium. Wound areas were imaged at 0, 24, and 48 h to assess closure.

2.11.4 Transwell Invasion Assay

Matrigel‑coated Transwell inserts were seeded with 2 \times 10⁵ cells in serum‑free medium. After 2–3 days, invading cells were fixed.

2.12 Cell Viability and Lipid Peroxidation Analysis

Cell viability was evaluated using the CellTiter-Glo 3D assay by measuring luminescence after treatment with ferroptosis inducers and inhibitors. Lipid peroxidation was assessed using the BODIPY C11 probe, with the ratio of oxidized to reduced fluorescence indicating peroxidation levels. Glioblastoma cells were treated for 24–72 hours. Luminescence and fluorescence were measured using microplate readers. Higher luminescence indicates higher viability, while a higher oxidized/reduced fluorescence ratio indicates greater lipid peroxidation.

2.13 U-87MG Xenografts and In Vivo Treatment

Regarding the animal housing environment, experimental animals are housed in a dedicated animal facility with stringent environmental controls. This facility complies with the CV grade standard. The indoor temperature is consistently maintained at 22 \pm 2 °C, a range that has been extensively validated. This temperature range effectively minimizes physiological variations caused by fluctuations, thereby ensuring the stability and reliability of experimental results. Humidity is controlled at 50% \pm 5%, which helps prevent respiratory diseases and ensures animal health. The lighting environment was kept on a 12-hour light/dark cycle to simulate natural circadian rhythms and maintain the animals’ biological regularity. In the xenograft experiments, 1 \times 10⁶ U-87MG glioblastoma cells, transfected with siRNA or control siRNA (sicon), were mixed with 100 µL Matrigel (1:1) and implanted subcutaneously into the right flank of 5–6-week-old male nude mice. Tumor growth was monitored every 2–3 days by caliper measurements and volume calculated using V = L \times W² \times 0.52. Mice were assigned to treatment groups once tumors reached 80–100 mm³. They were administered JKE-1674 (25 mg/kg) orally every other day, and atezolizumab (10 mg/kg) intravenously twice weekly, for a duration of 3 weeks. Throughout the treatment period, peripheral blood was collected weekly for monitoring liver function, specifically alanine transaminase (ALT) and aspartate aminotransferase (AST) levels, to assess potential drug toxicity. Following the completion of the 3-week treatment phase, the mice were euthanized via cervical dislocation, and their tumor tissues were collected for subsequent analysis, including histological evaluation, protein analysis, and RNA extraction. Our in vivo studies strictly follow the “ARRIVE and PREPARE” guidelines, which promise to improve the scientific nature and transparency of the study while ensuring ethical treatment of the animals.

2.14 Statistical Analysis

For normally distributed variables, an independent Student’s t-test was used to assess statistical significance, whereas the Mann-Whitney U test (also known as the Wilcoxon rank-sum test) was used for non-normally distributed data.

3. Results

3.1 Identification of GBM-Related Genes Through eQTL MR Analysis

We analyzed eQTL data from blood samples for 19,942 genes to explore associations with GBM. After filtering, 15,306 eQTLs and 52,674 SNPs were selected for MR analysis. Cross-referencing with GBM GWAS data identified 49,823 SNPs with F-statistics >10. A total of 15,035 genes were included in MR analysis. Using inverse-variance weighted (IVW) method, 250 genes showed significant associations with GBM risk (147 increased risk, 103 reduced risk) (Table 1, Supplementary Tables 3–7).

3.2 Differential Analysis and Identification of Candidate GBM Genes

To further screen for GBM-related genes, a differential analysis was conducted on mRNA data derived from TCGA cancerous tissues and adjacent normal tissues. Using a false discovery rate (FDR) threshold of less than 0.05 and a logFC value greater than 0.585 as screening criteria, a total of 9942 differentially expressed genes were identified (Supplementary Table 8). Heatmap and volcano plot analyses revealed significant differences in gene expression between the two groups, with fewer upregulated genes observed in the tumor group compared to the adjacent normal group (Fig. 1A,B). By intersecting the risk gene sets obtained from both TCGA and MR analyses, we identified 50 genes that were significantly associated with GBM. Among these, 39 genes were classified as high-risk genes for GBM, while 11 genes were identified as protective genes for GBM (Fig. 1C,D, Supplementary Table 9).

3.3 LASSO Regression Risk Model

Univariate Cox regression analysis conducted on the 50 candidate genes indicated that 11 genes exhibited a statistically significant link to the prognosis of GBM patients (Fig. 2A). Subsequent MR analysis indicated that the eQTLs of these 11 genes may have a significant impact on GBM prognosis (Fig. 2B). The Circos plot visually illustrates the chromosomal locations of these 11 genes (Fig. 2C). Cross-validation curve analysis and LASSO coefficient path analysis determined four risk genes (FSTL1, FXYD5, RRAS, and RNF216P1) as variables in the model and provided the model coefficients for each gene (Fig. 2D,E, Supplementary Table 10).

3.4 Validation of the Risk Model

TCGA tumor samples were divided randomly into training and testing cohorts. Four key genes were significantly higher in the high-risk group than in the low-risk group. High-risk individuals showed higher mortality and shorter survival times (Fig. 3A–C). Survival analysis confirmed that high-risk patients had notably lower survival rates (Fig. 3D–F). Data validation analysis was also performed using the GSE4421 dataset. Unfortunately, the RNA expression data in this dataset does not include RNF216P1, so the results from this analysis cannot be used (Supplementary Fig. 1). These findings were consistent when validated across all TCGA tumor samples, underscoring the importance of these four genes in GBM prognosis.

3.5 Prognostic Prediction Using the Risk Model

In conjunction with age and gender, was significantly associated with GBM prognosis, exhibiting a hazard ratio (HR) greater than 1.5. This suggests that an elevated riskScore is correlated with a poorer prognosis (Fig. 4A,B). Receiver operating characteristic (ROC) curves and the concordance index (C-index) further demonstrated that the riskScore possessed high specificity and sensitivity for predicting GBM prognosis (Fig. 4C,D, Supplementary Fig. 2). A more sophisticated predictive model integrating age, gender, and the riskScore was subsequently developed to enhance prediction accuracy (Fig. 4E). This model exhibited robust accuracy in forecasting 2- and 3-year survival rates, albeit with slightly diminished performance for 1-year survival prediction (Fig. 4F). Principal component analysis (PCA) effectively differentiated between high-risk and low-risk groups, particularly through the application of the four identified risk-associated genes (Fig. 4G–I).

3.6 Differential Signaling Pathways Between High and Low-Risk Groups

These genes were categorized into two groups: those more highly expressed in the high-risk group and those more prevalent in the low-risk group. GO and KEGG pathway enrichment analyses were then conducted for each gene group. Results revealed in the high-risk group, that gene sets, including those linked to cytokine activity, chemokine activity, and immune cell migration were significantly enriched. Conversely, in the low-risk group, gene sets associated with ion transport and neurotransmitter regulation, such as modulation of chemical synaptic transmission and regulation of metal ion transport, were mainly enriched (Fig. 5A–D, Supplementary Tables 11–13). Consistent with these findings, KEGG enrichment analysis mirrored the GO results, demonstrating high-risk group was related to cytokines, chemokines, and inflammation, such as interleukin-17 (IL-17) and tumor necrosis factor (TNF) signaling (Fig. 5E,F, Supplementary Table 14). The low-risk group was significantly associated with neuroactive ligand-receptor interactions and calcium signaling pathways, which are involved in ion transport and neurotransmitter regulation. The low-risk group exhibited increased expression of genes related to ion transport and neurotransmitter regulation (Fig. 5G,H, Supplementary Table 15). These findings indicate that the riskScore effectively distinguishes major gene expression patterns in GBM patients.

3.7 Gene Mutations and TMB in High and Low-Risk Groups

Next, we conducted a comprehensive examination of genetic mutation characteristics in both high- and low-risk groups based on TCGA data. The findings revealed notable variations in the mutation rates of the wild-type human p53 (TP53), epidermal growth factor receptor (EGFR), and phosphatase and tensin homolog (PTEN) across the two groups. In particular, in the high-risk group, mutations in PTEN were more prevalent, but mutation frequencies of TP53 and EGFR were more prevalent in low-risk group (Fig. 6A,B). Co-mutation analysis further revealed notable patterns of co-occurrence and mutual exclusivity. For instance, PTEN and TP53 co-occurred more frequently in the high-risk group, while TP53 co-occurred with isocitrate dehydrogenase 1 (IDH1) and alpha-thalassemia mental retardation X-linked (ATRX) in both groups. Conversely, TP53 and EGFR exhibited mutual exclusivity. Notably, the high-risk group exhibited a larger number of concurrently occurring gene pairs (16 pairs) in comparison to the low-risk group (9 pairs) (Fig. 6C,D). variant allele frequency (VAF) analysis revealed that the high-risk group had generally lower average VAF values than the low-risk group. Among the high-risk group, PTEN, TP53, and erythrocytic 1 (SPTA1) were the top genes based on VAF. In contrast, the low-risk group showed ATRX, retinoblastoma protein 1 (RB1), and TP53 as the most prominent genes. Importantly, the average VAF of TP53 differed significantly between the two groups (Fig. 6E,F).

We further explored the impact of tumor mutational burden (TMB) on GBM and found that patients with higher TMB levels had notably better survival outcomes than those with lower TMB levels (Fig. 6H). Nevertheless, no significant difference in TMB was observed (Fig. 6G). The patients were divided into four groups based on TMB levels and riskScores. The low TMB but high-risk group had the worst survival prognosis. The high TMB and low-risk group showed the better quality of life (Fig. 6I). Collectively, these findings indicate that the riskScore effectively distinguishes gene mutation patterns and characteristics among the different groups.

3.8 Tumor Immune Microenvironment in High and Low-Risk Groups

We utilized the R programming language to assess immune infiltration levels in the samples through multiple immune cell estimation approaches, such as Estimate, ssGSEA, MCPcounter, Timer, Epic, CIBERSORT, and Quantiseq. Out of these seven techniques, the high-risk group demonstrated markedly increased infiltration levels of T cells (specifically CD4+ and CD8+ T cells) when compared to the low-risk group. Furthermore, in nearly all methods except CIBERSORT, the high-risk group displayed significantly greater infiltration of monocytes and tumor-associated fibroblasts. (Fig. 7A–G, Supplementary Table 16). The high-risk group also exhibited higher stromal and immune scores overall. Additionally, scores associated with immune inflammatory responses—such as those involving major histocompatibility complex (MHC), interferon (IFN), antigen-presenting cell (APC), and T cell-related gene sets—were considerably elevated in this group (Fig. 7D,H, Supplementary Table 16). Using the ssGSEA method, the infiltration levels of myeloid-derived suppressor cells (MDSCs) were found to be significantly higher in the high-risk group compared to the low-risk group (Fig. 7G, Supplementary Table 16). Examination of immune checkpoint expression revealed that several key checkpoints, such as cytotoxic T lymphocyte-associated antigen-4 (CTLA4), Programmed cell death protein 1‌ (PDCD1), Programmed Cell Death 1 Ligand 1 (CD274), and Recombinant Hepatitis A Virus Cellular Receptor 2 (HAVCR2), were significantly overexpressed in the high-risk group (Fig. 7I,J). Collectively, these results highlight a clear distinction in the tumor immune microenvironment between the high-risk and low-risk groups.

3.9 Prediction of Immunotherapy Response in Different Groups

Study used the TIDE platform to judge immunotherapy response immune cell infiltration in GBM samples. Higher infiltration of CAFs and increased expression of IFN, CD8, and programmed cell death ligand 1 (PD-L1) in high-risk group. The high-risk group also had higher immune dysfunction scores, indicating more significant infiltration of tumor-associated immune cells. Unlike previous findings, the high-risk group had lower infiltration of MDSCs and M2 macrophages. However, there was no significant difference in predicted immunotherapy response between the two groups (Supplementary Fig. 3, Supplementary Table 17).

3.10 Prediction of Chemotherapy Drug Sensitivity in High- and Low-Risk Groups

The study used the R package “pRRophetic” to estimate TCGA GBM patients’ sensitivity to various chemotherapeutic agents. A significant inversion was found in dasatinib, lapatinib, and temsirolimus, suggesting higher riskScores correlate with better responses to these drugs. Conversely, a positive correlation was observed between riskScore and the IC₅₀ values of elacridar, lestaurtinib, and tubastatin A, indicating lower riskScores may lead to greater sensitivity to these drugs (Supplementary Fig. 4). A systematic analysis of the link between riskScore and chemotherapy drug sensitivity allowed us to identify differing drug sensitivity among patient groups with varying riskScores.

3.11 Relationship Between RiskScore Genes and Tumor Immune Microenvironment

The boxplots of gene expression demonstrated that the mRNA expression levels of the four risk-associated genes were significantly elevated in tumor tissues relative to adjacent normal tissues (Fig. 8A). Subsequently, we investigated the correlation between the abundance of tumor-infiltrating immune cells and the expression of these four risk genes to elucidate their roles in immune cell infiltration. Across six distinct immune cell infiltration scoring methodologies, FXYD5 exhibited a significant positive correlation with macrophage infiltration levels. Furthermore, FXYD5 demonstrated a positive relationship CD4+ and CD8+ T cells across the majority of the scoring methods employed. In the Epic and MCPcounter models, FXYD5 was found to be positively linked with CAFs. Additionally, using the ssGSEA, MCPcounter, and CIBERSORT approaches, a positive association was observed between FXYD5 and B cell infiltration levels (Fig. 8B–H). In a similar manner, RRAS expression levels were related to the presence of CAFs, T cells, and macrophages, within the microenvironment. However, RNF216P1 expression showed a bad relationship with multiple tumor-infiltrating immune cells, such as macrophages, CAFs, and T cells, as indicated by six different scoring systems (Fig. 8B–H). These results suggest that riskScore serves as a significant indicator for evaluating immune cell infiltration in GBM patients.

3.12 DCTD and RRAS Promote GBM Cell Proliferation, Migration, and Invasion

The GBM samples were categorized into two groups—low expression and high expression—according to the median expression values of the four risk-associated genes, after which survival analysis was conducted. The findings revealed that higher expression levels of DCTD and RRAS were linked to worse prognostic outcomes (Fig. 9A,B), whereas the remaining two genes did not exhibit a statistically significant difference between the two expression groups (Fig. 9C,D). Based on these results, we selected DCTD and RRAS for further experiments. We transfected GBM U-87MG cells with siRNAs targeting DCTD and RRAS.

We transfected GBM U-87MG cells with siRNAs targeting DCTD and RRAS. Knockdown of DCTD or RRAS led to a reduction in their respective protein levels (Fig. 10A,B). MTT, colony formation, wound healing, and invasion assays demonstrated that knockdown of DCTD or RRAS reduced U-87MG cell proliferation, migration, and invasion (Fig. 10C–J, Supplementary Fig. 5). Using a CDX model, we found that tumors with DCTD or RRAS knockdown were significantly smaller than those in the control group (Fig. 10K–N), We have also added a new cell line, U251, and obtained similar results (Supplementary Fig. 6) indicating that knockdown of DCTD or RRAS inhibits GBM CDX model tumor growth.

3.13 DCTD and RRAS Regulate Key Ferroptosis Proteins and Affect U-87MG Phenotype

Using GEPIA for correlation analysis, a significant positive correlation was found between RRAS and GPX4, a key ferroptosis gene (Supplementary Fig. 7A). Although the correlation between DCTD and GPX4 was not statistically significant, the p-value was close to 0.05 with a positive correlation coefficient (Supplementary Fig. 7B). Knockdown of DCTD or RRAS led to decreased expression of ferroptosis inhibitory proteins SLC7A11 and GPX4, and increased expression of the ferroptosis-promoting protein ACSL4 (Fig. 11A–D, Supplementary Fig. 5A–D). Colony formation, wound healing, and invasion assays showed that cells with si-DCTD or si-RRAS and exogenous ACSL4 had lower proliferation, migration, and invasion capabilities compared to cells with si-DCTD or si-RRAS alone. Conversely, overexpression of DCTD or RRAS enhanced these abilities, and knocking down ACSL4 in DCTD or RRAS overexpression cells increased proliferation, migration, and invasion (Fig. 11E–P, Supplementary Fig. 5G–I). This indicates that DCTD and RRAS influence the ferroptosis signaling pathway and affect U-87MG cell phenotype. Meanwhile, we investigated the mechanism by which DCTD regulates ferroptosis. Initially, using q-PCR, we demonstrated that DCTD knockout did not induce changes in the transcriptional level of ACSL4 (Supplementary Fig. 8A). Subsequently, we hypothesized whether DCTD modulates ACSL4 expression by influencing ubiquitination-related enzymes. To this end, we predicted E3 ubiquitin ligases regulating ACSL4 via the ubibrowser website and identified autocrine motility factor receptor (AMFR) as a regulator of ACSL4 (Supplementary Fig. 8B). Furthermore, through analysis on gepia2.cancer, we observed a positive correlation between DCTD and AMFR (Supplementary Fig. 8C). Next, we transfected si-RRAS and si-DCTD into the U251 cell line and assessed the protein expression levels of RRAS, DCTD, and AMFR using Western blotting. Our results revealed that RRAS and DCTD regulate the protein level of AMFR, with a positive correlation among them (Supplementary Fig. 8D,E). Finally, we designed an experiment to transfect an AMFR overexpression plasmid following DCTD knockout and examined the ubiquitination level of ACSL4. We found that the ubiquitination level of ACSL4 decreased significantly upon DCTD knockout. However, when the AMFR overexpression plasmid was transfected subsequent to DCTD knockout, the ubiquitination level of ACSL4 was restored (Supplementary Fig. 8F).

3.14 DCTD and RRAS Inhibit Ferroptosis in GBM Cells

To validate the effect of ferroptosis on GBM cells, we treated cells with the ferroptosis activator RSL3 and the ferroptosis inhibitors ferrostatin-1 and deferoxamine (iron chelator). Activation of the ferroptosis pathway significantly reduced cell viability in three GBM cell lines, and this effect was reversed by ferroptosis inhibitors (Fig. 12A–C), indicating that activation of the ferroptosis pathway decreases GBM cell viability. In U-87MG cells with DCTD or RRAS knockdown or overexpression, RSL3 treatment at varying concentrations for 24 hours showed that the IC₅₀ of RSL3 decreased with DCTD or RRAS knockdown and increased with overexpression. Cells treated with RSL3 showed reduced viability in si-DCTD, si-RRAS, OE-DCTD, and OE-RRAS groups compared to controls. Additionally, lipid peroxidation levels were higher in the DCTD or RRAS knockdown groups and lower in the overexpression groups after RSL3 treatment (Fig. 12D–O), suggesting that knockdown of DCTD or RRAS enhances ferroptosis in U-87MG cells treated with RSL3. We detected the intracellular reactive oxygen species (ROS) levels after the knockout and overexpression of RRAS and DCTD using a flow cytometer. We found that the expression level of ROS increased significantly after the knockout of RRAS and DCTD, while it decreased significantly after their overexpression (Supplementary Fig. 9).

3.15 JKE-1647 Combined With PD-L1 Inhibitor Inhibits Subcutaneous Tumor Growth

Given that JKE-1647 is more easily metabolized in vivo and induces ferroptosis [18], we chose JKE-1647 as the ferroptosis inducer for in vivo experiments. Since DCTD and RRAS are positively correlated with immune cell infiltration, we selected the PD-L1 inhibitor atezolizumab to verify the efficacy of immunotherapy for GBM. We established a subcutaneous tumor model with three groups: si-con + JKE-1647, si-con + JKE-1647 + atezolizumab, and si-DCTD/si-RRAS + JKE-1647 + atezolizumab. Consistent with in vitro results, ferroptosis activation significantly inhibited subcutaneous tumor growth, and the combination of JKE-1647 and atezolizumab had a stronger inhibitory effect than JKE-1647 alone. As expected, knockdown of DCTD or RRAS combined with JKE-1647 and atezolizumab more effectively inhibited GBM subcutaneous tumor growth (Fig. 13A–C, F–H), indicating that DCTD or RRAS knockdown promotes ferroptosis activation. Plasma ALT and AST levels showed no significant differences among the groups, suggesting that JKE-1647 combined with atezolizumab does not cause noticeable hepatotoxicity (Fig. 13D,E,I,J). We added the murine glioma cell line GL261 to establish a subcutaneous tumor model, and obtained the same conclusion (Supplementary Fig. 10).

4. Discussion

Numerous studies have demonstrated that eQTL analysis can effectively identify the association between SNP genotypes and gene dysregulation in GBM [7, 8, 9]. A meta-analysis of SNP data from glioma-associated loci identified 12 candidate susceptibility SNPs and investigated 12 glioma-associated genes that had not previously been implicated in GBM [1]. Another study identified rs3761124 as a functional variant located on 20q13.33, which is associated with the risk of glioma and GBM. This variant regulates the expression of stathmin-3 (STMN3) and other potential genes across various cellular contexts [7]. Cancer eQTL analysis can specifically focus on differentially regulated genes in GBM by utilizing fold-change expression levels, thereby highlighting the enrichment of genes associated with neuronal functions. While individual variations may potentially impact the results, the fold-change-based screening of probe sets uncovers gene dysregulation patterns that are characteristic of GBM [2]. A systematic evaluation of the genetic impact on tumor gene expression revealed 13 cis-eQTL target genes associated with GBM survival, underscoring the substantial genetic influence on tumor gene expression [8]. eQTL analysis not only investigates the impact of copy number variations on single-gene expression but also explores large-scale genomic associations across entire gene networks [9]. These results suggest that eQTL-related analyses can reveal the characteristics of gene regulation in GBM. In this study, we identified 250 genes associated with GBM by conducting MR analysis using blood eQTL data and GBM GWAS data.

Next, leveraging the TCGA dataset, we conducted a differential gene expression analysis comparing normal and tumor tissues to identify genes with significantly altered mRNA expression levels. Following this, we applied LASSO regression to construct a prognostic model, which determined the risk coefficients of four key genes: FSTL1, FXYD5, RRAS, and RNF216P1. This allowed us to calculate the riskScore for each sample. Notably, FSTL1, a member of the secreted protein acidic and rich in cysteine (SPARC) family, is up-regulated in glioblastoma, and its co-expression with p53 correlates with poor survival in GBM patients. In our prognostic model, it is associated with the high-risk population and a poor prognosis [19]. FSTL1 plays a pivotal role in temozolomide resistance by modulating discointeracting protein 2 homolog A (DIP2A) protein localization, H3K9 acetylation, and O6-methylguanine methyltransferase (MGMT) transcription [20]. Elevated FSTL1 expression inhibits the nuclear translocation of DIP2A, leading to increased MGMT expression and consequently enhancing glioblastoma resistance to temozolomide [21]. FXYD5, part of the FXYD ion-transport regulator family, has unclear precise roles in glioma, but its significance in our prognostic model points to a strong association with glioma prognosis [22]. At least half of the survey datasets demonstrate a marked increase in the expression levels of multiple ion channels, including FXYD5, in GBM [23]. A study has demonstrated that alterations in the RRAS gene are predominantly observed in gemistocytic astrocytomas and secondary glioblastomas, rather than fibrillary astrocytomas and primary GBM [24]. Furthermore, patients harboring RRAS deletions exhibit significantly reduced survival rates [25]. The existing literature has not documented the relationship between RNF216P1 and GBM. To date, RNF216P1 has only been identified in hepatocellular carcinoma, where it may be involved in an EZH2-related ceRNA pathway [26]. FSTL1, FXYD5, and RRAS have been investigated to a certain extent in GBM; however, research specifically addressing their roles is limited [27, 28]. We need a more detailed exploration of their associations with the GBM immune microenvironment.

In GBM, Tumor-associated macrophages (TAMs) consist of macrophages, making up the largest proportion of immune cells in the tumor [29]. Their function in GBM exemplifies a prevalent phenomenon observed in the TME [30, 31]. Unlike primary brain tumors, metastatic brain tumors lack microglia at their core and are instead populated by macrophages [18, 32, 33, 34], with recurrent GBM being predominantly macrophage-driven [35]. TAM-derived Oncostatin M (OSM), its receptor OSMR, leukemia inhibitory factor (LIF), and its receptor LIFR activate the STAT3 signaling pathway. This activation can induce a mesenchymal (MES)-like transformation in GBM cells both in vivo and in vitro [36, 37]. M2 macrophages sustain their immunosuppressive phenotype via integrin \alphav\beta5 (ITG\alphav\beta5), with osteopontin (OPN) serving as its principal ligand. The absence of OPN reduces M2 macrophage infiltration, enhances CD8+ T cell activity against GBM cells, and prolongs survival. ITG\alphav\beta3 promotes M2 macrophage polarization and abnormal blood vessel formation via the Src-PI3K-YAP pathway. TAMs play a key role in GBM progression. Our risk model effectively predicts TAM infiltration, especially M2 subtype levels [38, 39]. HAVCR2 promotes tumor-supportive macrophage differentiation through IL-6 signaling and correlates with increased macrophage infiltration in high-risk GBM [40].

Compared to TAMs, which constitute 50% of all viable GBM cells, MDSCs exert a greater suppressive effect on the TME than TAMs and Tregs [41, 42]. MDSCs are classified into two main categories, one of which is monocytic myeloid-derived suppressor cells (M-MDSCs), which resemble monocytes, which resemble neutrophils. M-MDSCs primarily tend to polarize into TAMs and rapidly exert immunosuppressive effects, while polymorphonuclear myeloid-derived suppressor cells (PMN-MDSCs) mainly induce long-term immune tolerance [43]. The chemokines, including C-X-C chemokine receptor type 4/C-X-C motif chemokine 12 (CXCR4-CXCL12), CXCR2-CXCL5/8, and CCR2-CCL2, play important roles in M-MDSC mobilization, with CCL2 and CCL15 being especially relevant in colorectal cancer [44, 45]. PMN-MDSC recruitment is mainly mediated by CXCR1-CXCL8, CXCR2-CXCL8, CCR5-CCL5, CXCL6, and CXCL12, along with factors such as CCL2, CCL3, and hypoxia [46, 47]. Genes highly expressed in the high-risk group are significantly enriched in chemokine-related pathways compared to those in the low-risk group [48, 49, 50, 51]. The elevated levels of MDSCs in the high-risk group indicate that riskScore is a key marker for distinguishing MDSC infiltration in GBM.

CD8+ T cells lose functionality in the immunosuppressive environment of glioblastoma. This decline is caused by gene expression changes, epigenetic modifications, and immunosuppressive cells like Tregs, TAMs, and MDSCs. Despite higher CD8+ T cell infiltration, high-risk patients still have poorer clinical outcomes [48, 49]. This may be attributed to the high levels of infiltrating Tregs, TAMs, and MDSCs, as well as the inhibitory effects of co-inhibitory immune checkpoints on CD8+ T cells [50].

To better characterize this risk model, we analyzed the link between riskScore and tumor mutations in depth. TP53, EGFR, and PTEN showed significant differences in mutation frequencies between high-risk and low-risk groups. Notably, the high-risk group had more co-mutations, indicating a more complex mutational profile. While TMB did not differ significantly between the groups, combining riskScore with TMB in subgroup analyses effectively distinguished patient survival outcomes, highlighting riskScore’s potential prognostic value. Additionally, examining the relationships between the four key risk genes and the TME revealed that FXYD5 and RRAS correlated positively with T cells, macrophages, and CAFs, whereas RNF216P1 expression correlated negatively with the infiltration of these immune cells.

Based on the pRRophetic “package predictions”, we confirmed that High-risk patients had lower IC₅₀ to Dasatinib, Lapatinib, and Temsirolimus, whereas low-risk patients showed lower IC₅₀ to Elesclomol, Lisitinib, and Tubastatin A. This systematic analysis helps clarify differences in chemotherapy drug sensitivity among GBM patients based on riskScore. Given that clinical trials are currently underway evaluating Lapatinib in combination with radiotherapy and temozolomide (RT/TMZ) for GBM treatment [51], riskScore may potentially be valuable in predicting the effectiveness of Lapatinib for GBM.

Our study elucidates the pivotal roles of DCTD and RRAS in the progression of GBM. By employing siRNA-mediated knockdown of DCTD and RRAS, we observed a marked reduction in the proliferation, migration, and invasion capabilities of U-87MG and U251 cells, as well as a suppression of tumor growth in cell-derived xenograft (CDX) models. The result showed DCTD and RRAS can accelerate GBM progression by promoting these biological processes. Further investigation revealed that DCTD and RRAS regulate the expression of key proteins associated with ferroptosis. Although the correlation between DCTD and the ferroptosis-related gene GPX4 did not achieve statistical significance, it approached significance, indicating a potential association. Following the knockdown of DCTD or RRAS, we detected a decrease in the expression of the ferroptosis-inhibitory protein SLC7A11 and its downstream effector GPX4, alongside an increase in the expression of the ferroptosis-promoting protein ACSL4. Collectively, these results imply that DCTD and RRAS may modulate the phenotypic characteristics of GBM cells by influencing the ferroptosis signaling pathway. By applying ferroptosis activators (RSL3) and inhibitors (Ferrostatin-1, Deferoxamine), we validated the impact of ferroptosis on GBM cell viability. The results demonstrated that the knockdown of DCTD or RRAS significantly reduced the IC₅₀ of RSL3, whereas their overexpression substantially elevated the IC₅₀, further underscoring the critical roles of DCTD and RRAS in ferroptosis regulation. When ferroptosis activators were employed in combination, the groups with DCTD and RRAS knockdown displayed augmented ferroptosis effects. In vivo experiments revealed that the combination of JKE-1647 with the PD-L1 inhibitor atezolizumab effectively suppressed the growth of subcutaneously implanted tumors. Moreover, the integration of DCTD or RRAS knockdown with this combined therapy exhibited even more pronounced inhibitory effects without inducing observable liver toxicity. These findings indicate that DCTD and RRAS not only play pivotal roles in GBM cell proliferation and metastasis but also modulate treatment sensitivity via ferroptosis regulation, potentially providing novel therapeutic targets for GBM in the future.

While our study provides valuable insights into GBM risk profiling and treatment sensitivity, there are several limitations to consider. Firstly, the study’s focus on a limited set of key genes might overlook other important biomarkers and their interactions. Secondly, since the ssGSEA method only measures the relative levels of immune cell infiltration, rather than their absolute abundance, the interpretation of our results has certain limitations. Further experimental validation is required to confirm our conclusions. Additionally, although we identified differential gene expression patterns and their associations with immune cell infiltration, the causal relationships and mechanisms remain unclear. The predictive accuracy of the riskScore model may vary across different patient populations and treatment contexts, and our findings are constrained by the data quality and sample sizes available. Lastly, while predicting drug sensitivity is promising, experimental validation in clinical settings is necessary to confirm these predictions and assess their practical utility in personalized GBM therapy.

5. Conclusion

This study integrated multi-omics data and experimental validation to investigate GBM pathogenesis and potential therapies. Using MR analysis of blood eQTL and GBM GWAS data, we identified 250 GBM-associated genes, and constructed a prognostic risk model based on four key genes (FSTL1, FXYD5, RRAS, RNF216P1). This model effectively stratified patients into high- and low-risk groups, with high-risk patients showing poorer survival, distinct TME features, and differential drug sensitivity. High-risk groups were enriched in immune-related pathways with increased infiltration of T cells, macrophages, and cancer-associated fibroblasts, along with upregulated immune checkpoints. Low-risk groups were linked to ion transport and neurotransmitter regulation pathways. Mutation analyses revealed differential patterns of TP53, EGFR, and PTEN mutations between groups. Experimental validation confirmed that DCTD and RRAS promote GBM cell proliferation, migration, invasion, and in vivo tumor growth by inhibiting ferroptosis via regulating SLC7A11, GPX4, and ACSL4. Combining the ferroptosis inducer JKE-1647 with the PD-L1 inhibitor atezolizumab—enhanced by DCTD/RRAS knockdown—significantly suppressed tumor growth without hepatotoxicity. These findings highlight the risk model’s prognostic value and identify DCTD/RRAS as potential targets for GBM therapy, supporting personalized treatment development.

Availability of Data and Materials

The datasets used and analyzed during the current study are available from the corresponding author on reasonable request.

References

[1] Yang TH, Kon M, Hung JH, Delisi C. Combinations of newly confirmed Glioma-Associated loci link regions on chromosomes 1 and 9 to increased disease risk. BMC Medical Genomics. 2011; 4: 63. https://doi.org/10.1186/1755-8794-4-63.

[2] Shpak M, Hall AW, Goldberg MM, Derryberry DZ, Ni Y, Iyer VR, et al. An eQTL analysis of the human glioblastoma multiforme genome. Genomics. 2014; 103: 252–263. https://doi.org/10.1016/j.ygeno.2014.02.005.

[3] Louis DN, Perry A, Wesseling P, Brat DJ, Cree IA, Figarella-Branger D, et al. The 2021 WHO Classification of Tumors of the Central Nervous System: a summary. Neuro-oncology. 2021; 23: 1231–1251. https://doi.org/10.1093/neuonc/noab106.

[4] De Boeck A, Ahn BY, D’Mello C, Lun X, Menon SV, Alshehri MM, et al. Glioma-derived IL-33 orchestrates an inflammatory brain tumor microenvironment that accelerates glioma progression. Nature Communications. 2020; 11: 4997. https://doi.org/10.1038/s41467-020-18569-4.

[5] Weller M, van den Bent M, Preusser M, Le Rhun E, Tonn JC, Minniti G, et al. EANO guidelines on the diagnosis and treatment of diffuse gliomas of adulthood. Nature Reviews. Clinical Oncology. 2021; 18: 170–186. https://doi.org/10.1038/s41571-020-00447-z.

[6] Qazi MA, Vora P, Venugopal C, Sidhu SS, Moffat J, Swanton C, et al. Intratumoral heterogeneity: pathways to treatment resistance and relapse in human glioblastoma. Annals of Oncology: Official Journal of the European Society for Medical Oncology. 2017; 28: 1448–1456. https://doi.org/10.1093/annonc/mdx169.

[7] Ali MW, Patro CPK, Zhu JJ, Dampier CH, Plummer SJ, Kuscu C, et al. A functional variant on 20q13.33 related to glioma risk alters enhancer activity and modulates expression of multiple genes. Human Mutation. 2021; 42: 77–88. https://doi.org/10.1002/humu.24134.

[8] Chen QR, Hu Y, Yan C, Buetow K, Meerzaman D. Systematic genetic analysis identifies Cis-eQTL target genes associated with glioblastoma patient survival. PloS One. 2014; 9: e105393. https://doi.org/10.1371/journal.pone.0105393.

[9] Wuchty S, Vazquez A, Bozdag S, Bauer PO. Genome-wide associations of signaling pathways in glioblastoma multiforme. BMC Medical Genomics. 2013; 6: 11. https://doi.org/10.1186/1755-8794-6-11.

[10] Mellman I, Chen DS, Powles T, Turley SJ. The cancer-immunity cycle: Indication, genotype, and immunotype. Immunity. 2023; 56: 2188–2205. https://doi.org/10.1016/j.immuni.2023.09.011.

[11] Zhou W, Ke SQ, Huang Z, Flavahan W, Fang X, Paul J, et al. Periostin secreted by glioblastoma stem cells recruits M2 tumour-associated macrophages and promotes malignant growth. Nature Cell Biology. 2015; 17: 170–182. https://doi.org/10.1038/ncb3090.

[12] Chen Z, Feng X, Herting CJ, Garcia VA, Nie K, Pong WW, et al. Cellular and Molecular Identity of Tumor-Associated Macrophages in Glioblastoma. Cancer Research. 2017; 77: 2266–2278. https://doi.org/10.1158/0008-5472.CAN-16-2310.

[13] Amoozgar Z, Kloepper J, Ren J, Tay RE, Kazer SW, Kiner E, et al. Targeting Treg cells with GITR activation alleviates resistance to immunotherapy in murine glioblastomas. Nature Communications. 2021; 12: 2582. https://doi.org/10.1038/s41467-021-22885-8.

[14] Otvos B, Silver DJ, Mulkearns-Hubert EE, Alvarado AG, Turaga SM, Sorensen MD, et al. Cancer Stem Cell-Secreted Macrophage Migration Inhibitory Factor Stimulates Myeloid Derived Suppressor Cell Function and Facilitates Glioblastoma Immune Evasion. Stem Cells (Dayton, Ohio). 2016; 34: 2026–2039. https://doi.org/10.1002/stem.2393.

[15] Raychaudhuri B, Rayman P, Ireland J, Ko J, Rini B, Borden EC, et al. Myeloid-derived suppressor cell accumulation and function in patients with newly diagnosed glioblastoma. Neuro-oncology. 2011; 13: 591–599. https://doi.org/10.1093/neuonc/nor042.

[16] Chongsathidkiet P, Jackson C, Koyama S, Loebel F, Cui X, Farber SH, et al. Sequestration of T cells in bone marrow in the setting of glioblastoma and other intracranial tumors. Nature Medicine. 2018; 24: 1459–1468. https://doi.org/10.1038/s41591-018-0135-2.

[17] Wang Y, Wu N, Li J, Liang J, Zhou D, Cao Q, et al. The interplay between autophagy and ferroptosis presents a novel conceptual therapeutic framework for neuroendocrine prostate cancer. Pharmacological Research. 2024; 203: 107162. https://doi.org/10.1016/j.phrs.2024.107162.

[18] Klemm F, Maas RR, Bowman RL, Kornete M, Soukup K, Nassiri S, et al. Interrogation of the Microenvironmental Landscape in Brain Tumors Reveals Disease-Specific Alterations of Immune Cells. Cell. 2020; 181: 1643–1660.e17. https://doi.org/10.1016/j.cell.2020.05.007.

[19] Nie E, Miao F, Jin X, Wu W, Zhou X, Zeng A, et al. Fstl1/DIP2A/MGMT signaling pathway plays important roles in temozolomide resistance in glioblastoma. Oncogene. 2019; 38: 2706–2721. https://doi.org/10.1038/s41388-018-0596-2.

[20] Zhou F, Tao J, Gou H, Liu S, Yu D, Zhang J, et al. FSTL1 sustains glioma stem cell stemness and promotes immunosuppressive macrophage polarization in glioblastoma. Cancer Letters. 2024; 611: 217400. https://doi.org/10.1016/j.canlet.2024.217400.

[21] Zhang XN, Yang KD, Chen C, He ZC, Wang QH, Feng H, et al. Pericytes augment glioblastoma cell resistance to temozolomide through CCL5-CCR5 paracrine signaling. Cell Research. 2021; 31: 1072–1087. https://doi.org/10.1038/s41422-021-00528-3.

[22] Yu A, Faiq N, Green S, Lai A, Green R, Hu J, et al. Report of safety of pulse dosing of lapatinib with temozolomide and radiation therapy for newly-diagnosed glioblastoma in a pilot phase II study. Journal of Neuro-oncology. 2017; 134: 357–362. https://doi.org/10.1007/s11060-017-2533-6.

[23] Zhou Y, Wu W, Bi H, Yang D, Zhang C. Glioblastoma precision therapy: From the bench to the clinic. Cancer Letters. 2020; 475: 79–91. https://doi.org/10.1016/j.canlet.2020.01.027.

[24] Jin X, Nie E, Zhou X, Zeng A, Yu T, Zhi T, et al. Fstl1 Promotes Glioma Growth Through the BMP4/Smad1/5/8 Signaling Pathway. Cellular Physiology and Biochemistry: International Journal of Experimental Cellular Physiology, Biochemistry, and Pharmacology. 2017; 44: 1616–1628. https://doi.org/10.1159/000485759.

[25] Shi Z, Chen Q, Li C, Wang L, Qian X, Jiang C, et al. MiR-124 governs glioma growth and angiogenesis and enhances chemosensitivity by targeting R-Ras and N-Ras. Neuro-oncology. 2014; 16: 1341–1353. https://doi.org/10.1093/neuonc/nou084.

[26] Gutierrez-Erlandsson S, Herrero-Vidal P, Fernandez-Alfara M, Hernandez-Garcia S, Gonzalo-Flores S, Mudarra-Rubio A, et al. R-RAS2 overexpression in tumors of the human central nervous system. Molecular Cancer. 2013; 12: 127. https://doi.org/10.1186/1476-4598-12-127.

[27] Weller M, Rhun E. How did lomustine become standard of care in recurrent glioblastoma? Cancer Treatment Reviews. 2020; 87: 102029. https://doi.org/10.1016/j.ctrv.2020.102029.

[28] Biasiotta A, D’Arcangelo D, Passarelli F, Nicodemi EM, Facchiano A. Ion channels expression and function are strongly modified in solid tumors and vascular malformations. Journal of Translational Medicine. 2016; 14: 285. https://doi.org/10.1186/s12967-016-1038-y.

[29] Ohta T, Kim YH, Oh JE, Satomi K, Nonoguchi N, Keyvani K, et al. Alterations of the RRAS and ERCC1 genes at 19q13 in gemistocytic astrocytomas. Journal of Neuropathology and Experimental Neurology. 2014; 73: 908–915. https://doi.org/10.1097/NEN.0000000000000110.

[30] Bai Y, Li L, Li J, Lu X. Association analysis of FXYD5 with prognosis and immunological characteristics across pan-cancer. Heliyon. 2024; 10: e30727. https://doi.org/10.1016/j.heliyon.2024.e30727.

[31] Zhao H, Liu H, Kang W, Zhan C, Man Y, Qu T. Analysis on EZH2: mechanism identification of related CeRNA and its immunoassay in hepatocellular carcinoma. BMC Medical Genomics. 2023; 16: 201. https://doi.org/10.1186/s12920-023-01594-9.

[32] van Vlerken-Ysla L, Tyurina YY, Kagan VE, Gabrilovich DI. Functional states of myeloid cells in cancer. Cancer Cell. 2023; 41: 490–504. https://doi.org/10.1016/j.ccell.2023.02.009.

[33] Chen P, Hsu WH, Han J, Xia Y, DePinho RA. Cancer Stemness Meets Immunity: From Mechanism to Therapy. Cell Reports. 2021; 34: 108597. https://doi.org/10.1016/j.celrep.2020.108597.

[34] Friebel E, Kapolou K, Unger S, Núñez NG, Utz S, Rushing EJ, et al. Single-Cell Mapping of Human Brain Cancer Reveals Tumor-Specific Instruction of Tissue-Invading Leukocytes. Cell. 2020; 181: 1626–1642.e20. https://doi.org/10.1016/j.cell.2020.04.055.

[35] Pombo Antunes AR, Scheyltjens I, Lodi F, Messiaen J, Antoranz A, Duerinck J, et al. Single-cell profiling of myeloid cells in glioblastoma across species and disease stage reveals macrophage competition and specialization. Nature Neuroscience. 2021; 24: 595–610. https://doi.org/10.1038/s41593-020-00789-y.

[36] Chen M, Ren R, Lin W, Xiang L, Zhao Z, Shao B. Exploring the oncostatin M (OSM) feed-forward signaling of glioblastoma via STAT3 in pan-cancer analysis. Cancer Cell International. 2021; 21: 565. https://doi.org/10.1186/s12935-021-02260-9.

[37] Wei J, Marisetty A, Schrand B, Gabrusiewicz K, Hashimoto Y, Ott M, et al. Osteopontin mediates glioblastoma-associated macrophage infiltration and is a potential therapeutic target. The Journal of Clinical Investigation. 2019; 129: 137–149. https://doi.org/10.1172/JCI121266.

[38] Cui X, Morales RTT, Qian W, Wang H, Gagner JP, Dolgalev I, et al. Hacking macrophage-associated immunosuppression for regulating glioblastoma angiogenesis. Biomaterials. 2018; 161: 164–178. https://doi.org/10.1016/j.biomaterials.2018.01.053.

[39] Kim JE, Patel MA, Mangraviti A, Kim ES, Theodros D, Velarde E, et al. Combination Therapy with Anti-PD-1, Anti-TIM-3, and Focal Radiation Results in Regression of Murine Gliomas. Clinical Cancer Research: an Official Journal of the American Association for Cancer Research. 2017; 23: 124–136. https://doi.org/10.1158/1078-0432.CCR-15-1535.

[40] Pang L, Khan F, Heimberger AB, Chen P. Mechanism and therapeutic potential of tumor-immune symbiosis in glioblastoma. Trends in Cancer. 2022; 8: 839–854. https://doi.org/10.1016/j.trecan.2022.04.010.

[41] González-Tablas Pimenta M, Otero Á Arandia Guzman DA, Pascual-Argente D, Ruíz Martín L, Sousa-Casasnovas P, et al. Tumor cell and immune cell profiles in primary human glioblastoma: Impact on patient outcome. Brain Pathology (Zurich, Switzerland). 2021; 31: 365–380. https://doi.org/10.1111/bpa.12927.

[42] Kumar V, Patel S, Tcyganov E, Gabrilovich DI. The Nature of Myeloid-Derived Suppressor Cells in the Tumor Microenvironment. Trends in Immunology. 2016; 37: 208–220. https://doi.org/10.1016/j.it.2016.01.004.

[43] Qian BZ, Li J, Zhang H, Kitamura T, Zhang J, Campion LR, et al. CCL2 recruits inflammatory monocytes to facilitate breast-tumour metastasis. Nature. 2011; 475: 222–225. https://doi.org/10.1038/nature10138.

[44] Ozga AJ, Chow MT, Luster AD. Chemokines and the immune response to cancer. Immunity. 2021; 54: 859–874. https://doi.org/10.1016/j.immuni.2021.01.012.

[45] Chun E, Lavoie S, Michaud M, Gallini CA, Kim J, Soucy G, et al. CCL2 Promotes Colorectal Carcinogenesis by Enhancing Polymorphonuclear Myeloid-Derived Suppressor Cell Population and Function. Cell Reports. 2015; 12: 244–257. https://doi.org/10.1016/j.celrep.2015.06.024.

[46] Inamoto S, Itatani Y, Yamamoto T, Minamiguchi S, Hirai H, Iwamoto M, et al. Loss of SMAD4 Promotes Colorectal Cancer Progression by Accumulation of Myeloid-Derived Suppressor Cells through the CCL15-CCR1 Chemokine Axis. Clinical Cancer Research: an Official Journal of the American Association for Cancer Research. 2016; 22: 492–501. https://doi.org/10.1158/1078-0432.CCR-15-0726.

[47] Sceneay J, Chow MT, Chen A, Halse HM, Wong CSF, Andrews DM, et al. Primary tumor hypoxia recruits CD11b+/Ly6Cmed/Ly6G+ immune suppressor cells and compromises NK cell cytotoxicity in the premetastatic niche. Cancer Research. 2012; 72: 3906–3911. https://doi.org/10.1158/0008-5472.CAN-11-3873.

[48] Reichel CA, Puhr-Westerheide D, Zuchtriegel G, Uhl B, Berberich N, Zahler S, et al. C-C motif chemokine CCL3 and canonical neutrophil attractants promote neutrophil extravasation through common and distinct mechanisms. Blood. 2012; 120: 880–890. https://doi.org/10.1182/blood-2012-01-402164.

[49] Veglia F, Perego M, Gabrilovich D. Myeloid-derived suppressor cells coming of age. Nature Immunology. 2018; 19: 108–119. https://doi.org/10.1038/s41590-017-0022-x.

[50] Kamran N, Chandran M, Lowenstein PR, Castro MG. Immature myeloid cells in the tumor microenvironment: Implications for immunotherapy. Clinical Immunology (Orlando, Fla.). 2018; 189: 34–42. https://doi.org/10.1016/j.clim.2016.10.008.

[51] Janjua TI, Rewatkar P, Ahmed-Cox A, Saeed I, Mansfeld FM, Kulshreshtha R, et al. Frontiers in the treatment of glioblastoma: Past, present and emerging. Advanced Drug Delivery Reviews. 2021; 171: 108–138. https://doi.org/10.1016/j.addr.2021.01.012.

Funding

Natural Science Foundation of Fujian Province(2024J08172)

Scientific Research Foundation for the Introduction of Talent of the First Affiliated Hospital of Fujian Medical University(YJRC3991)