Abstract

   The tumor immune microenvironment (TIME) significantly influences
   cancer progression and treatment. This study sought to uncover novel
   TIME-related glioma biomarkers to advance antitumor immunotherapies by
   integrating data from sequencing of bulk RNA as well as scRNA.
   Immunologic and epithelial-mesenchymal transition (EMT) characteristics
   were used to classify glioma patients into two immune subtypes (ISs)
   and two EMT subtypes (ESs). Patients in IS1 and ES1, characterized by
   high immune infiltration and low stemness scores, exhibited poor
   clinical outcomes and limited responsiveness to immunotherapy. A new
   risk signature was developed using 16 genes and validated in
   independent glioma cohorts. Among these, HAVCR2, IL18, LAGLS9, and
   PTPN6 emerged as hub genes, with IL18 identified as a potential
   independent indicator. The upregulation of IL18 in high-grade gliomas
   and U-251 MG cells aligned with bioinformatics analysis. These insights
   deepen the understanding of TIME-related mechanisms in glioma and
   highlight potential therapeutic targets, offering a theoretical
   foundation for effective antitumor immunotherapies in glioma.

   Keywords: Molecular subtypes, Immunologic and epithelial-mesenchymal
   transition gene sets, Tumor immune microenvironment, Immunotherapy

1. Introduction

   Gliomas exist as the most common brain tumor that arises from cells in
   the glia, significantly contributing to morbidity and mortality
   [[27]1]. Glioblastoma (GBM) stands as the most frequent variant of
   malignant glioma, comprising roughly 80 % of primary malignant brain
   tumors. It has a grim prognosis, with under two years of median overall
   survival (OS), significantly diminishing the quality of life [[28]2].
   Despite advances in understanding glioma biology and the development of
   targeted therapies in preclinical and clinical trials, the incidence of
   recurrence, cognitive deficiencies, disability, and mortality
   associated with glioma remains high [[29]3]. Most glioma treatments
   with notable survival benefits, such as radiation and chemotherapy,
   rely on nonspecific targeting of proliferating cells and are often
   compromised by glioma cell invasion and immune evasion [[30]3,[31]4].
   Therefore, it is crucial to identify unknown factors or further
   investigate known factors that influence the tumor microenvironment
   (TME) to promote glioma progression, malignant transformation, and drug
   resistance. This approach is vital for advancing the assessment and
   customization of therapeutic options for gliomas.

   The TME, which is the site of tumor emergence, is significantly
   implicated in its development [[32]5,[33]6]. TME primarily consists of
   the extracellular matrix, various cellular components such as
   fibroblasts and immune components, along with vascular networks that
   are sourced from adjoining tissues [[34][5], [35][6], [36][7]]. With
   the participation of immune-related genes (IRGs) and
   epithelial-mesenchymal transition (EMT)-related genes (ERGs), the TME
   functions crucially in tumor progression and treatment resistance,
   influenced by extracellular matrices, cytokines, growth signals, and
   other elements, thereby affecting cancer treatment [[37]8].

   Compared to brain tissue in the normal state, which exhibits a low to
   moderate immune response, the tumor immune microenvironment (TIME)
   significantly contributes to the aggressive advancement and therapeutic
   resistance in glioma [[38]9]. A suppressive TIME can severely impair
   GBM treatment [[39]10]. Remodeling TIME using a liposomal honokiol and
   disulfiram/copper codelivery system (CDX-LIPO) specifically targeting
   the brain can trigger autophagy and immune-mediated cell death in tumor
   cells. This process activates macrophages and dendritic cells
   infiltrating tumors, as well as the priming of T and NK (natural
   killer) cells, leading to an immunity response against tumors and tumor
   shrinkage [[40]10]. Recently, various IRGs have been identified, which
   can be categorized as low-risk and high-risk [[41][11], [42][12],
   [43][13]], suggesting that IRGs play differential and complex roles in
   the emergence, advancement, and metastatic spread of gliomas.

   Increasing evidence suggests EMT facilitates glioma invasion and
   progression. EMT, a multifactor-regulated process, involves various
   induction agents, families of transcription factors, and an extensive
   range of genes involved in signaling pathways [[44]14]. Enhanced EMT is
   associated with poor survival, whereas the EMT in the initial stages
   iscomparatively less malignant in GBM [[45]15], suggesting that EMT
   assessment could predict GBM prognosis and aid in the advancement of
   innovative therapeutic strategies and prognostic markers [[46]15].
   Additionally, the elevated expression of EMT genes in GBMs and
   pilocytic astrocytomas is indicative of mural cells proliferation
   within newly formed blood vessels [[47]16].

   Six collagen genes have been identified as regulators of both the
   glioma immunosuppressive tumor microenvironment and the EMT process,
   namely COL1A1, COL1A2, COL3A1, COL4A1, COL4A2, and COL5A2 [[48]17].
   This suggests that IRGs and ERGs may synergistically influence glioma
   behaviors such as invasion, migration, and metastasis. Cohort analysis
   of The Cancer Genome Atlas (TCGA) uncovered two immune subtypes (ISs)
   and two EMT subtypes (ESs) in patients with glioma based on IRG and ERG
   sets. This resulted in the elucidation of key immune- and EMT-related
   genes (IERGs), with sixteen chosen to construct a novel risk prediction
   model for glioma patients. Among these, four IERGs (HAVCR2, IL18,
   LAGLS9, and PTPN6) emerged as hub genes, showing high expression in M1
   macrophages and monocytes. Furthermore, the upregulation of IL18 in GBM
   indicates its potential as an independent clinical marker for patients
   with high-grade glioma. This integrated analysis of IRGs and ERGs in
   gliomas establishes a theoretical framework aimed at developing and
   implementing more effective antitumor immunotherapies.

