ABSTRACT
The pathogenesis of pheochromocytoma and paraganglioma (PCPG) catecholamine-producing tumors is exceedingly complicated. Here, we sought to identify important genes affecting the prognosis and survival rate of patients suffering from PCPG. We analyzed 95 samples obtained from two microarray data series, GSE19422 and GSE60459, from the Gene Expression Omnibus (GEO) repository. First, differentially expressed genes (DEGs) were identified by comparing 87 PCPG tumor samples and eight normal adrenal tissue samples using R language. The GEO2R tool and Venn diagram software were applied to the Database for Annotation, Visualization and Integrated Discovery (DAVID) to analyze Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways and Gene Ontology (GO). We further employed Cytoscape with the Molecular Complex Detection (MCODE) tool to make protein–protein interactions visible for the Search Tool for Retrieval of Interacting Genes (STRING). These procedures resulted in 30 candidate DEGs, which were subjected to Kaplan–Meier analysis and validated by Gene Expression Profiling Interactive Analysis (GEPIA) to determine their influence on overall survival rate. Finally, we identified ALDH3A2 and AKR1B1, two genes in the glycerolipid metabolism pathway, as being particularly enriched in PCPG tumors and correlated with T and B tumor-infiltrating immune cells. Our results suggest that these two DEGs are closely associated with the prognosis of malignant PCPG tumors.
INTRODUCTION
Pheochromocytoma and paraganglioma (PCPG) tumors originate from sympathetic nerve-derived adrenal medulla cells, as well as extra-adrenal thoracic and abdominal accessory ganglia, and show high vascularization and catecholamine secretion (Manger, 2006; Dahia, 2014). Patients with PCPG display symptoms including severe hypertension and sweating, which can lead to serious complications such as heart disease or cerebral hemorrhage. PCPG is a rare neuroendocrine tumor and about 10–20% of PCPG is metastatic, which results in an approximate 50% reduction in survival rate (Ghosal et al., 2020). The survival rate of patients with malignant PCPG remains low due to poor prognosis and lack of adequate treatment options. Currently, diagnosis of PCPG is difficult and malignant tumors are only detected when a metastasis is discovered. Since PCPG primary tumors may cause serious metabolic disorder that contributes to decreased survival rate (Martucci and Pacak, 2014), further investigations to identify novel prognostic molecular biomarkers to aid diagnosis, develop adequate therapeutics and understand the mechanism of PCPG pathogenesis are therefore of critical importance.
Previous studies have shown that some long noncoding RNAs (lncRNAs) are overexpressed in pan-cancers including PCPG (Posa et al., 2016; Wang et al., 2020) and are thus potential biomarkers. However, factors such as low abundance of circulating lncRNAs in body fluids, cost of analysis and low throughput mitigate against using lncRNAs as useful biomarkers for PCPG. To more efficiently explore the pathogenesis, diagnosis and therapy of cancers, newer techniques including gene chips have been widely utilized to detect differentially expressed genes (DEGs) between malignant tumors or sarcomas and normal tissues (Severgnini et al., 2006). These techniques have enriched the discovery of novel cancer biomarkers, which can be stored in databases. Consequently, these data will be applicable in the discovery of reliable clues for deeper research on cancers in the future. Bioinformatic tools have revealed new biomarkers of PCPG (Burnichon et al., 2011; Crona et al., 2013, 2015), proving that the extensive use of bioinformatic techniques can produce valuable information on the pathogenesis of cancers.
In this paper, we screened two PCPG data series, GSE19422 and GSE60459, from the Gene Expression Omnibus (GEO), a public gene microarray repository. We utilized the online GEO2R and R language to analyze differentially regulated genes in each series. Next, we employed Venn diagram software to detect common DEGs in both series. We also applied the Database for Annotation, Visualization and Integrated Discovery (DAVID) to check the molecular function (MF), cellular component (CC), biological process (BP) and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways of DEGs. Finally, we confirmed key genes and their edges from the DEGs using a protein–protein interaction (PPI) network online tool as well as the Cytoscape MCODE (Molecular Complex Detection) tool. To test whether these key genes are applicable to the prognosis of PCPG, we re-analyzed them using the Kaplan-Meier Plotter online database to show the relationship between their expression and patient survival rate. The key genes highlighted by Kaplan–Meier analysis were imported into Gene Expression Profiling Interactive Analysis (GEPIA) to compare their expression between PCPG tumors and normal adrenal tissue These manipulations resulted in eight DEGs, whose potential relationship with PCPG was analyzed using KEGG pathway enrichment. This yielded two DEGs, ALDH3A2 and AKR1B1, that are highly enriched in glycerolipid metabolism and that correlate closely with T and B lymphocytes of tumor-infiltrating immune cells (TIICs). We thus applied bioinformatic methods to identify two sensitive and effective biomarkers in PCPG patients that may provide useful clues regarding the mechanism of pathogenesis, diagnosis and treatment of PCPG tumors.
MATERIALS AND METHODS
Gene expression microarray
We downloaded the gene expression profiles in GSE19422 and GSE60459 for PCPG tumors and normal adrenal tissue from the GEO (www.ncbi.nlm.nih.gov/geo). We analyzed 95 samples in total, comprising 84 PCPG tumor samples and five normal adrenal tissue samples from GSE19422, based on the GPL6480 platform, and three PCPG tumor samples and three normal adrenal tissue samples from GSE60459, based on GPL13607. Thus, eight adrenal tissues were used as control and 87 tumor samples were used individually to identify DEG levels.
Differential expression analysis
We applied R language and GEO2R online tools to analyze DEGs between PCPG tumor specimens and normal adrenal specimens, with the P-value cut-off set to < 0.05 and |logFC| set to > 2. The DEGs from these two data sets were then analyzed by Venn software to identify common genes DEGs with logFC > 2 were regarded as up-regulated, while those with logFC < –2 were considered as down-regulated genes (Feng et al., 2019).
Gene Ontology and pathway enrichment analyses of DEGs
Gene Ontology (GO) analysis was performed to identify the functional roles of DEGs and proteins in (The Gene Ontology Consortium, 2000). KEGG, which is widely used in analyzing functional proteins and genes, was applied to define the nodes of molecular networks (Kanehisa et al., 2017).
DAVID is used to integrate DNA, RNA and protein functions (Zou et al., 2019). DAVID was used to describe the enrichment of MF, BP and CC.
PPI network and module analysis
PPI networks provide information about how proteins interact with other proteins to carry out their cellular functions. We used the STRING database to make a comprehensive annotation and analysis of the functional relationship of genes (Szklarczyk et al., 2017). Cytoscape was applied to visualize the PPI network. In addition, the MCODE tool was applied to screen every module of the interaction (degree cutoff = 2, max. depth = 100, k-core = 2, and node score cutoff = 0.2) (Feng et al., 2019).
Kaplan–Meier survival analysis
The prognostic impacts of the genes were checked with the Kaplan-Meier Plotter (https://kmplot.com/analysis), a public web tool for confirming the influence of a large group of genes on cancer survival whose sources include the GEO, EGA and TCGA databases. We also applied GEPIA, a tool for delivering overall functional results (Tang et al., 2017). Expressed RNA sequence data for 9,736 tumors and 8,587 normal samples, including 182 PCPG tumor patients and three normal people, from the TCGA and GTEx projects are served by GEPIA.
Tumor-infiltrating immune cells analysis
The prognostic results for the selected genes were checked using the Tumor Immune Estimation Resource (TIMER) database tool, which contains abundant information on TIICs and gene expression levels for 23 types of cancers from the TCGA database. We analyzed the correlation of eight screened significant survival-related key genes’ expression with PCPG TIICs (B cells, CD4+ T cells, CD8+ T cells, neutrophils, macrophages and dendritic cells). Correlation of tumor purity level with key gene expression was also considered. We also analyzed the correlation of different somatic copy number alterations (SCNAs) of the eight key genes with tumor immune lymphocytes in PCPG tumors using the SCNA module to assess whether the mutation of these genes will influence the PCPG tumor microenvironment. Genomic Identification of Significant Targets in Cancer (GISTIC 2.0) was utilized to define SCNAs that contained deep deletion (–2), arm-level deletion (–1), diploid/normal (0), arm-level gain (1) and high amplification (2). The infiltration levels of tumor immune lymphocytes were compared between PCPG tumors and normal tissue using the two-sided Wilcoxon rank-sum test (Kärjä et al., 2005).
RESULTS
Identification of common DEGs in PCPG tumors
In total, we identified 635 DEGs from GSE19422 and 1,402 from GSE60459 using the GEO2R online tool, and their expression was obtained by the R software package with |logFC| > 2 and P < 0.05 as cut-off points. We found 343 up-regulated genes and 292 down-regulated genes in GSE19422, and 474 up-regulated genes and 928 down-regulated genes in GSE60459 (Fig. 1A). Common DEGs in these two datasets were identified by Venn diagram software. As demonstrated in Fig. 1B, the Venn diagram revealed 151 common DEGs. In all, there were 30 up-regulated genes (logFC > 2) and 121 down-regulated genes (logFC < –2) in PCPG tumor tissues (Table 1).
Table 1. List of common differentially expressed genes from GSE19422 and GSE60459 data series containing 30 up-regulated and 121 down-regulated genes in PCPG tumor tissues
| DEGs | Gene names |
|---|
| Up-regulated | TPD52 SYT1///FRY PTPRN2 INSM1 SLC8A2 SLC8A3 NEFM RGS4 KCNH2 GPR19 MYCN GALNT18 PTPRN SCHIP1 SCG3 GABRB3 CHRNA3 SVOP CXXC4 GNG4 LINGO1 PSIP1 PLCD4 ATCAY PCSK2 PPP1R17 MAB21L1 SCN9A GLCCI1 |
| Down-regulated | MT1B PPP1R14A REEP6 RGN KCNK3 PTPRE GJA1 VEPH1 AKR1B10 RMDN2 ALDH3A2 DOCK8 APOE GATA6 MRAP ATP4A CYP21A2 SYT12 EFNA4 SLC6A6 TNXB TMEM200A IRS1 BAMBI CHST6 MT3 TMEM217 SLITRK4 CYP11B2 SCD5 SLC47A1 TBXA2R MCOLN2 IL1RL1 MGST1 MLIP KCNA4 ZNF275 SYTL5 LAMB3 AOC3 CYP3A5 FBLN1 SLC16A9 ABCC3 FGFR2 HSD3B1 FDX1 SULT2A1 P2RY2 REPS2 GALM CTH CALB2 ECHDC3 KLHL4 RND3 CPB1 TPD52L1 DHCR7 MT1X SOAT1 ALAS1 CCDC8 ANGPTL1 ICK VDR C4B DIRAS3 ASB4 PRKAR2B GRAMD1B CSDC2 GIPC2 MT1E RARRES2 CYP3A7 FGF9 SLC37A2 IGDCC4 EPHX1 ADGRL3 PPIF RBM47 RBBP7 MOCOS TBX3 AKR1B1 GSTA3 SNCG AMDHD1 POLE4 CYP11B1 MSMO1 MT1F FERMT1 PDGFRA IDH1 CCBE1 PGM5 CLDN1 GSTA5 TCEA3 MARCO US DDR2 TMEM61 LY6D C10orf82 SLIT2 BAALC MAPK4 AOX1 IL1R1 AMHR2 CELA2A APOC1 SLC51A AS3MT PSD3 WDR63 |
We conducted functional annotation of the 151 up- and down-regulated DEGs using the DAVID tool. In the BP category, up-regulated DEGs were enriched remarkably in the regulation of cell communication, while down-regulated DEGs were enriched in the cellular response to zinc ion, steroid metabolic process, C21-steroid hormone biosynthetic process, cholesterol metabolic and oxidation-reduction process, cellular response to cadmium ion, negative regulation of growth, and sterol metabolic process, among others (Table 2). CC category genes revealed a detectable enrichment of up-regulated DEGs in the regulation of synaptic vesicle membrane, cell junction, and integral component of membrane, while down-regulated DEGs were enriched in endoplasmic reticulum membrane, organelle membrane, integral component of plasma membrane, plasma membrane and perinuclear region of cytoplasm.
Table 2. Gene Ontology analysis of DEGs in PCPG
| Expression | Category | Term | Count | % | P-value | FDR |
|---|
| Up-regulated | GOTERM_BP_DIRECT | GO:0007154~cell communication | 2 | 0.041356493 | 2.90E-02 | 2.69E+01 |
| GOTERM_CC_DIRECT | GO:0030672~synaptic vesicle membrane | 3 | 0.062034739 | 4.54E-04 | 3.99E-01 |
| GOTERM_CC_DIRECT | GO:0030054~cell junction | 4 | 0.082712986 | 2.15E-03 | 1.88E+00 |
| GOTERM_CC_DIRECT | GO:0016021~integral component of membrane | 11 | 0.227460711 | 4.93E-02 | 3.59E+01 |
| GOTERM_MF_DIRECT | GO:0005432~calcium: sodium antiporter activity | 2 | 0.041356493 | 5.87E-03 | 4.80E+00 |
| Down-regulated | GOTERM_BP_DIRECT | GO:0071294~cellular response to zinc ion | 5 | 0.029936535 | 5.82E-06 | 9.04E-03 |
| GOTERM_BP_DIRECT | GO:0008202~steroid metabolic process | 6 | 0.035923841 | 7.94E-06 | 1.23E-02 |
| GOTERM_BP_DIRECT | GO:0006700~C21-steroid hormone biosynthetic
process | 4 | 0.023949228 | 2.11E-05 | 3.28E-02 |
| GOTERM_BP_DIRECT | GO:0008203~cholesterol metabolic process | 6 | 0.035923841 | 7.57E-05 | 1.18E-01 |
| GOTERM_BP_DIRECT | GO:0055114~oxidation-reduction process | 14 | 0.083822297 | 1.09E-04 | 1.69E-01 |
| GOTERM_BP_DIRECT | GO:0071276~cellular response to cadmium ion | 4 | 0.023949228 | 1.65E-04 | 2.56E-01 |
| GOTERM_BP_DIRECT | GO:0045926~negative regulation of growth | 4 | 0.023949228 | 2.33E-04 | 3.61E-01 |
| GOTERM_BP_DIRECT | GO:0006705~mineralocorticoid biosynthetic
process | 3 | 0.017961921 | 2.44E-04 | 3.78E-01 |
| GOTERM_BP_DIRECT | GO:0016125~sterol metabolic process | 4 | 0.023949228 | 3.16E-04 | 4.90E-01 |
| GOTERM_BP_DIRECT | GO:0006704~glucocorticoid biosynthetic process | 3 | 0.017961921 | 1.43E-03 | 2.20E+00 |
| GOTERM_BP_DIRECT | GO:0006805~xenobiotic metabolic process | 5 | 0.029936535 | 1.60E-03 | 2.46E+00 |
| GOTERM_BP_DIRECT | GO:0006749~glutathione metabolic process | 4 | 0.023949228 | 5.60E-03 | 8.35E+00 |
| GOTERM_BP_DIRECT | GO:0060045~positive regulation of cardiac
muscle cell proliferation | 3 | 0.017961921 | 8.71E-03 | 1.27E+01 |
| GOTERM_BP_DIRECT | GO:0033344~cholesterol efflux | 3 | 0.017961921 | 1.12E-02 | 1.60E+01 |
| GOTERM_BP_DIRECT | GO:0014066~regulation of phosphatidylinositol
3-kinase signaling | 4 | 0.023949228 | 1.39E-02 | 1.95E+01 |
| GOTERM_BP_DIRECT | GO:0050873~brown fat cell differentiation | 3 | 0.017961921 | 1.79E-02 | 2.45E+01 |
| GOTERM_BP_DIRECT | GO:0034651~cortisol biosynthetic process | 2 | 0.011974614 | 1.92E-02 | 2.60E+01 |
| GOTERM_BP_DIRECT | GO:0034447~very-low-density lipoprotein particle
clearance | 2 | 0.011974614 | 1.92E-02 | 2.60E+01 |
| GOTERM_BP_DIRECT | GO:0032342~aldosterone biosynthetic process | 2 | 0.011974614 | 1.92E-02 | 2.60E+01 |
| GOTERM_BP_DIRECT | GO:0046854~phosphatidylinositol phosphorylation | 4 | 0.023949228 | 2.27E-02 | 3.00E+01 |
| GOTERM_BP_DIRECT | GO:0007155~cell adhesion | 8 | 0.047898455 | 2.87E-02 | 3.64E+01 |
| GOTERM_BP_DIRECT | GO:0048015~phosphatidylinositol-mediated
signaling | 4 | 0.023949228 | 3.10E-02 | 3.86E+01 |
| GOTERM_BP_DIRECT | GO:0044550~secondary metabolite biosynthetic
process | 2 | 0.011974614 | 3.18E-02 | 3.94E+01 |
| GOTERM_BP_DIRECT | GO:0001701~in utero embryonic development | 5 | 0.029936535 | 3.26E-02 | 4.03E+01 |
| GOTERM_BP_DIRECT | GO:0032870~cellular response to hormone
stimulus | 3 | 0.017961921 | 3.39E-02 | 4.15E+01 |
| GOTERM_BP_DIRECT | GO:0034382~chylomicron remnant clearance | 2 | 0.011974614 | 3.80E-02 | 4.52E+01 |
| GOTERM_BP_DIRECT | GO:0002933~lipid hydroxylation | 2 | 0.011974614 | 3.80E-02 | 4.52E+01 |
| GOTERM_BP_DIRECT | GO:0036092~phosphatidylinositol-3-phosphate
biosynthetic process | 3 | 0.017961921 | 3.96E-02 | 4.66E+01 |
| GOTERM_BP_DIRECT | GO:0030308~negative regulation of cell growth | 4 | 0.023949228 | 4.32E-02 | 4.96E+01 |
| GOTERM_BP_DIRECT | GO:0048146~positive regulation of fibroblast
proliferation | 3 | 0.017961921 | 4.72E-02 | 5.28E+01 |
| GOTERM_CC_DIRECT | GO:0005789~endoplasmic reticulum membrane | 15 | 0.089809604 | 1.37E-03 | 1.61E+00 |
| GOTERM_CC_DIRECT | GO:0031090~organelle membrane | 4 | 0.023949228 | 1.89E-02 | 2.03E+01 |
| GOTERM_CC_DIRECT | GO:0005887~integral component of plasma
membrane | 17 | 0.101784217 | 2.02E-02 | 2.15E+01 |
| GOTERM_CC_DIRECT | GO:0005886~plasma membrane | 37 | 0.221530356 | 2.87E-02 | 2.93E+01 |
| GOTERM_CC_DIRECT | GO:0048471~perinuclear region of cytoplasm | 9 | 0.053885762 | 4.87E-02 | 4.47E+01 |
| GOTERM_MF_DIRECT | GO:0005506~iron ion binding | 9 | 0.053885762 | 5.79E-06 | 7.71E-03 |
| GOTERM_MF_DIRECT | GO:0016705~oxidoreductase activity, acting on
paired donors, with incorporation or reduction of
molecular oxygen | 4 | 0.023949228 | 5.65E-03 | 7.27E+00 |
| GOTERM_MF_DIRECT | GO:0004497~monooxygenase activity | 4 | 0.023949228 | 5.93E-03 | 7.62E+00 |
| GOTERM_MF_DIRECT | GO:0046934~phosphatidylinositol-4,5-bisphosphate
3-kinase activity | 4 | 0.023949228 | 7.14E-03 | 9.10E+00 |
| GOTERM_MF_DIRECT | GO:0016491~oxidoreductase activity | 6 | 0.035923841 | 8.91E-03 | 1.12E+01 |
| GOTERM_MF_DIRECT | GO:0020037~heme binding | 5 | 0.029936535 | 1.13E-02 | 1.40E+01 |
| GOTERM_MF_DIRECT | GO:0008395~steroid hydroxylase activity | 3 | 0.017961921 | 1.17E-02 | 1.45E+01 |
| GOTERM_MF_DIRECT | GO:0047783~corticosterone 18-monooxygenase
activity | 2 | 0.011974614 | 1.26E-02 | 1.56E+01 |
| GOTERM_MF_DIRECT | GO:0004507~steroid 11-beta-monooxygenase
activity | 2 | 0.011974614 | 1.26E-02 | 1.56E+01 |
| GOTERM_MF_DIRECT | GO:0042803~protein homodimerization activity | 11 | 0.065860376 | 1.74E-02 | 2.09E+01 |
| GOTERM_MF_DIRECT | GO:0008201~heparin binding | 5 | 0.029936535 | 1.89E-02 | 2.24E+01 |
| GOTERM_MF_DIRECT | GO:0004364~glutathione transferase activity | 3 | 0.017961921 | 2.07E-02 | 2.43E+01 |
| GOTERM_MF_DIRECT | GO:0004714~transmembrane receptor protein
tyrosine kinase activity | 3 | 0.017961921 | 2.41E-02 | 2.78E+01 |
| GOTERM_MF_DIRECT | GO:0016303~1-phosphatidylinositol-3-kinase
activity | 3 | 0.017961921 | 3.04E-02 | 3.37E+01 |
| GOTERM_MF_DIRECT | GO:0005088~Ras guanyl-nucleotide
exchange factor activity | 4 | 0.023949228 | 3.67E-02 | 3.92E+01 |
| GOTERM_MF_DIRECT | GO:0060228~phosphatidylcholine-sterol
O-acyltransferase activator activity | 2 | 0.011974614 | 3.74E-02 | 3.98E+01 |
| GOTERM_MF_DIRECT | GO:0046870~cadmium ion binding | 2 | 0.011974614 | 3.74E-02 | 3.98E+01 |
| GOTERM_MF_DIRECT | GO:0005509~calcium ion binding | 10 | 0.059873069 | 3.83E-02 | 4.06E+01 |
| GOTERM_MF_DIRECT | GO:0004908~interleukin-1 receptor activity | 2 | 0.011974614 | 4.35E-02 | 4.47E+01 |
In the MF category, up-regulated genes were enriched in calcium: sodium antiporter activity and down-regulated genes in iron ion binding, oxidoreductase activity, monooxygenase activity, phosphatidylinositol-4,5-bisphosphate 3-kinase activity, oxidoreductase activity, heme binding, steroid hydroxylase activity, corticosterone 18-monooxygenase activity, steroid 11-β-monooxygenase activity, protein homodimerization activity, heparin binding, glutathione transferase activity, etc. (Table 2).
KEGG pathway analysis indicated that DEGs in PCPG tumors were enriched in, for example, metabolic pathways, chemical carcinogenesis, steroid hormone biosynthesis, metabolism of xenobiotics by cytochrome P450, mineral absorption, drug metabolism - cytochrome P450 and glutathione metabolism (Fig. 2 and Supplementary Table S1).
Human PPI network analysis
We subjected the 151 DEGs to PPI multiple protein identifier analysis, which demonstrated 93 nodes and 139 edges (Fig. 3A). Fifty-eight DEGs were not incorporated into the PPI network analysis. Furthermore, the Cytoscape MCODE tool was applied for deeper analysis, which yielded 30 key genes among the 93 nodes (Fig. 3B).
Kaplan-Meier Plotter and GEPIA analysis of key genes
We used the Kaplan-Meier Plotter website to check the association with patient survival rate of these 30 key DEGs in PCPG tumors, and found that eight genes have a significant influence on survival. In contrast, the other 22 genes had no significant impact on survival (P > 0.05) (Table 3 and Fig. 4).
Table 3. Prognostic information for the 30 candidate genes
| Category | Genes |
|---|
| Genes with significant influence on survival (P < 0.05) | PCSK2 SCG3 ALDH3A2 AKR1B1 KCNH2 P2RY2 TBXA2R GNG4 |
| Genes without significant influence on survival (P > 0.05) | EPHX1 CYP3A5 MGST1 SULT2A1 HSD3B1 CYP21A2 FDX1 GSTA5 CYP11B2 CYP11B1 GSTA3 PTPRN2 PTPRN GALM AKR1B10 SCN9A IDH1 KCNK3 KCNA4 PDGFRA MYCN GJA1 |
Next, we applied GEPIA to compare the expression levels of these eight genes between 182 PCPG tumor patients and 3 normal people. Interestingly, we found that four of the eight were highly expressed in PCPG tumors. In comparison, the other four genes displayed low expression in PCPG tumors relative to normal adrenal tissue samples (Table 4 and Fig. 5).
Table 4. Validation of eight genes by GEPIA
| Category | Genes |
|---|
| Genes with high expression in PCPG (P < 0.05) | PCSK2 SCG3 TBXA2R GNG4 |
| Genes with low expression in PCPG (P < 0.05) | ALDH3A2 AKR1B1 KCNH2 P2RY2 |
To identify potential pathways involving these eight genes, we utilized KEGG pathway enrichment to analyze them with the DAVID tool (Table 5 and Fig. 6). We found that two key genes, ALDH3A2 and AKR1B1, were enriched in the glycerolipid metabolism pathway, but there was no significant KEGG pathway association for the other six genes. Since AKR1B1 and ALDH3A2 may influence the survival rate of PCPG cancer patients via glycerolipid metabolism, they therefore have greater potential than the other six genes to be biomarkers for PCPG tumors.
Table 5. Re-analysis of eight selected genes by KEGG pathway enrichment
| Category | Term | Count | % | P-value | Genes |
|---|
| KEGG_PATHWAY | myd00561: Glycerolipid metabolism | 2 | 0.16 | 0.029 | AKR1B1, ALDH3A2 |
TIICs are also an independent predictor of survival rate in many cancers. Therefore, we investigated whether the key genes revealed by our KEGG analysis correlate with PCPG-infiltrating lymphocytes using the TIMER tool and TCGA database. PCPG tumor purity was significantly positively correlated with KCNH2 expression but negatively correlated with TBXA2R expression. For the infiltrating lymphocytes, B cell infiltration level was predominantly positively correlated with expression of the ALDH3A2 and AKR1B1 genes, which were also enriched in the glycerolipid metabolism pathway. CD4+ T cell and dendritic cell infiltration were significantly positively correlated with TBXA2R expression. For the infiltration by macrophages, AKR1B1 and KCNH2 genes showed significant positive correlation. Intriguingly, all eight of the key DEGs showed weak or no significant correlations with neutrophil and CD8+ T cells in PCPG tumors (Fig. 7). For the correlation between SCNA and abundance of immune infiltrates in PCPG tumors, the infiltration by CD8+ T cells was increased by arm-level deletion of P2RY2 and decreased by arm-level deletion of ALDH3A2, arm-level gain of GNG4 and high amplification of TBXA2R. Further, the infiltration by CD4+ T cells was increased by arm-level deletion of PCSK2 and P2RY2. Infiltration by macrophages was increased by arm-level deletion of P2RY2 and SCG3, and decreased by high amplification of GNG4 (Fig. 8).
DISCUSSION
Aldehyde dehydrogenase 3 family member A2 (ALDH3A2), also known as fatty aldehyde dehydrogenase, has been shown to play an important role in lipid metabolism, including oxidation of fatty aldehyde, and lipid detoxification (Rizzo and Craft, 1991; González-Fernández et al., 2016). The ALDH3A2 gene consists of 10 exons that span approximately 30.5 kb and are located 50–85 kb from ALDH3 (Rogers et al., 1997). There are 19 different isozymes in the human ALDH family with similar functions (Marcato et al., 2011). They are biomarkers in cancer stem cells that show not only tumor features of aggressiveness and potential prognosis (Ma and Allan, 2011; Marcato et al., 2011), but also vital effects on the self-protection, differentiation and expansion of these cancer stem cells (Ma and Allan, 2011). Mutations in the ALDH3A2 gene in humans result in lipid metabolism errors leading to Sjögren–Larsson syndrome (Rizzo and Craft, 1991; Carney et al., 2004; Rizzo and Carney, 2005; Didona et al., 2007). There have been reports implicating ALDH3A2 especially in ovarian and prostate cancers. These reports indicate that ALDH3A2 is highly expressed in ovarian tumors and primary prostate cancers relative to normal tissues (van den Hoogen et al., 2010; Saw et al., 2012). On the other hand, ALDH3A2 was reported to be down-regulated in breast cancer, lung adenocarcinoma, lung squamous cell carcinoma, esophageal squamous cell carcinoma and HNSC, compared to normal tissues (Chang et al., 2018). These observations indicate that ALDH3A2 is an important gene in the progression of cancers.
Human aldose reductase (AKR1B1) plays an important role in glucose metabolism and lipid-aldehyde-mediated cell regulation (Barski et al., 2008). AKR1B1 also influences cell differentiation and growth by protecting against lipid peroxidation and steroidogenesis (Aida et al., 2000; Lefrançois-Martinez et al., 2004; Baba et al., 2009). AKR1B1 has been proposed to play varying roles in different types of cancers and is mostly related to abundant cell processes in tumors (Saraswat et al., 2006). Tammali et al. (2009) reported that AKR1B1 plays a role in colon cancer progression. Further, AKR1B1 was shown to be highly expressed in basal-like breast cancer, which is one of the breast cancers with the poorest prognosis (Wu et al., 2017), as well as in the development and progression of liver and lung cancers (Nakarai et al., 2015). Conversely, in prostate or endometrial cancers, AKR1B1 expression was significantly down-regulated relative to normal tissues (Laffin and Petrash, 2012; Hevir et al., 2013). Also, outer cortical fiber cells that expressed a high level of AKR1B1 appeared to show a dramatic change called the epithelial-to-mesenchymal transition (Yadav et al., 2011; Zablocki et al., 2011). AKR1B1 has also been shown to inhibit the migration of inflammatory cells (Yadav et al., 2007) and, hence, could increase the survival of tumor cells. AKR1B1 can also aid in the survival of tumor tissue because of its involvement in the regulation of epidermal growth factor, fibroblast growth factor, hypoxia-inducible factor-1α and vascular endothelial growth factor, which leads to angiogenesis (Tammali et al., 2011a, 2011b).
In addition, the tumor immune microenvironment is a factor that influences prognosis in patients with malignant tumors (Binnewies et al., 2018). Therefore, it has great potential to be a prognostic and effective biomarker for cancers (Taube et al., 2018). We report here that ALDH3A2 and AKR1B1 were associated with T cell and B cell infiltration in PCPG tumors. Although the infiltration by T and B lymphocytes in cancers may improve patient survival (Ruffini et al., 2009), studies in prostate cancer have suggested that the lack of TIICs predicted a high risk of progression and recurrence of tumors (Vesalainen et al., 1994). In contrast, intracapsular carcinoma correlated with weak expression of TIICs, while perineural invasion and capsular invasion correlated with strong expression of TIICs in prostate cancer (Vesalainen et al., 1994; Kärjä et al., 2005). The regulation of DEGs in TIICs is also an independent predictor of the prognosis of cancers in patients.
However, whether ALDH3A2 and AKR1B1 can serve as a prognostic marker for PCPG has not been investigated, and the correlation of both genes with tumor immune microenvironment in PCPG tumors has not been previously reported. Here, using bioinformatics analysis, we revealed ALDH3A2 and AKR1B1 as potential biomarkers, being down-regulated in PCPG tumors relative to normal adrenal tissues. The low expression of ALDH3A2 caused better prognosis, but the low expression of AKR1B1 resulted in poorer prognosis. Results from the Kaplan-Meier Plotter online tool showed that low expression of both genes could play key roles in the progression of PCPG and infiltration of immune cells into the tumor microenvironment. However, this needs further verification through a series of experiments. Nonetheless, data reported in this paper should provide useful information regarding potential biomarkers as well as insights into the biological mechanisms of malignant PCPG. As a limitation, some of the DEGs we screened may be from immune cells, but this should have little effect on the results because, in PCPG patients, not only the tumor cells but also the “complexity” of the PCPG tumor influences patient survival rates. We are now undertaking experiments to confirm our gene expression data in tumor cells and TIICs, to explore the mechanism by which AKR1B1 and ALDH3A2 influence the survival rate of PCPG patients.
CONCLUSIONS
Previous studies have confirmed that AKR1B1 and ALDH3A2 play important roles in the growth and development of diverse cancer types, but little research has been conducted on the relationship between AKR1B1, ALDH3A2 and PCPG tumors. Therefore, this paper reveals potentially useful biomarkers of PCPG tumors as well as new research directions for PCPG.
ACKNOWLEDGMENTS
We thank Zhenyu Dang (Fudan University), Tianchen Zhou (University of Iowa) and Jinrui Wen (Imperial College London) for their support and encouragement. This research was funded by the National Natural Science Foundation of China, grant number 81871577.
REFERENCES
- Aida, K., Ikegishi, Y., Chen, J., Tawata, M., Ito, S., Maeda, S., and Onaya, T. (2000) Disruption of aldose reductase gene (Akr1b1) causes defect in urinary concentrating ability and divalent cation homeostasis. Biochem. Biophys. Res. Commun. 277, 281–286.
- Baba, S. P., Wetzelberger, K., Hoetker, J. D., and Bhatnagar, A. (2009) Posttranslational glutathiolation of aldose reductase (AKR1B1): a possible mechanism of protein recovery from S-nitrosylation. Chem. Biol. Interact. 178, 250–258.
- Barski, O. A., Tipparaju, S. M., and Bhatnagar, A. (2008) The aldo-keto reductase superfamily and its role in drug metabolism and detoxification. Drug Metab. Rev. 40, 553–624.
- Binnewies, M., Roberts, E. W., Kersten, K., Chan, V., Fearon, D. F., Merad, M., Coussens, L. M., Gabrilovich, D. I., Ostrand-Rosenberg, S., Hedrick, C. C., et al. (2018) Understanding the tumor immune microenvironment (TIME) for effective therapy. Nat. Med. 24, 541–550.
- Burnichon, N., Vescovo, L., Amar, L., Libé, R., de Reynies, A., Venisse, A., Jouanno, E., Laurendeau, I., Parfait, B., Bertherat, J., et al. (2011) Integrative genomic analysis reveals somatic mutations in pheochromocytoma and paraganglioma. Hum. Mol. Genet. 20, 3974–3985.
- Carney, G., Wei, S., and Rizzo, W. B. (2004) Sjögren-Larsson syndrome: seven novel mutations in the fatty aldehyde dehydrogenase gene ALDH3A2. Hum. Mutat. 24, 186.
- Chang, P. M.-H., Chen, C.-H., Yeh, C.-C., Lu, H.-J., Liu, T.-T., Chen, M.-H., Liu, C.-Y., Wu, A. T. H., Yang, M.-H., Tai, S.-K., et al. (2018) Transcriptome analysis and prognosis of ALDH isoforms in human cancer. Sci. Rep. 8, 2713.
- Crona, J., Ljungström, V., Welin, S., Walz, M. K., Hellman, P., and Björklund, P. (2015) Bioinformatic challenges in clinical diagnostic application of targeted next generation sequencing: experience from pheochromocytoma. PLoS One 10, e0133210.
- Crona, J., Verdugo, A. D., Granberg, D., Welin, S., Stålberg, P., Hellman, P., and Björklund, P. (2013) Next-generation sequencing in the clinical genetic screening of patients with pheochromocytoma and paraganglioma. Endocr. Connect. 2, 104–111.
- Dahia, P. L. M. (2014) Pheochromocytoma and paraganglioma pathogenesis: learning from genetic heterogeneity. Nat. Rev. Cancer 14, 108–119.
- Didona, B., Codispoti, A., Bertini, E., Rizzo, W. B., Carney, G., Zambruno, G., Dionisi-Vici, C., Paradisi, M., Pedicelli, C., Melino, G., et al. (2007) Novel and recurrent ALDH3A2 mutations in Italian patients with Sjögren-Larsson syndrome. J. Hum. Genet. 52, 865–870.
- Feng, H., Gu, Z.-Y., Li, Q., Liu, Q.-H., Yang, X.-Y., and Zhang, J.-J. (2019) Identification of significant genes with poor prognosis in ovarian cancer via bioinformatical analysis. J. Ovarian Res. 12, 35.
- Ghosal, S., Das, S., Pang, Y., Gonzales, M. K., Huynh, T.-T., Yang, Y., Taieb, D., Crona, J., Shankavaram, U. T., and Pacak, K. (2020) Long intergenic noncoding RNA profiles of pheochromocytoma and paraganglioma: a novel prognostic biomarker. Int. J. Cancer 146, 2326–2335.
- González-Fernández, R., Hernández, J., Martín-Vasallo, P., Puopolo, M., Palumbo, A., and Ávila, J. (2016) Expression levels of the oxidative stress response gene ALDH3A2 in granulosa-lutein cells are related to female age and infertility diagnosis. Reprod. Sci. 23, 604–609.
- Hevir, N., Sinkovec, J., and Lanišnik Rižner, T. (2013) Decreased levels of AKR1B1 and AKR1B10 in cancerous endometrium compared to adjacent non-cancerous tissue. Chem. Biol. Interact. 202, 226–233.
- Kanehisa, M., Furumichi, M., Tanabe, M., Sato, Y., and Morishima, K. (2017) KEGG: new perspectives on genomes, pathways, diseases and drugs. Nucleic Acids Res. 45, D353–D361.
- Kärjä, V., Aaltomaa, S., Lipponen, P., Isotalo, T., Talja, M., and Mokka, R. (2005) Tumour-infiltrating lymphocytes: a prognostic factor of PSA-free survival in patients with local prostate carcinoma treated by radical prostatectomy. Anticancer Res. 25, 4435–4438.
- Laffin, B., and Petrash, J. M. (2012) Expression of the aldo-ketoreductases AKR1B1 and AKR1B10 in human cancers. Front. Pharmacol. 3, 104.
- Lefrançois-Martinez, A.-M., Bertherat, J., Val, P., Tournaire, C., Gallo-Payet, N., Hyndman, D., Veyssiére, G., Bertagna, X., Jean, C., and Martinez, A. (2004) Decreased expression of cyclic adenosine monophosphate-regulated aldose reductase (AKR1B1) is associated with malignancy in human sporadic adrenocortical tumors. J. Clin. Endocrinol. Metab. 89, 3010–3019.
- Ma, I., and Allan, A. L. (2011) The role of human aldehyde dehydrogenase in normal and cancer stem cells. Stem Cell Rev. Rep. 7, 292–306.
- Manger, W. M. (2006) An overview of pheochromocytoma: history, current concepts, vagaries, and diagnostic challenges. Ann. N. Y. Acad. Sci. 1073, 1–20.
- Marcato, P., Dean, C. A., Pan, D., Araslanova, R., Gillis, M., Joshi, M., Helyer, L., Pan, L., Leidal, A., Gujar, S., et al. (2011) Aldehyde dehydrogenase activity of breast cancer stem cells is primarily due to isoform ALDH1A3 and its expression is predictive of metastasis. Stem Cells 29, 32–45.
- Martucci, V. L., and Pacak, K. (2014) Pheochromocytoma and paraganglioma: diagnosis, genetics, management, and treatment. Curr. Probl. Cancer 38, 7–41.
- Nakarai, C., Osawa, K., Akiyama, M., Matsubara, N., Ikeuchi, H., Yamano, T., Hirota, S., Tomita, N., Usami, M., and Kido, Y. (2015) Expression of AKR1C3 and CNN3 as markers for detection of lymph node metastases in colorectal cancer. Clin. Exp. Med. 15, 333–341.
- Posa, I., Carvalho, S., Tavares, J., and Grosso, A. R. (2016) A pan-cancer analysis of MYC-PVT1 reveals CNV-unmediated deregulation and poor prognosis in renal carcinoma. Oncotarget 7, 47033–47041.
- Rizzo, W. B., and Carney, G. (2005) Sjögren-Larsson syndrome: diversity of mutations and polymorphisms in the fatty aldehyde dehydrogenase gene (ALDH3A2). Hum. Mutat. 26, 1–10.
- Rizzo, W. B., and Craft, D. A. (1991) Sjögren-Larsson syndrome. Deficient activity of the fatty aldehyde dehydrogenase component of fatty alcohol:NAD+ oxidoreductase in cultured fibroblasts. J. Clin. Invest. 88, 1643–1648.
- Rogers, G. R., Markova, N. G., De Laurenzi, V., Rizzo, W. B., and Compton, J. G. (1997) Genomic organization and expression of the human fatty aldehyde dehydrogenase gene (FALDH). Genomics 39, 127–135.
- Ruffini, E., Asioli, S., Filosso, P. L., Lyberis, P., Bruna, M. C., Macri, L., Daniele, L., and Oliaro, A. (2009) Clinical significance of tumor-infiltrating lymphocytes in lung neoplasms. Ann. Thorac. Surg. 87, 365–372.
- Saraswat, M., Mrudula, T., Kumar, P. U., Suneetha, A., Rao, T. S., Srinivasulu, M., and Reddy, G. B. (2006) Overexpression of aldose reductase in human cancer tissues. Med. Sci. Monit. 12, CR525–CR529.
- Saw, Y.-T., Yang, J., Ng, S.-K., Liu, S., Singh, S., Singh, M., Welch, W. R., Tsuda, H., Fong, W.-P., Thompson, D., et al. (2012) Characterization of aldehyde dehydrogenase isozymes in ovarian cancer tissues and sphere cultures. BMC Cancer 12, 329.
- Severgnini, M., Bicciato, S., Mangano, E., Scarlatti, F., Mezzelani, A., Mattioli, M., Ghidoni, R., Peano, C., Bonnal, R., Viti, F., et al. (2006) Strategies for comparing gene expression profiles from different microarray platforms: application to a case-control experiment. Anal. Biochem. 353, 43–56.
- Szklarczyk, D., Morris, J. H., Cook, H., Kuhn, M., Wyder, S., Simonovic, M., Santos, A., Doncheva, N. T., Roth, A., Bork, P., et al. (2017) The STRING database in 2017: quality-controlled protein-protein association networks, made broadly accessible. Nucleic Acids Res. 45, D362–D368.
- Tammali, R., Reddy, A. B. M., Ramana, K. V., Petrash, J. M., and Srivastava, S. K. (2009) Aldose reductase deficiency in mice prevents azoxymethane-induced colonic preneoplastic aberrant crypt foci formation. Carcinogenesis 30, 799–807.
- Tammali, R., Reddy, A. B. M., Saxena, A., Rychahou, P. G., Evers, B. M., Qiu, S., Awasthi, S., Ramana, K. V., and Srivastava, S. K. (2011a) Inhibition of aldose reductase prevents colon cancer metastasis. Carcinogenesis 32, 1259–1267.
- Tammali, R., Saxena, A., Srivastava, S. K., and Ramana, K. V. (2011b) Aldose reductase nhibition prevents hypoxia-induced increase in hypoxia-inducible factor-1α (HIF-1α) and vascular endothelial growth factor (VEGF) by regulating 26 S proteasome-mediated protein degradation in human colon cancer cells. J. Biol. Chem. 286, 24089–24100.
- Tang, Z., Li, C., Kang, B., Gao, G., Li, C., and Zhang, Z. (2017) GEPIA: a web server for cancer and normal gene expression profiling and interactive analyses. Nucleic Acids Res. 45, W98–W102.
- Taube, J. M., Galon, J., Sholl, L. M., Rodig, S. J., Cottrell, T. R., Giraldo, N. A., Baras, A. S., Patel, S. S., Anders, R. A., Rimm, D. L., et al. (2018) Implications of the tumor immune microenvironment for staging and therapeutics. Mod. Pathol. 31, 214–234.
- The Gene Ontology Consortium (2000) Gene ontology: tool for the unification of biology. Nat. Genet. 25, 25–29.
- van den Hoogen, C., van der Horst, G., Cheung, H., Buijs, J. T., Lippitt, J. M., Guzmán-Ramírez, N., Hamdy, F. C., Eaton, C. L., Thalmann, G. N., Cecchini, M. G., et al. (2010) High aldehyde dehydrogenase activity identifies tumor-initiating and metastasis-initiating cells in human prostate cancer. Cancer Res. 70, 5163–5173.
- Vesalainen, S., Lipponen, P., Talja, M., and Syrjänen, K. (1994) Histological grade, perineural infiltration, tumour-infiltrating lymphocytes and apoptosis as determinants of long-term prognosis in prostatic adenocarcinoma. Eur. J. Cancer 30A, 1797–1803.
- Wang, Z., Li, Y., Zhong, Y., Wang, Y., and Peng, M. (2020) Comprehensive analysis of aberrantly expressed competitive endogenous RNA network and identification of prognostic biomarkers in pheochromocytoma and paraganglioma. Onco Targets Ther. 13, 11377–11395.
- Wu, X., Li, X., Fu, Q., Cao, Q., Chen, X., Wang, M., Yu, J., Long, J., Yao, J., Liu, H., et al. (2017) AKR1B1 promotes basal-like breast cancer progression by a positive feedback loop that activates the EMT program. J. Exp. Med. 214, 1065–1079.
- Yadav, U. C. S., Shoeb, M., Srivastava, S. K., and Ramana, K. V. (2011) Aldose reductase deficiency protects from autoimmune- and endotoxin-induced uveitis in mice. Invest. Ophthalmol. Vis. Sci. 52, 8076–8085.
- Yadav, U. C. S., Srivastava, S. K., and Ramana, K. V. (2007) Aldose reductase inhibition prevents endotoxin-induced uveitis in rats. Invest. Ophthalmol. Vis. Sci. 48, 4634–4642.
- Zablocki, G. J., Ruzycki, P. A., Overturf, M. A., Palla, S., Reddy, G. B., and Petrash, J. M. (2011) Aldose reductase-mediated induction of epithelium-to-mesenchymal transition (EMT) in lens. Chem. Biol. Interact. 191, 351–356.
- Zou, B., Li, J., Xu, K., Liu, J.-L., Yuan, D.-Y., Meng, Z., and Zhang, B. (2019) Identification of key candidate genes and pathways in oral squamous cell carcinoma by integrated Bioinformatics analysis. Exp. Ther. Med. 17, 4089–4099.
