Evolution of Castration-Resistant Prostate Cancer in ctDNA during Sequential Androgen Receptor Pathway Inhibition

This article is featured in Highlights of This Issue, p. 4457

Metastatic prostate cancer is a lethal disease, but can be transiently controlled by inhibition of androgen receptor (AR) signaling (1). This inhibition is achieved clinically via suppression of androgen synthesis or direct antagonism of the AR ligand binding domain (LBD). To overcome these treatments, prostate cancers develop augmented AR signaling through mutations altering AR ligand binding affinity, AR amplification, or LBD truncation (2–4).

Several potent AR signaling inhibitors have clinical activity in AR-augmented prostate cancer (5), but cross-resistance is common (6, 7). Adaptive mechanisms can drive renewed resistance (8, 9), but most posttreatment tumors retain AR signaling, suggesting that the pathway remains an oncogenic driver (10). The extent to which the AR gene and the wider somatic genome continue to evolve and drive acquired resistance over sequential lines of AR-targeted therapy remains unknown, largely because of the difficulty of serially biopsying a treatment-standardized cohort of patients.

Recently, we reported primary outcomes from a phase II trial (NCT02125357) of sequential AR inhibitor treatment in metastatic prostate cancer that has already progressed on conventional androgen deprivation therapy (a clinical state described as metastatic castration-resistant prostate cancer, mCRPC; ref. 7). We analyzed pretreatment plasma cell-free DNA (cfDNA) samples collected from the 202 trial patients and demonstrated that circulating tumor DNA (ctDNA) fraction and key genomic alterations were independently prognostic for time to progression (TTP; ref. 11). Our prior work was limited to one time point, had short follow-up, and did not address acquired resistance. Here, we leverage this clinical trial population to present the first comprehensive analysis of longitudinal cfDNA collections (n = 458) from a large standardized mCRPC cohort. We show that most somatic driver alterations are consistently detected in ctDNA over sequential clinical progressions. However, temporal shifts in cancer cell populations represented in ctDNA were common, and the AR locus was highly plastic, evolving toward aggressive genotypes associated with increased AR activation and ligand independence. Collectively, this work provides insight into real-time evolution of the metastatic niche during treatment and proposes a genomic explanation for resistance to the most commonly used mCRPC therapies.

Blood samples were collected as part of an open-label randomized phase II crossover trial (NCT02125357) conducted in six centers in British Columbia, Canada (7). In the published trial, patients with newly diagnosed mCRPC were randomly assigned 1:1 to receive abiraterone 1,000 mg orally daily + prednisone 5 mg twice daily followed by enzalutamide 160 mg orally daily (arm A), or the opposite sequence of enzalutamide followed by abiraterone + prednisone (arm B). Patients were to crossover from the first treatment assigned (first-line treatment) to the second treatment (second-line treatment) at the first occurrence of either confirmed PSA progression, radiographic progression as assessed by CT scan (RECIST 1.1 criteria) and/or bone scan (PCWG2 criteria; ref. 12), unequivocal clinical disease progression, or unacceptable toxicity. The study was designed and conducted in accordance with Good Clinical Practice and the Declaration of Helsinki. All patients provided written informed consent, including for collection and analysis of peripheral whole blood samples. Approval for this study was granted by the University of British Columbia Research Ethics Board (Number H14-00738 and H18-00944). Patients were required to have histologically proven prostatic adenocarcinoma with metastatic disease on CT scan, MRI, and/or bone scan, and a rising PSA (PSA progression per PCWG2 criteria; ref. 12) with castrate levels of testosterone (≤1.7 nmol/L) with ongoing medical castration or previous bilateral orchiectomy. Prior use of CYP17A1 inhibitors or novel AR inhibitors was prohibited, whereas prior use of docetaxel for castration-sensitive disease was allowed. Exclusion criteria included contraindications to abiraterone and enzalutamide per the manufacturer's label and Eastern Cooperative Oncology Group performance status of > 2. The presence of visceral metastasis and pain requiring opioid analgesia were allowed. Patients were required to meet the same eligibility criteria to crossover to second-line treatment.

PSA response was defined as ≥30% PSA decline from baseline confirmed on repeat measurement ≥28 days later. Time to PSA progression (on first- or second-line therapy) was defined as the time from start of therapy to confirmed PSA progression. Confirmed PSA progression was defined as an increase of 2 μg/L and 25% from nadir confirmed by subsequent rising PSA ≥ 28 days later. For patients with no PSA decline, PSA progression was defined as an increase of 2 μg/L and 25% from baseline after ≥12 weeks of treatment. Patients without PSA progression were right-censored at the last treatment date. Overall survival was defined as the time from start of first-line therapy to time of death from any cause, or last follow-up (censored).

TTP (on first- or second-line therapy) was defined as the time from start of therapy to confirmed PSA progression, radiographic progression (PCWG2 criteria), clinical progression, or death from prostate cancer, whichever occurred first. If no events occurred, patients were right-censored at the last treatment date.

All Cox regression analyses, associated confidence intervals (CI), and Kaplan–Meier curves were calculated using R (version 3.6.0) with the “survival” package (version 2.44.1.1). CIs for PSA response rates, Pearson's χ² tests, rank-sum tests, and Fisher exact tests were calculated using Julia (version 1.1.0) with the “HypothesisTests” package (version 0.8.0).

