Comprehensive Analysis of Aberrantly Expressed Competitive Endogenous RNA Network and Identification of Prognostic Biomarkers in Pheochromocytoma and Paraganglioma

Introduction

Pheochromocytoma and paraganglioma (PCPG) is a neuroendocrine tumor originating from chromaffin cells and sympathetic chains, characterized by the production of catecholamines.¹ Previous studies have shown that mutations in susceptibility genes of PCPG provide a basis for explaining the occurrence and development of PCPG.² Patients with PCPG are typically characterized by severe hypertension, irritability, sweating, and other symptoms, potentially leading to serious complications such as cerebral hemorrhage and heart disease. At present, the diagnosis of PCPG mainly depends on medical history, physical examination, biochemical tests, and imaging examinations. Nevertheless, it is difficult to diagnose and perform effective preoperative preparation in the absence of catecholamine secretion. In addition, malignant PCPG can only be diagnosed when a metastasis is found, resulting in a poor prognosis. Novel molecular markers are helpful to clearly distinguish between benign and malignant PCPG.^3,4 In order to confirm the diagnosis and characterization of PCPG, it is therefore necessary to explore the mechanisms of tumorigenesis and metastasis in PCPG.

Long non-coding RNA (lncRNA) is a non-coding RNA greater than 200 nucleotides in length. In recent years, studies have found that lncRNA is involved in many biological processes and plays an important role in tumorigenesis, progression, and metastasis.⁵ A few studies have examined lncRNA in PCPG.^6,7 However, the specific role of other lncRNAs in PCPG is unclear. With the introduction of the competitive endogenous RNA (ceRNA) hypothesis and the understanding of the lncRNA–miRNA–mRNA network regulation mechanism, it has been shown that lncRNA directly silences mRNA expression and indirectly acts as a “sponge” of microRNA (miRNA) to regulate mRNA. As a post-transcriptional product, lncRNA performs as an miRNA sponge to fully bind to miRNA, reducing the binding of miRNA to mRNA, which results in decreased mRNA silencing and increased mRNA expression.⁸ Through the study of ceRNA regulatory networks, we can further understand the role of lncRNA in PCPG and identify prognostic biomarkers for PCPG.

In this study, we collected differentially expressed lncRNA (DElncRNA), differentially expressed miRNA (DEmiRNA), and differentially expressed mRNA (DEmRNA) in PCPG tissues compared with adjacent adrenal gland medulla tissues, based on The Cancer Genome Atlas (TCGA) database. We found differentially expressed mRNAs involved in hormone secretion, metabolism signaling pathways, and tumorigenesis biological behaviors. Through integration analysis, we constructed an lncRNA–miRNA–mRNA network for ceRNA analysis in PCPG. Subsequently, we obtained two lncRNA–miRNA–mRNA axes in PCPG and performed Kaplan–Meier survival analysis to identify prognostic biomarkers. We further validated the mRNA and protein expression of prognostic biomarkers by the Gene Expression Omnibus (GEO) database and immunohistochemistry (IHC) staining, separately.

Methods

RNA-Seq, miRNA Data, and Clinical Data from TCGA-PCPG