2. Results

2.1. Identification of novel ISs and ESs of gliomas and analysis of their
survival and different immune features

   The study framework is detailed in [49]Supplementary Fig. 1. Consensus
   clustering without supervision (k = 2) was applied to 663 gliomas from
   the TCGA cohort, focusing on immunologic and EMT gene signatures. This
   analysis categorized patients with glioma into two ISs (IS1 and IS2)
   and two ESs (ES1 and ES2) ([50]Fig. 1A–H). Based on Kaplan-Meier
   survival analysis, individuals in the IS1 and ES1 subtypes had better
   OS compared to those in IS2 and ES2 ([51]Fig. 1B–I). Gene set
   enrichment analysis (GESA) of HALLMARK pathways revealed significant
   differential enrichment in five pathways—interferon-gamma response,
   IL2-STAT5 signaling, apoptosis, hypoxia, and glycolysis—in IS1
   ([52]Fig. 1C) and ES1 ([53]Fig. 1J) based on NES, P-value, and FDR
   value.

Fig. 1.

   [54]Fig. 1
   [55]Open in a new tab

   Unsupervised consensus clustering of gliomas based on immunologic and
   EMT gene signatures in the TCGA cohort. (A and H) Heatmap depicting
   consensus clustering matrix (k = 2) for 1959 immunologic genes (A) and
   1263 EMT genes (H) across 663 glioma samples. (B and I) Kaplan-Meier
   survival analysis of patients with gliomas exhibiting different ISs (B)
   and ESs (I). (C and J) GSEA for pathways in HALLMARK terms of different
   ISs (C) and ESs (J). (D and K) Potential ICB response predicted with
   TIDE algorithm in different ISs (D) and ESs (K). (E and L) Stemness
   score based on OCLR algorithm in different ISs (E) and ESs (L). (F-G
   and M-N) TMB and MSI score in different ISs (F-G) and ESs (M-N).

   The TIDE algorithm, used to predict potential immunotherapeutic
   responses, showed significantly higher TIDE scores in IS1 and ES1
   compared to IS2 and ES2 ([56]Fig. 1D–K), indicating that IS1 and ES1
   patients had poorer responses to immune checkpoint blockade (ICB)
   immunotherapy and experienced decreased survival following therapies.
   The OCLR algorithm assessed stemness scores, showing higher scores in
   IS2 and ES2 than in IS1 and ES1 ([57]Fig. 1E–L), which correlates with
   greater tumor differentiation in IS2 and ES2 patients. Tumor mutational
   burden (TMB) and microsatellite instability (MSI) scores, critical for
   predicting immune efficacy, were higher in IS1 and ES1 compared to IS2
   and ES2 ([58]Fig. 1F–M) and lower in MSI score ([59]Fig. 1G–N).
   Interestingly, IS1 and ES1 exhibited similar survival prognoses and
   clinicopathological features.