Whole blood was collected in either 4 × 6 mL EDTA tubes (2014–2016) or 2 × 9 mL Streck Cell-Free DNA BCT tubes (2016–2019). Samples collected in EDTA tubes were kept on ice or at 4°C until processing, which occurred within 1 hour. For EDTA tubes, whole blood was aliquoted first and remaining samples were then centrifuged at 1,600 rpm for 10 minutes at 4°C, after which plasma was transferred to a new tube and spun for 10 additional minutes at 1,600 rpm. Samples collected in Streck tubes were maintained at room temperature prior to and during processing. Streck tubes were centrifuged at 1,600 rcf for 15 minutes, buffy coat was aliquoted, and plasma was transferred to a new tube and spun for 10 additional minutes at 5,000 rcf (or maximum attainable speed).

Plasma aliquots were stored at −80°C prior to DNA extraction. cfDNA was extracted from up to 6 mL of plasma with the QIAGEN Circulating Nucleic Acids kit. We included a 1-hour lysis incubation at 60°C, and cfDNA was eluted in 60 μL water. After extraction, cfDNA was quantified with the Qubit 2.0 Fluorometer and Qubit dsDNA HS Assay Kit, or the Quantus Fluorometer and QuantiFluor ONE dsDNA system. Buffy coat fraction aliquots from centrifuged whole blood samples (Streck) and whole blood aliquots (EDTA) were stored at −80°C prior to DNA extraction. White blood cell (WBC) DNA was extracted from the buffy coat fraction using the QIAGEN DNeasy Blood and Tissue kit or Promega Maxwell RSC Blood DNA kit and Maxwell RSC system, as per the manufacturer's instructions. WBC DNA was quantified with a NanoDrop spectrophotometer following extraction.

Sequencing libraries were prepared from a minimum of 10 ng of plasma DNA input per sample. WBC DNA was sheared into ∼180 bp fragments (with a Covaris focused ultrasonicator or via enzymatic digestion) prior to A-tailing, end repair, Illumina-compatible adapter ligation, and PCR amplification. Library quantification was carried out with a NanoDrop spectrophotometer, and each library was run on an ethidium bromide gel to confirm success. Purified sample libraries were multiplexed to obtain single pools with a combined mass of 1 to 2 μg. For targeted sequencing with the 72-gene or AR locus panels, library pools were hybridized to custom-designed Roche NimbleGen SeqCap EZ Choice capture probe sets. For whole-exome sequencing (WES), library pools were hybridized to the SeqCap EZ MedExome panel. SeqCap EZ HyperCap protocol was followed for hybridization and subsequent wash, recovery, and amplification of the captured regions. Agencourt AMPure beads enabled purification of library pools prior to quantitation with the Qubit 2.0 or Quantus Fluorometer. Pools were diluted and sequenced on an Illumina MiSeq (V3 600 cycle kit) or HiSeq 2500 (V4 250 cycle kit) machine. All deidentified targeted and whole-exome DNA sequencing data have been deposited in the European Genome-phenome Archive (EGA) database under the accession code EGAS00001005318 and are available under standard EGA-controlled release.

For both targeted sequencing and WES data, paired-end reads were aligned against an hg38 reference genome using Bowtie-2.3.0 (13). PCR and optical duplicate reads were marked using samblaster-0.1.24 (14) and omitted from analysis. Adapters in read 3′ ends were trimmed in paired-end mode using cutadapt-1.11 (15). Low-quality read tails (smoothed baseq < 30) were trimmed using an in-house algorithm. Per-base read coverages in target regions were quantified using bedtools-2.25.0 (16). cfDNA and WBC sample identities were verified based on SNP genotype profiles.

In targeted sequencing data (72-gene panel), somatic mutations (single-nucleotide variants and indels) were called in cfDNA samples by searching for variants with at least 8 supporting reads and a mutant allele fraction (MAF) of at least 0.5%. The allele fraction was required to be at least 20 times higher than the average allele fraction of the same allele in all WBC samples (n = 202), and at least 3 times higher than the allele fraction in the patient's WBC sample. The paired WBC sample had to have at least 20 reads covering the position. The comparisons to patient-matched and cohort-wide WBC samples were necessary to help eliminate (from somatic mutations calls) germline SNPs, stereotypical artifacts arising from sequencing and/or alignment error, as well as mutations of plausible clonal-hematopoietic origin. For base substitutions, the average mapping quality of mutation supporting reads was required to be at least 10, and the average distance of the mutant allele from the nearest read end was required to be at least 15 bases. Our mutation analysis pipeline has been previously validated in matched tissue cfDNA cohorts and in dilution series experiments (11, 17). Protein-level consequences of variants were predicted using ANNOVAR, and all somatic variants were inspected with the Integrative Genomics Viewer for alignment errors.

Somatic mutations were called in WES data using the same pipeline but requiring an allele fraction ≥ 5% to account for the lower sequencing depth. The allele fraction was required to be at least 25 times higher than the average allele fraction across all WBC samples, and at least 10 times higher than the allele fraction in the paired WBC sample.

Germline variants were called in targeted 72-gene panel WBC samples by searching for variants with an alternate allele fraction of at least 15%, and at least 8 supporting reads. Germline variants with population allele frequency 0.5% or higher in the KAVIAR or ExAC databases were discarded. Protein-level consequences of variants were predicted using ANNOVAR. Variant pathogenicity was assessed using ClinVar annotations (18) and based on our previously published protocols (19).

Differences in ctDNA fraction and sequencing depth between serial samples can affect sensitivity of mutation detection. To reduce the likelihood of false negatives in comparing mutation presence/absence across samples, mutations independently called in one cfDNA sample were considered not amenable for detection in a later cfDNA timepoint if the expected allele fraction MAF_earlier * (C_later / C_earlier) was below 2%, where C_later and C_earlier were ctDNA fractions of the later and earlier samples, respectively. In addition, if the probability of observing eight or more mutant reads in the later sample (based on coverage and expected allele fraction) was less than 90%, the mutation was considered not amenable for detection.

