Intratumor heterogeneity and clonal evolution revealed in castration-resistant prostate cancer by longitudinal genomic analysis.

Intratumor heterogeneity is a key driver for local relapse and treatment failure. Thus, using multifocal prostate cancer as a model to investigate tumor inter-clonal relationships and tumor evolution could aid in our understanding of drug resistance. Previous studies discovered genomic alterations by comparing hormone-sensitive prostate cancer (HSPC) with castration-resistant prostate cancer (CRPC) in large cohorts. However, most studies did not sequentially sample tumors from the same patient. In our study, we performed whole-exome sequencing (WES) on 14 specimens from five locally relapsed patients before and after androgen-deprivation therapy. We described the landscape of genomic alterations before and after treatment and identified critical driver events that could have contributed to the evolution of CRPC. In addition to confirming known cancer genes such as TP53 and CDK12, we also identified new candidate genes that may play a role in the progression of prostate cancer, including MYO15A, CHD6 and LZTR1. At copy number alteration (CNA) level, gain of 8q24.13-8q24.3 was observed in 60% of patients and was the most commonly altered locus in both HSPC and CRPC tumors. Finally, utilizing phylogenetic reconstruction, we explored the clonal progression pattern from HSPC to CRPC in each patient. Our findings highlight the complex and heterogeneous mechanisms underlying the development of drug resistance, and underscore the potential value of monitoring tumor clonal architectures during disease progression in a clinical setting.

Introduction

Genomic instability is the engine for tumorigenesis [1]. Characterized by uncontrolled growth and cell division, tumors expand and accumulate heterogeneous genomic alterations. This intra-tumor heterogeneity is one of the main culprits for local relapse and treatment failure [2]. Therefore, deciphering the clonal evolutional path of drug-resistant clones could provide insights for biomarker discovery and development of personalized treatment [3].

Prostate cancer is the most common cancer and the second leading cause of cancer death among men in the United States [4] and is markedly heterogeneous in terms of pathology and clinical presentation. Multiple tumor foci are commonly detected within a single prostate, underscoring the scale of heterogeneity in the tumor of the same patient [5]. Androgen-deprivation therapy (ADT) is fundamental in the treatment of the majority of prostate cancer patients. However, 60–70% patients will progress to castration-resistant prostate cancer (CRPC) within 2,3 years following the initiation of ADT. After progression to CRPC, the estimated mean survival of men is 9–36 months [6]. Although previous studies reported main driver genes such as AR amplifications, TP53 and RB1 deletions in CRPC samples, they mainly analyzed the metastatic lesions, but did not perform sequential analysis of tumors from the same patient. Furthermore, most previous studies focused on distant metastasis, leaving a significant knowledge gap in characterizing the transition from local HSPC to local CRPC in the same patient. In this study, we investigated genomic changes associated with clonal evolution in locally relapsed patients. We performed whole-exome sequencing (WES) on multiple tumor foci both at the treatment-naïve and castration-resistant stages within the same individuals. Comparison of mutational profiles within primary tumors and between primary tumors and drug-resistant tumors revealed extensive intratumor heterogeneity. Clonal evolution analysis revealed diverse patterns of how tumors progress from HSPC to CRPC, providing insights for more precise therapeutic intervention.

Materials and methods

Study design and sample information

A total of 5 CRPC samples and 9 matched treatment naïve tumor samples were obtained from prostate FFPE biopsy in 5 patients. Multiple tumor foci were profiled from all patients in HSPC, but two of the 5 cases (Patient-L and Patient-Y) were selected for only one biopsy to perform sequencing, owing to low tumor content. The blood samples were collected from each patient at the time of initial diagnosis as a germline control. Patients understood the procedures and implications of the study and provided written informed consent. Sample collection for this study was approved by the ethical committee at Changhai Hospital (CHEC2019-012). All the HSPC patients received ADT therapy in accordance with the guideline [7]. When necessary, transurethral resection of the prostate (TURP) was performed to relieve the symptoms of dysuria necessary. CRPC was defined as castrate levels of serum testosterone < 50 ng/dL or 1.7 nmol/L plus one of the biochemical progression or radiological progression according to EAU guidelines [7]. Follow-up database PC-FollowTM [8] was used to deposit all patient clinical information.

