Identification of vital prognostic genes related to tumor microenvironment in pheochromocytoma and paraganglioma based on weighted gene co-expression network analysis

Pheochromocytoma and paraganglioma (PCPG) is a rare neuroendocrine tumor. This study aims to identify vital prognostic genes which were associated with PCPG tumor microenvironment (TME). We downloaded transcriptome data of PCPG from TCGA database and calculated the immune scores and stromal scores by using the ESTIMATE algorithm. DEGs related to TMB were then identified. We conducted WGCNA to further extract the TME-related modules. GO, KEGG pathway analysis, and PPI network were performed. Survival analysis was conducted to identify the hub genes associated with the prognosis of PCPG. A total of 150 PCPG samples were included in this study. We obtained 1507 and 2067 DEGs based on immune scores and stromal scores, respectively. WGCNA analysis identified the red module and brown module were correlated with immune sores while the turquoise module and red module were significantly associated with stromal scores. Functional enrichments analysis revealed that 307 TME-related genes were correlated with the inflammation or immune response. Survival analysis showed that three TME-relate genes (ADGRE1, CCL18, and LILRA6) were associated with PCPG prognosis. These three hub genes including ADGRE1, CCL18, and LILRA6 might be involved in the progression of PCPG and could serve as potential biomarkers and novel therapeutic targets.


INTRODUCTION
Pheochromocytoma and paraganglioma (PCPG) is a rare neuroendocrine tumor that produces catecholamines (CA) and originates from adrenal medulla or extra-adrenal ganglia [1]. It refers to pheochromocytoma as the neoplasm arises from the adrenal gland, whereas it is termed paraganglioma when the tumor originates from extra-adrenal tissues. Tumoral secretion of catecholamines accounts for characteristic clinical symptoms particularly episodic or sustained hypertension, palpitations, headache, and diaphoresis. Also, PCPG can lead to disorders in insulin metabolism, cardiovascular morbidity, and even mortality [2]. Most PCPGs are benign in their clinical presentation, however, have the potential risk of malignant transformation [3,4]. Patients with metastatic PCPG have a 5-year survival rate ranging from 40-77% [3]. At present, surgery is the main way to treat nonmetastatic PCPG, yet there has been no standardized therapeutic regimen for metastatic PCPG [4]. Further investigation into the molecular mechanism using AGING bioinformatics methods would provide novel insight into diagnosis and prognosis of PCPG.
Tumor microenvironment (TME) is defined as the cellular environment where the tumor cells are located, and is composed of different cell types including immune cells, stromal cells and endothelial cells etc [5]. In recent years, TME has been regarded as a critical factor in tumor progression and metastasis. However, the impact of TME on PCPG remains unclear. Immune cells and stromal cells are the major components of non-tumor cells in TME, and play an essential role in tumor diagnosing and prognostic predicting. Accumulated evidence has shown that tumor gene expression profile can quantify the immune infiltrating landscape of tumors.
The ESTIMATE (Estimation of STromal and Immune cells in MAlignant Tumor tissues using Expression data) algorithm proposed by Yoshihara et al. [6] has been applied to estimate the proportion of infiltrating stromal and immune cells in malignancy by using the expression profile. Previous studies have proved that the ESTIMATE algorithm was effective for calculating stromal scores, immune scores and tumor purity, thus providing a useful method for predicting the prognosis of patients using gene expression data [6,7]. However, there were no studies available exploring tumor microenvironment of PCPG utilizing ESTIMATE algorithm or high throughput data.
The weighted gene co-expression network analysis (WGCNA) has been widely used in biology and medical research and offers an effective approach to detect a group of genes with similar expression patterns as well as their relevant biological processes and pathways [8,9]. In this way, modules are defined as the clusters of highly correlated genes.
In the present study, we extracted the expression data of PCPG cohorts from TCGA database and then applied the ESTIMATE algorithm and the WGCNA method to perform a comprehensive analysis of tumor microenvironment-related genes for the first time. Our results indicated that the tumor microenvironment associated genes could affect the clinical prognosis of PCPG patients, suggesting that it might provide a basis for development of new prognostic biomarkers and therapeutic targets.