For cfDNA samples sequenced with the 72-gene panel, ctDNA fractions were estimated based on the allele fractions of autosomal somatic mutations in nonamplified regions. MAF in autosomes is elevated when a mutation is accompanied with loss of the other allele (i.e., LOH). In this situation, the MAF and ctDNA fraction are related as ctDNA fraction = 2 / (1 / MAF + 1). Because deletions are not detectable at low ctDNA fractions and LOH is very common in frequently mutated genes such as TP53, we conservatively assumed all mutations to have associated LOH. To account for sampling noise, we modeled the mutant read count as arising from a binomial distribution, and conservatively calculated what the true MAF would be if the highest observed MAF was a 95% quantile outlier. After calculating a ctDNA fraction estimate for each somatic mutation, the highest estimate was adopted as the estimate for the sample, as this mutation was the most likely to be truncal to the metastatic lineage. A diploid tumor model was assumed, because our targeted panel lacks sufficient genome breadth to reliably identify whole genome duplication events. These three assumptions are potential sources of error in calculating ctDNA purity from targeted data. Nonetheless, here and in our prior prostate cancer work (11), we have shown that ctDNA fraction estimates derived from targeted 72-gene panel data are highly correlated with orthogonal approaches based on whole-exome copy-number profiles. Our methodology for approximating ctDNA fraction has also been applied to targeted cfDNA sequencing data from bladder cancer (20).

For samples with WES, ctDNA fractions were estimated based on genome-wide copy-number and heterozygous SNP allele fractions (Supplementary Figs. S1 and S2). Specifically, we divided the genome into segments based on coverage logratio discontinuities identified in any sample from a patient (see Patients and Methods section “Copy-number segmentation”). We then tested combinations of ctDNA fraction and tumor copy-number levels to calculate expected autosomal logratios LR = log₂[(F * C + (1 – F) * 2) / 2] (where F = ctDNA fraction and C = number of mutant allele copies), and selected the model that minimized discrepancy between the expected and observed exome-wide autosomal logratios. Candidate models were rejected if any regions had a negative copy number, or if the model suggested a zero copy status for a large genomic region (which is highly biologically improbable). The best model was selected manually to reduce error arising from aneuploidy and potential polyclonality captured by ctDNA; prior experience with automated tools has suggested low accuracy for cfDNA whole-exome data. We corroborated our WES ctDNA fraction estimates using ASCAT (version 2.5.2; ref. 21) and achieved a Pearson correlation of 0.88 between the manual and automated estimates for the 48% of samples where ASCAT was able to find a model.

To model the accuracy of our ctDNA fraction estimates, we compared estimates from the mutation- and copy-number–based procedures using a Bland–Altman plot, and found that a model C = C_true * 2^{N(0, 0.12)} + N(0, 0.06) composed of multiplicative noise ∼ 2^{N(0, 0.12)} and additive noise ∼ N(0, 0.06) represented the data well (Supplementary Fig. S3).

Read counts inside panel target regions were calculated using bedtools-2.25.0 (16). A set of 138 reference cfDNA patient samples with no detectable somatic mutations or amplifications from targeted sequencing were used to generate an aggregate reference coverage profile. To do this, we first calculated coverage logratios L_k = log₂(R_k,s / R_k,1) for each reference cfDNA sample against an index sample (arbitrarily selected as the first reference sample), where R_k,s represented the read count of capture region k in sample s. We then fitted a Loess regression curve to a two-dimensional scatter plot of (G_k, L_k) pairs—where G_k was the GC fraction of region k—and subtracted the Loess curve from the logratio values to correct for GC bias and differences in global sequencing depth. By inverting the aforementioned formula, we then converted each sample's GC-corrected logratios back to read counts. Finally, we constructed the aggregate reference coverage profile by calculating the median positional coverage of all 138 GC-corrected reference samples. The cfDNA samples used as negative controls were sequenced with the same panel, and to a similar depth.

To calculate coverage logratios for patient cfDNA samples, each sample's coverage profile was compared against the reference profile, and the resulting logratios were GC corrected using Loess regression. A gene's coverage logratio was calculated as the median logratio of all target regions inside the gene.

For samples sequenced with the 72-gene panel, gene deletions and gains were called at coverage logratios of ≤ –0.3 and ≥ 0.3, respectively. These thresholds were determined by examining the distribution of coverage logratios and heterozygous SNP allele fractions in negative control cfDNA samples from healthy volunteers and patients with ctDNA-negative prostate cancer, as previously described (11).

A gene's copy number in the cancer cells contributing to ctDNA was calculated as P * [1 + (2^L − 1) / F], where L = the gene's coverage logratio, P = the gene's copy number in normal cells, and F = ctDNA fraction. To derive a CI for a gene's copy number, we assumed that the gene's coverage logratio was normally distributed as L ∼ N(μ, σ²), and used cfDNA samples sequenced with both the 72-gene and whole-exome panel to estimate σ² = Var(L_72G − L_WES) / 2, where L_72G and L_WES were coverage logratios measured with the two panels. For simplicity, this calculation assumed that coverage logratios from the two panels were identically distributed.

To identify temporal copy-number changes between same-patient samples sequenced with the whole-exome or AR panels, we used a segmentation procedure incorporating information from all consecutive samples. Each chromosome was represented as a sequence of coverage logratios L₁, …, L_n ordered by genomic position. To identify a breakpoint a, we searched for integers a and b defining a region L_a, …, L_b such that the statistical difference relative to region L₁, …, L_a-1 was maximized. The difference was defined as Z = (μ_1..a − μ_a..b) / sqrt(1 / a + 1 / (b − a)), where μ_a..b = mean (L_a, …, L_b). For each chromosome, the breakpoint position a that produced the highest Z score in any sample from the patient was identified, and the chromosome was split into two segments if the Z score exceeded a specified threshold. This process was repeated recursively for each resulting segment.

