RNA Binding Protein-Mediated Immunosuppressive Tumor Microenvironment Contributes to Poor Prognosis in High-Grade Gliomas

The abnormal expression of RNA binding proteins (RBPs) in tumors can regulate the functions of immune genes and affect immune microenvironment. The characterization of immune inltration and its regulatory mechanism from the perspective of RBPs and RBP-regulated immune genes in high-grade gliomas were comprehensively studied by bioinformatics. Prognosis-related RBPs and associated immune genes were obtained from the cancer genome atlas database (TCGA). By clustering RBPs, the potential relationship between the different clusters of RBPs expression and immune inltration was explored. Immune-related RBP constructs immunosuppressive microenvironment in glioma by regulating immune genes and immune-related genes. Immunosuppression was the predominant type of phenotype with poor prognosis, while immunoinammation was the predominant type of immunoinammatory phenotype with good prognosis. RBP can construct immunosuppressive microenvironment by regulating immunosuppressive cells and stromal activation-related factors through directly and indirectly ways. Stromal activation plays a bridging role in tumor immunosuppressive regulation. The immunophenotyp depends on different enrichment states of immune-related RBP, which can change dynamically with the change of tumor pathological grade and affect the prognosis of patients.Our study suggests glioma related RBPs can potentially establish an immunosuppressive microenvironment by regulating immune genes and immune-related genes, which can be directly affected by RBPs clusters with distinct prognosis. immunosuppressive immunosuppressive regulating stromal activation. The association between stromal activation markers and immunosuppressive cell markers suggests that stromal activation plays a brideing role in mediating immunosuppressive microenvironment in high-grade gliomas.


Introduction
Gliomas are commonly prevalent malignancies found in the central nervous system (CNS), accounting for 24% of all primary CNS tumors and 81% of primary malignant CNS tumors (1). In 2007, the WHO classi cation of CNS tumors divided gliomas into grades I-IV, with grade III and IV being under high grade gliomas. The ndings have also shown that patient prognosis was associated with clinical grade, where the higher grade predicted worse prognosis (2). Although glioma can be treated postoperatively with adjuvant chemoradiotherapy, the 5-year survival rate of patients with high-grade glioma remains less than 5% due to dose limitations associated with radiation therapy and radiotherapy resistance (3).
Therefore, an urgent need remains to seek new effective biomarkers or treatments modalities to improve outcomes. Recent advances in immunotherapy against malignancies have brought renewed hope for cancer patients (4,5).
Tumor immunotherapy is a novel therapeutic modality that primarily exploits the body's immune system to control and kill tumor cells by enhancing or restoring antitumor immunity. Moreover, compared with conventional chemoradiotherapy, cancer immunotherapy has demonstrated the advantages of less damage to healthy organs, lower probability of drug resistance, and effective clearance of the residual tumor cells (6). The combination of cancer immunotherapy with conventional surgery and chemoradiotherapy has been observed to signi cantly improve the overall survival of patients with malignant tumors (7,8). However, immunotherapy has been found to be relatively ineffective in a larger number of patients in the clinic. This has been attributed to the ndings that tumors are able to establish an immunosuppressive environment and thereby escape from the tumor elimination actions performed by host immune system (9). Tumor cells can express various inhibitory molecules and secrete tumor associated factors to create an immunosuppressive microenvironment, which can effectively mediate immune tolerance towards tumors, but the regulatory mechanisms are not clearly understood (10).
RNA binding proteins (RBPs) are proteins that can bind to a wide range of RNAs, accounting for 7.5% of the genes encoding proteins (11). These RNA binding proteins can determine the function of each transcript by regulating mRNA stability, RNA localization, processing, splicing, degradation, and other related functions (12). A number of previous studies have con rmed that RBPs play critical roles in the tumor development, progression, and tumor immunity across various human tumors (13,14). The regulation of immune microenvironment by RBPs in glioma has been previously reported (15). However, there are few studies on the validation of this regulation mode and its regulatory mechanism from the perspective of RBPs and RBP-regulated immune genes.
In this study, transcriptomic and clinical data of a high-grade glioma cohort was extracted from the TCGA database, and the prognosis-related RBPs and related immune genes were also obtained. By clustering RBPs and immune genes related to RBPs, the potential relationship between the different expression clusters of RBPs, RBPs-related immune genes and immune cell in ltration in gliomas was investigated, and the possible causes for the presence of different expression clusters and the mechanisms leading to immune alterations were analyzed. Meanwhile, a risk model was established by screening ve different genes in the prognostic RBPs of high-grade glioma for the clinical prediction. We comprehensively studied the characterization of immune in ltration and its regulatory mechanism from the perspective of RBPs and RBP-regulated immune genes by bioinformatics. The ow diagram has been displayed in

