Substantial interindividual and limited intraindividual genomic diversity among tumors from men with metastatic prostate cancer.

Tumor heterogeneity may reduce the efficacy of molecularly guided systemic therapy for cancers that have metastasized. To determine whether the genomic alterations in a single metastasis provide a reasonable assessment of the major oncogenic drivers of other dispersed metastases in an individual, we analyzed multiple tumors from men with disseminated prostate cancer through whole-exome sequencing, array comparative genomic hybridization (CGH) and RNA transcript profiling, and we compared the genomic diversity within and between individuals. In contrast to the substantial heterogeneity between men, there was limited diversity among metastases within an individual. The number of somatic mutations, the burden of genomic copy number alterations and aberrations in known oncogenic drivers were all highly concordant, as were metrics of androgen receptor (AR) activity and cell cycle activity. AR activity was inversely associated with cell proliferation, whereas the expression of Fanconi anemia (FA)-complex genes was correlated with elevated cell cycle progression, expression of the E2F transcription factor 1 (E2F1) and loss of retinoblastoma 1 (RB1). Men with somatic aberrations in FA-complex genes or in ATM serine/threonine kinase (ATM) exhibited significantly longer treatment-response durations to carboplatin than did men without defects in genes encoding DNA-repair proteins. Collectively, these data indicate that although exceptions exist, evaluating a single metastasis provides a reasonable assessment of the major oncogenic driver alterations that are present in disseminated tumors within an individual, and thus may be useful for selecting treatments on the basis of predicted molecular vulnerabilities.

Personalized approaches for cancer treatment rely on the premise that defined molecular alterations in tumors will confer susceptibility to specific therapeutics, and those malignancies without the alteration will not respond[1]. Consequently, identifying subsets of patients with and without a particular oncogenic driver is predicted to enhance efficacy and reduce futile treatment. This concept requires that heterogeneity exist between individuals with respect to the composition of therapy vulnerability characteristics and is supported by studies demonstrating extensive diversity in genomic aberrations and treatment responses between cancers with the same histological classification[2-8]. To achieve full success, precision oncology also requires limited molecular heterogeneity within an individual such that the targeted vulnerability is represented in most or all tumor cells.

While striking successes exploiting vital drivers of tumor cell growth have been achieved[9], deep molecular profiling has uncovered a substantial degree of intratumoral heterogeneity with divergent cancer clones existing within a primary tumor[10-12]. The extent of this diversity indicates the potential for misclassification using current biopsy-based strategies that may under-sample tumors, and provides clear pathways to rapid treatment resistance. To date, most studies assessing the extent of tumor heterogeneity have characterized primary tumors. However, localized carcinomas are generally managed with surgery or radiation, which are often curative. A major unmet need involves the systemic treatment of metastatic cancer, which usually involves multiple spatially distributed anatomic sites. There is limited information concerning the molecular diversity of the metastatic burden within an individual[13-15]. Such knowledge has considerable implications for the success of precision oncology.

Here, we analyzed the intra-individual molecular diversity of metastatic prostate cancer (mPC), a disease that accounts for more than 27,000 deaths annually in the US[16]. The mainstay of initial therapy for mPC, essentially unchanged for more than 50 years, relies on suppressing circulating levels of testosterone and consequently suppressing the androgen receptor (AR) signaling program in PC cells[17]. This therapy, conceptually one of the first precision medicine approaches for treating cancer, is highly effective with response rates exceeding 90% and durations of disease control lasting approximately 24 months. At the time of disease progression, a nearly universal event, metastatic castration resistant prostate cancers (mCRPC) exhibit a spectrum of recurrent molecular aberrations, a number of which are predicted to confer susceptibility to targeted therapeutics[18-21]. Our objective was to determine whether the molecular composition of a single metastatic focus at a particular point in the treatment continuum provides a reasonable assessment of the major oncogenic drivers of the other dispersed metastases in a patient. We used whole exome sequencing (WES), array CGH and transcript profiling to characterize the composition of metastatic tumors and evaluated the consistency and diversity of molecular aberrations within and between individuals.

RESULTS

Genomic landscapes of metastatic prostate cancers

We resected 176 primary or metastatic tumors from 63 men with mCPRC at the time of rapid autopsy and determined mutations by WES (n = 141 tumors from 56 men), copy number aberrations by array CGH (n = 149 tumors from 60 men), and transcript expression by microarray hybridization (n = 171 tumors from 63 men). A small subset of tumors could not be evaluated by all three methods due to limited tumor quantity. Histological assessments of each metastasis determined that the vast majority (n = 156) were adenocarcinomas while 20 tumors from 2 men were of small cell neuroendocrine histology. All men were treated with androgen deprivation therapy (ADT) and following disease progression, most also received additional treatments including at least one additional AR pathway-targeted agent and at least one systemic chemotherapy, most commonly docetaxel ().