To account for variable levels of coverage logratio noise between samples, the Z scores calculated for each sample were divided by a sample noise score N, defined as the average difference between consecutive logratios.

Somatic genomic profiles from targeted sequencing approaches such as our 72-gene panel do not enable robust characterization of subclonal cancer cell populations, as insufficient somatic mutations are captured. However, a shift in the relative fractions of cancer cell populations between two samples can be inferred even from sparse evidence. Here, we considered a population shift to have occurred between two cfDNA samples sequenced with the 72-gene panel if at least one somatic mutation pair showed a significant change in allele fraction relative to each other, or if at least one gene showed a significant change in its ctDNA fraction adjusted copy number.

To identify a change in the relative allele fraction of two mutations A and B between two cfDNA samples, we started with the null hypothesis p_1A / p_1B = p_2A / p_2B, where p_1A was the true allele fraction of mutation A in sample 1. This null hypothesis states that the relative allele fraction of mutations A and B did not change between samples. We then ran a million rounds of simulation generating allele fractions s_1A ∼ Bin(n_1A, o_1A) / n_1A, s_1B ∼ Bin(n_1B, o_1B) / n_1B, s_2A ∼ Bin(n_2A, cs_2A) / n_2A, and s_2B ∼ Bin(n_2B, cs_2B) / n_2B, where o_1A = observed allele fraction of mutation A in sample 1, n_1A = observed total read depth at the position of mutation A in sample 1, and c = observed change in average allele fraction between the samples (o_2A + o_2B) / (o_1A + o_1B). We then tested how often the observed relative allele fraction o_2A / o_2B was more extreme than the simulated s_2A / s_2B. This gave us the probability of the observed change in relative allele fraction occurring due to sampling noise. A population shift was called if this P value was less than 0.00001 for any mutation pair. Samples with ctDNA < 5% were excluded, based on a previous dilution series experiment where mutation detection sensitivity in this cohort approached 100% for samples with ctDNA ≥ 5% (11). If the probability parameter for any binomial distribution would have been zero, it was instead replaced with 1 / read depth, in accordance with Cromwell's rule.

To identify a change in the copy number of gene G between two cfDNA samples, we started with the null hypothesis that the gene had the same copy number (C₁ + C₂) / 2 in both samples. We then ran a million rounds of simulation where we perturbed the observed ctDNA fractions with noise (Supplementary Fig. S3), then calculated coverage logratios by inverting the C = P * [1 + (2^LR − 1) / F] formula, then perturbed the logratios with noise ∼ N(0, S * σ_G), converted back to copy numbers, and stored the absolute difference of the two simulated copy-number measurements. The amount of noise added to logratios was based on σ_G = standard deviation of the gene's coverage logratios in control cfDNA samples with zero ctDNA, and S = sample noise multiplier determined by comparing coverage logratios of adjacent genomic regions. A copy-number change was called if less than 1% of simulations (P < 0.01) yielded a difference stronger than the observed difference |C₁ - C₂|.

In pairs of samples with WES data, we used a different method to identify mutation-based evidence for a shift in cancer cell populations, because the number of tested hypotheses was too large when testing all mutation pairs. We considered a mutation to show evidence for population shift between two WES samples if the ctDNA fraction adjusted MAF was different (Bonferroni adjusted P < 0.00001) between two cfDNA timepoints according to Fisher exact test, and the MAF was higher than 10% in one sample. To adjust for ctDNA fraction in the two samples, the number of mutant reads in each timepoint was divided by C / A, where C = ctDNA fraction of the sample and A = average ctDNA fraction of the two samples. Unlike the method used for the 72-gene panel, this method relies on accurate knowledge of the ctDNA fraction for the two samples.

Somatic structural rearrangements were detected using Breakfast (github.com/annalam/breakfast), with the –merge-duplicates option. Briefly, breakfast extracts 30 bp anchors from the 5′ and 3′ end of unaligned reads and aligns them against the genome. If both anchors align to the genome, an optimal breakpoint inside the original read is sought such that the halves align to their respective anchor regions with the fewest mismatches. Split-read alignments supporting an identical junction structure were grouped together. Somatic rearrangements were required to be supported by at least four unique junction-spanning reads, and to have no supporting reads in any of the 202 WBC samples. All rearrangement junction sequences were inspected using UCSC BLAT, and were discarded if alternative genomic alignments in homologous regions were detected. A rearrangement was classified as a deletion if the 5′ anchor was located upstream of the 3′ anchor in the same strand of the same chromosome, and as a tandem duplication if the 5′ anchor was located downstream of the 3′ anchor. A rearrangement was classified as gene inactivating if it was a deletion or tandem duplication overlapping the gene's coding region, or an inversion or translocation anywhere between the transcription start site and stop codon.

A rearrangement was classified as an AR gene structural rearrangement (AR-GSR) if at least one breakpoint fell inside the AR gene body (chrX: 67,544,623–67,730,619 in hg38). An AR-GSR was considered LBD truncating if it was (Supplementary Fig. S4):

• 1) A deletion that fully or partially overlaps the LBD in AR exons 5 to 8, but does not overlap AR exons 1 to 3

• 2) A tandem duplication that fully overlaps AR exons 1 to 3, but does not fully overlap the LBD in AR exons 5 to 8

