Diversity in Androgen Receptor Action Among Treatment-naïve Prostate Cancers Is Reflected in Treatment Response Predictions and Molecular Subtypes.

BACKGROUND: Metastatic prostate cancer (CaP) treatments are evolving rapidly but without evidence-based biomarkers to predict responses, and to maximize remissions and survival.
OBJECTIVE: To determine the activity of androgen receptor (AR), the target for default first-line systemic treatment, in localized treatment-naïve CaP and its association with clinical risk factors, molecular markers, CaP subtypes, and predictors of treatment response. DESIGN SETTING AND PARTICIPANTS: We examined 452 bona fide AR target genes in clinical-grade expression profiles from 6532 such CaPs collected between 2013 and 2017 by US physicians ordering the Decipher RP test. Results were validated in three independent smaller cohorts (n = 73, 90, and 127) and clinical CaP AR ChIP-Seq data. Association with CaP differentiation and progression was analyzed in independent datasets. OUTCOME MEASUREMENTS AND STATISTICAL ANALYSIS: Unsupervised clustering of CaPs based on AR target gene expression was aligned with clinical variables, differentiation scores, molecular subtypes, and predictors of response to hormonal therapy, radiotherapy, and chemotherapy. AR target gene sets were analyzed via Gene Set Enrichment Analysis for differentiation and treatment resistance, Ingenuity Pathway Analysis for associated biology, and Cistrome for genomic AR binding site (ARBS) composition. RESULTS AND LIMITATIONS: Expression of eight AR target gene subsignatures gave rise to five CaP clusters, which were preferentially associated with CaP molecular subtypes, differentiation, and predictors of treatment response rather than with clinical variables. Subsignatures differed in contribution to CaP progression, luminal/basal differentiation, CaP biology, and ARBS composition. Validation in prospective trials and optimized quantitation are needed for clinical implementation.
CONCLUSIONS: Measurement of AR activity patterns in treatment-naïve CaP may serve as a first branch of an evidence-based decision tree to optimize personalized treatment plans. PATIENT
SUMMARY: Treatment options for metastatic prostate cancer are increasing without information needed to choose the right treatment for the right patient. We found variation in the behavior of the target for the default first-line therapy before treatment, which may help optimize treatment plans.

Introduction

The treatment landscape for metastatic prostate cancer (CaP) is evolving rapidly. Androgen deprivation therapy (ADT), the default first-line treatment, alone or in combination with radiation therapy (RT), induces a remission, but time to prostate-specific antigen (PSA) or clinical progression varies from several months to years [1], [2], [3]. Treatments for castration-resistant CaP (CRPC) include second-generation ADT or chemo-, immuno-, or radiotherapeutics [1], which provide short-lived remissions of variable duration [4], [5], [6], [7]. Recent trials of ADT combined with chemotherapy or second-generation ADT as first-line systemic therapy showed a survival advantage and delayed treatment resistance compared with ADT alone [2], [3], [8], [9], adding treatment options and questioning ADT monotherapy as the standard of care. Despite this expansion in therapies, each of which benefits patient subsets, the field lacks predictive biomarkers to facilitate treatment choices, design treatment plans that maximize remission for each consecutive therapy, and prolong survival of individual patients [10], [11].

Evidence-based decision-making tools that start with the distinct biological features of each unique patient’s CaP are needed if we are to optimize personalized treatment plans. The aim of this study was to determine whether the activity of a target for a specific CaP treatment prior to administration of this treatment predicts therapeutic response. We examined the spectrum of activity of the androgen receptor (AR), the target for ADT that is impacted also by RT [12], [13], in treatment-naïve localized CaP. Analyzing 6532 cases, we found that AR action varies widely and associates preferentially with the predictions of response to ADT, RT, or chemotherapy, and a CaP’s basal-luminal differentiation status and/or molecular subtype rather than with clinical variables. Distinct patterns of AR activity affected CaP biology and were reflected in AR binding site (ARBS) composition. Our studies indicate that a primary CaP’s pattern of AR activity may predict individual patients’ response to first-line therapy for metastatic disease and that this information can be developed into a first branch of a decision tree to optimize treatment plans.

Patients and methods

Study population

Deidentified transcriptome data from 6532 CaP patients were obtained from clinical use of the Decipher RP test between 2013 and 2017 ordered by US physicians. Patients treated with radical prostatectomy who had one or more adverse pathology features (ie, positive margins, extraprostatic extension, or seminal vesicle invasion) were eligible for the Decipher test.