Data
The data was obtained from The Cancer Genome Atlas Database (GDC; https://portal.gdc.cancer.gov/).
The data consisted of ve normal brain samples and 698 glioma (grade II-IV) samples (downloaded on May 25, 2021). In addition, two datasets (325 and 693, containing data and clinicopathological information) were downloaded from the Chinese Glioma Genome Atlas (CGGA) website (http://www.cgga.org.cn/) on April 3, 2021. However, for data downloaded from both the databases, only expression data and clinicopathological information of patients with grade III-IV gliomas were extracted. The genes for RBPs were obtained from the Eukaryotic RNA binding proteins (EuRBPDB) website (https://EuRBPDB syshospital.org), which included a total of 2949 RBPs (see Supplementary Table S1). 2483 immune genes were obtained from the Tumor Immune Estimation Resource (TIMER) Database (https://cistrome.shinyapps.io/timer/) (see Supplementary Table S2).

Exclusion criteria
The exclusion criteria used for this study was as follows: (a) glioma patients without overall survival (OS) information or with an OS<30 days (for reduced statistical bias); (b) patients lacking clinical information or expression data; (c) patients with gliomas graded as I-II. After initial data screening, a total of 403 samples (including 398 tumor and 5 normal samples) from the TCGA database, and 650 samples from the CGGA database were included in this study. The various clinicopathological and molecular characteristics of glioma patients have been shown in Table 1.

Data processing
The "limma" package in R language(16)was thereafter used to identify the levels of various differentially expressed RBPs between the normal brain and high-grade glioma tissues with the settings of log twofold change (log2FC)>0.5 and false discovery rate (FDR)<0.05. A total of 1289 RBPs were found to be differentially expressed in gliomas. Meanwhile, the expression data of 2483 immune genes were also extracted. The correlation analysis between RBPs and immune genes was carried out using "limma" package of R software. A total of 1120 immune related RBPs and 565 immune genes related to RBPs were obtained by setting following parameters: Correlation coe cient >0.5, P<0.001 (see Supplementary  Table S3). 308 immune genes were differentially expressed following the screening criteria of FDR < 0.05 and log2FC>0.5 in 565 immune genes.

Functional enrichment analysis
The biological functions of the various differentially expressed RBPs and related immune genes were observed using Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analyses. We used the "clusterpro ler" package (17) in the R language for enrichment of the various genes, and both P and FDR values less than 0.05 were considered as statistically signi cant.

Cluster analysis
Unsupervised clustering was applied to analyze the potential immune gene expression clusters associated with different RBPs using the "ConsensuClusterPlus" package with a 1000-fold replicate to ensure the stability of this clustering(18).
2.6 Gene Set Variation Analysis (GSVA) and functional annotation GSVA enrichment analysis was performed using the "GSVA" R package to observe the various differences in the biological processes in RBPs modi cation cluster (19). Thereafter, the gene set of "C2.cp.kegg.v6.2.symbols" downloaded from the Molecular Signatures Database (MSigDB) was used to conduct the GSVA analysis. An adjusted P<0.05 was considered as statistically signi cant. The "ClusterPro ler" R package was then used to functionally annotate the RBPs-associated genes with a cutoff value of FDR <0.05.

Estimation of cell in ltration in TME
CIBERSORT was then used to examine the relative fractions of the 22 in ltrating immune cell types in each tumor sample and to determine the variations among immune cells present in the different clusters (20). The codes were obtained from CIBERSORT website (https://cibersort.stanford.edu).
Thereafter, the immune score and stromal score of the tumor samples were calculated using the "estimate" R package.

Construction of prognostic model
Univariate Cox regression analysis was performed using the "survival" R package to assess the possible correlation between patients' OS and RBPs (P-value<0.001). Thereafter, signi cant candidate genes were then further screened using Lasso Cox regression. The different signi cant candidate genes were then screened using multivariate Cox analysis, and a prognostic risk score model for RBPs was constructed based on the signi cant candidate genes screened above, and a speci c risk score was then calculated to evaluate patient outcome. The risk score for each sample was de ned as follows: Risk score = β1* Exp1 + β2* Exp2 + βi * Expi where β represents coe cient, Exp represents expression level. The different patients with glioma were then divided into low-and high-risk groups according to the median value of the risk score. The survival differences between the two different groups were compared using the Kaplan-Meier method. Additionally, "Survival" receiver operating characteristic (ROC) package was used to validate the model, and to calculate the area under curve (AUC) to evaluate the e cacy of the immune gene prognostic model (21). In addition, the samples from the CGGA dataset were used to accurately assess the predictive power of this prognostic model. Finally, nomogram with calibration plots was carried out to predict the likelihood of OS using the "RMS" R package.

Validation of expression levels and prognostic signi cance
We extracted the data from the human protein atlas (HPA) online database (http://www.proteinatlas.org/) for detecting the expression of ve different RBPs at the translational level to validate the expression of these RBPs in gliomas(22).

Clinical sample testing
We collected 3 clinical samples of high-grade glioma for mRNA and protein detection of RBP, the samples included a 35-year-old male (grade III glioma, Patient 1), a 50-year-old female (grade IV glioma, Patient 2), and a 37-year-old male (grade IV glioma, Patient 3). The total protein extracted from tissue or cells were resolved on 10% SDS-PAGE gel (BioRad, USA), The proteins were therafter electrophoretically transferred to PVDF membranes (BioRad, USA). The blot was incubated with 5% skimmed milk and incubated with primary antibody at 4°C. The membranes were then incubated with a secondary anti-rabbit antibody for 2h, and the speci c protein bands were detected using the ChampChemi imaging system (BioRad, USA), and the protein expression level was normalized to GAPDH. The main antibodies used in the study: Anti-CLICL1 (BioSS, China), anti-CTIF (Thermo sher, USA), anti-EIF3L (BioSS, China), anti-PXN (BioSS, China), anti-TOP2A (PTGCN, China).

Statistical analysis
Kaplan-Meier curves and two-sided rank tests were employed to assess whether there was a signi cant difference in overall survival between the various populations. The categorical data was analyzed using χ2 or Fisher's exact test. The t-test or Wilcoxon test were then used for continuous data. For all statistical analyses, P<0.05 was considered as statistically signi cant. All statistical analyses in this study were done based on the R programming language (version 4.0.2; https://www.r-project.org/).

Ethics approval and consent to participate
All experimental protocols for this study was approved by the Institutional Ethics Committee of the

Possible impact of the different expression clusters of RBPs on glioma immune in ltration
Through univariate regression analysis of 1120 Immune-related RBPs screened from TCGA database, a total of 749 prognostic RBPs were obtained. Considering the impact of different expression clusters of RBP on its associated genes, we performed unpredicted clustering of expression clusters of 749 differential RBPs using the "ConsensusClusterPlus" package, which nally led to the identi cation of two different modi cation clusters (see Supplementary Fig. S1A). Cluster A comprised of 187 patients and Cluster B comprised of 211 patients ( Figure 2A). Indeed, Cluster A was found to enrich more patients whose tumors were of grade IV exhibited signi cantly higher mortality compared with Cluster B. On the contrary, Cluster B enriched greater number of patients with tumors classi ed as grade III, and less frequent mortality. Kaplan-Meier curves indicated that Cluster A conferred a worse prognosis among the two ( Figure 2E).
To explore the potential impact of different expression clusters of RBPs on glioma immunity, we did GSVA enrichment analysis. The analysis results showed that the various immune related pathways such as that of antigen processing and presentation, natural killer cell-mediated cytotoxicity, immune de ciency, allograft rejection, graft versus host were found to be enriched in the Cluster A. However, Cluster B was enriched for wingless-type MMTV integration site family (WNT) signaling, Hedgehog (Hh) signaling, calcium signaling, and phosphatidylinositol signaling ( Figure 2B). We hypothesized that the different pathways and regulators involved in immunosuppression and immune rejection were enriched in glioma by Cluster A leading to poor prognosis of patients.
Thereafter, heatmap of immune in ltration by CIBERSORT showed a signi cantly elevated expression of immunosuppressive related cell include regulated T cells (Tregs), M2 macrophages, naive CD4 T cells, non-activated NK cells, and Cells with immune function such as CD8 T cells, activated CD4 T memory cells, gamma delta cells, M0, M1 macrophages and neutrophils, stromal score, and immune score that could be associated with immunosuppression in Cluster A. This cluster was consistent with the previously reported immune phenotype characterized by immunosuppression and immune rejection 23 . However, the immune cells that in ltrated in Cluster B included activated NK cells, monocytes, activated dendritic cells, activated mast cells, eosinophils, stromal score, and exhibited a markedly decreased immune score. This cluster showed an innate immune in ltration ( Figure 2C). This was consistent with the previously reported immune in amed phenotype and was primarily characterized by adaptive immune cell in ltration and immune activation (23). These results suggest that immunoinvasive phenotypes of immunosuppression and immunorejection are present in the RBP clustering pattern with poor prognosis. In addition, immune related pathway enhancement was also accompanied by an increased expression of immune checkpoint molecules PD-L1 ( Figure 2F). This clustering cluster could effectively distinguish the two samples with different prognostic effects on patients ( Figure 2G).
To verify the potential impact of RBPs expression clusters on glioma immune in ltration, we extracted 308 different immune genes associated with 749 RBPs as described above (see Supplementary Table  S4). A total of 244 immune genes were found to be associated with prognosis in subsequent univariate regression analyses. Among these, 66 genes were observed to be protective factors with HR <1, while 178 genes were classi ed as risk factors with HR >1 (see Supplementary Table S5). Moreover, two distinct clusters of modi cation were identi ed by clustering on 244 genes (see Supplementary Fig. S1B). Cluster A included 213 patients, and Cluster B contained 185 patients. The 178 risk factors were thereafter enriched in Cluster A, corresponding to patients with tumor grade IV and a higher proportion of deaths. On the contrary, Cluster B enriched 66 different protective factors, corresponding to patients with tumor grade III and a signi cantly lower mortality rate ( Figure 2D). Kaplan-Meier curves indicated that Cluster A conferred a worse prognosis, and there was a better survival advantage observed for Cluster B ( Figure  2H).  Figure 2I). However, the pathways associated with highly expressed genes in Cluster B were predominantly enriched in T cell receptor signaling pathway, B cell receptor signaling pathway, Vascular endothelial growth factors (VEGFs) signaling pathway, mitogen-activated protein kinase (MAPK) signaling pathway and natural killer cell cytotoxic activity ( Figure 2J).
The subsequent heatmap of immune in ltration by CIBERSORT showed that Cluster A displayed signi cantly higher stromal score and immune score. In ltrated immunosuppressive related cell include naive CD4 T cells, non-activated NK cells, M2 macrophages, and Cells with immune function such as CD8 T cells, activated CD4 T memory cells, gamma delta T cells, M0, M1 macrophages and neutrophils exhibited an elevated expression in Cluster A. This cluster was consistent with the previously reported immune phenotype characterized by immunosuppression and immune rejection (23). However, Cluster B possessed lower stromal score and immune score. The plasma cells, activated NK cells, monocytes, activated mast cells and eosinophils were highly expressed in Cluster B ( Figure 3A), thus presenting with a phenotype characterized by adaptive immune cell in ltration and immune activation. This expression cluster of immune genes was almost consistent with that observed with RBPs.
To further clarify the correspondence between clusters of RBPs expression and those of relevant immune genes, we did enrichment heatmap of expression clusters of RBPs and immune gene expression clusters. It was noted from the correspondence of RBPs and immune gene expression clusters, immune gene expression Cluster A corresponded directly with the expression Cluster A of RBPs, which enriched more patients with tumor grade IV, and caused higher patient mortality. However, the Cluster B of immune gene expression corresponded to that of RBPs expression, which was primarily enriched in patients with tumor grade III and exhibited relatively low mortality ( Figure 3B). Additionally, the alluvial diagram also illustrated this interesting relationship (see Supplementary Fig. S1C). The relationship of protein interaction between RBPs and related immune genes revealed a complex regulatory relationship between them upon analysis ( Figure 3C). There was harmonious relationship between them (such as the link between RBP protective factor and immune gene protective factor, RBP risk factor and immune gene risk factor) and antagonistic relationships between them also (such as the link between RBP risk factors and RBP protective factors, and the link between immune gene risk factors and immune gene protective factors).

Regulation mechanism of different expression patterns of RBPS on immune invasion
To observe whether RBPs can activate in ammation-related factors to construct tumor microenvironment, we investigated the expression of various chemokines and cytokines characterizing three different gene clusters. Cytokines and chemokines selected from the published literature were extracted, among which TGRB1, Smad9, Twist1, CLDN3, TGFBR2, Acta2, COL4A1, Zeb1, and Vim were considered to associate with transcripts of the transforming growth factor (TGF) B/EMT pathway. PD-L1, CTLA-4, IDO1, LAG3, HAVCR2, PD-1, PD-L2, CD80, CD86, TIGIT, and TNFRSF9 have been reported to be involved in transcripts of immune checkpoints. TNF, IFNG, TBX2, GZMB, CD8A, PRF1, GZMA, CXCL9, and CXCL10 were associated with transcripts related to immune activation (24,25). The results clearly showed that these pathways that were involved primarily in immunosuppression were enriched in RBPs Cluster A and were differentially expressed compared with Cluster B (see Supplementary Fig. S1E). The diagram of protein interaction between RBPs and in ammatory pathway-related genes shows that glioma related RBPs can directly regulate in ammatory pathway-related genes ( Figure 3D). This result illustrated that RBPs could activate immune-related pathways to create tumor microenvironment in gliomas.
To verify the possible distribution pro le of immunosuppressive cells in immune gene expression clusters, we investigated the expression of different cell surface markers associated with immune regulation. Immune regulatory cells and their markers were obtained from the published literature, including regulatory T cells whose major surface markers include CD25, CCR7, FOXP3, CD4, CD127, CD152, and GITR; Regulatory B cells whose critical surface markers comprise CD19, CD20, CD24, CD27, and CD38; Tumor associated macrophages whose surface markers contain CD14, CD16, CD64, CD86, and CD163; and myeloid derived suppressor cells which can encompass several surface markers including that of CD11B, CD33, CD34, and HLA-DR(26). These factors were all signi cantly enriched in the expression Cluster A ( Figure 3G), and there were marked differences in the expression of these factors between clusters A and B ( Figure 3F), thereby illustrating that there was a distinct immunosuppressive phenomenon established in the Cluster A of immune genes.
It has been suggested that stromal activation mediates immunosuppression in tumors(23). To con rm this phenomenon, markers of stromal activation were extracted from previous studies, comprising 141 factors(27) (see Supplementary Table S6). According to the distribution of these factors in the immune gene cluster pattern, stromal activation-related factors were all enriched in the immune gene cluster A pattern, suggest that immunosuppression in high-grade gliomas is associated with tumor stromal activation ( Figure 3G). From the protein interaction between RBPs, immunosuppressive cell markers and stromal activation markers, we believe that RBP regulates both immunosuppressive genes and stromal activation-related factors to constitute the immunosuppressive microenvironment in gliomas. The protein interaction between immunosuppressive cell markers and stromal activation markers was stronger than the regulation of immunosuppressive genes by RBPs, that implies the bridging role of stromal activation in tumor immunosuppressive regulation ( Figure 3E).

The speci c causes for the regularity of RBP clustering pattern
In order to explore the speci c reasons of RBP clustering rules, the lasso regression analysis was applied on 749 genes associated with the prognosis (Figure 3I-J), which resulted in 32 factors with a signi cant prognostic impact. Among them, 14 factors were found to be protective with HR <1 and 18 were risk factors with HR >1 ( Figure 3K). Based on the degree of enrichment of the 32 factors in both the clusters, the risk factors were enriched in Cluster A, while the protective factors were enriched in Cluster B and depicted a low risk ( Figure 3H). thereby indicating that the enrichment status of protective and risk factors could potentially determine the expression of clusters and the survival prognosis of patients.
In addition, we further investigated that whether the status of the enrichment of risk and protective factors could signi cantly change as a result of changes in tumor grade. We included grade II-III glioma samples to contrast with grade III-IV glioma samples. The results showed that the factors of RBPs with prognostic impact were markedly decreased in grade II-III glioma samples compared with grade III-IV glioma samples, and the factors of RBPs were only found to be partially the same between the two groups (see Supplementary Fig. S1F). Among the 32 factors that were found to be signi cantly prognostic in grade III-IV gliomas, only 24 of them remained prognostic in grade II-III glioma samples, and the remaining eight factors were found to be statistically insigni cant (see Supplementary Fig. S1G). There were no signi cant changes in the conversion of the various protective factors to risk factors or vice versa for any of the 24 factors with consistent effects in grade II-III, and in grade III-IV samples (see Supplementary Fig. S1H). However, for the enrichment of 24 factors in the RBPs clustering of grade II-III gliomas (see Supplementary Fig. S1D), the risk factors were further enriched in the Cluster B of poor OS prognosis, which directly corresponded to the samples with glioma grade III. However, the different protective factors were enriched in Cluster A with a good prognosis, which was consistent to the samples with glioma grade II (see Supplementary Fig. S1H). The results clearly demonstrated that the status of the enrichment of risk factors and protective factors could signi cantly change with the tumor grade, and the protective factors could always be enriched in the tumor grade with better prognosis.

Prognostic model
We next constructed a novel prognostic model and to validate it, we proceeded to download the relevant glioma transcriptomic data from CGGA and extracted only a total of 650 expression matrices of the expression data of grades 3-4 with complete clinical information. After extracting the RDB details from transcriptome data, the differential genes were obtained by setting the parameters such that P-value was below 0.01 in both KM and COX methods, and the difference in 5-year survival between the high and low risk groups was > 20%. A total of 395 differential genes were screened and a total of 390 prognosis related RBPs were identi ed by univariate regression analysis. A total of 190 genes were thus identi ed by intersecting the prognostic genes in the TCGA and those in the CGGA databases respectively. We screened 17 different genes with the best potential prognostic signi cance among 190 genes using the lasso regression algorithm. Among them, ARPP21, CHD3, CTIF, SNRPN, and EIF3L acted as protective factors with HR < 1, while ASPM, BRCA1, CLIC1, DDX60L, ERI1, KIF4A, MBD2, MLEC, NCAPD2, PXN, REXO2, and TOP2A were found to be risk factors with HR >1 ( Figure 4A). The expression pro les of these 17 RBPs in the tumors have been shown in Figure 4B. The forest plots of the prognostic impact of these 17 factors in the CGGA database and their expression pro les have been shown in Supplementary Fig.  S1I-J.
We further analyzed the relationship between the expression levels of 17 factors in the prognostic model and immune subtypes. Among the risk factors, REXO2, MBD2, ERI1, KIF4A and TOP2A were all highly expressed in the C5 immune subtype, which corresponded to the immune silent type, and was characterized by M2-type macrophages as the dominant in ltrating cells. CLIC1 was highly expressed in the C4 immune subtype, which corresponded to the lymphocyte depletion type. DDX60L was highly expressed in the C6 immune subtype, which was classi ed as TGF-B dominant type. Both C4 and C5 immune subtypes are common in gliomas. Among the protective factors, ARPP21 and CTIF were highly expressed in C1 immune subtype, which corresponded to tissue healing type and had Th2-biased acquired immune invasion. CHD3 is highly expressed in the C2 immune subtype, which corresponds to IFN-γ dominant type ( Figure 4C). The rest of the factors including ASPM, BRCA1, MLEC, NCAPD2, PXN, EIF3L and SNRPN in immune express no obvious difference between different subtypes. The immune subtypes belonging to these 17 factors were consistent with the immune types corresponding to our previous RBP clustering model.
We identi ed ve distinct prognostic factors by using multiple stepwise Cox regression analysis. CTIF and EIF3L were protective factors with HR <1 in the TCGA ( Figure 4D) and CGGA ( Figure 4E) databases, while CLIC1, PXN, and TOP2A were considered as potential risk factors with HR >1. These were consistent with their attributes in the model. Kaplan-Meier survival curves in the TCGA and CGGA datasets further indicated that high expressions of CTIF and EIF3L were favorable for patient survival, while high expressions of CLIC1, PXN, and TOP2A resulted in poor outcomes ( Figure 5J-K). We calculated the risk score for each patient according to the following formula: Risk score=(0.0060*ExpCLIC1) + (-0.0679* Exp CTIF) + (-0.0139*Exp EIF3L) + (0.0381* Exp PXN) + (0.0192* Exp TOP2A) After dividing the patients into different high-and low-risk groups according to the median value of the risk score in the cohort, the Kaplan-Meier survival curve indicated that patients in the high-risk group displayed a worse prognosis ( Figure 4F). In the analysis of the relationship between the high-low risk group and the immune subtype, we found that the high-risk group in the line-up was highly expressed in the C5 type of the immune subtype, while the low-risk group had no obvious enrichment characteristics in the C1-C6 immune subtype ( Figure 5D). In addition, the high-risk group had higher immune and matrix scores ( Figure 5F-G). This result also con rms our previous inference that RBP can construct an immunosuppressive microenvironment through matrix activation, leading to poor prognosis of patients. The AUC of ROC curve in the RBP score system was observed to be 0.893 ( Figure 4G), thereby indicating a nice diagnostic performance. In addition, the risk curves and gene expression heatmaps of the high-risk and low-risk groups were clearly distinguished by the risk scores derived from the ve RBPs as shown in Figure 4H.
To further validate the predictive ability of prognostic mode for clinical predictive power, we used the CGGA dataset. The risk scores for patients in the CGGA cohort were calculated using the same formula, and patients were then segregated into high-and low-risk groups based on the obtained median value. The results showed that glioma patients in the high-risk group displayed a worse OS and survival rates compared with those in the low-risk group ( Figure 5A). The AUC of ROC curve was 0.771 ( Figure 5B), thus clearly illustrating that this RBP signature possessed a stable predictive ability in glioma patients.
Furthermore, the risk curves and gene expression heatmaps of the high-risk and low-risk groups were markedly distinguished based on the risk scores constructed using the ve RBPs as shown in Figure 5C.
To assess the independent prognostic value of the RBP signature, we next performed univariate and multivariate Cox analyses in the TCGA cohort. TCGA dataset results indicated that the RBP risk strati cation method acted not only an OS-related prognostic factor in glioma, which was validated as an independent prognostic factor for glioma patients in each independent and combined cohort ( Figure 5H-I). Among them, ndings of univariate analysis indicated that HR was 1.2949 (1.2380-1.3544), P= 1.86E-29, and multivariate analysis showed that HR was 1.1564 (1.0754-1.2434), P= 8.71E-05. The similar results were obtained during the validation of the CGGA dataset (see Supplementary Fig. S2A-B). A nomogram model was built for clinically predicting the prognosis of high-grade LGG patients using the 5signature RBPs ( Figure 5E).

Clinical validation and drug response characteristics of RBP-related prognostic models
Immunohistochemical ndings for CLIC1, CTIF, EIF3L, PXN, and TOP2A were obtained from the Human Protein Atlas database (https://www.proteinatlas.org/). Moreover, compared with the normal brain glial cells, CLIC1 and EIF3L exhibited slightly elevated expression, whereas PXN and TOP2A showed strong high expression, and CTIF showed decreased expression in high-grade glioma ( Figure 6A). Additionally, clinical high-grade glioma samples (one was grade III and two were grade IV) further con rmed the mRNA (see Supplementary Fig. S2C-G) and protein expression level of the above factors in high-grade glioma ( Figure 6B-F). The mRNA and protein expression level of CLIC1, PXN and TOP2A showed high expression, EIF3L showed decreased expression in high-grade glioma. CTIF only detected high expression of mRNA in tumor, but not protein expression was detected. Finally, we predicted the drug sensitivity of 5 factors in the line-up, PXN was more sensitive to Dimethylaminoparthenolide, Carmustine, Oxaliplatin, Imexon and Ifosfamide (COR <-0.45, P <0.001) ( Figure 6G-K), whereas EIF3L was more sensitive to Chelerythrine (COR =0.481, P <0.001) ( Figure 6L). No signi cant correlation was found between CLIC1,CTIF and TOP2A and the included drugs (COR >0.45 or COR <-0.45).

Discussion
Tumor tissue is a complex entity composed of several tumor cells, stromal cells, immune and in ammatory cells as well as other non-cellular components. In this microenvironment (tumor microenvironment, TME) in order to facilitate tumor proliferation, multiple cytokines secreted by tumor cells, stromal cells and immune cells display complicated relationships such as those of mutual cooperation, inhibition, and antagonism with each other, and further constitute a cytokine network (28, 29). Tumors can exploit cytokine networks to exert growth promoting effects through a variety of different mechanisms including that of proin ammation, immunoediting as well as immune escape (30).
In this study, we investigated the immune regulation model and mechanism of immune-related RBPs in high-grade gliomas. By clustering 749 RBPs factors affecting prognosis were divided into A and B clustering modes, we found that patients with higher grade of tumor in mode A had poorer prognosis, while patients with lower grade of tumor in mode B had better prognosis. Correspond to the two different RBPs clustering, two expression clusters possessed distinctly different TME cell in ltration characteristics were generated. Among them, Cluster A showed a high stromal score and immune score, and the in ltrated immune cells including Tregs that can contribute to immunosuppression, M2 macrophages, non-functional naive B cells and naive CD4 T cells, non-activated NK cells, but also functional cells such as CD8-T cells, gamma delta T cells, M0 and M1 macrophages. This cluster embodied a distinct phenotype which was characterized by immunosuppression and mixed immune rejection. On the contrary, the Cluster B re ected the immune in amed phenotype, and the tumor present contained a higher abundance of immune cells, such as activated NK cells, monocytes, activated dendritic cells, activated mast cells, eosinophils, while stromal score and immune score were relatively signi cantly decreased compared with Cluster A. Both expression clusters were identical to the phenotypes that have been previously reported (15,23).
Moreover, adopting a different strategy from the previous studies, we used a dichotomous approach for RBPs clustering because the gliomas included in this study only contained those of grades III-IV.
Therefore, there are only two types of immune invasion phenotypes corresponding to the RBPs clustering.
While Cluster A simultaneously exhibited features common to both immunode ciency and immune rejection, we hypothesized that immunosuppression played a more dominant role in there. We con rmed this hypothesis from the perspective of immune genes. Compared with other studies with 3 or more clusters, we speculated that with increasing clusters, a signi cantly higher number of expression clusters that were predominantly associated with tumor immune disorders might appear in the expression clusters distributed in the intermediate clusters. The expression cluster with the best prognosis may then only express the type that can effectively activate a robust autoimmune response, whereas the expression cluster with the worst prognosis may be immunosuppressive in nature, as shown by Zhang (23).
We also conducted the mechanistic studies to gain insight into the regulation of immune responses by RBPs, and found that the protein interaction between RBP and related immune genes re ects a complex network regulatory relationship. There are harmonious relationship between them (such as the link between RBP protective factor and immune gene protective factor, RBP risk factor and immune gene risk factor) and antagonistic relationships between them also (such as the link between RBP risk factors and RBP protective factors, and the link between immune gene risk factors and immune gene protective factors). Furthermore, RBPS can regulate immunosuppression by directly and indirectly pattern. Risk factors of RBPs can directly regulate the immunosuppressive cell markers and stromal activation markers. There were also correlations between immunosuppressive cell markers and stromal activation markers, and the number of such correlations was signi cantly greater than the correlations between RBPs and immunosuppressive cell markers. These results indicated that RBPs could not only directly regulate immunosuppressive cells, but also indirectly regulate immunosuppressive cells by regulating stromal activation. The association between stromal activation markers and immunosuppressive cell markers suggests that stromal activation plays a brideing role in mediating immunosuppressive microenvironment in high-grade gliomas.
We also explored the pattern for the different expression clusters of RBPs. We performed an enrichment analysis by acquiring 32 different factors representative in 749 RBPs with signi cant prognostic impact. It was found that the different expression clusters of the risk and protective factors can effectively determine the overall prognosis. This conclusion was further corroborated by our enrichment using 244 RBPs associated immune genes. In addition, we also found the status of enrichment of risk and protective factors was observed to vary with tumor grade, and protective factors were always found to be enriched in tumor grade with better prognosis when more grade glioma samples were added.
Finally, we selected few hub RBPs based on univariate Cox regression analysis, Lasso analysis and multivariate Cox regression analysis. A total of ve different RBPs including CLIC1, CTIF, EIF3L, PXN, and TOP2A were identi ed as prognostically relevant hub RBPs. Among them, EIF3L and CTIF served as protective factors with hazard ratios less than 1, while CLIC1, PXN and TOP2A were found to be risk factors with hazard ratios greater than 1. The relationship between model risk score and immune subtype, immune score and matrix score further con rmed our previous inference that RBP could construct immunosuppressive microenvironment through matrix activation, leading to poor prognosis of patients.
EIF3L is an important member of eukaryotic initiation factor 3 (eIF3) family, and it plays a central role in the initiation of eukaryotic translation. A number of previous studies have con rmed EIF3L acts as a protective factor in glioma. EIF3L has been reported to be more highly expressed in low-grade gliomas relative to high-grade gliomas and can be associated with better overall survival (31). CTIF is a pivotal factor involved in the rst round of translation driven by the nuclear cap binding protein complex, which can primarily contribute to the quality control during protein and mRNA synthesis. The role of this factor in glioma remains poorly studied. There have been prior studies on liver cancer(32) and gastric cancer(33) that have conclusively established a tumor promoting role of CLIC1 in tumors. CLIC1 has been also found to be highly expressed in glioblastoma and can lead to reduced patient survival, and CLIC1 silencing can signi cantly decrease proliferation, clonogenicity and tumorigenic capacity of glioblastoma stem cells (34). PXN has been reported to be closely associated with the occurrence and development of a variety of tumors, such as liver cancer (35), cervical carcinoma(36), and gastric cancer (37). PXN was found to substantially enhance the tumor cell migration and invasion in glioblastoma(38). TOP2A is a critical enzyme that controls and alters the topological state of DNA during transcription, and it has been found to play a pro-tumorigenic role in glioblastoma (39), liver cancer (40), bladder cancer(41), as well as prostate cancer(42).

Conclusions
Glioma related RBPs can potentially establish an immunosuppressive microenvironment by regulating immune genes and immune-related genes, which can be directly affected by RBPs clusters with distinct prognosis.

Tables
Due to technical limitations, Table 1 is only available as a download in the Supplemental Files section. Figure 1 Flow chart of this study