• 3) A tandem duplication that starts downstream of AR exon 3 and ends before the stop codon in exon 8

• 4) An inversion or translocation with at least one breakpoint located between the end of AR exon 3 and stop codon in exon 8, with a junction flank on the centromeric side containing AR exons 1 to 3

To determine which AR-GSRs showed evidence for an increase in their allele fraction between consecutive cfDNA samples, we calculated the number of unique rearrangement supporting reads R₁ and R₂ in each sample, and divided them by the median AR sequencing depth A₁ and A₂ in those samples to account for differences in global sequencing depth between samples.

In our published trial, we randomized 202 patients with mCRPC to either abiraterone (an inhibitor of androgen synthesis) or enzalutamide (an AR LBD antagonist), followed by crossover to the other AR signaling inhibitor at progression (7). Patients provided serial plasma cfDNA samples before protocol treatment (baseline; representing clinically progressing mCRPC) and at progression on first- and second-line therapy with either inhibitor (Fig. 1A; Supplementary Table S1; Supplementary Fig. S5). In total, 458 samples were collected and subjected to initial targeted sequencing (median read depth 724x) with a panel encompassing 72 prostate cancer genes to estimate ctDNA quantity and to inform on DNA-level temporal changes in driver genes (Supplementary Tables S2 and S3).

Somatic mutations with an allele fraction above 1% were detected in 119 of 202 (59%) baseline samples, 82 of 158 (52%) first progression samples, and 67 of 98 (68%) second progression samples (Supplementary Table S4). Using the highest allele fraction mutation in nonamplified regions, estimated ctDNA fractions (as a proportion of total cfDNA) ranged between undetectable and 89%, and were strongly correlated with clinical markers of tumor burden (Supplementary Table S5; Supplementary Fig. S6). Select samples (mainly those with high ctDNA fractions, but including some low ctDNA fraction controls) were additionally subjected to targeted sequencing of the entire AR gene and flanking regions (n = 175; median unique read depth 1,093x) and/or WES (n = 134; median depth 144x; Fig. 1A; Supplementary Table S3). Technical concordance was high between targeted sequencing and WES (Fig. 1B; Supplementary Fig. S7). ctDNA fractions estimated based on whole-exome somatic profiles were consistent with mutation-based estimates from targeted sequencing (Fig. 1C; Supplementary Fig. S8).

ctDNA fractions between consecutive cfDNA collections were correlated (Fig. 1D). The correlation was stronger when timepoints were separated by < 180 days (the median collection interval; Pearson r = 0.75 vs. 0.49). Only 21 of 63 (33%) patients with multiple samples and ctDNA fraction < 2% at baseline had ctDNA ≥ 2% in a subsequent progression sample. Changes in ctDNA fraction between collections were correlated with changes in serum PSA, a clinical surrogate for depth of therapy response (ref. 22; Fig. 1E; Supplementary Fig. S9). Eleven cfDNA samples were collected too early (before progression) due to clinical circumstances, and expectedly had lower ctDNA fractions than all other samples (median undetectable vs. 3.9%, P = 0.008, rank-sum test). Collectively, these results suggest that repeat ctDNA testing for genomic characterization may have limited utility in patients with consistently low tumor burden.

Patients with baseline ctDNA ≥ 5% (the cohort median) exhibited shorter overall survival (median 16.6 months vs. not reached; HR, 4.57; 95% CI, 3.01–6.92; P < 0.001) and time from start of first-line therapy to second progression (median 7.8 vs. 19.1 months; HR, 2.63; 95% CI, 1.88–3.67; P < 0.001) than patients with low or undetectable ctDNA. This association was quantitative, with higher ctDNA fractions associated with progressively worse overall survival and TTP (Fig. 1F; Supplementary Fig. S10). Importantly, baseline ctDNA fraction achieved higher prognostic accuracy than a multivariate model combining all routine clinical prognostic factors (Supplementary Fig. S11). In our trial, almost no patients benefited from second-line abiraterone, but 51% of patients experienced a PSA response of 30% or greater to second-line enzalutamide (Supplementary Fig. S5; ref. 7). Two characteristics showed the strongest association with TTP on second-line enzalutamide: ctDNA fraction at crossover and TTP on first-line abiraterone (Fig. 1G and H; Supplementary Table S6).

To compare driver gene alteration status over time, we focused on 164 same-patient sample pairs (across 66 patients) with sufficient ctDNA for robust genotyping (ctDNA ≥ 5%; see Patients and Methods for justification). In these pairs, ctDNA fraction adjusted MAFs and copy numbers were strongly correlated (Fig. 2A and B). We observed high temporal concordance for the copy number, rearrangement, and mutation status of key prostate cancer driver genes associated with tumor initiation and development, including SPOP, FOXA1, PTEN, and TP53 (Fig. 2C and D; Supplementary Table S7). Across all regions captured via targeted sequencing, 78% of somatic mutations detected in an earlier sample were independently redetected in the subsequent sample. The rate of redetection increased to 91% after excluding discordances explained by insufficient sequencing depth and/or reduced ctDNA fraction in the later sample (Supplementary Fig. S12). Mutations estimated to be present in over 50% of ctDNA-shedding cells were redetected more often than less dominant mutations (93% vs. 76%, P < 0.0001, Fisher exact test). Together, these data demonstrate that pretreatment and postprogression ctDNA populations share a large common reservoir of cancer driver alterations.