As reported in previous studies of mCRPC, we found recurrent aberrations involving the AR, ERG, TP53, RB1, SPOP, CHD1, ZBTB16/PLZF and regions of chromosome 8q gain (including the MYC locus) and 8p loss (). We also identified recurrent mutations in FOXA2, which encodes a transcription factor expressed in tumors exhibiting neuroendocrine features ()[22,23]. Metastases from five men exhibited hypermutated genomes with complex structural aberrations in MSH2 and MSH6 mismatch repair genes[24]. All metastases and the primary tumors from these men were hypermutated, suggesting that mismatch repair deficiency occurred early in the genesis of these carcinomas.

Overall, the molecular landscapes of the lethal PCs exhibited a subset of aberrations observed recurrently across men, and also combinations of events unique to tumors originating within an individual. Excluding men with hyper-mutated tumors, the average number of protein coding mutations averaged 44 events per tumor, and the overall percentage of tumor genomes affected by structural changes, termed the DNA copy number alteration (CNA) burden[25], averaged 38%, though there were wide ranges in these metrics. Of interest most tumors from a given man exhibited similar numbers of mutations and CNA burden, and there were no differences in overall mutation abundance, or aberrations in any specific genes that associated with specific visceral or bone metastatic sites ().

AR Activity is Maintained in the Majority of CRPC Metastases

