Increased transcriptional and metabolic capacity for lipid metabolism in the peripheral zone of the prostate may underpin its increased susceptibility to cancer

The human prostate gland comprises three distinct anatomical glandular zones, namely the peripheral, central and transitional zones. Although prostate cancer can arise throughout the prostate, it is more frequent in the peripheral zone. In contrast, hyperplasia occurs most frequently in the transitional zone. In this paper, we test the hypothesis that peripheral and transitional zones have distinct metabolic adaptations that may underlie their different inherent predispositions to cancer and hyperplasia. In order to do this, we undertook RNA sequencing and high-throughput metabolic analyses of non-cancerous tissue from the peripheral and transitional zones of patients undergoing prostatectomy. Integrated analysis of RNAseq and metabolomic data revealed that transcription of genes involved in lipid biosynthesis is higher in the peripheral zone, which was mirrored by an increase in fatty acid metabolites, such as lysolipids. The peripheral zone also exhibited increased fatty acid catabolic activity and contained higher level of neurotransmitters. Such increased capacity for de novo lipogenesis and fatty acid oxidation, which is characteristic of prostate cancer, can potentially provide a permissive growth environment within the peripheral zone for cancer growth and also transmit a metabolic growth advantage to newly emerging clones themselves. This lipo-rich priming may explain the observed susceptibility of the peripheral zone to oncogenesis.


INTRODUCTION
Prostate cancer remains the second most commonly diagnosed cancer in males worldwide (15% of all male cancers diagnosed in 2012) [1]. Transcriptomic and metabolomic approaches are being used for the development of new biomarkers for both the early detection of prostate cancer, and to distinguish indolent from aggressive forms of this disease [2][3][4][5], but there have been relatively few attempts to integrate these approaches into a single analysis. The prostate gland has several important metabolic adaptations, notably an inhibition of tricarboxylic acid (TCA) cycle activity in epithelial secretory cells resulting in export of citrate into the seminal fluid [6][7][8], and the ability of prostate cells to accumulate high amounts of cholesterol, which has been

Research Paper
Oncotarget 84903 www.impactjournals.com/oncotarget attributed to enhanced uptake and synthesis of cholesterol [9]. These metabolic processes, amongst others, are perturbed within prostate tumours.
Among the most prominent changes observed during prostate cancer progression are those occurring in lipid-associated metabolic pathways. Several changes to key gene regulators of lipid metabolism, such as fatty acid synthase (FASN) and alpha-methylacyl-CoA racemase (α-AMACR) have been described in prostate cancer [10][11]. Recent molecular and metabolic profiling of prostate cancer also identifies lipid metabolism as a key pathway that undergoes metabolic reprogramming during prostate cancer development [12][13][14]. These changes include an upregulation of genes and metabolites responsible for de novo lipid biosynthesis and fatty acid β-oxidation, and are likely to facilitate tumour development by meeting energy demand for increased cell proliferation and also by supplying important metabolic intermediates for new cell growth.
The occurrence of cancer is not equally distributed across the three glandular zones of the prostate gland; the peripheral (PZ), central, and transitional (TZ) zones [15]. These zones are located within the same prostatic tissue but in different areas of the prostate and are not easily distinguishable by anatomical or histological means other than their location relative to the urethra. The majority (~70%) of diagnosed prostate cancer is found in the PZ [16,17].
In this study we comprehensively measure whole genome expression by next-generation sequencing and combine this with high-throughput metabolic profiling of the PZ and TZ from non-cancerous prostate tissue. Our aim was to identify the molecular and metabolic differences between the two zones that could underlie their differential susceptibility to prostate cancer. Using comprehensive analyses and novel omics integrative approaches we clearly demonstrate, for the first time, that the PZ is innately capable of increased lipid metabolism compared to the TZ, which may predispose the PZ to carcinogenesis and tumour growth. Our study significantly increases our understanding of prostate biology and highlights the unique metabolic nature of the prostate gland.