Because rapid autopsy studies have revealed interlesion genomic heterogeneity in mCRPC (23), we hypothesized that treatment may drive temporal changes in the proportions of genomically distinct cancer cell populations represented in ctDNA. Therefore, we searched for quantitative differences in somatic MAFs and gene copy numbers in the 164 sample pairs with ctDNA ≥ 5% (see Patients and Methods). Using targeted sequencing alone (72-gene panel), 36% of patients (24/66) displayed evidence for a temporal shift in detected cancer cell populations (Supplementary Table S8). MAF profiles changed significantly in 21 patients (Fig. 3A; Supplementary Table S9), and copy-number profiles shifted in 17 patients (with an overlap of 12 patients; Fig. 3B, Supplementary Fig. S13).

To validate across a broader genomic context, we analyzed serial cfDNA WES data from 36 patients, focusing on patient-matched sample pairs where high ctDNA fraction (≥30%) enabled more sensitive comparisons. Among patients with evidence for population shifts from initial targeted sequencing, exome-wide changes in mutation and copy-number profiles were confirmed in all evaluable patients (Fig. 3C–E; Supplementary Figs. S14 and S15). In one extreme case (patient 180), we observed an almost complete switch in ctDNA populations between baseline and second progression, and a mixture of the two distantly related populations in the first progression sample (Fig. 3D). Whole-exome analysis supported the absence of ctDNA changes in patients classified as temporally similar based on targeted sequencing data (Supplementary Figs. S16–S18). One patient (026) was excluded from analysis because his mCRPC baseline and “progression” samples shared no mutations, and he was diagnosed with concurrent metastatic bladder cancer (Supplementary Fig. S19).

With the exception of the AR (addressed below), the only recurrent posttreatment changes detected by our 72-gene panel were WNT pathway changes (n = 3 patients), PI3K pathway alterations (n = 3), and cell-cycle gene changes (n = 4, including two patients that acquired a CCND1 amplification at progression; Supplementary Fig. S20).

Missense mutations in the AR LBD result in differential ligand affinity, and prior work has suggested genotype switching in response to enzalutamide or abiraterone (24, 25). Accordingly, we observed high variation between collection timepoints for somatic AR LBD mutations (Supplementary Fig. S21). We detected a marked increase in AR L702H mutations previously linked to abiraterone resistance through glucocorticoid agonism (glucocorticoids are administered concurrently with abiraterone; Fig. 4A; refs. 24, 26). Interestingly, although 6 patients acquired AR L702H mutations after exposure to abiraterone plus glucocorticoids, 2 patients (035 and 163) acquired L702H mutations after enzalutamide, without accompanying glucocorticoids (Supplementary Fig. S22). Both patients harbored additional AR LBD mutations, and patient 163 harbored a mismatch repair defect and hypermutation (27). The AR F877 L mutation linked with enzalutamide resistance in model systems was not identified in any samples (28, 29).

Mutations in AR amino acids 742 to 743 were found in 6 patients, all of whom had previously received bicalutamide, a first-generation AR signaling inhibitor. These mutations are a known resistance mechanism to bicalutamide, and showed decreased frequency after enzalutamide and/or abiraterone treatment (Fig. 4A), consistent with preclinical evidence of enzalutamide activity against these mutants (30).

Somatic AR amplifications were detected in 62 of 119 (52%) patients at baseline and 74 of 149 (50%) progression samples (excluding samples with undetectable ctDNA; Supplementary Fig. S21). Given the switch in AR mutation genotypes, we hypothesized that AR amplifications detected before and after treatment may constitute different events in some patients. Nine of 66 patients with sufficient ctDNA in two or more samples to quantify AR copy number displayed a significant temporal increase in AR copy number per average ctDNA-shedding tumor cell (Fig. 4B). Three patients displayed a significant decrease in AR copy number. Consistent with these observations, across all evaluable patients, mean AR copy number increased between the three collection timepoints (7.8 vs. 10.1 vs. 12.3 copies; Fig. 4C). Because an increased AR copy number can result in elevated AR protein expression and greater ligand sensitivity (31), populations with higher AR copy number may have greater fitness during treatment. Accordingly, higher baseline AR copy number was associated with shorter overall survival, even in multivariate models incorporating ctDNA fraction (Fig. 4D; Supplementary Table S6).

To explore the genomic context of AR copy-number increase, we analyzed a 30 Mbp window around AR in 175 cfDNA samples (see Patients and Methods). Consistent with our initial observation of temporally increasing AR gene copy number, a spatial profile of average chromosome X copy-number change between consecutive timepoints revealed a focal amplification centered at the AR and its upstream region (Fig. 4E). All evaluable patients with AR copy-number change also displayed a change in chromosome X amplicon structure (Fig. 4F; Supplementary Fig. S23). The observed changes in amplicon boundaries suggest that AR copy-number increases are due to treatment-induced selection for a distinct cell population with a different pattern of AR gene copy gain. Interestingly, an additional 9 patients displayed a change in AR amplicon structure without significant change in bulk AR copy number, consistent with the notion that patients harbor multiple similar AR locus genotypes.

An enhancer region 640 kb upstream of AR was recently reported to contribute to increased AR expression through focal tandem duplications in patients with mCRPC (32, 33). In our cohort, most AR amplifications also encompassed the enhancer (Fig. 4G). Indeed, 15 patients exhibited AR locus copy gains where the amplification peak began at the AR enhancer (Supplementary Fig. S24). Only 7 patients showed an AR amplification that excluded the enhancer. Fourteen patients showed evidence for a focal amplification of the AR enhancer. Together, these data suggest that potent AR signaling inhibition precipitates population restructuring that strongly favors lineages with high-level AR gene and enhancer amplification.