Evaluation of immune/stromal scores and DEGs screening
A total of 150 cases of PCPG from TCGA database were included in the present study and 3 normal cases were used. Based on the ESTIMATE algorithm, the stromal scores were ranged from -2002.755 to 1445.792, and the immune scores were ranged from -1503.612 to 2140.047. We divided all PCPG cases into high stromal score group and low stromal score group according to median value. Analogously, all PCPG patients were also classified into high immune score group and low immune score group according to median value. Based on stromal scores, 2067 DEGs related to stromal scores (including 1965 up-regulated DEGs and 102 down-regulated DEGs) were screened. Besides, 1507 DEGs related to immune sores (including 1442 up-regulated DEGs and 341 down-regulated DEGs) were obtained based on immune score ( Figure  1A, 1B). The Venn plots showed that there were a total of 1281 co-upregulation genes and 34 codownregulation genes ( Figure 1C, 1D).

Construction of weighted gene co-expression modules
The "WGCNA" package in R was applied to put the DEGs with similar expression patterns into modules by average linkage clustering. The best soft-thresholding parameter was set at 8 (scale-free R 2 = 0.85) to guarantee a scale-free network (Figure 2A-2D). A sample dendrogram was then constructed based on the similarity between the samples and the clinical characteristics of each sample ( Figure 3C, 3D). Finally, 7 gene co-expression modules of immune scores-related genes and 5 gene co-expression modules of stromal scores-related genes were identified, respectively ( Figure 3A, 3B). Module-trait relationship analysis showed that Module Eigengene (ME) of the red module and brown module were highly correlated with immune sores while the turquoise module and red module were significantly associated with stromal scores (Figures 2E-2H, 4D, 4E). In this way, we obtained 509 immune sores-related genes and 677 stromal scores-related genes from those key modules. 307 intersection genes were selected for further analysis ( Figure 4F).

Functional enrichment analyses
Functional enrichment analyses were conducted to investigate the biological processes and pathways relevant to these 307 genes using Metascape databases. The top 20 enriched biological processes and pathways were demonstrated in Figure 5 and Table 1, including regulation of cell activation, myeloid leukocyte activation, adaptive immune response, leukocyte migration, regulation of cytokine production, response to bacterium, Staphylococcus aureus infection, etc. The results showed that biological processes and pathways were mainly related to the inflammation or immune response.

PPI network, GO and KEGG analysis
All intersection genes were mapped into STRING database and the PPI networks were built using Cytoscape software. The plug-in MCODE screened four top core modules, which were shown in Figure 5A-5D. These PPI networks of core modules related to tumor microenvironment might be of great importance for PCPG development and progression. Enrichment analyses of GO and KEGG were conducted to assess the functions of the genes in the four identified modules, respectively. The most enriched GO and KEGG terms preserved in each module were displayed in Tables 2, 3, respectively. For the first module (module 1), positive regulation of cytokine production and cytokine receptor activity were the most significant enrichment in biological process and molecular functions, respectively. In the meantime, KEGG analysis revealed that genes in this module were enriched in complement and coagulation cascades pathway and Osteoclast differentiation pathway. For the second module (module 2), cellular response to biotic AGING stimulus, secretory granule membrane and cytokine receptor binding were the most significant enrichment in biological process, molecular function and cellular component, respectively. For the third module (module 3), the GO analysis indicated that these genes were mainly involved in extracellular matrix structural constituent. Protein digestion and absorption, ECMreceptor interaction, and PI3K-Akt signaling pathway were the most significantly enriched pathways. The genes in the fourth module (module 4) were primarily enriched in positive regulation of leukocyte cell-cell adhesion, Osteoclast differentiation pathway and Th17 cell differentiation pathway.

Survival analysis of intersection genes
Additional survival analysis was conducted on the 307 intersection genes to evaluate their effects on the survival of PCPG. Finally, three hub genes (ADGRE1, CCL18, and LILRA6) ( Figure 4) were identified to be associated with overall survival time (p < 0.05). We found that the high levels of CCL18 and LILRA6 were significantly correlated with poor survival outcomes, while the high levels of ADGRE1 was associated with longer overall survival.

