Identification of ITPR1 as a Hub Gene of Group 3 Medulloblastoma and Coregulated Genes with Potential Prognostic Values

The Group 3 Medulloblastoma (Grp3-MB) is an aggressive molecular subtype with a high incidence of metastasis and deaths. In this study, were used an RNA sequencing data (RNA-Seq) from a Brazilian cohort of MBs to identify hub genes associated with the metastatic risk. Data validation were performed by using multiple large datasets from MBs (GSE85217, GSE37418, and EGAS00001001953). DESeq2 package in R software was used to identify the differentially expressed genes (DEGs) in our RNA-Seq data. The DEGs data were accessed to construct the modules/graphs of co-expression and to identify hub genes through Cytoscape platform. The coregulated genes were enriched by the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway, and the protein–protein interaction (PPI) network was visualized by Cytoscape. The Kaplan–Meier plotter and ROC curves were used to validate the diagnostic and prognostic values of specific biomarkers identified through this model. We identified that inositol 1,4,5-trisphosphate receptor type 1 (ITPR1) as a downregulated hub gene, with a high diagnostic accuracy to Grp3-MBs and associated with tumor metastasis. In addition, we identified genes significantly correlated with ITPR1 that were associated with metastasis in Grp3-MB (ATP1A2, MTTL7A, and RGL1) and worst overall survival in MBs (ANTXR1 and RGL1). Our findings suggest that the ITPR1 hub gene is potentially involved in the metastatic process for Grp3-MB. Our data also provide evidence of targets that may serve as prognostic predictors and/or regulators for the metastatic process that maybe explored for further research of individualized therapy to Grp3-MBs.


Introduction
Medulloblastoma (MB) is the most common malignant embryonal brain tumor in the pediatric population (Northcott et al. 2019). The current consensus established by the World Health Organization (WHO) recognizes four main molecular and clinico-pathological distinct MB subtypes: winglessactivated (WNT-MB), sonic hedgehog-activated (SHH-MB), and non-WNT/non-SHH MBs: Group 3 (Grp3-MB) and Group 4 (Grp4-MB) (Louis et al. 2021). Grp3-MBs represent approximately 25% of all MB cases, being the deadliest across molecular subgroups and frequently (15-20%) showing metastasis at diagnosis in infants and older children (Cavalli et al. 2017;Juraschka and Taylor 2019).
Next-generation sequencing (NGS) techniques such as RNA sequencing (RNA-Seq) and robust bioinformatics analysis have been used over the years to uncover potential biomarkers in MBs (Guo et al. 2020b). Although integrative genomic studies have clarified several mechanisms on MBs tumorigenesis, (Juraschka and Taylor 2019;Menyhárt and Győrffy 2019), it is of note the paucity of genes with clinical relevance that ultimately could be used as metastatic biomarkers and predictors for prognostication, improving therapeutic strategies for Grp3-MB.
In this study, we identified the inositol 1,4,5-trisphosphate receptor 1 (ITPR1) gene, a regulator of calcium signaling (Ca 2+ ), homeostasis, and autophagy in cancer (Gambardella et al. 2020;Sneyers et al. 2020) as a hub gene for Grp3-MBs. In addition, we showed that ITPR1 is specifically downregulated with a high diagnostic accuracy to Grp3-MB in our and multiple independent cohorts of MB. Interestingly, the low expression of ITPR1 was associated with the presence of metastasis in MBs. In addition, we identified 17 genes significantly coregulated with ITPR1, highlighting the ATP1A2, MTTL7A, and RGL1 which demonstrated association with metastasis in Grp3-MB, and ANTXR1 and RGL1 correlated to worse overall survival (OS) in MBs. Our study provides for the first-time insights about the association of ITPR1 as a prognostic biomarker and a potential therapeutic target to Grp3-MBs.