AR proteins with a truncated LBD can facilitate target gene transcription independent of ligand availability, and can drive resistance to approved AR signaling inhibitors in preclinical models (34). LBD-truncated AR transcripts can be produced via aberrant splicing, but recent work suggests that genomic breakpoints falling within the AR intragenic region can also give rise to truncated transcripts (3, 35). Here, we identified AR-GSRs in 69 of 175 (39%) cfDNA samples subjected to targeted sequencing of the entire AR locus (Supplementary Table S10). All detected AR-GSRs were unique to a single patient. Of the 170 distinct AR-GSRs identified across all timepoints, 36% were intragenic, 50% juxtaposed AR with another chromosome X region, and 14% were interchromosomal. Copy-number profiles inside AR were consistent with detected AR-GSR breakpoints: among 8 patients with clear evidence for an intragenic AR copy-number change, a concordant AR-GSR was identified in all except one (077, where the expected breakpoint was inside a LINE element and therefore unresolvable; Supplementary Fig. S25). These data are consistent with the AR being one of the most frequently rearranged genes in mCRPC (4, 35).

Based on structural analysis (see Patients and Methods; Supplementary Fig. S4), 34 unique AR-GSRs had potential to produce a truncated AR isoform lacking part or all of the LBD (exons 5–8) while preserving the DNA binding domain (exons 1–3; Supplementary Table S10). The number of detected AR-GSRs was strongly correlated with AR copy number (Fig. 5A). cfDNA samples with amplified AR were 9 times more likely to carry AR-GSRs (61% vs. 7%, P < 0.001, Fisher exact test), including LBD truncating rearrangements (26% vs. 3%, P < 0.001, Fisher exact test).

Recent studies have linked LBD truncating AR-GSRs to de novo treatment resistance (11, 35, 36), but it is unknown whether AR signaling inhibitors actively select for these genotypes in responding patients. Here, the fraction of patients carrying AR LBD truncating rearrangements increased successively with progression on AR signaling inhibitors (11% at baseline vs. 20% at first progression vs. 25% at second progression; Fig. 5B). Comparing coverage-normalized read counts at consecutive timepoints, we observed that although 72% of all AR-GSRs increased in allele fraction over time, LBD truncating rearrangements were under stronger positive selection than other AR-GSRs (P = 0.003, rank-sum test; Fig. 5C). Eight patients displayed new LBD truncating rearrangements in progression samples, with no supporting evidence at baseline despite a median ctDNA fraction of 41% and median sequencing depth of 1,388x at the breakpoints. Conversely, no LBD truncating rearrangements decreased to undetectable levels over time, compared with 13% of other AR-GSRs (P = 0.0036, Fisher exact test; Fig. 5C).

Intriguingly, we observed an enrichment for intragenic deletions removing AR exons 5 to 7 while preserving the final exon 8 (Fig. 5D and E). In patients 071 and 110, multiple independent deletions matching this pattern were detected solely in progression samples. Because deletions of exons 5 to 7 juxtapose frame-incompatible exons 4 and 8 (which if spliced together would yield a frameshifted AR isoform), their apparent positive selection may suggest a benefit for preservation of the 3′-untranslated region (and polyadenylation site) or adjacent regulons.

Combining results from all analysis approaches, treatment-induced genomic changes were detected in 53% of evaluable patients (Supplementary Table S8). The vast majority (86%) of these patients exhibited AR locus changes, suggesting that differential selective pressure on distinct AR genotypes may be driving the overall population shifts (Fig. 6A). In support of this concept, 7 patients exhibited temporal changes only at the AR locus. Patients with primary resistance against first-line treatment (TTP < 3 months) did not show systematic increases in AR copy number (Fig. 6B), suggesting that most changes were occurring in the context of disease that was initially sensitive to treatment. Consistent with this, patients with evidence for temporal changes in ctDNA populations were more likely to have experienced a deep PSA response to intervening therapy than other patients (median PSA response 76% vs. 28%, P = 0.028; Fig. 6C). Indeed, all 8 patients with greater than 90% PSA decline and ctDNA fraction ≥ 10% in serial cfDNA collections showed population shifts.

This study provides insights into mCRPC genomic evolution in living patients as their disease responds and progresses through standard-of-care treatment with abiraterone and enzalutamide. Importantly, observed shifts in the predominant plasma ctDNA profile and AR genotype over sequential clinical progressions propose a Darwinian explanation for acquired resistance to these widely used AR signaling inhibitors.

Somatic profiles changed over time in most patients where AR signaling inhibitors elicited strong responses. Differences in the relative abundance of mutations, copy numbers, and rearrangements as a proportion of total ctDNA supported a model in which the metastatic niche harbors several phylogenetically related competing populations. This is consistent with genomic analyses of metastatic lesions obtained at autopsy, which show interlesion heterogeneity at the time of patient death (23).

Very few temporal changes were detected in TP53, FOXA1, or SPOP, supportive of prior work indicating that somatic alterations in these genes are early events truncal to the metastatic lineage (23, 37, 38). However, we observed a striking convergence of genomic changes upon AR. Compared with baseline mCRPC, the temporal changes detected here in AR after progression on abiraterone and/or enzalutamide reflect a shift toward even higher levels of AR augmentation. Although the concept of an AR genotype that changes during mCRPC treatment is well accepted for AR LBD mutations (39), previous efforts to understand the clinical relevance of AR amplification have not explored whether amplifications detected at different timepoints represent the same event (24, 25, 40). Here, the mean AR copy number in ctDNA increased by 50% after progression on abiraterone and enzalutamide. Furthermore, drastic temporal differences in AR locus structure within individual patients suggested that stronger AR amplification after AR-targeted therapy often arose through selection of a different tumor population relative to baseline mCRPC. This treatment-associated pressure toward increased AR copy number may help explain why most mCRPC with AR amplification initially respond to potent AR inhibitors, and is consistent with observations that AR amplification associates with poor prognosis (41, 42).

