Signature Panel of 11 Methylated mRNAs and 3 Methylated lncRNAs for Prediction of Recurrence-Free Survival in Prostate Cancer Patients

Abstract

Background

Radical prostatectomy is the main treatment for prostate cancer (PCa), a common cancer type among men. Recurrence frequently occurs in a proportion of patients. Therefore, there is a great need to early screen those patients to specifically schedule adjuvant therapy to improve the recurrence-free survival (RFS) rate. This study aims to develop a biomarker to predict RFS for patients with PCa based on the data of methylation, an important heritable contributor to carcinogenesis.

Methods

Methylation expression data of PCa patients were downloaded from The Cancer Genome Atlas (TCGA), Gene Expression Omnibus database (GSE26126), and the European Bioinformatics Institute (E-MTAB-6131). The stable co-methylation modules were identified by weighted gene co-expression network analysis. The genes in modules were overlapped with differentially methylated RNAs (DMRs) screened by MetaDE package in three datasets, which were used to screen the prognostic genes using least absolute shrinkage and selection operator analyses. The prognostic performance of the prognostic signature was assessed by survival curve analysis.

Results

Five co-methylation modules were considered preserved in three datasets. A total of 192 genes in these 5 modules were overlapped with 985 DMRs, from which a signature panel of 11 methylated messenger RNAs and 3 methylated long non-coding RNAs was identified. This signature panel could independently predict the 5-year RFS of PCa patients, with an area under the receiver operating characteristic curve (AUC) of 0.969 for the training TCGA dataset and 0.811 for the testing E-MTAB-6131 dataset, both of which were higher than the predictive accuracy of Gleason score (AUC = 0.689). Also, the patients with the same Gleason score (6–7 or 8–10) could be further divided into the high-risk group and the low-risk group.

Conclusion

These results suggest that our prognostic model may be a promising biomarker for clinical prediction of RFS in PCa patients.

Keywords:

• prostate cancer
• recurrence-free survival
• methylation
• prognostic signature

Introduction

Prostate cancer (PCa) is one of the most common cancers among men. It was estimated that there were 248,530 new cases and 34,130 deaths in the United States in 2021.¹ Radical prostatectomy is considered as the first-line treatment option for PCa patients, but recurrence (biochemical, BCR; or clinical) frequently occurs in 10–40% of patients after curative surgery,^2–4 which leads to the cumulative 5-year recurrence-free survival (RFS) rate of only about 65%.^3,5 Clinically, prostate-specific antigen (PSA) level,⁶ Gleason score,⁷ and tumor, node, metastasis (TNM) staging⁸ are widely used to predict tumor recurrence. Nevertheless, their prediction accuracy remains unsatisfactory (<70%).⁹ Therefore, there is a great need to identify more effective biomarkers for early screening patients who will possess poor RFS and then specifically scheduling postoperative radiotherapy and chemotherapy for them to reduce the probability of recurrence.¹⁰

Recent accumulating evidence pinpoints that epigenetic modifications are heritable contributors to the carcinogenesis of PCa;^11,12 among them, DNA methylation is one of the most important epigenetic modifications of the genome,¹³ indicating methylated genes may represent potential biomarkers for the prediction of RFS in PCa patients. This hypothesis has been confirmed in some studies as following. Protocadherin (PCDH)-10 and PCDH17 methylation were detected in patients with PCa, but not in controls.^14,15 The high methylation levels of PCDH10 and PCDH17 were significantly associated with poor BCR-free survival.^14,15 A meta-analysis of seven studies showed that the 5-year BCR-free survival for patients with a high methylation status in the paired-like homeodomain transcription factor 2 (PITX2) gene was significantly lower than that for patients with a low methylation status (71% vs 90%).¹⁶ Xu et al reported that patients with a lower methylation level of the long interspersed nucleotide elements (LINE-1) gene [hazard ratios (HR) = 3.34, 95% confidence interval (CI) = 1.32–8.45] and a higher methylation level of DNA repeats D4Z4 gene [HR = 4.12, 95% CI = 1.32–12.86] exhibited markedly increased risks of BCR and significantly shorter BCR-free survival time.¹⁷ Stott-Miller et al found the presence of promoter hypermethylation of the homeobox D3 (HOXD3) gene in patients with PCa recurrence compared with those without recurrence. The median time for RFS was shorter in the high methylation group than that in the HOXD3 low methylation group.¹⁸ Analysis of The Cancer Genome Atlas (TCGA) data demonstrated that the methylation status in the gene of programmed death-1 receptor (PDCD1)¹⁹ or solute carrier organic anion transporter family member 4C1 (SLCO4C1)²⁰ was an independent prognostic biomarker for poor RFS in patients with PCa. Combination with a 52-gene methylation signature was also indicated to improve the prediction power of standard clinical-pathological parameters for RFS [the area under the receiver operating characteristic (ROC) curve (AUC): 0.78 vs 0.73].²¹ However, effective epigenetic biomarkers for RFS prediction remain lacking in PCa.

