Metabolic heterogeneity signature of primary treatment-naïve prostate cancer

To avoid over- or under-treatment of primary prostate tumours, there is a critical need for molecular signatures to discriminate indolent from aggressive, lethal disease. Reprogrammed energy metabolism is an important hallmark of cancer, and abnormal metabolic characteristics of cancers have been implicated as potential diagnostic/prognostic signatures. While genomic and transcriptomic heterogeneity of prostate cancer is well documented and associated with tumour progression, less is known about metabolic heterogeneity of the disease. Using a panel of high fidelity patient-derived xenograft (PDX) models derived from hormone-naïve prostate cancer, we demonstrated heterogeneity of expression of genes involved in cellular energetics and macromolecular biosynthesis. Such heterogeneity was also observed in clinical, treatment-naïve prostate cancers by analyzing the transcriptome sequencing data. Importantly, a metabolic gene signature of increased one-carbon metabolism or decreased proline degradation was identified to be associated with significantly decreased biochemical disease-free patient survival. These results suggest that metabolic heterogeneity of hormone-naïve prostate cancer is of biological and clinical importance and motivate further studies to determine the heterogeneity in metabolic flux in the disease that may lead to identification of new signatures for tumour/patient stratification and the development of new strategies and targets for therapy of prostate cancer.


INTRODUCTION
Prostate cancer is a clinically heterogeneous disease covering a wide variety of phenotypes, ranging from indolent disease to aggressive lethal disease. At present, there is a lack of methods for reliably identifying such phenotypes, rendering use of first-line androgen ablation therapy a challenge with a risk of over-or under-treatment [1]. As such, there is an urgent need for signatures to stratify indolent from aggressive prostate cancer which could be developed through identification of differences in energy metabolisms of prostate cancer subtypes [2].
Reprogramming energy metabolism is an important hallmark of cancer [3,4] and abnormal metabolic characteristics of cancers have been implicated as potential diagnostic/prognostic biomarkers [5][6][7][8][9]. Thus in contrast to normal resting cells whose energy production is mainly based on mitochondrial oxidative

Research Paper
Oncotarget 25929 www.impactjournals.com/oncotarget phosphorylation, the majority of cancer cells rely on aerobic glycolysis coupled to increased glucose uptake and lactic acid production/secretion (the Warburg effect) [3,7]. The augmented glucose consumption can be observed using ( 18 F)-fluoro-2-deoxyglucose positron emission tomography (FDG-PET) [10,11]. However, prostate cancers are characterized by a lack of increased glycolysis rendering FDG-PET less effective in imaging this malignancy [12], (although a small group of tumours do utilize aerobic glycolysis [9]). Instead, prostate cancers have abnormal metabolisms for citrate and choline [13][14][15] and are characterized by increased levels of cholinecontaining compounds (ChoCCs). These compounds have been linked to increased cell proliferation and survival [16]. Tracing metabolites with ( 11 C)-and ( 18 F)labeled choline derivative-based PET and PET/CT scans is increasingly used in the clinic. Together, these glucose and choline-based imaging studies demonstrated the metabolic heterogeneity of clinical prostate cancer [17,18]. However, it is not clear whether prostate cancer subtypes are metabolically different, or whether metabolic heterogeneity (intra-tumoural and inter-tumoural) can affect tumour response to treatment. In view of this, a better understanding of metabolic heterogeneity of prostate cancer at a molecular level will provide mechanistic insights into disease progression and may lead to identification of novel diagnostic/prognostic signatures for stratification of patients' prostate cancers.
To focus on the metabolic heterogeneity of prostate cancer cells and avoid complication from the stromal compartment, we used our unique panel of transplantable, patient-derived xenograft (PDX) models, derived from hormone-naïve prostate cancer tissues (www. livingtumorlab.com) [19,20]. Using these hormonenaïve prostate cancer PDX models and clinical samples, we, for the first time, demonstrated the heterogeneity of expression of genes involved in cellular energetics and macromolecule biosynthesis and identified an inter-tumour metabolic heterogeneity signature in clinical samples that is associated with patients' prognosis.

