An immune and epithelial–mesenchymal transition-related risk model and immunotherapy strategy for grade II and III gliomas

Grade II and III gliomas are heterogeneous and aggressive diseases. More efficient prognosis models and treatment methods are needed. This study aims to construct a new risk model and propose a new strategy for grade II and III gliomas. The data were downloaded from The Cancer Genome Atlas (TCGA), the Gene Expression Omnibus (GEO), gene set enrichment analysis (GSEA), and the EMTome website for analysis. The Human Cell Landscape website and the Genomics of Drug Sensitivity in Cancer website were used for single-cell analysis and drug susceptibility analysis. Gene set enrichment analysis, gene function enrichment analysis, univariate and multivariate Cox regression analyses, Pearson’s correlation analysis, log-rank test, Kaplan–Meier survival analysis, and ROC analysis were performed. We constructed an immune-related prognostic model associated with the isocitrate dehydrogenase 1 (IDH1) mutation status. By analyzing the immune microenvironment of patients with different risk scores, we found that high-risk patients were more likely to have an inflammatory immune microenvironment and a higher programmed death ligand-1 (PD-L1) expression level. Epithelial–mesenchymal transition (EMT)-related gene sets were significantly enriched in the high-risk group, and the epithelial–mesenchymal transition phenotype was associated with a decrease in CD8+ T cells and an increase in M2 macrophages. Transforming growth factor-β (TGF-β) signaling was the most important signaling in inducing epithelial–mesenchymal transition, and TGFB1/TGFBR1 was correlated with an increase in CD8+ T cytopenia and M2 macrophages. Survival analysis showed that simultaneous low expression of TGFBR1 and PD-L1 had better survival results. Through single-cell analysis, we found that TGFB1 is closely related to microglia and macrophages, especially M2 macrophages. Finally, we discussed the sensitivity of TGFB1 inhibitors in gliomas using cell line susceptibility data. These results demonstrated a potential immunotherapy strategy in combination with the TGFB1/TGFBR1 inhibitor and PD-1/PD-L1 inhibitor for grade II and III gliomas.


Introduction
Grade II and III gliomas are the most common primary brain tumors and proved to be with substantial heterogeneity in terms of pathological features and clinical outcomes (Louis et al., 2021). In order to distinguish the different pathological features of patients and develop individualized treatment strategies, glioma-related biomarkers have been identified. The IDH1 was reported to mutate frequently in gliomas (Yan et al., 2009). It had been proved that gliomas with the IDH1 mutation were more sensitive to chemotherapy and radiotherapy, resulting in a better prognosis (Sanson et al., 2009). Other biomarkers, for example, codeletion of 1p and 19q (1p/19q) (Smith et al., 2000), Capicua (CIC) transcriptional repressor mutation (Gleize et al., 2015), loss of chromosome 9p, mutation of phosphatidylinositol-4,5bisphosphate 3-kinase catalytic subunit alpha (PIK3CA) and phosphoinositide-3-kinase regulatory subunit 1 (PIK3R1) (Draaisma et al., 2015), and deletion of cyclin-dependent kinase inhibitor 2 A (CDKN2A) (Reis et al., 2015), were confirmed to have prognostic value and important for rational selection of surgery, radiotherapy, and chemotherapy treatment. However, over 50% of grade II and III gliomas eventually develop into highly aggressive gliomas, indicating the need for a more efficient prognosis model and treatment methods (Zhu et al., 2021).
In the past decades, there has been little progress in the treatment of gliomas. Although immunotherapy successfully promoted the treatment results of other cancer types and was the major research direction for gliomas, limited progression had been made in immunotherapy treatment of gliomas (Yang et al., 2022). In CheckMate-498, a randomized clinical trial, nivolumab combined with bevacizumab and nivolumab combined with chemoradiotherapy in newly diagnosed glioma patients with O 6methylguanine DNA methyltransferase (MGMT) promoter unmethylation were both ineffective (Reardon et al., 2020). Comparing with IDH-wild gliomas, the IDH-mutant gliomas have significantly low tumor-infiltrating lymphocytes and PD-L1 expression (Berghoff et al., 2017). Numerous studies had demonstrated that gliomas were infiltrated by immune cells that made up to 30% of a tumor's mass (Kaminska et al., 2021). Such extensive accumulation of innate immune cells in gliomas might be misleading as these events did not reflect the effective anti-tumor immunity. This phenomenon enlightened us that tumor immune infiltration in gliomas might be accompanied by other pathological processes that promote tumor progression.
As mentioned previously, the construction of a new risk model and improving the efficacy of immunotherapy are essential for the treatment of gliomas. In this study, we attempted to construct an immune-related risk model and proposed a feasible strategy for improving the immunotherapy efficacy for gliomas.