Our purpose in this study was to develop and validate a novel methylation signature panel to identify patients at a high risk of poor RFS using the TCGA dataset and data collected from Gene Expression Omnibus (GEO) or European Bioinformatics Institute (EMBL-EBI) databases. The superior prognostic performance of this signature panel to clinical parameters was evaluated comprehensively, including AUC, concordance index (C-index), and stratification analyses. Compared with the study of Geybels et al,²¹ our methylation signature not only included protein-encoding messenger RNAs (mRNAs), but also long non-coding RNAs (lncRNAs). Previous studies suggested that the prognostic power of methylated lncRNAs-²² or lncRNA-mRNA-based²³ prognostic signature was higher than that of mRNA alone. Accordingly, our signature may be more helpful for predicting prognosis and guiding the individualized treatment of PCa patients.

Materials and Methods

Data Collection and Preprocessing

Matched methylation (platform, Illumina Human Methylation 450), level 3 mRNA-seq (platform, Illumina HiSeq 2000 RNA Sequencing), and corresponding clinical data of PCa patients were obtained from the TCGA database (https://portal.gdc.cancer.gov/) on October 25, 2019, in which 51 samples were normal controls and 495 were PCa (375 samples reported the recurrence status, including 47 from patients with recurrence and 328 from patients without recurrence). Furthermore, methylation datasets of PCa were also downloaded from GEO (http://www.ncbi.nlm.nih.gov/geo/) or EMBL-EBI database (https://www.ebi.ac.uk). A total of 83 samples reporting the status of recurrence (recurrence, n = 17; non-recurrence, n = 68) were included in the GSE26126 dataset (platform, Illumina HumanMethylation27 BeadChip),²⁴ which was used for the validation of modules; there were 108 samples (recurrence, n = 15; non-recurrence, n = 93) in the E-MTAB-6131 dataset (platform, Illumina Human Methylation 450), which was used for validations of both module and survival results.

The annotations of lncRNAs and mRNAs in each dataset were performed using the HUGO Gene Nomenclature Committee (http://www.genenames.org/).²⁵ RNAs with a median expression value of 0 were deleted. The overlapped lncRNAs and mRNAs in all three datasets were used for the following analyses.

Identification of Co-Methylation Modules

Based on all shared methylated RNAs in TCGA (training), GSE26126 (testing), and E-MTAB-6131 (testing) datasets, co-methylation networks were constructed using the weighted gene co-expression network analysis (WGCNA) software (v1.61; https://cran.r-project.org/web/packages/WGCNA/index.html).²⁶ Briefly, the verboseScatterplot function was conducted to explore the correlations in the methylation level and the connectivity of all RNAs between any two datasets to confirm their comparability. Based on the criterion of scale-free topology, the pickSoftThreshold function was used to select an appropriate soft-thresholding power (β) to construct the weighted adjacency matrix which was subsequently transformed into a topological overlap matrix (TOM), a measure for the correlations between the methylation levels of two genes. The hierarchical clustering dendrogram was obtained based on the TOM-based dissimilarity. The DynamicTreeCut function²⁷ was applied to identify modules with cutHeight = 0.995 and minSize = 100. The modulePreservation function²⁸ was carried out to assess the preservation of identified modules between the training set and two testing sets, with Z-score >5 as the statistical threshold. In addition, moduleTraitCor and moduleTraitPvalue algorithms were chosen to calculate the correlation between module eigengenes and clinical traits in the TCGA dataset.

Identification of Differentially Methylated RNAs (DMRs) and Expressed RNAs (DERs)

The DMRs, including differentially methylated lncRNAs (DMLs) and differentially methylated genes (DMGs) between recurrence and non-recurrence PCa samples, were identified using the MetaDE.ES function in the MetaDE package (v1.0.5, https://cran.r-project.org/web/packages/MetaDE/). Briefly, the heterogeneity was assessed across three datasets (TCGA, GSE26126, and E-MTAB-6131) using tau² statistic and Chi-square-based Q-test. Genes with tau² = 0 and Qpval >0.05 were considered to be homogeneous and used for the differential analysis. The gene expression difference was determined by the MetaDE.pvalue algorithm, with a false discovery rate (FDR) <0.05 selected as the significance threshold. The expression consistency of DMLs and DMGs in three datasets was detected by calculating the log₂FC(fold change). The heatmap.sig.genes function in the MetaDE package was used to plot the heatmap of DMRs in three datasets.

The differentially expressed RNAs (DERs) between recurrence and non-recurrence PCa samples or between PCa and normal controls were identified using limma package in R (v3.34.7; https://bioconductor.org/packages/release/bioc/html/limma.html).²⁹ A p-value of <0.05 was considered as the statistical threshold.

Development and Validation of a Prognostic Scoring Model Based on DMRs

A Venn diagram (http://bioinformatics.psb.ugent.be/webtools/Venn/) was constructed to screen the overlap between co-methylation module RNAs and DMRs. The shared genes were used for survival analysis. Based on the survival information of patients in the TCGA dataset, univariate Cox regression analysis in the “survival” package of R (v2.41–1; http://bioconductor.org/packages/survivalr/) was used to evaluate the association between the methylation levels of DMRs and RFS. Significant DMRs with a log-rank p < 0.05 were selected for multivariate Cox regression analysis. The signature identified by multivariate analysis was further set as input for an L1 penalized (Least Absolute Selection and Shrinkage Operator, LASSO) Cox-proportional hazard model analysis (penalized package, v0.9–5; http://bioconductor.org/packages/penalized/)^30,31 to obtain an optimal signature panel for prognosis prediction. The prognostic risk scoring model was constructed based on the methylation levels of prognostic RNAs (Methy_DMRs) and their LASSO coefficients (β_DMRs):

Risk score = β_DML1 × Exp_DML1 … + β_DMLn × Exp_DMLn + β_DMG1 × Exp_DMG1 + ….β_DMGn × Exp_DMGn

The risk score was calculated for each patient in the TCGA dataset. Patients were divided into a low-risk group (risk score below the median value) and a high-risk group (risk score above the median value). Kaplan–Meier survival curve analysis and the Log rank test were used to estimate the RFS time of two risk groups. ROC analysis with calculation of AUC was used to evaluate the prognostic ability of the risk scoring model. The prognostic robustness of the risk scoring system was validated in the testing set (E-MTAB-6131).

Furthermore, univariate and multivariate Cox regression and stratification analyses were conducted to estimate the independent prognostic ability of the risk scoring model from clinical characteristics, with p < 0.05 as the statistical significance. The superior predictive efficiency of the risk score for RFS compared with other clinical characteristics was determined by time-dependent ROC curve and C-index analyses.

Function Enrichment Analyses of the Prognostic RNAs

To explore the function of prognostic mRNAs, public web servers gProfiler (http://biit.cs.ut.ee/gprofiler/gost)³² and The Database for Annotation, Visualization and Integrated Discovery (DAVID; v6.8; https://david.ncifcrf.gov) were searched. Data sources included Gene Ontology (GO) terms, Kyoto Encyclopedia of Genes and Genomes (KEGG), Reactome, and Wiki pathways. A p-value <0.05 was selected as the statistical threshold.

The functions of lncRNAs or some mRNAs were predicted according to their co-expressed mRNAs. The cor.test function (https://stat.ethz.ch/R-manual/R-devel/library/stats/html/cor.test.html) in R was used to calculate the Pearson correlation coefficients (PCC) between prognostic lncRNAs and stable module mRNAs or between module mRNAs using expression and methylation data, respectively. The co-expression pairs with the PCC > 0.4 were used to construct the co-expression network which was visualized in the Cytoscape software (v3.6.1; www.cytoscape.org/).

Results

Identification of Co-Methylation Modules

A total of 9745 RNAs were found to be in interaction with three datasets (TCGA, GSE26126, and E-MTAB-6131). Thus, methylation levels of these RNAs were collected for the WGCNA analysis to mine PCa-associated co-methylation modules. The methylation levels of these RNAs were positively correlated between any two datasets (TCGA-GSE26126, cor = 0.84, p < 1e-200; TCGA-E-MTAB-6131, cor = 0.86, p < 1e-200; GSE26126-E-MTAB-6131, cor = 0.9, p < 1e-200), indicating a good comparability between these datasets. The soft-threshold power of β was selected as 8 to construct a scale-free network (scale-free R² = 0.9, Figure 1A; mean connectivity = 1, Figure 1B). Genes in the TCGA dataset were clustered into 10 co-methylation modules represented by branches with different colors in the dendrogram (black, blue, brown, green, grey, magenta, pink, red, turquoise, and yellow) (Figure 1C; Table 1). These co-methylation modules were also formed in the analysis of GSE26126 (Figure 1D) and E-MTAB-6131 (Figure 1E) datasets using the same manner as the TCGA dataset. Among these 10 modules, the black, blue, brown, yellow, and turquoise modules were considered to be preserved (Z-score >5; Table 1). These five preserved modules were also proved to be significantly associated with crucial clinical characteristics of PCa patients, such as recurrence (turquoise, blue: negatively associated; black, brown, yellow: positively associated) (Figure 1F). These findings suggest that the genes (including 39 lncRNAs and 2531 mRNAs) in these five modules may be PCa recurrence-associated.

Table 1 Co-Methylation Modules Identified by WGCNA

ID	Color	Module Size	Number of mRNAs	Number of lncRNAs	Preservation Z-score
Module 1	Black	206	206	0	6.6439
Module 2	Blue	374	362	12	8.5168
Module 3	Brown	302	295	7	7.8463
Module 4	Green	234	234	0	0.6804
Module 5	Grey	1902	1868	34	2.9263
Module 6	Magenta	124	124	0	2.4207
Module 7	Pink	203	203	0	4.6150
Module 8	Red	233	224	9	4.2880
Module 9	Turquoise	1438	1427	11	20.7257
Module 10	Yellow	250	241	9	7.4672

Note: Bold indicated the preserved modules identified in three datasets.

Abbreviations: WGCNA, weighted gene co-expression network analysis; mRNAs, messenger RNAs; lncRNAs, long non-coding RNAs.

Figure 1 WGCNA analysis. (A) soft-threshold power β selected when the R² reached 0.9 for the first time; (B) the mean connectivity corresponding to different power β values; C-E, clustering dendrogram of co-methylation modules from TCGA (C), GSE26126 (D) and EMBL-EBI (E); (F) the correlations between modules and clinical traits.

Identification of DMRs in PCa Recurrence Samples

Using the metaDE method, 985 DMRs (including 888 DMGs and 7 DMLs) were identified in PCa recurrence samples compared with non-recurrence tissues, including 577 hypomethylated and 318 hypermethylated (Figure 2A). A total of 192 DMRs (including 5 lncRNAs and 187 mRNAs) were shared with the genes in the five preserved modules (Figure 2B), consisting of 16 (all were mRNAs) in the black module, 30 (including two lncRNAs and 28 mRNAs) in the blue module, 25 (all were mRNAs) in the brown module, 101 (including three lncRNAs and 98 mRNAs) in the turquoise module and 20 (all were mRNAs) in the yellow module (Figure 2C). These findings indicate that these 192 DMRs may represent epigenetic biomarkers for PCa recurrence and were used for the following survival analysis.

Figure 2 Identification of differentially methylated genes. (A) Heat map of differentially methylated RNAs in three datasets; (B) Venn diagram to obtain the overlap between differentially methylated RNAs and module genes; (C) the module distribution of the overlapped genes.

Identification of an Epigenetic Signature Panel and Development of a Risk Scoring Model for RFS Prediction in PCa Patients

Univariate Cox regression analysis identified the expressions of 42 DMRs (including four DMLs and 38 DMGs) were significantly associated with RFS of PCa patients (log-rank p < 0.05). Multivariable Cox regression model demonstrated that fifteen of them (including four DMLs and 11 DMGs) were independent prognostic factors. LASSO modelling further filtered 14 DMRs as the optimal prognostic panel [DMLs: MEG3 (maternally expressed 3), DSCR9 (Down syndrome critical region 9), HCP5; DMGs: MMP7 (matrix metallopeptidase 7), SLCO3A1 (solute carrier organic anion transporter family member 3A1), KCNF1 (potassium voltage-gated channel modifier subfamily F member 1), RFXAP (regulatory factor X associated protein), NTRK3 (neurotrophic receptor tyrosine kinase 3), NAV1 (neuron navigator 1), HOXA13 (homeobox A13), HAS2 (hyaluronan synthase 2), CBX2 (chromobox 2), HIST1H2AJ (H2A clustered histone 14), SNX4 (sorting nexin 4)] (Table 2). Six prognostic genes (MEG3, MMP7, SLCO3A1, KCNF1, RFXAP, and NTRK3) had positive LASSO coefficient and HR > 1, suggesting patients with the high methylation levels of these genes may have poor RFS, while the other eight genes (DSCR9, HCP5, NAV1, HOXA13, HAS2, CBX2, HIST1H2AJ, and SNX4) had negative LASSO coefficient and HR < 1, implying the high methylation levels may play a protective role against recurrence and caused death (Table 2).

Table 2 The Optimal Methylation Signature Panel for Prognosis Prediction

Symbol	Type	Module	Methylation Level	Univariate Cox Regression Analysis			LASSO Coefficient
				HR	95% CI	P-value
MEG3	lncRNA	Blue	Downregulated	1.180	1.034–4.021	1.19E-02	0.2901
MMP7	mRNA	Turquoise	Upregulated	1.631	1.045–7.561	2.80E-03	0.4701
SLCO3A1	mRNA	Turquoise	Upregulated	5.852	1.018–7.363	4.80E-02	0.4693
KCNF1	mRNA	Turquoise	Upregulated	2.260	1.203–3.129	8.70E-03	0.3800
RFXAP	mRNA	Black	Upregulated	1.693	1.034–4.108	7.10E-03	0.3384
NTRK3	mRNA	Turquoise	Upregulated	4.930	1.260–6.929	2.20E-02	0.2656
DSCR9	lncRNA	Turquoise	Upregulated	0.541	0.169–0.815	4.41E-02	−0.0112
HCP5	lncRNA	Turquoise	Downregulated	0.367	0.095–0.423	1.05E-02	−0.1128
NAV1	mRNA	Turquoise	Downregulated	0.117	0.0184–0.739	2.30E-02	−0.1493
HOXA13	mRNA	Turquoise	Downregulated	0.224	0.0509–0.987	4.80E-02	−0.2229
HAS2	mRNA	Turquoise	Downregulated	0.170	0.0308–0.938	4.20E-02	−0.5077
CBX2	mRNA	Turquoise	Downregulated	0.011	0.000507–0.216	3.10E-03	−0.6602
HIST1H2AJ	mRNA	Turquoise	Downregulated	0.002	0.000844–0.503	2.70E-02	−1.0755
SNX4	mRNA	Yellow	Downregulated	0.755	0.138–0.941	2.30E-02	−2.3318

Abbreviations: HR, hazard ratio; CI, confidence interval; LASSO, least absolute shrinkage, and selection operator.

A risk scoring model was established based on the methylation levels of the above 14 genes and their corresponding LASSO coefficients (Table 2). Using the median risk score in each dataset as the cut-off, patients were assigned to the low-risk group and the high-risk group. Obviously, more patients in the high-risk group developed recurrence (44/187 vs 3/188, p < 0.001). Kaplan–Meier curves showed that the high-risk group had a significantly shorter RFS rate compared with the low-risk group (TCGA: HR = 18.72, 95% CI = 5.800–60.43, p = 3.521e-13, Figure 3A; E-MTAB-6131: HR = 3.455, 95% CI = 1.089–10.90, p = 2.417e-02, Figure 3B). ROC curve analysis showed that AUC of the training cohort was 0.987, 0.947, and 0.969 in predicting 1-, 3- and 5-year RFS, respectively (Figure 3C); while it was 0.851, 0.789, and 0.811 for the testing dataset in predicting 1-, 3- and 5-year RFS, respectively (Figure 3D). These findings indicate that our epigenetic signature panel had high accuracy in predicting RFS for patients with PCa.

Figure 3 The prediction performance of our 14-DMRs-based risk score system for recurrence-free survival. (A) Kaplan-Meier survival curve analysis of the TCGA dataset; (B) Kaplan-Meier survival curve analysis of the EMBL-EBI dataset; (C) receiver operator characteristic (ROC) curve analysis of the TCGA dataset; (D) receiver operator characteristic curve analysis of the EMBL-EBI dataset.

Abbreviations: HR, hazard ratio; AUC, area under the ROC curve.

To explore whether our epigenetic signature panel was independent of clinical pathological characteristics for RFS prediction, univariate and multivariate Cox regression analyses were performed. As a result, Pathologic_T, Radiation therapy, Targeted molecular therapy, Gleason score (Figure 4A), PSA, and the risk score status exhibited significant associations with RFS in univariate Cox regression analysis, but only Gleason score and the risk score status remained significant in multivariate Cox regression analysis (Table 3), revealing these two variables were independent prognostic factors. To investigate whether the risk score was superior to Gleason score for RFS prediction, stratification, time-dependent ROC curve, and C-index analyses were then conducted. Stratification analysis showed that patients with the same Gleason score [(6–7) (Figure 4B) or (8–10) (Figure 4C)] could be further divided into the high-risk group and the low-risk group according to their risk score. Similarly, AUC (0.959 vs 0.689) and C-index (0.885 vs 0.723) of Gleason score were lower than those of the risk score (Figure 4D). These findings implied an improved prognostic power of our epigenetic signature panel. Thereby, our signature panel can be integrated with the Gleason score model to obtain better prognostic effects in clinic. This theory is verified in Figure 4D, in which AUC and C-index for the combined model were higher than those of Gleason score (AUC = 0.984 vs 0.689; C-index = 0.898 vs 0.723) and risk score (AUC = 0.984 vs 0.959; C-index = 0.898 vs 0.885) alone, respectively.

Table 3 Univariate and Multivariate Cox Regression of Clinical Features and Risk Score

Variables	TCGA (N=375)	Univariate Analysis			Multivariate Analysis
		HR	95% CI	P-value	HR	95% CI	P-value
Age (years, mean ± SD)	60.85 ± 6.84	1.023	0.979–1.068	3.02E-01	-	-	-
Pathologic_M (M0/M1/-)	354/2/19	6.115	0.893–16.35	8.49E-01	-	-	-
Pathologic_N (N0/N1/-)	273/54/48	1.691	0.853–3.353	1.51E-01	-	-	-
Pathologic_T (T2/T3/T4/-)	144/219/6/6	2.865	1.586–5.177	4.22E-04	0.976	0.457–2.085	9.51E-01
Radiation therapy (Yes/No/-)	50/317/8	2.189	1.080–4.438	4.40E-02	1.006	0.424–2.388	9.88E-01
Targeted molecular therapy (Yes/No/-)	324/42/9	2.522	1.276–4.986	5.85E-03	0.761	0.325–1.776	5.26E-01
Gleason score (6/7/8/9/10)	37/181/55/98/4	2.57	1.858–3.554	4.89E-10	2.196	1.451–3.322	1.99E-04
Prostate-specific antigen	1.51 ± 3.17	1.057	1.020–1.095	4.43E-04	1.016	0.978–1.055	4.16E-01
Risk score status (High/Low)	187/188	18.72	5.800–60.43	3.52E-13	12.67	3.845–41.82	3.02E-05

Note: Bold indicated the statistical results met the significance threshold of p-value < 0.05.

Abbreviations: SD, standard deviation; M, metastasis; N, node; T, tumor; HR, hazard ratio; CI, confidence interval; TCGA, The Cancer Genome Atlas.

Figure 4 The superiority of our 14-DMRs-based risk score system to clinical indicators. (A) Kaplan-Meier survival curve analysis to show the association of Gleason score with recurrence-free survival; (B) stratification analysis for Gleason score (6–7) using the risk score; (C) stratification analysis for Gleason score (8–10) using the risk score; (D) time-dependent ROC curve analysis constructed according to various models.

Abbreviations: HR, hazard ratio; AUC, area under the receiver operator characteristic curve; C-index, concordance index.

Function Analysis of Prognostic Genes

To obtain the possible functions of prognostic genes, 888 DMGs were used as the input to search the gProfiler and DAVID databases. The results showed the prognostic genes were mainly involved in GO:0042127~regulation of cell population proliferation (NTRK3), GO:0006915~apoptotic process (NTRK3), GO:0022407~ regulation of cell–cell adhesion (HAS2) (Table S1), GO:0045893~positive regulation of transcription, DNA-templated (RFXAP), GO:0030335~positive regulation of cell migration (NTRK3, HAS2), GO:0043065~positive regulation of apoptotic process (NTRK3), GO:0000122~negative regulation of transcription from RNA polymerase II promoter (CBX2), and hsa04310:Wnt signaling pathway (MMP7) (Table S2). All these biological processes were carcinogenesis-related, further confirming their importance for PCa.

A total of 148 co-expression pairs were obtained between prognostic lncRNAs and module mRNAs according to their mRNA expression (Figure 5), from which MEG3 was shown to co-express with RFXAP (cor = 0.78); DSCR9 could co-express with CBX2 (cor = 0.61). Thus, the functions of lncRNAs MEG3 and DSCR9 may be associated with the functions of RFXAP and CBX2, respectively. Furthermore, the co-expression methylation relationships were calculated for genes in each module, from which HCP5 were found to interact with HIST1H2AJ (cor = 0.31); while HIST1H2AJ may interact with CBX2 (cor = 0.03). Thus, the functions of HCP5 and HIST1H2AJ may also be similar to CBX2.

Figure 5 The co-expression network between differentially methylated lncRNAs and mRNAs identified unstable modules. The color indicated the corresponding module.

Validation of the RNA Expression Levels of Prognostic Genes

A total of 115 DERs were identified between recurrence and non-recurrence PCa samples. Among them, prognostic MEG3, NTRK3, NAV1, CBX2, and SNX4 were found to be DERs and their RNA expression levels were opposite to their methylation levels (Table 4). Furthermore, 169 DERs were obtained between PCa and normal controls, among which SLCO3A1, KCNF1, RFXAP, NTRK3 (downregulated), CBX2, HIST1H2AJ, and SNX4 (upregulated) were opposite to their methylation levels (Table 4). These findings revealed that mRNA expressions of these genes may be driven by their methylation levels.

Table 4 The RNA Expression Levels of Prognostic Genes

Symbol	Type	Recurrence vs Non-Recurrence			Tumor vs Control
		Log2FC	P-value	FDR	Log2FC	P-value	FDR
MEG3	lncRNA	0.166016	0.0313	0.2537	−0.35998	1.04E-03	0.008423
MMP7	mRNA	−0.00127	0.122	0.990023	0.007234	1.15E-01	0.931766
SLCO3A1	mRNA	0.07484	0.0618	0.500229	−0.84417	1.74E-16	1.41E-15
KCNF1	mRNA	−0.09793	0.0817	0.661785	−1.66242	2.98E-06	2.41E-05
RFXAP	mRNA	0.032439	0.0552	0.447327	−0.05774	6.40E-03	0.051835
NTRK3	mRNA	−0.36891	0.00523	0.042399	−1.04021	1.15E-07	9.35E-07
DSCR9	lncRNA	0.372407	0.0162	0.131543	0.832172	3.69E-11	2.99E-10
HCP5	lncRNA	0.019159	0.0984	0.796839	−0.07392	4.06E-02	0.32906
NAV1	mRNA	0.247454	0.0029	0.02352	−0.13033	1.24E-02	0.100329
HOXA13	mRNA	0.007909	0.108	0.873183	−0.19936	3.49E-06	2.83E-05
HAS2	mRNA	0.140339	0.068	0.551047	−0.87177	1.51E-04	0.001222
CBX2	mRNA	0.613398	2.81E-05	0.000228	1.083018	1.36E-16	1.10E-15
HIST1H2AJ	mRNA	−0.06935	0.103	0.830599	1.24546	2.18E-03	0.017627
SNX4	mRNA	0.02881	0.0218	0.1768	0.088201	6.30E-04	0.005108

Notes: Bold indicate the genes with p-value < 0.05 and their mRNA expression trend was opposite to their methylation levels in table 2.

Abbreviations: FC, fold change; FDR, false discovery rate.

Discussion

In the present study, we first mined PCa recurrence-related co-methylation modules using the WGCNA method and then used the intersection between module genes and recurrent-associated DMRs to construct the prognostic signature. Thus, the performance of our signature panel for RFS prediction may be better than that of the 52-gene methylation signature reported by Geybels et al²¹ which only used DMRs between patients with high (8–10) and low (≤6) Gleason scores. In line with this hypothesis, we developed a signature panel of 11 methylated mRNAs (MMP7, SLCO3A1, KCNF1, RFXAP, NTRK3, NAV1, HOXA13, HAS2, CBX2, HIST1H2AJ, SNX4) and 3 methylated lncRNAs (MEG3, DSCR9, and HCP5). This signature was shown to independently predict the 5-year RFS of PCa patients, with AUC of 0.969 for the training TCGA dataset and AUC of 0.811 for the testing E-MTAB-6131, both of which were higher than the predictive accuracy of the combined model of 52-gene methylation signature and clinical features (AUC = 0.78) in the study of Geybels et al.²¹ Similar to the reports by Zeng et al,²³ our lncRNA-mRNA-based methylation signature was demonstrated to have a higher prognostic power than the lncRNA- (AUC = 0.959 vs 0.743; C-index = 0.885 vs 0.723) and mRNA-based signature (AUC =0.959 vs 0.93; C-index = 0.885 vs 0.864). The superiority of our methylation signature to clinical features for RFS prediction was also evidence: TNM and PSA were even not identified to be associated with RFS in the multivariate analysis; patients with the same Gleason score (especially those with Gleason score lower than 8 who were usually considered to have a favorable prognosis in clinic³³) could also be further divided into the high-risk group and the low-risk group. These results suggest that our prognostic model may be a promising biomarker for clinical prediction of RFS in PCa patients.

Of these 14 prognostic mRNAs or lncRNAs, the methylation of 3 RNAs (MEG3, SLCO3A1, and NTRK3) was previously demonstrated to be associated with the prognosis of cancer patients: Zhang et al reported that hypermethylation of MEG3 in plasma was associated with worse RFS and overall survival (OS) in cervical cancer patients.³⁴ Gao et al identified that hypermethylation of MEG3 promoter in retinoblastoma tissues was highly associated with poor survival of retinoblastoma patients.³⁵ Li et al observed that the high methylation rate of MEG3 indicated poor prognosis of breast cancer patients (HR: 2.14).³⁶ DNA hypermethylation in the SLCO3A1 promoter region was detected in a small subset of colorectal cancer patients and in HCT116 and Caco-2 cell lines.³⁷ Higher methylation level of SLCO3A1 (also known as OATP3A1) was associated with shorter post-treatment survival of patients with chronic lymphocytic leukaemia.³⁸ Univariate Cox regression and Kaplan–Meier survival analyses showed that high methylation located in NTRK3 gene was significantly associated with poor prognosis in patients with esophageal squamous cell carcinoma (HR = 1.79).³⁹ In agreement with these studies, we also found that patients with a high methylation level of MEG3 (HR = 1.18), SLCO3A1 (HR = 5.852) and NTRK3 (HR = 4.930) were at a higher risk of having unfavorable RFS.

Usually, high CpG island methylation is associated with silencing of tumor suppressor genes, while low methylation results in elevated expression of proto-oncogenes, which ultimately influence the malignant characteristics of tumor cells. This theory has also been validated for some prognostic genes in our risk scoring model: hypermethylation of MEG3 promoter was highly associated with low MEG3 expression in retinoblastoma tissues.³⁵ The use of DNA methyltransferase inhibitor 5-aza-2-deoxycytidine (5-Aza-dC) significantly increased the expression level of MEG3 in retinoblastoma³⁵ and esophageal cancer.⁴⁰ Hypermethylated MEG3 significantly reduced cancer cell proliferation,^35,41 invasiveness,⁴⁰ and promoted apoptosis³⁵ by depressing MEG3 expression and then activating the activity of the Wnt/β-catenin pathway. The force of 5-Aza-dC³⁶ or pcDNA3.1-MEG3^36,42 slowed down cell viability and migration of breast cancer or PCa cells in vitro. The transcript of SLCO3A1 was also observed to be increased following treatment with 5-Aza-dC.³⁸ NTRK3 methylation silenced NTRK3 expression which induced neoplastic transformation in vitro and tumor growth in vivo; reconstitution of NTRK3 induced apoptosis in colorectal cancers.⁴³ DNA methylation profiling analysis and validation experiments showed that NAV1 was significantly hypomethylated in breast cancers;⁴⁴ while overexpression of members of the neuron navigator gene family was reported to promote invasion and predict poor prognosis.^45,46 The breast tumours that overexpressed CBX2 had a clear reduction in DNA methylation.⁴⁷ Elevated CBX2 expression was correlated with poor clinical outcome in breast cancer⁴⁷ and PCa cohorts.⁴⁸ Depletion of CBX2 abrogated cell viability and induced caspase 3-mediated apoptosis in metastatic PCa cell lines.⁴⁸ In line with these studies, we found a negative association between the methylation level and expression level of MEG3, NTRK3, NAV1, and CBX2 in PCa samples. High methylation of NAV1 and CBX2 was also identified to exert a protective role against recurrence in favor of RFS in PCa patients.

Although the expression of some genes was not directly implied to be methylation-related previously, their roles in cancer progression may indirectly explain the resultant functions of methylation: Bayesian network analysis of GEO microarray datasets found low expressed KCNF1 was a predictor of the risk for site-specific metastasis of breast cancer.⁴⁹ RFXAP expression was relatively lower in pancreatic adenocarcinoma tissues and pancreatic cancer-cell lines than that in normal pancreatic tissues or normal pancreatic ductal epithelial cell line. Its low expression level was positively correlated with high tumor stage and poor survival.⁵⁰ RFXAP silencing was also proved to enhance cell viability of pancreatic adenocarcinoma cells.⁵⁰ Histone family gene HIST1H2AJ was upregulated in gynecological tumors.⁵¹ Similar to these studies, we also identified that KCNF1 and RFXAP were downregulated, while HIST1H2AJ was upregulated in PCa samples compared with controls.

Our dry data analysis revealed that the mRNA expressions of HCP5, HOXA13 and HAS2 were not significant in PCa samples. However, some wet experiments obtained significant results. Hu et al found that HCP5 was highly expressed in PCa tissues.⁵² High expression of HCP5 was positively correlated with the metastasis status of PCa patients.⁵² Knockdown of HCP5 inhibited the proliferation, colony formation and induced apoptosis of PCa cells.⁵² Dong et al reported that HOXA13 expression was sharply increased in carcinoma tissues and independently associated with poor prognosis of PCa patients. Forced expression of HOXA13 promoted the proliferation, migration, and invasion, whereas inhibited the apoptosis of PCa cells.⁵³ HAS2 is a gene encoding the hyaluronan synthase, which is a structural component of extracellular matrices and thus plays important roles in cell proliferation and motility. Inhibition of HAS2 by antisense⁵⁴ or 4-Methylumbelliferone⁵⁵ reduced the growth rate and suppressed invasion and chemotactic motility of PCa cells. The high mRNA expression level was opposite to the low methylation level identified in our study, indirectly indicating that these genes may also be methylation-driven in PCa.

Of the three prognostic lncRNAs, only the functional mechanisms of HCP5 and MEG3 were previously explored in PCa. Proto-oncogenic HCP5 was shown to act as a competing endogenous RNA (ceRNA) to sponge microRNA (miR)‑4656 and then led to the upregulation of miR‑4656 target gene cell migration inducing hyaluronidase 1 (CEMIP),⁵² MEG3 functioned as a ceRNA to modulate miR-9-5p/Quaking protein 5 (QKI-5) axis⁵⁶ or directly bound to EZH2.⁴² However, their functions remain unclear. In this study, we predicted HCP5 and MEG3 may, respectively, be co-expressed with HIST1H2AJ and RFXAP, which may represent novel mechanisms to explain the roles of HCP5 and MEG3 in PCa recurrence. Except the study by Yegnasubramanian et al that showed DSCR9 was hypermethylated in PCa compared to control cells,⁵⁷ no scholars attempted to explore its function mechanism. In this study, we newly predicted that there was an interaction between DSCR9 and CBX2.

Some limitations should be acknowledged in this study. First, this is a study to develop and validate the predictive power of our methylation signature panel for RFS prediction using three independent retrospective data downloaded from publicly available dataset. There was unavoidable bias in these three datasets (such as patient selection), which led to slight differences in the predictive results (ie AUCs in the training cohort were higher than those in the testing dataset; the separation between the high-risk and low-risk categories is better in the training set than in the testing set). In the future, several cohorts of patients in our hospital should be prospectively collected to reconfirm the methylation level of signature genes and their prediction ability for RFS (especially for patients with a subdivided Gleason score of 3+4, 4+3, 4+5, 5+4; or specific postoperative treatments; these were not available in the retrospective TCGA data). Second, the associations between methylation and mRNA expression of all our prognostic genes were only estimated by their expression trend and previous studies in other cancers. Methylation inhibitor 5-aza-dC should be used to treat PCa cells to directly verify whether the expression of these signature genes was methylation-mediated and their roles in the progression of PCa. Third, the co-expression relationship between HCP5/MEG3/DSCR9 and their co-expressed mRNAs should be investigated by immunoprecipitation or biotin-labeled RNA pull-down assays.

Conclusion

In summary, we identified and developed a novel 14-DMRs-based risk score system for predicting RFS in PCa patients. This risk scoring model could independently and accurately classify patients into high-risk and low-risk groups and its prognostic power was higher than the clinical Gleason score. Thus, the signature panel may potentially serve as a promising biomarker for guiding individualized treatment of PCa patients. Although methylation of some genes was predicted to involve PCa recurrence by changing their mRNA expression levels, further in vitro and in vivo studies are needed to confirm our hypothesis.

Data Sharing Statement

All data were downloaded from TCGA (https://portal.gdc.cancer.gov/), GEO (GSE26126; http://www. ncbi.nlm.nih.gov/geo/) and EMBL-EBI (E-MTAB-6131; https://www.ebi.ac.uk).

Ethics Approval and Consent to Participate

No experiments on humans and animals were used in this study. Therefore, no ethical approval is waived.

Disclosure

The authors declare that they have no conflicts of interest in this work.

Additional information

Funding

This work was supported by Science and Technology Innovating and Popularizing Project of Guangdong (No.2019A141405007), Medical Science and Technology Research Fund Project of Guangdong (No. A201926, A2019560, and C2019108).