Reconstitution of the ERG Gene Expression Network Reveals New Biomarkers and Therapeutic Targets in ERG Positive Prostate Tumors.

BACKGROUND: Despite a growing number of studies evaluating cancer of prostate (CaP) specific gene alterations, oncogenic activation of the ETS Related Gene (ERG) by gene fusions remains the most validated cancer gene alteration in CaP. Prevalent gene fusions have been described between the ERG gene and promoter upstream sequences of androgen-inducible genes, predominantly TMPRSS2 (transmembrane protease serine 2). Despite the extensive evaluations of ERG genomic rearrangements, fusion transcripts and the ERG oncoprotein, the prognostic value of ERG remains to be better understood. Using gene expression dataset from matched prostate tumor and normal epithelial cells from an 80 GeneChip experiment examining 40 tumors and their matching normal pairs in 40 patients with known ERG status, we conducted a cancer signaling-focused functional analysis of prostatic carcinoma representing moderate and aggressive cancers stratified by ERG expression.
RESULTS: In the present study of matched pairs of laser capture microdissected normal epithelial cells and well-to-moderately differentiated tumor epithelial cells with known ERG gene expression status from 20 patients with localized prostate cancer, we have discovered novel ERG associated biochemical networks.
CONCLUSIONS: Using causal network reconstruction methods, we have identified three major signaling pathways related to MAPK/PI3K cascade that may indeed contribute synergistically to the ERG dependent tumor development. Moreover, the key components of these pathways have potential as biomarkers and therapeutic target for ERG positive prostate tumors.

Background

Cancer of the prostate (CaP) is the second leading cause of cancer related death in men in the USA with estimated 238,590 new cases and 29,720 deaths in 2013 1. Although significant progress has been made in the early detection and treatment of the organ confined CaP, better understanding of the cancer biology is essential in addressing current challenges related to the development of more precise diagnostic and prognostic biomarkers and rational therapy targets. 2 Recent investigations have been providing new insights into common genomic alterations in CaP. 3 The driver genetic alterations and molecular mechanisms contributing to CaP onset and progression need to be better understood. 4-8 Despite a growing number of studies evaluating CaP specific gene alterations, oncogenic activation of the ETS Related Gene (ERG) by gene fusions remains the most validated cancer gene alteration in CaP. 9-11

Evaluations of the CaP epithelial transcriptome revealed overexpression of the ERG mRNA in prostate cancer cells. 12-14 Further, our quantitative evaluations of ERG mRNA expression in laser capture micro-dissected (LCM) matched benign and tumor cells of 114 patients, identified ERG overexpression in over two thirds of patients. 14 Prevalent gene fusions have been described between the ERG gene (or other members of the ETS Family of transcription factors) and promoter upstream sequences of androgen-inducible genes, predominantly TMPRSS2 (transmembrane protease serine 2) and other androgen inducible genes. 15-17 ERG knock-down in TMPRSS2-ERG positive CaP cell culture and xenograft models revealed significantly reduced tumor cell growth. 18 Further, altered ERG expression in cell culture models of CaP resulted in differentiation-associated or invasion-associated gene expression. 18-20

With the development of ERG monoclonal antibodies, recent reports provided first insights into the detection of the protein products of the ERG gene fusions in CaP. 21, 22 The ERG oncoprotein expression in CaP is highly tumor cell specific (>99%) and concordant (>95%) with the ERG gene fusion status. 23-26 Despite the extensive evaluations of ERG genomic rearrangements, fusion transcripts and now the ERG oncoprotein, the prognostic value of ERG remains to be better understood. 11, 27 Cooperation of ERG with other CaP-associated alterations, primarily the loss of PTEN (phosphatase and tensin homolog) has shown association with the disease progression. 28-33

