Patient-Derived Organoids from Multiple Colorectal Cancer Liver Metastases Reveal Moderate Intra-patient Pharmacotranscriptomic Heterogeneity

Drug sensitivity and gene expression profiling of PDOs from resected CRLMs

An overview of the ex vivo pharmacogenomic platform for analysis of PDOs derived from resected CRLMs is shown in Fig. 1A. Tumor cells were isolated from multiple CRLMs (n = 75 lesions) from twenty-nine patients and cultured ex vivo. PDOs were successfully established for 39 lesions (lesion-wise success rate 52%) from twenty-two patients [patient-wise success rate 76%; 95% binomial proportion confidence interval (CI): 58%–88%; Supplementary Table S1], including five distinct metastatic lesions from one patient, three lesions from four patients, two lesions from another five patients, and one lesion from twelve patients (Fig. 1B). There was a moderate correlation between the number of lesions attempted to culture per patient and the numbers of successfully established PDOs (r = 0.39, P = 0.035, Pearson's correlation test, Fig. 1B), indicating that successful culture was not only dependent on the number of sampled tumors per patient, but also on the biological traits of the sampled cancer cells. Most PDOs had an epithelial architecture and presented with a lumen enclosed by a single epithelial cell layer (Fig. 1C left and Supplementary Fig. S1). Epithelial colorectal cancer tissue origin was assessed by immunostaining of the proteins CDX2, CK7, and CK20 (Fig. 1C right and Supplementary Fig. S2). The PDOs exhibited large differences in growth rates and the time from tumor resection to completed drug screening varied from 3.7 to 25 weeks (median = 9.2 weeks).

• Download figure
• Open in new tab
• Download powerpoint

Figure 1.

A living biobank of colorectal cancer liver metastases (CRLM). A, Included patients and study overview. B, PDO growth success rate and the number of lesions. C, PDO growth morphologies and epithelial markers. CDX2, caudal type homeobox 2 transcription factor; CK7, cytokeratin 7; CK20, cytokeratin 20. D, Representativeness of the PDO cohort. Unsupervised principal component analysis of gene expression data highlighting lesions from patients with successful PDO cultures (black, n = 31 tumors from 14 patients) in the context of the in-house consecutive cohort of surgically treated patients with colorectal liver metastases (n = 290 tumors and n = 174 patients). E, Principal component analysis plot showing correspondence between gene expression from individual PDOs and their respective tumors. Multiple-lesion CRLMs and corresponding PDOs from five patients are indicated with colors and numbers inside the plot, whereas single-lesion CRLMs and corresponding PDOs are shown in the legend. Distributions of Spearman's correlation coefficients are from matched and non-matched PDOs and tumors. For twenty-one PDOs analyzed, corresponding tumor tissue was not analyzed. p, patient; PC, principal component.

Representativeness of the PDOs for the resected CRLMs was evaluated by gene expression profiling. Principal component analysis (PCA) of a subset of the corresponding tissue samples showed that the tumors grown as PDOs were representative of a large in-house consecutive cohort of resected CRLMs (n = 290 tumors; 174 patients, 1–15 tumors per patient, unpublished data; Fig. 1D). Furthermore, gene expression profiles of these PDOs and their original tumors showed an overall fair clustering (Fig. 1E). The mean Spearman's correlation coefficient between matched PDOs and tumors was significantly higher than the mean between unmatched PDOs and tumors (difference in mean Spearman's correlation = 0.27, P = 1.3 × 10⁻⁹). Principal component (PC) 2 of the tumor samples (but not PC1) correlated with a gene expression–based enrichment score for liver-specific genes (Supplementary Fig. S3), showing that infiltration of non-malignant cells in the tumor microenvironment has an important impact on gene expression signals in the tissue samples, and this is a likely explanation for deviations between some of the tumors and the corresponding PDO.

All PDOs were evaluated for sensitivities to thirty-eight clinically relevant drugs (Supplementary Table S2) and sensitivities were calculated as DSS values (34), representing the area over the dose–response curve across nine different concentrations of each drug. Replicate screens for each PDO showed a high technical reproducibility, with a median Pearson's correlation coefficient across samples of 0.98 (range 0.87–0.99, Supplementary Fig. S4A and S4B). Notably, the reproducibility varied among the different drugs and two drugs (ipatasertib and linsitinib) had a particularly high variability between technical replicates (SD = 5.5 and 6.4, respectively; Supplementary Fig. S4C). Replicate screens of ten PDOs performed from 1 to 17 months apart and following a freeze/thaw cycle showed strong biological reproducibility, with Pearson's correlations ranging from 0.81 to 0.93 (Supplementary Fig. S4D). Also, drugs with a similar mechanism of action clustered together in hierarchical clustering based on Pearson's correlations (Supplementary Fig. S5). Comparisons of the three drugs targeting EGFR (afatinib, cetuximab, and erlotinib) showed Pearson's correlation coefficients ranging from 0.86 to 0.94, and the two drugs targeting MEK (binimetinib and trametinib) had a correlation of 0.80. Resistance to EGFR inhibition was confirmed in PDOs with RAS mutations (KRAS mut = 14, NRAS mut = 3, RAS wt = 22; Fig. 2A; ref. 41), and the median DSS for cetuximab in RAS-mutated PDOs was 5.0 (range, 0–19), compared with 21 (range, 5.2–29) for wild-type PDOs (P = 0.002, Wilcoxon rank-sum test). Pharmacological associations were also found for a few targeted drugs with different mechanisms of action, such as between idasanutlin (MDM2/TP53i) and palbociclib (CDKi; Pearson's correlation 0.61) or alisertib (Aurorai; Pearson's correlation 0.51; Supplementary Fig. S5).

