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

Background Long non-coding RNA (lncRNA) functions as a competitive endogenous RNA (ceRNA) and plays an important role in the biological processes underlying tumorigenesis. However, studies describing the function of lncRNA in pheochromocytoma and paraganglioma (PCPG) remain largely unknown. Our study aims to construct a regulatory ceRNA network and explore prognostic biomarkers for PCPG through a comprehensive analysis. Methods PCPG data from The Cancer Genome Atlas (TCGA) were utilized to obtain differentially expressed lncRNAs (DElncRNAs), microRNAs (DEmiRNAs), and mRNAs (DEmRNAs). Kaplan–Meier analysis was used to detect prognostic biomarkers and Cytoscape was utilized to construct a regulatory network of ceRNA. Potential lncRNA–miRNA–mRNA axes were inferred by correlation analysis. GO and KEGG pathways were constructed using “clusterProfiler” and “DOSE” R-packages. Immunohistochemistry (IHC) staining was performed to validate differential protein expression levels of genes in the axes. Finally, the GSE19422 dataset and Pan-Cancer data were applied to validate the expression pattern and survival status of mRNAs, respectively. Results A total of 334 DElncRNAs, 116 DEmiRNAs, and 3496 DEmRNAs were identified and mainly enriched in hormone secretion, metabolism signaling, metastatic and proliferative pathways. Among these differentially expressed genes, 16 mRNAs, six lncRNAs, and two miRNAs were associated with overall survival of patients with PCPG and sequentially enrolled in the ceRNA network. Two lncRNA–miRNA–mRNA regulatory axes were predicted: AP001486.2/hsa-miR-195-5p/RCAN3 and AP006333.2/hsa-miR-34a-5p/PTPRJ. The GSE19422 dataset and IHC analysis validated that mRNA and protein levels of RCAN3 and PTPRJ were upregulated in PCPG tissues compared with adjacent adrenal gland medulla tissues. Pan-Cancer data showed that the upregulated expression of RCAN3 and PTPRJ was associated with favorable overall survival and disease-free survival. Conclusion A regulatory lncRNA–miRNA–mRNA ceRNA network was successfully constructed and 24 prognostic biomarkers were identified for PCPG patients. These findings may contribute toward a better understanding of the biological mechanism of tumorigenesis and enable further evaluation of the prognosis of patients with PCPG.


Introduction
Pheochromocytoma and paraganglioma (PCPG) is a neuroendocrine tumor originating from chromaffin cells and sympathetic chains, characterized by the production of catecholamines. 1 Previous studies have shown that mutations in susceptibility genes of PCPG provide a basis for explaining the occurrence and development of PCPG. 2 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. 5 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. 8 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. 9 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. Lowexpression 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. submit your manuscript | www.dovepress.com

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. 9 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

DovePress
OncoTargets and Therapy 2020:13 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. 10 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. 11 Gene symbols were annotated with the GPL6480 annotation file.
The mRNA expression of prognostic biomarkers involved in the lncRNA-miRNA-mRNA axis was plotted 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 antihuman 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.

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  DovePress functions of these differentially expressed mRNAs in PCPG are related to the metabolism and hormone secretion signaling pathways, there were some metastasisassociated 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".

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).

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 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).

Discussion
Previous studies have indicated that the mechanism of posttranscriptional 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. 13 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. 14 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

DovePress
OncoTargets and Therapy 2020:13 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, 14 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) 15 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). 10 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. 15 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, 16 hypoxia, 17 and chemical hypoxia in PC-12 cells, 18 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. 19 More recently, another group identified C9orf147 and BSN-AS2 as independent prognostic lncRNAs for PCPG. 14 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 20 and miR-483-5p. 21 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 22 and decreased migration and invasion of colorectal carcinoma cells by targeting smoothened. 23 miR-148b-3p was also associated with the osteogenic differentiation capability of BM-MSCs. 24 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 14 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. 27 A number of mRNAs have been found as prognostic biomarkers in cancers. 28 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, 29 but also plays a vital role in suppression of cancer. 30 PTPRJ, also named CD148, could negatively regulate T-cell activation and play a role in feedback inhibition of T-cell immune responses, 31 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, 32 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.