2.2. TIME features of distinct TCGA cohort subtypes

   To evaluate the association of glioma subtypes within immune scores,
   the ESTIMATE algorithm evaluated stromal and immune cell infiltration,
   along with tumor purity, across different ISs and ESs. This analysis
   generated four scores: immune, stromal, ESTIMATE, and tumor purity. The
   results indicated IS1 and ES1 showed increased immune scores
   ([60]Supplemental Figs. 2A and E), stromal scores ([61]Supplemental
   Figs. 2B and F), and ESTIMATE scores ([62]Supplemental Figs. 2C and G),
   along with lower tumor purity ([63]Supplemental Figs. 2D and H)
   compared to IS2 and ES2. Higher tumor purity in IS2 and ES2 suggested
   fewer infiltrated immune and stromal cells in these gliomas. The
   infiltration profiles of 22 distinct immune cell types were assessed by
   the CIBERSORT algorithm across various ISs and ESs, and the results
   were displayed using heatmaps ([64]Supplemental Figs. 2I and K). Clear
   differences were noted in various populations of immune cells in the
   comparison of IS1 and IS2, and between ES1 and ES2. Notably, the
   composition of B cell plasma, resting T cells CD4^+ memory, monocytes,
   and M2 macrophages varied significantly between IS1 and IS2 and between
   ES1 and ES2 ([65]Supplemental Figs. 2J and L).

2.3. Linkage between ISs and ESs with ICPs, ICD modulators, and HLA family

   Considering the critical roles of immune checkpoint proteins (ICPs),
   immunogenic cell death (ICD) modulators, and HLA family genes in cancer
   immunity and the effectiveness of vaccines against mRNA, distinct
   expression of these elements was assessed in different ISs and ESs. The
   levels of expression for 25 ICD modulators, 47 ICPs, and 24 HLA family
   genes were investigated across different ISs and ESs. The analysis
   revealed that 19 ICD modulators, 39 ICPs, and 23 HLA-related genes were
   differentially expressed in these subtypes ([66]Supplemental Fig. 3).
   For instance, 18 ICD modulators were significantly upregulated in IS1
   and ES1 gliomas. Additionally, 38 ICPs, excluding ADORA2A, CD27, IDO2,
   KIR3DL1, TIGIT, TNFRSF25, TNFSF18, TNFSF9, and VSIR, showed significant
   upregulation in IS1 and ES1 gliomas ([67]Supplemental Fig. 3). All
   HLA-related genes, except HLA.J, exhibited significantly higher
   expression levels in IS1 and ES1 gliomas ([68]Supplemental Fig. 3).
   These findings suggest that IS1 and ES1 gliomas share similar
   expression profiles of ICPs, ICD modulators, and HLA family genes,
   indicating that patients with IS1 and ES1 gliomas may benefit from
   similar mRNA vaccine therapies.

2.4. Exploring Time-related modules by WGCNA

   To recognize IERGs, WGCNA on the gliomas-TCGA cohort constructed
   co-expression networks and pinpointed gene modules linked with immune
   scores. Setting the soft threshold to 12 ([69]Fig. 2A), the analysis
   identified 20 key TIME-related gene modules, each color-coded ([70]Fig.
   2B). The eigengene adjacency heatmap illustrated correlations among
   these modules ([71]Fig. 2C). Further exploration revealed their
   associations with clinical factors, such as grade, age, gender, and
   various immune-related parameters ([72]Fig. 2D). The brown module,
   showing the strongest links to stromal, immune and ESTIMATE score, was
   selected to conduct deeper analysis ([73]Fig. 2D). Scatter plots
   demonstrated significant positive associations of the brown module
   membership with immune score ([74]Fig. 2E), stromal score ([75]Fig.
   2F), ESTIMATE score ([76]Fig. 2G), and tumor purity ([77]Fig. 2H).

Fig. 2.

   [78]Fig. 2
   [79]Open in a new tab

   WGCNA screening for TIME-related module in the TCGA cohort. (A)
   Analysis of the scale-free topology fit index and mean connectivity for
   various soft-threshold powers (β), selecting a power of 12. (B)
   Dendrogram of all genes clustered based on a dissimilarity measure
   (1-TOM). (C) Correlation heatmap of different modules. (D)
   Determination of module-trait relationships in patients with glioma.
   (E-H) Scatter plots of module eigengene significance for immune score
   (E), stromal score (F), ESTIMATE score (G), and tumor purity (H) in the
   brown module. (For interpretation of the references to color in this