• Download figure
• Open in new tab
• Download powerpoint

Figure 2.

Pharmacological patterns of PDOs from colorectal cancer liver metastases. A, Pharmacogenomic relationship between RAS mutation status and response to EGFR inhibition and MEK inhibition. B, Relative drug sensitivities across the PDO cohort. Heatmap visualizes drug-wise median-centered DSS values, with higher (red) and lower (blue) values, indicating higher and lower sensitivity, respectively. Samples were ordered by hierarchical clustering using Pearson's correlation distance and complete linkage. Drugs were manually grouped according to mechanism of action and response patterns. Missing values are indicated in gray. The boxplot below shows the distribution of DSS values per PDO, with 25^th, 50^th, and 75^th percentiles. For pairs of highly correlated drugs, only one was included in the calculations (see Materials and Methods). C, Association of the TP53 gene expression signature and TP53 mutation status with idasanutlin sensitivity. The P values are from Wilcoxon rank sum tests. Δ DSS indicates difference in medians. For multi-lesion samples, the mean DSS of the lesions was used as input. Abbreviations: DSS, drug sensitivity score; p, patient; ssGSEA, single-sample gene set enrichment analysis.

Most CRLMs have strong ex vivo vulnerabilities to selected clinical drugs

A heatmap of sensitivities to the thirty-eight evaluated drugs (drug-wise median-centered DSS values; Fig. 2B) showed that the PDOs formed clusters based primarily on the sensitivity to EGFR and/or MDM2 inhibition, as well as to SN-38 and nucleosides. The difference in sensitivity to MDM2 inhibition was defined by TP53 activity (Fig. 2C), and the four PDOs (from three patients) with the most “TP53 wild-type–like” phenotype (based on a gene expression signature of the five most upregulated genes in TP53 wild-type compared with mutated primary colorectal cancers) had outlier sensitivity to idasanutlin (P = 0.0071, Wilcoxon rank-sum test). TP53 mutation analyses confirmed wild-type status for two of the sensitive PDOs, but identified heterozygous mutations in the two PDOs with the strongest sensitivity (both from patient 18), and additionally indicated wild-type status in a PDO with only intermediate sensitivity (TP53 wild-type versus mutated: P = 0.012, Wilcoxon rank-sum test; Fig. 2C; Supplementary Fig. S6). Notably, the third PDO from patient 18 was completely resistant to MDM2 inhibition by idasanutlin.

In contrast with EGFR inhibition, sensitivity to MEK inhibition was not associated with RAS mutation status (P = 0.97 for binimetinib and P = 1.0 for trametinib, Wilcoxon rank sum test; Fig. 2A) and two RAS-mutated lesions from each of two patients (patient 3, KRAS mutated and patient 10, NRAS mutated) showed apparent intra-patient heterogeneity in MEK inhibitor sensitivities (patient 3: ΔDSS = 7 for trametinib and 6 for binimetinib; patient 10: ΔDSS = 9 for trametinib and 4 for binimetinib). Notably, the RAS wild-type PDO derived from patient 1 was resistant to both EGFR and MEK inhibitors. This PDO harbored the rare BRAF^D594G mutation and was also resistant to the BRAF^V600 inhibitor encorafenib (Fig. 2B).

For approved drugs other than those targeting EGFR or MEK, we observed mostly moderate differences in drug sensitivity among the PDOs (Fig. 2B). The activity of 5-FU and the multi-kinase inhibitor regorafenib was generally low, although we observed fair differential drug responses below the reported maximum blood plasma concentrations (C_max; Supplementary Fig. S7). Encorafenib showed no activity, consistent with wild-type BRAF codon 600 status in all PDOs. No/limited activity was also found for oxaliplatin and TAS-102, but this was likely due to too low drug concentrations. In contrast, the active metabolite of irinotecan, SN-38, showed high activity in nearly all the PDOs with a moderate to strong differential activity. Fifteen of the twenty-two patients were treated with neoadjuvant combination therapies (5-FU, oxaliplatin, irinotecan, bevacizumab, cetuximab, and panitumumab combinations) and the PDOs showed ex vivo sensitivities to 5-FU, SN-38, and EGFR inhibitors corresponding with the tumor responses according to RECIST, implying recapitulation of pharmacological responses in the pre-clinical models (oxaliplatin was not compared due to low activity at the tested concentrations; Supplementary Table S3).