DNA extraction and quantification

For all 5 cases, tumor DNA was extracted from formalin-fixed paraffin embedded (FFPE) biopsy samples and patient-matched normal DNA was extracted from mononuclear cells from peripheral blood, respectively, using the QIAamp DNA FFPE Tissue Kit (Qiagen) and the DNA Blood Mini Kit (Qiagen). DNA quality and yield were measured and accessed using a Qubit fluorometer and the Qubit dsDNA HS Assay Kit (Invitrogen). All of the experiments were performed according to the manufacturer's protocol.

WES library generation and sequencing

Whole-exome capture libraries were constructed from 100 ng per sample as input material for the DNA library preparations. Illumina TruSeq™ Sample Prep kits were used for library preparations and NimbleGen SeqCap EZ Human Exome Library was used for exome capture before sequencing. According to the manufacturer's instructions, genomic DNA was fragmented by sonication to a median size of 350 bp. Then, genomic DNA fragments were end-repaired, ligated with the Illumina sequencing adapters, and amplified. Finally, DNA libraries were subjected to WES using the Illumina HiSeq X TEN platform (2 × 150 bp paired-end reads).

Sequencing data analysis

Paired-end reads in FastQ format were quality checked by FastQC and processed to high-quality using Trimmomatic v0.33 [9] to remove adapters and perform trimming. The reads were aligned to the reference human genome (build hg19), using the Burrows-Wheeler Aligner (BWA-MEM) v0.7.15 [10] (parameters -M -t 24). BAM files were processed and sorted from SAM files using Samtools v1.8 and performed deduplication using Picard v1.130 [11]. Genome Analysis Toolkit (GATK) v3.2 (https://www.broadinstitute.org/gatk/) was used to adjust and refine the alignments via Indel Realignment and Base Quality Score Recalibration.

Somatic mutation and copy number alterations calling

Somatic mutations including SNV and small scale insertions/deletions (INDELs) were detected using Varscan2 v2.3.9 in the paired mutation calling mode [12]. The annotation of the somatic mutations was performed with ANNOVAR tool [13]. Mutations with minor allele frequency > 5% in the 1000 Genomes cohort(http://www.internationalgenome.org), the ExAC resource (http://exac.broadinstitute.org), and the NHLBI GO Exome Sequencing Project were removed, but all COSMIC variants were retained. Mutations within the blacklist that we compiled according to the MutSigCV paper were also filtered and removed [14]. Additionally, variants in segmental duplications (UCSC hg19. genomicSuperDups) or repetitive elements (RepeatMasker http://www.repeatmasker.org/) were removed. Finally, all of the mutations after filtering were reviewed manually on the Integrated Genomics Viewer (IGV). FACETS [15] algorithm was utilized to detect somatic allele-specific copy number alterations to determine the tumor purity of the samples. Regions with total copy number = 1 were defined as a loss event, and total copy number = 0 were defined as a homozygous deletion event. Regions with total copy number = 3–8 were defined as a gain event, and total copy number > 8 were defined as an amplification event based on field standard [16].

Cancer cell fraction and clone evolution analysis

We conducted PyClone [17], an algorithm to infer the clonality by the Bayesian clustering method, to infer accurate clonal structures across primary and CRPC tumors. The somatic mutations with depth higher than 40 × were used as input to PyClone. Each mutation input contained the major and minor copy number variations determined from FACETS, and tumor purity was also used to estimate the cancer cell fraction (CCF). The reference and alternate read-depth of each mutation was calculated from the BAM files using Bam-readcount v0.8.0 (https://github.com/genome/bam-readcount) with minimum base quality 20. PyClone beta-binomial model with the “parental_copy_number” and “–tumor_contents” option was run for 50,000 iterations for each case. PyClone clusters were discarded if they contained two or fewer variants. Variant clusters identified by PyClone were imported into ClonEvol [18], a tool for phylogenetic tree inference from CCF clusters, to track and visualize each tumor's clonal evolution. CCF clusters were also visualized in R using the fishplot package [19].

Driver events annotating and mapping onto the clonal evolution trees

To evaluate whether the potential driver events play a significant role in the clone evolution of CRPC, we compiled two driver event lists for somatic mutations and CNA, respectively. Driver gene lists for mutation were selected from genes that overlapped with Network of Cancer Genes (NCG6.0 http://ncg.kcl.ac.uk/), or frequently mutated in CPGEA cohort (Chinese Prostate Cancer Genome and Epigenome Atlas) or highlighted as the significantly mutated genes (SMGs) in prostate cancer representing 13 Western cohorts as described previously [20]. Driver gene mutations were labeled in clonal evolution trees according to the results of ClonEvol. Driver events for copy number alterations were selected from regions identified as recurrently amplified and deleted in a unified “Pan-Cancer” analysis [21], or primary prostate cancer cohorts such as CPGEA [20], TCGA [22], and MSKCC [23], and CRPC cohorts such as SU2C-PCF Dream Team [24]. Driver CNA events were labeled to the clonal clusters by inferring temporal ordering of SNAs and CNAs using Canopy, and by matching the presence/absence status of the copy number events with that of the clusters across samples. Finally, we used Canopy [25], a statistical framework and computational procedure for identifying the subpopulations, to help label the CNA events and verify the inference of the clonal evolution tree.

Results

Patient cohort and clinical information

We selected 5 patients to comprise our small cohort to investigate genomic mutations associated with the transition from HSPC to CRPC following ADT therapy. Clinical information, including treatment history and duration of therapy, and pathological characteristics are summarized in Supplementary Table 1. All tumor samples for WES were collected by prostate biopsy, and sample sources were included in Supplementary Table 2. Briefly, at the initial diagnosis, TURP was performed for patients with urinary symptoms. Multifocal biopsy was performed to harvest the specimens. Samples of high tumor content (>70%, curated by uropathologists and represented the proportion of tumor in the prostate biopsy tissue) were selected for WES. Sites of HSPC tumor biopsy were shown in Figs. 3, 4 and 5. In adherence to the EAU guideline on prostate cancer [7], all patients received the ADT treatment until progressed to CRPC and continued to receive it afterwards (Figs. 1A, 6). The median time between initiation of ADT and CRPC was 141 days (30–224 days) (Supplementary Table 1) (Fig. 1A). Guided by MRI-imaging, one active focus from each patient was biopsied at the CRPC stage by an experienced urologist. After additional quality control, multifocal biopsies from three patients, single biopsies from two patients at the HSPC stage, as well as all five CRPC biopsies underwent WES. The average coverage of WES for tumor and control blood samples ranged from 144 × to 255 × and 78 × to 124 ×, respectively (Supplementary Table 2).

Fig. 3

Clone evolution of case Patient-Y and Patient-L. Clone evolution for the two cases in which only one sample was collected in HSPC. The left panel shows the sites of primary tumor biopsies and the pathological results and Gleason score (Gleason S), and the red dot means biopsy was positive. Fishplots for two cases are shown in the middle panel. The most recent common ancestor (MRCA) clone is marked in gray. The percentage represents the cancer cell fraction (CCF) of the subclone. Give an example to explain fishplot: MRCA in Patient-Y contained somatic alterations TP53, AXIN2, SOX9, 8q gain and 17q gain probably relevant for prostate carcinogenesis. One subclone marked in blue within the MRCA clone evolved to become the dominant clone in CRPC by acquiring additional alterations, including mutations in DNMT3B and Xq12 amplification harbored AR gene. One subclone marked in red disappeared in CRPC due to the clonal selection, which may be sensitive to ADT. The right panel illustrates the phylogenetic relationships by clone evolution tree of the case. The length of branches connecting clones is proportional to the number of mutations that the clone or subclone contained. Identified driver events are marked on the tree and fishplot. CRPC seeding subclones are highlighted in blue.

Fig. 4

Clone evolution from HSPC to CRPC of Patient-C. Clone evolution for Patient-C with multifocal biopsy in HSPC. In the uppermost panel, Bell plots represent the clonal dynamics over time of the individual tumor, and the sphere of cells shows the clonal admixture or subclonal population of individual samples. The percentage range in bell plots indicates the 95% confidence interval of the CCF in the samples. The bottom panels are organized similar to Fig. 3. Fishplots of patients’ multifocal biopsies in HSPC were manually generated based on the bell plots and subclonal population. The clonal evolution tree with branch length scaled by the log2 ratio of the number of variants in each clone. The potential driver variants are highlighted.

Fig. 5

Clone evolution from HSPC to CRPC of Patient-B and Patient-N. Clone evolution for patients collected multifocal biopsies in HSPC. Legend same as in Fig. 4.

Fig. 1

Biopsy, treatment timelines, and mutation landscape of the cohort. (A) Biopsy and treatment timelines for the five patients. We performed whole-exome sequencing on tumor DNA from treatment naiive biopsy samples and CRPC biopsy samples. Green vertical lines and words indicate time points of biopsies, and the red vertical lines indicate the time of CRPC diagnosis. The blue boxes indicate the periods from ADT to the diagnosis of CRPC, with the treatment duration recorded as days. (B), (C), (D), (E) An overview of somatic alterations detected in matched HSPC and CRPC tumors across 5 cases. Each column represents an individual tumor. Panel (A) shows distribution of tumor mutation burden. Each subsequent panel displays the coding SNVs and Indels mutations, and the number of segments for copy number alterations. The bottom panel (E) shows the proportion of clonal and subclonal mutations.

Fig. 6

Clinical course of each patient. Clinical treatment course and PSA responses of each patient. The duration of androgen-deprivation therapy are denoted as colored areas, and the time of sampling as arrows. The top panel of each case showed the number of somatic mutations.

The landscape of somatic mutations and copy number alterations

Changes in the mutational landscape from primary tumor to resistant tumor can provide critical information regarding the mechanisms involved in tumor evolution and provide potential treatment targets. Towards this end, we applied stringent single nucleotide variant (SNV) calling pipelines to our datasets. The average mutation burden was 1.07 per Mb for HSPC samples and 1.02 per Mb for CRPC samples (Fig. 1B). In total, an average of 52 (30–80) SNVs and 5 (3–6) indels were detected in the nine treatment-naïve biopsies. The five matched CRPC prostate biopsies presented an average of 55 (7–145) somatic SNVs and 4 (3–6) somatic indels (Fig. 1C). At the CNA level, we identified an average of 23 (5–35) gains and 4 (0–17) deletions in the nine treatment-naïve biopsies. In the five matched CRPC prostate biopsies, an average of 27 (4–59) gains or amplifications and 8 (0–29) deletions were detected (Fig. 1D). Based on clonal analysis [26,27], the proportion of clonal variants that we were able to assign was 56.58% in the HSPC, whereas in CRPC, we estimated that 67.24% variants were clonal variants (Fig. 1E).

We also selected 10 mutations of Patient-Y for conducting validation experiments using Sanger sequencing (Supplementary Fig. 7). These mutations existed in both HSPC and CRPC samples, and the majority of genes are novel for prostate cancer studies. Finally, eight out of ten mutations have been experimentally validated.

To understand the intratumor heterogeneity of HSPC and how it is reflected in the evolution and establishment of CRPC, we further investigated the non-synonymous somatic mutations of all cases. We observed that tumors from the same patient always shared some common mutated genes, indicating the tumor's clonal origin despite the heterogeneous mutation patterns across foci and time (Fig. 2A). On average we detected 23 (5–36) mutated genes shared between the HSPC and CRPC samples from the same patient, suggesting that the resistant tumors evolved directly from primary tumors. In all but one case, we found that the number of mutated genes decreased from HSPC to CRPC, suggesting that ADT treatment was effective in eliminating most tumor subclones, although new mutations did arise after ADT. To define the driver events for individual patients, we combined Network of Cancer Genes (http://ncg.kcl.ac.uk/) and known prostate cancer genes (frequently mutated in previously analyzed prostate cancer cohorts as potential driver genes [20]. In each patient, 7 (6–8) driver genes were identified on average, including CDK12 in Patient-B, CHD6 in patient C, FANCA gene in patient-N, ZMYM3 in patient L, and TP53 in patient Y (Fig. 2A).

Fig. 2

Driver events and regional distribution of alterations. (A) The regional distribution of nonsynonymous mutations in five patients. The Venn diagram depicts the number of mutated genes, and the heat map indicates the presence of a mutation (color) or its absence (gray) in the individual tumor. Right showed the gene names of mutations and driver genes (red). (B) Overview of somatic copy number alterations identified as driver events across five patients. Alterations in CRPC tumors are indicated in a darker shade and those detected in treatment naive biopsy samples in a lighter shade. Right showed the gene names within driver events. (For interpretation of the references to color in this figure legend, the reader is referred to the web version of this article.)

Furthermore, we discovered novel genes that may be associated with prostate cancer. MYO15A, which is associated with increased genomic instability, was found to harbor a frame shift deletion in Patient-N [28]. CHD6 and LZTR1 were mutated in all samples from Patient-C. CHD6 is an oncogene whose activation is associated with oxidative stress in tumor microenvironments, including in colorectal, uterine, gastric, lung and pancreatic cancers [29]. Non-sense mutation of LZTR1 has been reported to activate the RAS-dependent pathways which lead to dysregulation of growth and differentiation in leukemia [30]. To prove that these three genes may play a role in the progression of castration resistance in prostate cancer, we provide a comparison between the mutation frequency profiles in the primary tumors and CRPC samples using the previously published cohorts (Supplementary Fig. 8). We observed that the mutation frequencies of the three genes in metastatic CRPC (880 samples) and neuroendocrine cancer (54 samples) were significantly higher than those in the primary samples (1400 samples). The CHD6 gene is marked as a significantly mutated gene in the M/DFCI cohort [31], and LZTR1 is a known tumor suppressor gene. Using the RNA expression data from the TCGA Western primary cohort (550 samples) and Chinese primary cohort CPGEA (268 samples), we found that the expression differences of the three genes between tumor tissues and normal tissues were statistically significant in the Chinese population (P < 0.01). In summary, these results provided additional evidence that these new candidate genes may play a role in prostate cancer tumorigenesis and/or metastasis, which also are the direction we need to perform experiments to verify in the future.

We also examined tumor heterogeneity at the copy number level (Supplementary Fig. 3). MYC gene gain was the most frequent clonal copy number variations found in three patients (Fig. 2B). Gain of 17q25.2 which harbored an oncogene, SEPT9, was also detected in two patients (Patient-N and Patient-Y).

CRPC-specific somatic alterations (SNVs and CNAs) were found in 4 patients, but not in Patient-C. Across all five patients, two gained AR amplification in CRPC after ADT. Focal amplifications of AR were restricted to castration-resistant tumors, consistent with the known mechanism of drug resistance in CRPC [18]. In Patient-Y, a subclone with a missense mutation of DNMT3B was only present in the CRPC sample. The mutation pattern of Patient-B suggests that his CRPC sample acquired a “hyper-mutator” phenotype with 33 CRPC-specific mutated genes. Among these, STYK1 was a known cancer driver gene, which encodes a putative kinase and is known to be overexpressed in CRPC [32]. The mutation in STYK1 may induce tumor cell invasion which may contribute to CRPC [32]. Intriguingly, 18 out of 69 mutated genes were typical Cys2/His2-type (C2H2-type) zinc finger (ZF) family genes in Patient-B, and 10 were private to his CRPC sample. Upon manual review of these mutations in IGV, we found that 68% mutations were in the Cys2His2 ZF domains. Consistent with previous studies, the ZF gene mutations can occur in the same individual and influence DNA-binding activity [33]. We integrated the previous cohort and found that CRPC samples were enriched with more ZF gene mutations compared with HSPC (Supplementary Fig. 4). That said, the number of ZF mutations was positively correlated with the total number of mutations in the corresponding tumor samples. These data suggest that ZF family genes might have important roles in the progression of prostate cancer, however, elucidating the mechanism behind these phenomena will require further experimentation.

Clonal evolutionary history in cases of Patient-Y and Patient-L

For Patient-Y and Patient-L, we analyzed a single biopsy each from their treatment-naïve tumor and their CRPC. We applied the PyClone [17] method to analyze the clonal population structure in each tumor, and the ClonEvol [18] method to order and visualize inferred clonal evolution patterns (Fig. 3). In both patients, the predominant subclone after treatment was private to CRPC and harbored driver mutations not detected in the treatment-naïve tumors.

Patient-Y was diagnosed with small cell carcinoma. The primary tumor harbored a clonal TP53 mutation and 8q24 gain encompassing MYC at the initial diagnosis (Fig. 3). We did not detect the deletion of RB1 in Patient-Y (Supplementary Fig. 5), although a previous study reported that concurrent loss of RB1 and mutation or deletion of TP53 was present in 53.3% of castration-resistant neuroendocrine histology [34]. After ADT treatment, a subclone encompassing AR amplification and DNMT3B mutation appeared, whereas a clone characterized by CIC and ITCH mutation was no longer detectable. DNA methyl-transferase enzyme (DNMT3B) is a key enzyme that catalyzes CpG methylation and is frequently mutated in many cancers [35]. Consistent with the recent findings in castration-resistant prostate metastases, somatic mutations of DNMT3B are associated with a novel epigenomic subtype in metastatic CRPC, underscoring DNMT3B’s important role in tumor progression [36].

In Patient-L, the castration resistance was associated with expansion of the subclonal 11q deletion which encompassed FLI1, a prostate cancer gene reported in a recent Chinese primary cohort [20]. The predominant clone at diagnosis had SPOP, NRAS, and ZMYM3 mutations, and experienced sharply decline following treatment, indicating that this clone may be sensitive to ADT.

In Patient-Y and Patient-L, the predominant CRPC subclones were undetectable in primary tumor. However, we cannot rule out the possibility of existence of such resistant subclones in the primary tumor due to intra-tumor heterogeneity. Indeed, different loci of the same prostate cancer may exhibit distinct Gleason scores and tumor clonal composition. To mitigate sampling bias, we conducted multifocal sampling in HSPC in the remaining three patients to more precisely delineate their tumor clonal evolution.

Clonal evolutionary history of cases with multifocal biopsy in HSPC

We profiled two, three, and two biopsies from the HSPC tumors of Patients B, C, and N, respectively (Figs. 4, 5), along with their CRPC tumors. We constructed phylogenetic trees for each multifocal prostate cancer case based on their somatic mutations.

Similar to the findings from Patients Y and L, the CRPC of Patient-B appeared not to be derived from either HSPC biopsies, rather it was seeded directly from the MRCA (most recent common ancestor) which had mutations in known prostate cancer genes, including CDK12 and MYC. His CRPC biopsy obtained a new mutation in STYK1.

In both Patient-C and Patient-N, we found that their CRPC were derived from subclones in the primary tumor after experiencing additional driver events (Figs. 4, 5). In the case of Patient-C, the three primary biopsies had distinct molecular signatures, even though two of them had the same Gleason score (5 + 4). Phylogenetic analysis revealed that the CRPC was most likely derived from subclone PC-1, which was characterized by the loss of STK11, an early event in PC-1 and a potential driver of castration resistance [37]. The predominant subclones from PC-2 and PC-3 were effectively eliminated by ADT treatment.

For Patient-N, PN-1 was a dominant subclone at initial diagnosis and characterized by FANCA mutation and 17q25.2 gain. PN-1 likely gave rise to CRPC, which further acquired a copy gain of AR. FANCA is a tumor suppressor gene involved in DNA damage repair (DDR), and the mutation c. 3482C>T (p.T1161M) has been reported in prostate cancer and recorded in COSMIC (COSM705286). Consistent with the observation that the aberrations in DDR pathway genes were relatively unchanged between hormone-naïve tumor and CRPC, the “FANCA mutation” subclone might not be sensitive to ADT and eventually evolved into CRPC [38].

Driver events in CRPC and clinical course

In three out of the five cases, we found that the dominant subclones after drug treatment were private to CRPC. These newly established subclones likely acquired adaptive driver events during clonal evolution. Evidently, the driver event AR amplification is a biomarker for primary tumor developing into CRPC and is considered a key mechanism responsible for castration resistance and cancer relapse (Figs. 3, 5). In the other two cases, Patient-C and Patient-N, we detected that subclones expanded from the primary tumor to seed CRPC, potentially indicating that these pre-existing subclones may be castration-resistant [39]. Notably, in Patient-N, the driver event 17q25.2 gain in the pre-existing resistant subclones was identified in previous studies such as Pan-cancer [21], and primary prostate cancer cohorts [20] (Supplementary Table 3). The 17q25.2 gain was also detected in the founding MRCA clone from the primary tumor of Patient-Y, suggesting that it may play a significant role in the progression of prostate cancer.

How are clinical features correlated with clonal evolution? We found that PSA (Prostate-Specific Antigen) may not accurately reflect the molecular dynamics of subclonal progression in some patients (Fig. 6). Thus, tracking tumor evolution could provide critical information about tumor progression and drug resistance that standard clinical tests such as PSA failed to provide.

Discussion

The evolutionary history of prostate cancer, especially clonal origin and metastatic spread can be elucidated by analyses of subclonal architecture [40], [41], [42], [43]. Gundem et al. collected multiple metastases from different anatomic sites and clarified the complex pattern of prostate cancer metastasis, such as polyclonal seeding and seeding of tumor clones between metastases [40]. Based on this discovery, Woodcock et al. described the heterogeneity of routes to metastasis in more detail, and found that the heterogeneity associated with metastasis is often acquired within the prostate [44]. However, these studies usually performed spatial sampling from multiple metastases and did not focus on the alterations directly associated with treatment. Previous studies have demonstrated that the clonal evolution constructed by sequential sampling may yield useful information. For example, Haffner et al. reported a case in which the lethal clone arose from a relatively low-grade primary tumor focus, which demonstrated that the Gleason score alone might not be sufficient to determine the lethal clone of prostate cancer [45]. In this study, we aimed to dissect the clonal evolution from HSPC to CRPC by analyzing the mutation landscape of heterogeneous tumor samples from five individual patients. With the exception of one patient, most cases exhibited fewer mutations in CRPC compared to HSPC, suggesting that the pretreatment tumor had a higher degree of genetic diversity [26]. On one hand, ADT treatment was able to eliminate such subclones. On the other hand, this genetic diversity likely empowered tumor cells to undergo the evolutionary bottleneck introduced by ADT treatment, a finding consistent with results from the TRACERx consortium [26]. In addition, the clonal architecture of HSPC and CRPC specimen in individual patients further suggested that they are derived from a common ancestor, supporting the monoclonal origin of lethal CRPC [24,41].

We determined likely driver events underlying the transition from HSPC to CRPC. Gene mutations associated with the transition included well-known prostate cancer genes such as TP53, SPOP and CDK12, as well as novel potential tumor-related genes such as MYO15A, CHD6 and LZTR1. At the CNA level, AR gain and amplification was found in emerging subclones in CRPC in two patients, highlighting AR's driver role in prostate cancer development and specifically in the progression to CRPC [24,46]. We also observed the most frequently altered locus was the clonal gain of chromosome 8q24.13-8q24.3, which validated a previous study in which 8q gain was found to play a crucial role in the pathogenesis and progression of prostate tumor [47]. Interestingly, 19p deletion and 17q25.2 gain persisted in the clones and subclones from HSPC to CRPC and did not get eliminated by treatment, suggesting that mutations associated with these chromosomal abnormalities can drive tumor progression and are insensitive to ADT. We also reported subclones eliminated by the ADT treatment, such as “SPG20-driven” subclone in Patient-N, and “SPOP-driven” subclone in Patient-L. SPOP mutation is an earlier event in prostate tumorigenesis and is much more frequent in non-castration diseases compared with CRPC. Patients with SPOP mutations have a higher response rate to abiraterone, therefore SPOP mutations are indicative of a favorable prognosis [48]. Our results support that tumors with SPOP mutation also respond well to ADT treatment. Importantly, tracking of the dynamics of subclones or mutational events could provide a novel and more effective indicator of ADT treatment response.

Several studies previously investigated the genomic alterations associated with CRPC and resistance to ADT therapy [49]. Grasso et al. performed exome-seq on 50 CRPC patients and 11 treatment-naïve, localized prostate tumors and identified a diverse range of potential driver somatic variations [24]. Some significantly mutated genes are replicated in our relatively small cohort. Similarly, Mateo et al. sequenced 61 patients with matched hormone-naive and metastatic CRPC biopsies, and illustrated the differences in genomic abnormalities between the nonlethal primary prostatic cancers and the metastatic CRPC [38]. However, these studies did not directly investigate tumor evolution in the same patients. In contrast, despite the small cohort size, our study focused on closing the gap in the transition from HSPC to CRPC in the same patient. Our findings are clinically relevant and can offer insights to enable more informed clinical decision-making in the era of precision medicine. For instance, our preliminary findings suggest that the information of the clonal or subclonal changes during treatment could be helpful in monitoring and diagnosis of the progression of CRPC when clinical indicators may not be informative. Indeed, previous studies using phylogenetic circulating tumor DNA (ctDNA) profiling have demonstrated value in investigating the subclonal architectures and delineating genomic evolution of tumors under therapy [50].

Our study had several limitations. First, we were limited by sample size. In addition, we were only able to conduct multifocal profiling for three patients. Thus, our findings need to be replicated or validated in much larger cohorts. Second, our genomic profiling is based on whole-exome sequencing, which has limited power both in identifying cancer driver events and in reconstructing clonal evolutionary history. Prostate cancer is characterized by high genome instability and relatively low mutation load. Thus, clonal analysis based on DNA rearrangement holds the promise to better define the evolutional trajectory from primary to CRPC [51]. Finally, cancer evolution is characterized not only by genetic mutations, but also epigenome and transcriptomic changes. New multiple-omics data and integrated analysis will likely provide more accurate description of the evolutionary history of prostate cancer.

In summary, our work discovered putative driver events underlying the evolution of CRPC. Once validated, these molecular events can serve as biomarkers for disease progression. We highlight that the origin of CRPC and evolution process from primary tumor to relapsed tumor is not directly predictable from the Gleason score of the primary tumor, but can be monitored through molecular profiling and clonal analysis. The dynamic changes in tumor heterogeneity could be exploited for monitoring and management of CRPC in future clinical settings.

Data availability

The WES data in FastQ format are available at the Genome Sequence Archive (GSA) for Human at the BIG Data Center (http://bigd.big.ac.cn/gsa-human/), Beijing Institute of Genomics with the accession number PRJCA004439.

CRediT authorship contribution statement

Wenhui Zhang: Methodology, Writing – original draft. Tao Wang: Software, Visualization. Yan Wang: Writing – review & editing. Feng Zhu: Data curation. Haoqing Shi: Writing – original draft. Jili Zhang: Software. Ziwei Wang: Visualization. Min Qu: Resources. Huaru Zhang: Visualization. Tianyi Wang: Validation. Yuping Qian: Validation. Jinjian Yang: Project administration, Investigation. Xu Gao: Supervision, Resources. Jing Li: Conceptualization, Funding acquisition, Methodology.

Declaration of Competing Interest

The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.