Patient demographics
Clinical characteristics of patients, including the demographic and tumour stage details, are shown in Table 1. Their average age was 62 ± 6.7 years, and the average PSA preoperatively was 8.7 ± 4.2 μg/L. All patients had an American Society of Anaesthesiologists (ASA) co-morbidity score of 2, which reflects the average general health status of men in this age group.

The two prostatic zones are transcriptionally distinct
We undertook RNA sequencing analysis of noncancerous prostate tissue derived from the peripheral and transitional zone of nine patients (seven from PZ and seven from TZ, Table 1) in order to, first, determine whether the two prostate zones differ significantly in their transcriptional profiles, and secondly to identify the genes responsible for any such differences. Unsupervised hierarchical clustering on the normalised expression levels revealed that the two zones are transcriptionally distinct ( Figure 1A). In addition, the transcriptional profiles of the respective zones from the prostate specimens were more similar to each other, than were the two zones from the same prostate gland, as evidenced by the distinct clustering of zones rather than clustering of individual patients. Further analysis identified a total of 2,764 transcripts that were differentially expressed between the two zones after multiple testing corrections (FDR q < 0.05). Of these 1,252 were higher in the peripheral zone compared to the transitional zone ( Figure 1B).
We then sought to identify whether there was an enrichment of a specific cell type within each zone. The expression of the stromal marker (vimentin, VIM), the basal markers (keratin (KRT)5 and KRT14), and the luminal cell marker (KRT8) were extracted from the RNAseq data ( Figure 2A). VIM expression was higher in the TZ samples by 1.8-fold (q < 0.001), whereas expression of KRT8 was lower by 2-fold (q <0.005). These reflect the zonal cell distribution, characterised by increased stromal composition in the transitional zone compared to the peripheral zone.
In order to independently validate the zonal origin of our samples we used an 11-gene signature previously shown to classify zone of origin of prostate cancer samples [18]. Principal component analysis (PCA) separated accurately the two zones using the 11-gene signature (Supplementary Figure 1), confirming our sampling, and also further supporting this diagnostic gene signature. Moreover, the comparative expression of the eleven genes between the two zones from our study correlated accurately with the direction of change previously reported [18] (Figure 2B).

Transcription of genes involved in fatty acid metabolism is higher in the peripheral zone
In order to identify whether the distinct gene signatures between the non-cancerous PZ and TZ are associated with particular biological processes, we performed functional analysis to identify pathways that may be overrepresented in each zone. We found that the genes that were relatively more highly expressed in the non-cancerous PZ (q-value < 0.05, 1252 genes) were Oncotarget 84904 www.impactjournals.com/oncotarget enriched in functions related to fatty acid metabolism and steroid biosynthesis ( Table 2). These pathways remained significant when over-representation analysis was undertaken on more stringent expression lists (q-value < 0.01, 833 genes) and on less stringent expression lists (q-value < 0.1, 1500 genes). These included the major metabolic lipid enzymes involved in de novo lipogenesis, FASN, acetyl-CoA carboxylase isoform alpha (ACACA), acyl-CoA synthetase long-chain (ACSL)1 and ACSL3, which were upregulated by 2.8-, 1.4-, 2.0-, and 1.9fold respectively ( Table 3). Several genes involved in cholesterol and steroid biosynthesis were also increased, including farnesyl-diphosphate farnesyltransferase 1 (FDFT1) and 24-dehydrocholesterol reductase (DHCR24), which were upregulated by 1.8-and 2.4-fold respectively ( Figure 3, Table 3).
In contrast, genes that were relatively more highly expressed in the non-cancerous TZ belonged to a variety of pathways predominantly involved in signal transduction and immune regulation (Supplementary Table 1) and were enriched in Gene Ontology functions of extracellular matrix organisation and immune response and signal transduction (data not shown).

