Phytanoyl-CoA 2-Hydroxylase-Interacting Protein-Like Gene Is a Therapeutic Target Gene for Glioblastoma Multiforme.

Glioblastoma multiforme (GBM) is the most common primary CNS cancer and has a poor prognosis. This study searched for significant genes and the mechanisms involved in GBM. We used the Gene Expression Omnibus (GEO) to test the WHO normal and IV glioma database, used R tool to identify the significant gene, and finally, combined these with The Cancer Genome Atlas (TCGA) to verify the significant genes. Subsequently, we explored the biological mechanisms involved. Phytanoyl-CoA 2-hydroxylase-interacting protein-like gene (PHYHIPL) is downregulated in grade IV glioma (GBM). The downregulation of PHYHIPL in GBM is accompanied by poor overall survival in the TCGA database, which indicates that PHYHIPL is a protection gene in GBM development. Bioinformatics analysis shows that the poor prognosis with downregulated PHYHIPL may be the result of the TNF signaling pathway and the IL-17 signaling pathway, but good prognosis accompanied by upregulated PHYHIPL may be the result of retrograde endocannabinoid signaling and the cAMP signaling pathway. Protein-protein interactions (PPI) net indicated that PHYHIPL may play a vital role in cell metabolism, and we hypothesize that the downregulation mechanism may be the result of mutations of the ß-catenin gene and the endogenous siRNA, as shown in previous studies. PHYHIPL may be a target gene for the treatment and prognosis of GBM.

Background