RNA-Sequencing from a Data Generated In-House Cohort
In a previous study of our group, 17 pediatric MB samples from Brazilian patients (clinical data are summarized in Supplementary Table 1) were sequenced and deposited in the GEO database under there Cord number GSE181293. The RNA-Seq data of our in-house cohort was used to identify the differentially expressed genes (DEGs) and also to gain insights on hub genes for Grp3-MBs as follows.

Co-expression Network and Hub Gene Selection
DESeq2 package in R software (version 3.4.3) was used to identify the DEGs into our RNA-Seq data (Love et al. 2014). The counts expression data (for Grp3-MB) were normalized and used as input for network construction in CEMiTool (followed script standard) (Russo et al. 2018) which DEGs with similar expression patterns were grouped into modules via the average linkage hierarchical clustering. Within the networks of co-expression, the genes that have the highest number of connections and/or interactions were identified, using the connectivity module measured by the absolute value of the Pearson correlation coefficient > 0.5, and the graphs were generated to represent interacting genes (hub genes) through the extension called GeneMania (Montojo et al. 2010) in the Cytoscape 3.7.1 platform (https:// cytos cape. org/). We screened genes closely linked in the modules and assessed the correlation of these genes expression value with presence of metastasis using the clinical data of the pediatric MBs GSE85217 (n = 643) (Cavalli et al. 2017) that were downloaded from R2 Platform (Genomics Analysis and Visualization Platform -http:// r2. amc. nl). We used the Student's t test (non-parametric Mann-Whitney test), and a p value ≤ 0.05 was considered as statistically significant. Thus, the ITPR1 hub gene was retained for further analysis.

Evaluation of ITPR1 Expression and Identification of Coregulated Genes in Multiple Pediatric MBs Retrospective Cohorts
Different datasets were downloaded from R2 Platform (Genomics Analysis and Visualization Platform -http:// r2. amc. nl) to evaluate ITPR1 expression levels. Two out of three datasets were derived of microarray data: GSE85217 (n = 643) (Cavalli et al. 2017); GSE37418 (n = 73) (Robinson et al. 2012); and one of RNA-Seq data: EGAS00001001953 (n = 186) (Northcott et al. 2017). ITPR1-coregulated genes were identified using common elements of all four MB datasets with a p value < 0.05 and a Pearson's correlation coefficient > 0.5. Further, a Venn diagram was performed to overlap coregulated genes with ITPR1 in all cohorts (data not shown). To better visualize the profile expression of ITPR1-coregulated genes, we generated a heatmap using the ComplexHeatmap package (Gu et al. 2016) with pre-set metrics as Pearson's correlation distance parameter and the algorithm WardD2.

Protein-Protein Interaction (PPI) Network
The protein-protein interactions (PPI) among the common ITPR1-coregulated genes were evaluated using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http:// www. strin gdb. org/) (Szklarczyk et al. 2015). Validated interactions with a combined score > 0.9 were selected as significant (Bozhilova et al. 2019) Then, integration of PPI networks was constructed using the Cytoscape software (https:// cytos cape. org) (Shannon 2003) with Gene-MANIA plugin to identify biological network integration for gene prioritization (Warde-Farley et al. 2010).

KEGG Pathway Enrichment analysis
We conducted Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis to better understand the putative functions of ITPR1-coregulated genes using the Enrichr database (https:// maaya nlab. cloud/ Enric hr/) (Kanehisa 2000). A p value of < 0.05 was considered to represent a statistically significant enrichment of DEGs in pathway.

Survival and Metastasis Analysis
Kaplan-Meier curves and the long-rank tests were performed to evaluate overall survival (OS) according to the expression of the ITPR1-coregulated genes using the GSE85217 dataset (Cavalli et al. 2017). This dataset was preferred to this analysis due to the higher number of MB cases available to be analyzed, along with better clinical data associated to genetic dataset. Patients were divided into two groups (low vs. high) according to the median expression of each gene. First, we considered the data from all subtypes of MBs samples primary (n = 343) and metastasis (n = 167). Then, we used only data of Grp3-MBs primary samples (n = 62) and metastasis (n = 43) information.

Statistical analysis
We performed statistical analysis comparing expressions among two different groups by Student's t test (non-parametric Mann-Whitney test). Analyses between more the two groups were used One-way-ANOVA non-parametric followed Tukey post test. The receiver operating characteristic (ROC) curve was applied to evaluate the potential value of ITPR1 expression levels in discriminating molecular MBs subtypes using the GSE85217 and to examine the diagnostic capacity of ITPR1 expression for distinguishing MBs and normal tissue using expression data of EGAS00001001953 and GSE3526 datasets (Roth et al. 2006). The accuracy was determined by the area under the curve (AUC). Graphs were generated with GraphPad Prism 8.0 (GraphPad Software, San Diego, CA, USA), considering p value < 0.05 as statistically significant.

RNA-Seq Analysis Reveals the ITPR1 as a Hub Gene in Our Pediatric Cohort of MB
In order to identify key genes in Grp3-MB tumor, we first identified DEGs in Grp3-MB compared to other molecular subtypes (WNT-MB, SHH-MB, and Grp4-MB) in our cohort. Next, co-expression network analysis identified a total of five specific key modules related to Grp3-MB (Supplementary Fig. 1). In one out of these modules, we identified ITPR1 as a hub gene (Fig. 1). Previously, the ITPR1 was found downregulated in Grp3-MBs highligthing a potential therapeutic target (Castillo-Rodríguez et al. 2018). Thus, we Fig. 1 Gene networks of coexpression module displaying the most connected genes (hub genes). Analysis of coexpression reveal ITPR1 as hub gene in Group 3-MBs. The size of the node is proportional to its degree. The genes with co-expression are represented in blue, and the genes with interactions are written in red. The module was constructed by CEMiTool package selected the ITPR1 as an interesting gene related with tumor biology for subsequent analysis given that ITPR1 into Grp3-MBs carcinogenesis are poorly understood to date.

The Low Expression of ITPR1 Gene Is Associated with Metastasis and Represents a Potential Biomarker in Grp3/Grp4-MBs
In our cohort, we showed that low expression gene of ITPR1 in Grp3-MBs when compared to other MBs subtypes (Fig. 2a). Thus, the expression profile of ITPR1 was also evaluated in three additional independent databases GSE85217 (Fig. 2b), GSE37418 (Fig. 2c), and EGAS00001001953 (Fig. 2d) of MBs. Interestingly, it was possible to validate a specific low expression profile of ITPR1 in Grp3-MBs or Grp4-MBs compared to WNT and SHH-MBs.
Once we observed that ITPR1 is downregulated in both Grp3 and Grp4-MBs, we evaluated its potential as a biomarker of these two specific subgroups. ROC analysis . The profile of ITPR1 gene expression was significantly observed downregulated in Grp3/Grp4-MBs when compared with other molecular subgroups of MB. ANOVA were conducted using the Tukey multiple comparisons post test to assess the statistical significance of the differences. * indicates (p < 0.001). e ROC curve analysis for prediction of G3 and G4 MBs among pediatric patients by ITPR1 level. f Low levels of ITPR1 are associated with metastasis in MB samples primary (n = 343) and metastasis (n = 167) (p < 0 .05) showed that the AUC was 0.846, suggesting a high diagnostic accuracy in the discrimination between Grp3/Grp4-MBs from non-Grp3/Grp4-MBs (Fig. 2e). At a threshold of 6.245 of ITPR1 expression, the optimal sensitivity and specificity were 82.1% and 77.4%, respectively. In addition, ITPR1 expression presented a significant diagnostic value for separating MB samples from healthy cerebellum (cut-off value of 9.151 normalized expression, AUC 0.999, 95% CI 0.997-1.000, sensitivity 99.5%, specificity 88.9%) as illustrated in (Supplementary Fig. 2). Remarkably, the lower expression of ITPR1 was also associated with the presence of metastasis in MB (Fig. 2f). Together, these results suggest that the low expression of ITPR1 gene is a potential biomarker of Grp3/Grp4-MBs and is potentially related to metastasis process in this setting.

ITPR1-Coregulated Genes Potentially Regulate Pathways Related to Metastasis in MBs
To gain further insights into the putative functions of the hub gene ITPR1, we identified genes significantly coregulated with ITPR1 expression exploring data of  (n = 186). Red indicates that the expression of genes is relatively upregulated, and blue indicates that expression of genes is relatively downregulated. The heatmap was generated using the ComplexHeatmap package and Pearson distance as metric. Clustering using the Ward.D2 algorithm three MB public databases (GSE85217, GSE37418, and EGAS00001001953) and our RNA-Seq data. A total of 17 genes were shown to be commonly coregulated with ITPR1 in all four datasets, which roles previously described in cancer are summarized in Supplementary Table 2. Interestingly, the expression profile of these 17 ITPR1-coregulated genes in four datasets converges to a low expression of these genes to Grp3-MBs in the most of samples, as shown in Fig. 3a-d. The PPI network and KEGG enrichment analysis revealed that the seventeen ITPR1-coregulated genes present genetic and co-expression interactions (Fig. 4a) and are significantly associated with the signaling pathways regulating pluripotency of stem cells, regulation of actin cytoskeleton, Ras signaling pathway, and cell tight junction (Fig. 4b). The genes that are involved into each correlated pathway are described in Supplementary Table 3. Altogether, our results indicate an important regulatory network between ITPR1 and coregulated genes that potentially are involved into a signaling net associated with metastatic processes in MBs.

ITPR1-Coregulated Genes Are Associated with Metastasis and Poor Prognosis in MBs
To evaluate the prognostic value and clinical relevance of the 17 ITPR1-coregulated genes, we performed association and survival analysis using the GSE85217 dataset (Cavalli et al. 2017). Interestingly, the low expression of all the 17 genes was significantly associated with the occurrence of metastasis when analyzed in all MBs datasets (p < 0.001) (Supplementary Fig. 3). In addition, we observed that the expression of three out of 17 common-genes, the ATP1A2, MTTL7A, and RGL1 were also significantly associated with presence of metastasis specifically in Grp3-MBs patients (p < 0.001) (Fig. 5a-c). Regarding prognostic value of the 17 ITPR1coregulated genes, lower expression levels of ANTXR1 and RGL1 were significantly associated with worse OS in MBs (p = 0.0037 and p = 0.029, respectively, Fig. 5d-e). Collectively, our findings suggest that ITPR1 and coregulated genes are potential biomarkers for of poor prognosis and metastasis in MBs. Fig. 4 a PPI analysis reveals interactions potential of ITPR1 and its coregulated genes. Circles represent genes and the results within the circle represent the structure of proteins. The ITPR1-RGL1-ANTXR1 axis regulation indicated interaction of genes associated with worse prognosis. ITPR1-ATP1A2-MTTL7A-RGLI suggests a putative gene interaction coregulated with metastatic risk potential in Group 3-MBs. Lines colors represent evidence of the interaction of proteins between genes: genetic interactions (line in blue) and co-expression (line in red). b KEGG pathway enrichment analysis in MBs. Abbreviations: KEGG, Kyoto Encyclopedia of Genes and Genomes. The significant enriched KEGG terms in MBs were based on their functions. Enrichment was done by Enrichr, using KEGG databases. The grayscale bars represent the p value, increased from the bottom up

Discussion
"Omics" studies (e.g., exome/proteomics/single-cell approaches) and integrative large NGS data analysis are valuable tools to identify novel biomarkers in cancer. Although many studies regarding biomarkers for MBs have been published (Huang et al. 2020), a sharper identification of hub genes and their interactions can be very helpful to unravel molecular mechanism associated with prognosis and metastasis, particularly for Grp3-MBs.
In the current study, we demonstrated that ITPR1 is a hub gene downregulated in Grp3-MBs and associated to metastasis. The ITPR1 gene is highly expressed in the Purkinje cells in the central nervous system (SNC), and the modulation of intracellular calcium homeostasis is the main biological function of this gene (Foskett et al. 2007;Nakanishi et al. 1991). Recently, the ITPR1 has been reported to present a complex role in cancer dynamics by controlling cell survival and cell death via regulation of apoptosis and autophagy (Parys et al. 2021). This is especially relevant, given that programmed cell death through autophagy plays crucial role in progression and metastasis in cancer (Su et al. 2015) and the low expression of autophagy-related genes is associated with poor prognosis and tumor progression (Wang et al. 2015a;Xu et al. 2019). Based on these reports and in our results, we suggest that the downregulation of ITPR1 is involved in the crosstalk between reduced autophagy and increased metastatic mechanisms and may play an important role towards a malignant phenotype in Grp3-MBs. Yet, further studies are needed to test and to explore this hypothesis.
Previously, the ITPR1 was found downregulated in Grp3/ Grp4-MBs highlighting a novel and promising therapeutic target. Meanwhile, this data was not explored the underlying molecular mechanisms of the ITPR1 into Grp3/Grp4-MBs carcinogenesis (Castillo-Rodríguez et al. 2018). We screened the gene expression profile of ITPR1 among our and others three pediatric MBs dataset: GSE85217 (Cavalli et al. 2017); GSE37418 (Robinson et al. 2012); and EGAS00001001953 (Northcott et al. 2017), and the low expression of ITPR1 gene was also found on Grp3/Grp4-MBs when compared to WNT/SHH-MBs. Of note, our ROC analysis showed that low expression of ITPR1 gene is a potential predictive biomarker for segregating Grp3/Grp4-MBs subtypes.
Furthermore, by screening the degree of nodes in the PPI network and KEGG enrichment analysis, we showed that overlapping of these ITPR1-coregulated genes present genetic and co-expression interactions and were significantly associated with the signaling pathways regulating pluripotency of stem cells, regulation of actin cytoskeleton, Ras signaling pathway and tight junction. The Ras signaling was Kaplan-Meier survival plots showing that low-expression of the d ANTXR1 and e RGL1 reveals statistical significance associated with poor survival in MBs patients p = 0.0037 and 0.029, respectively. Red, low-risk group; blue, high-risk group. X-axis, time (months); Y-axis, overall survival probability reported as regulator key element of several aspects of normal cell growth and malignant transformation. Further, it was also associated with metastatic dissemination and poor survival in breast cancer (Wright et al. 2015). Surprisingly, we showed that ITPR1-coregulated genes are also involved in relevant hallmarks of cancer as metastasis, cell invasion, migration, differentiation, and proliferation, as described in Supplementary Table 2. Intriguingly, we also observed that in most of the samples of pediatric MB tissues from all four cohorts, 17 ITPR1-coregulated genes were also found downregulated in the majority of samples of Grp3-MBs. Together, these ITPR1-coregulated genes are involved in pathways associated to the metastatic markers of Grp3-MBs.
Finally, we also found that low gene expression of ATP1A2, MTTL7A, and RGL1 were significantly associated with the occurrence of metastasis, particularly for Grp3-MBs. In addition, we demonstrated that low gene expression of RGL1 was associated with worse OS in MBs. Studies addressing the role of these genes in tumorigenic processes also are very scarce. Interestingly, the low gene expression of ATP1A2 was also found in breast cancer as a predictive and a prognostic factor (Bogdanov et al. 2017). Yet, the association between the downregulation of RGL1, metastasis, and worse OS on MBs was neither described nor explored to date.

Conclusion
In summary, this study highlights ITPR1 as a hub gene for Grp3-MB. Using an integrative bioinformatics approach, ITPR1 is potentially involved in the metastatic process of MBs. Nevertheless, further investigations are needed to fully verify and validate the function of this gene in MB setting and to balance their role in the underlying mechanisms of Grp3-MB metastization.