Metabolites involved in de novo lipid biosynthesis are higher in the peripheral zone
PCA analysis of metabolic profiles from the PZ and TZ of 18 men (a total of 36 samples) did not clearly differentiate the two zones (Supplementary Figure 2A). Subsequent partial least square-discriminant analysis (PLS-DA) gave clearer separation but this was not statistically significant (data not shown). This is likely due to the influence of other factors, notably inter-individual differences.
We then used N-way ANOVA to identify the metabolites that were significantly different between the two prostate zones, having factored in the structure of the data (multiple measurements per patient). This identified 50 metabolites that were different between the zones (p < 0.05), of which 44 were higher in the PZ and 6 higher in the TZ (Table 4).
We performed enrichment analysis of the metabolites that were upregulated in each zone to give us an insight into the pathways that are involved. Similar to the transcriptional data, steroid hormone biosynthesis and fatty acid biosynthesis were among the    [18]. The comparative expression between the two zones mirrors that reported previously.
Oncotarget 84906 www.impactjournals.com/oncotarget most overrepresented pathways in the PZ, compared to phenylalanine biosynthesis and the pentose phosphate pathway, both of which were over-represented in the TZ ( Figure 4, Supplementary Table 2).

Integration of transcriptional and metabolomics data
Network-based integrative analysis of the combined RNAseq and metabolomics data further supported the strong differences in fatty acid and steroid metabolism between the two zones identified by the individual analyses. In particular, over-representation analysis of biological processes finds lipid and fatty acid metabolism pathways amongst the most significant between the PZ and TZ with a high proportion of differentially active molecules ( Figure 5, Supplementary Table 3). In addition, integrative analysis of the combined datasets performed better in identifying lipid metabolism as a key pathway differentiating the PZ and TZ, compared to RNAseq data alone (rank 16 compared to rank 50; Supplementary Tables  3 and 4). Genes and metabolites differentially expressed between the PZ and TZ are forming distinct clusters in the human molecular interactome network closely related with lipid metabolism (Supplementary Figure 3).