Validation of expression levels of hub genes
We further validated the differential expression of ADGRE1, CCL18, and LILRA6 between PCPG tissue and normal tissue in TCGA database and ULCAN database. The results showed that only LILRA6 were differently expressed between tumor and normal tissues. The promoter methylation level of LILRA6 was also significantly lower tumor tissue compared with normal tissue. Besides, pan-cancers analysis showed that LILRA6 was not specific to PCPG. However, there were few molecular markers of PCPG at present and LILRA6 might be a vital marker in PCPG immune microenvironment. (Figure 6).

DISCUSSION
Increasing evidence indicates that the tumor microenvironment plays a critical role in tumor progression and metastasis [5,10]. Tumor development is highly related to the TME, and any alterations of tumor microenvironment may influence the clinical outcomes of malignancies. However, the effect of TME differs in different types of cancer. The correlation between TME and the PCPG prognosis remains poorly understood. In this study, we calculated the stromal scores and immune scores using the ESTIMATE algorithm to estimate the level of infiltrating stromal and immune cells in PCPG. A total of 2067 DEGs related to stromal scores and 1507 DEGs related to immune sores were obtained from comparison of lowand high-level groups. Then, the WCGNA analysis was conducted and 307 genes related to PCPG tumor microenvironment were identified. Finally, three tumor AGING microenvironment-related genes (ADGRE1, CCL18, and LILRA6) significantly associated with PCPG prognosis were obtained by survival analysis.
Functional enrichment analysis was conducted and PPI networks were built to further investigate the biological functions. The results revealed that the biological processes and pathways were mainly related to the inflammation or immune response. Likewise, the GO and KEGG analyses showed that the top GO terms and pathways enriched in the four core modules were also significantly associated with the inflammation and immune response. According to the previous studies, it was generally acknowledged that systemic inflam-mation played a central role in tumorigenesis and cancer progression. Several reports have shown that fever of unknown origin and systemic inflammatory syndrome were associated with IL-6 producing pheochromocytoma, which tended to have a larger tumor volume as well as an elevated risk of pheochromocytoma multisystem crisis [11,12]. Furthermore, the inflammation has long been considered to be closely associated with the pathogenesis of RCC [13]. A variety of systemic inflammatory indexes have been identified to exhibit a predictive value for RCC, including NLR, lymphocyteto-monocyte ratio and platelet count [14]. Karin et al. [15,16] reported that chronic inflammation caused by AGING prolonged infection with a bacterium, parasite, or virus was a main driving force in cancer development. It was found that the persistent inflammatory microenvironment set by HBV and HCV virus infections induced Hepatocellular carcinoma development [17][18][19]. Similarly, Chronic Helicobacter pylori infection increased the likelihood of mucosa-associated lymphoid tissue cancer and gastric cancer development. Cytokines were the key signaling molecules of communication between cells in the inflammatory tumor-microenvironment. The cytokines released during the persistent infection tended to induce several molecular signaling cascades, ultimately promoting neoplastic processes [20,21].

AGING
CCL18 is a C-C chemokine mainly secreted by M2 type tumor associated macrophages which acts mainly by binding to its corresponding chemokine receptor CCR8 [22]. Emerging evidence indicates that CCL18 serves numerous functions not only closely related to immune and inflammatory modulation, but also involved in cancer progression [23]. Wang et al. [24] reported that CCL18 could promote tumor angiogenesis, repress anticancer immune reaction and reshape TME, thus, leading to malignant progression in diverse human cancers. Moreover, another study revealed that CCL18 affected the replicative ability of tumor cells by promoting cell  transformation and altering the number of tumor cell chromosomes [25]. High expression CCL18 was tightly related to poor prognosis in various tumor, including ovarian cancer [26], pancreatic ductal adenocarcinoma [27] gallbladder carcinoma [28] and gastric cancer [29] etc. Additionally, CCL18 was expressed not only in cells of the immune system, but also in tumor cells. Liu et al. [22] reported for the first time that CCL18 was upregulated in bladder cancer (BC) cells, which further promote cell migration, invasion and epithelialmesenchymal transition (EMT). CCL18 overexpression has also been revealed to be associated with the proliferation, invasiveness, and angiogenesis of oral cancer cells, as well as the TNM stage [30]. However, there was no study investigating the role in PCPG. This study revealed that the high expression of CCL18 was associated with poor prognosis and might be a promising candidate gene affecting the occurrence and development of PCPG.

