Integrative Analyses and Verification of the Expression and Prognostic Significance for RCN1 in Glioblastoma Multiforme

Glioblastoma multiform is a lethal primary brain tumor derived from astrocytic, with a poor prognosis in adults. Reticulocalbin-1 (RCN1) is a calcium-binding protein, dysregulation of which contributes to tumorigenesis and progression in various cancers. The present study aimed to identify the impact of RCN1 on the outcomes of patients with Glioblastoma multiforme (GBM). The study applied two public databases to require RNA sequencing data of Glioblastoma multiform samples with clinical data for the construction of a training set and a validation set, respectively. We used bioinformatic analyses to determine that RCN1 could be an independent factor for the overall survival of Glioblastoma multiform patients. In the training set, the study constructed a predictive prognostic model based on the combination of RCN1 with various clinical parameters for overall survival at 0.5-, 1.0-, and 1.5-years, as well as developed a nomogram, which was further validated by validation set. Pathways analyses indicated that RCN1 was involved in KEAS and MYC pathways and apoptosis. In vitro experiments indicated that RCN1 promoted cell invasion of Glioblastoma multiform cells. These results illustrated the prognostic role of RCN1 for overall survival in Glioblastoma multiform patients, indicated the promotion of RCN1 in cell invasion, and suggested the probability of RCN1 as a potential targeted molecule for treatment in Glioblastoma multiform.