A panel of PDX models of hormone-naïve prostate cancers demonstrates heterogeneity of expression of genes involved in energy metabolism and recapitulates metabolic heterogeneity of clinical samples
A panel of 11 well-characterized LTL (Living Tumour Laboratory) prostate cancer PDX-models that are derived from hormone-naïve prostate cancer tissues, presents a unique opportunity for focusing on metabolic alterations in prostate cancer cells by analyzing human genes (while stromal genes are of mouse origin) in vivo [19]. To investigate whether these hormone-naïve PDX models present with molecular metabolic heterogeneity, transcriptome microarray data were generated and cluster analysis was performed using mean pathway scores (derived from per gene z-scores normalized to cancer mean and standard deviation, Supplementary Table 2). As shown in Figure 2, the metabolic pathway scores demonstrated variation among these models. For example, four tumours had increased pathway scores for glycolysis, while seven other tumours had decreased scores for this pathway. Heterogeneity was evident for all pathways.
To investigate whether such heterogeneity of metabolic gene expression observed in the PDX models represents the clinical situation, we analyzed transcriptome sequencing data of treatment-naïve, primary prostate cancer specimens of the VPC cohort (n = 14). The VPC prostate cancer per-gene z-scores were normalized to cancer mean and standard deviation, to generate mean pathway scores. The cluster analysis combining PDX models and VPC clinical samples showed that (in Figure  3), all the PDX models were clustered into different arms as observed for the clinical samples, suggesting that this set of PDX models recapitulates the heterogeneous metabolic nature of the clinical samples. In addition, the matched PDX/parental tissue pairs LTL-331/VPCT24 and LTL418/VPCT23, clustered closely together, as seen in dendrogram arms A and B respectively, illustrating the fidelity of the PDX models (arrow brackets). Interestingly, the LTL-313 xenograft lines, which had been derived from needle biopsies at five different foci of one patient's primary prostate cancer (asterisk), show a diverse heterogeneity signature in each focus. The observed differences between the five cores from which these models originate suggest the existence of diverging events that are continually associated with the development of intra-tumoural metabolic heterogeneity of prostate cancer. In summary, we documented molecular evidence of cancer metabolic heterogeneity in PDX models and clinical hormone-naïve prostate cancer samples. Importantly, the metabolic heterogeneity discovered in PDX models was validated in clinical prostate cancer samples.

Change in expression of genes involved in metabolism is consistent with the metabolomic changes in prostate cancer
To determine whether the gene expression involved in metabolic pathways is consistent with the metabolic properties of prostate cancer cells, we first analyzed the transcriptomic differences, between benign and treatmentnaïve prostate tumours, of genes involved in metabolic pathways and we then explored the correlation with clinical metabolomics. A comparison of the MSKCC dataset [21] consisting of 28 benign prostate tissues and Oncotarget 25930 www.impactjournals.com/oncotarget 112 prostate adenocarcinomas from treatment-naïve patients, showed significantly increased expression of genes involved in choline/phosphocholine generation (p < 0.0002) and in proline synthesis (p < 0.0029) ( Table  2) in tumours compared to benign tissue. These results reflect the findings of several recent clinical metabolomics studies reporting increased choline metabolites and increased proline [22][23][24][25] in prostate cancer compared to benign tissue.
Interestingly, expression of genes involved in pyruvate to acetyl-CoA formation was significantly decreased (p < 0.0092). As expected, there were no significant differences in the expression of glycolysis pathway genes (i.e. upstream of the PDH complex) between benign tissue and cancer (Table 2), which is consistent with previous reports and reflects the clinical finding that FDG-PET is less effective in imaging prostate malignancy [15] than other cancers.