AGING
LILR, also known as the leukocyte Ig-like receptors, are a family of innate immune receptors that finely balance the functions of immune system and dictate their response to infected, stressed, and an aggressive tumor behavior [31]. Members of the LILR family played a vital role in various immune response processes such as cytokine production, DC maturation, and co-stimulatory molecules expression [32,33]. In a previous study, Jones et al [34]. reported that the paired activating and inhibitory LILRB3 and LILRA6 receptors could shape the local inflammatory responses in epithelial tumors by interacting with a ligand associated with the expression of cytokeratin 8. The balance of signals transmitted by paired LILRB3/LILRA6 receptors might therefore determine whether an immune response was triggered to the necrotic tumor cell, thereby leading to a wider immune response within the tumor microenvironment. In our study, the upregulation of LILRA6 was associated with poor prognosis in PCPG. Therefore, the development of specific LILRA6 inhibitors may provide an attractive target for anti-tumor immunotherapy in further study.
ADGRE1, named EMR1 in the past, is a member of adhesion G protein-coupled receptor family and has been regarded as a specific marker for eosinophils in humans [35]. However, the exact function of ADGRE1 is still unknown. Qi et al. [36] reported that ADGRE1, a down-regulated hub gene, was associated with the alpha fetoprotein (AFP) level in hepatocellular carcinoma (HCC) and may have crucial roles in HCC progression. Waddell et al [37] suggested the hypothesis that ADGRE1 could evolved in the process of immune selection and pathogen recognition, thus, participating in the host defense. Moreover, a previous study showed that afucosylated chimeric antibodies directed against AGING  ADGRE1 provided a promising target for immune cell ablation strategies [38]. In this study, we found that ADGRE1 was associated with longer overall survival in PCPG which indicated a protective role in PCPG biogenesis. However, the mechanism of ADGRE1 in PCPG development remains unknown and further research is required to be investigated.
In our study, we found that high expression levels of CCL18 and LILRA6 were significantly correlated with AGING poor survival, while the high levels of ADGRE1 were associated with longer overall survival. Interestingly, the results showed that all of the three hub genes belong to the 1281 upregulated genes. That is to say, the stromal and immune cells seem to be playing a dual role, either as the promoter or suppressor, during PCPG development and progression. To date, several studies have described the dual effect of the tumor stroma in the tumor-host interaction [39,40]. The stroma cells preferentially served an antitumor role in the tumor microenvironment. With ongoing tumor growth, stromal cells have been shown to promote growth and invasive behavior of tumor cells by secretion of growth factors, neovascularisation and facilitating recycling of anaerobic metabolic products [41]. Along with stromal cells, immune cells infiltrating the tumor tissue were also associated positively or negatively with tumor progression [42]. It was well known that CD8+ T cells were essential for immune defense and cytotoxic antitumor immunity. By contrast, CD4+ regulatory T cells suppressed anti-tumor immunity, thereby limiting the long-term efficacy of cancer immunotherapy [43]. Previous studies have reported that an imbalance of immune response may lead to oncogenesis and cancer progression [42]. Therefore, the immunomodulation may play multiple roles in the progression of PCPG and the specific mechanism remains an area of intense study.
There were several limitations in this study. First of all, due to the rather low incidence of PCPG, there were rare additional dataset of PCPG in several most commonly used database. We mainly explored the potential functions and hub genes of tumor immune microenvironment of PCPG based on WGCNA and TCGA database. Secondly, this was a retrospective study and the sample size was limited. Further prospective study is required. Finally, we just showed that these DEGs were correlated with ESTIMATE stromal or immune scores; however, we did not clarify whether these genes were expressed in immune or stromal cells. Whether these genes were expressed in immune or stromal cells or tumor cells required further cytological experiment.