Gene expression profiling, bioinformatics, and statistical analyses

Evaluation of AR activity in primary CaPs has mostly involved quantifying expression of small AR target gene signatures and/or smaller patient cohorts without considering its association with CaP progression. To assess AR activity more comprehensively, we examined gene expression profiles from 6532 localized treatment-naive specimens captured in the CLIA-certified Decipher GRID microarray platform [14] for the expression patterns of 452 bona fide AR target genes. We validated these genes both to be androgen-responsive and to contain ARBS(s) to which AR binds in an androgen-induced manner [15]. Demographic, pathological, and clinical patient information is shown in Supplementary Table 1. Clinical variables used to predict the risk of CaP progression (Gleason grade group, tumor-node-metastasis stage, and presurgical PSA) [1] were available for most patients. A total of 5000 cases had available molecular markers and gene signatures that were derived from gene expression data and were relevant to CaP therapy response. Molecular markers that were assessed included ERG and ETS fusion status (ERG+, ETS+), SPINK1 expression status (SPINK+), or the absence of these three markers (triple negative cases, triple–) [14]. Scores for gene signatures that reflect basal or luminal cell lineages such as the breast cancer–derived but CaP-validated PAM50 classifier (basal, luminal A, and luminal B subtypes [16]), or a signature that characterizes benign basal and luminal prostate epithelial cells (basal-like and luminal-like subtypes [17]) were considered. Furthermore, we incorporated gene signatures that classify prostate adenocarcinoma versus small-cell carcinoma, reflect genomic instability, or predict response to ADT (ADT response signature [ARS aka ADT-RS]) [18], RT (postoperative radiation therapy outcomes score [PORTOS]) [19], docetaxel, or dasatinib. Behavior of the 452 AR target genes was compared with the AR activity score from a smaller signature of nine AR target genes [20]. The incidence of molecular features was consistent with literature: 40.9% of cases were ERG+, 10.1% were ETS+, and 2.5% showed small-cell carcinoma features (Table 1).

Table 1

Distribution of molecular features studied among CaP clustersa

Molecular feature	All clusters(%)	Cluster 1 (%)	Cluster 2 (%)	Cluster 3 (%)	Cluster 4 (%)	Cluster 5 (%)
Molecular subtype
ERG+	40.9	31.63	44.11	1.85	88.19	39.56
ETS+	10.12	13.6	7.93	4.16	9.92	13.63
SPINK+	10.96	10.43	16.04	22.43	0.12	10.98
Triple–	38.02	44.33	37.43	71.56	1.77	35.51
Basal/luminal like
Basal like	29.46	85.71	27.01	2.89	4.01	27.65
Luminal like	70.54	14.28	72.99	97.1	95.98	72.35
Pam50 classifier
Basal	42	74.94	70.14	14.45	10.86	33.88
Luminal A	27.1	22.78	9.09	29.83	33.29	39.88
Luminal B	30.9	2.26	20.76	55.72	55.84	26.24
Small cell classifier
Adenocarcinoma	97.5	94.67	99.2	98.27	98.7	96.65
Small cell carcinoma	2.5	5.33	0.8	1.73	1.3	3.35
ADT response (ARS)
Lower	19.06	29.36	13.72	19.65	12.51	20.64
Average	49.6	50.91	51.6	46.7	44.04	52.57
Higher	31.34	19.72	34.67	33.64	43.44	26.79
AR-A (AROS)
Lower	11.92	42.85	9.27	1.27	0.59	8.02
Average	87.88	57.14	90.73	98.38	99.4	91.98
Higher	0.2	0	0	0.35	0	0
RT response (PORTOS)
Average	81.46	63.38	84.58	89.83	94.33	77.02
Higher	18.54	36.62	15.42	10.17	5.67	22.98
Docetaxel sensitivity
Lower	2.08	5.22	3.21	0.46	0	1.4
Average	90.28	94.44	95.98	86.47	79.1	92.36
Higher	7.64	0.34	0.8	13.06	20.9	6.23
Dasatinib sensitivity
Lower	0.84	0	0	1.73	2.72	0.31
Average	95.18	89.12	91.26	98.15	97.28	99.38
Higher	3.98	10.88	8.73	0.12	0	0.31