Gene expression of a metabolic signature demonstrates heterogeneity in clinical samples and clinical outcome of prostate cancer patients
To further validate our finding of metabolic heterogeneity of treatment-naïve prostate tumours and explore its clinical significance, we analyzed the correlation of our heterogeneous metabolic signature  Z-scores were normalized to all tumours (mean/standard deviation) for each gene within a pathway. Average z-scores were then used to generate pathway scores. Red, upregulated; blue, downregulated; white, no change from tumour average. A, B, and b1 represent various arms of the dendogram.
Oncotarget 25932 www.impactjournals.com/oncotarget with clinical outcome in the MSKCC dataset [21]. As shown in Figure 4, heterogeneity of metabolic pathways between tumours is evident. Dendrogram arms I and II indicate two broadly different groups of tumours. Group II shows further heterogeneity as indicated by dendrogram arms (i and ii). The MSKCC validation cohort showed significant heterogeneity (Wilcoxon Rank test, p < 0.05) of metabolic pathways involved in cholesterogenesis, choline metabolism, gluconeogenesis, glutaminolysis, glycolysis, ketogenesis, lipogenesis, oxidative phosphorylation involving complexes I, IV and V, the TCA cycle and onecarbon metabolism (Supplementary Table 4). Pathways with less heterogeneity also occurred, including fatty acid activation and ketolysis. Interestingly, proline synthesis and proline degradation pathways showed opposite directions of expression. The proline synthesis pathway was consistently upregulated and proline degradation pathway was downregulated compared to benign tissue controls, which together supports metabolomics results [17] of increased proline content in aggressive prostate cancer. This is also consistent with results in the previous section (Table 2) where we showed that the proline synthesis pathway score was significantly increased in cancer compared to benign tissue (p < 0.0029).
To determine the clinical relevance of the observed metabolic heterogeneity we carried out survival curve analyses (5 year-overall survival as well as biochemical recurrence-free survival) for the MSKCC cohort (Supplementary Table 5). Patients with elevated expression levels of genes involved in one- The prostate cancer cohort of the Vancouver Prostate Centre (VPC) consisted of 14 untreated primary tumours and 3 matched benign prostate tissues [44]. Abbreviations: T, tumour; B, Benign; Dx, diagnosis; PSA, prostate specific antigen; nrn, never reached nadir; N, no; Y, yes; LTFU, lost to follow-up. There were no differences in stromal content across the cohort based on comparisons of desmin and vimentin expression (fibroblasts) [44]. www.impactjournals.com/oncotarget carbon metabolism, compared to controls ( > 2 Standard Deviations, Supplementary Table 6), showed significantly poorer biochemical-free survival compared to patients with lower expression of these genes (p < 0.01) ( Figure  5A-5B). Additionally, patients with low expression levels of genes involved in proline degradation, compared to controls (Supplementary Table 6), had significantly poorer biochemical-free survival compared to patients with higher expression levels ( Figure 5C-5D). Therefore, increased one-carbon metabolism or decreased proline degradation in these patients was associated with significantly decreased survival, suggesting these metabolic changes form clinically important metabolic signatures. There was a likely-correlation (Spearman) of one carbon metabolism and proline degradation pathway scores with Pathology stage (p = 0.06 and p = 0.08 respectively). This suggests that increased one carbon metabolism or decreased proline degradation (therefore increased proline levels) is observed in more aggressive tumours compared to less aggressive tumours. Such increased one carbon metabolism or decreased proline degradation signature was also observed in 50% or 75% (respectively) of PDX models with metastatic ability suggesting a predictive value of aggressiveness.

