Elevated Hexose-6-Phosphate Dehydrogenase Regulated by OSMR-AS1/hsa-miR-516b-5p Axis Correlates with Poor Prognosis and Dendritic Cells Infiltration of Glioblastoma

Objective Glioblastoma (GBM), a type of malignant glioma, is the most aggressive type of brain tumor and is associated with high mortality. Hexose-6-phosphate dehydrogenase (H6PD) has been detected in multiple tumors and is involved in tumor initiation and progression. However, the specific role and mechanism of H6PD in GBM remain unclear. Methods We performed pan-cancer analysis of expression and prognosis of H6PD in GBM using the Genotype-Tissue Expression Project (GTEx) and The Cancer Genome Atlas (TCGA). Subsequently, noncoding RNAs regulating H6PD expression were obtained by comprehensive analysis, including gene expression, prognosis, correlation, and immune infiltration. Finally, tumor immune infiltrates related to H6PD and survival were performed. Results Higher expression of H6PD was statistically significantly associated with an unfavorable outcome in GBM. Downregulation of hsa-miR-124-3p and hsa-miR-516b-5p in GBM was detected from GSE90603. Subsequently, OSMR-AS1 was observed in the regulation of H6PD via hsa-miR-516b-5p. Moreover, higher H6PD expression significantly correlated with immune infiltration of dendritic cells, immune checkpoint expression, and biomarkers of dendritic cells. Conclusions The OSMR-AS1/ miR-516b-5p axis was identified as the highest-potential upstream ncRNA-related pathway of H6PD in GBM. Furthermore, the present findings demonstrated that H6PD blockading might possess antitumor roles via regulating dendritic cell infiltration and immune checkpoint expression.


Introduction
Glioma tissue derived from the nerve epithelium is a diffuse, highly invasive cerebral tumor with a poor prognosis. Glioblastoma (GBM), a type of malignant glioma, is the most aggressive variety of brain tumor and is associated with high mortality. GBM patients typically have short survival times (approximately 1 year) and an extremely poor prognosis, with a 5-year overall survival (OS) rate of less than 4% [1]. The prognosis of GBM patients remains dismal despite multiple therapies such as radical surgery, chemotherapy, and radiotherapy [2,3]. Thus, detailed underlying molecular mechanisms and prognostic predictors in GBM patients are urgently needed.
Hexose-6-phosphate dehydrogenase (H6PD) is involved in the downregulation of nicotinamide adenine dinucleotide phosphate (NADPH) in the endoplasmic reticulum lumen via converting glucose 6-phosphate to 6-phosphogluconate [4,5]. The H6PD gene is upregulated in multiple malignancies, and early studies reported that H6PD led to cancer cell growth [6,7]. Knockdown of H6PD reduced the proliferation and migration of breast cancer cells [5]. Downregulation of H6PD with short interfering RNAs in the murine colon and breast cancer cells correlated with reduced glucose uptake and decreased proliferation [8]. Silencing of H6PD contributed to reduced invasion and migration of gallbladder cancers [7]. Furthermore, increased H6PD plays a crucial role in purine metabolic reprogramming and the progression of glioma stem cells [9].
Nevertheless, the prognosis and regulatory mechanisms of H6PD in GBM remain scarcely explored. Moreover, the potential function of H6PD in regulating tumor immune cell infiltration in GBM is incompletely understood. In this study, we performed survival analysis for H6PD in a pan-cancer study. Next, the potential upstream (noncoding RNA) pathways of H6PD, including long noncoding RNAs (lncRNAs) and microRNAs (miRNAs), are also explored in GBM. Finally, we explored the correlation between H6PD expression and immune infiltration in GBM. Collectively, our study suggests that elevated H6PD regulated by OSMR-AS1/hsa-miR-516b-5p correlates with tumor immune infiltration and poor prognosis in GBM patients.