RNA sequencing data
The IDH1 somatic mutation status for 500 samples, gene expression data for 525 samples, and the corresponding clinical datasheets for 515 samples were obtained from TCGA website (https://portal.gdc.cancer.gov/). Among these grade II and III glioma samples, 493 samples with RNA sequencing data and IDH1 mutation information were subjected to subsequent analyses. Log2 (x+1) normalization was performed for all gene expression data. Rows and columns with more than 50% missing values were removed. The study report fully met TCGA publication requirements.

Construction and validation of the risk model
GSEA was performed to determine how the immunological pathways and corresponding immune genes differ between IDH1-wild (IDH1-wt) (n = 117) and IDH1 mutation (IDH1mut) (n = 376) samples in TCGA cohort. An annotated gene set file (c7. immunesigdb.v7.4. symbols.gmt; downloaded from the Molecular Signatures Database) was selected for use as the reference gene set. The expression profiles of the top 50 genes expressed in the IDH1-wt and IDH1-mut groups were analyzed via univariate Cox regression analysis. In this analysis, genes were regarded as significant at p < 0.05. For the significant genes, least absolute shrinkage and selection operator (LASSO) Cox analysis was performed by using the glmnet R package. Then, the nine candidate genes were analyzed by multivariate Cox regression analysis based on progression-free survival (PFS). Finally, two independent prognostic factors for PFS were analyzed by multivariate Cox regression analysis based on overall survival (OS) to construct the risk model (Risk Score i Coefficient of(i)*Expressionofgene(i)). The log-rank test and Kaplan-Meier survival analysis were used to assess the predictive ability of the prognostic model. The maxstat R package was applied to determine the best cutoff value, Frontiers in Genetics frontiersin.org 02 and then patients were classified into low-risk and high-risk groups. We performed ROC analysis using the R package pROC.

Correlation of immunophenotype with the risk model
Using the deconvo_CIBERSOR method of the IOBR R package, the immune cell infiltration score was calculated. The R package ESTIMATE was used to calculate the immune scores for each patient based on gene expression. Pearson's correlation coefficients for risk scores and immune infiltration were calculated. The rank sum test was used to compare differences in immune cells and immune-related molecules between highrisk and low-risk groups. p < 0.05 was considered significant.

Functional enrichment analysis
The STRING website (https://cn.string-db.org/) was used to construct a 32-gene functional protein association network. We used the Kyoto Encyclopedia of Genes and Genomes Application Programming Interface (KEGG API) to obtain genetic annotations for the latest KEGG pathway. Enrichment analysis was performed using the R package clusterProfiler.

Epithelial-mesenchymal transition (EMT) gene set-related analysis
The 61 EMT gene sets were downloaded from the EMTome website (http://emtome.org/) (Vasaikar et al., 2021). GSEA was performed between the high-risk and low-risk groups. The enrichment scores of the EMT gene sets were analyzed by univariate Cox analysis through the EMTome website. The gene set with the highest prognostic significance and the highest GSEA enrichment score was selected for gene set variation analysis (GSVA) to obtain the enrichment scores of samples. Pearson's correlation coefficients for enrichment scores, risk scores, immune cells, and immune-related molecules were calculated. The log-rank test and Kaplan-Meier survival analysis were applied to assess the predictive ability of the enrichment score. By reviewing the literature, we identified five EMT phenotype-related signaling pathways and downloaded corresponding gene sets from the GSEA website (http://www.gsea-msigdb.org/), namely, REACTOME_ SIGNALING_BY_TGFB_FAMILY_MEMBERS, GOBP_ CANONICAL_WNT_SIGNALING_PATHWAY, BIOCARTA_ RAS_PATHWAY, GOBP_NOTCH_SIGNALING_PATHWAY, and GOBP_PHOSPHATIDYLINOSITOL_3_KINASE_ SIGNALING (Dongre and Weinberg, 2019). We divided the samples into two groups according to the enrichment score, with 263 samples in the high-enrichment score group and 262 samples in the low-enrichment score group, and performed GSEA to obtain the most significant enrichment signaling pathway.
2.7 Effects of transforming growth factorβ (TGF-β)-associated molecules, PD-L1 (CD274), and CTLA-4 on immune cell infiltration and prognosis Pearson's correlation coefficients for enrichment scores, TGFB1, TGFBR1, TGFB2, and TGFB3 were calculated, and the correlation matrix was plotted. Multivariate Cox regression analysis was performed on TGFB1, TGFBR1, TGFB2, TGFB3, PD-L1, and cytotoxic T-lymphocyteassociated protein 4 (CTLA-4). For TGFB1, TGFBR1, CD274, CTLA-4, and immune cells, we calculated the Pearson's correlation coefficient and plotted the correlation scatterplot. The log-rank test and Kaplan-Meier survival analysis were used to assess the predictive ability of TGFBR1 and PD-L1. The maxstat R package was applied to determine the best cutoff value, and then patients were classified into the low-expression group and the highexpression group. The samples were then divided into four groups (TGFBR1-H+CD274-L, TGFBR1-H+CD274-H, TGFBR1-L+CD274-L, and TGFBR1-l+CD274-H) for survival analysis, and the differences in immune cell infiltration between groups were compared by the rank sum test, and violin charts were plotted.
2.8 Single-cell analysis was performed to determine cell localization of TGF-β signaling-related molecules This analysis was conducted through the Human Cell Landscape website (https://db.cngb.org/HCL/index.html) (Han et al., 2020). Platform creators analyzed >700,000 single cells from >50 human tissues (2-4 replicates per tissue in general) and cultures. Through the brain section of the Gallery module, we can acquire the single-cell data matrix and analysis results related to brain tissue.

Anti-TGFB1 drug susceptibility analysis
We adopted the Genomics of Drug Sensitivity in Cancer website (https://www.cancerrxgene.org/) and the EMTome website (http://emtome.org/) for drug sensitivity analysis. The two websites provided drug susceptibility data on the TGFB1 inhibitor LY2109761 in different cell lines, as well as online analysis tools. Using online tools from both websites, we performed drug susceptibility analysis.
Frontiers in Genetics frontiersin.org 2.10 Data processing platform All the data processing was performed on Sangerbox (http://vip. sangerbox.com/home.html), a powerful platform based on R (Shen et al., 2022), including RNA sequencing data normalization, merging of datasets, removing batch effects, gene set enrichment analysis, gene set variation analysis, gene function enrichment analysis, univariate and multivariate Cox regression analyses, LASSO-Cox regression analysis, Pearson's correlation analysis, log-rank test, Kaplan-Meier survival analysis, ROC analysis, the immune cell infiltration score and immune infiltration score calculation, and KEGG enrichment analysis.

Construction and validation of a risk model associated with IDH1 and immune status
We applied GSEA between IDH1-wt (n = 117) and IDH1mut (n = 376) groups using c7. immunesigdb.v7.4. symbols.gmt as the reference gene set. All 3919 immunerelated gene sets with FDR < 0.05 were enriched in the IDH1-wt group, suggesting that the IDH1-wt group was more correlated with immune response than the IDH1-mut group. The top 50 genes enriched in the two groups were used for univariate Cox regression analysis. The results showed that all 100 genes were significantly associated with prognosis ( Figure S1). These 100 genes were then analyzed by LASSO Cox, and nine candidate genes were obtained for risk modeling ( Figures 1A-B). In order to better predict PFS, we performed multivariate analysis on PFS for nine genes and obtained two genes with independent prognostic value for PFS, NOG, and IGFBP2 ( Figure 1C). Through multivariate analysis of NOG and IGFBP2 in OS, we obtained the expression coefficients of these two genes and constructed a risk model ( Figure 1D). Risk score = -0.4253 * expression of NOG + 0.3954 *expression of IGFBP2.
We calculated the risk scores of all 493 TCGA samples and performed survival analysis, plotting K-M curves and ROC curves. The results showed that the risk model was a good predictor of patients' OS (p = 1.3e-24, HR = 5.52) and PFS (p = 3.2e-13, HR = 3.20) (Figures 2A-D). Clinically, the Frontiers in Genetics frontiersin.org IDH1 mutation status was used to determine patient prognosis, and we explored whether our risk model could further stratify the IDH1 mutation status. In both the IDH1-mut and IDH1-wt groups, a high-risk score indicates a worse prognosis (Figures 2E-F), whereas the IDH1 status in the low-risk group does not affect prognosis ( Figure 2G). Multivariate analysis of IDH1 and risk score suggested that risk score is an independent prognostic factor relative to IDH1 ( Figure 2N). GSE107850 (including 195 samples) and GSE43388 (including 43 samples) were used to validate the risk prediction capabilities of risk models. Among them, the study endpoint of the GSE107850 sample was PFS after radiotherapy (RT) and temozolomide (TMZ) treatment, and the study endpoint of the GSE43388 sample was OS, so we used these two sets of sample data to verify the predictive efficacy of risk models on PFS and OS, respectively. Furthermore, we used the GSE107850 dataset to analyze the predictive power of risk models for response to radiotherapy and chemotherapy treatment. As shown in Figures 2H-M, the risk model can predict not only the patient's OS (p = 3.2e-6, HR = 7.89) and PFS (p = 8.8e-8, HR = 2.84) but also the patient's response to radiation (p = 7.7e-6, HR = 3.45) and chemotherapy (p = 1.9e-3, HR = 2.39).

EMT-related signaling pathways are closely related to the risk model
In order to explore the signaling pathways related to the risk model, we first used the STRING website to construct a gene regulatory network of 32 genes for the risk model ( Figure 4A), and then we performed KEGG enrichment analysis on this regulatory network ( Figure 4B). Some signaling pathways that are closely related to EMT had been enriched, including TGF-β signaling pathway (Hao et al., 2019), signaling pathways regulating pluripotency of stem cells (Mani et al., 2008;Shibue and Weinberg, 2017), Hippo signaling pathway (Cordenonsi et al., 2011), PI3K-Akt signaling pathway, and mTOR signaling pathway (Song et al., 2014). This suggested that the EMT phenotype may be also formed at the same time as the tumor formed an inflammatory immune microenvironment.
3.4 EMT phenotypes are significantly associated with M2 macrophages and are regulated by the TGF-β signaling pathway The EMTome website provided 61 gene sets that are associated with the EMT phenotype, and we downloaded all of them for GSEA between the two risk groups. At the same time, we used the online tool of the EMTome website to perform univariate Cox regression analysis of each gene set enrichment score. The top gene sets are listed in Table 1; Table 2. PMID29212455: Wang_et_al. 2017 was the gene set with the most significant prognostic value and highest normalized enrichment score (NES), and we plotted the GSEA enrichment curve for this gene set ( Figure 5A).
In order to determine the main signaling pathways that induce EMT phenotypes, we reviewed the relevant literature, identified five candidate signaling pathways, downloaded the gene sets of each pathway through the GSEA website, and   Note: EMT gene sets were acquired from the EMTome website, and GSEA was performed between the high-risk and low-risk groups. FDR <0.05 was considered significantly different in gene set enrichment between the two groups.

Compared with PD-L1 and CTLA-4, TGFB1/TGFBR1 is much associated with the immunosuppressive microenvironment
In this section, we explored specific molecules of the TGF-β signaling pathway, and we selected TGFB1, TGFBR1, TGFB2, and TGFB3 as research subjects because the drugs targeting these molecules are currently in the clinical research stage (Tauriello et al., 2022), and exploring them will make our research conclusions more likely to guide clinical practice. First, we analyzed the correlation between four candidate molecules and EMT enrichment scores, and the results showed that all four molecules were significantly correlated with EMT enrichment scores, among which TGFBR1 had the highest correlation (p < 0.001, r = 0.55) ( Figure 6A). We then performed multivariate Cox regression analysis on four candidate molecules together with PD-L1 and CTLA-4, and the results suggested that PD-L1, TGFB2, and TGFBR1 were independent prognostic factors ( Figure 6B). Based on the previous conclusions, we selected TGFB1/TGFBR1 for subsequent analysis.

Simultaneous blocking of TGFB1/ TGFBR1 and PD-L1 might significantly improve survival
We calculated the optimal cutoff value for TGFBR1 and PD-L1 using the R package maxstat (maximally selected rank statistics with several p-value approximations, version: 0.7-25). We then grouped the samples according to cutoff values and performed survival analysis. Patients with higher PD-L1 and TGFBR1 expression had shorter survival (p = 1.4e-9, HR = 2.92; p = 1.0e-7, HR = 2.58) (Figures 8A-B). Patients with both high expression of TGFBR1 and PD-L1 had the worst prognosis, and those with low expression of both TGFBR1 and PD-L1 had the best prognosis, while those with high expression of one of the two molecules had a moderate prognosis ( Figure 8C; Supplementary  Table S1). The group with the best prognosis tended to have more CD4 + T naïve cells and CD8 + T-cell infiltration (although no statistically significant difference was achieved) and less  Figures 8D-H). These results suggested that simultaneous blocking of TGFB1/ TGFBR1 and PD-L1 is more likely to confer a survival benefit than blocking PD-L1 alone.

TGFB1 is closely related to microglia and macrophages in brain tissues
We analyzed the cellular localization of TGF-β-related molecules in brain tissue through Human Cell Landscape, a single-cell analysis website. In four single-cell samples of brain tissue, we found that TGFB1 and TGFBI (transforming growth factor-beta induced, a protein induced by TGFB1) were mainly expressed in microglia and macrophages, especially M2 macrophages. The results of the clustering analysis with marker genes are shown in Figures 9A-D and Table 3.

Gliomas may be one of the tumors most sensitive to the TGFB1 inhibitor
LY2109761 is an inhibitor of TGFB1 and is widely used in in vitro studies of various tumors. Through the Genomics of Drug Sensitivity in Cancer website and the EMTome website, we analyzed the susceptibility results of LY2109761 in various tumors. Figure 10A shows the IC 50 values of LY2109761 in different LGG cell lines, and the IC 50 values of LY2109761 in all available tumor cell lines are shown in Figure 10B. After standardizing the susceptibility data of all tumor cell lines, we found that LGG ranked sixth in drug sensitivity among 29 tumor species ( Figure 10C). This suggested that TGFB1/ TGFBR1 inhibitors may be sensitive in gliomas.

Discussion
Previous studies had demonstrated that the IDH1 mutation status has a significant effect on the immune microenvironment of diffuse glioma, and IDH1 wild-type patients generally have higher lymphocyte infiltration and PD-L1 expression (Yan et al., 2009;Berghoff et al., 2017). Also, the IDH1 mutation status has a significant correlation with prognosis (Sanson et al., 2009). Therefore, we can utilize these features of IDH1 to construct a risk model related to the inflammatory immune microenvironment and survival. Through GSEA, we found that the 3,919 immune-related gene sets with significant differences between the two groups were enriched in the Frontiers in Genetics frontiersin.org IDH1-wild group, suggesting that the IDH1-wild type is more related to immune response, which is consistent with the conclusions mentioned previously. According to the gene rank list of GSEA, top 50 genes in each group were selected for univariate analysis, and the results showed that these 100 genes had significant prognostic value. By performing LASSO Cox and multivariate Cox regression analyses on 100 genes, we constructed a two-gene risk model consisting of IGFBP2 and NOG. IGFBP2 is a pleiotropic oncogene and plays a role in the occurrence and development of a variety of tumors (Brouwer-Visser and Huang, 2015). It has been confirmed that IGFBP2 can induce tumor epithelial-mesenchymal transformation and metastasis through the NF-κB signaling pathway (Gao et al., 2016). NOG is a natural inhibitor of bone morphogenetic protein (BMP), especially BMP2 and BMP4, which are the members of the TGF-β family. BMPs suppressed the tumorigenic function of human glioma-initiating cells by inducing cell differentiation, cell cycle arrest, and apoptosis (Bao et al., 2013). Several reports have shown that BMP4 is expressed in low-grade gliomas and that it serves as a favorable prognostic marker in gliomas (Bao et al., 2013;Nayak et al., 2020;Zhou et al., 2020). In addition, BMP4 was able to abolish cancer stem cell populations in human cancers, including malignant gliomas Piccirillo and Vescovi, 2006). Interestingly, as an inhibitor of BMPs, the expression of NOG and BMP2/4 was positively correlated. The high-risk group had lower levels of BMPs, consistent with the studies mentioned previously. We verified that the risk model has good prognostic value in TCGA cohort and two GEO cohorts and can further stratify the IDH1-mut and IDH1-wt groups. The multivariate Cox regression analyses further confirmed the independent prognostic value of the  Note: Single-cell differential gene expression analysis for four brain tissues was performed on the Human Cell Landscape website; p < 0.05 was considered statistically significant.
Frontiers in Genetics frontiersin.org risk model. Thus, this model might be used as a risk indicator in clinical practice.
To further investigate the relationship between the risk model and the immune microenvironment, we analyzed the correlation between immune infiltration scores, immunerelated molecules, and immune cell infiltration with immunerelated prognostic models. The stromal score, immune score, and ESTIMATE score were positively correlated with the risk score. The high-risk group had higher expression of immunosuppressive molecules and also had higher expression of immunostimulant molecules. Tregs were elevated in the highrisk group. However, B naïve cells, mast resting cells, CD4+ T memory resting cells, CD8+ T cells, neutrophils, and M0 and M1 macrophages were also elevated in the high-risk group. Similar conclusions had been reached in other people's studies. Berghoff et al. (2017) found significantly higher levels of PD-1-positive tumor-infiltrating lymphocytes and PD-L1 expression in IDH-wild-type gliomas than IDH-mutant gliomas. Liu et al. (2020) demonstrated higher CTLA-4 expression in higher-grade IDH-wild-type tumors than lowergrade IDH-mutant tumors. We would summarize these phenotypes as "hot" tumor microenvironments. "Hot" tumors and high expression of PD-L1 were known as hallmarks of sensitivity to immunotherapy (Xiong et al., 2018;  Frontiers in Genetics frontiersin.org 2022), while inexplicably, the efficiency of immunotherapy was limited in gliomas, regardless of the IDH1 status and tumor grades (Blumenthal et al., 2016;Bouffet et al., 2016;Johanns et al., 2016;Lukas et al., 2018;Reardon et al., 2020). We hypothesized that the patient's inflammatory immune microenvironment is accompanied by other malignant phenotypes. We then established a 32-gene functional protein association network for NOG and IGFBP2 on the STRING website. These 32 genes were analyzed for gene function enrichment using the KEGG database. The TGF-β signaling pathway, Hippo signaling pathway, PI3K-AKT-mTOR pathway, and signaling pathways regulating pluripotency of stem cells were found to be closely linked to the risk model. It had been reported that stem cell properties could be acquired by tumor cells through EMT. Induction of EMT in immortalized human mammary epithelial cells was sufficient to induce the expression of stem cell markers, enhance self-renewal, and increase the number of tumor-initiating cells (Mani et al., 2008;Shibue and Weinberg, 2017). TGF-β had been thought to be the most important factor inducing EMT via the classic Smad and non-Smad pathways (Lamouille et al., 2014;Hao et al., 2019). The Hippo pathway and PI3K-AKT-mTOR pathway had also been proved to be related to EMT (Cordenonsi et al., 2011;Song et al., 2014). It had been reported that there exists a strong correlation between EMT and immune activation. Further analysis demonstrated high expression of immune checkpoints and other druggable immune targets such as PD-1, PD-L1, CTLA-4, OX40L, and PD-L2 in patients with the EMT phenotype (Mak et al., 2016). Therefore, we speculated that the inflammatory immune microenvironment of gliomas is accompanied by the EMT phenotype.
In subsequent analyses, we found that EMT enrichment scores were significantly positively correlated with risk scores, M2 macrophage infiltration, Tregs, PD-L1, and CTLA-4 expression and negatively correlated with CD8 + T-cell infiltration. Although Mak et al. (2016) discovered that the EMT phenotype is always accompanied with immune activation and higher expression of immune checkpoint molecules and declared that immunotherapy might have potential, the reality was far from that (10). The inflammatory immune microenvironment is accompanied by EMT, which in turn induces immunosuppression against the inflammatory immune microenvironment. Reversing the EMT phenotype might be necessary for immunotherapy treatment.
Dongre and Weinberg (2019) reviewed the main signaling pathways that induce EMT, including the TGF-β signaling pathway, the WNT signaling pathway, the Notch signaling pathway, the PI3K signaling pathway, and the RAS signaling pathway. Through the enrichment analysis of the aforementioned pathway, we found that the TGF-β signaling pathway is the most important. TGFB1/TGFBR1, TGFB2, and TGFB3 were elected for further investigation. Among these four molecules, TGFBR1 not only had the highest correlation with EMT enrichment scores but also was an independent prognostic factor relative to PD-L1. Immunocorrelation analysis showed that TGFB1/ TGFBR1 had more powerful immunosuppressive properties than PD-L1 and CTLA-4, especially in inducing M2 macrophage infiltration and CD8 + T-cell depletion. According to the different expressions of TGFBR1 and PD-L1, we performed survival analysis in groups, and the results showed that patients with both TGFBR1 and PD-L1 expression had obvious survival advantages, and the high expression of either molecule led to poor prognosis. This indicates that the combined inhibition of TGFB1/ TGFBR1 and PD-1/PD-L1 has a good clinical application prospect.
Through single-cell analysis, we further determined that TGFB1 and TGFBI are mainly derived from microglia and M2 macrophages. As resident macrophages of the central nervous system (CNS), microglia are associated with diverse functions essential to the developing and adult brain during homeostasis and disease (Borst et al., 2021). Microglia-derived TAM (tumor-associated macrophages) increased angiogenesis and suppressed T-cell proliferation. Depletion of TAM provides survival advantages and delays recurrence when combined with standard-of-care treatment such as irradiation (Akkari et al., 2020). Numerous studies have demonstrated that gliomas are infiltrated by immune cells that make up to 30% of a tumor's mass (Nduom et al., 2015). The predominant population consists of glioma-associated microglia and macrophages, and their numbers inversely correlate with patients' survival (Gieryng et al., 2017). We speculated that by synthesizing and secreting TGFB1, microglia and M2 macrophages simultaneously induced EMT phenotype and immunosuppression.
Finally, we explored the relative drug sensitivity of the TGFB1 inhibitor in glioma cell lines through the Genomics of Drug Sensitivity in Cancer website and the EMTome website. Although these studies are in vitro experiments, the relative sensitivity between different tumor species can still give us some hints that gliomas have relatively good sensitivity relative to most tumors.
In this study, we constructed an immune-related prognostic model associated with the IDH1 mutation status. This model enables further risk stratification of patients with different IDH1 mutation states. By analyzing the immune microenvironment of patients with different risk scores, we found that high-risk patients were more likely to have an inflammatory immune microenvironment and a higher PD-L1 expression level, although clinical studies showed that patients with different IDH1 mutation states did not benefit from PD-1/ PD-L1 inhibitors. We speculated that there were other malignant phenotypes that accompanied the inflammatory immune microenvironment, so we performed KEGG analysis on the risk model gene and found that it may be closely related to the EMT phenotype. This hypothesis was confirmed because EMT-related gene sets were significantly enriched in the highrisk group. Subsequently, we found that the EMT phenotype was associated with a decrease in CD8 + T cells and an increase in M2 macrophages, which is different from the risk model. By analyzing the main signaling pathways that induce the EMT phenotype, we found that TGF-β was the most important one in gliomas, and TGFB1/TGFBR1 showed stronger immunosuppressive properties than PD-L1 and CTLA-4, especially in inducing an increase in CD8 + T cytopenia and M2 macrophages. It is clinically instructive that simultaneous low expression of TGFBR1 and PD-L1 has obvious survival advantages over other expression modes. Through single-cell analysis, we also found that TGFB1 is closely related to microglia and macrophages, especially M2 macrophages, which can explain why the increase in TGFB1/TGFBR1 expression is accompanied by a significant increase in M2 macrophages. Finally, we discussed the sensitivity of TGFB1 inhibitors in gliomas using cell line susceptibility data. From these analyses, we demonstrated a viable clinical strategy in combination with TGFB1/TGFBR1 inhibitors and PD-1/PD-L1 inhibitors for the treatment of high-risk gliomas.

Data availability statement
The original contributions presented in the study are included in the article/Supplementary Material; further inquiries can be directed to the corresponding author.