R-packaged GDCRNATools were utilized in this study.⁹ RNA sequencing (RNA-seq), miRNA, and clinical features of TCGA-PCPG data were downloaded from TCGA (https://portal.gdc.cancer.gov/). The RNA-seq and miRNA data from different patients were integrated separately into two files. Duplicated samples were removed from the RNA-seq matrix and miRNA matrix and the count data of duplicated genes were averaged. Raw count data were normalized by the trimmed mean of M values (TMM) method implemented in edgeR and further transformed by the voom method provided in limma. Low-expression genes (logCPM <1 in more than half of the samples) were filtered out. Differentially expressed genes were subsequently calculated by the “edgeR” package by comparing PCPG with normal adrenal gland tissue.

Identification of Differentially Expressed lncRNA, miRNA, and mRNA

All long non-coding genes and protein coding genes were separately reported by setting the “gene Type” argument in “gdcDEReport”. Gene symbols and biotypes based on the Ensembl 38.90 annotation were reported in the output. DElncRNA and DEmRNA were identified by inclusion criteria. The inclusion criteria were “|log2(fold change) | >2, FDR <0.05”. The same procedure was performed with miRNA data to find DEmiRNA. Volcano maps and heatmaps of lncRNA, miRNA, and mRNA were then constructed.

GO and KEGG Pathway Analysis

A functional enrichment analysis of differentially expressed mRNAs was performed via the R/Bioconductor packages “clusterProfiler” and “DOSE”. Gene Ontology (GO) and the Kyoto Gene and Genomic Encyclopedia (KEGG) were analyzed and the bubble maps of the GO and KEGG pathways were drawn. After generating the ceRNA network, mRNAs in the ceRNA network were analyzed for GO and KEGG enrichment analysis via R packages.

Survival Analysis

Survival analysis of differentially expressed lncRNA, miRNA, and mRNA in PCPG was performed with GDCRNATools. RNAs with p-value <0.05 were identified as prognostic biomarkers.

Expression Level of Prognostic Biomarkers Between Metastatic and Non-Metastatic PCPG Groups

Based on the clinical file downloaded by the GDCRNATools package, 168 PCPG patients with metastatic status were enrolled in this study and the clinical data and normalized expression value of prognostic markers were extracted. The unpaired t-test was used for statistical analysis. GraphPad version 8.3.1 was used to display the results.

Construction of lncRNA–miRNA–mRNA ceRNA Regulatory Network

Based on RNA-seq and miRNA data, lncRNA–miRNA interactions and mRNA–miRNA interactions were analyzed to search for shared miRNAs. There are three criteria for determining the ceRNA between lncRNA and mRNA. First, lncRNA and mRNA share a certain number of miRNAs. Second, the correlation between lncRNA and mRNA expression is positive. Third, those shared miRNAs should have similar effects to regulate the expression levels of lncRNA and mRNA.⁹ The miRcode database was used to predict lncRNA–miRNA interactions, and mirTarBase database was used to predict mRNA–miRNA interactions. The correlation coefficient was measured to assess the linear association between lncRNA and mRNA. Modulation of similarity and sensitivity correlations was used to calculate the possible regulatory effects of consensus miRNAs on lncRNA and mRNA. hyperPValue <0.01 and corPValue <0.01 were set for analysis. Prognostic lncRNA, differentially expressed miRNA, and prognostic mRNA were selected to generate the “node” and “edge” files. Then, we imported the node and edge data into Cytoscape to build a ceRNA network.

Prediction of lncRNA–miRNA–mRNA Axis and Correlations Between PCPG Genotypes and the lncRNA–miRNA–mRNA Axis

It is well known that shared miRNAs are capable of binding lncRNAs and mRNAs. The more miRNAs bind to the lncRNAs, the fewer miRNAs will bind to the mRNAs. High mRNA expression is associated with more lncRNAs and fewer miRNAs. According to this trend, the lncRNA–miRNA–mRNA axis was identified from the ceRNA network. The correlations between lncRNAs and mRNAs were further calculated to confirm the regulatory axis.

Based on 17 susceptibility genes, patients with PCPG were divided into four genotypes: the pseudohypoxia group, the kinase signaling group, the Wnt signaling group, and the disease-modifying gene group.¹⁰ The expression of susceptibility genes and lncRNA–miRNA–mRNA axis genes in PCPG was extracted and Pearson correlation analysis was used to analyze the correlations between PCPG genotypes and lncRNA–miRNA–mRNA axis.

Validation of mRNA Expression and Survival in Pan-Cancer

A gene expression microarray data (GSE19422) was downloaded from Gene Expression Omnibus (GEO) database. The dataset contained 84 primary pheochromocytoma (PCC)/paraganglioma tumors and six normal adrenal gland tissues, and the platform of GPL6480-Agilent-014850 Whole Human Genome Microarray 4x44K G4112F (Probe Name version) was used for the analysis.¹¹ Gene symbols were annotated with the GPL6480 annotation file.

The mRNA expression of prognostic biomarkers involved in the lncRNA–miRNA–mRNA axis was plotted using TCGA-PCPG database and GSE19422 dataset. Overall survival and disease-free survival of these prognostic biomarkers in Pan-Cancer (excluding PCPG) were analyzed using the Gene Expression Profiling Interactive Analysis (GEPIA) database (http://gepia.cancer-pku.cn/).¹²

Immunohistochemistry (IHC)

Fifteen pairs of PCPG and adjacent adrenal gland tissues were collected from the Department of Urology, Second Xiangya Hospital of Central South University. At least two pathologists evaluated the HE staining of PCPG and adjacent adrenal gland tissues. Written informed consent was obtained from all patients and this study was approved by the Clinical Research Ethics Committee of the Second Xiangya Hospital of Central South University.

Patient samples harvested by surgical operation were fixed in neutral formalin and embedded in paraffin. Tissue slides were prepared and dewaxed by baking in an oven. The antigen retrieval was then carried out in a microwave oven with sodium citrate buffer pH 6.0 for 15 minutes, then the endogenous peroxidase was blocked by incubation in 3% hydrogen peroxide for 25 minutes, and then 3% BSA was used to block for 30 minutes at room temperature. Then, the sample was incubated with rabbit anti-human RCAN3 polyclonal antibody (1:200; Bioss, Beijing, China, #bs-11542R) or rabbit anti-human PTPRJ polyclonal antibody (1:200; Proteintech, Wuhan, China, #55123-1-ap) at 4°C overnight, followed by incubation with the secondary antibody, goat anti-rabbit antibody (1:200; Servicebio, Wuhan, China, #G23301), at room temperature for 50 minutes. Substrate and hematoxylin staining were used, and the slides were imaged by a microscope. The region of interest in the adjacent adrenal gland was the medulla. The IHC profiler was used to evaluate IHC score by ImageJ. The calculated scores were defined as: high positive = score 4, positive = score 3, low positive = score 2, and negative = score 1. Statistical analysis was performed using the paired t-test. The two groups were taken to have significant differences when p<0.05.

Results

Identification of Differentially Expressed lncRNA, miRNA, and mRNA

For a better understanding of the procedure used in this study, the workflow is presented in Figure 1. Altogether, 178 PCPG samples and three normal adrenal gland samples were included in the study. In total, 3496 mRNAs, 334 lncRNAs, and 116 miRNAs were obtained and met the inclusion criteria by comparing PCPG with normal adrenal gland tissue. The volcano plots and heatmaps are shown in Figure 2. The top 25 upregulated and top 25 downregulated mRNAs are shown in Table 1. Differentially expressed lncRNAs and differentially expressed miRNAs are shown in Tables 2, 3 and , respectively.

Table 1 Top 25 Downregulated and Top 25 Upregulated mRNA

Table 2 Top 25 Downregulated and Top 25 Upregulated lncRNA

Table 3 Top 25 Downregulated and Top 25 Upregulated miRNA

Table 4 Prognostic Biomarkers in Pheochromocytoma and Paraganglioma

Figure 1 Workflow of the study.

Figure 2 Differentially expressed miRNAs, mRNAs, and lncRNAs in tumor and normal adrenal gland tissues. (A) Volcano plots of miRNAs, mRNAs, and lncRNAs. Red: upregulated in tumor tissues. Green: upregulated in normal adrenal gland tissues. (B) Heatmap of differentially expressed lncRNAs (334 lncRNAs), differentially expressed mRNAs (3496 genes), and differentially expressed miRNAs (116 miRNAs).

GO and KEGG Pathways Involved in PCPG

To investigate the function of differentially expressed mRNAs, analysis of the GO and KEGG pathways was performed. The main enriched biological processes were “synapse organization” (Figure 3A). The main enriched cellular components were “cell body” (Figure 3B). The main enriched molecular functions were “tubulin binding” (Figure 3C) (Supplementary Table 1). Differentially expressed genes were enriched in 59 KEGG signaling pathways (Supplementary Table 2) and the top 10 KEGG pathways are displayed in Figure 3D. Although the major functions of these differentially expressed mRNAs in PCPG are related to the metabolism and hormone secretion signaling pathways, there were some metastasis-associated pathways, such as “ECM-receptor interaction”, “adherens junction”, and “cell adhesion molecules”; pathways associated with proliferation, such as “MAPK signaling pathway”; and cancer metabolism-associated pathways, such as “choline metabolism in cancer”, “central carbon metabolism in cancer”, and “AMPK signaling pathway”.

Figure 3 GO and KEGG pathways enriched in tumor tissues. (A) The top 10 enriched biological processes of differentially expressed genes. (B) The top 10 enriched cellular components of differentially expressed genes. (C) The top 10 enriched molecular functions of differentially expressed genes. (D) Bubble plot showing the top 10 enriched pathways.

Identification of Prognostic Biomarkers

To identify the prognostic biomarkers in PCPG, Kaplan–Meier curves were plotted for each differentially expressed mRNA, miRNA, and lncRNA. In total, 16 mRNAs, six lncRNAs, and eight miRNAs were associated with the prognosis of PCPG patients (Table 4). For lncRNAs, the Kaplan–Meier curves of overall survival indicated that high expression of AP006333.2, AP001486.2, and AL035071.1 was associated with good prognosis, whereas high expression of AL133163.3, AC156455.1, and AC008969.1 was associated with poor prognosis (Figure 4). For miRNAs, the Kaplan–Meier curves of overall survival in this ceRNA network revealed that high expression of hsa-miR-148b-3p and hsa-miR-338-3p was associated with good prognosis (Figure 5). For mRNAs, the Kaplan–Meier curves of overall survival showed that high expression of VPS37B, TTC38, RCAN3, ORAI2, NTNG1, NKX3-1, EPHA2, ENOSF1, and ACAA1 was associated with good prognosis, whereas high expression of PTPRJ, MED12, MCM10, KIAA1958, IREB2, CCNO, and ANKRD36 was associated with poor prognosis (Figure 6).

Figure 4 Kaplan–Meier survival curves for six differentially expressed long non-coding RNAs associated with overall survival in pheochromocytoma and paraganglioma. Six DElncRNAs are presented (p<0.05): AC008969.1, AC156455.1, AL035071.1, AL133163.3, AP001486.2, and AP006333.2. The initial pathological diagnosis of patients was set as day 0. Horizontal axis represents overall survival time (days) and vertical axis represents survival probability.

Figure 5 Kaplan–Meier survival curves for two differentially expressed miRNAs associated with overall survival in pheochromocytoma and paraganglioma. Two DEmiRNAs are presented (p<0.05): hsa-miR-148b-3p and hsa-miR-338-3p. The initial pathological diagnosis of patients was set as day 0. Horizontal axis represents overall survival time (days) and vertical axis represents survival probability.

Figure 6 Kaplan–Meier survival curves for 16 differentially expressed mRNAs associated with overall survival in pheochromocytoma and paraganglioma. Sixteen DEmRNAs are presented (p<0.05): ACAA1, ANKRD36, CCNO, PTPRJ, ENOSF1, IREB2, KIAA1958, RCAN3, MCM10, MED12, NKX3-1, TTC38, NTNG1, ORAI2, VPS378, and EPHA2. The initial pathological diagnosis of patients was set as day 0. Horizontal axis represents overall survival time (days) and vertical axis represents survival probability.

Different Expression of Prognostic Biomarkers Between Metastatic and Non-Metastatic PCPG Patients

By investigating the difference in expression of these prognostic biomarkers between metastatic patients and non-metastatic patients, we found that expression of hsa-miR-130b-3p and IREB2 was significantly higher in metastatic patients than in non-metastatic patients. In contrast, AP006333.2, hsa-miR-137, ACAA1, and NTNG1 were significantly more downregulated in metastatic patients compared with non-metastatic patients (Figure 7).

Figure 7 Different expression of prognostic biomarkers between metastatic and non-metastatic PCPG patients. (A) Expression of lncRNA between metastatic and non-metastatic PCPG patients. (B) Expression of miRNA between metastatic and non-metastatic PCPG patients. (C) Expression of mRNA between metastatic and non-metastatic PCPG patients.

Regulatory ceRNA Network and Two lncRNA–miRNA–mRNA Axes in PCPG

By exploring the interaction between prognostic lncRNA and differentially expressed miRNA, and the interaction between differentially expressed miRNA and prognostic mRNA, we found 20 shared miRNAs. The ceRNA regulatory network was based on 16 mRNAs, six lncRNAs, and 20 shared miRNAs (Figure 8A). Based on the difference of expression and correlation between lncRNA and mRNA, two lncRNA–miRNA–mRNA axes were identified: AP001486.2/hsa-miR-195-5p/RCAN3 and AP006333.2/hsa-miR-34a-5p/PTPRJ (Figure 8B). There were statistically significant positive correlations between AP001486.2 and RCAN3 (cor=0.46, p=1.4e-10) and between AP006333.2 and PTPRJ (cor=0.53, p=2.2e-14) (Figure 8C). There were significant negative correlations between AP001486.2 and hsa-miR-195-5p (cor=−0.26, p=0.00057), RCAN3 and hsa-miR-195-5p (cor=−0.19, p=0.01), AP006333.2 and hsa-miR-34a-5p (cor=−0.26, p=0.00057), and PTPRJ I and hsa-miR-34a-5p (cor=−0.23, p=0.0023).

Figure 8 ceRNA network and regulatory axes. (A) A regulatory lncRNA–miRNA–mRNA ceRNA network in PCPG patients. Red: downregulated. Green: upregulated. Diamond shape: lncRNA. Rounded rectangle: miRNA. Triangle: mRNA. ceRNA: competitive endogenous RNA; PCPG: pheochromocytoma and paraganglioma. (B) Two lncRNA–miRNA–mRNA axes in the ceRNA network. (C) The expression of AP001486.2 was positively correlated with RCAN3 (R=0.46, p=1.4×10−10), whereas the expression of AP006333.2 was positively correlated with PTPRJ (R=0.53, p=2.2×10−14). The expression of hsa-miR-195-5p was negatively correlated with RCAN3 (R=−0.19, p=0.01) and AP001486.2 (R=−0.26, p=0.00057), whereas the expression of hsa-miR-34a-5p was negatively correlated with PTPRJ (R=−0.23, p=0.0023) and AP006333.2 (R=−0.26, p=0.00057).

Correlations Between PCPG Genotypes and lncRNA–miRNA–mRNA Axes

In order to investigate the correlations between susceptibility genes and lncRNA–miRNA–mRNA axes, we extracted the expression value of 20 susceptibility genes and the lncRNA–miRNA–mRNA axes. We found that RET highly correlated with AP006333.2 in PCPG samples (R=0.66, p<2.2e-16) (Figure 9).

Figure 9 Correlations between PCPG genotypes and lncRNA–miRNA–mRNA axes. (A) Correlation heatmap of susceptibility genes and lncRNA–miRNA–mRNA axes genes. (B) Correlation plot between RET and AP006333.2.

Enrichment Analysis Predicting the ceRNA Network

To predict the function of the ceRNA network, 16 mRNAs were involved in enrichment analysis of GO and KEGG pathways. The top 5 enriched biological processes (protein kinase B signaling, post-anal tail morphogenesis, positive regulation of cell junction assembly, negative regulation of protein kinase B signaling, and somitogenesis) are shown in Figure 10A. The top 5 enriched cellular components (ruffle membrane, specific granule, leading edge membrane, ruffle, and anchored component of synaptic membrane) are shown in Figure 10B. The top 5 enriched molecular functions (hydro-lyase activity, carbon-oxygen lyase activity, beta-catenin binding, nuclear hormone receptor binding, and lyase activity) are shown in Figure 10C. The top 5 enriched KEGG pathways (axon guidance, alpha-linolenic acid metabolism, biosynthesis of unsaturated fatty acids, fructose and mannose metabolism, and fatty acid degradation) are shown in Figure 10D. These results show that mRNAs in the ceRNA network were enriched in metabolism-related pathways.

Figure 10 Enriched GO and KEGG pathways in the ceRNA network. (A) The top 5 enriched biological processes of the ceRNA network. (B) The top 5 enriched cellular components of the ceRNA network. (C) The top 5 enriched molecular functions of the ceRNA network. (D) Bubble plot showing the top 5 enriched pathways of the ceRNA network.

Validation of Expression Levels of RCAN3 and PTPRJ and Prognostic Value in Pan-Cancer

Based on TCGA database, we found that RCAN3 and PTRPRJ expression levels in PCPG were higher than in normal adrenal gland tissue (RCAN3: p=0.0021; PTPRJ: p=0.00031). To validate the mRNA expression of RCAN3 and PTPRJ, we found a similar expression pattern of RCAN3 and PTPRJ in pheochromocytoma/paraganglioma primary tumors based on the GSE19422 dataset (RCAN3: p=0.00016; PTPRJ: p=0.014) (Figure 11). The results of IHC further indicated higher protein expression levels of RCAN3 and PTPRJ in pheochromocytoma/paraganglioma primary tumors compared with adjacent adrenal gland medulla tissues (RCAN3: p=0.0406, PTPRJ: p=0.0192) (Figure 12). To further identify the prognostic value of RCAN3 and PTPRJ, Kaplan–Meier curves of overall survival and disease-free survival in Pan-Cancer were plotted, and the high expression levels of RCAN3 and PTPRJ were associated with a favorable prognosis (Figure 13).

Figure 11 PTPRJ and RCAN3 mRNA expression in TCGA and GSE19422 dataset. (A) PTPRJ and RCAN3 mRNA expression levels in PCPG were compared to normal adrenal gland tissue samples based on TCGA-PCPG database. (B) PTPRJ and RCAN3 mRNA expression levels in PCPG were compared to normal adrenal gland tissue levels based on the GSE19422 dataset.

Figure 12 Immunohistochemistry (IHC) analysis of protein expression of RCAN3 and PTPRJ in adjacent adrenal gland and PCPG tissues. (A) IHC results of RCAN3 and PTPRJ protein expression in adjacent adrenal gland and PCPG tissues. (B) Histogram of protein expression of RCAN3 and PTPRJ indicating that they were both upregulated in PCPG tissue compared with adjacent adrenal gland tissues. Scale bar, 50 µm.

Figure 13 Overall survival and disease-free survival of RCAN3 and PTPRJ mRNA expression in Pan-Cancer.

Discussion

Previous studies have indicated that the mechanism of post-transcriptional regulation is complicated. Many factors, such as mRNA, miRNA, lncRNA, and circular RNA (circRNA), are involved. In our study, we integrally analyzed RNA-seq data and miRNA data of TCGA-PCPG and constructed a regulatory lncRNA–miRNA–mRNA network. Based on the ceRNA network, six lncRNAs, two miRNAs, and 16 mRNAs were identified as prognostic biomarkers in PCPG.

Several studies have identified differentially expressed genes of PCPG, which were different from those in our study. One previous study compared malignant PCPG with the benign samples from TCGA and COMETE databases to investigate malignancy-induced differentially expressed genes. It found 11 signaling pathways that were related to metabolic alteration and cancer signaling.¹³ Among these signaling pathways, dopamine signaling and cAMP signaling were also reported in our study. Another study used 183 pheochromocytoma samples and three control samples to identify differentially expressed RNAs that were annotated by conversion through the Ensembl file (Homo sapiens GRCh38.87.chr.gtf). In that study, the threshold was adjusted by p-value <0.01 and |log2 fold change (FC)| >2.0.¹⁴ Circadian entrainment, cholinergic synapse, and aldosterone synthesis and secretion were considered as the paramount pathways in their study. In our study, we included 178 PCPG samples and three normal adrenal gland samples from TCGA database and used the threshold of p<0.05 and |log2(fold change) | >2 to identify differentially expressed RNAs (lncRNA, miRNA, and mRNA) which were annotated within Ensembl 90. These identified differentially expressed RNAs were enriched in some KEGG signaling pathways that were related to metabolism alteration, hormone changes produced by the adrenal gland and tumorigenesis biological behavior. Some KEGG signaling pathways overlapped with study by Liang et al,¹⁴ such as glutamatergic synapse, circadian entrainment, cholinergic synapse, and aldosterone synthesis and secretion. Furthermore, some metastasis-associated pathways (such as ECM–receptor interaction, adherens junction, and cell adhesion molecules), proliferation-associated pathway (such as the MAPK signaling pathway), and cancer metabolism-associated pathways (such as choline metabolism in cancer, central carbon metabolism in cancer, and the AMPK signaling pathway) were mentioned in our study.

Functional, pathogenesis-based classification deepens our understanding of PCPG. The mutation of susceptibility genes dysregulates biological functions in PCPG. Based on this functional classification of susceptibility genes, PCPG patients can be divided into two clusters (cluster 1 and cluster 2)¹⁵ and further divided into four major groups: the pseudohypoxia group (VHL, EGLN1, EPAS1, SDHx, SDHAF2, FH), the kinase signaling group (RET, HRAS, TMEM127, MAX), the Wnt signaling group (CSDE1, MAML3), and the disease-modifying gene group (ATRX, KMT2D, SETD2, TERT, TP53).¹⁰ In our study, we found SDHA, SDHB, and RET as differentially expressed genes and the GO analysis revealed that the tricarboxylic acid cycle enzyme complex is involved in cellular components. This suggested that most of the patients with PCPG in TCGA database were in the pseudohypoxia group, in accordance with the previous study’s results.¹⁵

Studies on lncRNA in PCPG are emerging and several lncRNAs have been found to play a vital role in tumorigenesis and malignancy. Limited studies found that SOX2OT, EFNA3, and MALAT-1 were involved in injury induced by dihydrogen peroxide,¹⁶ hypoxia,¹⁷ and chemical hypoxia in PC-12 cells,¹⁸ respectively. A previous study indicated that LINC00472, RP11-455b3.1, and RP4-806M20.3 were associated with a poor clinical outcome in PCPG and the three lncRNAs served as prognostic biomarkers.¹⁹ More recently, another group identified C9orf147 and BSN-AS2 as independent prognostic lncRNAs for PCPG.¹⁴ In the present study, we identified that high expression of AC008969.1, AC156455.1, and AL133163.3 predicted a good prognosis in PCPG, while high expression of AP006333.2, AP001486.2, and AL035071.1 was associated with a poor prognosis. These six lncRNAs had not been reported in previous studies and needed to be further studied for their function and underlying mechanisms in PCPG.

Many studies have reported that miRNAs act as prognostic biomarkers in PCPG, such as miR-183²⁰ and miR-483-5p.²¹ In our study, we found that high expression of two miRNAs, hsa-miR-148b-3p and hsa-miR-338-3p, was associated with a good prognosis in PCPG. Previous studies reported that hsa-miR-338-3p enhanced the risk of esophageal squamous cell carcinoma²² and decreased migration and invasion of colorectal carcinoma cells by targeting smoothened.²³ miR-148b-3p was also associated with the osteogenic differentiation capability of BM-MSCs.²⁴ These studies suggested various roles of miR-338-3p and miR-148b-3p in different diseases. In addition, we predicted that hsa-miR-195-5p and hsa-miR-34a-5p were involved in two lncRNA–miRNA–mRNA axes. hsa-miR-195 was also involved in a ceRNA network constructed by Liang et al¹⁴ and hsa-miR-195-5p was an anti-oncogene to inhibit cell proliferation and angiogenesis.^25,26 hsa-miR-34a-5p was reported to induce apoptosis and cell-cycle arrest, and suppressed tumor invasion and migration.²⁷

A number of mRNAs have been found as prognostic biomarkers in cancers.²⁸ In our study, mRNAs in the ceRNA network were enriched in metabolism-related signaling pathways. Two genes (RCAN3, PTPRJ) are important in two predicted lncRNA–miRNA–mRNA axes. RCAN3 was a favorable prognostic biomarker while PTPRJ was an unfavorable marker in PCPG. RCAN3 not only may be involved in immune response as an endogenous regulator of the calcineurin–nuclear factor of activated T cells (NFAT) signaling pathway,²⁹ but also plays a vital role in suppression of cancer.³⁰ PTPRJ, also named CD148, could negatively regulate T-cell activation and play a role in feedback inhibition of T-cell immune responses,³¹ which may result in a poor prognosis in PCPG. However, in malignant cancers, PRPTJ acted as a tumor suppressor to inhibit the proliferation of cancer,³² suggesting a potential different mechanism by which PRPTJ serves as a favorable biomarker in Pan-Cancer.

There are some limitations of this study. First, the ceRNA network produced by our inclusion and exclusion criteria was not static and the updated databases may lead to more interactions of lncRNA–miRNA and mRNA–miRNA and expend the ceRNA network. Second, more cohorts of PCPG are required to validate novel prognostic biomarkers.

In summary, we successfully constructed an lncRNA–miRNA–mRNA regulatory ceRNA network in PCPG and identified 16 mRNAs, six lncRNAs, and two miRNAs as prognostic biomarkers. AP001486.2/hsa-miR-195-5p/RCAN3 and AP006333.2/hsa-miR-34a-5p/PTPRJ axes were predicted. These findings could contribute toward a better understanding of the biological mechanism of tumorigenesis of PCPG.