Identication 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 Brazinian 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, 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 coexpression and to identify hub genes through Cytoscape platform. The co-regulated 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 specic biomarkers identied through this model. We identied that Inositol 1,4,5-triphosphate 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 identied genes signicantly 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 ndings 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: Wingless-activated (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 clari ed 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 metastastatic biomarkers and predictors for prognostication, improving therapeutic strategies for Grp3-MB.
In this study, we identi ed 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 speci cally 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 identi ed 17 genes signi cantly 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 rst-time insights about the association of ITPR1 as a prognostic biomarker and a potential therapeutic target to Grp3-MBs.

Materials And Methods
RNA-Sequencing from a data generated in-house cohort In a previous study of our group, 17 pediatric MB samples from Brazilian patients 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 coexpression were identi ed the genes that have the highest number of connections and/or interactions, using the connectivity module measured by the absolute value of the Pearson correlation coe cient > 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://cytoscape.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 signi cant. Thus, the ITPR1 hub gene was retained for further analysis.
Evaluation of ITPR1 expression and identi cation of co-regulated 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 identi ed using common elements of all four MB datasets with a p-value <0.05 and a Pearson's correlation coe cient > 0.5. Further, a Venn diagram was performed to overlap coregulated genes with ITPR1 in all cohorts (data not shown). To better visualize the pro le expression of ITPR1-coregulated genes we generated a heatmap using the Complex Heatmap 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 commom ITPR1-coregulated genes were evaluated using the Search Tool for the Retrieval of Interacting Genes (STRING) database (http://www.stringdb.org/) (Szklarczyk et al. 2015). Validated interactions with a combined score >0.9 were selected as signi cant (Bozhilova et al. 2019) Then, integration of PPI networks was constructed using the Cytoscape software (https://cytoscape.org) (Shannon et al. 2003) with GeneMANIA plugin to identify biological network integration for gene prioritization (Warde-Farley et al. 2010).

KEGG pathways 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://maayanlab.cloud/Enrichr/) (Kanehisa and Goto, 2000). A p-value of < 0.05 was considered to represent a statistically signi cant 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 analysed, 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). Analyzes between more the two groups were used One-way-ANOVA non-parametric followed Tukey post test. Receiver operating characteristic (ROC) curves were used to evaluate the discrimination of ITPR1 expression levels into molecular MBs subtypes. 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 signi cant.

Results
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 rst identi ed DEGs in Grp3-MB compared to other molecular subtypes (WNT-MB, SHH-MB and Grp4-MB) in our cohort. Next, co-expression network analysis identi ed a total of ve speci c key modules related to Grp3-MB ( Supplementary Fig. 1). In one out of these modules, we identi ed ITPR1 as a hub gene (Fig. 1). Previously, the ITPR1 was found downregulated in Grp3-MBs highlighting a novel and promising therapeutic target (Beckmann et al. 2019;Castillo-Rodríguez et al. 2018). Thus, we 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 pro le of ITPR1 was also evaluated in three additional independent databases GSE85217 ( Once we observed that ITPR1 is downregulated in both Grp3 and Grp4-MBs, we evaluted its potential as a biomarker of these two speci c subgroups. ROC analysis 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). Remarkably, the lower expression of ITPR1 was also assoaciated 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 identi ed genes signi cantly coregulated with ITPR1 expression exploring data of three MBs public databases (GSE85217, GSE37418, 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 1. Interestingly, the expression pro le 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 revelead that the seventeen ITPR1coregulated genes present genetic and co-expression interactions (Fig. 4a) and are signi cantly associated with the signaling pathways regulating pluripotency of stem cells, regulation of actin cytoskeleton, Ras signaling pathway and cell tigh-junction (Fig. 4b). The genes that are involved into each correlated pathway are described in Supplementary Table 2. Altogether, our results indicate an important regulatory network between ITPR1 and co-regulated 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 signi cantly associated with the ocurrance of metastasis when analyzed in all MBs datasets (p<0.001) Supplementary Fig. 2). In addition, we observed that the expression of three out of 17 common-genes, the ATP1A2, MTTL7A and RGL1 were also signi cantly associated with presence of metastasis speci cally in Grp3-MBs patients (p<0.001) (Fig. 5a-c). Regarding prognostic value of the 17 ITPR1-coregulated genes, lower expression levels of ANTXR1 and RGL1 were signi cantly associated with worse OS in MBs (p=0.0037 and p=0.029, respectively- Fig. 5d-e). Collectively, our ndings suggest that ITPR1 and coregulated genes are potential biomarkes for of poor prognosis and metastasis in MBs. Discussion a "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 identi cation 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 (Beckmann et al. 2019;Castillo-Rodríguez et al. 2018). We screened the gene expression pro le 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 signi cantly associated with the signaling pathways regulating pluripotency of stem cells, regulation of actin cytoskeleton, Ras signaling pathway and tigh junction. The Ras Signaling was 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, diferention and proliferation, as described in Table 1. Intriguingly, we also observed that in most of the samples of pediatric MBs 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 signi cantly associated with the ocurrance 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 were neither described, nor explored to date.

Conclusion
In summary, this study highlights ITPR1 as a hub gene for Grp3-MB. Using an integrative bioinformatic 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.

Ethical Statement
Ethics approval and consent to participate: This study was approved by the HC/FMRP-USP Research Ethics Committee (Approval number: CAAE: 62604816.4.000.5440) and informed consent was obtained from all patients included.
Consent for publication: Not applicable.
Availability of data and materials: The data and other items supporting the results of the study will be made available upon reasonable request. Figure 1 Gene networks of co-expression module displaying the most connected genes (hub genes). Analysis of co-expression 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 Hierarchical unsupervised clustering of primary MB datasets. Samples were clustered according molecular subgroups of MBs: SHH (blue), WNT (purple), Group 3 (pink) and Group 4 (green). Heatmaps summarizing a downregulated pro le of expression genes positively coregulated with ITPR1 (FGL2,  ALDH1A, ATP1A2, FGFR2, CEND1, GR1A4, METTL7A, NAV2, RGL1, JAM2, PCGF5, ZCCH24, ANTXR1, WWTR1, ZFP36L1, RILPL2, MSN) in (a) our cohort (n=17); (b) GSE85217(n=643), (c) GSE37418(n=73) and (d) EGAS00001001953 (n=186). Red indicates that the expression of genes is relatively upregulated, 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 Figure 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-ANTXR1axis-  Analysis of ITPR1-gene expression with biological and clinical importance in MBs samples by GSE85217: