Investigation of the molecular mechanisms underlying postoperative recurrence in prostate cancer by gene expression profiling.

The present study aimed to identify potential genes associated with prostate cancer (PCa) recurrence following radical prostatectomy (RP) in order to improve the prediction of the prognosis of patients with PCa. The GSE25136 microarray dataset, including 39 recurrent and 40 non-recurrent PCa samples, was downloaded from the Gene Expression Omnibus database. Differentially-expressed genes (DEGs) were identified using limma packages, and the pheatmap package was used to present the DEGs screened using a hierarchical cluster analysis. Furthermore, gene ontology functional enrichment analysis was used to predict the potential functions of the DEGs. Subsequently, Kyoto Encyclopedia of Genes and Genomes pathway enrichment analyses were performed to analyze pathway enrichment of DEGs in the regulatory network. Lastly, a protein-protein interaction (PPI) network of the DEGs was constructed using Cytoscape software to understand the interactions between these DEGs. A total of 708 DEGs were identified in the recurrent and non-recurrent PCa samples. Functional annotation revealed that these DEGs were primarily involved in cell adhesion, negative regulation of growth, and the cyclic adenosine monophosphate and mitogen-activated protein kinase (MAPK) signaling pathways. Furthermore, five key genes, including cluster of differentiation 22, insulin-like growth factor-1, inhibin β A subunit, MAPK kinase 5 and receptor tyrosine kinase like orphan receptor 1, were identified through PPI network analysis. The results of the present study have provided novel ideas for predicting the prognosis of patients with PCa following RP.

Introduction

Prostate cancer (PCa), as one of the most common men's malignancy in America, is the second leading cause of cancer-related death in men (1). Although more than 80% of PCa was diagnosed as localized disease and commonly treated by radical prostatectomy (RP), postoperative recurrence occurred in about 15% of patients within 5 years and up to 40% within 10 years (2). Recurrence of localized PCa following treatment can lead to lethal metastatic castration-resistant PCa. Various biomarkers have been reported for PCa recurrence surveillance, including preoperative prostate specific antigen (PSA) value, Gleason score, lymph node invasion and others, but not cancer-specific and inaccurate (3). Therefore, more efforts should be devoted for identifying disease specific markers of PCa recurrence that can better directly offer practical aid to drug treatment and lead to improved survival and reductions in morbidity.

Although the mechanism underlying PCa is not yet completely understood, multiple genes to help predict PCa risk have been proposed by considerable researches. Brian R. Hu et al (4) reported that AXIN2 expression could not only predict PCa recurrence, but also promoted tumor growth and metastasis in vivo and vitro. Hao et al (5) found that XPO6 expression was elevated in PCa and maybe a potential prognostic biomarker for PCa recurrence. Additionally, some other targets from blood and (or) urine have been examined and identified, including KLK2-KLK3 SNP rs2735839, 17p12 SNP rs4054823 and Eotaxin-1 (6,7). However, few of these profiles have been adopted in the clinic after RP to predict recurrence PCa. Therefore, there is still a need for novel tumor biomarkers that can help improve prediction of prostate cancer recurrence upon clinical variables.