ADT = androgen deprivation therapy; AR = androgen receptor; AR-A = AR activity; AROS = AR output signature; ARS = ADT response signature; CaP = prostate cancer; PORTOS = postoperative radiation therapy outcomes score; RT = radiation therapy.

Data for each molecular feature are expressed as the percentage of total CaP cases showing that feature’s score or activity.

CaP subtypes, molecular markers, and gene signatures were analyzed using Gene Set Enrichment Analysis (GSEA), Ingenuity Pathway Analysis (IPA), and Cistrome tools in GRID data and multiple independent clinical CaP gene expression and AR-ChIP-Seq datasets, as detailed in the Supplementary material.

Results

AR action partitions based on molecular markers and predicted treatment response

Of the 452 AR target genes, 450 mapped to the Decipher GRID data (Supplementary Table 2). Unsupervised clustering based on the expression patterns of these 450 genes yielded five similarly sized clusters (clusters 1–5; Fig. 1). CaP clustering showed modest associations with some clinical variables used for risk stratification, for example, slightly more GGG5 cases in cluster 5 (Table 2). For other clinical variables, no or moderate associations were seen, for example, slightly younger patients and more extraprostatic extension in cluster 4. More pronounced associations were found between clusters and scores for predicted drug sensitivities, cellular differentiation, and molecular CaP subtype. For each cluster, one or more molecular features showed a considerable shift from the average obtained from all CaPs (Table 1). Cluster 1 cases differed most from average with regard to cell differentiation, with a higher percentage of basal-like, PAM50 basal, and small-cell carcinoma CaPs (doubled to 5.33% from 2.5%) and the lowest PAM50 luminal B scores, and lowest predicted response to ADT and docetaxel, yet the highest response to RT. A small subset from clusters 1 and 2 also showed the highest sensitivity to dasatinib. Although more luminal like, cluster 2 cases still had a high percentage of PAM50 basal CaPs, albeit with a higher fraction of luminal B subtypes and the lowest percentage of small-cell carcinoma. Cluster 3 showed the highest percentage of triple cases (71.6%), with very few ERG+ (1.8%) and ETS+ (4.16%) CaPs, yet the lowest basal-like score. Cluster 4 contained a high percentage (88%) of ERG+ cases, very few SPINK+ cases, the highest percentage of luminal B subtypes, highest ADT responses, and lowest RT responses, but highest docetaxel sensitivity. Cluster 5 was closest to average, although it contained the highest percentage of luminal A CaPs and somewhat lower sensitivity to dasatinib. These CaPs represent early-stage disease; no significant genomic instability was observed. The AROS-1 score, which summarizes the expression of nine target genes [20], differed markedly between clusters, with the lowest score in cluster 1 and highest in cluster 4.

Fig. 1

AR action partitions based on molecular markers and predicted treatment response. Deidentified Decipher GRID transcriptome data obtained from clinical use of the Decipher RP test for 6532 individual CaP cases were subjected to hierarchical clustering for 452 AR target genes, of which 450 mapped to GRID gene expression data. Heatmap shows clustering of CaP cases based on AR target gene expression. Five major CaP clusters are marked by colored boxes on the dendrogram at the top of the heatmap: cluster 1 in red, cluster 2 in blue, cluster 3 in black, cluster 4 in yellow, and cluster 5 in orange. Clusters 1–5 are labeled also at the bottom of the heatmap, where the annotations a, b, and c denote subclusters within each major cluster. Eight AR target gene sets the expression of which contributes to clustering are marked at the right of the heatmap. Panels at the top of the heatmap show annotation tracks for the molecular subtypes and gene expression signatures, and presurgical clinical variables. Legends at the bottom of the heatmap indicate scores and activity levels of gene signatures, subtyping, and parameters that were assessed. Color key on the left top corner indicates row z scores. ADT = androgen deprivation therapy; AR = androgen receptor; AR-A = AR activity score derived from nine AR target genes; ARS = ADT response signature; CaP = prostate cancer; GG = grade group; GGG = Gleason grade group; Lum = luminal; PORTOS = postoperative radiation therapy outcomes score; PSA = presurgical serum prostate-specific antigen levels (in ng/mL); pStage = pathological stage; RP = radical prostatectomy; RT = radiation therapy.

Table 2

Distribution of clinical variables studied among CaP clusters

