Identification of Key Pathways and Genes in the Orai2 Mediated Classical and Mesenchymal Subtype of Glioblastoma by Bioinformatic Analyses

Background Ca2+ release-activated Ca2+ channels (CRAC) are the main Ca2+ entry pathway regulating intracellular Ca2+ concentration in a variety of cancer types. Orai2 is the main pore-forming subunit of CRAC channels in central neurons. To explore the role of Orai2 in glioblastoma (GBM), we investigated the key pathways and genes in Orai2-mediated GBM by bioinformatic analyses. Methods Via The Cancer Genome Atlas (TCGA), French, Sun, and Gene Expression Omnibus (GEO) (GDS3885) datasets, we collected 1231 cases with RNA-seq data and analyzed the functional annotation of Orai2 by gene ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway analysis. Univariate and multivariate survival analyses were applied to 823 patients with survival data. Results We discovered that Orai2 was markedly upregulated in GBM compared to normal brain samples and lower-grade gliomas (LGG). Survival analysis found that higher expression of Orai2 was independently associated with a worse prognosis of patients with the classical and mesenchymal subtypes of GBM. Simultaneously, Orai2 expression was higher in tumors of the classical and mesenchymal subtypes than other subtypes and was significantly correlated with classical- and mesenchymal-related genes. GO and KEGG pathway analysis revealed that genes significantly correlated with Orai2 were involved in the JNK pathway. Through screening transcriptomic data, we found a strong association between Orai2 and apoptosis, stemness, and an epithelial-mesenchymal transition- (EMT-) like phenotype. Conclusion As a prognostic factor, Orai2 is obviously activated in the classical and mesenchymal subtypes of GBM and promotes glioma cell self-renewal, apoptosis, and EMT-like by the JNK pathway. These findings indicate that Orai2 could be a candidate prognostic and therapeutic target, especially for the classical and mesenchymal subtypes of GBM.