CONCLUSIONS
In summary, by using a comprehensive bioinformatics analysis, we identified three tumor microenvironment- AGING related genes (ADGRE1, CCL18, and LILRA6) which were closely associated with the prognosis of PCPG.
These key genes may act as potential biomarkers and novel therapeutic targets. Further studies are required to analysis and validate the function of the identified genes and pathways in vitro and in vivo.

Data collection and preprocessing
The overall workflow of the present study was shown in Figure 7. We downloaded the level 3 gene expression profile, including 150 PCPG samples and 3 normal samples, from the cancer genome atlas (TCGA) database (https://portal.gdc.cancer.gov/). The immune scores and stromal scores were then calculated by using the ESTIMATE algorithm accordingly [6].

Differential gene expression analysis
Enrolled samples were divided into low-and high-level groups according to the median value of scores from the ESTIMATE algorithm. To reveal the correlations between the gene expression profiles and the stromal scores, "limma" package was applied to identify differentially expressed genes (DEGs). The selection criteria for DEGs were as follows: log |fold change (FC)| > 1 and adj. p <0.05. Heatmaps in this study were generated using "pheatmap" R package as described previously [44]. Likewise, we performed the same analysis procedure in the immune scores group.

Weighted gene co-expression network analysis (WGCNA)
After screening the differentially expressed genes (DEGs), WGCNA was conducted using the "WGCNA" package in R to identify the co-expression modules and significant genes that have similar patterns of expression to each other [45,8]. The dynamic tree cut method was employed to identify modules. Modules were named using different colors. Highly similar modules were identified by cluster analysis and the merged cutting height of 0.25 was set. Pearson correlation coefficients between modules and clinical traits were calculated, and modules significantly correlated with immune/stromal scores were selected (P < 0.05). The intersection genes in those key modules were then subjected to further analysis and the Venn diagrams were generated with the VennDiagram R package.

Functional enrichment analysis
Metascape (http://metascape.org/) is an online database providing a comprehensive set of functional annotation tools to determine the biological relevance behind large list of genes [46]. We uploaded the intersection genes to perform GO analysis and pathway enrichment analysis. P-value <0.05 was-considered statistically significant.

Construction of PPI network, GO analysis, and KEGG
All intersection genes were mapped to the Search Tool for the Retrieval of Interacting Genes (STRING) database, which is an online database for constructing protein-protein interaction (PPI) networks [47]. The MCODE plug-in in the Cytoscape software was used to screen the core modules [48]. To reveal the potential functions of significant genes in candidate modules, Gene Ontology (GO) term and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analyses were then performed using R software, with pvalue <0.05 as the cutoff value.

Survival analysis
Survival analysis was performed using the survival R package. The Kaplan-Meier plots was drawn to assess the relationship between intersect genes and overall survival, and the log-rank test was applied to test the significance of the difference between the two. P < 0.05 was considered to be significant.

Validation of expression levels of hub genes
UALCAN database (http://ualcan.path.uab.edu/) was used to validate the different expression levels, promoter methylation level of the hub genes between PCPG tissues and normal tissues. Additionally, the expression levels of hub genes were also compared with other tumor types to explore whether these genes were specific to PCPG.

Data availability statement
All data generated or analyzed during the present study was downloaded from the cancer genome atlas (TCGA) database and UALCAN database, and could be available freely from the corresponding author (Ning Xu).

AUTHOR CONTRIBUTIONS
Ning Xu and Yong Wei designed the work; Zhi-Bin Ke, Fei Lin, Hang Chen, and Xuan Tao carried out the experiments and analysed the data with the guidance of Ning Xu; Chun-Xian Chen, Dong-Ning Chen, Xiong-Lin Sun prepared the manuscript. All authors drafted or critically revised the manuscript for important intellectual content and approved the final version of the manuscript.

CONFLICTS OF INTEREST
The authors declare that they have no conflicts of interest.