Clinical variable	All clusters (%)	Cluster 1 (%)	Cluster 2 (%)	Cluster 3 (%)	Cluster 4 (%)	Cluster 5 (%)
Gleason grade group
1	7.65	9.31	9.45	5.90	7.09	6.47
2	43.15	42.79	47.86	37.23	47.28	40.56
3	27.84	24.29	24.51	31.68	30.97	28.55
4	8.93	8.74	6.95	13.53	7.92	8.35
5	12.43	14.87	11.23	11.68	6.74	16.07
PSA (ng/ml)
<10	76.04	73.96	79.74	70.33	78.63	76.63
10–20	17.77	19.13	13.34	21.97	17.45	17.75
>20	6.18	6.90	6.92	7.71	3.92	5.62
p stage
pT2	44.86	43.81	50.00	50.54	37.97	41.87
pT3a	34.50	29.85	33.50	34.93	41.35	33.74
pT3b	17.77	22.94	13.20	12.12	18.42	21.54
pT4	2.86	3.40	3.30	2.40	2.26	2.85
Age
1st quartile	60.7	61.1	60.5	62.3	59.7	60.6
Median	65.8	65.7	65.9	66.6	65.4	65.7
3rd quartile	69.4	69.6	69.6	69.8	68.8	69.2
Extraprostatic extension
Yes	51.42	50.89	47.97	46.50	58.86	53.17
No	48.58	49.11	52.03	53.50	41.14	46.83
Seminal vesicle involvement
Yes	18.30	23.30	13.75	12.40	18.67	22.58
No	81.70	76.70	86.25	87.60	81.33	77.42
Lymph node involvement
Yes	4.80	5.48	4.12	4.13	4.08	5.84
No	95.20	94.52	95.88	95.87	95.92	94.16
Bladder neck invasion
Yes	11.80	12.84	14.41	8.93	10.00	12.24
No	88.20	87.16	85.59	91.07	90.00	87.76
Surgical margins
Yes	50.85	46.35	50.50	52.44	53.93	51.11
No	49.15	53.65	49.50	47.56	46.07	48.89

CaP = prostate cancer; GGG = Gleason grade group; PSA = prostate-specific antigen.

Except for age, data for each clinical variable are expressed as the percentage of total CaP cases showing that feature’s score or activity. GGG, PSA, and p stage data were derived from 4996, 2926, and 4715 cases, respectively. Information on age is expressed in years and derived for 4853 patients. Information on extraprostatic extension, seminal vesicle involvement, lymph node involvement, bladder neck invasion, and surgical margins was derived from 4809, 4787, 4329, 1144, and 4909 patients, respectively.

Patterns of AR action, and not global AR activity, underlie CaP clustering

CaP clustering was not caused by differences, up or down, in the expression of the entire 450-gene signature. Analysis of Pearson correlations underlying the hierarchical clustering revealed that dissimilarities in the expression levels of eight AR target gene subsignatures (range: 32–90 genes, color coded in Fig. 1; Supplementary Table 2) produced the clusters. For instance, the expression of the black-coded signature was high in cluster 1 (which has a higher percentage of basal/basal-like carcinoma), lower in clusters 2–4, and intermediate in cluster 5, whereas the purple-coded gene set showed the highest expression in ERG+ cluster 4, but lower expression in clusters 1 and 2, and intermediate expression in clusters 3 and 5. Expression for the green-coded subsignature was lowest in cluster 4 and more enriched in others, particularly cluster 3. The gray-coded gene set showed higher expression in parts of clusters 1 and 2, and much lower expression in parts of cluster 4 and 5. Other subsets showed neither extremely high nor low expression among clusters (Fig. 1).

Differential CaP behavior of the eight gene sets was confirmed via GSEA. Similar to the 452 AR target gene signature [15], two gene sets (green and blue) were depleted in metastatic CRPC compared with localized treatment-naïve CaP [21], whereas the black-, purple-, orange-, and coral-coded subsignatures showed no difference, and the pink and gray ones were enriched in CRPC, suggesting the latter to be more relevant to CRPC progression. Comparing data sets derived from luminal and basal cells, luminal enrichment was seen for blue-, green-, coral-, orange-, and purple-coded gene sets, but no differences were found for black, pink, and gray signatures, indicating more basal features combined with luminal features in the latter (Fig. 2A and Supplementary Fig. 1).

Fig. 2