DISCUSSION
Whereas the genomic and transcriptomic heterogeneity of prostate cancer is well documented [26], less is known about metabolic heterogeneity of this malignancy and its relevance to malignant progression. Studies of the heterogeneity of prostate cancer, in particular hormone-naïve cancers, have been impeded by a lack of clinically relevant experimental models. Commonly used prostate cancer xenograft models based on cultured human prostate cancer cells do not accurately recapitulate the tumour heterogeneity, nor the cancer-stroma architecture of the original cancer specimens from which the cell lines were derived [19,20]. Importantly, the majority of prostate cancer cell lines and xenograft models used were derived from advanced metastatic castrationresistant cancer, rather than untreated hormone-naïve cancer. In this study, we demonstrated the heterogeneous expression of cancer cell genes involved in energy metabolisms in 11 patient-derived prostate cancer PDX models derived from hormone-naïve tumours. These models belong to a previously produced panel of prostate cancer PDX models [19,20] and to our knowledge currently form the largest collection of hormone-naïve prostate cancer PDX models in the field. Furthermore, the metabolic heterogeneity signature in these models suggests, for the first time, that a panel of PDX models can recapitulate heterogeneous metabolic signatures observed in clinical samples. The metabolic program that correlated with poor prognosis in the MSKCC dataset was associated with metastasis in PDX models. As such, these models provide a valuable platform for studying metabolic heterogeneity.
Using RNA sequencing strategies and microarrays to investigate heterogeneity of expression of genes involved in cellular energetics and macromolecular biosynthesis pathways, the present study documents intertumour metabolic expression heterogeneity of several pathways/sub-pathways in PDX models and in two clinical cohorts of, treatment-naïve prostate cancers as summarized in Figure 6. Consistent with these findings, heterogeneity of several metabolic pathways has been Figure 6: Summary of heterogeneity of expression of genes related to core metabolic pathways involved in cellular energetics and macromolecule biosynthesis in hormone-naive prostate cancer. The one-carbon metabolism pathway score was significantly increased (red) and the proline degradation pathway score was decreased (blue) compared to benign and both were associated with decreased biochemical-free survival in the MSKCC subset. It has been reported that the one-carbon metabolism and proline metabolism pathways interconnect with glycolysis and pentose phosphate pathway to provide ATP, and NADH and NADPH (mediators of electron transfer for energy, reductive biosynthesis and redox defense) [8,9,42] Oncotarget 25938 www.impactjournals.com/oncotarget observed in a meta-analysis of datasets of 22 different tumour types [27]. Importantly, we identified metabolic gene signatures of increased one-carbon metabolism or decreased proline degradation associated with significantly decreased biochemical disease-free survival of patients. One limitation of this study is the lack of metabolomics validation. Further integrated analysis of transcriptomic and metabolomics data and functional studies will provide more insights regarding the significance of metabolic heterogeneity in prostate cancer.
It is widely reported that cancer cells can direct glucose metabolism toward aerobic glycolysis [28]. However, consistent with our own results, emerging evidence suggests that prostate cancer cells can also direct glucose towards another ATP-generating pathway involving 3-phosphoglycerate conversion to serine through phosphoglycerate dehydrogenase (PHGDH), resulting in activation of one-carbon metabolism [29][30][31]. One-carbon metabolism facilitates NADPH generation [32] and the anabolic synthesis of amino acids, proteins, nucleotides and phospholipids, and has a role in the maintenance of redox balance and methylation reactions involved in post-translational modifications [33,34] and so may be a potential driver of oncogenesis. Indeed, methotrexate, an inhibitor of the one-carbon cycle, is used as an anticancer agent and preclinical studies are evaluating several other enzymes in this pathway [33]. Importantly, however, methotrexate has recently been associated with accelerated progression of a previously indolent prostate cancer [35]. An alternative therapeutic, metformin, has also been shown to indirectly inhibit one-carbon metabolism [36] and is being tested in preclinical trials for prostate cancer [37][38][39]. Several agents targeting metabolic enzymes exist for other metabolic syndromes that could be re-purposed for prostate cancer [37].
Proline metabolism promotes cancer cell proliferation and energy production [40]. The cycling of proline catabolism and proline synthesis acts as a redox shuttle between the mitochondria and cytosol to maintain redox homeostasis [41,42]. Furthermore, increased proline biosynthesis contributes to cancer cell growth by providing NAD+/NADP+ molecules to the pentose phosphate pathway and to glycolysis [40]. Here we show that a decreased pathway score for proline degradation is associated with significantly decreased biochemical disease-free survival of prostate cancer patients. Concomitantly, the expression of the proline synthesis pathway genes was consistently increased in prostate tumours implying a strong propensity of these tumours to activate the proline synthesis cycle. In support of our findings, LC/GC-MS metabolomic profiling identified increased proline levels with increased progression of prostate cancer [22][23][24][25].
Increased proline synthesis has also been reported as the major metabolic shift in metastatic breast cancer cells, in melanoma cell lines compared to melanocytes and in ovarian cancer cells (OVCAR3) compared to ovarian cancer stem cells [40]. Furthermore, a recent study of 133 metabolic genes identified five genes whose siRNAinduced knockdown resulted in strong growth-inhibitory effects on breast cancer cells. Included were PYCR1 involved in proline synthesis and the serine biosynthesis gene, PDGDH, involved in one-carbon metabolism [43].
In conclusion, we have, for the first time, demonstrated heterogeneity of expression of genes related to metabolic pathways in both hormone-naïve prostate cancer PDX models and clinical samples and provided insight into the metabolic heterogeneity of hormone-naïve prostate cancer associated with patient clinical outcome. These results not only suggest the potential of a metabolic heterogeneity signature as a biomarker for patient risk stratification, but also motivate further work to determine whether the observed differences in expression of genes involved in metabolic pathways translate into differences in metabolic flux. A better understanding of metabolic heterogeneity of hormone-naïve prostate cancer will lead to the identification of new signatures for patient stratification and the development of new strategies and targets for therapy of prostate cancer.

Patient-derived xenograft (PDX) modelsdiscovery cohort
Transplantable xenograft lines of patients' prostate cancer tissues were developed by grafting of fresh cancer tissue (obtained from needle biopsies), and serial transplantation, into the subrenal capsule (SRC) graft site of male NOD/SCID mice supplemented with testosterone as previously described [19]. Animal care and experiments were carried out in accordance with the guidelines of the Canadian Council on Animal Care.