Amplification or mutation of the AR gene occurred in 63% of men (). Such aberrations were extremely rare in untreated primary PCs (Fisher's exact test; P < 0.0001) (). Using a gene expression signature of AR transcriptional output[26], we found that 88% of patients had tumors with robust AR activity, despite extensive prior treatment to suppress AR function (). To evaluate this more fully, we assessed AR expression and AR output across multiple metastases from individual men. AR activity scores () and AR expression () diverged across individuals, but were highly concordant when comparing metastasis within individuals. AR expression was positively associated with AR locus amplification (p ≤ 0.05) but did not differ between lymph node, visceral or bone metastatic sites (). AR genomic aberrations were generally concordant within individuals: all tumors were consistent in having no AR aberrations in 8 men, and for 24 men all tumors were found to have AR copy gain or mutation. In 9 men there was discordance with one or more tumors exhibiting an AR copy gain or mutation and one or more tumors with no aberration. We identified one man with evidence of convergent evolution where several metastases had AR mutation and other metastases had AR amplification, each potentially contributing to treatment resistance. In two men with small cell carcinomas confirmed by histology and immunohistochemical stains showing chromogranin expression, high neuroendocrine (NE) and low AR activity signatures were concordant across all metastases ().

Though AR expression was usually associated with AR activity, we identified outliers with high AR activity but with absent or limited expression of the AR itself (). Previous reports indicate that other nuclear hormone receptors such as the glucocorticoid receptor (GR;NR3C1), progesterone receptor (PGR), or estrogen receptor (ESR1) could activate canonical AR target genes[27-30]. In a small minority of tumors, GR;NR3C1 and ESR1 were expressed in the absence of AR (), and thus could account for the maintenance of at least a component of the AR transcriptome.

Limited genomic diversity within patients

The substantial concordance of AR- or NE-associated gene expression across metastases within an individual patient prompted a further assessment of intra-individual molecular diversity. Presumptive driver aberrations identified in an index tumor from each man were generally shared in all tumors. For example, patient 99-091 had AR copy gain, high AR activity, TP53 mutation, and lacked an ERG rearrangement (). The other 4 evaluable tumors from 99-091 shared these and other features. CNAs were also highly concordant (). At the level of mutations, a subset of non-synonymous single nucleotide variants (NSVs) were shared, though each tumor had private mutations, of which few were predicted to be deleterious (). To further assess the potential for these metastasis-private NSVs to influence pathogenesis, we examined WES data from 333 primary PCs from TCGA[31] and 150 metastatic CRPC tumors from the PCF/SU2C study[20]. Of 51 genes with NSVs exclusive to an individual metastasis from patient 99-091, only two, CSMD1 and NCOR2, were altered at frequencies exceeding 5% in either the TCGA or SU2C tumors. Similar results were observed across multiple tumors from other men (). Collectively, these data suggest that the vast majority of mutations unique to an individual tumor are not drivers of tumor behavior.

Using an integrated analysis that incorporated CNAs, mutations, and gene expression, we found that all tumors from an individual closely congregated with other tumors from the same individual, and were distinct from tumors in other men (). Unsupervised clustering based on CNA profiles grouped tumors from the same individual together with near 100% accuracy (). Of note, where the primary tumor removed at the time of rapid autopsy was analyzed, it grouped with the metastatic tumors from the same individual, rather than with primary tumors from different men.

We next evaluated the intra-individual consistency of a recurrent driver aberration, the TMPRSS2-ERG rearrangement, present in approximately 40% of primary PCs[32]. FISH using TMPRSS2 and ERG probes on 53 tumors from 13 men, found 100% concordance, all positive or all negative, between all metastasis within an individual (). Using array CGH to identify chromosomal loss between the TMRPSS2 and ERG genes, 32 of 34 men including the 13 men with FISH results had 100% concordance across all metastases (94%) with divergence in two men. For one, 00-029, FISH demonstrated all tumors were positive for a TMPRSS2-ERG rearrangement.

Clustering tumors based on expression profiles also grouped tumors by patient of origin, excepting bone metastasis (). This may reflect inherent tumor biology, microenvironment effects, the overall lower tumor cell content of bone metastasis, or technical factors associated with processing for transcript quantitation. While we and others have reported differences in the expression of specific genes between bone and visceral sites[33], we did not identify mutations or CNAs associated with metastasis location.

We next determined the overall molecular similarity and diversity of tumors using a statistical framework. We focused on the subset of 38 patients where multiple tumors per individual were successfully sequenced. Of the 123 tumors from this cohort, 5389 genes harbored a somatic NSV of which 2843 were from hypermutated tumors, leaving 2546 NSVs derived from 112 tumors from 34 patients. Collectively, there were 984 regions with defined CNAs in these tumors. We compared the CNAs in each tumor with every other tumor and determined that the median percent of corresponding aberrations that differed between tumors from the same individual was 6.7% compared to 22.1% from different individuals. The median inter- versus intra-individual tumor difference was 15.7% (95% CI 15.5-15.9%; P < 0.001; 2-sided 2-sample Wilcoxon rank sum test). Heatmaps with complete linkage dendrograms using copy number aberrations () or mutations () clustered tumors according to patient of origin.

Cell cycle progression inversely correlates with AR activity

The concordance of AR activity or NE activity in tumors derived from the same patient prompted further analyses to determine if a relationship existed between these programs and a key metric of treatment resistance: tumor cell proliferation. We evaluated a 31 gene signature of cell cycle progression (CCP) shown to associate with PC mortality when assessed in primary tumors[34]. A wide-range of CCP activity was evident across the CRPC tumors (). We confirmed that CCP activity was associated with cell proliferation determined by Ki67 immunohistochemistry on the same tumor samples (r = 0.48; P < 0.005)(). CCP activity was generally concordant between tumors within an individual compared to tumors from different individuals (). There were modest positive associations between CCP activity and mutation burden (r = 0.21; P = 0.023) ().

Contrary to expectations, AR activity was inversely associated with CCP activity (r = −0.33; P < 0.001) (), a finding confirmed in an independent dataset comprising mCRPC tumors ()[19]. AR activity was also inversely associated with E2F1, a key regulator of cell proliferation (P < 0.001) (). To determine the functional relevance of this correlation, we compared the proliferation of LNCaPWT cells expressing endogenous AR versus LNCaPAR cells engineered to overexpress AR. LNCaPAR proliferation was significantly greater than LNCaPWT cells at very low AR ligand concentrations of 10−11 M R1881 (P = 0.002) () that approximate castrate serum testosterone concentrations achieved with androgen deprivation therapy[35]. However, at higher ligand levels (10−9 M), approximating concentrations in a non-castrate male, the growth of LNCaPWT cells substantially exceeded LNCaPAR cells whose proliferation was suppressed by the higher androgen concentration (). Consistent with these results, after the introduction of AR into AR-null PC3 cells 10−9 M R1881 also significantly inhibited growth (P = 0.02) (). These findings support the concept that in addition to oncogenic roles, the AR may have activities that suppress aspects of tumor progression [36,37].

Cell cycle progression correlates with FA gene expression

We next sought to identify genes with expression levels positively associated with elevated cell proliferation. We computed genome-wide correlations of CCP score to each gene and found that elevated expression of Fanconi Anemia (FA) complex components including FANCA, FANCI, FANCD2, BRCA1 and BRCA2 associated with those tumors exhibiting high CCP activity (). Genes comprising the FA pathway were expressed at low levels in untreated primary PCs, which also had low CCP scores (). For several FA genes, including FANCA, BRCA2 and PALB2 we found significant relationships between copy number and transcript levels though most tumors did not have alterations in genes encoding FA components (). A subset of FA genes are positively regulated by E2F transcription factors, primarily E2F1[38,39]. We found a significant relationship between high CCP activity and loss of the retinoblastoma (RB1) tumor suppressor locus (P < 0.01) (), which also associated with E2F1 expression and remained significant after removing the highly proliferative NE metastases[40]. Tumors with high E2F1 expression also had elevated CCP scores (r = 0.8; P < 0.001) and elevated activity of a signature comprising 15 FA pathway genes (r = 0.78; P < 0.001) ().

To determine if FA components can influence cell proliferation, we suppressed expression of individual FA genes using siRNAs in LNCaP, 22Rv1, VCaP and PC3 cells. Suppression of each FA gene reduced the proliferation of each cell line (). Consistent with the known role of the FA complex in facilitating DNA repair, suppression of individual FA complex genes significantly increased phosphorylation of H2AX (γ–H2AX), a marker of DNA double strand breaks, in LNCaP cells following exposure to the genotoxic drug carboplatin (). Of the men comprising this rapid autopsy cohort, 20 received carboplatin chemotherapy during their clinical course. Those with a somatic aberration in a DNA repair pathway gene, defined as a homozygous deleterious event in any FA pathway component, or heterozygous inactivating event in the ATM gene, events recently shown to associate with responses to the PARP inhibitor olaparib[41], had significantly greater durations of response to carboplatin chemotherapy as determined by time on drug treatment (P = 0.02) ().

DISCUSSION

Comprehensive molecular assessments of PC have identified numerous recurrent genomic and epigenomic aberrations in genes implicated in driving cancer development, of which a substantial fraction have clinical implications for directing treatment[12,19,20,42-45]. These findings provide rationale for tumor sequencing to identify facets of tumor biology that prioritize and individualize therapeutics. However, deep sequencing of tumor genomes has also identified mutational heterogeneity, particularly within primary tumors. Branching evolution and the coexistence of multiple cancer lineages occur in a spectrum of carcinomas which have implications with respect to accurately sampling tumors, and for evolutionary mechanisms that may rapidly promote treatment resistance[11,12,46].

In the present study of mCRPC, we determined that constellations of genomic events were unique to a given patient from which a tumor was derived, demonstrating substantial inter-patient heterogeneity. In contrast, the molecular aberrations identified in multiple spatially-distinct metastases derived from a single individual and sampled at a single point in time were highly concordant, though a limited number of unique events specific to any particular tumor were evident of which very few represent known oncogenic drivers. We found that overall mutation burden, CNA burden, AR activity signature, and cell cycle output were generally concordant when comparing tumors within an individual. Our results provide further evidence for monoclonal seeding of metastatic sites, though our methods precluded assessments of their polyclonal composition or inter-mixing of metastases as reported by others[15]. It is interesting to consider that the ability of tumor clones from one site to seed other metastatic foci and intermix would result in a degree of homogenization, leading to the practical consequence that sampling one site has the potential to interrogate the extent of mutational heterogeneity present within an individual. Further studies are required to ascertain the influence of intra-tumoral heterogeneity with respect to treatment responses of tumors with variable frequencies of clones that possess druggable alterations.

Though the analysis of a single metastatic tumor site clearly does underestimate the total burden of molecular aberrations found in the totality of all metastases, the vast majority of drivers and actionable features are represented as common roots across all tumors or result from convergent evolution conferred by therapeutic pressure[15], and most molecular differences between metastases do not appear to influence tumor phenotypes. Though there are exceptions, these findings lead to the practical conclusion that clinical decision-making based on a biopsy from a single metastatic site is reasonable. This assumption may not be relevant when sampling a primary tumor where greater molecular diversity appears to occur[12], and must be tempered as only three men in our study were treated with the potent AR pathway antagonists abiraterone or enzalutamide that could induce diverse resistance mechanisms that include AR amplification, AR mutation, and expression of AR splice variants that may differ across metastases within an individual. In some situations, simply determining the activity status of a particular pathway will be sufficient for directing therapy, whereas in other instances the mechanism by which the pathway is activated, for example a specific AR ligand binding site mutation, will have treatment implications.

Our findings are supported by prior studies reporting substantially greater similarities in the genomic and epigenomic composition between metastases relative to the more expansive diversity found within regions of primary tumors[10-12,45-49]. In addition to the severe bottlenecks imposed by the metastatic process, the time scales involved in the evolutionary stages of cellular transformation to invasion and metastasis may underlie the disparity in molecular heterogeneity within primary tumors and between metastases. It has been estimated that decades of chronological time are required for cells with initiating events to acquire sufficient additional mutations that endow subclones with characteristics enabling invasion and metastatic spread[10]. In contrast, intervals between the acquisition of metastatic capability and the ability to grow in distant organs span only months to a few years[10]. Indeed, prostatic intraepithelial neoplasia and primary PCs are readily evident in men aged 30-40 years, decades before metastatic PC is commonly diagnosed[50]. These years of intra-organ tumor growth allow for substantial evolutionary diversity, clonal expansion, clonal mixing, and also for the genesis of completely independent primary tumors.

In the present study, an additional bottleneck limiting tumor diversity may be treatment pressures that likely eliminate sensitive tumor clones and allow the expansion of resistant clones. Extreme diversity across metastases would likely result in the rapid expansion of one or a few metastases and regression or stasis of others. Our finding that most metastatic sites within a given individual exhibited concordant proliferation rates suggests that there are limited differences across metastases with respect to treatment responses to androgen deprivation therapy or taxane chemotherapy. This finding may reflect the fact that most successful therapies in mPC are directed toward the AR which functions as a core lineage oncogene[51] or else toward fundamental cellular processes involved in cell division such as microtubule dynamics[6]. In this context, the observation that FA pathway activity is associated with the subset of tumors with high proliferation rates, likely reflecting important roles in replication stress[52], coupled with the high frequency of somatic events involving DNA repair proteins supports targeting this pathway in selected patients[41].

ONLINE METHODS

Patient Enrollment and Tissue Acquisition

All procedures involving human subjects were approved by the Institutional Review Board (IRB) of the University of Washington and of the Fred Hutchinson Cancer Research Center. Samples were obtained from patients who died of metastatic CRPC and who signed written informed consent for a rapid autopsy performed within 6 hours of death, under the aegis of the Prostate Cancer Donor Program at the University of Washington[53]. The inclusion criteria for the Prostate Cancer Donor Program included a diagnosis of prostate cancer, informed consent to perform an autopsy and informed consent to use biospecimens for biomedical research. Visceral metastases were identified at a gross level. Biopsies of all bone sites, except ribs (a dyke cutter was used for rib biopsies) were obtained using a drill with an attached trephine (11 mm diameter coring drill bit). All specimens were either flash frozen and embedded in Optimal Cutting Temperature compound (OCT; Tissue-Tek) or formalin fixed and embedded in paraffin respectively (bone metastases were decalcified in 10% formic acid before paraffin embedding). This study included 176 localized and metastatic tumors from 63 men with CRPC and a matched benign tissue sample from each individual to facilitate the identification of somatic mutations. For the molecular studies, the investigators were blinded as to sample identities and samples were randomly assigned to exome sequencing runs and CGH analyses. Samples were included for the molecular studies if histological analysis estimated >50% neoplastic cells, DNA quantity ≥ 3 μg and the presence of high molecular weight DNA determined by gel electrophoresis.

Transcript Profiling

A total of 176 tumors from 63 patients were collected for profiling by expression microarray. For five tumors, RNA was of insufficient RNA quality or quantity leaving 171 tumors from 63 patients that were successfully profiled. A Leica CM3050S cryostat was used to cut 8 micron sections, which were collected on PEN Membrane Frame Slides (Life Technologies) and immediately fixed in 95% ethanol. Sections were briefly stained (5–10 seconds) with Hematoxylin solu tion according to Mayer (Sigma-Aldrich) then dehydrated in 100% ethanol. Laser capture microdissection was performed with an Arcturus Veritas instrument, using both UV cutting and infrared capture lasers, with CapSure® Macro LCM Caps (Life Technologies) collecting 5,000-10,000 tumor cells per sample. Digital photographs were taken of tissue sections before, during, and after LCM and assessed by a pathologist to confirm the tumor content. Following LCM, captured cells were lysed in Arcturus RNA Extraction Buffer. RNA was isolated using the Arcturus PicoPure RNA Isolation Kit and DNAse treated using the Qiagen RNase-Free DNase Set. The bone metastases were sampled using a 1 mm diameter tissue punch of the frozen OCT-embedded tumors in a −20°C cryostat. The sample was obtained from the region of the block where there was tumor based upon a section of an adjacent decalcified FFPE block. RNA was isolated from the tissue cores using the RNeasy Plus Micro Kit (Qiagen Inc.). Subsequently, total RNA from all tumors was amplified for two rounds using the Ambion MessageAmp aRNA Kit. Probe labeling and hybridization to Agilent 44K whole human genome expression oligonucleotide microarray slides (Agilent Technologies, Inc., Santa Clara, CA) was performed following the manufacturer's suggested protocols and fluorescent array images were collected using the Agilent DNA microarray scanner G2565BA. Agilent Feature Extraction software was used to grid and extract the data. Data was loess normalized within arrays (normexp background correction with offset 50) and quantile normalized between arrays in R using the Limma Bioconductor package by batch. Control probes were removed, duplicate probes averaged, and spots flagged by Agilent Feature Extraction software as being foreground feature nonuniformity or population outliers were assigned a value of “NA”. An additional normalization step was applied to remove systematic batch effects due to date of RNA-amplification by application of the ComBat function within the sva Bioconductor package to the log2 Cy3 signal intensities. Expression microarray data are deposited in the Gene Expression Omnibus database under the accession number GSE74685

Copy Number Profiling

A total of 176 tumors from 63 patients were collected for profiling by microarray comparative genomic hybridization (CGH). For 27 tumors, DNA was of insufficient quality or quantity leaving 149 tumors from 60 patients that were successfully profiled. Genomic DNA was isolated from OCT sections of localized and metastatic tissues with ≥50% tumor content, 1mm punches from frozen bone metastases as well as patient-matched normal tissues using the Qiagen QiaAmp DNA Micro kit incorporating RNAse-treatment. DNA quality was assessed by agarose gel electrophoresis. One microgram of high molecular weight genomic DNA from each sample was labeled by random priming using the Agilent Genomic DNA Enzymatic Labeling Kit (Cy3-dUTP). A pool of reference normal DNA (Promega) was labeled with Cy5-dUTP. Cy3 and Cy5 probes were combined and hybridized to Agilent 2×400K SurePrint G3 CGH Microarrays and washed following the manufacturer's specifications. Fluorescent array images were collected using the Agilent DNA microarray scanner G2505C and Agilent Feature Extraction software. Control probes were removed, duplicate probes averaged, and log2 ratios defined as gProcessed-Signal/rProcessedSignal. Recurrent germline copy number variations (CNVs) were identified across benign tissues uninvolved with cancer obtained from 56 patients where DNA of sufficient quality and quantity was obtained for this study, as well as the Grasso, et. al copy number dataset of 28 benign samples, and a list of CNV regions defined by the Broad Institute and markers within these regions were removed from further analysis. The log2 ratios were corrected for the GC content wave effect using 1Mb windows for genome GC content. Outliers were smoothed and the normalized log2 ratios were segmented using the Circular Binary Segmentation (CBS) algorithm implemented with the DNAcopy Bioconductor package using default parameters and undo.splits = “sdundo” (undo.sd=1). The mode of the histogram of segment means was set to zero. To determine significantly recurrent regions of SCNA, the GISTIC 2.0 version 6.2 GenePattern module was applied to the segmented data using amplifications and deletions thresholds of 0.2, join.segment.size of 4 probes, a copy-ratio cap of 2.5, arm-level peel-off enabled, broad length cutoff of 0.7 chromosome arms, and confidence level of 95%. CGH array data are deposited in the Gene Expression Omnibus database under the accession number GSE74685.

Clustering

Unsupervised clustering of expression and CN profiles was performed in R on the 5000 most variable genes for each dataset using Pearson correlation as the distance metric and average-linkage hierarchical clustering. Clusters were visualized using the ape R package. Circular plots of expression and CN profiles by chromosome location were created with the OmicCircos R package version 1.8.0.

Statistical Analyses

Pearson's correlation coefficient was used to study the relationships between variables shown in scatterplots using the cor.test function in R. RB1 inactivation status was determined as homozygous if there was a two copy loss or 1 copy loss accompanied by mutation. Single copy loss or mutation only was categorized as heterozygous loss. The CCP activity in each of these groups was compared to the RB1 wild-type group using both pairwise t-tests and pairwise Wilcoxon rank sum tests, using the pairwise.t.test and pairwise.wilcox.test functions in R. The pairwise t-test and pairwise Wilcoxon rank sum test were also used to compare the number of nonsynonymous mutations and percent genome altered in different metastatic tissue sites, and all copy number vs. expression comparisons. Assessments of the proportions of bone or soft tissue metastases with homozygous PTEN cn loss or high gain of AR were performed by Fisher's Exact test using GraphPad Prism version 6.02, GraphPad Software, La Jolla California USA. Duration on carboplatin treatment was compared for percentages of patients harboring a DNA repair defect with those who did not by Kaplan-Meier plot with logrank test using GraphPad Prism version 6.02. Patients with a DNA repair defect were defined as those with any combination of the following criteria: homozygous or heterozygous loss in ATM by copy number loss and/or mutation, homozygous loss of one of the following 15 Fanconi anemia associated genes (FANCA, FANCB, FANCC, FANCD2, FANCE, FANCF, FANCG, FANCI, BRIP1, FANCL, FANCM, PALB2, SLX4, BRCA1) or germline mutation in BRCA2 and/or ATM.

Exome sequencing

A total of 176 tumors from 63 patients were collected for whole exome sequencing. For 35 tumors, DNA was of insufficient quality or quantity leaving 141 tumors from 56 patients that were successfully sequenced using the Nimblegen V2 or V3 platforms as previously described[18] with the following modification: in a subset of tumors, individually barcoded libraries were pooled in pairs prior to capture. Genomic DNA was isolated as described above. Sequencing was per formed using the Illumina Hiseq 2000 with either 50 bp or 100 bp paired end sequences. Reads were mapped to the human reference genome sequence (hg19) with bwa v0.7.1[54]. After remov al of PCR duplicate pairs we performed local realignment around indels using the Genome Analysis Toolkit (GATK)[47]. We subsequently called mutations using the Mutect software package with the following parameters: “--minimum_normal_allele_fraction 0.02 -- max_alt_alleles_in_normal_count 12 --intervals poscont.list --fraction_contamination 0.02”. To remove common polymorphisms and enrich for likely somatic mutations, we imposed a number of additional requirements, including requiring variants to be observed in at least 10% of reads at a position and removing variants present within a modified database of SNPs (dbSNP v137) that had first been stripped of all COSMIC variants. We investigated mutations for significance with MutSigCV using standard parameters, and inspected mutations in the top 50 significant genes manually using the Integrated Genomics Viewer (IGV) to remove sequence artifacts. The data are available for visualization and analysis in the cBioPortal for Cancer Genomics at: http://www.cbioportal.org/study.do?cancer_study_id=prad_fhcrc”.

Immunohistochemistry

One hundred and fifty-five metastatic PCa tissues (including 73 soft tissue metastases and 82 bone metastases) from a total of 50 autopsy patients (up to 4 sites per patient) were embedded in paraffin. A human tissue microarray (TMA) was constructed with duplicate 1 mm cores using the aforementioned tissues. Due to insufficient tissue, not all tumors subjected to molecular profiling were included on the TMA. Formalin-fixed paraffin-embedded TMA sections (5 μm) were deparaffinized and rehydrated with 3 changes of xylene and graded ethanol. Antigen retrieval was performed with 10 mM citrate buffer (pH 6.0) in a high pressure cooker for 30 min. Endogenous peroxidase and avidin/biotin were blocked respectively (Vector Laboratories Inc.). Sections were then blocked with 5% normal goat-horse-chicken serum at room temperature for 1h, and incubated with primary antibody at 4°C overnight. After washing 3 times with 1X PBS, slides were incubated with biotinylated secondary antibody (Vector Laboratories Inc.), followed by ABC reagent (Vector Laboratories Inc.) and stable DAB (Invitrogen Corp.). All sections were lightly counterstained with hematoxylin and mounted with Cytoseal XYL (Richard Allan Scientific). Mouse or rabbit IgG were used as appropriate for negative controls. Antibodies sources and dilutions: ER beta was from Biogenex cat # PU385-UP and used at 1:100; ER alpha was from Abcam, cat# ab173467, used at 1:500; GR/NR3C1 was from Santa Cruz, cat# sc-8992, used at 1:75; Androgen Receptor was from Biogenex, cat# F39.4.1, used at 1:60; Chromogranin A was from Dako, cat# DAK-A3, used at 1:100, and Ki-67 was from Dako, cat# M7240, used at 1:100.

Assessments of Inter- and Intra-Individual Tumor Molecular Heterogeneity

Mutation and copy number aberrations for 141 tumor samples from 56 patients that had successful sequencing and CGH analyses performed were initially examined separately. We identified 5389 mutations across tumor samples. Based on visualization of individual profiles, we flagged 5 patients as hypermutated; these patients together contributed 2843 (53%) unique mutations. To summarize heterogeneity across samples from non-hypermutated patients, we removed these patients and all mutations that were present only for these patients. This left 2546 mutations from 129 tumor samples from 51 patients. Next, we identified 984 regions with defined chromosome copy gain or loss across tumor samples from all 56 patients. We collapsed copy numbers with 1 or 2 gains into a single category (and similarly for 1 or 2 losses). The resulting mutation and copy number aberrations were then combined to create aggregate tumor samples. Finally, to compare heterogeneity across samples from the same individual vs between individuals, we restricted attention to individuals with at least 2 samples. This left a total of 3530 mutations or copy number aberrations from 133 tumors from 39 patients.

We quantified heterogeneity between two tumor samples using a “diversity index,” defined as the proportion of corresponding aberrations that differed between samples. For statistical inference, we randomly drew 10,000 pairs of tumor samples from the same individual and, separately, drew another 10,000 random pairs of tumor samples from different patients. We summarized distributions of the resulting diversity index using kernel density estimates and tested for equality of distributions using a 2-sided 2-sample Wilcoxon rank sum test.

As an alternative visualization of heterogeneity across tumor samples, we created heatmaps of the same diversity index calculated for all possible pairs of tumor samples. The heatmaps were sorted based on results of a complete linkage cluster analysis using the diversity index as the Euclidean distance between samples.

For the 3D plot of tumor molecular relationships (e.g. Figure 2C), we calculated a sample similarity score between the primary tumor of each patient and all other tumors using (a) single nucleotide alterations, (b) copy number alterations, and (c) gene expression. Similarity in single nucleotide alterations was measured by the number of impacted genes that are shared between two tumors. Similarity in copy number alterations was calculated as the inner (dot) product of the vectors of GISTIC-score per-gene for each tumor. Similarity in expression was measured by Pearson correlation across all genes.

Figure 2

Relationships between AR activity and the expression of AR and other nuclear hormone receptors

(a) The relationship between AR transcript expression plotted as mean centered log2 ratio to the AR activity score for each tumor. Blue points represent adenocarcinomas where AR activity levels correlate with AR level; black points are neuroendocrine tumors; red circles are tumors with high AR activity and relatively low AR expression. There was a positive overall correlation between AR transcript levels and AR activity score (r = 0.74, P < 0.001) using Pearson's correlation coefficient,

(b) Transcript levels of nuclear hormone receptors for the 15 tumors with high AR activity and low AR expression. NR3C1, glucocorticoid receptor; PGR, progesterone receptor; ESR1, estrogen receptor alpha, ESR2, estrogen receptor beta.

(c) Immunohistochemical assessment of AR (NR3C4), GR (NR3C1), PSA (KLK3) and chromogranin (CHGA) protein in mCRPC tumors from three men with different AR expression and AR activity relationships: 99-091 expresses AR and the AR regulated protein PSA while 03-139 lacks AR expression but expresses GR (NR3C1) and PSA (KLK3). Black asterisk indicates cells with positive staining. Black scale bar is 200 micrometers.

We generated one plot for each patient. Each coordinate axis in a plot represents one of the above three similarity measures. Plots were interactively rotated in 3D to obtain the clearest viewing angle. In the example plot of Figure 2C, the vertex at the top marks the location of an index tumor. The red plot points mark samples from the same patient as the index tumor. The blue plot points mark all other tumor samples.

Gene Manipulation and Cell Growth Assessments

LNCaP, PC3, VCaP and 22Rv1 cells were obtained from ATCC and used within six passages of growth after receipt or validated using STR (DNA Diagnostics) and confirmed to be mycoplasma negative. Cell lines were grown in RPMI 1640 supplemented with 10% fetal bovine serum (FBS). siRNAs were obtained from either Sigma (FANC F, L, M and BRIP1) or Qiagen (all others and the Allstars negative control). Target sequences are shown below. Cells were transfected with the siRNAs using Liopofectamine RNAiMAX (Invitrogen) following the manufacturer's protocol. Briefly, for gene expression analyses, 250,000 cells were plated in a 6 well dish with 50nM RNAi and 5μl of RNAiMAX in Optimem. After 8 hours the media was replaced with RPMI 1640 with 10%FBS. After 5 days RNA was extracted using the RNeasy kit (Qiagen). cDNA was made using SuperScript II (Invitrogen) and oligo dT primers. Gene expression was determined using qRT-PCR. Each assay comprised three biological replicates. Reactions were performed using Power SybrGreen on an Applied Biosystems (Foster City, CA) 7900HT Fast Real-Time PCR System. Cycling conditions included a 10 minute incubation at 95°C followed by 40 cycles of 95°C for 15 seconds with a one minute extension at 60°C. Disassociation curves were generated using the machines default conditions. Expression was normalized to RPL13A expression. Gene specific primer sequences are shown below.

Growth assays were performed using CellTiter 96 Aqueous One Solution (Promega). For each assay/condition, 10,000 cells per well were transfected as above using 10nM RNAi and 0.2μl of RNAiMAX in 100ul of Optimem. After 8 hours 100μl of RPMI 1640 with 20% FBS was added to each well. Each day, 10μl of CellTiter 96 Aqueous One Solution was added to designated wells, incubated at 37°C for 2 hours and the absorbance was measured at 490nM.

LNCaP-AR cells were the kind gift of Dr. Charles Sawyers (Memorial Sloan Kettering Cancer Center). PC3-AR were created by cloning the androgen receptor from a normal prostate cDNA library obtained from Clontech (Mountain View, CA) into the pBABE retroviral expression vector. After virus production and infection cells were selected using puromycin and clones selected and screened for AR expression. AR expression was determined by immunoblot using an antibody against AR (GTX62599, Genetex, Irvine, CA).

AR expressing cells were assayed for their response to varying concentrations of the synthetic androgen R1881. Five thousand cells were plated in charcoal stripped growth media for one day prior to exposure to androgen. Cell growth was determined as above using CellTiter 96 Aqueous One Solution. Percent cell growth was normalized to day 0 or to vector controls. Percent growth in each group was compared to the control using two-sample t-tests. Variances in each group were compared and not found to be significantly different by F test.

DNA Damage and γ-H2AX Assessment

Cells were transfected in 96 well plates with siRNAs as described above. Three days after transfection half of the plates were dosed with 50 μM of Carboplatin (Sigma) or vehicle. After 16 hours of exposure to Carboplatin or vehicle, the level of phosphorylated γ-H2AX was determined using the H2A.X (pSer139) Human In-Cell ELISA Kit (Abcam). The manufactures procedures were followed precisely. PCR primer sequences and siRNA sequences are in the Supplemental Tables 5 and 6.