AR target gene sets differentially involved in CaP progression and differentiation. (A) GSEA analyses of eight AR target gene sets between localized treatment-naïve CaP and metastatic CRPC (GSE32269; top panel) and beween benign basal epithelial cells compared with benign luminal epithelial prostate cells (GSE67070; bottom panel). The dashed red line indicates FDR < 0.25, which is considered significant. Red bars indicate significant positive enrichment, blue bars indicate significant negative enrichment, and gray bars indicate no significant enrichment. GSEA plots are shown in Supplementary Figure 1. (B) Eight AR target gene sets were subjected to GSEA analysis as above to compare their behavior in ERG+ versus ERG– CaPs. Top panel: gene expression study (GSE84042) in which ERG status was verified via whole-exome sequencing. Middle and bottom panels: gene expression study (GSE120741 and GSE21032) in which 30% of CaP cases with the highest ERG gene expression were considered ERG+ and 30% of CaP cases with the lowest ERG expression were considered ERG–. GSEA enrichment plots are shown in Supplementary Figure 2. (C) Luminal and basal gene signatures from Decipher GRID were subjected to GSEA analysis as in Figure 1B to compare their behavior in ERG + versus ERG– CaPs. GSEA enrichment plots are shown in Supplementary Figure 3. AR = androgen receptor; CaP = prostate cancer; CRPC = castration-resistant CaP; FDR = false discovery rate; GSEA = Gene Set Enrichment Analysis.

We then used GSEA to validate AR target gene set expression patterns in three independent treatment-naïve localized CaP cohorts that were analyzed using diverse profiling platforms [22], [23], [24]. The considerably smaller number of cases (n = 73, 90, and 127) [22], [23], [24] permitted more stringent analyses of key features such as ERG+ status, a readily evaluated marker present in a considerable portion of CaPs [25]. For the three datasets, the purple-coded gene set (high expression in ERG+ cluster 4) was enriched, whereas the green-coded set (depleted in cluster 4) was depleted in ERG+ cases. The second and third study also revealed gray gene set depletion in ERG+ cases, reminiscent of its lower expression in part of cluster 4 cases (Fig. 2B and Supplementary Fig. 2). GSEA for the basal and luminal signatures in ERG+ and ERG– cases verified enrichment of luminal phenotype in ERG+ cases in all three studies, again validating the GRID association between ERG positivity and luminal features (cluster 4; Fig. 2C and Supplementary Fig. 3). To further characterize CaP’s heterogeneity in AR action, we verified cluster-specific expression and activity of genes, such as SPOP, PTEN, TP53, CHD1, and MYC, whose somatic alterations are used to subtype CaP and impact AR action [25]. We noted differences in their expression or activity among clusters, such as lower CHD1 expression in cluster 2 or near absence of mutant p53 activity in cluster 1, potentially contributing to differential AR output in these cases [26], [27]. Similarly, mutual exclusivity of SPOP mutations in ERG+ cases [25] was linked to cluster 4–specific AR action (Supplementary Fig. 4 and Supplementary Table 3).

AR target gene expression patterns impact CaP biology and reflect ARBS composition

We asked whether variable expression of the eight AR target gene sets affects CaP biology. First, we conducted IPA on individual AR gene sets, which revealed preferential enrichment in upstream regulators, canonical pathways, and diseases and functions (Supplementary Tables 4–6 and Supplementary Fig. 5–7) for each set. For instance, the pink gene set was enriched in catecholamine biosynthesis, whereas the black and grey signatures overlapped with the pluripotency of human embryonic stem cells. We performed similar IPA on reconstituted AR target gene set expression patterns among clusters (Fig. 3A). To account for the above-mentioned subtle variations in intracluster behavior of some gene sets, each cluster was subdivided into two or three subclusters (Fig. 1). IPA returned different numbers of entries per (sub)cluster as enriched or depleted, which ranged from 34 to 96 for upstream regulators, 219 to 292 for canonical pathways, and 9 to 105 for diseases and functions (Supplementary Tables 7–10). Subclusters within a cluster behaved more similarly; more pronounced differences were seen among the five principal clusters. Overall, differential combinations for the expression levels of the eight AR target gene sets observed in the (sub)clustering sufficed to induce shifts in CaP biology (Fig. 3B–D and Supplementary Tables 11–13). For instance, specific AR target gene set expression levels in basal cluster 1 enriched for cell functions related to neuronal development and cancer drug resistance by drug efflux. In contrast, more luminal clusters such as cluster 4 preferentially returned other cell biology, for example, lipid synthesis. IPA also suggested similarities between basal cluster 1 and more luminal cluster 5. Both were, for example, specifically enriched in cytokine (eg, CCL2) and growth factor signaling (eg, BMPs). IPA prediction of decreased SRC action in cluster 5 recalled its low dasatinib sensitivity in those CaPs. IPA reported different activities for transcription factors (TFs) such as p53, which is linked to context-dependent AR action [15], among clusters: it was predicted to be down in cluster 4, which showed the lowest percentage of cases with response signature for wild-type p53 action. In addition, differential involvement of AR and related nuclear receptors that can partially substitute for AR function was seen between clusters. Canonical pathway analyses confirmed variable androgen signaling, estrogen receptor action, and steroid metabolism among (sub)clusters. The Analysis Match feature, which considers results for the three IPA analyses, confirmed that subclusters within the same cluster were most similar and clusters differed considerably, with clusters 1 and 5 appearing most similar (Fig. 3E). IPA results verified that variable expression patterns of AR target gene subsets alter CaP cell biology.

