PARGP1, a specific enhancer RNA associated with biochemical recurrence of prostate cancer

Abstract

Background: Enhancer RNA (eRNA) is a non-coding RNA transcribed by enhancers and has a crucial role in controlling gene transcription. However, the potential functions of eRNA in prostate cancer (PCa) are still not fully understood.

Methods: Clinical information and expression data of eRNAs and their target genes for PCa were downloaded from TCGA. Each eRNA in PCa samples was divided into high- and low-expression groups. The Kaplan–Meier method was used to identify the eRNAs significantly associated with the overall survival of PCa. Their correlations with age, prostate-specific antigen (PSA), and tumor stage were analyzed using Spearman coefficients. Finally, Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEEG) analyzes were used to predict the possible biological functions and pathways associated with eRNAs.

Results: We identified 12 eRNAs associated with PCa survival, of which PARGP1 was the most clinically significant. It was associated with biochemical recurrence in PCa patients and was also a valuable complement to distinguish between PSA 0–4 ng/ml and 4–10 ng/ml. GO, and KEEG analysis showed that herpes simplex virus-1 infection was the most biologically significant enriched process.

Conclusions: PARGP1 is associated with biochemical recurrence in prostate cancer and contributes to assessing the outcome.

KEYWORDS:

• PARGP1
• enhancer RNA
• prostate cancer
• prognosis
• biochemical recurrence

Introduction

Prostate cancer (PCa) is one of the most common malignant tumors worldwide and ranks first among malignant tumors in men in the United States (Miller et al. 2019). In China, the morbidity of PCa ranks seventh among male cancers (Liu et al. 2019). According to the Global Burden of Disease 2017, there were 144,887 cases of PCa and 51,718 deaths in China (Pang et al. 2016). When PCa grows rapidly or metastasizes, the disease worsens, and the late metastasis and mortality rates are higher (Sebesta and Anderson 2017; Abdollah et al. 2018). Prostate-specific antigen (PSA) is a traditional screening parameter. However, PSA has low specificity, especially when the total PSA value is 4–10 ng/ml; therefore, there is an urgent need to identify new reliable diagnostic and prognostic biomarkers or therapeutic targets for PCa (Nguyen-Nielsen and Borre 2016; Eapen et al. 2019; Eskra et al. 2019).

Enhancer RNAs (eRNAs) are transcribed and produced by enhancers, 500–5 kb in length (Zhang et al. 2019). These regulatory molecules stimulate the activity of critical epigenetic regulatory enzymes, promoting gene expression. Numerous studies showed that eRNAs could be used as new targets for tumor therapy (Puc et al. 2015; Zhao et al. 2016; Bose and Berger 2017).

There is a paucity of research on eRNA in PCa. Therefore, in this study, we screened eRNAs and the target genes related to outcomes PCa using The Cancer Genome Atlas (TCGA) and performed clinical correlation analysis and functional enrichment to reveal their potential biological functions and signal pathways related to PCa. Our ultimate aim is to improve outcome evaluation and personalized treatment.

Materials and methods

Data download and collation