Introduction
Gliomas are the most common primary intracranial tumors. Glioblastoma multiforme (GBM) is the most aggressive form (WHO grade IV) with high morbidity among different types of gliomas [1]. The standard treatment methods start with maximal surgical resection, followed by radiotherapy with concomitant and adjuvant temozolomide. However, most patients will die within 1 to 2 years. A median progressionfree survival of 6.2 to 7.5 months and a median overall survival of 14.6 to 16.7 months have been reported in a clinical trial [2]. One of the reasons for this high mortality is the therapeutic resistance of GBM, mainly due to its diffuse growth, genetic heterogeneity, and microvascular proliferation [3].
The discovery of glioma stem cells (GSCs) is an important event in glioma research, which has profoundly affected our study of GBM patterns, tumor heterogeneity, tumor cell and microenvironment relationships, and treatment strategies [4]. Glioma stem cells, also known as glioma-initiating cells (GICs), are tumor cells with self-renewal and infinite differentiation potential and have a high degree of invasion, migration and tumorigenicity. GSCs maintain the vitality of tumor cell populations through self-renewal and immortalization, and the ability to move and migrate makes it possible to spread tumor cells. GSCs can also be dormant for a long time and have multiple drug-resistant molecules that are insensitive to the external physical and chemical factors that kill tumor cells and are the drivers of tumor progression and recurrence [5].
Epithelial-mesenchymal transition (EMT) is a process in which epithelial properties are impaired and interstitial properties are enhanced [6]. The EMT in GBM has gradually attracted increased attention in recent years. However, what occurs during the progression of GBM is not complete EMT but a transformation process similar to epithelialmesenchymal transition, called EMT-like, which refers to a process in which the epithelial phenotype is less pronounced and more biased toward the mesenchymal phenotype, mainly manifested by downregulation of epithelial markers such as E-cadherin and increased expression of interstitial markers such as N-cadherin and fibronectin, which leads to a decrease in cell-to-cell adhesion and loss of cell polarity, ultimately enhancing tumor invasion and migration [7]. Therefore, as one of the most challenging malignancies worldwide, effective therapies for GBM are needed.
Since Bailey and Cushing first proposed the systematic classification of neuroepithelial tissue tumors in 1926, the concept of histogenesis has dominated CNS tumor classification for almost a century. The histology-based WHO classification and grading system serves as a "gold standard," playing an important role in the diagnosis and treatment of CNS tumors. In addition, the pathological diagnostic criteria of gliomas, especially GBMs, which have caused artificial heterogeneity and complexities in investigations, did not have useful effects on clinical management. In 2016, the WHO CNS tumor histology classification (revised version 4) integrated the histological phenotype and gene phenotype for the first time. IDH, 1p/19q, TP53, pTERT, ATRX, and other indicators were included in the integrated diagnostic criteria for gliomas. Since then, the diagnosis and treatment of gliomas have entered the era of molecular typing [8]. However, glioma diagnosis with the integration of histological and gene phenotype information still needs further improvement. Verhaak et al. analyzed the expression profiles of 202 GBMs in The Cancer Genome Atlas (TCGA) and classified GBMs into four subtypes: proneural, neural, classical, and mesenchymal. Each subtype differs greatly in terms of its cellular features, genetic contexts, and signaling pathways involved. In TCGA classification of GBMs, different molecular subtypes do not share the same response and survival benefit to traditional treatment [9]. The molecular typing of GBM makes the inclusion and grouping of target patients in clinical trials more homogenous and balanced. This system also helps to identify important molecular targets in each subtype to design feasible therapeutic strategies.
The c-Jun NH2-terminus kinases (JNKs) are a family of stress-activated serine threonine protein kinases of the MAPK pathway. Mitogen-activated protein kinases (MAPKs) are serine-threonine protein kinases that regulate various cellular activities including proliferation, differentiation, apoptosis, survival, inflammation, and innate immunity. Compromised MAPK signaling pathways have also been found to play key roles in the pathogenesis of various diseases including cancer and neurodegenerative disorders. Three distinct genes, known as JNK1, JNK2, and JNK3, encode ten splice variants. Four JNK1 isoforms, four JNK2 isoforms, and two JNK3 isoforms have been identified. JNK1 and JNK2 are expressed in most tissues, whereas JNK3 is expressed predominantly in neuronal cells [10][11][12]. JNK is involved in the regulation of various signaling pathways and plays an important role in cell proliferation, apoptosis, and differentiation [13]. Therefore, JNK has been implicated in a variety of human diseases, including tumors. Increasing evidence from animal and clinical studies suggests a critical role for aberrant activation of JNK in tumor development, and JNK is drawing increasing attention as a promising target of antitumor therapy in the near future [14,15]. Recent studies have shown that the JNK signaling pathway plays an important role in the malignant progression of GBMs [16].
Ca 2+ release-activated Ca 2+ (CRAC) channels mediate Ca 2+ influx in many cell types. These channels are formed by the tetraspanning plasma membrane proteins, named Orai1, Orai2, and Orai3. These ORAI proteins mediate Ca 2+ influx by store-operated Ca 2+ entry (SOCE) [17]. CRAC channels are hexameric complexes composed of individual or potentially multiple Orai homologues. The transmembrane domains are highly conserved between all three Orai homologues [18], and all the Orai homologues can function as Ca 2+ channels when they are overexpressed. The properties of ectopically expressed Orai1, Orai2, and Orai3 channels are similar to those of endogenous CRAC channels, including activation by Ca 2+ store depletion, high Ca 2+ selectivity, inward rectification, and Ca 2+ -dependent inactivation (CDI). However, the three ORAI homologues differ in some of their channel properties, including fast and slow CDI and their sensitivity to pharmacological inhibitors such as 2aminoethoxydiphenyl borate [19].
Recent results of our research team showed that Orai1 is upregulated in glioma specimens and glioma cell lines compared with nontumor control brain tissue, and the regulatory degree is positively correlated with the grades of gliomas. Further analysis found that the use of the inhibitor SKF96365 or RNAi downregulation of Orai1 significantly inhibited the invasion and migration of glioma cells and the epithelial-mesenchymal transition-(EMT-) like process.
In vivo experiments showed that Orai1 can inhibit C6 colloidal downregulation of invasion and migration of tumor cells in a rat brain [20]. Motiani et al. found significant upregulation of Orai1 in primary GBM cell lines, and Orai1 knockdown significantly inhibited GBM cell proliferation and invasiveness [21]. However, as a homolog of Orai1, the role of Orai2 in gliomas has not yet been elucidated. In the present study, we aimed to evaluate the Orai2-related key genes and pathways by analyzing RNA-seq data across microarray chips including a total of 1231 glioma samples, and to the best of our knowledge, this is the largest and most comprehensive study of the main role and internal mechanism of Orai2 in GBM.