Fig. 3

AR target gene expression patterns impact CaP biology. (A) Overview of Ingenuity Pathway Analysis (IPA) strategy. Based on the gene expression patterns observed in the Decipher GRID heatmap shown in Figure 1, for each CaP subcluster, the expression level of each AR target gene set was assigned a value ranging from –4 to +4. The combined expression patterns for the eight AR target gene sets (resulting in the 450 gene signature) per subcluster were subjected to IPA. Each column in Figure 3A represents a subcluster. (B) IPA upstream regulator analysis of gene expression pattern for each column/subcluster from Figure 3A. Each column represents a subcluster, and each row represents one upstream regulator returned as significantly inhibited (blue) or activated (red) in IPA analysis. The legend below Figure 3B specifies classes of upstream regulators returned. Supplementary Table 11 provides a more detailed view of the heatmap shown in Figure 3B, with each row annotated. Supplementary Table 8 provides specifics on IPA results used to generate the heatmap in Figure 3B. (C) IPA canonical pathway analysis performed as detailed for Figure 3B. Only differentially enriched entries are shown. Each column represents a subcluster, and each row represents one canonical pathway returned as significantly enriched in IPA analysis. Supplementary Table 12 provides a more detailed view of Figure 3C heatmap, with each row annotated. Supplementary Table 9 provides details of IPA results used to generate the heatmap. (D) Results for IPA diseases and functions analysis, which was performed as described for Figure 3B. Each column represents a subcluster, and each row represents one disease and function entry returned as significantly inhibited (blue) or activated (red) in IPA analysis. Only differentially enriched entries are shown. Supplementary Table 13 provides a more detailed view of the Figure 3D heatmap, with each row annotated. Supplementary Table 10 provides details on IPA results used to generate the heatmap. (E) Heatmap representing results for IPA analysis match studies. Legend on the right specifies the overall z scores. AR = androgen receptor; CaP = prostate cancer.

To verify that ARBS composition contributes to intercluster AR target gene expression heterogeneity, we performed Cistrome analysis of genomic regions encompassing ARBSs in the eight target gene sets. Five to 93 consensus transcription factor binding sites (TFBSs) were identified per gene set, with no correlation with gene set size and little overlap (Fig. 4A and 4B, and Supplementary Table 14). Gene set–specific TFBSs are bound by functionally diverse TFs that are relevant to the CaP clusters in which the TFBSs are over-represented. For instance, the black-coded gene set ARBSs showed preferential enrichment for TFBSs bound by factors that regulate embryogenesis and development, such as SOX8 and ATF4, consistent with more basal and stem cell qualities of cluster 1. However, ARBSs in the purple-coded gene set, highly expressed in ERG+ cluster 4, were enriched for TFBSs for ERG and related factors such as FLI1. These results indicate that differential composition of ARBSs underlies variable AR target gene expression among clusters and revealed CaP roles for TFs not yet studied extensively, such as Nkx1-1 and MYOG, which we confirmed to be expressed in GRID data (Supplementary Fig. 8).

Fig. 4