To explore more meaningful molecular biomarker for predicting the prostate cancer prognosis, technologies with high-throughput screen was implied to identify the genes. Microarray data GSE 25136 with 39 recurrent and 40 non-recurrent PCa was published and analyzed by Stephenson et al (8) via leave-one-out-cross-validation (LOOCV) approach and the results showed that Etoposide-induced 2.4 mRNA (EI24) and mitogen-activated protein kinase kinase kinase kinase 4 (MAP4K4) were the most highly overexpressed genes and erythrocyte membrane protein band 4.9 (EPB49) was the most highly underexpressed gene in recurrent tumors compared with primary PCa and may be the potential biomarker. Subsequently, Sun and Goodison (9) conducted a more advanced computational algorithm to analyze the Microarray data GSE25136 and acquire more accurate biomarkers for predicting the prognosis of PCa. With technological development, bioinformatics has been a mainstream tool to analyze the microarray data. In the present study, microarray data GSE25136 ( was employed to identify differentially expressed genes (DEGs) between PCa and PCa recurrence samples with Limma package in R language. Furthermore, gene ontology (GO) and pathway enrichment analysis was performed to screen the DEGs. Lastly, PPI networks of DEGs was constructed by Cytosacpe mapping software and hub genes was identified by the STRING database. Therefore, it is better for us to further understand the molecular mechanisms of PCa.

Materials and methods

Microarray data

The gene expression profiles of GSE25136 were downloaded from the GEO database. GSE25136 based on Affymetrix GPL96 platform (Affymetrix Human Genome U133A Array), was submitted by Sun and Goodison (9) and updated on Jul 01, 2016. The GSE25136 dataset contained 79 PCa samples treated by radical prostatectomy (RP) in 1993 and 1999, including 39 recurrent and 40 non-recurrent PCa samples. When serum level of PSA consecutively increased at least 3 times post operation, the patients were classified as disease recurrence; non-recurrent patients with an undetectable PSA (<0.05 ng/ml) for at least 5 years after RP were identified. The clinical characteristics of all 79 patients has been completely described by Stephenson et al (8). In briefly, Median PSA level and Prostatectomy Gleason sum of patients with recurrence were higher than those in non-recurrenct PCa patients, and the number of patients with Extracapsular extension, positive surgical margins, seminal vesicle invasion were greater in recurrent PCa group compared with non-recurrent PCa group.

Identification of DEGs

The raw data files used for the analysis included cell files (Affymertix platform). The data was preprocessed by R biocondutor RMA Packages, and DEGs were identified by limma packages in recurrent PCa compared with non-recurrent PCa samples. DEGs were identified with a change fold and defined a P-value cutoff of <0.05 to be statistically significant. Hierarchical clustering analysis was applied to categorize the data into two groups that had similar expression patterns. Heatmap was performed by the pheatmap package analysis with joint between-within distances. Expression values from multiple clones or probe sets mapping to the same Unigene Cluster ID were averaged.

Gene ontology (GO) analysis of DEGs

The Database for Annotation, Visualization and Integrated Discovery (DAVID; https://david.ncifcrf.gov/) provides a comprehensive set of novel and powerful tools for assigning biological meaning to a set of genes (10). The false recovery rate (FDR) <0.05 was used as the cut-off criterion for GO functional enrichment analysis by DAVID.

Pathway enrichment analysis of DEGs in the regulatory network

KEGG (http://www.genome.jp/) is acknowledge base for systematic analysis of gene functions, linking genomic information with higher-order functional information (11). P<0.05 was used as the cut-off criterion for the Kyoto Encyclopedia of Genes and Genomes pathway enrichment analysis using DAVID.

Integration of protein-protein interaction (PPI) network and module analysis

The search Tool for the Retrieval of Interacting Genes (STRING) database is an online tool designed to evaluate PPI information. STRING (version 10.0) covers 5,214,234 proteins from 1,133 organisms. To evaluate the interactive relationships between the DEGs, the DEGs were mapped to STRING, and only experimentally validated interactions with a combined score >0.4 were selected as significant. Then, the Cytoscape software was used to construct the PPI network. The plug-in Molecular Complex Detection (MCODE) was used to screen the modules of PPI network in Cytoscape with a degree cutoff=2, node score cutoff=0.2, k-core=2, max depth from seed=100. The criteria were set as follows: MCODE scores >4 and number of nodes >4. P<0.05 was considered to have significant differences.

Results

After data, including 39 recurrent PCa samples and 40 non-recurrent PCa samples, was downloaded from GEO database and preprocessed, 708 DEGs, including 212 up genes and 496 down genes were identified using limma packages on the basis of the cut-off criteria (P<0.05 and fold control (FC) ≥1.4 criteria) in recurrent samples compared with non-recurrent samples. Subsequently, DEGs were performed hierarchical clustering analysis and it can accurately classify the prostate samples as recurrent PCa tissues and non-recurrent PCa tissues (Fig. 1: left, recurrent PCa; right, non-recurrent PCa). Additionally, Top up 50 DEGs and down 50 DEGs expression with most significant was shown in Fig. 1 (P<0.05).

Figure 1.

Identification of DEGs between recurrent and non-recurrent prostate cancer tissues by Hierarchical cluster analysis. (Left): Non-recurrent PCa group, (Right): Recurrent PCa group. Each row represents a single gene; each column represents a tissue sample. The gradual color change from green to red represents the changing process from dowregulation to upregulation.

GO functional enrichment analysis

In order to gain further insight into the function of the identified DEGs, we uploaded DEGs to the online biological classification software DAVID to identify typical GO categories. Go analysis showed that DEGs were significantly enriched in biological processes (BP), including cell adhesion, negative regulation of growth, extracellular matrix organization, negative regulation of cell migration, apoptotic signaling pathway (Table I). For cell components, DEGs were enriched in focal adhesion, extracellular exosome, cell-cell adherens junction, proteinaceous extracellular matrix (Table I). Finally, Go molecular function analysis showed that DEGs were enriched in protein binding, protein homodimerization activity, cadherin binding involved in cell-cell adhesion, insulin receptor binding (Table I).

Table I.

Gene ontology analysis of differentially expressed genes associated with PCa recurrence.

Category	Term/gene function	Gene count	%	P-value
GOTERM_BP_DIRECT	GO:0007155~cell adhesion	38	6.551724138	3.82E-07
GOTERM_BP_DIRECT	GO:0045926~negative regulation of growth	8	1.379310345	1.33E-06
GOTERM_BP_DIRECT	GO:0030198~extracellular matrix organization	22	3.793103448	1.72E-06
GOTERM_BP_DIRECT	GO:0030336~negative regulation of cell migration	12	2.06896552	2.57E-04
GOTERM_BP_DIRECT	GO:0097190~apoptotic signaling pathway	10	1.72413793	4.71E-04
GOTERM_CC_DIRECT	GO:0005925~focal adhesion	42	7.241379	8.83E-12
GOTERM_CC_DIRECT	GO:0070062~extracellular exosome	149	25.68966	1.26E-11
GOTERM_CC_DIRECT	GO:0005913~cell-cell adherens junction	23	3.965517241	4.87E-04
GOTERM_CC_DIRECT	GO:0005578~proteinaceous extracellular matrix	23	3.965517	3.33E-05
GOTERM_MF_DIRECT	GO:0005515~protein binding	372	64.13793103	3.67E-15
GOTERM_MF_DIRECT	GO:0042803~protein homodimerization activity	46	7.931034483	2.34E-05
GOTERM_MF_DIRECT	GO:0098641~cadherin binding involved in cell-cell adhesion	20	3.448275862	0.002756127
GOTERM_MF_DIRECT	GO:0005158~insulin receptor binding	6	1.034482759	0.0028733

BP, biological process; CC, cell component; MF, molecular function.

KEGG pathway analysis

We employed KEGG pathway enrichment analysis to identify the most significantly enriched pathways of the DEGs. 10 biological pathways which significant enriched with DEGs including cAMP signaling pathway, MAPK signaling pathway, Adherensjunction, Calcium signaling pathway, Pathways in cancer, Proteoglycans in cancer, Transcriptional misregulation in cancer, Leukocyte transendothelial migration, Focal adhesion and Ras signaling pathway (Table II).

Table II.

KEGG pathway analysis of differentially expressed genes associated with PCa recurrence.

Pathway ID	Name	Gene count	%	P-value	Genes
hsa04024	cAMP signaling pathway	23	0.021929197	5.74E-05	PLD1,VAV3,PTGER3,PDE3B, GRIA3, PDE4D, GABBR2, CNGB1, BDNF, HTR1A,ATP2B4, NPY, GRIA2, PDE4A, RAC1, CREB3L2, RYR2, GNAS, CAMK2B, ADCY10, FSHB, GLP1R, NFATC1
hsa04010	MAPK signaling pathway	22	0.020975754	0.004317491	FGFR2, FGF8, FGF7, CACNA1I, TAOK3, MAPK11, MECOM, FLNC, FLNB, CDC42, MAP4K4, CASP3, BDNF, RPS6KA4, PAK2, SOS1, RAC1, CACNA1G, EGF, DUSP7, MAP2K5, NFATC1
hsa04520	Adherens junction	9	0.00858099	0.012675942	ACTB, CDC42, TCF7, CSNK2A1, BAIAP2, RAC1, SSX2IP, PTPN1, ACTN3
hsa04020	Calcium signaling pathway	16	0.015255094	0.012798992	SLC8A2, PTGER3, SPHK2, SPHK1, CACNA1I, VDAC1, GNAL, ATP2B4, ATP2A3, RYR3, PDE1A, CACNA1G, RYR2, GNAS, CAMK2B, ADRA1D
hsa05200	Pathways in cancer	27	0.025742971	0.025195199	FGFR2, FGF8, TCF7, PTGER3, CTBP2, FGF7, RXRB, ITGA2, IGF1, STAT1, MECOM, CDC42, CASP3, LAMB2, HDAC2, CXCR4, SOS1, RAC1, MDM2, NKX3-1, GNAS, GNB3, GNG3, RARB, EGF, WNT6, APC
hsa05205	Proteoglycans in cancer	16	0.015255094	0.031484838	ACTB, PPP1R12B, ITGA2, IGF1, MAPK11, FLNC, FLNB, KDR, CDC42, CASP3, HPSE, SOS1, RAC1, MDM2, CAMK2B, WNT6
hsa05202	Transcriptional misregulation in cancer	14	0.013348207	0.035045912	SUPT3H, FLT1, RXRB, IGF1, PAX5, GRIA3, GZMB, HDAC2, REL, LYL1, MDM2, PBX1, IGFBP3, CDK14
hsa04670	Leukocyte transendothelial migration	11	0.010487877	0.035804774	ACTB, CDC42, ICAM1, NOX3, VAV3, CXCR4, CLDN5, NOX1, RAC1, MAPK11, ACTN3
hsa04510	Focal adhesion	16	0.015255094	0.039374107	ACTB, CDC42, FLT1, VAV3, LAMB2, PAK2, SOS1, PPP1R12B, RAC1, IGF1, ITGA2, ACTN3, FLNC, EGF, FLNB, KDR
hsa04014	Ras signaling pathway	17	0.016208537	0.042199384	FGFR2, PLD1, FGF8, FGF7, FLT1, IGF1, BRAP, KDR, CDC42, PAK2, REL, SOS1, TEK, RAC1, GNG3, GNB3, EGF

FDR<0.05.

Construction of the PPI network

Cytosacpe mapping software was employed to construct the PPI network of DEGs. A total of 663 nodes and 8,871 edges were analyzed using plug-ins MCODE. The top 5 significant modules with MCODE scores >4 and nodes >4 in whole network were screened by analysis in the STRING database, and the hub gene in each cluster, also called the seed, was identified by on the basis of the highest modules scoring in the cluster including Insulin-like growth factor-1 (IGF-1) (Fig. 2A), mitogen-activated protein kinase kinase 5 (MAP2K5) (Fig. 2B), Receptor tyrosine kinase like orphan receptor 1 (ROR1) (Fig. 2C), Inhibin beta A (INHBA) (Fig. 2D), and differentiation-22 (CD22) (Fig. 2E). All clusters are named after the hub gene name, and of these clusters, IGF-1 modules showed a highest MCODE scores in whole network, with 20.606. Additionally, IGF-1 modules consisted of 34 nodes and 340 edges; MAP2K5 modules consisted of 43 nodes and 306 edges; ROR1 modules consisted of 71 nodes and 464 edges; INHBA modules consisted of 51 nodes and 119 edges; CD22 modules consisted of 38 nodes and 74 edges.

Figure 2.

PPI sub-network of hubgenes. (A) IGF-1 modules with 34 nodes and 340 edges; (B) MAP2K5 modules with 43 nodes and 306 edges; (C) ROR1 modules with 71 nodes and 464 edges; (D) INHBA modules with 51 nodes and 119 edges; (E) CD22 modules with 38 nodes and 74 edges in whole network.

Discussion

PCa is the fourth leading global cause of human malignancies worldwide, and is a product of mutation in genomics including cumulative genetic, epigenetic, somatic, and endocrine aberrations (12). Of the differential expression of genes caused by various mutations, some specific genes are positively or negatively associated with therapy resistance and poor outcomes in PCa. The wide application of microarray and high through put sequencing has made it possible to identify the more appropriate genes to predict the prognosis of PCa after RP from thousands of genes in human genome level (12). In the present study, we extracted the data from GSE25136 and 708 DEGs between recurrent PCa samples and non-recurrent PCa samples using bioinformatics analysis were screened out. Functional annotation showed that these DEGs were mainly involved in cell adhesion, focal adhesion, protein binding, cAMP signaling pathway and MAPK signaling pathway. In addition, to better understand the interaction of these DEGs, a PPI network was constructed and we identified four key genes, including CD22, IGF-1, INHBA, MAP2K5 and ROR1, that can provide new ideas for predict the prognosis in PCa following RP.

The GO term analysis showed that these DEGs were mainly involved cell adhesion, focal adhesion, and protein binding. In addition, cAMP and MAPK signaling pathway were shown to participate int PCa recurrence by KEGG pathway analysis. Classical signal transduction pathway and cAMP signaling pathway have been extensively studied in the context of carcinogenesis by regulating cellular growth and proliferation. cAMP-dependent protein kinase (PKA), as a critical mediator of cAMP signaling pathway, has been demonstrated that it is overexpressed in PCa and has been examined as a potential biomarker for predicting the outcome of PCa patients (13). Androgens are required for the initiation and the development of PCa via stimulating the AR signaling pathway, and androgen ablation therapies, such as chemical or surgical castration, have become a standard against PCa (14). There is a highly relevant cross-talk between cAMP and AR signaling pathway in PCa progression, because not only cAMP and PKA activation may result in the stimulation of AR but androgens can also regulate the activity of PKA (15). In addition, cyclic nucleotide phosphodiesterases (PDEs) are involved in the metabolism of cAMP by regulating its degradation (16). PDE4D, as a kind of PDEs, is highly expressed in PCa and has been implicated to promote PCa progression (16). Furthermore, members of the PDE4D subfamily are classified as long, short and super-short. PDE4D7, as a long isoform member, is downregulated in androgen-independent prostate cancer cells compared with androgen sensitive prostate cancer cells, and inhibited its growth by compartmentalising cAMP (17). R Böttcher et al (18) also showed that it was up-regulated in localized PCa samples compared with the normal adjacent prostate tissues, while its expression diminished with emergence of Castration resistant prostate cancer (CRPC). Other PDE4D isoform composition, such as PDE4D5 and PDE4D9, was also upregulated in PCa and played an import role in PCa progression (19). MAPK signaling pathway also played a vital role in regulating cellular behaviors in response to extracellular stimuli. Dysregulation of p38 MAPK, as a main subgroup of MAPK signaling pathway, are associated with tumor stages and poor survival of PCa patients (20). The emergence of Castration resistant prostate cancer (CRPC) caused by certain co-activators or through MAPK signaling pathway activities, which lead to the overexpression of anti-apoptotic genes and survival of the cancer cells, thus increasing the PCa related death (21).

Finally, the PPI network with DEGs was constructed and the hub genes exhibiting the highest degree of connectivity were identified, including CD22, IGF-1, INHBA, MAP2K5 and ROR1. IGF-1 was identified as one of the most DEGs in the recurrent PCa samples. IGF-1, also known as somatomedin 1, is a mitogen that plays a key role in regulating various cell biological behavior, including cell proliferation, differentiation, and apoptosis via endocrine, paracrine and autocrine mechanisms (22). IGF-1 binds to the insulin-like growth factor 1 recep-tor (IGF-1R), which is a tyrosine kinase receptor, and initiates a cascade of downstream signal transduction pathways (23). Results published recent studies evaluated the role of IGF-1 in PCa and showed that higher circulating IGF-1 levels were consistently associated with increased risk of PCa in epidemio-logic studies (24). IGF-1, which is synthesized locally in an autocrine or paracrine manner by PCa cells, may stimulate PCa growth and development (25). Then, IGF has been a pivotal target gene for PCa therapy. Magnolol has been demonstrated that it served as a novel anti-PCa agenet via regulating the expression of IGF-1 in vitro (26). Apigenin effectively suppressed PCa cells growth and metastasis in TRAMP mice by attenuating IGF-I signaling (27). In addition, IGF-1 genotypes and haplotypes were associated with wore survival of PCa patients with bone metastasis (28). Mitogen-activated protein kinase kinase 5 (MAP2K5), also known as MEK5, was overexpressed in PCa, which was associated with tumor metastases and unfavourable survival outcome of PCa patients (29).

Receptor tyrosine kinase like orphan receptor 1 (ROR1), also known as neurotrophic tyrosine kinase, receptor-related 1 (NTRKR1), is a transmembrane protein belonging to the receptor tyrosine kinase (RTK) family. Down-regulation of ROR1 inhibited human colorectal cancer cell growth and promoted apoptosis (30). ROR1 has been shown to be overexpressed in several solid tumors and its unique expression by malignant cells surface is a target for novel therapeutics, especially monoclonal antibodies (mAbs) for the treatment of cancer (31). Inhibin beta A, also known as INHBA, is a subunit of both activin and inhibin, two closely related glycoproteins with opposing biological effects. INHBA is overexpressed in various cancers, including gastric cancer, rothelial carcinoma of the urinary bladder and upper tract and colorectal cancer, demonstrating its association with poor prognosis of patients (32–34). Cluster of differentiation-22 (CD22), as a molecule belonging to the SIGLEC family of lectins, is a transmembrane glycoprotein expressed by mature B cells. Recently, Tuscano and colleagues reported that CD22, a hallmark marker on B lymphocytes, was expressed on lung cancer cells and might serve as a new target for therapy (35). However, Pop et al (36) reported that the surface of lung cancer cells did not detect CD22 expression, and cannot be killed by anti-CD22 immunotoxins, which have not previously been directly associated with initiation and progression of PCa, according to the present results.

Although data from GSE25136 have been analyzed by previous authors, and some DEGs have been identified, the uniqueness of the present study is that Limma package in R language, as one of the most fashionable Algorithmic Language now, was applied to analyze the data from GSE25136 and different DEGs compared with previous report was identified. Furthermore, DEGs in PCa recurrence related BP and signaling pathway was screened out, which may help us better understand the potential mechanism of PCa recurrence. Additionally, a PPI network of DEGs was constructed and 5 hub nodes with higher degrees were identified and could be used to predict the prognosis of PCa.

In conclusion, the results of this study provide a comprehensive bioinformatics analysis of DEGs to increase the understading of the mechanism underlying PCa recurrence. The study showed that CD22, IGF-1, INHBA, MAP2K5 and ROR1 may be pivotal for participating in PCa recurrences. However, these functions need to be confirmed by further molecular biological experiments.