TCGA Data Download, Process, and Analysis
Pan-cancer gene expression data were acquired from The Cancer Genome Atlas (TCGA) database (https://tcga-data.nci.nih.gov/tcga/, accessed on 15 February 2022). We analyzed the differential expression of H6PD using the R package [10]. We also analyzed weighted Pearson correlations and p-values using the R package "weights" (https://CRAN.R-project.org/package=weights); p-values less than 0.05 were considered statistically significant.

GEPIA Database Analysis
GEPIA (http://gepia.cancer-pku.cn/detail.php, accessed on 15 February 2022) is a web-based tool for gene expression and correlation analysis based on TCGA and data from the Genotype-Tissue Expression Project (GTEx) [11]. GEPIA was used to analyze the expression of H6PD and long noncoding RNA (lncRNA) in different types of cancers. Survival analysis for H6PD in pan-cancer studies, including overall survival (OS) and disease-free survival (DFS), was conducted using GEPIA. In addition, prognostic values of the candidate lncRNAs were assessed by GEPIA and the R package. PDCD1, CD274, CTLA4, SIGLEC15, HAVCR2, TIGIT, LAG3, and PDCD1LG2 were selected as major immune checkpoints. The associations between H6PD and immune checkpoints in GBM were investigated using data from the GEPIA database.

ENCORI Database Analysis
ENCORI (http://starbase.sysu.edu.cn/, accessed on 17 February 2022) is a web server for analyzing RNA-related studies [12]. Candidate miRNAs were obtained using EN-CORI. Several target prediction programs obtained predictive miRNAs of H6PD, including RNA22, miRmap, PITA, microT, miRanda, PicTar, and TargetScan. In addition, parameters for degradome data (low stringency), clip data (medium stringency), and pan-cancer type (one cancer type) were set. Only predicted miRNAs acquired in at least three programs were considered candidate miRNAs of H6PD. The correlation between miRNAs and H6PD in GBM was analyzed by using ENCORI. miRNAs negatively correlated with H6PD were Brain Sci. 2022, 12, 1012 3 of 12 selected for subsequent survival analysis. In addition, candidate lncRNAs that might bind to candidate miRNAs were obtained from ENCORI.

Differential Expression Analysis by microRNA-seq
The miRNA sequencing of glioblastoma is publicly available from the GEO under accession number GSE90603. The reads were mapped to the human genome (hg38). The resulting matrix of reading counts was then loaded into RStudio (R version 3.6.3, R core Team, 2020, R Foundation for Statistical Computing, Vienna, Austria). The differential expression analysis was performed using the R/Bioconductor package DESeq (version 1.16.0, R Foundation for Statistical Computing, Vienna, Austria).

GEO2R Analysis
GEO2R (http://www.ncbi.nlm.nih.gov/geo/geo2r, accessed on 21 February 2022), a web tool available from the NCBI, was applied to identify differentially expressed genes between GBM and normal samples. GEO2R carried out comparisons using the GEOquery and Limma R packages from the Bioconductor project to identify the differentially expressed miRNAs in GBM and normal brain samples. The statistical significance of miRNA differential expression was determined by assigning |Log2 fold change| (|log2FC|) values and adjusted p-values. An adjusted p-value < 0.05 and the absolute|log2FC| > 1 were set as the cut-off criteria.

Prediction of lncRNA and ceRNA Network Construction
Based on the ceRNA hypothesis [13,14], we analyzed the positive correlation between H6PD and targeted lncRNAs. We constructed a lncRNA-miRNA-mRNA interaction network for H6PD to understand the gene regulation of H6PD.

UCSC Xena Database and Kaplan-Meier Plotter Analysis
The UCSC Xena database (http://xena.ucsc.edu/, accessed on 22 February 2022) is a visualization and analysis web tool for correlations between genomic and phenotypic variables. The database contains TCGA data. The database provided information on survival outcomes and gene expression. The expression and survival curves of lncRNAs were acquired by combining the GEPIA database and the "survival" package-derived R Project.

TIMER Database Analysis
TIMER (https://cistrome.shinyapps.io/timer/, accessed on 21 February 2022) is a classical and authoritative database for evaluating immune infiltrates and clinical implications in diverse cancer types. TIMER was used to generate the correlation of H6PD expression with immune infiltrates in GBM. The survival module investigated the association between clinical outcomes and immune infiltrates.

Immune Infiltration Analysis
Tumor immune infiltrates were obtained using a single GSEA (ssGSEA) method with a GSVA R package based on TCGA datasets [15]. The correlation between H6PD and eight immune cell types was determined using the Spearman correlation test. Graphs were obtained using the ggplot2 R package (version 3.6.3, R core Team, 2020, Vienna, Austria). The GEPIA database was utilized to identify the relationships between H6PD and biomarkers of immune cells.

UACLAN
UALCAN is designed to provide easy access to a publicly available database; it is a comprehensive and interactive online tool that can be used to analyze cancer OMICS data [16]. Graphs and plots describing survival information of lncRNAs in GBM patients were obtained from UALCAN.

The Human Protein Atlas (THPA) Analysis
THPA was used to assess the H6PD expression of GBM and normal brain tissues based on immunohistochemistry.

Statistical Analysis
H6PD expression analysis was performed with TIMER, GEPIA, THPA, and R projects using the "ggplot2" package. Survival analysis was performed with GEPIA, TIMER, ENCORI, and R projects. Expression of miRNAs was performed with R projects. Expression of lncRNAs was obtained from GEPIA. The association between H6PD and immune checkpoints in GBM was evaluated using GEPIA. Kaplan-Meier survival analysis and Cox regression were also utilized to assess the prognostic value of H6PD expression. A value of p < 0.05 was considered statistically significant in all analyses.

Expression and Survival Analysis for H6PD in Pan-Cancer Studies
Differences in H6PD were detected in 14 cancers based on TCGA cancer and standard data ( Figure 1A

Analysis of Candidate miRNAs and lncRNAs Binding to H6PD
According to the ceRNA hypothesis [13,14], there is a positive correlation between lncRNA and H6PD, while a negative correlation exists between miRNA and H6PD. We obtained three candidate miRNAs that could competitively bind to H6PD: hsa-miR-124-3p, hsa-miR-516b-5p, and hsa-miR-506-3p ( Figure 3A). Downregulation of hsa-miR-124-3p and hsa-miR-516b-5p was detected from GSE90603, while downregulation of hsa-miR-506-3p was not. Subsequently, 27 candidate lncRNAs related to hsa-miR-124-3p were upregulated in GBM ( Figure 3B-D). Statistical significance of DLEU1 for predicting the prognosis of GBM was detected, while statistical significance was not found for the other 26 lncRNAs (Figure 3E,F). OS and DFS analyses demonstrated that GBM patients with higher DLEU1 expression had a good clinical prognosis. Thirty-five candidate lncRNAs related to hsa-miR-516b-5p were upregulated in GBM ( Figure 3G-J). Statistical significance of

Analysis of Candidate miRNAs and lncRNAs Binding to H6PD
According to the ceRNA hypothesis [13,14], there is a positive correlation between lncRNA and H6PD, while a negative correlation exists between miRNA and H6PD. We obtained three candidate miRNAs that could competitively bind to H6PD: hsa-miR-124-3p, hsa-miR-516b-5p, and hsa-miR-506-3p ( Figure 3A). Downregulation of hsa-miR-124-3p and hsa-miR-516b-5p was detected from GSE90603, while downregulation of hsa-miR-506-3p was not. Subsequently, 27 candidate lncRNAs related to hsa-miR-124-3p were upregulated in GBM ( Figure 3B-D). Statistical significance of DLEU1 for predicting the prognosis of Brain Sci. 2022, 12, 1012 6 of 12 GBM was detected, while statistical significance was not found for the other 26 lncRNAs ( Figure 3E,F). OS and DFS analyses demonstrated that GBM patients with higher DLEU1 expression had a good clinical prognosis. Thirty-five candidate lncRNAs related to hsa-miR-516b-5p were upregulated in GBM ( Figure 3G-J). Statistical significance of OSMR-AS1 for predicting the prognosis of GBM was detected, while statistical significance was not found for the other thirty-four lncRNAs ( Figure 3K,L). GBM patients with higher OSMR-AS1 expression were negatively associated with OS and RFS. It was found that OSMR-AS1 was positively correlated with H6PD, while DLEU1 was negatively correlated with H6PD ( Figure 3M,N). Considering expression analysis, survival analysis, and correlation analysis, OSMR-AS1 might be the highest-potential upstream lncRNA of the hsa-miR-516b-5p/H6PD axis in GBM.

Association of Immune Infiltration, Overall Survival, Immune Checkpoints, and H6PD Level in GBM
Infiltration levels of various immune cells under different copy numbers of H6PD in GBM were provided. Significant changes in immune cell infiltration levels were observed under various copy numbers of H6PD in GBM, including B cells, CD8+ T cells, CD4+ T Figure 3. Expression, survival analysis, correlation of candidate miRNAs, lncRNAs, and H6PD. To ascertain whether some lncRNAs modulated H6PD, the upstream miRNAs bound to H6PD were predicted, and finally, three miRNAs, namely hsa-miR-124-3p, hsa-miR-516b-5p, and hsa-miR-506-3p, were found (A). Downregulation of hsa-miR-124-3p and hsa-miR-516b-5p was detected from GSE90603, while downregulation of hsa-miR-506-3p was not. Twenty-seven candidate lncRNAs related to hsa-miR-124-3p were upregulated in GBM (B-D). DLEU1 was significantly associated with prognosis in GBM patients, while the other 26 lncRNAs were not (E,F). OS and DFS analyses demonstrated that GBM patients with higher DLEU1 expression had a good clinical prognosis. Thirty-five candidate lncRNAs related to hsa-miR-516b-5p were upregulated in GBM (G-J). Statistical significance of OSMR-AS1 for predicting prognosis of GBM was detected, while statistical significance of the other 34 lncRNAs was not (K,L). GBM patients with higher OSMR-AS1 expression were negatively associated with OS and RFS. DLEU1 was negatively correlated with H6PD (M), while OSMR-AS1 was positively correlated with H6PD (N). (ns: no significance; * p-value < 0.05; ** p-value < 0.01; *** p value < 0.001).

Association of Immune Infiltration, Overall Survival, Immune Checkpoints, and H6PD Level in GBM
Infiltration levels of various immune cells under different copy numbers of H6PD in GBM were provided. Significant changes in immune cell infiltration levels were observed under various copy numbers of H6PD in GBM, including B cells, CD8+ T cells, CD4+ T cells, and neutrophils. We found a significant positive correlation between H6PD expression and dendritic cell infiltration ( Figure 4C). Higher dendritic cell infiltration and expression of H6PD were significantly positively associated with poor outcomes (Figure 4D,E). A correlation was found between biomarkers of dendritic cell infiltration and H6PD in GBM. CD8A, HLA-DPB1, TGAX, NRP1, SIRPA, THBD, and XCR1 were positively associated with H6PD expression (Figure 5 A-G). The association of H6PD with biomarkers of dendritic cells in GBM is summarized in Table 1. Higher expression of eight checkpoints in GBM was detected ( Figure 6A). The correlations between H6PD and immune checkpoints in GBM were assessed ( Figure 6B-I). H6PD expression was significantly positively correlated with CD274, CTLA4, LAG3, PDCD1LG2, TIGHT, and SIGLEX15.

Discussion
Malignant glioma is the most common primary brain tumor type and has a poor prognosis. According to previous studies, the median survival time for GBM patients after treatment is approximately one year [17]. Despite improvement in surgery and chemoradiotherapy for GBM, the clinical prognosis of GBM patients is extremely poor. Accumulating evidence indicates that H6PD plays an essential role in tumor initiation, progression, and drug resistance in multiple human cancers, including gallbladder, prostate, and breast cancer [6,7]. However, the function, molecular mechanism, and immune infiltration of H6PD in GBM remain unclear.
The study conducted a pan-cancer analysis of H6PD's expression using TCGA and GTEx data. H6PD was upregulated in GBM and downregulated in CESC, PCPG, and UCEC, suggesting that H6PD may play an oncogenic role in these four cancers. Although CECS, PCPG, and UCEC are different cancer types and come from different tissues, they share the same result from H6PD downregulation. H6PD controls cancer cell proliferation and migration through pleiotropic effects, unfolded protein response, calcium homeostasis, and redox balance, which may result in a downregulation in CECS, PCPG, and UCEC . Future studies assessing gene expression might provide necessary insight into the biological mechanism. Low H6PD expression was detected in normal brain tissues and low-grade glioma. Nevertheless, abundant H6PD expression was found in GBM and related to a poor prognosis. H6PD was not statistically significant in predicting the prognosis of patients with other cancer types. Previous studies also indicated that upregulation of H6PD plays an essential role in purine metabolic reprogramming and glioma stem cell progression [9].
The ceRNA hypothesis is a novel mechanism for RNA interactions. Based on the ceRNA hypothesis, lncRNA acts as a ceRNA to interact with mRNA by competing with the corresponding miRNA [13]. The ceRNA hypothesis has provided a novel perspective that ceRNA-miRNA-mRNA interactions (or lncRNAs acting as ceRNAs) can regulate many cancers' biological processes for cancer therapy [27,28]. Thirty-five candidate lncRNAs related to hsa-miR-516b-5p were upregulated in GBM. It was found that OSMR-AS1 works as a prognosticator for GBM in the present study. GBM patients with higher expression of OSMR-AS1 possessed poor OS and RFS. A positive correlation between OSMR-AS1 and H6PD was detected. OSMR-AS1 was found to be a differentially expressed lncRNA and an independent prognosis factor for patients with hepatocellular carcinoma [29]. OSMR-AS1 was detected as a survival-related lncRNA in glioblastoma using TCGA databases [30]. Moreover, OSMR-AS1 was a lncRNA related to pro-neural to mesenchymal transition and was positively associated with poor prognosis in glioma [31].
Immune checkpoint blockading promotes immune infiltration and antitumor immune response [32]. Immune infiltration can contribute to tumor progression by altering the tumor microenvironment and extracellular matrix [14]. The relationship between H6PD and immune response remains unclear. Our work suggested that significant changes in immune infiltration with various copy numbers of H6PD, including B cells, CD8+ T cells, CD4+ T cells, and neutrophils, were observed in GBM. H6PD expression was significantly positively associated with dendritic cell infiltration. Moreover, higher immune cell infiltration of dendritic cells and higher expression of H6PD were significantly positively associated with poor outcomes. Subsequently, we found that biomarkers of dendritic cells, including CD8A, HLA-DPB1, TGAX, NRP1, SIRPA, THBD, and XCR1, were positively associated with H6PD expression. Additionally, we investigated the association between H6PD and immune checkpoints. Higher expression of eight checkpoints in GBM was detected ( Figure 6A). The correlation between H6PD and immune checkpoints in GBM was assessed ( Figure 6B-I). H6PD expression was significantly positively correlated with CD274, CTLA4, LAG3, PDCD1LG2, TIGHT, and SIGLEX15, indicating that targeting H6PD might increase the efficacy of immunotherapy in GBM.

Conclusions
H6PD was highly expressed and correlated with an unfavorable prognosis in GBM. The OSMR-AS1/miR-516b-5p axis was identified as the highest-potential upstream ncRNArelated pathway of H6PD in GBM. Furthermore, the present findings demonstrated that H6PD blockading might possess antitumor roles via regulating dendritic cell infiltration and immune checkpoint expression. However, basic experimental studies and clinical trials need to be carried out in the future to further validate these bioinformatic analysis results for GBM.

Conflicts of Interest:
The authors declare that they have no conflict of interest.