AR target gene expression patterns reflect ARBS composition. (A) Overview of Cistrome analysis strategy. For each AR gene set shown in Figure 1, AR binding peaks were extended by 1 kb 5′ and 1 kb 3′, and subjected to Cistrome analysis. (B) Heatmap summarizing clustering of over-represented transcription factor binding sites in ARBSs in the AR target gene sets. Rows represent AR target gene sets; columns represent binding sites returned as significantly enriched in Cistrome analysis (blue). Supplementary Table 14 provides a more detailed view of the heatmap shown in Figure 4B, with each column annotated. (C) Heatmap showing distribution of AR-ChIP-Seq peaks that overlap with ARBSs in the eight AR target gene sets in 88 clinical CaPs. Red implies that ARBS peak is present; blue indicates that ARBS peak is absent. Rows represent ARBSs; columns represent CaPs with AR-ChIP-Seq data available. (D) Heatmap for gene expression for the eight gene sets in CaP samples for which AR-ChIP data are available in Figure 4C. Rows represent AR target gene expression; columns represent CaPs with RNA-Seq data available. AR = androgen receptor; ARBS = AR binding site; CaP = prostate cancer; TFBS = transcription factor binding site.

Finally, we verified whether heterogeneous AR target gene expression is reflected in ARBS patterns using AR-ChIP-Seq data from 88 treatment-naïve CaPs [23]. While a majority of 452 AR target genes (>50%) overlapped with the AR ChIP-Seq peaks in a majority (>50%) of CaPs (Supplementary Fig. 9A and 9B), distribution of AR ChIP-Seq peaks from the 88 CaPs was heterogeneous across the eight gene sets (Fig. 4C). Furthermore, RNA-Seq data from the same CaPs on which AR-ChIP-Seq was done (Fig. 4D) showed heterogeneity in AR target gene expression over the eight gene sets. Together, these data validate the heterogeneity of AR cistrome in clinical CaP.

Discussion

Novel CaP drugs along with innovative trial designs and genomic profiling are rapidly transforming the metastatic CaP treatment landscape. However, we lack unbiased methods to decide on the right treatment for the right patient at the right time so that each consecutive therapy induces maximal remission and prolongs overall survival with optimal quality of life.

We reasoned that assessing treatment target activity before treatment could guide such decisions. Our first-in-field comprehensive analysis of the activity of AR, the target for the default first-line systemic treatment for non–organ-confined CaP [1], in thousands of treatment-naïve localized CaPs from four independent cohorts supports this possibility. Previously unrecognized heterogeneity in AR transcriptional output was reflected in the expression patterns of eight distinct AR target gene sets, but did not correlate with overall down- or upregulation of the overarching AR target gene signature. Mechanistically, this argues against simple correlation of AR target gene expression with AR level, but is consistent with molecular diversity in AR transactivation mediated by interaction with transcriptional regulators and ARBS sequence composition [15], which was supported by our IPA and Cistrome analyses. Plasticity in AR transcriptional action is well recognized during CaP progression, being caused by AR signaling axis adaptations to overcome selective pressure by ADT. Our findings here are the first to uncover the spectrum of baseline heterogeneity in AR action in CaP before systemic treatment, which is not fully captured by smaller AR target signatures, by simple quantification of AR target gene expression levels into an AR activity score, or by assessing a single AR target gene such as PSA or PSMA (Table 1 and Supplementary Fig. 10). Up- or downregulation of specific AR gene sets was mostly associated with treatment response predictions, but less so with clinical risk stratification. Treatment response findings likely reflect that AR differentially controls CaP biological processes. An AR target gene subset regulates DNA damage repair [12], [13], providing a rationale for differential association with RT response. Similarly, preferential silencing of specific AR target genes occurs after docetaxel treatment [28]. Dasatinib inhibits SRC, which induces context-dependent AR cistrome composition [29], and was tested for CaP treatment [30]; our data suggest that specific AR activity patterns may identify patients who will benefit most. That over/underexpression of select gene sets, rather than whole AR action, relates more strongly with the response to ADT was somewhat unexpected, although reminiscent of the association of molecular markers studied here with fractions of AR action and CaP treatment responses: ERG affects specific subsets of AR target genes [31] and CaP progression after ADT [32], as reflected in ERG + cluster 4, which shows the highest response to ADT. Luminal B CaPs are more likely to respond to ADT [16], consistent with congruent luminal B and ADT response scores (compare clusters 1 and 4).

