Genomic and transcriptomic alterations during metastasis are considered to affect clinical outcome of colorectal cancers, but detailed clinical implications of metastatic alterations are not fully uncovered. We aimed to investigate the effect of metastatic evolution on in vivo treatment outcome, and identify genomic and transcriptomic alterations associated with drug responsiveness.
We developed and analyzed patient-derived xenograft (PDX) models from 35 patients with colorectal cancer including 5 patients with multiple organ metastases (MOMs). We performed whole-exome, DNA methylation, and RNA sequencing for patient and PDX tumors. With samples from patients with MOMs, we conducted phylogenetic and subclonal analysis and in vivo drug efficacy test on the corresponding PDX models.
Phylogenetic analysis using mutation, expression, and DNA methylation data in patients with MOMs showed that mutational alterations were closely connected with transcriptomic and epigenomic changes during the tumor evolution. Subclonal analysis revealed that initial primary tumors with larger number of subclones exhibited more dynamic changes in subclonal architecture according to metastasis, and loco-regional and distant metastases occurred in a parallel or independent fashion. The PDX models from MOMs demonstrated therapeutic heterogeneity for targeted treatment, due to subclonal acquisition of additional mutations or transcriptomic activation of bypass signaling pathway during tumor evolution.
This study demonstrated in vivo therapeutic heterogeneity of colorectal cancers using PDX models, and suggests that acquired subclonal alterations in mutations or gene expression profiles during tumor metastatic processes can be associated with the development of drug resistance and therapeutic heterogeneity of colorectal cancers.
Colorectal cancers exhibit various genomic and transcriptomic changes during metastasis because of genomic instable nature of cancer cells. However, detailed clinical implications of metastasis-associated alterations have not been fully understood. Herein, we developed and analyzed patient-derived xenograft models from patients with multiple organ metastases. Phylogenetic analysis showed that mutational alterations were closely connected with transcriptomic and epigenomic changes during the tumors' metastatic process. Primary tumors with heterogeneous subclonal architecture showed more dynamic subclonal changes during metastasis, and subclonal acquisition of additional mutations or transcriptomic activation of bypass signaling pathway were responsible for the development of treatment resistance. Especially, gain of ERBB2 L755S mutation and activation of TGFβ signaling pathway were associated with the resistance to an ERBB2 inhibitor, lapatinib and a PI3K inhibitor, BYL719, respectively. Thus, therapeutic heterogeneity from genome and transcriptome dynamics during metastasis should be considered for the treatment of patients with colorectal cancer with metastasis.
Tumor metastasis is a major cause of cancer-related deaths, and, to date, appropriate treatment modalities are very limited in most types of cancers (1, 2). To develop the strategies for preventing or targeting tumor metastasis, understanding the nature and biology of tumor metastasis is inevitable. Genomic instable nature of cancer cells can result in genetic diversity of cells within a given tumor [intratumoral heterogeneity (ITH); ref. 3], and genetic and epigenetic alterations during metastasis process has been suggested to contribute to tumor progression, treatment resistance, and survival outcome (4–7). Recent advances of high-throughput genomic analysis technology provided landscape of genomic alterations specific to metastasized tumors (8), and comparative studies between primary and matched metastasized tumors have described the evolutionary characteristics of metastatic tumors (9–11). However, detailed clinical implications of genetic and epigenetic alterations during metastasis remain limited for most cancer types.
Patient-derived xenografts (PDXs) are cancer models that are established by the transfer of patient tumor tissues into immunodeficient mice. PDXs retain several key characteristics of patients' tumors including histology, genomic alterations, and transcriptomic signatures (12–14). Previous studies have also suggested that PDX models exhibit drug responsiveness comparable to that of the primary tumors in the patients (15–17), and overcome many limitations of conventional in vitro cell line models and cell line xenograft models (13, 18). Most critically, PDX tumors have the capacity to retain the ITH of the original patient tumors and therefore can be applied to the study of clonal selection and tumor evolution (13, 19).
Colorectal cancer is the third most prevalent cancer and the fourth common cause of cancer-related death worldwide (20), and a cancer type reported to exhibit significant ITH (21, 22). In this study, we examined the clonal dynamics and tumor evolution of colorectal cancer during metastasis by analyzing the genome, transcriptome, and epigenome of patient and PDX tumors from patients with colorectal cancer with multiple organ metastases (MOMs), and experimentally demonstrated that acquired genomic and transcriptomic heterogeneities during tumor evolution can be associated with the heterogeneity of drug responsiveness.
Materials and Methods
Colorectal cancer patient sample collection and generation of PDX models
Colorectal cancer tissue samples from primary and metastasis sites and blood samples were obtained from individuals who underwent colectomies at Gil Medical Center from 2014 to 2015. All samples were obtained with informed consent at the Gil Medical Center, and the study was approved by the institutional review board in accordance with the Declaration of Helsinki. For PDX models, surgically resected tissues were minced into pieces approximately ∼2 mm in size and injected subcutaneously in the flanks of 6-week-old NOD/SCID/IL2γ-receptor null (NSG) female mice (The Jackson Laboratory). When a tumor volume reached > 700 to 1,000 mm3, the mouse was sacrificed and tumor tissues were extracted and cryopreserved in liquid nitrogen and stored at −80 °C for generating future PDX passages. Mice were cared for according to institutional guidelines of the Institutional Animal Care and Use Committee (IACUC) of Seoul National University (No. 14-0016-C0A0).
Exome sequencing processing and variant calling
DNA reads were aligned to the merged references for human GRCh19 and mouse mm10 genome versions using BWA-MEM (23). Sorting and marking duplicate reads were performed using Picard tools. Additional processing was conducted to remove mouse sequences using an in-house bioinformatics pipeline, followed by the IndelRealignment, BaseQualityScoreRecalibration, PrintReads commands as Best Practices recommendations of the Genome Analysis Tool Kit (GATK; ref. 24).
After processing the BAM files, Mutect (25) and IndelGenotyper were used for somatic mutation calling, followed by ANNOVAR for functional annotations. Variants with at least 8 read depths, or 4 alternate allele depths and a minimum genotype quality of 20 were maintained. For somatic indels, a strand bias Phred-scale P value of greater than 20, by Fisher exact test, was discarded. Common variants found in dbSNP142 were also removed. Finally, exonic and splicing variants were selected based on the RefSeq database and having population frequencies lower than 0.01, based on The Exome Aggregation Consortium, 1000 Genomes Project, and NHLBI ESP6500.
The pairwise distance of SNPs were calculated based on 4,784 rare nonsilent somatic mutations in tumor tissues from the patient No. 21. The maximum composite likelihood substitution model, which is an accurate alignment-free estimator of the number of substitutions per site based on the lengths of exact matches between pairs of sequences, were applied using MEGA (version 7.0) program (26).
Copy-number alteration analysis using whole exome sequencing data
Copy-number alterations (CNAs) were identified in whole exome sequencing (WES) data based on the read per kilobase per million mapped reads (RPKM) value of exonic regions from CONIFER (27). The logarithm of germline blood and sample data were used for further CNA analysis. Somatic CNAs were segmented and called with the Bioconductor package “DNAcopy.”
To identify CNAs emerging during the generation of PDXs, we compared copy-number profiles between patient and derived PDX tumors. First, a comparative value was obtained for each probe region by subtracting the log2-transformed copy-number value of patient tumor from that of PDX tumor. The circular binary segmentation analysis was then applied to relative copy-number data with the “DNAcopy” R package. If |segment value| > 0.3 and size > 5 Mb, segment calls were retained. Second, segment calls from relative data and each sample were summarized to gene-level calls. A discordant CNA call between the patient and PDX tumor was considered an actual PDX-specific CNA, and “false” discordant CNA calls, which is called in both samples but showed a high relative value because of the high purity of the PDX tumor, were eliminated. The number of discordant CNA calls was divided by the total number of genes (excluding genes with a neutral copy number call in both datasets).
Calculation of heterogeneity score and analysis of clonal architecture
The subclonal genomic structure of the patient and PDX tumors was inferred using mutant allele fractions and copy-number data from WES. Based on the Bayesian clustering method, SciClone (28) and PyClone (29) estimated the clonal population structure of each sample. In particular, PyClone was used for the estimation of cellular prevalence of each cluster of each sample with respect to normal cell contamination and copy-number changes. PyClone was run for 10,000 iterations with a burn in period of 1,000 iterations using a beta binomial parameter of 1,000. We plotted clusters with at least 5 mutations and selected the representative genes for each cluster if mutated genes were included in our cancer gene list or showed high mutation frequencies in cBioPortal colorectal cancer studies (http://www.cbioportal.org). Of the primary patient–PDX tumor pairs, the number of input mutations in 2 patients (patient No. 1 and No. 2) was too small or large to obtain the results in PyClone analysis.
RNA-seq data processing and expression analysis
RNA-seq reads from each whole transcriptome sequencing (WTS) experiments were aligned to the same reference as for WES using STAR aligner (30). The following procedures were performed as the Best Practices workflow for RNA-seq using GATK. Gene expression levels were quantified using deduplicated BAM files by fragments per kilobase of exon per million mapped reads (FPKM) using HTSeq-count (31) based on the Homo sapiens GRCh37 Ensemble v65. FPKM values were normalized and log transformed with edgeR. We also used pairwise average-linkage clustering for Pearson distance measurement of common clustering patterns. The DESeq2 algorithm (32) was used to detect metastatic tumor-specific gene expression patterns when compared with the primary tumor of each patient.
Methyl-capture sequencing was performed using the SureSelectXT Human Methyl-Seq Target Enrichment System, which utilized a capture-then-bisulfite-convert approach, as directed by the manufacturer. Genomic DNA (3 μg) of PDX tumors was sheared to a median fragment size of 150 to 200 bp using a Corvaris E210. Libraries were generated and hybridized using SureSelectXT Methyl-Seq Kit (Agilent). Bisulfite conversion was performed using EZ DNA Methylation Gold Kit (Zymo Research). The enriched and bisulfite-converted libraries were indexed by PCR amplification as directed in the protocol and purified using AMPure XP beads (Beckman Coulter). After the quality and quantity of the library sample was assessed by Qubits dsDNA HS Assay Kit (Thermo Fisher Scientific) and the Bioanalyzer DNA 1000 Kit (Agilent), the libraries were sequenced with HiSeq instruments.
The raw sequence data were trimmed for potential adapter sequences using trimgalore and aligned to combined reference using Bismark program with bowtie1. Deduplication and extraction of methylation was also performed with Bismark. Because the capture-then-bisulfite-convert approach captures the top strand of the DNA, only reads that aligned to the original top strand were considered for calling cytosine methylation. Only CpGs with at least 7× coverage were kept for comparisons. The methylation ratio was the number of methylated Cs divided by the total of methylated and unmethylated Cs. Average methylation level of each CpG island and CpG shores were calculated annotated against RefSeq gene database. For 13,199 genes, expression of each gene was correlated with multiple methylation probes per gene and the region that is most anticorrelated with expression data was selected for a given gene.
Phylogenetic analysis using sequencing data
Phylogenetic trees were constructed separately using (i) point mutations, (ii) gene expression, and (iii) methylation from each patient and corresponding PDX models with multiple metastatic tumors using the neighbor joining method of the MEGA7 program (26). Rare functional mutations were used for WES phylogeny tree analysis.
Generation of cancer-related gene list
Single sample gene set enrichment analysis projection on a collection of hallmark gene sets
An unsupervised gene enrichment method that calculates a separate enrichment score (ES) to gene set for each individual sample, independent of sample phenotype labeling, single sample gene set enrichment analysis (ssGSEA), was performed on expression and methylation data. We scored signatures of hallmark processes of the Molecular Signatures Database (MSigDB) through the ssGSEA Projection module of GenePattern, ES was obtained and further analyzed by transforming into a z-score by subtracting the mean of the ESs assigned to all other gene-sets and by dividing the result to their SD. A positive ES denotes a significant overlap of the signature gene set with groups of genes at the top of the ranked list, whereas a negative ES denotes a significant overlap of the signature gene set with groups of genes at the bottom of the ranked list.
Histological analysis and IHC
After resection, tumor tissue samples were fixed in 10% neutral buffered formalin and embedded in paraffin. Hematoxylin-Eosin staining (H&E) was performed according to standard protocols. For IHC, tumor tissue samples obtained from patients and PDX mice were stained using antibodies including CEA (Cell Marque), CK7 (Dako), CK20 (Dako), Ki-67 (Dako), E-cadherin (Zymed), Vimentin (Zymed), CD3 (Dako), and CD31 (Cell Marque). Formalin fixed and paraffin embedded tissue sections of tumor samples were immunohistochemically stained for expression of those antibodies using a BenchMArk ULTRA automatic immunostaining device (Ventana Medical Systems) with the OptiView DAB IHC Detection Kit (Ventana Medical Systems) according to the manufacturer's instructions.
In vivo pharmacologic studies
For PDX mice, drug treatments began after tumors reached approximately 200 mm3. For standard treatments, mice were randomly divided into 2 treatment groups, consisting of 5 mice in each group: (i) vehicle only, (ii) 5-FU (ApexBio, 50 mg/kg, weekly) + oxaliplatin (ApexBio, 5 mg/kg, weekly). For targeted treatments, mice were randomly divided into 4 treatment groups consisting of 5 mice in each group: (i) vehicle only, (ii) lapatinib (ChemieTek, 30 mg/kg, twice a day), (iii) trametinib (ApexBio, 2 mg/kg, daily), and (iv) BYL719 (ChemieTek, 25 mg/kg, daily). The vehicle for 5-FU and oxaliplatin was 5% (w/v) glucose in water (Sigma). The vehicle for lapatinib was 0.5% (v/v) methylcellulose (Sigma) and 0.5% (v/v) Tween 80 (Sigma) in PBS. The vehicle for trametinib was 0.5% (v/v) hydroxypropyl methylcellulose (HPMC; Sigma) and 0.2% (v/v) Tween 80 in PBS. The vehicle for BYL719 was 10% (v/v) ethanol (Merck), 30% (v/v) polyethylene glycol (Sigma), and 60% (v/v) Phosal 50 PG (Lipoid). 5-FU and oxaliplatin was administered via intraperitoneal injection whereas lapatinib, trametinib, and BYL719 were all administered via oral injection for 21 to 26 days.
The data in this study have been submitted to the National Center for Biotechnology Information (NCBI) BioProject (http://www.ncbi.nlm.nih.gov/bioproject) under accession number PRJNA400542.
Generation of genomically defined PDXs from patients with colorectal cancer
We generated a PDX cohort of colorectal cancers and analyzed the mutational profiles of the original tumors and the derived PDX models. We implanted a total of 95 tumor tissues from 42 patients with colorectal cancer (detailed clinical information in Supplementary Table S1), and successfully established 72 serially-transplantable PDXs from 35 patients (Fig. 1A; Supplementary Table S2). Anatomical location of primary tumors, tumor stage, microsatellite instability (MSI) status, and organ sites of metastatic tumors were not associated with PDX engraftment rates (Supplementary Table S3).
Given that a genomically well-defined PDX cohort is a valuable tool for assessing drug efficacy and understanding drug resistance mechanisms, we analyzed the genomic profiles of the original patient and PDX tumors using WES. We sequenced 76 of the original tumor samples from 29 patients (29 primary and 47 metastatic tumors; Fig. 1A; Supplementary Table S2), and 55 PDX samples from 24 patients (21 primary and 34 metastatic tumors; Fig. 1A; Supplementary Table S2). The most frequently mutated cancer genes in our PDX tumors paralleled those seen in colorectal cancers in The Cancer Genome Atlas (TCGA) database (Fig. 1B; ref. 35); with the exception of an increased frequency of PIK3CA mutations in our PDX tumors. These suggest that overall, our colorectal cancer tumors and derived PDX models retain mutation profiles similar to that of other colorectal cancer cohorts. In addition, we performed RNA sequencing for 40 PDX samples from 10 patients (9 primary and 31 metastatic tumors) and DNA methyl-capture sequencing for 18 PDX samples from 3 patients with MOMs (3 primary and 15 metastatic tumors; Fig. 1A; Supplementary Table S2).
In mutational analysis, the allele frequencies (AF) of mutations were correlated between these patient and PDX tumors (average Pearson correlation coefficient = 0.60), but AF usually increased in PDX tumors consistent with the notion that tumor cells are enriched in PDX samples (Fig. 1C; Supplementary Fig. S1A). Somatic CNAs (SCNA) profiling using WES data showed that SCNAs were also generally correlated between patient tumors and corresponding PDX tumors (average Pearson correlation coefficient = 0.57; Fig. 1C; Supplementary Fig. S1B). We also estimated the clonal architecture using PyClone (29) and found that each tumor had 2 to 6 major subclones, which were well conserved and enriched in PDX tumors (average changes in clonal prevalence: 0.286; Fig. 1C; Supplementary Fig. S1C; Supplementary Table S4).
In histologic analyses, the PDX tumors retained comparable histologic architecture of the primary patient tumors, and cancer cells exhibited similar expression patterns with respect to carcinoembryonic antigen (CEA), cytokeratin 7 (CK7), CK20, Ki-67, and E-cadherin staining (Fig. 1C; Supplementary Fig. S2). When evaluating the tumor microenvironment, stromal vimentin expressions of the primary patient tumors were retained in half of the derived PDX tumors (Fig. 1C; Supplementary Fig. S2), and vessels of human origin (CD31-positive cells) were replaced with mouse vessels (Fig. 1C; Supplementary Fig. S2). In addition, CD3+ lymphocytes were enriched in PDX tissues, suggesting that human lymphocytes are expanding in immunocompromised states of early PDX passages (Supplementary Fig. S2). Taken together, our observations showed that our PDX models retain genomic and histologic characteristics of the primary patient tumors, but lose some characteristics associated with tumor microenvironment.
Analysis of genomic alterations during tumor metastasis
In our colorectal cancer cohort, we have 5 patients with MOMs, and successfully generated PDX models from the primary tumors of these 5 patients and 25 of their metastases (patient No. 10, No. 14, No. 21, No. 34, and No. 35 in Supplementary Table S2). In data analysis from WES on the original patient samples and WES/RNA-seq/methyl-capture sequencing of the derived PDX samples (unfortunately the RNA quality from the patient samples precluded RNA-seq), phylogenetic analyses revealed branched evolution for the primary tumors and derived PDXs for all 5 patients (Fig. 2), consistent with the notion that the evolutionary processes were well conserved in the PDX models. Furthermore, the evolutionary processes, determined by somatic mutations, mRNA expressions, and DNA methylation, exhibited similar patterns in the primary patient tumors and the derived PDX samples (Fig. 2), indicating that genomic changes during the cancer evolution process were closely connected with transcriptomic and epigenomic changes.
In tumors with MOMs, the somatic mutations were divided into 3 categories: truncal (mutations present in all tumor regions), branched (mutations present in at least 2, but not all regions), and private mutations (mutations present in 1 region). In the 5 patients analyzed, the average rates of truncal, branched, and private mutations were 30.5% (range, 16.1–39.1%), 16.3% (range, 9.6–25.9%), and 53.2% (range, 36.8–63.7%), respectively (Supplementary Table S5). Of the 608 cancer-related genes with known importance (the selection criteria are described in Materials and Methods), mutations in KRAS (3/5), PIK3CA (3/5), TP53 (2/5), and APC (2/5) were detected as truncal mutations in 2 or more patients (Fig. 2), suggesting that mutations of these genes play a critical role in tumor initiation and early development of colorectal cancers.
Comparison of genetic changes between metastasis and PDX generation
PDX models provide unique tumor microenvironment different to human, raising the possibility that increased genomic instability occurs. To estimate the de novo genetic changes during PDX generation in mice, we compared the changes of mutations and CNAs during PDX generations with those during tumor evolution. In patient No. 21, primary tumor (T75) was metastasized to liver (T74, T79, T80, and T91) and regional lymph nodes (T81, T82, and T191), and we generated 4 different PDX models by engrafting T75 (Fig. 3A). Analysis of WES data showed that estimated differences between primary tumor and PDXs (range, 0.07–0.12, mean: 0.09) were lower than those between primary and metastasized tumors, especially tumors metastasized to liver (range, 0.20–0.23, mean: 0.22; Fig. 3B and C). These data suggest that genetic alterations during the PDX generation were not significantly augmented compared with the metastasis within the patient.
Clonal architecture dynamics during tumor metastasis
Next, we analyzed the clonal architectures of the PDX tumors from 5 patients with MOMs, because clonal architecture of patient tumors was conserved in PDX tumors (Fig. 1C; Supplementary Fig. S1C). The SciClone analysis, modeling both clonal and subclonal mutation clusters (28), revealed that the mutations in the tumor samples from each individual could be organized into 2 to 8 clusters (Fig. 4). All metastasized tumors showed multiple subclones, and these subclones were usually found in primary tumors (Fig. 4), suggesting that metastasis occurs via polyclonal seeding, not a single cell. Patients having primary tumors with a small number of subclones (2 to 3 subclones; patient No. 10 and No. 35) showed little change in subclonal architecture during metastasis (Fig. 4). However, patients having primary tumors with a large number of subclones (5 to 7 subclone; patient No. 14, No. 21, and No. 34) exhibited dynamic changes in subclonal architecture during metastasis (Fig. 4), and these changes were largely dependent on the organs where the cancer had metastasized (Fig. 4). Tumors metastasized to lymph nodes showed similar subclonal architecture with primary tumors, but tumors metastasized to ovary or liver showed enrichment of specific subclones showing little ratios in primary tumors (Fig. 4). For example, C4 subclone in patient No. 14 was enriched in tumors metastasized to ovary, and C3 suclone in patient No. 21 and C4 subclone in patient No. 34 were significantly enriched in tumors metastasized to liver (Fig. 4), suggesting that loco-regional metastasis to lymph nodes and distant metastasis to liver or ovary probably occurs in a parallel or independent fashion and distant organs provide unique tumor microenvironment for distinct clonal selection.
Transcriptomic and epigenomic changes during tumor metastasis
We also analyzed the RNA-seq data of PDX tumors from patients with MOMs. Clustering analysis and principal component analysis (PCA) showed that samples from the same patient were more adjacently clustered compared with metastases to the same organ (Supplementary Fig. S3). When we investigated enriched hallmark gene sets (36) among the metastasized tumors using single sample gene set enrichment analysis (ssGSEA), gene sets of “UV_response,” “DNA_repair,” “Apoptosis,” “Myc_ targets,” and “Protein secretion” were enriched in metastasized tumors in 4 or more patients (Fig. 5). Interestingly, gene sets of “Epithelial_mesenchymal_transition” and “Hypoxia,” which were known to promote metastasis (37, 38), were depleted in metastasized tumors in 5 or more patients (Fig. 5). It is probable that metastasized tumors are exposed to more oxygen compared with primary tumor sites, but detailed investigation was needed in further studies.
We also performed ssGSEA using DNA methyl-capture sequencing data and investigated hallmark gene sets enriched in metastasized tumors. However, there is little gene set enriched or depleted in metastasized tumors of all tested 3 patients with MOMs (Fig. 5). The promoter methylation in gene set of “Cholesterol_homeostasis” was enriched in metastasized tumors of all 3 patients, and the mRNA expression levels of this gene set were inversely correlated with methylation status (Fig. 5). Although cholesterol was reported to promote metastasis (39), gene set of “Cholesterol_homeostasis” was downregulated in metastasized tumors, and detailed molecular mechanisms need to be investigated.
Tumor heterogeneity analysis of patient with colorectal cancer with MSI-H status
Among the 5 patients with MOMs, we focused on 1 patient exhibiting the most divergent mutational profiles due to high levels of microsatellite instability (patient No. 21). This case showed prominent differences in mutation and expression profiles between 2 distinct groups of tumors; primary/LN group and liver-metastasized group (Fig. 4; Supplementary Fig. S4).
In subclone analysis, the C1 subclone was found in all tumors (Fig. 4), and possesses several truncal mutations including CTNNB1, ERBB3, NF1, BRCA2, SYK mutations (Fig. 4; Supplementary Fig. S4A). In addition, copy neutral loss of heterozygosity (cnLOH) in chromosome 3p, which contains CTNNB1 gene, were also found as truncal alterations (Supplementary Fig. S4B). The differences between primary/LN group and liver-metastasized groups were predominantly based on the different prevalence of 2 subclones, C2 and C3 (Fig. 4). In the primary/LN group, the most enriched subclone was a C2 subclone, which included mutations such as TSC2 and MSH6 mutations (Fig. 4; Supplementary Fig. S4A). Interestingly, convergent bi-allelic inactivation of TP53 gene was observed in metastatic tumors in LN group (Supplementary Fig. S5). However, the liver-metastasized group showed enrichment of the C3 subclone, which harbored the TP53 R273C mutation (Fig. 4; Supplementary Fig. S4A). Comparison of copy-number (CN) values for chromosome 7 (chr7) between the 2 groups showed that primary/LN group showed significant increase in CN alterations on chr7 (Supplementary Fig. S4B; P = 0.012). Notably, copy-number changes for EGFR showed significant differences between the primary/LN group compared with the primary/liver group (P = 0.004), and the pathway-level expression analysis also showed that EGFR signaling in cancer was upregulated in the primary/LN group (Supplementary Fig. S4C). In addition, expression of genes in “Signaling by Wnt” increased in primary/liver group (Supplementary Fig. S4C), which was compatible with decrease of promoter methylation in “Wnt_beta_catenin_signaling” pathway (Supplementary Fig. S4D).
For the X201 samples (which were derived from recurrence in the liver after adjuvant chemotherapy with FOLFIRI and cetuximab), we identified a subclone that was not evident in the pre-chemotherapy samples (C6; Fig. 4). In these treatment-resistant subclones, potential driver mutations included mutations of EGFR S464L, was previously reported to be related to the cetuximab resistance (40).
Association of genomic and transcriptomic signatures with therapeutic heterogeneity as delineated by PDX models
Next, we investigated the responsiveness of targeted treatment based on genomic profiling in the PDX models from patient No. 21 (Supplementary Table S6). The ERBB3 G284R mutation was reported to activate the ERBB2 signaling pathway, which can be targeted by the ERBB2 inhibitor, lapatinib (41). NF1 mutations result in activation of the Ras-MEK-ERK signaling pathway, which can be targeted by a MEK inhibitor, trametinib (42), and PIK3CA E545A mutations, found in X80, can be targeted by PI3K inhibitors such as BYL719 (Supplementary Table S6; refs. 43, 44). When the normal growth curves of these PDX tumors were examined, the tumor with a PIK3CA E545A activation mutation (X80) showed the fastest growth (Supplementary Fig. S6).
The combination treatment of 5-florouricil and oxaliplatin (5-FU/Oxa), which is a standard nontargeted treatment for colorectal cancer, showed a significant tumor inhibitory effect in 3 out of the 5 samples tested (Fig. 6A; Supplementary Fig. S7A). Lapatinib treatment resulted in significant tumor inhibition in 3 samples (X75, X79, X80; Fig. 6B; Supplementary Fig. S7B). However, 1 sample with ERBB2 L755S mutation (X91) was refractory to lapatinib treatment (Fig. 6B; Supplementary Fig. S7B), which was explained by the presence of the lapatinib-resistant mutation in ERBB2 gene (ERBB2 L755S; Supplementary Table S6; ref. 45). We validated this mutation with droplet digital PCR (ddPCR), and found the ERBB2 L755S mutation to be a private mutation (i.e., found in only 1 sample; Supplementary Fig. S8). Trametinib treatments showed effective inhibition on tumor growth in 4 sample out of 5 samples (Fig. 6C; Supplementary Fig. S7B). Unexpectedly, BYL719 was effective only in the primary tumor (Fig. 6D; Supplementary Fig. S7B), which was not anticipated by its PIK3CA mutation status.
To understand the factor(s) associated with BYL719 responsiveness, we compared genomic and transcriptomic differences between responsive and nonresponsive tumors. In mutational profiling, mutations in TP53 (R273C and A159V) were detected only in metastatic tumors (Supplementary Table S6). However, overexpression of mutant TP53 did not appear to significantly contribute to resistance to BYL719, as compared with the wild-type (Supplementary Fig. S9). From expression profiling, the TGFβ signaling pathway was enriched in the metastatic tumors (Supplementary Fig. S10), and treatment of TGFβ1 increased the IC50 for BYL719 (Fig. 6E), suggesting that activation of TGFβ signaling pathway is associated with resistance to BYL719. Taken together, in vivo experiments using our PDX models verified that subtle differences in the genome and transcriptome could largely affect the therapeutic responsiveness.
During metastasis, genetic and epigenetic alterations are generated due to the genomic instability of cancer cells, and are dynamically reshaped via evolutionary process and adaptation to environmental changes. The tumor evolution process is generally delineated as a branched evolution model, based on phylogenetic reconstructive analyses (10, 11, 46–48). The branched evolution model emphasizes the importance of targeting truncal alterations. However, resistance for targeted drug may also occur due to subclonal alterations, at the level of gene mutations or gene expression as shown in our data (Fig. 6).
It has been ambiguous whether the subclonal alterations associated with drug resistance exist as minor subclones in the original tumor or emerge during the evolutionary process. The Big Bang model of tumor growth suggests that a single expansion produces diverse subclones where clonal and subclonal alterations can develop early in cancer growth (49). A study using a complex DNA barcode labeling system showed that, during the development of drug resistance caused by ITH, the majority of resistant clones were part of preexisting subpopulations (50). Our analysis of clonal architecture of metastatic tumors also demonstrated that most subclones in metastatic tumors were found in primary tumors (Fig. 4), suggesting that genomic divergence during metastasis occurs usually by evolutionary reshaping of preexisting subclones. However, our ddPCR results also supports the notion that some drug-resistant alterations could appear in the middle of evolutionary process from de novo alterations (Supplementary Fig. S8), and the branched or private alterations play critical roles in the development of drug resistance. In our study, 1 metastatic tumor with the ERBB2 L755S mutation showed the resistance to lapatinib, and transcriptional activation of the TGFβ signaling pathway increases the drug resistance to BYL719 (Fig. 6). Therefore, for the treatment of cancer patients with MOMs, genomic and transcriptomic heterogeneities should be considered for pertinent treatment, especially in patients with colorectal cancer having primary tumors with a large number of subclones, because subclonal architecture changes dynamically in these patients.
Our analyses provide several applicable points for the targeted treatment of patients with CRC with MOMs. Considering the heterogeneous nature of primary and metastatic tumors, targeting the truncal alterations would be the best way for the treatment of patients with MOMs, because all primary and metastatic tumors possess the truncal alterations. In our study, targeting NF1 mutations, which were truncal alterations, with trametinib reduced tumor growth in all tested samples (Fig. 6C). However, at the same time, the incidence of resistance due to subclonal alterations needs to be carefully examined. Because multiple biopsies are usually difficult to perform, 1 alternative way is to analyze the genomic profiles of circulating tumor cells or cell-free DNA, which represent the summation of genomic alterations from all cancer cells. Detection of resistance-related alterations in liquid biopsy can guide to treat the patients with other targeted drugs.
Detailed analysis of subclonal architecture of tumors from patients with colorectal cancer with MOMs suggested some implications about mode of tumor metastasis. Several subclones were found simultaneously in primary and metastatic tumors (Fig. 4), which support polyclonal metastasis rather than metastasis via a single cell dissemination. And, initial number of subclones in primary tumors is associated with diversity in subclonal architecture according to metastasis, and primary tumors harboring more subclones showed more diverse subclonal architecture among metastatic tumors, probably due to higher chance of clonal selection in metastasized tumor microenvironment. In addition, loco-regional metastasis to lymph nodes and distant metastasis to liver or ovary exhibited distinctive subclonal architectures, suggesting that loco-regional and distant metastasis probably occurs in parallel or independent ways.
PDX models retain several characteristics of patient tumors and showed drug responsiveness compatible with the patient's original tumors (15–17). Our comparison of the patient's primary tumors with the derived PDX tumors indicated that the derived PDXs appear to maintain the characteristics of the patient original tumor on the basis of mutational profile and signatures, SCNAs, tumor clonal architectures and histology (Fig. 1). PDX models also retained the characteristics related with evolutionary processes in terms of somatic mutation, gene expression, and DNA methylation (Fig. 2). In addition, genetic changes during PDX generation were similar or lower than genetic changes during metastasis processes in patients (Fig. 3). However, PDX models showed some changes in its tumor microenvironment, such as a decrease in extracellular proteins and vessels of human origin (Supplementary Fig. S2), which was probably due to the replacement of these cells with mouse stromal cells. In addition, analysis of the mutation AFs and clonal architecture in the PDX models exhibited an enrichment of cancer cells having driver mutations or drug responsiveness-associated mutations, compared with the patient samples. For example, the prevalence of PIK3CA mutations in our PDX cohort was higher than our patient cohort (Fig. 1B), which suggests that the subclones with PIK3CA mutations (preexist or de novo) had growth advantages in the PDX environment (Supplementary Fig. S6). Therefore, PDX models show accordant prediction of drug responsiveness in patients by augmenting the portion of subclones with driver mutations or drug responsiveness-determining mutations.
Our study demonstrated that development of genomic and transcriptomic alterations during the metastatic process plays critical roles in determining drug responsiveness, and, as a result, are highly associated with clinical outcome of patients. Furthermore, PDX models are valuable tools for studying the correlation between tumor evolution and therapeutic heterogeneity. More studies of tumor heterogeneity during metastasis with respect to therapeutic effectiveness would pave the ways for prediction and improvement of clinical outcome in patients with colorectal cancer.
Disclosure of Potential Conflicts of Interest
C. Lee is a consultant/advisory board member for Nabsys Inc., Genome's Company, Bionano Genomics, and Seven Bridges. No potential conflicts of interest were disclosed by the other authors.
Conception and design: S.-Y. Cho, C. Lee, W.-S. Lee, H. Park, J.-I. Kim
Development of methodology: S.-Y. Cho, H. Park, J.-I. Kim
Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S.-Y. Cho, J. Chae, D. Na, W. Kang, A. Lee, S. Min, J. Kang, J. Lee, C.O. Sung, H. Park
Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S.-Y. Cho, J. Chae, D. Na, W. Kang, C.O. Sung, J.H. Chuang, J.-I. Kim
Writing, review, and/or revision of the manuscript: S.-Y. Cho, J. Chae, C. Lee, J.-I. Kim
Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): B. Choi
Study supervision: S.-Y. Cho, C. Lee, W.-S. Lee, J.-I. Kim
This work was supported by Gil Hospital Research Grant (Grant No. GCU-2015-5092); the Korean Healthcare Technology R&D project through the Korean Health Industry Development Institute (KHIDI) funded by the Ministry of Health & Welfare, Republic of Korea (Grant No. HI13C2148); the International Research & Development Program of the National Research Foundation of Korea (NRF) funded by the Ministry of Science, ICT & Future Planning (Grant No. 2015K1A4A3047851); the National Research Foundation of Korea (NRF) grant funded by the Korea government (MSIP; Grant No. 2016R1E1A1A01943469 and 2017R1C1B2002183); the Basic Science Research Program through the National Research Foundation of Korea (NRF) funded by the Ministry of Education (Grant No. 2017R1D1A1B03035453); and the National Cancer Institute of the NIH (Grant No. P30CA034196). C. Lee is a distinguished Ewha Womans University Professor supported in part by the Ewha Womans University Research grant of 2018.
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.