Glioma is the most common type of primary intracranial tumor [1]. According to the WHO classification of tumors of the central nervous system, there are 4 grades [2], and these grades have a potential value for the determination of prognosis. Through molecular and technical developments, many important and clinically valuable genes have been identified, and epidermal growth factor receptor (EGFR) has been shown to be a regulatory gene for GBM [3]. Gene of phosphate and tension homology deleted on chromosome ten (PTEN) is a tumor-suppressor gene, which can also protect humans from tumor invasion [4]. The mutation of the TP53 gene is a typical mechanism in 18- to 45-year-old glioma patients [5]. Cell proliferation markers also play a vital role in the development of glioma, and a study has demonstrated that Ki67 labeling index and prognosis have statistical significance [6]. Although there are many certifications in the marker genes, they either have a positive or a negative function in glioma. The survive time of grade II glioma tumor is more than 5 years, but patients with high-grade glioma still have a poor clinical prognosis. The grade III glioma tumors survive 2–3 years, while the grade IV glioma tumors survive less than 1 year. Thus, it is necessary to find the marker genes in WHO high-grade glioma tumors. This study searched for the marker gene in grade IV glioma tumors (glioblastoma multiforme (GBM).

There are many methods used to find marker genes. High-throughput gene expression is a general method used to explore the marker genes. The GEO project was initiated in response to the growing demand for a public repository for high-throughput gene expression data. According to shared science data around the world, big data mining is an effective tool for exploring the mechanism of disease development, which is bioinformatics analysis. The Cancer Genome Atlas (TCGA) database comprises 33 types of cancers, and contributes to improving the prevention, diagnosis, and treatment of cancer. We used the GEO to find the marker genes, and then used the TCGA database to verify them. The work flow of the present study is shown in Figure 1.

Figure 1

The workflow of the study.

In this study, we selected a GEO that provides normal brain and grade IV glioma tumors expression profiling by the array. We found that 47 cases were diagnosed as normal brain (n=13) and GBM (n=34) of any histologic type in initial surgical treatment and fresh frozen material was obtained from an array of dataset GSE50161. The immune marker gene has demonstrated that specific immunotherapeutic approaches may be most effective for each tumor type [7], but there is still a large amount of significant information that needs further analysis.

Material and Methods

The data were downloaded from the TCGA and GEO, according to the TCGA policy, a sample informed consent document that TCGA suggests tissue source sites use when collecting specimens from prospective project participants. This document ensures that patients considering donating tissue specimens to TCGA recognize the risks and benefits of participation and understand the nature of their inclusion in the project. All of the 47 GSE50161 patients were invited to participate in this Institutional Review Board-approved study (COMIRB 95-500 and 09-0906) [7].

First, we selected the array of dataset GSE50161 for analysis, and the array of dataset GSE50161 was downloaded by the GEOquery R package [8]. The dataset contains normal brain (n=13) and GBM (n=34) tissues.

Second, we found the different expression protein coding Genes (DEGs) of GBM in GSE4412, (1) Data pre-processing: The raw data were preprocessed by the GEOquery R package [8]. (2) The gene annotation was by hgu133plus2.db R package [9]. The protein-coding genes were downloaded from Ensembl (). (3) DEG analysis. The DEGs in samples were analyzed by the Limma package [10]. We plotted a heat map and the volcano map of DEGs. It is clear the PHYHIPL is a notable DEG. Its logFC=−1.917168, AveExpr=10.60324, t=−4.061331, p.value=1.870214e-04, adj.p.value=9.084314e-04, so we used the TCGA to verify the marker gene.

Third, the overall survival of the marker gene in TCGA-GBM. TCGA portal () is a web tool used to analyze the prognosis of the marker gene, which was provided by Shutan Xu, Shaying Zhao Lab, University of Georgia. We chose the TCGA-GBM to analyze. The threshold was: the marker gene for gene name, mean ±SD for dividing patient, and choosing all for subtype.

Fourth, we used the GDCRNATools [11] to download the TCGA-GBM RNA-seq data; there are 159 samples in it, including 154 primary tumor samples and 5 normal samples. (1) We plotted PHYHIPL gene expression in 2 different groups, (2) we divided the 154 primary tumor samples into high- and low-expression group according to the mean expression of PHYHIPL gene, and used Limma [10] to explore the DEGs of the 2 group, clusterProfiler R package [12] for further bio-information exploration, the gene ontology (GO), and the Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway.

Finally, we used bioinformatics to explore the PPI net of the marker gene and used the published articles to explore the down-regeneration mechanism. (1) The exploration of the mechanism. We used String to analyze the PPI [13], which has 9’643’763 proteins from 2031 organisms and 1’380’838’440 interactions. It is a powerful web tool for PPI. (2) Although there are few related articles, we still found some significant value of the downregulation mechanism.

Results

DEG analysis

We found that there were 18 764 genes and 2625 different expression genes (DEGs), including 1203 upregulated and 1422 downregulated DEGs. The results are shown by the volcano map in Figure 2, and we showed the DEGs with a heat map in Figure 3. PHYHIPL gene has not been reported before and it was a notable DEG. Therefore, we selected the PHYHIPL gene to explore the expression and prognosis in TCGA-GBM.

Figure 2

The volcano plot by R tool of GSE 50161. The black point is the unchanged genes; the green point is the downregulated genes; the red is upregulated genes. The threshold was adjusted P-value <0.001 and |log2 fold change (FC)|=1.5 was set.

Figure 3

The heat map of GSE50161. The blue is the GBM; the red is normal brain tissues. The threshold was adjusted P-value <0.001 and |log2 fold change (FC)|=1.5 was set.

PHYHIPL gene overall survival in TCGA-GBM

The threshold was set as previously stated. Figure 4 shows that patients benefited from high expression of PHYHIPL, with the log-rank P=0.038.

Figure 4

Overall survival. The red is the PHYHIPL high-expression group; the black is the PHYHIPL low-expression group. Log-rank p=0.038.

PHYHIPL gene expression in TCGA-GBM

Figure 5 shows that PHYHIPL gene is expressed less in primary GBM than in normal tissue. t=−8.410, p<0.001, CI=(−2.743, −0.903).

Figure 5

PHYHIPL expression in TCGA-GBM. The primary tumor mean expression is 3873.513, the solid tissue normal mean expression is 9486.2, the t test was used in the comparison, the t=−8.410, p≤0.001.

PHYHIPL gene in primary GBM

We selected the primary GBM for further analysis and there are 154 primary tumor samples. The mean expression of PHYHIPL gene was 3873.513, so we divided them into 2 groups according to the mean expression: the high-expression sample was 59 and the low-expression sample was 95. DEG analysis also used the Limma package [10]. There were 16 158 genes in it, and the threshold was an adjusted p-value=0.05 and |log2 fold change (FC)|=1 was set. There were 499 upregulated genes and 170 downregulated genes. The volcano map was shown in Figure 6. Subsequently, we used the clusterProfiler R package [12] to explore the downregulated genes and upregulated genes in GO and KEGG pathways (Figures 7, 8).

Figure 6

The volcano plot of TCGA-GBM. The black point is the unchanged genes; the green point is the downregulated genes; the red is upregulated genes.

Figure 7

The downregulated DEGs in the GO and KEGG pathways. The threshold was set at minimal size of genes annotated by Ontology term for testing=10, the cutoff value of p-value is 0.05, q value cutoff is 0.2. The size of blots represents the enriched count numbers, the color of blots is the adjusted p-value. It only shows the top 10 enriched GO and the top 5 enriched KEGG pathway. (A) Cellular Component (CC); (B) Biological Process (BP); (C) Molecular Function (MF); (D) KEGG pathway.

Figure 8

The upregulation DEGs gene ontology and KEGG pathway. The threshold was set at minimal size of genes annotated by Ontology term for testing=10, the cutoff value of p-value is 0.05, q value cutoff is 0.2. The size of blots represents the enriched count numbers, the color of blots is the adjusted p-value. It only shows the top 10 enriched GO and KEGG pathways. (A) Cellular Component (CC); (B) Biological Process (BP); (C) Molecular Function (MF); (D) KEGG pathway.

The GO is constructed of Cellular Component (CC), Biological Process (BP), and Molecular Function (MF). The pictures only show the top 10 Clustering results. CC enrichment suggests that downregulated DEGs were significantly enriched in the extracellular matrix, proteinaceous extracellular matrix, endoplasmic reticulum lumen, and collagen trimer. BP enrichment suggests that downregulated DEGs were significantly enriched in extracellular structure organization, extracellular matrix organization, and skeletal system development. MF was significantly enriched in receptor-ligand activity, enzyme inhibitor activity, cytokine activity, and growth factor binding. The KEGG enrichment suggests that downregulated DEGs were significantly enriched in cytokine-cytokine receptor interaction, the TNF signaling pathway, the IL-17 signaling pathway, transcriptional misregulation in cancer, and AGE-RAGE signaling pathway in diabetic complications. CC enrichment suggests that upregulated DEGs were significantly enriched in postsynapse, synaptic membrane, and neuronal cell body. BP enrichment suggests that upregulated DEGs were significantly enriched in the regulation of membrane potential, modulation of chemical synaptic transmission, and the regulation of trans-synaptic signaling. MF was significantly enriched in channel activity, passive transmembrane transporter activity, and substrate-specific channel activity. The KEGG enrichment suggests that upregulated DEGs were significantly enriched in retrograde endocannabinoid signaling, morphine addiction, the cAMP signaling pathway, and GABAergic synapse.

PPI net

There were 10 proteins in PPI net, including StAR-related lipid transfer domain containing 7(STARD7), outer mitochondrial membrane lipid metabolism regulator (OPA3), enoyl-CoA hydratase 1 (ECH1), dehydrogenase/reductase 4 (DHRS4), peptidase, mitochondrial processing beta subunit (PMPCB), armadillo repeat containing 4 (ARMC4), transmembrane and coiled-coil domains 6 (TMCO6), phosphodiesterase 9A (PDE9A), complement C1q binding protein (C1QBP), and chromosome 1 open reading frame, human C21orf33 (C21orf33). The average node degree was 3.09 and the PPI enriched p-value is 0.0331 (Figure 9).

Figure 9

The PPI net of PHYHIPL by String. It shows STARD7, OPA3, ECH1, DHRS4, PMPCB, ARMC4, TMCO6, PDE9A, C1QBP, and C21orf33. The average node degree was 3.09 and the PPI enriched p-value was 0.0331.

Discussion

GBM is associated with poor prognosis and is the most common primary CNS cancer. The mechanism is not clear. It contains epigenetic alterations [14,15], mutation [16], copy numbers variations [17], DNA damage [18], microenvironment [19], and extracellular matrix [20]. The downregulated DEGs were significantly enriched in the extracellular matrix and the KEGG enrichment suggests that downregulated DEGs were significantly enriched in the TNF signaling pathway and the IL-17 signaling pathway. The KEGG enrichment suggests that upregulated DEGs were significantly enriched in retrograde endocannabinoid signaling and the cAMP signaling pathway.

The TNF signaling pathway plays a vital role in apoptosis, inflammation, and immunity. It is an important mediator of tumor progress [21]. EGFR gene amplification and mutation are common in GBM, but EGFR inhibitor is not active in GBM, mainly because the TNF-JNK-Axl-ERK signaling pathway mediates the resistance [22]. IL-17 is a pro-inflammatory cytokine that sustains a chronic inflammatory microenvironment that favors tumor formation. The IL-17R-Act1-TRAF4-MEKK3-ERK5 signaling pathway is vital for IL-17 to mediate keratinocyte proliferation and tumor formation [23]. Endocannabinoid signaling plays an antiproliferative role in GBM [24], and forskolin and 3-isobutyl-1-methylxanthine can inhibit the proliferation and proapoptotic activity in glioma stem cells through activating the cAMP signaling pathway [25].

We also found significant value in the PPI. STARD7 is a lipid transfer protein that is a mediator for phosphatidylcholine to mitochondria [26]. Want/β-catenin signaling can enhance the STARD7 expression [27], but in the adrenocortical adenomas, the PHYHIPL was downregulated according to the mutations of the β-catenin gene [28]. OPA3 had a protective function in mitochondria [29]. ECH1 plays a vital role in lymphatic metabolism [30]. DHRS4 is a kind of peroxisomal member [31]. PMPCB is the mediator for the majority of mitochondrial proteins [32]. ARMC4 is a new gene that has been found in peritoneal metastatic gastric carcinoma [33]. PDE9A is expressed highly in the brain, and it is a cGMP-specific phosphodiesterase [34]. C1QBP is a mitochondrial protein [35] that plays a vital role in promoting cell survival. Gene downregulation may be the result of miRNA and siRNA. It was reported that endogenous siRNA derived from RMRP results in off-target effects in humans, and PHYHIPL was identified as the most promising candidate [36].

Conclusions

PHYHIPL is expressed in the brain (RPKM57.9) (Gene ID: 84457). To the best of our knowledge, PHYHIPL has not been previously reported in GBM, so there is sparse related information on the function of PHYHIPL.

Bioinformatics analysis showed that the KEGG enrichment suggests that downregulated DEGs are significantly enriched in the TNF signaling pathway and the IL-17 signaling pathway, and they are the oncology pathways in GBM. The KEGG enrichment analysis suggests that upregulated DEGs are significantly enriched in retrograde endocannabinoid signaling and the cAMP signaling pathway, and they are the tumor-suppressor pathways in GBM. GBM prognosis may be connected to the pathways.

The PPI net showed that PHYHIPL may be a metabolism gene. Thus, we speculate that the PHYHIPL is related to GBM cell mitochondria, which regulates cell metabolism. It is well known that cell metabolism plays a vital role in tumor cell development. The decrease in PHYHIPL gene expression is due to mutations of the β-catenin gene and the endogenous siRNA in GBM. The mechanism should be explored in the future, which will contribute to the treatment and prognosis of GBM.

Data availability (where applicable)

The datasets analyzed are available in the Cancer Genome Atlas data portal and the gene expression omnibus database (GSE50161).