Gene expression analysis of xenografts using RNA microarray analysis
RNA microarray analysis of PDX tissue was performed as previously described [19]. Data are available at GEO accession number GSE41193. Total RNA samples were prepared and processed using Agilent's One-Color Microarray-Based Gene Expression Analysis Low Input Quick Amp Labelling v6.0. An input of 100 ng of total RNA was used to generate cyanine-3-labeled cRNA. Samples were hybridized on Agilent SurePrint G3 Human GE 8X60K Microarray (Design ID, 028004). Arrays were scanned with the Agilent DNA Microarray Scanner at 3-μm scan resolution and data were processed with Agilent Feature Extraction 11.0.1.1. The processed signal was quantile normalized with Agilent GeneSpring 12.0. www.impactjournals.com/oncotarget

Validation cohorts
The prostate cancer cohort of the Vancouver Prostate Centre (VPC) consisted of 14 untreated primary tumours and 3 matched benign prostate tissues [44] as described in Table 1. Gleason scores ranged from 7-10 and samples consisted of ≥ 70% of malignant glands. The MSKCC cohort consisted of 112 untreated primary prostate tumours and 28 benign prostate tissues [21]. Data are available at GEO accession number GSE21032.

Gene expression analysis of tumor sections using RNA sequencing
VPC clinical tumour sections were processed and RNA sequencing was performed at the BCCA Michael Smith Genome Sciences Centre according to standard protocol as previously described [44]. Briefly, RNAseq data reads were first mapped onto the hg19 human reference genome and exon-exon junctions by spliceaware aligner STAR [45], using known gene model annotation from Ensembl Release 75. Reads with an unmapped mate or multi-mapped location were filtered out using Bam Tools [46]. Data are available at GEO accession number GSE55016.
Using aligned RNA-seq reads, gene expression profiles for each sample were calculated on the basis of gene annotation from Ensembl Release 75. Only reads that were unique to one gene and exactly corresponded with the structure of the gene were counted for the corresponding genes by using tool HTSeq [47]. In order to eliminate the variance of sequencing depth among samples, the raw read counts were normalized by R package DESeq [48].

Metabolic pathway score
Twenty-three pathways analyzed in this study include those of choline metabolism [13,14] as well as energy production and macromolecule biosynthesis, as previously curated for cancer cells [6,27,49,50] outlined in Figure 1 (and detailed in Supplementary Table 1). Assigning an 'upregulation' or 'downregulation' status to a metabolic pathway requires the contribution of all genes in the pathway. To this end, we selected well-established genes for each pathway, excluding upstream effectors.
In the PDX models, gene expression z-scores were derived from the Mean/SD of all 11 xenograft lines. The pathway score was then calculated from the average z-scores of genes within each pathway/sub-pathway. To compare the VPC cohort with the PDX models, z-scores for VPC gene expression were derived from the Mean/SD of the 14 VPC cancer samples. We obtained enough matched patient tissue for a paired comparison for 2 of the PDX models and original patient tissues (VPCT24:LTL331 and VPCT23:LTL418).
To validate metabolic heterogeneity in a larger cohort of primary untreated prostate cancer specimens (MSKCC), and to ascertain potential clinical predictive value, prostate cancer gene expression was normalized to benign gene expression (Mean/SD) to generate per gene z-scores for the genes expressed in 23 metabolic pathways/sub-pathways ( Figure 1).

Statistical analyses
The pathway scores calculated from the average z-scores of genes within a pathway provide information about the direction of the gene/pathway regulation. In the MSKCC cohort, Wilcoxon rank tests (p < 0.05) were applied to determine if pathway scores were significantly up-or down-regulated in tumours compared to benign tissue. Pathways containing fewer than 6 genes were not tested by the Wilcoxon rank test. The MSKCC cohort was also subjected to t-test (p < 0.05) using log2 gene expression (cancer vs benign) and survival curve analysis using pathway z-scores normalized to benign (Kaplan Meier plots) performed in Graphpad 6 Prism (significance determined using the Mantel-Cox test and the Gehan-Breslow test). Cox regression analysis and Spearman correlation were also undertaken to determine if metabolic pathway scores correlated with clinical parameters. Hierarchical clustering was used as previously described [51] using the R language (http://cran.r-project.org/).

CONFLICTS OF INTEREST
The authors declare that they have no conflicting interests.

GRANT SUPPORT
This study was supported in part by the Canadian

Author contributions
DL and SLE were involved in the study design, data analysis, data interpretation, literature search and generation of figures; SQ carried out data analysis; SQ, NN were involved in generation of figures; NN, SYCC and PWG performed literature searches; HX, AMH, SYCC carried out data collection; FM and RHB performed