Studies have also suggested prognostic implication of the relative abundance of specific TMPRSS-ERG splice variants (full length vs. truncated) in CaP progression. 20, 34 Although, a considerable number of studies have addressed the clinical correlations of the oncogenic activation of ERG, critical pathways affected by ERG remain to be better defined. Since the analysis of pairwise gene interaction of fusion gene networks, including ERG, indicated universal principles in human neoplasia 35, bioinformatic approaches have been used to address the ERG-regulated gene expression signature in CaP. Studies on ERG-associated gene expression patterns revealed signatures of epigenetic reprogramming and elevated expression of WNT-associated and attenuated cell death pathways 36 Association of epigenetic alterations as well as the presence of WNT signature in ERG positive human prostate tumors have been confirmed by several studies 37-39 A functional interface of ERG and C-MYC (v-myc myelocytomatosis viral oncogene homolog (avian)) has also been noted in the regulation prostate epithelial differentiation markers 18 Integrative molecular concept analysis addressing the CaP-associated gene expression signatures have revealed ERG as the most upregulated feature in prostatic adenocarcinoma. These data also highlighted the enrichment of differentially expressed networks of the protein biosynthesis pathway and androgen signaling in the transition from benign to pre-invasive and from localized to metastatic prostate cancer 40

Using gene expression dataset from matched prostate tumor and normal epithelial cells from 40 patients 14, 18, 41; we conducted a cancer signaling-focused functional analysis of prostatic carcinoma representing moderate and aggressive cancers stratified by ERG expression. Our findings suggest the activation of MAPK (mitogen activated kinase-like protein) and PI3K (Phosphoinositide 3-kinase) cascades by synergistic actions of several genes (FGF2 (fibroblast growth factor 2), Angiotensin II, Neprilysin, AR (androgen receptor), E2F1 (E2F transcription factor 1), Tcf (T-cell specific factors), and DAN (neuroblastoma, suppression of tumorigenicity 1)) via three major distinctive signaling pathways.

Methods

Selection of Cases, LCM isolation of Tumor and Matching Benign Cells, Detection of ERG