INTRODUCTION
Glioblastoma multiforme (GBM) represents the most prevalent brain cancer in adults and has a dismal prognosis and poor quality of life (Omuro and DeAngelis, 2013). The current treatment strategies for GBM are maximum surgical resection followed by a combination of chemotherapy and radiotherapy (Gilbert et al., 2013;Ostrom et al., 2015). Even with the advancement in therapeutic options over recent decades, recent studies have demonstrated that the median survival of GBM patients is 16.6 months, which decreases after 2 years, with a survival rate of only 34% (Gilbert et al., 2013). Several studies have illustrated that some omics markers within tumors could impact patients' survival, like the status of Isocitrate dehydrogenase 1/2 (IDH1/2) mutation, glioma-CpG island methylator phenotype (G-CIMP), methylation of O-6methylguanine-DNA methyltransferase (MGMT), and codeletion for chromosome 1p and 19q (1p/19q codeletion) (Hartmann et al., 2010;Wick et al., 2013;Hainfellner et al., 2014;Louis et al., 2016). In addition, further improvements have been made with subspecialized care, improved resection methods precisely targeted radiotherapy, and early systemic salvage therapies (Jayamanne et al., 2018). However, patients with GBM still have a poor prognosis due to the GBM's aggressive behavior, rapid progression, and frequent recurrence (Soffietti et al., 2014). Thus, it is imperative to search for a novel biomarker with good prediction for the prognostic signature of GBM via various methods, to explore the molecular mechanisms for precisely targeted treatments in GBM.
Reticulocalbin-1 (RCN1), a calcium-binding protein, contains six conserved regions and is located in the endoplasmic reticulum (Ozawa and Muramatsu, 1993), which regulates calciumdependent activities combined with reticulocalbin 2 (RCN2) (Nakakido et al., 2016). Dysregulation of RCN1 protein has been reported in multifarious diseases, including cancer, cardiovascular, and neuromuscular diseases (Grzeskowiak et al., 2003;Liu et al., 1997;Zhang et al., 2006). It has been found that RCN1 is involved in breast cancer (Nakakido et al., 2016), colorectal cancer (Nimmrich et al., 2000), liver cancer , kidney cancer (Giribaldi et al., 2013), and non-small cell lung cancer . In addition, it is reported that down-regulation of RCN1 facilitates apoptosis and necroptosis in prostate cancer cells (Liu et al., 2018). Overall, the above findings have revealed that dysregulation of RCN1 could contribute to tumorigenesis and progression. However, the relationship between RCN1 and prognosis, nor its biological functions in GBM, have been completely investigated.
In our study, GBM samples in the Cancer Genome Atlas (TCGA) database were enrolled in the training set, and cases in the Chinese Glioma Genome Atlas (CGGA) were used for the external validation set, to assess the RCN1-related prognostic signature in GBM. Gene Expression Profiling and Interactive Analyses (GEPIA2) (Tang et al., 2017) (http://gepia2.cancer-pku. cn/#index) were used to profile the tissue-wise expression of RCN1 in GBM. In addition, the next version of Cell type Identification By Estimating Relative Subsets Of RNA Transcripts (CIBERSORTx, (Newman et al., 2019) https:// cibersortx.stanford.edu/) was used to illustrate the abundances of the infiltration immune cells correlating with expression of RCN1 in GBM, and the single sample gene search enrichment analyses (ssGSEA) (Barbie et al., 2009) was utilized to determine potential pathways of RCN1 involved (Figure 1).

Data Acquisition
The RNA sequencing data (level 3) of GBM with corresponding clinical data were downloaded from the TCGA database to assess the prognostic impact of RCN1-signature on patients with GBM, and were enrolled into the training set. Then, GBM samples from the Chinese Glioma Genome Atlas (CGGA) dataset were employed as an external validation set for further verification.

Differential Analysis of Reticulocalbin-1 Expression in Glioblastoma Multiforme
Gene Expression Profiling and Interactive Analyses (GEPIA) (Tang et al., 2017) is an online resource for gene expression analysis, and GEPIA2 (http://gepia2.cancer-pku.cn/#index) (Tang et al., 2019) was an updated version of GEPIA, containing 163 GBM cases of TCGA and 207 normal brain samples of the Genome Tissues Expression database (GETx) (GTEx Consortium, 2015). We performed a differential analysis of RCN1 expression through GEPIA2 and used a boxplot to show the results with log2 of transcript count per million [log 2 (TPM + 1)] presenting the expression level of RCN1.

Survival Analysis
As described previously (Ogłuszka et al., 2019), we used survMisc package in R to divide GBM samples into high and low RCN1 groups based on the optimum cutoff. For the exploration of the association between the expression of RCN1 and survival, Kaplan-Meier (KM) analysis with log-rank test was conducted by using survival package.

Evaluation of Reticulocalbin-1-Based Prognostic Value
To explore the effect of the variation of RCN1 expression on OS in GBM patients, we estimated the relative risk and log2-based transformation for survival status, followed by fitting smooth line using a locally estimated scatterplot smoothing (LOESS) method (Zheng et al., 2020). We then compared the predictive effect of RCN1 with different prognostic factors by receiver operating characteristic (ROC) curve (Zhou et al., 2019). We regarded the RCN1 expression as a single continuous covariate for further regression analysis based on the result of LOESS analysis. We adjusted for clinical features, which changed the regression coefficients of RCN1 by more than 10%, or of p < 0.1 through the univariate analysis (Kernan et al., 2000). Next, adjust I model was determined after adjusting confounders, while all clinical factors were enrolled as adjust II model. Then, univariate and multivariate analyses using Cox proportional hazard regression model were performed. We also applied LOESS method to visually assess the relationship between RCN1 expression and OS in adjust I and adjust II models, respectively. Interaction test and stratified analysis (Soria et al., 2015) were also carried out to assess the differential prognostic value of RCN1 in accordance with a different model. A two-tailed p < 0.05 was considered to be statistically significant.

Construction and Comparison of Three Prognostic Models
After determining the prognostic features, five clinical prognostic factors in GBM (age, gender, IDH status, chemotherapy, as well as radiotherapy) were enrolled to construct the prognostic model as model 1. Model 2 (the expression of RCN1) and model 3 (model 1 + model 2) were compared with model 1 to estimate the robustness of different prognostic features to GBM prognosis, respectively. We used discrimination, calibration, and model improvement capability to evaluate the different models. The discrimination of these models was estimated by ROC curve, concordance index (C-index) (Harrell et al., 1996), as well as decision curve analysis (DCA) (Vickers and Elkin, 2006;Vickers et al., 2008;Rousson and Zumbrunn, 2011;Kerr et al., 2016). Three corresponding models were constructed with coxph function in survival package, and then the risk value of patients in each model was calculated with predict function, and thus the ROC curves were constructed by timeROC (Heagerty et al., 2000;Blanche et al., 2013) with those risk values. Three models were firstly built with cph function, and then their C-index at different times with the prediction error curve was calculated by cindex function in pec package. The stdca package was applied to visualization for DCA. The calibration curves of different models were completed via caoplot function in pec package. Notably, we used the bootstrap method with 1,000 resamples both for analyses of discrimination and calibration. Moreover, the improvement capability of the model was assessed through the net reclassification improvement (NRI) and the integrated discrimination improvement (IDI) by using IDI. INF.OUT function in survIDINRI package (Pencina et al., 2008). After the best model was determined, regplot package was employed to construct the diagram of the nomogram.

External Validation of the Prognostic Signature
The prognostic capability of the RCN1-based signature was externally validated in the validation set. The gene expression information with corresponding clinical factors of GBM was obtained from the CGGA database. Similarly, we divided GBM patients into high-and low-RCN1 groups by the optimum cut-off value, and then performed the KM survival curve. Afterward, using predict function, the time-independent ROC curves, C-index, and calibration curves were employed to assess the accuracy of prognostic model RCN1 signature-based. mixed cell population, by using gene expression data (Newman et al., 2019). We applied CIBERSORTx to analyze the infiltration level of 22 immune cells in high-and low-RCN1 groups, using samples from TCGA and CGGA datasets, respectively. After enabling batch correction, performing the "Bulk mode", and selecting the quantile normalization algorithm, the results were represented with an absolute score for the proportion of 22 immune cell subsets of GBM samples. Consecutively, the samples with p < 0.05 were retained after repeating the crossover operation 500 times (Ali et al., 2016). Wilcoxon rank-sum test was applied to identify the differences between the two groups.

ssGSEA
The ssGSEA was used to identify the differentially enriched hallmarks for a single sample (Barbie et al., 2009). To identify key pathways related to RCN1, we chose to focus on 50 hallmark gene sets, which were designed to highlight gene sets contained in the Molecular Signatures Database (MSigDB) (Subramanian et al., 2005). Gene symbol profiles for homo sapiens were downloaded from MSigDB database (Liberzon et al., 2015). Then we estimated the degree of each hallmark's ssGSEA profile in two groups, using the gsva package, both in the training and validation cohort. Next, by limma package, differential analysis was performed; and |t| > 1 or adjusted p < 0.05 were considered as statistically significant.

Cell Culture
Human GBM cell lines U87 and A172 were purchased from American Type Culture Collection (ATCC, Manassas, Virginia, United States), which were authenticated with a short tandem repeat. Cells were set on a humidified incubator with 5% CO 2 at 37°C, as well as cultured in Dulbecco's Modified Eagle's Medium (GIBCO, Billings, MT, United States) added with 10% fetal bovine serum (GIBCO).

Small Interfering RNA Transfection
The small interfering RNA of RCN1 (si-RCN1) sequences and the corresponding negative control were designed and purchased from RiboBio (Guangzhou, China). For transient silencing, A172 and U87 GBM cell lines were transfected with negative control or si-RCN1 by Lipofectamine TM 3,000 Reagent (Invitrogen, United States) according to the manufacturer's instruction. After 48 h, cells were harvested for quantitative real-time-polymerase chain reaction (qRT-PCR) analysis. Target sequences for transient silencing were as follows: si-RCN1-1: GAAGCTAACTAAAGA GGAA; si-RCN1-2: CCAGGCATCTGGTATATGA; negative control siRNA was obtained from RiboBio (Guangzhou, China).

qRT-PCR
By using ReverTra Ace ® qPCR RT Master Mix with gDNA Remover (TOYOBO, Shanghai, China), total RNA was extracted and then was used to synthesize the first complementary DNA (cDNA) strand according to the manufacturer's protocol. The qRT-PCR reaction was carried out to estimate the RNA levels and ACTIN was used as the internal reference. The primers used for qRT-PCR were as follows: RCN1 forward 5′-AAGGAGAGGCTAGGGAAGATT-3′ and reverse 5′-ATCCAGGTTTTCAGCTCCTCA-3'; ACTIN forward 5′-CACCATTGGCAATGAGCGGTTC-3′ and reverse 5′-AGG TCTTTGCGGATGTCCACGT -3'. The relative normalized expression of the target genes was compared with that of ACTIN, and the mRNA expression of each gene was calculated with the 2 −ΔΔCt method .

Cell Invasion Assays
Cell invasion was measured through wound healing and transwell migration assays following the manufacturers' instructions. In brief, cells were plated in 6-well plates and cultured at 37°C with 5% CO 2 . After 24 h, at which cells were reached on 80% confluence, we used a sterile 10 μL disposable serological pipette to make a straight-line scratch, and then cells were harvested after 48 h. Images of the scratch width were taken using an inverted microscope (Olympus IX73 Inverted Microscope, Olympus, Beijing China) at 0 and 48 h after the scratch, and then calculated by ImageJ software (version 1.52, National Institutes of Health, United States).
As for the transwell migration assay, it was performed using the Boyden chamber with a gelatin-coated polycarbonate filter with an 8-μm pore size (Neuro Probe, Gaithersburg, MD, United States). Cells were added to the upper chamber in 24well plates at a density of 5.0 × 10 4 cells per well, and the lower chamber was filled with 800 μL 10% FBS for 24 h. Cells in the transwell chamber were fixed with 4% paraformaldehyde for 15 min, stained with 0.1% crystal violet for 30 min, and then observed by a BDS500 Inverted Biological Microscope (Chongqing Optec Instrument Co., Ltd., China).

Statistical Analysis
All the data were presented as mean ± standard deviation and the statistical analyses were performed by IBM SPSS Statistics 25.0 software (International Business Machines Corporation, United States). The student's t-test was performed to evaluate the significant difference between the two groups.

Patients Characteristics
A total of 361 GBM samples (153 patients from TCGA as the training cohort and 208 patients from CGGA as the validation cohort) were obtained in our study, as shown in Table 1.

Reticulocalbin-1 Was Elevated in Glioblastoma Multiforme and May Act as an Oncogene
As shown in Figure 2A, we found that the expression of RCN1 was higher in the GBM samples (n 163) than in the normal brain tissues (n 207), by using the GEPIA2 tool. In the TCGA database, as the near-linear correlation between the variation of RCN1 expression and OS revealed through LOESS ( Figure 2C), RCN1 expression was considered as a single continuous variation for further analysis. A total of 153 samples were clustered into the high-(n 97) or low-RCN1 group (n 56) by the optimal cut-off value (4.144). Patients with higher RCN1 expression had a worse OS than those with a low one in GBM (p 0.001) ( Figure 2B).

Reticulocalbin-1 is an Independent Prognostic Signature for Glioblastoma Multiforme
Considering the interference of the confounding factors, identifying and then adjusting for potential confounding factors was conducted. We firstly found that RCN1 may be an independent prognostic signature when compared with other signatures (Supplementary Figure S1). In the training cohort, we then identified variants (age and radiotherapy) to be adjusted and then enrolled these clinical features into the adjusted I model. The adjustments for age, gender, IDH status, chemotherapy, and radiotherapy were included in the adjusted II model. Both non-adjusted and two adjusted models were analyzed by using the Cox regression analysis to further investigate whether RCN1 could estimate OS independently ( Figure 3). As shown in  Table 2). Furthermore, after adjusting for five predominant clinical features (age, gender, IDH status, chemotherapy, and radiotherapy), RCN1 independently predicted prognosis in the training cohort (HR 1.67, 95% CI 1.13∼2.45, p 0.009) ( Table 2). In the same way, we also found the near-linear correlation between the variation of RCN1 expression and OS both in adjust I model ( Figure 2D) and adjust II model ( Figure 2E), thus, we enrolled RCN1 as a single continuous variation for further analysis. In addition, subgroup analysis showed that there were no statistical differences neither in the non-adjusted model nor adjusted I model, except for chemotherapy (P interaction 0.0019 in the non-adjusted model, P interaction 0.0118 in the adjusted I model, respectively) ( Table 3), revealing that RCN1 might be an independent prognostic factor for OS in patients with GBM.

Construction and Evaluation of Three Prognostic Models
The clinical features and RCN1 were enrolled to construct the prognostic models for GBM. We firstly built three prognostic models (model 1: five clinical variants, model 2: the expression of RCN1, and model 3: model 1 + model 2) and then evaluated them. Model 3 had a higher area of under curve (AUC), better C-index, and lower prediction error compared with model 2 and model 1 (Figures 4A-C). DCA showed that the net benefit of  model 3 in 0.5 and 1 year is better than the other two models, but there was no significant difference in 1.5 years ( Figures 4D-F). It was found that the calibration of model 3 was better than that of model 1 and model 2 in 0.5 and 1 year, while the calibration of the three models was poor in 1.5 years (model 2 was better than others) ( Figures 4G-I). As for model improvement capability, when model 1 was considered as the reference, the NRI and IDI of model 3 were both positive, in which the NRI of 1.5 years was increased by 16.7% (p 0.064) and the IDI of 1.5 years was increased by 3.7% (p 0.034); on the contrary, the NRI and IDI of model 2 were both negative, where the IDI of 0.5 year was decreased by 11.4% (p 0.0016) meanwhile the NRI was decreased by 27.4% (p 0.040) ( Table 4). From the above results, it can be determined that model 3 had good discrimination and calibration in the prediction of OS. Therefore, we developed a nomogram in accordance with model 3 to assess OS at 0.5-, 1.0-, and 1.5-years in the TCGA dataset, in which each signature was assigned points according to its risk contribution to OS ( Figure 5).

External Validation for Nomogram in Chinese Glioma Genome Atlas
To assess whether an RCN1 signature-based model can similarly play a prognostic value in different populations, a total of 208 GBM samples from the CGGA database as an external validation cohort were used to assess its prediction performance. According to the optimum cutoff (5.521), high-(n 147) and low-(n 61) RCN1 groups were determined. Consistent with the findings in the training cohort, the Kaplan-Meier curve revealed patients with high RCN1 represented a worse OS than those with low expression (p 0.0047) ( Figure 6A). Moreover, the timedependent ROC curves were performed and the AUCs of the 0.5-, 1.0-, and 1.5-year survival for the constructed nomogram in the training cohort were 0.737, 0.673, and 0.694, respectively, in the validation cohort ( Figure 6B). The nomogram shared a C-index of more than 0.6 at different times ( Figure 6C) and a relatively low prediction error ( Figure 6D). Finally, the calibration curves for this nomogram in the validation cohort at 0.5 year was poor, while those at 1 year and 1.5 years were good ( Figures 6E-G).

Differential Abundances of Infiltrative Immune Cells
By using the CIBERSORTx algorithm, the relative proportions of 22 immune cells between two groups in GBM were obtained. In the TCGA dataset, there were three types of infiltrative immune cells with significant difference at different groups, whereas in CGGA there was fourteen, in which T cells CD4 memory resting, eosinophils, and macrophages M0 were differentially expressed in two data at the same time (Figure 7). With more detail as shown by bar plots in Figure 7A, in the training cohort, the infiltration level of eosinophils was significantly higher in the low-risk group, whereas the infiltration level of macrophages M0 and T cells CD4 memory resting was significantly higher in the high-risk group. In the validation cohort, the infiltration levels of B cells naive, dendritic cells resting, mast cells activated, macrophages M1, eosinophils, T cells CD4 memory resting, neutrophils, NK cells activated, and monocytes were significantly higher in the low-risk group, whereas the B cells memory, macrophages (M0, M2), T cells regulatory (Tregs), and plasma cells were significantly higher in the high-risk group ( Figure 7B).

Pathway Enrichment Analysis
The ssGSEA was used to identify signaling pathways RCN1involved in GBM, and then demonstrate significant differences in the enrichment of MSigDB hallmark gene set in the TCGA and   CGGA databases, respectively. The results indicated that there was no significant pathway screened in the TCGA database ( Figure 8A), whereas the KRAS-signaling-DN pathway was significantly involved in the low-RCN1 group and some pathways, including reactive oxygen species, MYC targets V2 and V1, apoptosis, and DNA repair pathways, in the high-RCN1 group in the CGGA database ( Figure 8B).

Knockdown of Reticulocalbin-1Inhibited Cell Invasion in Glioblastoma Multiforme
Finally, in order to elucidate the effect of RCN1on GBM cell invasion, we conducted a series of morphological and molecular biological experiments. The si-RCN1-1 and si-RCN1-2 could effectively reduce endogenous RCN1 mRNA expression by mRNA levels in both U87 and A172 cells (Figures 9A,B). Later, we confirmed that si-RCN1 could decrease cell invasion ( Figures 9C-F).

DISCUSSION
GBM is the most aggressive brain tumor. The prognosis of patients with GBM remains poor, although treatment strategies, including maximum surgical resection, radiation, and chemotherapy, have been conducted. The role and mechanism of biomarkers in GBM tumorigenesis are very important for the development of treatment of patients with GBM (Sasmita et al., 2018). Recently, ample molecular markers have been identified, which could provide new insight regarding GBM formation and progression. The discoveries and roles of molecular biomarkers, GBM-specific microRNAs, and GBM stem cells were delineated (Sasmita et al., 2018;Hassn Mesrati et al., 2020). There are some genetic mutation features, epigenetic modification, and some molecular alterations in situ. As highlighted by Alireza Mansouri and his colleagues, the methylation of MGMT promoter has been identified to provide better outcome prediction when GBM patients receive temozolomide chemotherapy (Mansouri et al., 2019). Meanwhile, the amplification of CCND2 suggests better clinical outcomes in IDH-mutant patients while the elevated total copy number variation, co-amplification of CDK4/MDM2, amplification of PDGFRA, and CDKN2A show worse association (Zheng et al., 2013;Sasmita et al., 2018;Mirchia et al., 2019;Louis et al., 2021).
Large-scale studies further demonstrated that IDH-mutant GBMs tend to have a higher concentration of 2-hydroxyglutarate (an oncometabolite produced as a result of the IDH1/2 mutation), which may promote infiltration via upregulation of HIF-1α and VEGF(K. Mirchia and Richardson, 2020;Sasmita et al., 2018). For IDH-wide type GBM patients, the mutation of BRAF may benefit them while the co-alteration of EGFR/PTEN/CDKN2A and mutation of PIK3CA and H3K27M correlates with worse clinical outcomes (Mirchia and Richardson, 2020;Umehara et al., 2019). Circulating biomarkers have also been developed in GBM because of their non-invasive potential (Jelski and Mroczko, 2021;Kefayat et al., 2021), but the sensitivity and specificity are still problems to be solved (Müller Bark et al., 2020;Raza et al., 2020;Jones et al., 2021). Additionally, there are It could not even be proven the release of potential biomarkers were exclusive from tumor cells or not. According to a systematic review, about 133 distinct biomarkers were identified from 1853 Frontiers in Molecular Biosciences | www.frontiersin.org October 2021 | Volume 8 | Article 736947 10 patients and evaluated based on level of evidence (IA∼IVD) using an adapted framework from the National Comprehensive Cancer Network guideline (Raza et al., 2020). Nevertheless, few can reach level IA and further refinements are needed. Therefore, it is urgent to identify new biomarkers that can illustrate the indepth molecular mechanism, enhance the diagnosis, respond to treatment, and provide prognostic prediction of GBM. Here, we identified the potential of RCN in GBM. RCN1, as a calciumbinding protein located in the lumen of the endoplasmic reticulum, containing six conserved regions with similarity to a high-affinity calcium-binding motif, the EF-hand. Recently, several studies have revealed that RCN1 acts as an oncogene involved in tumor progression. In our study, we found that the expression of RCN1 is higher in GBM samples than normal brain tissues by utilizing the GEPIA2 online tool, which is consistent with previous findings that high RCN1 expression is present in various malignancies (Amatschek et al., 2004). Recently, RCN1 has been demonstrated to be a prognostic marker in various cancers, such as non-small cell lung cancer  and renal cell carcinoma (Giribaldi et al., 2013). However, the function and the underlying mechanism of RCN1 involved in GBM remains a vague notion. Based on the observation in previous studies Giribaldi et al., 2013;Liu et al., 2018), we hypothesized that RCN1 could be a promising prognostic marker or as a response predictor for targeted therapy in GBM. Interestingly, we confirmed our hypothesis through TCGA and validated them in CGGA and experiments. At the same time, the mechanisms behind RCN1 were also explored and the identified pathways were consistent with previous studies (Fukasawa et al., 2021;Pucci et al., 2021;Sighel et al., 2021;Xue et al., 2021).
The convenient access to the public database allows us for the application of large-scale gene expression profiling and database mining for potential correlation between genes and overall survival of a variety of malignancies including GBM (Aldape et al., 2015;Zhou et al., 2021). In the present study, RCN1 was identified as an independent factor for the prognosis of GBM patients in TCGA, according to the results of univariate and multivariate analysis, with or without adjustment for confounders on the survival difference as described in a previous study (Erturk and Tas, 2017). By identifying the confounding factors, we further revealed the independent prognostic role of RCN1 in GBM patients. Since several clinical characteristics have been identified to influence the outcomes of GBM patients, interaction may exist (Brankovic et al., 2019) and it is plausible to perform subgroup analysis to estimate the interaction between RCN1 and selective clinical characteristics (Gönen, 2003). Here we innovatively identified confounding factors to increase the reliability of the results. Discrimination and calibration are the most commonly used indicators in evaluating prediction models. However, a systematic review found that 63% of the studies reported the discrimination information of prediction models, but only 36% of the studies reported the calibration information (Wessler et al., 2015); we reported both discrimination and calibration.
Recently, endoplasmic reticulum (ER) stress has been reported to have a considerable impact on cell growth, proliferation, metastasis, invasion, angiogenesis, and chemoradiotherapy resistance in various cancers. Accumulated studies have demonstrated that ER stress governs multiple pro-tumoural attributes in the cancer cell mainly via reprogramming the function of immune cells. The tumor microenvironment can be shaped because of the functional impact of ER stress responses in endothelial cells, cancer-associated fibroblasts, and other stromal cells during cancer progression (Chen and Cubillos-Ruiz, 2021). ER stress and downstream autophagy in the regulation of cell fate may function in temozolomide treatment and they have the potential to be therapeutic targets in GBM (He et al., 2019). Besides, it was revealed that the ER stress-related genes-based risk model can serve as a prognostic factor to predict the outcome for patients and be correlated with immune and inflammation responses in glioma . Among these genes, RCN1, as an ER-resident calciumbinding protein, is verified as one of ER stress-related genes, of which the depletion causes the ER stress-induced cells' apoptosis in various cancers (Huang et al., 2020;Liu et al., 2018;Xu et al., 2017). In the early study, it was demonstrated that RCN1 was identified as a genuine phagocytosis ligand to stimulate microglial phagocytosis of apoptotic neurons that were subsequently targeted by phagosomes (Ding et al., 2015). It was also observed that RCN1 interacts with SEC63p to activate protein translocation and quality control pathways in the ER (Honoré, 2009). As mentioned above, RCN1 proteins can not only have intrinsic roles in tumor development but also serve as a regulator involved in immune-related activities.
Considering the powerful role of tumor microenvironment cells contributing significantly to prognosis, we further investigated the immune cell infiltration between the highand low-risk groups in the TCGA and CGGA datasets, respectively. It was revealed that the infiltration levels of T cells CD4 memory resting, eosinophils, and macrophages M0 were differentially changed in both TCGA and CGGA cohorts. Furthermore, we found that macrophages M0 was significantly higher in the high-risk group both in the TCGA and CGGA datasets, representing a certain relationship between them. Furthermore, our results also suggested that high RCN1 may have a certain biological function in GBM. Therefore, to further investigate the role RCN1 played in GBM, we performed in vitro experiments and verified that the downregulation of RCN1 inhibited the cell invasion in GBM cell lines. The results of functional experiments are consistent with the above mentioned in our study, and further confirmed the critical role of RCN1 in GBM.
Several limitations should be noted in our study. Firstly, although we selected clinical samples from two different databases for mutual validation, verification in a larger sample is needed in the future. In addition, the variables involved in our prognostic model are easy to obtain, which is undoubtedly very convenient for clinical application, but it also needs to be further confirmed in clinical practice. Thirdly, although we verified the results of rigorous data mining, more animal or clinical exploration is still urgently needed. In addition, as for the absence of 1p19q characterization in TCGA, the status of 1p19q codeletion was not included for analyses. Finally, considering a lack of verification from clinical samples, the effectiveness of the RCN1-signature in GBM patients needs to be further verified by well-designed investigations in clinic. And how RCN1 regulates ER stress biology still needs to be further illustrated by experiments.
Overall, we found that high RCN1 has a poor OS for GBM patients, and confirmed RCN1 as an independent prognostic factor and developed a prognostic predictive model based on RCN1, which performed well in the prediction of OS for GBM patients.

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 authors.

AUTHOR CONTRIBUTIONS
WL performed the validation experiment and wrote the manuscript. HC collected the data and analysed in bioinformatics ways. BL supported the bioinformatics analyses and wrote the manuscript. CO analysed the statistics results. JX, QY, and MZ conceptualize the study and revised the original manuscript. JX also designed the validation experiment in vitro.