Samples and Patients.
Transcriptome and clinical data of glioma patients were retrieved from The Cancer Genome Atlas (TCGA) dataset (http://cancergenome.nih.gov/); French and Sun datasets were available on R2 Genomics Analysis and Visualization Platform (http://r2.amc.nl or http://r2platform.com) and Microarray data of GEO (http://www.ncbi.nlm.nih.gov/geo); GSE23806 was enrolled in this study. GSE23806 profile included microarray data from 17 cases of neurospheres, 11 primary tumors (corresponding tumor tissue of glioblastoma patients), 27 cases of glioblastoma stem-like cell lines, and 36 cases of conventional glioma cell lines. The mRNA expression level of 539 GBM samples was measured through AffyU133a microarray platform. Meanwhile, 530 lower-grade glioma (WHO II and WHO III) and 172 GBM patients were measured through an Illumina HiSeq microarray platform. All of these data are available from the TCGA database. Both the French (284 cases included) and Sun (180 cases included) were acquired from R2 Genomics Analysis and Visualization Platform. Cluster analyses were performed and visualized using Cluster/Java Treeview. Patients were divided into the Orai2 low and Orai2 high groups by the mean expression levels. ChIPBase v2.0 (http://rna.sysu.edu.cn/chipbase/) was used to perform Pearson's correlation tests. The Retrieval of Interacting Genes/Proteins (STRING database v10.5) resource was utilized for PPI network analysis and prediction of protein-protein interactions.

Identification of DEGs.
Morpheus, an R-associated web application, was applied to filtrate DEGs among four subtypes in GBM. The P ≤ 0:01 and | log FC| ⩾ 1 were considered the cutoff criterion. All results of DEGs were downloaded in text format, with hierarchical clustering analysis being conducted later in Morpheus (https://software .broadinstitute.org/morpheus/).

GO and Pathway Enrichment Analysis of DEGs.
The online tool, Database for Annotation, Visualization and Integrated Discovery (DAVID, https://david.ncifcrf.gov/) provided comprehensive information for the list of genes by GO and KEGG pathway analyses. In addition, GO enrichment analysis included three different aspects: biological process (BP), molecular function (MF), and cellular component (CC) [22]. KEGG enrichment analysis was associated with genomic information's functional interpretation and practical application [23]. P < 0:05 was considered to indicate a statistically significant difference.

Statistical Analysis.
In this study, the prognostic value of Orai2 was estimated by Kaplan-Meier analysis and Cox proportional hazard model analysis using SPSS statistical software (version 19). Other statistical computations and drawing of figures were performed with several packages (ggplot2, pheatmap, and circlize) in the statistical software environment R, version 3.3.2 (http://www.r-project.org). For all statistical methods, P < 0:05 was considered a significant difference.

Results
3.1. Survival Analysis of Orai2 Expression Status. TCGA datasets retrieved from the GEPIA (Gene Expression Profiling Interactive Analysis) showed that among the several types of tumors that are common in humans, the tumors with the highest expression of Orai2 are glioblastoma multiforme (GBM), diffuse large B-cell lymphoma (DLBCL), acute myeloid leukemia (AML), pheochromocytoma and paraganglioma (PCPG), and thyroid carcinoma (THCA) (Figure 1(a)). We analyzed Orai2 mRNA expression in the French and Sun Brain datasets from the R2 analysis and visualization platform. The mRNA expression results showed that Orai2 was overexpressed in GBM samples compared with normal brain tissues (Student's t-test, P < 0:05) (Figure 1(b)). Compared with lower-grade glioma (WHO II and WHO III), Orai2 expression was significantly increased in GBM from the TCGA dataset (Student's t-test, P < 0:001) (Figure 1(c)). When we classified the data according to glioma tissue type in the French and TCGA datasets, we found that the expression level of Orai2 in GBMs was also higher than that of astrocytomas, oligoastrocytomas, and oligodendrogliomas. There was no significant difference in the expression of Orai2 in astrocytomas, oligoastrocytomas, and oligodendrogliomas (Figures 1(d) and 1(e)). To further confirm the Orai2 expression results in GBM, we analyzed data in TCGA. The results revealed that compared with other subtypes, Orai2 was significantly upregulated in the classical and mesenchymal subtypes of GBMs (Figure 1(f)). Accordingly, the Kaplan-Meier survival analysis revealed that patients with high levels of Orai2 had shorter overall survival (OS) times than patients with low levels of Orai2 in the French and TCGA datasets (Figure 1(g)). We evaluated the prognostic value of Orai2 by Kaplan-Meier survival curve analysis in four subtypes of GBM. The results showed that patients with higher Orai2 expression in the classical and mesenchymal subtypes had a shorter overall survival, while the other two subtypes had no prognostic value (Figure 1(h)).

Orai2 Expression Characteristics in Classical and
Mesenchymal Subtypes of GBMs. By the analysis of the TCGA database, Verhaak et al. found that the classical type has astrocyte characteristics and expresses neural precursors and stem cell markers (NES, Notch pathway, and Sonic hedgehog pathway), and its tumor cells are highly proliferative; the mesenchymal type has astrocyte features, which mainly express mesenchymal markers (CHI3L1, CD44, and vimentin) and astrocyte markers (GFAP, Sox9, and CD44), the high activity of their markers often suggests epithelialmesenchymal transition (EMT), and necrotic and inflammatory reactions often occur in tumors; classical and mesenchymal subtypes of GBM patients have a poor prognosis [9,24]. Orai2 RNA levels were strongly positively associated with  CESC  CHOL  COAD  DLBC  ESCA  GBM  HNSC  KICH  KIRC  KIRP  LAML  LGG  LIHC  LUAD  LUSC  MESO  OV  PAAD  PCPG  PRAD  READ  SARC  SKCM  STAD  TGCT  THCA  THYM  UCEC   classical subgroup markers (NOTCH1, NES, NOTCH1, RELA, and OLIG2) (Figures 2(a) and 2(c)) and mesenchymal subgroup markers (TGFB1, VIM, TIMP1, CD44, and CHI3L1) (Figures 2(b) and 2(d)). These results indicate that Orai2 is strongly associated with the classical and mesenchymal phenotypes in human GBMs.

Orai2 Clinical Expression Status.
In the TCGA dataset (n = 539), high expression of Orai2 was selectively associated with elder ages at first-time diagnosis and shorter overall survival years, days to progression, and days to recurrence. Furthermore, higher Orai2 expression was obviously related to shorter overall survival years and temozolomide (TMZ) chemoradiation and long-term therapy with TMZ in classical patients. Meanwhile, higher Orai2 expression was distinctly associated with shorter overall survival years and GBM recurrence in the mesenchymal samples (Supplementary Table 1). The above results suggest that Orai2 may be involved in the recurrence of GBMs and could be a new diagnostic and clinical prognostic biomarker, especially in the classical and mesenchymal types of GBMs.

Identification and Hierarchical Clustering Analysis of
DEGs. First, a total of 7193, 6056, 6635, and 6606 DEGs were obtained separately from TCGA datasets between classical or mesenchymal and other subtypes of GBM with the criteria of P < 0:05 and | log FC| ⩾ 1:0. Hierarchical clustering analysis was conducted through Morpheus, a web-based online tool, with the series matrix data of the DEGs. The heat map is shown in Figures 3(a)-3(e) (the top 100 genes). Next, 3404 and 3030 genes were acquired from classical and mesenchymal GBMs by the mean Orai2 expression (Figures 3(c) and 3(f)). Finally, we summarized the abovementioned hub genes that overlapped among groups to identify DEGs, yielding a total of 185 genes (Figures 3(g)-3(i)).  Table 2. We found that Orai2 was mainly involved in the regulation of signaling, nervous system development, regulation of the JNK cascade, glycoprotein biosynthetic process, regulation of transport, regulation of the stress-activated MAPK cascade, negative regulation of phosphorylation, regulation of I-kappaB kinase/NF-kappaB signaling, and positive regulation of receptor-mediated endocytosis (Figure 4(a)).

Orai2-Related
To investigate whether Orai2 expression is correlated with JNK pathway-related molecule expression, we analyzed data in TCGA and found that Orai2 was correlated with JNK pathway-related molecules in the classical and mesenchymal subtypes of GBM samples (Figure 4(b)). To further investigate the downstream mechanism of Orai2 in regulating the JNK pathway, we mapped the interactome of Orai2 with GSC self-renewal, apoptosis process, and EMTassociated proteins by using STRING (https://string-db.org) (Figure 4(c)). Furthermore, we evaluated the association between the JNK pathway and GSC self-renewal, apoptosis, and EMT from TCGA transcriptomic data. A strong correlation between Orai2 and GSC self-renewal, apoptosis, and EMT was observed (Figures 4(d) and 4(e)).
3.6. Association of Orai2 and Self-Renewal, Apoptosis, and EMT. One GEO database (GSE23806) was used to analyze the expression of Orai2 in glioblastoma stem-like cell lines, conventional glioma cell lines, neurospheres, and primary tumors (corresponding tumor tissue of glioblastoma patients). We found that Orai2 expression was relatively high in glioblastoma stem-like cell lines and neurospheres compared to conventional glioma cell lines and primary tumors ( Figure 5(a)). At the same time, we utilized another independent GEO microarray dataset (GDS3885) and found that Orai2 was closely related to stem-related genes that were significantly upregulated in GSCs relative to the matched nontumor stem cells (NSTCs) (Figure 5(b)). In addition, to test Orai2 as a biomarker and a regulator of glioma stem cells, we further found a relationship between Orai2 and GSC markers, such as PROM1, NES, SOX2, OLIG2, NOTCH1, SOX9, MET, NOTCH2, and DLL3. Orai2 was obviously To explore whether Orai2 is involved in EMT-like in classical and mesenchymal subtypes of GBM, we calculated the relationship between Orai2 and EMT-related markers by cluster analysis (Supplementary Fig. 1A, B) and Pearson's correlation analysis (Supplementary Fig. 1C, D). The results showed that orai2 was negatively correlated with epithelial marker CDH1 (E-cadherin) expression (CDH1 CLA ,  Fig. 1C, D). These results suggested that Orai2 is likely to be involved in the EMT-like process in GBM and plays a facilitating role. Finally, schematic of the Orai2/JNK signaling pathway is shown in Figure 5(f).

Discussion
In summary, this is the largest and most comprehensive study identifying the key pathways and genes in Orai2mediated GBM by bioinformatic analyses. We revealed that Orai2 expression was significantly higher in GBM (WHO IV) than lower-grade glioma (WHO II and WHO III). Patients with high expression of Orai2 have a poor prognosis compared with low expression of Orai2 in GBM patients. Furthermore, we found that the expression level of orai2 mRNA in the classical and mesenchymal types of GBM was significantly higher than that of other subtypes. Genes associated with Orai2 were mainly involved in GSC self-renewal, apoptosis, and EMT. In clinical terms, higher Orai2 expression was independently associated with a worse prognosis of patients with GBM, especially in the classical and mesenchymal types of GBM. Our group has previously studied and found the expression characteristics of orai1 in gliomas across all grades in 61 patients [20]. In this study, we collected and characterized orai2 expression levels in 1231 glioma samples from four large independent databases. We found that the mRNA levels of Orai2 were higher in GBM than in nonneoplastic brain tissues, indicating that Orai2 acted as a tumor promoter in GBM. Notably, our findings showed that Orai2 was expressed at relatively higher levels in the classical and mesenchymal subtypes. At the same time, we also found that Orai2 RNA levels were strongly positively associated with classical and mesenchymal subgroup markers. We concluded that Orai2 is strongly associated with the classical and mesenchymal phenotypes in human GBMs.
In this study, GO term enrichment and KEGG pathway analysis results revealed that Orai2 was involved in and regulated GSC self-renewal, apoptosis, and EMT. Simultaneously, according to our results, Orai2 was likely to regulate these biological processes through the JNK pathway. This hypothesis requires in vitro and in vivo experiments to further confirm it in the near future.
Some studies have shown that JNK, as a key molecule in the maintenance of GSC, may be a feasible molecular target for glioma [14,15]. Moreover, related studies have suggested that JNK plays a pivotal role in other malignant tumors, such as ovarian, pancreatic, and lung cancers. Okada et al. found that JNK inhibitors could prevent pancreatic tumor formation via eliminating cancer stem-like cells [25][26][27]. JNK has been extensively studied in the induction of apoptosis. Upon activation, phosphorylated JNK translocates to the nucleus where it phosphorylates and regulates transcription factors such as c-Jun, p53, and c-Myc, which are involved in the induction of cell apoptosis [28]. The alkylating anticancer drugs temozolomide (TMZ) and nimustine (ACNU) induce apoptosis in glioma cells at late times following treatment by activating JNK and c-Jun [29].
Interestingly, other studies have confirmed that the JNK inhibitor SP600125 was sufficient to inhibit cell viability, migration, and invasion of GBM cells through decreasing the phosphorylation levels of both c-Jun and Akt [30]. Combined inhibition of PI3K and JNK can inhibit the invasion and migration of glioblastoma and lamellipodia and membrane ruffle formation [31]. Likewise, Ueno et al. found that knocking down JNK in TMZ-resistant GBM cells could enhance invadopodia formation through the JNK-paxillin axis to accelerate ECM degradation [32].
However, to date, there have been few reports on the specific mechanism of action of Orai2 in glioblastoma, so our bioinformatics analysis of Orai2 is very meaningful. We have shown through bioinformatics analysis that orai2 may participate in the JNK pathway to exert antitumor effects, which provides direction for future research, and we will further test the relationship between these molecules. We hoped that these and subsequent related experiments could benefit GBM patients through targeting Orai2/JNK-related factors in the near future.

Data Availability
The data used to support the findings of this study are included within the article.

Conflicts of Interest
There is no potential conflict of interest.