Strong differential drug sensitivities among the PDOs were observed for the antimetabolites methotrexate and gemcitabine (Fig. 3A; Supplementary Figs. S8 and S9), as well as for a number of targeted drugs that have no decisive genomic markers, but that are currently tested in late-phase clinical trials in patients with mCRC and/or other cancer types. These targeted agents included the Aurora A kinase inhibitor alisertib (phase 3), the Bcl inhibitor navitoclax (phase 2), the HSP90 inhibitor luminespib (phase 2), the dual PI3K/mTOR inhibitor gedatolisib (phase 3), and the PLK1 inhibitor volasertib (phase 3; Fig. 3B; Supplementary Figs. S8 and S9). There were no clear responders to crizotinib (ALK inhibitor) or LGK974 (PORCN inhibitor), likely related to the low frequency in mCRC of the genomic aberrations associated with sensitivity to these agents. In accordance, no outlier expression profiles were observed for the corresponding target genes (Supplementary Fig. S10).

• Download figure
• Open in new tab
• Download powerpoint

Figure 3.

Drugs with strong differential activity and personalized therapy nomination. Approved drugs (A) and investigational drugs (B) in late-phase clinical studies with strong differential activity across the PDO cohort. C, Illustration of use of ex vivo drug response profiling to nominate potential therapies for individual patients; here, patient 29. Error bars in the vertical bar plot indicate variation in median-centered DSS values between the three screened PDOs from the patient. Abbreviations: cDSS, median-centered drug sensitivity score; C_max, maximum plasma drug concentration; E_max, maximum drug activity; moa, mechanism of action; p, patient.

Promising therapies could be nominated for eighteen (82%) of the patients based on the ranked relative DSS values and E_max (% viability) across all drugs (Table 1; Supplementary Fig. S11). The PDOs from three metastatic lesions from patient 29 (Fig. 3C) illustrate the several heuristic parameters used for drug nominations (unambiguous dose–response curves, E_max < 50% viability, median-centered DSS > 5, and low DSS variation among metastases from the same patient). PDOs from seven of the patients (patients 6, 7, 9, 11, 13, 20 and 21) were generally sensitive to several drugs, including standard of care therapies (Table 1). PDOs derived from patients 25 and 29 were sensitive to drugs approved for other indications and drugs in phase 2/3 clinical trials, and PDOs from patients 1, 18, and 28 were sensitive only to drugs in development (in phase 2/3 trials). No drug candidates could be nominated for four (18%) of the patients (patients 3, 10, 12 and 23) according to the defined criteria. However, based on absolute (in contrast with relative) drug activity, the three RAS wild-type PDOs from patient 12 were all sensitive to EGFR inhibition (Fig. 2B).

There was large variation among the PDOs in the overall level of drug sensitivities and the number of strong drug vulnerabilities (Fig. 4A). The two PDOs from patient 26 showed an overall strong sensitivity pattern [median DSS 8.2 and 7.2; interquartile range (IQR): 13 and 14 for t1 and t2, respectively; Fig. 4A], whereas the PDO from patient 8 was resistant to all drugs beyond the few vulnerabilities listed in Table 1 (median DSS 2.3; IQR: 0.3–8.9). This latter PDO also had mesenchymal-like gene expression characteristics that have been linked to drug resistance, including upregulated expression of gene sets related to “TGFβ,” “EMT,” “LGR5 stem cells,” and “WNT,” as well as downregulation of the “differentiation” and “CDX2” gene sets (Fig. 4B). Reduced CDX2 and CK20 protein expression further supported a poorly differentiated phenotype (Fig. 4B). To explore the potential for more efficient growth inhibition of this multi-drug–resistant and KRAS^G12V-mutated PDO, we designed a combination screen based on the strongest observed vulnerability: sensitivity to inhibition of the PI3K/AKT/mTOR pathway (Fig. 4C), which was associated with the apparent loss of PTEN expression (Fig. 4D). Combined mTOR or AKT inhibition with MEK inhibition may be synergistic in this genetic background (42, 43), and both combination screens with everolimus (mTOR) or ipatasertib (AKT) plus trametinib suggested a synergistic negative effect on PDO viability, as evaluated by a zero interaction potency model (Fig. 4E).

• Download figure
• Open in new tab
• Download powerpoint

Figure 4.