Our findings suggest that information on intra-CaP action of CaP treatment targets can be exploited to rationally design evidenced-based treatment plans that are tailored to individual patients’ CaP biology and maximize response. One obvious limitation is the unknown accuracy of treatment response prediction signatures; with the exception of the PORTOS, these still require validation in prospective clinical trials. Another important consideration is the comprehensiveness of the target readout, which should determine full activity of all current or future targets. AR regulates not only mRNA, but also lncRNA and circular RNA levels, and its cistrome is impacted by somatic alterations [26], [27], all of which may not be captured by standard RNA-based readouts. Potential influences of CaP heterogeneity, clonality, or multifocality cannot be neglected, nor the real possibility of a shift in target activity during disease progression (Supplementary Fig. 11). The latter is an important consideration for building second and subsequent branches, and leaves, once prior treatments fail. In addition, procurement of metastatic CaP tissues and cells may be more challenging than obtaining a diagnostic biopsy or radical prostatectomy sample.

Conclusions

Our findings provide important novel insights into AR action in treatment-naive CaP, may explain variability in ADT and RT responses, and may provide a first branch of an evidence-based decision tree to optimize personalized treatment plans earlier in CaP progression. Next steps involve validating and refining AR target gene signatures to be used, alone or in combination with somatic alterations, in correlative studies of relevant clinical trials [2], [8], [9] and archival tissues from exceptional (non)responders in studies on preclinical monotherapeutics or combination therapeutics using models that allow modulation of AR activity and molecular markers.

Hannelore V. Heemers, Elai Davicioni, Jean-Noel Billaud and Song Liu had full access to all the data in the study and take responsibility for the integrity of the data and the accuracy of the data analysis.

Study concept and design: Ben-Salem, Heemers.

Acquisition of data: Ben-Salem, Hu, Wang, Y. Liu, Alshalalfa, Zhao, Billaud, Davicioni, Heemers.

Analysis and interpretation of data: Ben-Salem, Hu, Alshalalfa, Billaud, Davicioni, S. Liu, Heemers.

Drafting of the manuscript: Billaud, S. Liu, Heemers.

Critical revision of the manuscript for important intellectual content: Ben-Salem, Hu, Alshalalfa, Zhao, Wang, Venkadakrishnan, Senapati, Kumari, D. Liu, Sboner, Barbieri, Feng, Billaud, Davicioni, S. Liu, Heemers.

Statistical analysis: Hu, Wang, Y. Liu, Alshalalfa, Zhao, Billaud, S. Liu.

Obtaining funding: None.

Administrative, technical, or material support: Ben-Salem, Hu, Wang, Y. Liu, Alshalalfa, Zhao, Venkadakrishnan, Senapati, Kumari, D. Liu, Sboner, Barbieri, Feng.

Supervision: Davicioni, S. Liu, Heemers.

Other: None.

Hannelore V. Heemers certifies that all conflicts of interest, including specific financial interests and relationships and affiliations relevant to the subject matter or materials discussed in the manuscript (eg, employment/affiliation, grants or funding, consultancies, honoraria, stock ownership or options, expert testimony, royalties, or patents filed, received, or pending), are the following: Yang Liu, Xin Zhao, and Elai Davicioni are employees at Decipher Biosciences. Mohammed Alshalalfa is a former employee at Decipher Biosciences. Christopher E. Barbieri is a coinventor of a patent issued to Weill Medical College of Cornell University on SPOP mutations in prostate cancer. Felix Feng has consulted for Astellas, Bayer, Blue Earth Diagnostics, Celgene, Clovis, Genentech, Jansen, Myovant, Roivant, and Sanofi; is a cofounder of PFS Genomics, a molecular diagnostics company focused on breast cancer; is listed on patents for the use of the PAM50 and PORTOS signatures in prostate cancer; and receives grant funding from Zenith Epigenetics. Jean-Noel Billaud is employed by Qiagen, Inc.

These studies were supported in part by / grant CA166440 (to Hannelore V. Heemers), NIH grant U24CA232979 (to Song Liu), the NIH WCM SPORE in Prostate Cancer P50CA211024-01 (to Christopher E. Barbieri), and Prostate Cancer Foundation Young Investigator Award (to Deli Liu).

The authors thank Dr. Cassandra Talerico of the Department of Cancer Biology at the Cleveland Clinic for review of the manuscript.