Using our previously reported gene expression dataset (GSE32448) from laser capture micro-dissected matched normal and tumor cells from 40 patients 18, 41; we performed evaluation of biochemical pathways in ERG positive and ERG negative tumors. Briefly, from a PSA (prostate specific antigen) screened patient group with no prior neoadjuvant treatment, 300 cases were examined under an approved protocol from the Institutional Review Board approval of Walter Reed National Military Medical Center (WRNMMC) and Uniformed Services University of the Health Sciences (USUHS). Following radical prostatectomy frozen tumor tissues were obtained by ex-vivo biopsy of palpable tumors representing mainly index tumors. Prostates were preserved as whole-mounted FFPE fixed specimens. Following the evaluation of all tumor foci in whole-mounted sections, tumor growth pattern (Gleason score) and also the glandular differentiation (differentiation) were recorded in both the FFPE and the matching frozen tissues. We have selected cases representing two extremes from the continuum of prostate cancers. In one group (n=20) cases were selected with Gleason sum 6 (n=11) and 7 (3+4, n=9) and tumors cells with well- to moderately differentiated (WMD) morphology. The laser capture procedure was performed on the single-cell level capturing only cells with well to moderately differentiated morphology from Gleason patter 3. In the second group (n=20) cases with Gleason sum of 8-9 were selected. Single cell resolution of laser capture was directed to cells with PD morphology. Cells were isolated by laser capture microdissection (LCM) (Arcturus Pixel 2) from OCT embedded and H&E stained frozen prostate sections (2000 laser shots per area as described earlier 14, 41-43 The morphological homogeneity was approximately 90% for cancer or benign epithelial cells. As a quality control for tumor or benign epithelial cells ERG, AMACR (alpha-methylacyl-CoA racemase) and PCA3 (prostate cancer antigen 3) transcripts were analyzed by QRT-PCR, confirming the morphological assessment 43 Presence or absence of ERG expression, TMPRSS2-ERG transcripts, and the expression of major ERG splice variants were derived from the previously published data 14, 34, 41 TMPRSS2-ERG fusion junctions (fusion “a, b and/or c”) at mRNA levels was detected in 18 of the 20 samples with WMD tumor cell morphology. TMPRSS2-ERG fusion junction was found only in 9 cases among the samples from PD tumor cells as previously shown 41.

Microarray Expression Profiling and Statistical Analysis

RNA extracted from the LCM benign and malignant epithelial cells and were hybridized to Affymetrix GeneChips (Human Genome U133 Plus 2.0 chip) and processed as decribed 14, 42 Expression profiles in CEL format containing probe intensities corresponding to 80 samples were [NCBI/GEO: GSE32448] processed using R (http://r-project.com/). Additional modules used for data processing, normalization, clustering were obtained from Bioconductor Project (http://www.bioconductor.org) 44. The initial processing of CEL files was performed using an updated GeneChip Robust Multiarray Averaging (GC-RMA) algorithm with new probe annotations 45, 46 for EntrezGene features provided by Affymetrix (http://brainarray.mbni. med.umich.edu/ Brainarray /DatabaseCustomCDF /genomic_curated_CDF.asp). For the purpose of this study, each set of intensities obtained from EntrezGene microarray analysis is referred to as an “experiment” file.

The data for present study thus consisted of experiments obtained for prostate tumor cells and matched benign prostate cells from the same cases. For each EntrezGene feature the change in gene expression level was calculated by two ways: matched fold changes and mean fold changes. For each matched pair of experiments the change in gene expression level was calculated as the difference between gene expression levels in tumor sample vs. normal sample. This approach takes into account variability in each individual's gene expression pattern. Mean fold changes were computed for each EntrezGene feature as difference between mean expression in case matched samples and mean expression in control samples for each defined comparison groups of samples.

The following datasets were subjected to clustering analysis: (a) full set of 40 cases with Fold-Changes, (b) a set of only WMD tumor cells (Gleason score ≤7), (c) a set of only poor differentiated tumor cells (Gleason score > 7).

Clustering Analysis

Clustering analysis of experiments in a fold change format was performed with Pvclust 47. This package allows for hierarchical clustering using various algorithms with Multiscale Bootstrap Resampling (MRB). The method is based on resampling of data, and applicable to a large class of problems including hierarchical clustering. In the MRB, we set the gene set size of bootstrap samples to several values. Let N be the original data size, and N' be that for bootstrap samples (N'= N*0.5; 0.6; 0.7; 0.8; 0.9; 1; 1.1; 1.2; 1.3 and 1.4). In our case, N = 18873 EntrezGene features, and N' = 9436, 11323, 13211, 15098, 16985, 18873, 20760, 22647, 24534 and 26422 random EntrezGene features. MRB uses 1000 bootstrap replicates with random N' data sizes. On the basis of constructed bootstrap samples MRB measures the accuracy of defined clusters as p-values, which ranges from 0 to 1. MRB yields Approximatelly Unbiased (AU) p-value for each cluster, a cluster is considered statistically significant if AU p-value < 0.05. High statistical approval of AU p-value usage during cluster analysis is described in 48. In our work we use standard R-package representation of these p-values, therefore statistically significant range of AU p-values is from 95 to 100 47.

Clustering analysis was performed with Single Linkage, Complete Linkage, Average Linkage and Centroid Linkage algorithms. To determine the distance between clusters, Euclid metrics was used.

Clustering analysis was performed on a total set of EntrezGene features for all experiments, as well as on partial sets 49 Specifically, 1000, 3000, and 5000 expression profiles with the biggest standard deviations between log2 transformed fold changes of expression level in tumor tissue compared to normal tissue were used for clustering analysis first, and the cluster stability was then tested on a full set of genes (all features found on the investigated Affymetrix platform).

Identification of differentially expressed genes

To identify differentially expressed genes between defined groups of samples, we used Bayesian moderated t-statistics coupled with false discovery rate (FDR) adjustment = 0.05. This test borrows information about gene-wise signal variances across all genes and uses empirical Bayes approach to derive posterior estimates of variance. Formally it shrinks variances towards a pooled estimate, thus increasing statistical power when number of samples is small. Moreover, this test is insensitive to missing data and doesn't require normality of the data 50

From the lists of differentially expressed genes, the genes (EntrezGene features) were selected based on expression levels: either downregulated or upregulated at least 1.5 mean fold change values.

Functional Analysis

We used the knowledge-based MetaCore 6.7 software package 51 to identify molecular processes that are affected by the gene's products of translation and transcription. MetaCore™ analyzes the global biological network by various algorithms, such as direct interactions, shortest paths, analyze network, transcription regulation, and auto-expand, expand by one interaction, and others. 51 We used transcription regulation algorithm to identify the most overconnected transcription factors in a set of differentially expressed genes in our experimental data and Canonical Pathway modeling to identify canonical pathways enriched with differentially expressed genes from our data.

Transcription regulation algorithm starts building a network by first constructing a smaller network that consists of the initial set of genes of interest and then by adding closest transcription factors (transcription factors that are connected to a gene in a set by “transcription regulation” edge). Then, a separate network is built around every transcription factor, using the Auto Expand algorithm (which gradually expands sub-networks around every transcription factor) while limiting the expansion to the objects from the initial list. The algorithm delivers a list of networks, one per transcription factor.

For every initial object the Canonical Pathway modeling algorithm searches the MetaCore™ Database for all canonical pathways that start at, end at, or pass through that initial object, and builds a network consisting of all objects and links from all such pathways. The set of networks thus obtained is then prioritized by the number of initial objects they contain. If two or more networks contain the same number of the initial objects, they are further prioritized by giving preference to those with the least total number of objects.

Statistically significant network modules were combined for easier visualization of interconnecting processes taking place in prostate tumors with the expression of ERG oncogene.

Enrichment Analysis

A total of four functional ontologies of biological processes were used: GO processes (www.geneontology.org) and three MetaCore ontologies - GeneGo Process Networks Ontology (163 networks of cell processes created from literature curated data), GeneGo Pathway Maps (163 maps of main cellular activation processes, also constructed from annually curated data), Diseases (by Biomarkers) (known biomarkers of various diseases). In all three MetaCore ontologies molecular processes are prioritized based on statistical significance. To identify common molecular-biological processes the Compare Experiments workflow was used. This workflow determines common processes based on the analysis of statistically significant biological processes and expression data for genes participating in these processes. Visualization of generated hypotheses in the form of networks was also performed in MetaCore 51.

Results

In order to evaluate the biological role of elevated ERG expression in prostate tumors with ERG activation, we investigated global expression profiles of prostatectomy specimens from PSA (prostate specific antigen)-screened patients, who had no prior androgen ablation treatment. ERG expression was evaluated in human prostate cancer samples by using Affymetrix oligonucleotide microarrays. The data was generated from quality controlled laser capture microdissected tumor cells (n=40) and matched benign epithelial cells (n=40) from OCT embedded frozen prostate specimens. The specimens were obtained from age and race (Caucasians) matched cases, selected from a cohort of over 300 patients. Table 1 summarizes ERG expression levels in each tumor - normal pair sample. Applying a 1.5X cutoff, ERG overexpression is prevalent in well-to-moderately differentiated tumors (tumors that are associated with lower, ≤7 Gleason score, Table 1): 15 out of 20 well-to-moderately differentiated (WMD) tumors had increased ERG levels vs 6 out of 20 PD samples (p-value = 0.04953, Chi-squared test).

Table 1

Tumor cell differentiation and ERG expression levels in tumor samples. ERG+: samples with overexpressed ERG; ERGn: samples with relative ERG expression changes not exceeding the 1.5X cutoff.

Patient#	Diff. status	ERG Fold Change value	ERG mRNA		Patient ID	Diff. status	ERG Fold Change value	ERG mRNA
17	WMD	16.3	ERG+		31	PD	8.9	ERG+
19	WMD	15.2	ERG+		28	PD	8.7	ERG+
18	WMD	12.9	ERG+		26	PD	8.6	ERG+
3	WMD	12.5	ERG+		24	PD	7.9	ERG+
15	WMD	10.2	ERG+		33	PD	7.6	ERG+
1	WMD	9.7	ERG+		37	PD	5.4	ERG+
8	WMD	9.0	ERG+		40	PD	1.5	ERGn
4	WMD	8.3	ERG+		29	PD	1	ERGn
5	WMD	8.1	ERG+		32	PD	1	ERGn
16	WMD	7.7	ERG+		25	PD	1	ERGn
10	WMD	7.2	ERG+		27	PD	1	ERGn
13	WMD	6.7	ERG+		21	PD	1	ERGn
14	WMD	5.9	ERG+		30	PD	1	ERGn
11	WMD	4.3	ERG+		36	PD	1	ERGn
6	WMD	2.6	ERG+		38	PD	1	ERGn
7	WMD	1.1	ERGn		34	PD	1	ERGn
2	WMD	1	ERGn		39	PD	1	ERGn
9	WMD	1	ERGn		35	PD	1	ERGn
12	WMD	1	ERGn		23	PD	1	ERGn
20	WMD	1	ERGn		22	PD	1	ERGn

WMD tumors associate with better overall clinical prognosis (see, for example, 2). As noted before 41, majority (75%) of WMD samples harbored ERG overexpression. Hierarchical clustering analysis of our dataset may help validate the observation that ERG overexpression is more frequently observed in WMD samples. Thus, clustering analysis was performed with Single Linkage, Complete Linkage, Average Linkage and Centroid Linkage algorithms with multiscale bootstrap resampling (MRB). To determine the distance between clusters, Euclid metrics was used. On the basis of constructed bootstrap samples, MRB measured the accuracy of defined clusters as p-values, which ranged from 0 to100. MRB yielded Approximatelly Unbiased (AU) p-value for each cluster; a cluster was considered statistically significant if AU p-value >95. For clustering analysis, we used a total set of genes, about 20,000 (please refer to Methods for explanation how these genes were determined), as well as sets of 1000, 3000 and 5000 genes with the largest standard deviations between log2 transformed fold changes in normal vs. tumor samples expression profiles.

In the analyzed sets two extremes (well-to-moderately and poorly differentiated (PD)) tumors were equally represented. Clustering of all 40 samples with either full or truncated gene set did not yield clusters that correlated with either an observed differentiation status or ERG expression feature (data not shown). In contrast, when we performed clustering analysis on a set of 20 WMD samples with either a full or any of the truncated set of genes we observed a statistically significant partitioning that correlated with ERG+ and ERG- groups (Figure 1). This correlation was not observed in PD samples likely reflecting progression-associated molecular complexity of PD tumor cells.

Figure 1

Clustering of 20 samples from well differentiated tumors with the whole set of differentially expressed gene (DEGs). ERGn samples are in orange frames, ERG+ samples are in purple frames. For both clusters Approximately Unbiased (AU) p-value=96. Values on the edges of the clustering are p-values (%). Red values are AU p-values, and green values are Bootstrap Probability values (Explanations about AU and BP values are found on the website http://www.is.titech.ac.jp/~shimo/prog/pvclust/ along with the description of pvclust R package that was used to perform clustering. In brief, BP values are computed by bootstrap resampling and AU values are computed by multiscale bootstrap resampling and thus are considered a better representation of an unbiased p-value for clusters.) Clusters strongly supported by the data have AU values of larger than 95% and are highlighted in red in the Figure.

Hierarchical clustering analysis confirmed that the set of genes expressed differentially in normal and cancerous tissue, partitioned samples into distinct clusters with ERG expression in WMD and not in PD prostate cancer samples. This result is consistent with several previously reported observations for this dataset 14, 41 indicating that elevated ERG expression in WMD samples suggests for the role of ERG oncogene in early stages of CaP onset.

We next explored the mechanism of ERG overexpression in the subset of WMD cancer samples by building a “causal” model, or network, using differentially expressed genes (DEGs) as the seed nodes and collection of some 400,000 protein-protein interactions in MetaBase™ (Thomson Reuters) as the source of edges 51 DEGs used for causal network construction were identified by comparing expression profiles of ERG-overexpressing, WMD tumor samples with normal tissues from matched patients (Figure 2). A total of 207 genes (for the list, see Addition file 1) with > 1.5 fold change and p-value of differential expression <0.05, were processed for further analysis. However, comparison between ERG positive or ERG negative tumors (regardless of differentiation status) without reference to matched normal epithelial cells did not yield significant enrichment of DEGs to perform further analysis (Additional file 2).

Figure 2

Flowchart of sample selection and DEGs identification for the construction of causal network.

Analysis of key “triggers” and “effectors” 52 in the set of 215 DEGs (“ERG+ dependent genes”) according to MetaCore™ database of protein interactions (Thomson Reuters) suggests the following scenario of progression of normal to a ERG overexpressing cancer cells in the early stages of prostate cancer development (Figure 3). In this context, AR (androgen receptor) drives the unscheduled expression of the TMPRSS2-ERG in CaP cells, resulting in the production of truncated ERG (32 amino terminal amino acids) oncoprotein. The ERG oncoprotein as a transcription factor may then increase the expression of its targets including the activation of C-MYC oncogene1 and TWIST1 (twist basic helix-loop-helix transcription factor 1) that is associated with increased proliferative potential (mean fold change = 3.53; p-value = 6.8e-05). These two transcription factors (TFs) act as coactivators for several TFs such as AR, E2F1 (E2F transcription factor 1), Tcf (Lef) and changing the expression levels of ERG dependent genes.

Figure 3

Possible scenario of preinvasive to invasive stages of prostate tumorigenesis.

Network reconstruction method allows us to create hypotheses about further cancer development events (Figure 3 and 4). For example, TWIST1 by cooperation with histone deacetylase HDAC5 (histone deacetylase 5) activates HIF1A (hypoxia inducible factor 1, alpha subunit) 53 transcription factor and induces overexpression of GHR (growth hormone receptor) and Tissue factor. Indeed, our data indicates (see above) that TWIST1 is significantly overexpressed in WMD, ERG+ samples. Increased amount of GHR and Tissue factor may then initiate downstream signal transduction by MAPK and PI3K cascades and cause activation of downstream transcriptional factors. This pathway is shown in bold purple line on Figure 4.

Figure 4

Detailed self-regulation scheme of cancer development in ERG+ well differentiated tumor samples. Thin green links indicate activation interactions. Thin red links indicate inhibition interactions. Bold purple links - key pathway from Fig.3. Bold green links - DAN dependent AR activation. Bold beige links - MAPK cascade activation pathways. - down regulated ERG+ dependent genes. - upregulated ERG+ dependent genes. - triggers of the most important pathways.

An additional activation of MAPK pathway can possibly be induced due to an increased intercellular concentration of FGF2 (fibroblast growth factor 2) and Angiotensin II. Higher level of FGF2 and Angiotensin II may result from decreased expression of metallopeptidase Neprilysin in tumor cells, observed in this study (Figure 4). In normal process, Neprilysin breaks down signaling peptides FGF2 and Angiotensin II and thus inhibits hyperactivation of MAPK pathway, and, specifically, angiogenesis 54 In addition it has been shown that Angiotensin II serves as an activator for EGF (epidermal growth factor), which is also capable of activating the MAPK pathway 55 (Figure 4, bold beige lines).

Another route of activation that can lead to enhanced AR activity starts with a downregulated oncogene suppressor, DAN (neuroblastoma, suppression of tumorigenicity 1). In normal cellular context this protein is excreted into intercellular space and inhibits signaling peptides of BMP (Bone morphogenic proteins) family. BMPs initially reported to be regulators of bone biology have also been shown to influence other tissues. BMP receptors can phosphorylate transcription factors of SMAD family 56 Activated transcription factors SMAD, in turn, form complexes with CBP (CREB-binding protein) 57 which then can bind with AR and can enhance AR activity 57, 58 (Figure 4, bold green lines).

Overall, Figures 3 and 4 show a model of a cancer self-regulatory circuit that was uncovered by analyzing DEGs in ERG overexpressing WMD prostate tumor cells. Interestingly, we can discern three major pathways (shown in bold green, purple and beige lines on Figure 4) that contribute to the tumor development in the studied samples. In addition to the top-level process of cancer development shown in Figure 4 our newtwork analysis identified specific processes that some of ERG+ dependent genes participate in, thus contributing to the tumorigenesis. First, several ERG+ dependent genes are related to angiogenesis. The key factor in this process is the transcription factor HIF1 (Hypoxia Induced Factor 1). HIF1 may be activated in ERG+ tumors by several complementary pathways: by active MAPK cascade; by covalent modification brought about by histone deacetylases HDAC1 (which in our dataset has increased expresion) or HDAC5 (mediated by TWIST1); and by direct interaction with Angiotensin II, which is turned on by AR or IL-8. These observations may also reflect normal functions ERG in vascular endothelial cells 59-61 Second, activation of a sub-set of phospholipases that regulate arachidonic acid synthesis can be expected in the well-to-moderately defferentiated ERG overexpressing tumors. Phospholipases PLA2G4C and PLA2G7 are found among upregulated ERG+ dependent genes, and in addition, their activity can be increased by MAPK cascade. An observed decreased SFN (14-3-3 sigma) expression (which in normal processes inhibits cell cycle kinase CDK) can lead to accelerated cell cycle and uncontrolled proliferation 62 Third, translation processes can be accelerated by increased expression levels of some ribosomal proteins, translation initiation factor eIF2 and translation elongation factor eEF2, all of which are present in ERG, as well as, in the ERG-downstream C-MYC activated genes.

According to the analyzed data two independent pathways can be envisioned that support ERG+ cancer cell survival - protein folding increase and secretion of surplus products of translation in intercellular space. In ERG+ tumors an increased expression of TAOK3 protein (serine/threonine protein kinase TAO3) was noted, which activates IRE3 (endoplasmic reticulum to nucleus signaling 3) kinase responsible for launching cell survival pathways during unfolded protein response. IRE3 activates transcription factor XBP1 (X-box binding protein 1) which regulates a number of proteins involved in protein folding in endoplasmic reticulum and elimination of protein surplus 63 Further, the increased expression of myosin-6 can cause the secretion of large amounts of newly synthesized proteins, peptides, glycoproteins in intercellular region since myosin-6 forms structures that transport newly formed vacuoles from Golgi apparatus to cell surface through actin filaments 64, 65 Rab-11A protein from GTPase superfamily regulates vacuole transport to cell surface 66 and is overexpressed in most ERG+ samples. Both of the described processes and activated AKT pathway can inhibit apoptosis.

In summary, the analysis of the data on differentially expressed genes in ERG positive, WMD prostate tumor samples with TMPRSS2-ERG fusion reveals the following signature. Our data suggests the activation of MAPK and PI3K cascades by synergistic actions of several genes (FGF2, Angiotensin II, Neprilysin, AR, E2F1, Tcf, and DAN via three major distinctive signaling pathways (Figure 4). The major finding of this study is that an “out of context” expression of ERG in prostate cancer cells may activate processes resembling angiogenesis that may reflect the normal endothelial function of ERG. Further findings include increased proliferation, suppressed apoptosis and increased protein synthesis likely contributing to the overall process through distinct signaling pathways. We thus envision that key kinases from the identified pathways may be considered as therapeutic targets for potential ERG+ tumors.

Discussion

Lack of molecular typing of prostate tumors presents significant challenges in selecting targeted treatment strategies. For patients with localized cancer and low initial PSA levels, expectant management might be a preferred strategy, however, these patients must be distinguished from patients with progressive cancer phenotypes who require aggressive treatment 67 TMPRSS2-ERG fusion, detected in 50 to 70% prostate cancers, has been a subject of extensive evaluations 68-73 In this study we studied WMD subset of prostate tumor samples that had elevated levels of ERG expression in a defined cohort of radical prostatectomy patients with no previous hormone therapy. This type of WMD tumors may present special treatment challenges since typically they are associated with low (≤7) Gleason scores and low PSA levels at initial diagnosis and as such may prompt expectant management of patients. Understanding specific molecular mechanisms of cancer development in these tumors and identification of possible biomarkers and/or treatment targets can help clinicians to stratify patients into different treatment categories.

In the present study, we have identified 215 differentially expressed genes by comparing gene expression profile between ERG positive and ERG negative tumors versus matched normal epithelial cell in the well differentiated group. This gene set was then used for the reconstruction of causal network relying on protein interactions manually curated from literature and algorithms from MetaCore™ (Thomson Reuters). Application of system biology methods to a large set of data allowed an insight to complex interactions and causal dependencies that exist between potential disease biomarkers, transcription factors, key enzymes and other functional proteins represented by the set of differentially expressed genes. 51 These causal relationships can then be easily visualized on a network (such as shown in Figure 4) in the form of distinct pathways, processes and signal transduction routes. Some of the key benefits of visualizing protein interactions on a network constructed from knowledge base of literature curated interactions using well-validated and widely accepted algorithms 51 include ability to assess dependencies, key regulators, auxillary metabolic enzymes, as well as being able to see fine details of individual pathways. For example, some important genes vital to a pathway might have not been detected in a gene expression experiment either due to an experimental design, or due to a lack of expression at the particular timepoint of the experiment, or simply due to statistical constraints placed on the experiment. Our network (Figure 4) shows that HIF1A is an important participant of AR/ERG activation pathway (bold purple lines), but HIF1A was not found among 215 differentially expressed genes in our study. Thus, HIF1A could have been missed if we were to only construct a gene signature from our data without a consideration of functional significance revealed by network reconstruction method.

The causal network we constructed by analyzing the gene expression profiles of well differentiated, ERG-positive tumor samples offers a working hypotheses of the ERG-positive cancer phenotype. It follows from the network analysis that three major pathways leading to cancer development and progression are at play (Figure 4). All three pathways are related to a well-documented and cancer-related MAPK/PI3K cascade, but our network offers finer details on molecular mechanism of the MAPK cascade activation. Androgen-receptor activated ERG→MYC, TWIST1→HIF1A pathway, FGF2→Angiotensin II pathway and AR-activating DAN→SMAD1 pathway work synergistically and interdependently to promote tumorigenesis. A scrutiny of the network along these pathways can also yield a list of potential treatment targets that may be considered. Three major cancer pathways are supplemented by “typical” cancer-associated processes of angiogenesis, cell proliferation, increased translation and apoptosis inhibition.

Table listing genes that are differential expressed in well differentiated samples of CaP with ERG-TMPRSS2 fusion in comparison to normal prostate samples of the same patients.

Click here for additional data file.

Table listing differential expressed genes counts for different comparisons of experimental groups. For valuable differential expressed gene lists enrichment analysis is shown for pathways and other gene sets.

Click here for additional data file.

Conclusions

In the present study of matched pairs of isolated normal epithelial cells and WMD tumor epithelial cells with known ERG gene expression status from 20 patients with localized prostate cancer, we have discovered novel ERG associated biochemical networks. Using causal network reconstruction methods, we have identified three major signaling pathways related to MAPK/PI3K cascade that may indeed contribute synergistically to the ERG dependent tumor development. Moreover, the key components of these pathways have potential as biomarkers and therapeutic target for ERG positive prostate tumors.