Pharmacogenomic relationships for PDOs. A, Distribution of drug activities in PDOs ranked by 75^th percentile. B, Gene set enrichment analysis of PDOs showing median-centered scores for fourteen pre-specified gene sets (see Materials and Methods). Samples and gene sets were clustered according to correlation distance and average linkage. Protein expressions for CDX2 (caudal type homeobox 2 transcription factor) and CK20 (cytokeratin 20) are shown below. C, Waterfall plot of median-centered DSS for patient 8. Higher values indicate higher activity relative to other PDOs. The four drugs with cross-dataset range in DSS <5 were excluded. D, PTEN expression in PDOs, showing particularly low expression levels for patients 8, 15, and 30. E, Drug synergy testing between MEK inhibition and mTOR or AKT inhibition in PDO with PTEN loss. cDSS, median-centered drug sensitivity score; DAB, 3,3′-diaminobenzidine; p, patient; ssGSEA, single-sample gene set enrichment analysis; ZIP, zero interaction potency. *The ipatasertib plus trametinib combination was performed in replicate due to the overall technical variability observed for ipatasertib (Supplementary Fig. S4C), showing comparable results (ZIP = 5, ZIP_max = 10).

Intra-patient inter-metastatic pharmacotranscriptomic heterogeneity is not pronounced

Comparisons of the pharmacological profiles of PDOs from multiple CRLMs from each of ten patients showed strong intra-patient correspondences in drug sensitivities to the majority of the tested drugs. The multiple PDOs largely clustered patient-wise, both in hierarchical clustering and PCA (Figs. 2B and 5A). By careful evaluation of all dose–response curves we observed convincing inter-metastatic heterogeneity (not associated with technical variability in replicate screens) in less than one tenth of all the 38 × 10 drug–patient comparisons, including in patient 3 (MEKi), patient 10 (MEKi and navitoclax), patient 18 (idasanutlin), patient 23 (EGFRi/MEKi) and patient 29 (LCL 161). The differential response to idasanutlin in patient 18 was associated with heterozygous versus homozygous TP53 mutation status in the sensitive versus resistant PDOs, but all PDOs from this patient had a TP53 wild-type–like phenotype by gene expression. Furthermore, although a gene expression signature of sensitivity to EGFR inhibition (23) was correlated with erlotinib sensitivity across all PDOs, this could not explain the inter-metastatic pharmacological heterogeneity to EGFRi/MEKi among five RAS/BRAF wild-type PDOs from patient 23 (Supplementary Fig. S12).

• Download figure
• Open in new tab
• Download powerpoint

Figure 5.

Correspondence between drug sensitivity and gene expression variation. A, Principal component analyses of drug sensitivity (y-axis) and gene expression profiles (x-axis) of PDOs. PDOs derived from multiple CRLMs from each of ten patients are color coded and indicated with patient numbers. All other PDOs are shown in gray. Two drugs with frequent NAs (linsitinib and ipatasertib), and the seven drugs with lowest cross-sample variance (vitamin C, oxaliplatin, encorafenib, disulfiram, crizotinib, TAS-102, and tretinoin) were excluded from PCA analyses of drug sensitivity. B, Distributions of Spearman's correlation coefficients for drug sensitivity and gene expression separately. Correlations for related and unrelated PDOs are indicated. C, Scatterplot of Spearman's correlation coefficients between pairs of PDOs at the drug sensitivity (horizontal axis) and gene expression (vertical axis) levels. Linsitinib and ipatasertib were excluded also from analyses in (B) and (C) due to frequent NAs. Each dot represents one comparison between two PDOs, and comparisons of PDOs from the same patient are color coded, whereas comparisons of unrelated PDOs are shown in gray. Note: the PCA function does missing value imputation (by mean of variable). Abbreviations: p, patient; PC, principal component.

Gene expression profiles were similarly compared to investigate the extent to which the variation in drug sensitivities was recapitulated at the mRNA level. Gene set analyses in PDOs confirmed that PDOs from the same patients clustered more closely in hierarchical clustering analyses than unrelated PDOs also at the gene expression level (Fig. 4B). Furthermore, PCA comparisons showed a significant correlation of PC1 values at the gene expression and drug sensitivity levels (Spearman's ρ 0.48, P = 0.002), largely with sample-wise clustering of PDOs in both dimensions (Fig. 5A). Intra-patient correspondence was further confirmed by comparisons of pairwise Spearman's correlation coefficients for either drug sensitivity or gene expression between different tumors from the same patient and between tumors from different patients (Fig. 5B). PDOs from the same patients had stronger correlations at both levels, indicating lower intra-patient compared with inter-patient heterogeneity. Pairwise comparisons of Spearman's correlation coefficients for DSS and gene expression between PDOs demonstrated that both unrelated and related PDOs were correlated at the two different data levels, although PDOs from different lesions from the same patient had the strongest correlation (Fig. 5C).