DISCUSSION
Although cancer is more likely to arise in the peripheral than the transitional zone of the prostate, there is currently little understanding of the origins of such a difference. As there is no clear significant anatomical difference, it is likely that such susceptibility is driven by molecular differences. Earlier reports using microarraybased genome-wide expression profiling between the two zones identified some genes but with limited functional interpretation [18][19][20]. In this report, we have employed whole-genome RNA sequencing and parallel highthroughput metabolite profiling of non-cancerous prostate tissue derived from the peripheral and transitional zones to comprehensively study the molecular and metabolic differences that underpin the increased susceptibility of the PZ to develop cancer. In addition, we have applied novel 'omics integration bioinformatics analysis to improve our predictions for the differential pathways that characterise each zone. With the use of genome sequencing and comprehensive metabolomics analysis of the two zones we have shown that the PZ has a unique transcriptional and metabolic signature that is characterised by increased capacity for both fatty acid and steroid biosynthesis. We have demonstrated that the PZ has an almost 3-fold increased expression of FASN, which encodes the enzyme that catalyzes the formation of long-chain fatty acids from acetyl-CoA. FASN gene expression is linked to the content of fatty acids in non-cancerous prostate tissue [21], and is consistent with the increased levels of phospholipids seen by metabolomics analysis. High FASN expression is associated with prostate cancer progression [22] and is linked to higher Gleason score both at the gene [23] and protein [24] level. Additionally, gene expression of ACACA, the rate-limiting enzyme in FA biosynthesis, was also higher in the PZ. Inhibition of ACACA by RNAi led to inhibition of LNCaP prostate cancer cell growth and subsequent de novo lipogenesis [25].
Increased de novo fatty acid and cholesterol synthesis is one of the features of prostate cancer cells [26,27] and acts as a source of lipids for membrane production but also of signalling molecules activating oncogenic pathways, such as the phosphoinositide 3-kinase (PI3K) pathways [28]. The bulk of cell membrane lipids are phospholipids (PLs), such as phosphatidylcholine (PC) and phosphatidylethanolamine (PE), and other lipids, such as sterols, sphingolipids, and lyso-PLs. Lysolipids in particular have been suggested as biomarkers for prostate cancer [29]. It is notable that our metabolomics analysis identified several such lipids to be higher in the PZ compared to TZ, which supports a higher de novo Oncotarget 84908 www.impactjournals.com/oncotarget lipogenesis activity in the PZ. More recently, increased lipogenesis in prostate cancer tumours was also found to be protective from endogenous and exogenous insults that lead to oxidative stress-induced cell death [30]. This could be another manner by which newly emerging cancerous clones in the PZ are evading elimination.
A significant proportion of the genes leading to cholesterol and steroid biosynthesis were higher in the PZ compared to the TZ. These included FDFT1, which is expressed at higher levels in human prostate cancer specimens and in aggressive cancers [31], and DHCR24, also known as seladin-1, an androgen regulated gene that is expressed significantly higher in prostate cancer and is positively related to T stage [32]. Circulating levels of cholesterol as well as enhanced cholesteryl ester accumulation in cellular lipid droplets are positively associated with prostate cancer and aggressiveness [33,34]. These results are consistent with reports that pharmaceutical interventions, e.g. statins that target cholesterol metabolism, are associated with reduced prostate cancer recurrence [35,36].
Increased FA biosynthetic activity in the PZ was coupled with increased FA catabolic capability. ACSL1 and ACSL3 were highly expressed in the PZ compared to the TZ by more than 2-fold. The protein encoded by the ACSL3 gene is an isozyme of the long-chain fatty-acid coenzyme A ligase family that converts free long-chain fatty acids into fatty acyl-CoA esters that then enter fatty acid β-oxidation for subsequent catabolism. ACSL3 is an androgen receptor (AR)-stimulated gene [37] and was recently identified in formalin-fixed prostate cancer tissue as a new 5′ translocation partner of ETS variant 1 (ETV1), an important gene in driving prostate cancer progression [38,39]. Interestingly, pristanoyl-CoA oxidase 3 (ACOX3), encoding an enzyme involved in peroxisomal β-oxidation, was also higher in the PZ. Consistent with the gene expression, higher levels of several acyl carnitines were found in the PZ. These medium and long-chain acyl carnitines contain FAs and are mediating their transport into mitochondria for β-oxidation. Acyl carnitines are sourced from diet as well as produced endogenously, and have been linked to prostate cancer and associated with aggressive forms of the disease [29,[40][41][42]. Increased mitochondrial and peroxisomal FA β-oxidation is a feature of prostate cancer progression and is likely a unique metabolic adaptation that provides ATP and acetyl-CoA in prostate cancer [13,43,44].    two zones (p < 0.05 for metabolites, q < 0.05 for genes). Node size corresponds to number of genes that belong to the GO term. Node color represents statistical significance: white -not significant, yet some genes present, in most cases significant before multiple testing correction; scale yellow to orange (from 5E-2 to 5E-7) -the darker the orange the more significant after multiple testing correction.
Oncotarget 84912 www.impactjournals.com/oncotarget Interestingly, we also found increased levels of three neurotransmitters in the PZ, namely N-acetylaspartyl-glutamate (NAAG), serotonin (5-HT), and acetylcholine (Ach). NAAG is a neurotransmitter that is found in abundance in the mammalian nervous system. It is catabolised to glutamate and N-acetyl-L-aspartate (NAA) by two groups of glutamate carboxypeptidases (GCP II and III). Both GCPs are membrane bound and have also been identified in the human prostate. GCPII in the prostate is called prostate specific membrane antigen (PSMA) [45] and although is expressed in both benign and cancerous prostate tissue, the degree of expression is higher in prostate cancer and increases with cancer severity [46]. Radiolabelled indium coupled with an anti-PSMA antibody has shown some promise in detecting prostate cancer recurrence after surgery (ProstaScint ® scan, Cytogen Corporation, US). Inhibitors of these GCPs have also shown anti-tumour properties in animal models of prostate cancer [47]. There are no studies comparing the immunohistochemical staining of PSMA in different zones of the prostate, however, the abundance of NAAG in non-cancerous PZ tissue could indicate the low activity of GCPII. McDunn and colleagues found the breakdown product NAA to be associated with organ-confined prostate cancer when compared to cancer that had extended beyond the capsule of the prostate, indicating more catabolism of NAAG with advanced disease [48].
Levels of 5-HT were also significantly higher in the PZ. 5-HT is a metabolite derived from tryptophan, and functions as a neurotransmitter but also has growth promoting properties in several cancer cell lines, including androgen resistant prostate cancer [49][50][51]. In hormone resistant prostate cancer, neuroendocrine (NE)-like cells emerge and are associated with poor prognosis. NE cells typically secrete 5-HT, amongst other neurotransmitters, that are thought to have an autocrine as well as a paracrine effect, promoting cell growth and migration [52]. Higher levels of 5-HT in the PZ therefore could favour the growth and migration of cancer. Finally, Ach was also found to be higher in the PZ, compared to the TZ. Ach is a neurotransmitter associated with the parasympathetic nervous system, but increased levels were recently linked to high-grade/advanced prostate cancer in both animal models and in human tissue [53]. Although increased neurotransmitter metabolites in the PZ could represent the abundance of nerve tissue in this part of the prostate compared to the TZ, their emerging role in prostate cancer development and progression could also indicate another metabolic feature of the PZ that could explain cancer predilection to this part of the gland. Further studies are required to confirm these findings.
One limitation of the study is that the non-cancerous samples were obtained from patients with confirmed prostate cancer, albeit the sampling site was histologically confirmed to be cancer-free. Such prostates could potentially have acquired an altered metabolic profile as part of the carcinogenesis process. If that was the case, we could expect such metabolic alterations to be universally acquired across the prostate, spanning both the peripheral and transitional zones. Our data show that there are major differences in the metabolic and transcriptional profiles of these two zones that have a direct relevance for prostate cancer initiation and growth.
One explanation for the observed transcriptional differences between the two zones could be the inherent differences in the cellular composition of the zones, namely the peripheral zone being epithelial cell-rich and the transitional zone being stromal cell-rich. Such compositional differences could be the origin of the molecular differences between the two zones but do not detract from the observation that increased transcriptional and metabolic capacity for lipid metabolism is a characteristic of the peripheral zone.
Recently, the genetic complexity of morphologically normal prostate tissue was uncovered using genomewide DNA sequencing and suggested that mutational processes, consistent with field effects, are underlying prostate carcinogenesis [54]. In our study, we extended this observation to include genome-wide RNA sequencing and high-throughput metabolic profiling of morphologically normal prostate tissue distant from cancer. We have focused on the apparent differences in the two major zones of the prostate, the peripheral and the transitional zones, and have identified the potential underlying mechanisms for the increased susceptibility of the PZ to prostate cancer. We propose that the PZ is constitutionally primed for metabolic reprogramming that is vital for prostate cancer progression and is characterised by increased de novo FA synthesis, FA catabolism and increased levels of neurotransmitters. Such 'priming' can benefit prostate cancer initiation and progression in two ways. First, by providing a permissive environment, full of vital cell growth intermediates, for novel cancer clones to arise and thrive, and secondly, by transmitting a metabolic advantage to the newly emerging clones themselves. In this way, the clones that arise in the PZ have an advantage over those arising in the TZ.

Patient selection and sampling
Eighteen patients undergoing radical prostatectomy for organ confined prostate cancer were consented to allow the use of prostate tissue for research via the Biorepository at the Norfolk and Norwich University Hospital (NNUH) under ethical approval granted from the Faculty of Medicine and Health Sciences Research Ethics Committee of the University of East Anglia (FMHS 20122014-37). Patients underwent endoscopic extraperitoneal radical prostatectomy. All procedures were conducted by a single surgeon. The prostate gland was removed from the Oncotarget 84913 www.impactjournals.com/oncotarget abdominal cavity immediately after resection and rapidly biopsied after extraction to avoid ischaemic artefacts using a standard core biopsy instrument. The extracted whole prostatectomy specimen was placed on a surgical table or equivalent and the subsequent steps were taken to collect the biopsy samples: (1) the apex and base of the gland were identified, the prostate was then cut transversely (axial section) half way along the gland, (2) using the midline and the urethra as guides, biopsies were taken from the peripheral and transitional zones avoiding obvious tumour sites where feasible, (3) after collection of biopsies the prostate was sent for standard diagnostic histopathological examination. One tissue core from each zone was placed in RNAlater and stored at −80°C for downstream RNA sequencing. Two tissue cores from each zone were also placed in 80% HPLC-grade methanol/water for 24 hours for downstream metabolite profiling as per manufacturer's guidelines (Metabolon Inc.). The methanol extracts were then frozen at -80 ºC and tissue cores were subsequently sent for histological analysis. As the cores destined for metabolite profiling and RNA sequencing were taken from directly adjacent sites we used the histological evaluation on the metabolite profiling cores to determine whether samples were cancerous or not. The cores were confirmed during histological assessment to consist of a mixture of epithelial and stromal cells.
Clinical characteristics of patients are shown in Table 1.

RNA sequencing and analysis
A total of 18 PZ and TZ tissues cores from nine patients were sent for RNA sequencing. Histology of the directly adjacent region later confirmed that 14 of these cores were of normal histology, and were used for analysis in the current study, and four of malignant histology. Prostate biopsy cores stored in RNAlater and weighing between 4 and 16 mg were homogenised with a QIAGEN TissueRuptor before total RNA was extracted with the QIAGEN RNeasy Mini kit. The resulting RNA was quality checked with an Agilent Bioanalyzer and samples with RIN values higher than 7 were further processed. Samples were subjected to ribodepletion with the Ribo-Zero Magnetic Gold rRNA Removal Kit (Illumina) before constructing Illumina barcoded TruSeq RNA libraries. Sequencing of 18 libraries, in pools of four, was performed over four and a half lanes of Illumina HiSeq 2500/2000 in High-Output mode using 100 bp paired-end reads. RNA-seq reads were first processed by removing Illumina adapters using Trim Galore! version 0.3.7 (Babraham Bioinformatics) and reads with Phred quality of basecalls higher than 20 and with a length of higher than 60 bp were carried forward. SortMeRNA version 1.9 [55] was used to filter any remaining ribosomal RNA from the adaptor and quality trimmed reads.
Reads were analysed using the TopHat-Cufflinks pipeline [56] aligned to the Ensembl GRCh37.75 reference genome. Only the 14 samples with non-cancerous histology were taken for further analysis. The Cuffnorm package was used to generate normalized expression levels for each sample. Genes with low expression (counts less than 100) were filtered out before performing unsupervised hierarchical clustering using Euclidean distance and an average linkage function to generate the hierarchical tree presented in Figure 1. We standardised expression values per gene across all samples and then used colour to visualise the variation in gene expression, where red denotes the highest and blue the lowest expression for that gene across all samples. The Cuffdiff package was also used for differential expression analysis of the PZ and TZ. Differentially expressed genes were identified as those that were different between the two zones at q<0.05 and the zone average FPKM values was higher than 1 in both zones. KEGG Pathways [57] that were overrepresented were identified using Database for Annotation, Visualization and Integrated Discovery v6.8 (DAVID; http://david.abcc.ncifcrf.gov/) [58]. Significant pathways were identified according to their over-representation (ORA) scores obtained from DAVID. ORA scores are defined as the Benjamini -Hochberg adjusted modified Fisher's exact test calculated from a 2x2 table test crossclassifying genes according to differential expression and pathway membership. RNA sequencing data have been deposited in ArrayExpress (E-MTAB-5021).

Metabolite profiling and analysis
A total of 72 samples from eighteen men undergoing radical prostatectomy were sent to Metabolon Inc. (USA) to undergo ultrahigh performance liquid chromatography -mass spectroscopy (LC-MS) and gas chromatography -mass spectroscopy (GC-MS) with a high resolution accurate mass (HRAM) platform as previously described [59,60]. These included two cores (one from each prostate lobe; left or right) taken from each zone (peripheral or transitional). Samples were sent in two batches. The first batch (n = 8 patients, 32 samples) returned 266 individual metabolites. The second batch (n = 10 patients, 40 samples) was processed after instrument upgrade (Metabolon ® , USA) using the exact techniques and standardisation procedures as with the first batch and yielded 413 metabolites. There were 160 common metabolites, amongst which 150 were endogenous, between the two datasets that were taken forward for analysis. As we had two cores for each zone per patient we only further processed those cores with confirmed benign histology. In the case of both cores being benign the mean was calculated for subsequent multivariate analyses. Principal component analysis (PCA), partial least square-discriminate analysis (PLS-DA), and pathway analysis was carried out with the open-source web-based MetaboAnalyst software v3.0 (http://www.metaboanalyst. ca) [61,62]. Data uploaded onto MetaboAnalyst were log 2 www.impactjournals.com/oncotarget transformed and autoscaled by dividing the metabolite values by the standard deviation. Over-representation (ORA) analysis and subsequent ORA scores were obtained from MetaboAnalyst by applying Fisher's exact test.
Individual metabolite differences between the two zones were examined with univariate methods by analysis of variance (N-way ANOVA) in Matlab (The Mathworks, Cambridge, UK). For each metabolite the average of the zone (TZ, PZ) was used as the response variable and the log 10 (to correct for departure from normality) of each measurement, the patient (18 factors), and the laterality of the core (left, right) as explanatory variables.

Integration of transcriptomics and metabolomics data
Metabolomics and transcriptomics analyses were integrated with the aid of a human molecular interaction network. We used PathwayCommons version 8 (www. pathwaycommons.org; accessed 7th of July 2016) to collect interactions between genes, proteins and small molecules that correspond to biochemical pathways, gene regulation pathways and protein-protein interactions in human cells. The gene expression and metabolite concentrations data were then used to discover the extent of the connectivity of nodes (molecules) involved in differentiating the two prostatic zones. In effect, we find one connected component comprising 723 nodes (genes and metabolites) including 362 differentially expressed genes directly associated with metabolites also found as significantly differentiated. All network analyses and visualisations were performed in Cytoscape 3.1 [63]. We tested overrepresentation of Gene Ontology terms with Bingo plugin [64] in Cytoscape by applying hypergeometric test (alpha = 0.05) to genes directly associated with metabolites (Supplementary Figure 3 and  Supplementary Table 2) and complete connected components of differentially expressed genes (Supplementary Table 3). Biological processes were selected as significant for q < 0.05 after applying correction for multiple testing with the Benjamini-Hochberg test.

ACKNOWLEDGMENTS AND FUNDING
This work was supported by a Prostate Cancer Foundation Challenge Award and an Institute Strategic Programme Grant in Food & Health by the UK Biotechnology and Biological Sciences Research Council (BB/J004545/1). Next-generation sequencing and library construction was delivered via the BBSRC National Capability in Genomics (BB/J010375/1) at the Earlham Institute (Norwich, UK) by members of the Platforms and Pipelines Group.