Our work promotes LBD truncating AR-GSRs as a major driver of clinical resistance to AR signaling inhibitors. Robust in vitro work has underscored the potential of these genomic variants to drive resistant phenotypes via expression of constitutively active AR protein isoforms (3, 35), but the alternative hypothesis that AR-GSRs represent a neutral side effect of AR amplification has never been ruled out. Here using serial sampling in a live patient population, we show that the frequency of LBD truncating AR-GSRs increases 2.5-fold after treatment with abiraterone/enzalutamide, and that other AR rearrangements are under weaker selection. Combined with mutation and copy-number observations, the consistent change toward more active AR paints a picture of powerful selective pressure toward hyperaugmented AR and emergence of lethal cross-resistance.

The emergence of somatic populations with heightened AR augmentation supports the concept that a large proportion of resistant disease retains some reliance upon the AR, even after potent pathway suppression. As such, renewed targeting of the AR signaling pathway may still have an antitumor effect, if therapeutic strategies are independent of the AR LBD. Nevertheless, it may be more beneficial to temporarily reduce the selective pressure on the AR. This concept is supported by a recent phase III clinical trial demonstrating that second-line taxane-based chemotherapy outperforms treatment with sequential AR signaling inhibitors (6). Correlative studies have also suggested that patients with AR copy gain in plasma cfDNA may respond better to chemotherapy than AR signaling inhibitors (43), supporting the general hypothesis that chemotherapy is less dependent on AR genotype.

The overall abundance of ctDNA is increasingly recognized as an impactful prognostic factor in metastatic cancers (44). We and others have reported on the potential for ctDNA fraction to add to clinical prognostic models in mCRPC (11, 45). Although the current study was not focused on biomarker discovery, we note that ctDNA fraction may help identify patients likely to be responsive to second-line enzalutamide. This concept will be tested in a biomarker-driven phase II trial (NCT04015622).

Our analysis of treatment-driven genomic changes was limited to patients with intermediate outcomes, as sequential collections were not available from patients whose disease responded exceptionally badly (i.e., rapid death) or exceptionally well (i.e., no progression). Among patients with sequential collections, we adjusted for temporal changes in ctDNA fraction, but were forced to exclude samples with low ctDNA fraction, introducing a potential source of bias. Furthermore, this study focused on known driver genes and coding regions, limiting our ability to identify novel changes in the noncoding genome (46), or comprehensively resolve and characterize individual subclonal populations. Future studies should utilize whole-genome sequencing on samples with high ctDNA fractions to inform on clonal structure.

M. Annala reports grants from Jane and Aatos Erkko Foundation and Academy of Finland during the conduct of the study. S. Taavitsainen reports grants from Jane and Aatos Erkko Foundation during the conduct of the study. D.J. Khalaf reports personal fees from Janssen outside the submitted work. D.L. Finch reports other support from Janssen, Astellas, and Bayer outside the submitted work. J. Vergidis reports personal fees from BMS and non-financial support from Bayer outside the submitted work. B.J. Eigl reports personal fees from AstraZeneca, Pfizer, Seagen, Johnson and Johnson, EMD Serono, and Merck outside the submitted work. C.K. Kollmansberger reports personal fees from Astellas, Janssen, Pfizer, BMS, Eisai, Merck, Ipsen, and Bayer outside the submitted work. M. Nykter reports grants from Jane and Aatos Erkko Foundation and Academy of Finland during the conduct of the study. K.N. Chi reports grants from Astellas and Janssen during the conduct of the study. K.N. Chi also reports grants and personal fees from AstraZeneca, Bayer, Novartis, Pfizer, Point Biopharma, Roche, and Sanofi, as well as personal fees from Daiichi Sankyo, Merck, and Bristol Myers Squibb outside the submitted work. A.W. Wyatt reports personal fees from AstraZeneca, Merck, and Astellas; grants and personal fees from Janssen; and grants from ESSA Pharma outside the submitted work. No disclosures were reported by the other authors.

M. Annala: Conceptualization, data curation, software, formal analysis, investigation, visualization, methodology, writing–original draft, writing–review and editing. S. Taavitsainen: Data curation, formal analysis, investigation, methodology. D.J. Khalaf: Data curation, formal analysis, investigation. G. Vandekerkhove: Investigation, methodology, writing–review and editing. K. Beja: Investigation, methodology. J. Sipola: Data curation, software, methodology. E.W. Warner: Investigation, methodology. C. Herberts: Investigation, methodology, writing–review and editing. A. Wong: Investigation. S. Fu: Investigation. D.L. Finch: Investigation. C.D. Oja: Investigation. J. Vergidis: Investigation. M. Zulfiqar: Investigation. B.J. Eigl: Investigation. C.K. Kollmansberger: Investigation. M. Nykter: Supervision, project administration. M.E. Gleave: Investigation, project administration. K.N. Chi: Conceptualization, resources, supervision, funding acquisition, investigation, writing–original draft, project administration, writing–review and editing. A.W. Wyatt: Conceptualization, resources, data curation, supervision, funding acquisition, methodology, writing–original draft, project administration, writing–review and editing.

This work was supported by the Canadian Institutes of Health Research, Canadian Cancer Society Research Institute, Prostate Cancer Canada, Movember Foundation, Prostate Cancer Foundation, Jane and Aatos Erkko Foundation, Academy of Finland, Terry Fox New Frontiers Program, BC Cancer Foundation, Janssen, and Astellas. No funding sources were involved in the design or execution of the study. The authors are thankful to all contributing patients and their families.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.