Fragments Per Kilobase Million (FPKM) files containing eRNA sequencing data, clinical data, and survival data for PCa and 32 other tumors from TCGA were downloaded from the UCSC Cancer Browser (https://xena.ucsc.edu/); the PreSTIGE algorithm was used to predict the target genes of the eRNAs, and downloaded the RNA-Seq of pca from the TCGA, screening the common genes to both for further analysis (Corradin et al. 2014; Vučićević et al. 2015). If several probes detected the same eRNA, the average value was taken as the eRNA expression value. The eRNA expression matrix was combined with patient survival data using the limma software package.

Screening survival-related eRNA

The eRNAs were divided into high and low expression groups based on the median value of expression. Kaplan–Meier survival curves were used to compare the high- and low-expression groups’ overall survival rates. The log-rank test was used to determine whether there were differences between groups for overall survival. Statistical significant differences were assumed when the p-value was < 0.05.

The Spearman method was used to identify eRNA and target gene pairs with correlation coefficients r > 0.04 and p < 0.001 in Pca (Gu et al. 2019). ERNAs associated with overall survival (p < 0.05) and their target gene levels (p < 0.001) were considered key eRNAs.

Clinical correlation analysis

The R software limma package was used to analyze the correlations between survival-related eRNA and target genes and biochemical recurrence status, age, Gleason score, pathological T stage, lymphatic metastasis status, and distant metastasis of PCa patients. The ggpubr package was used to draw images.

Functional annotation and enrichment analysis

To further understand the biological information of eRNA, we analyzed co-expression with other genes in PCa. The screening condition was that the correlation coefficient would be more than 0.4 and p < 0.001. The ‘ClusterProfiler’ package in R software was used for Gene Ontology (GO) functional enrichment analysis and Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment analysis.

Statistical analysis

The eRNAs were divided into high and low expression groups based on the median value of expression. The overall survival for the study population was generated using the Kaplan–Meier method, and differences in survival were validated using the log-rank test. The Spearman method was used to identify the correlation coefficient. All statistical analyzes were performed in R, version 4.0.0 (http://www.r-project.org/) and the accompanying software package.

Results

ERNAs related to PCa outcome

Based on RNA sequencing data from PCa patients downloaded from TCGA, there were 1580 eRNAs in PCa, 37 of which were significantly associated with overall survival (Kaplan–Meier log-rank test, p < 0.05; Table 1). The PreSTIGE algorithm was used to predict the target genes for these 37 eRNAs. When the levels of these 37 eRNAs were compared with the levels of their predicted target gene mRNAs, only 12 were found to be significantly correlated (Spearman’s rank correlation coefficient r > 0.4, p < 0.001, Table 1).

Table 1. List of overall survival associated eRNAs derived from enhancers.

			Correlation between eRNA and the Neighboring Target
Symbol	Kaplan–Meier survival analysis and Log-Rank p-Value	Predicted Target	p-value	Correlation coefficient r
ZNF518A	0.001	CCNJ	<0.001	0.644
PARGP1	0.003	AGAP4	<0.001	0.571
CDK6-AS1	0.005	CDK6	<0.001	0.459
AC087672.2	0.023	LY96	<0.001	0.446
AC025539.1	0.032	HS3ST1	<0.001	0.461
LINC02731	0.033	JAM3	<0.001	0.564
AL353764.1	0.035	NFIL3	<0.001	0.533
AC025183.1	0.040	IRX4	<0.001	0.769
AC009950.1	0.043	SP140	<0.001	0.453
AC009950.1	0.043	SP140L	<0.001	0.514
FZD10-AS1	0.046	FZD10	<0.001	0.788
SLC7A11-AS1	0.046	SLC7A11	<0.001	0.638

ERNA PARGP1 is the key eRNA in PCa

PARGP1 was found to be a critical eRNA in PCa. In PCa patients, overall survival was better in the PARGP1 low-expression group than in the high-expression group (Kaplan–Meier log-rank test, p = 0.003; Figure 1(a)). Levels of PARGP1 and the target gene AGAP4 mRNA were positively correlated (Figure 1(b), r = 0.571, p < 0.001). AGAP4 protein expression was significantly different between PCa patients and normal subjects (Figure 2).

Figure 1. Impact of eRNA PARGP1 on prostate adenocarcinoma (PRAD). (a) Kaplan–Meier survival curve for patients with PARGP1-high and PARGP1-low expression. (b) Scatterplot showing the association between PARGP1 and AGAP4 levels.

Figure 2. Expression of AGAP4 in tumor and normal prostate tissue.

The effect of PARGP1 on overall survival in the other 32 tumors was limited to prostate adenocarcinoma (PRAD), head and neck squamous cell carcinoma (HNSC), brain lower grade glioma (LGG), and pheochromocytoma and paraganglioma (PCPG) with high specificity (Table 2, Figure 3).

Figure 3. Impact of eRNA PARGP1 on other cancers. (a) Kaplan–Meier overall survival curve for patients with PARGP1-high and PARGP1-low expression on head and neck squamous cell carcinoma (HNSC). (b) Kaplan–Meier overall survival curve for patients with PARGP1-high and PARGP1-low expression in brain lower grade glioma (LGG). (c) Kaplan–Meier overall survival curve for patients with PARGP1-high and PARGP1-low expression on pheochromocytoma and paraganglioma (PCPG).

Table 2. Kaplan–Meier survival analysis and gene expression correlations for PARGP1 and AGAP4 across 33 cancer types from the Cancer Genome Atlas (TCGA) project.

Tumor Type			PARGP1 and AGAP4
Abbrev-iation	Detail	PARGP1Overall Survival Log-Rank P-Value	Correlation P-Value	Correlation Coefficient r
ACC	Adrenocortical carcinoma	0.073	<0.001	0.597
BLCA	Bladder Urothelial Carcinoma	0.308	<0.001	0.559
BRCA	Breast invasive carcinoma	0.715	<0.001	0.636
CESC	Cervical squamous cell carcinoma and endocervical adenocarcinoma	0.296	<0.001	0.590
CHOL	Cholangiocarcinoma	0.408	0.004	0.473
COAD	Colon adenocarcinoma	0.904	<0.001	0.688
DLBC	Lymphoid Neoplasm Diffuse Large B-cell Lymphoma	0.667	<0.001	0.698
ESCA	Esophageal carcinoma	0.262	<0.001	0.790
GBM	Glioblastoma multiforme	0.394	<0.001	0.665
HNSC	Head and Neck squamous cell carcinoma	≤0.001	<0.001	0.671
KICH	Kidney Chromophobe	0.191	0.001	0.390
KIRC	Kidney renal clear cell carcinoma	0.364	<0.001	0.766
KIRP	Kidney renal papillary cell carcinoma	0.746	<0.001	0.697
LAML	Acute Myeloid Leukemia	0.277	<0.001	0.634
LGG	Brain Lower Grade Glioma	0.003	<0.001	0.720
LIHC	Liver hepatocellular carcinoma	0.072	<0.001	0.677
LUAD	Lung adenocarcinoma	0.203	<0.001	0.731
LUSC	Lung squamous cell carcinoma	0.684	<0.001	0.749
MESO	Mesothelioma	0.357	<0.001	0.519
OV	Ovarian serous cystadenocarcinoma	0.171	<0.001	0.671
PAAD	Pancreatic adenocarcinoma	0.060	<0.001	0.666
PCPG	Pheochromocytoma and Paraganglioma	0.033	<0.001	0.437
PRAD	Prostate adenocarcinoma	0.003	<0.001	0.571
READ	Rectum adenocarcinoma	0.153	<0.001	0.645
SARC	Sarcoma	0.151	<0.001	0.677
SKCM	Skin Cutaneous Melanoma	0.586	<0.001	0.760
STAD	Stomach adenocarcinoma	0.888	<0.001	0.760
TGCT	Testicular Germ Cell Tumors	0.609	<0.001	0.467
THCA	Thyroid carcinoma	0.231	<0.001	0.544
THYM	Thymoma	0.699	<0.001	0.592
UCEC	Uterine Corpus Endometrial Carcinoma	0.202	<0.001	0.447
UCS	Uterine Carcinosarcoma	0.483	<0.001	0.605
UVM	Uveal Melanoma	0.746	<0.001	0.891

PARGP1 is related to the state of biochemical recurrence of PCa

There was no significant difference in PARGP1 expression among PCa patients of different ages, pathological T-stages, N-stages, or distant metastases (Figure 4(a–d), p > 0.05). However, there was a significant difference in PARGP1 expression between patients with PSA less than 4 ng/ml and those with 4–10 ng/ml (Figure 4(e), p < 0.05). The expression of PARGP1 was significantly higher in patients with biochemical recurrence than in those without (Figure 4(g), p < 0.001), and the expression level of the target gene AGAP4 showed consistent results (Figure 4(h), p < 0.001).

Figure 4. Gene expression related to prostate cancer. (a–d) Expression levels of PARGP1 were not related to age in patients with prostate cancer or to the TNM stage of the tumor. (e–f) The difference between the expression level of PARGP1 and the value of PSA. There was a significant difference in PARGP1 expression level for PSA between 0–4 ng/ml and 4–10 ng/ml. (g–h) The expression levels of PARGP1 and AGAP4 were related to biochemical recurrence of prostate cancer (p < 0.05).

Function and pathway enrichment analysis

The GO and KEGG pathways were analyzed using the ‘clusterProfiler’ package in R software to explore further the molecular mechanisms of the 892 genes co-expressed with PARGP1 on the occurrence and progression of PCa. The functions of these genes were focused on biological processes, molecular functions, and cellular composition (Figure 5(a)). The main pathway for KEGG enrichment was infection with herpes simplex virus 1 (HSV-1) (Figure 5(b)).

Figure 5. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEEG)enrichment analysis of co-expressed genes. The ordinate is the name of the metabolic pathway in GO and KEGG, and the abscissa is the number of genes annotated to the pathway.

Discussion

There is a growing body of evidence suggesting that enhancers have inherent promoter capabilities and can transcribe non-coding RNA (Kim et al. 2010). Enhancer transcription and its product, enhanced eRNA, have been widely studied in mouse and human cells, which possess several common characteristics. Mammalian eRNA is usually low in abundance, unspliced, and retained in the nucleus, and they are usually transcribed in both directions (Andersson et al. 2014; Mikhaylichenko et al. 2018). In human cancers, the importance of eRNAs in the activation of oncogenes or oncogenic signaling pathways has received increasing attention. For example, the activation of ESR1 can increase eRNA transcription in breast cancer (Li et al. 2016). Together, these findings suggest the critical role of eRNA in tumorigenesis and their potential effects in eRNA-targeted therapy (Léveillé et al. 2015).

Statistics show that the incidence of PCa is increasing in younger ages, with the incidence of PCa in people aged 15–40 years steadily increasing at a rate of 2% per year; men under 40 years of age are more likely to develop metastatic disease and have a higher mortality rate than older men (Bleyer et al. 2020; Sartor 2020). PSA has dramatically improved the detection of PCa as a traditional screening indicator; however, when PSA values are between 4–10 ng/ml, we need to do more specific tests to rule out other factors, including inflammation or catheterization.

In the present study, there were significant differences in expression levels of PARGP1 in patients with PSA values 0–4 ng/ml and 4–10 ng/ml, suggesting that PARGP1 can be used as a potential diagnostic marker to differentiate benign from malignant prostate lesions when the PSA is less than 10 ng/ml. Expression levels of PARGP1 in patients with biochemical recurrence were significantly higher than those without, and the target gene AGAP4 showed the same results. Eexpression of PARGP1 and AGAP4 may be associated with a significantly increased risk of biochemical recurrence. Based on these findings, we screened the functionally unannotated eRNA PARGP1 and found that it may be one of the essential eRNAs in PCa. Based on GO enrichment analysis, the predominant KEGG pathway was HSV-1 infection. HSV-1 is a neurotropic herpesvirus; the γ34.5 gene is a critical determinant of the neurovirulence of this virus (Taneja et al. 2001). Recent studies showed that UV-inactivated HSV-1 is involved in killing PCa cells by activating natural killer cells via Toll-like receptor 2. HSV-1 may in the future be a novel option for the treatment of PCa, especially neuroaggressive PCa (Kim et al. 2012; Vogel et al. 2014; Samudio et al. 2019).

Our study also has some limitations. Our analysis based on TCGA yielded PARGP1, the most critical eRNA in PCa, and AGAP4, the predicted target gene, which only provides a direction for subsequent PCa research and has not been validated in large samples. Nevertheless, based on the analysis of these large data sets, we found that PARGP1 was also associated with patient survival in PRAD, HNSC, LGG, and PCPG with high specificity, limiting the scope of the study and reducing unnecessary duplication of studies. PARGP1 may also be used as a potential diagnostic marker for the diagnosis of PCa and the detection of biochemical recurrence. Novel therapeutic approaches targeting PARGP1 expression may help induce protective immunity to treat these cancers effectively.

Conclusions

We analyzed the effect of PARGP1 on PCa outcome and found that it can be used as an ancillary means to distinguish between PSA 0–4 ng/ml and 4–10 ng/ml in patients with PCa. It is also related to the biochemical recurrence in patients with PCa. Detection of the PARGP1 expression level is critical for predicting patients’ outcome, diagnosis, and management with PCa.

Authors’ contributions

Project design and implementation were conceived by Weiguo Chen, and Feng Zhou, and Xi Zhang. Data collection, collation and analysis was performed by Qi Zhou, and Zhiyu Zhang, and Lu Ma. Manuscript drafting and editing were performed by Xiaojie Ang, and Zekun Xu. All authors read and approved the final manuscript.

Data sharing statement

Available in ‘figshare’ at https://doi.org/10.6084/m9.figshare.13498374.v3.

Ethical statement

The authors are accountable for all aspects of the work in ensuring that questions related to the accuracy or integrity of any part of the work are appropriately investigated and resolved.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Additional information

Funding

This study was supported by the Project of State Key Laboratory of Radiation Medicine and Protection, Soochow University [grant number GZK1202017], Gusu Medical Talent Foundation [grant number GSWS2020021].