Fresh tissues and paraffin-embedded wax blocks were consecutively derived from patients who underwent surgical resection of CRC at Taipei Veterans General Hospital, Taiwan from February 13, 2017 to October 13, 2018. All studies were approved by the Institutional Review Board of Taipei Veterans General Hospital, Taiwan as well as the Institutional Review Board of Academia Sinica, Taiwan. Written consent forms were obtained from all the patients. In two patients, non-tumorous tissues were not obtained. Paired primary tumor tissues from the colon along with liver metastases were obtained from 3 patients. The average size of the tumor samples used to prepare PDO is 4-5 mm in diameter. Non-tumorous tissue was minced, digested with collagenase and dispase, and processed as described (14). The crypts were collected, embedded in matrigel, and cultured in Basal culture medium supplemented with 50% Wnt conditioned medium, 20% R-Spondin conditioned medium, 10% Noggin conditioned medium, 1x B27, 1.25 mM n-Acetyl Cysteine, 10 mM Nicotinamide, 50 ng/ml human EGF, 10 nM Gastrin, 500 nM A83-01, 3 M SB202190, 10 nM Prostaglandin E2, and 100 mg/ml Primocin (In vivogen). Tumor tissue was minced, digested with Liberase, and processed as described previously (15). Tumor organoids were cultured in the same medium described above except without Wnt. Similarly, organoids were prepared from metastatic tumors from liver site in the same medium without Wnt, and adjacent non-tumorous liver tissues in media supplemented with 100 ng/ml FGF10 (Peprotech), 25 ng/ml HGF (Peprotech) and 10 μM Forskolin (Tocris) (16). In general, we established the organoids from the fresh tumor samples as passage 0. The organoids were passaged once in 6 days using TrypLE Express (Thermo Fisher Scientific) in 1:3 to 1:5 dilution based on cell density. After initial passages, the organoids (at passages 3 and 4) were frozen in a standard cell freezing media containing 10% dimethyl sulfoxide (DMSO) and stored at -80°C for 2 to 3 days and then shifted to liquid nitrogen tank for long term storage as organoid bio-bank.

In brief, we have established a bio-bank of 151 CRC-PDOs from 148 CRC patients, including 3 patients with matching liver metastatic tissues. The utilization of PDOs in this study is summarized. PDOs derived from advanced stages (III and IV) (n = 42) were subjected to in vitro drug sensitivity assay against oxaliplatin and divided into oxaliplatin-sensitive (OS) and -resistant (OR) groups. Four OR and 4 OS PDOs were used in xenograft tumor model study in mice for in vivo validation of the in vitro drug response. In addition, the drug responses in 17 PDOs derived from stage III (n = 12) and stage IV (n = 5) patients were further compared to the clinical outcomes of these patients after FOLFOX adjuvant therapy. To establish personalized medicine, 8 OS and 8 OR PDOs were subjected to RNA sequencing to build up molecular portraits of transcriptomes, followed by gene signature-based selection of candidate drugs.

Whole mount immunofluorescence staining of organoids was performed as described earlier (17). Organoids were washed in ice-cold wash buffer (DMEM/F12 with 10% FBS), and the matrigel domes disrupted. Organoids were pelleted down and fixed with 4% paraformaldehyde in PBS for 30 min, and permeabilized with 0.5% Triton X-100 in PBS for 15 min. Organoids were washed in immunofluorescence (IF) wash buffer (0.1% BSA, 0.2% Triton X-100, 0.05% Tween 20), incubated with blocking buffer (1% BSA in IF wash buffer) for 1 h, and incubated with primary antibody rabbit anti-EpCAM (AbCam), rabbit anti-Mucin 2 (Santa Cruz), mouse anti-CK20 (Santa Cruz), or mouse anti-Chr-A (Santa Cruz) in blocking buffer for 2 h, followed by incubation with secondary antibody anti-mouse Alexafluor Texas Red (Invitrogen) or anti-rabbit or mouse Alexaflour 488 (Invitrogen) in blocking buffer for 1 h. After mounting in SlowFade anti-fade with DAPI (ThermoFisher), organoids were observed and micro graphed using Carl Zeiss LSM510 laser scanning microscope imaging system with Zen analysis software (Zeiss). Pathological sections of formalin-fixed paraffin-embedded patient tissue, xenografts and organoids were processed, imaged and analysed in the Histopathology unit at Institute of Biomedical Sciences, Academia Sinica. For organoids histology, whole organoids were pelleted down and fixed as described above. Then organoids were suspended in 2% low melting agarose, followed by dehydration, paraffin embedding, sectioning and Hematoxylin Eosin (H&E) staining.

Drug sensitivity assay was performed in organoids as described with modifications (14). Organoid cultures were gently disrupted into single cell suspension with TrypLE Express, and 3000 cells were seeded in 7 μl per well of matrigel diluted with equal volume of medium in 96-well plates in four replicates and cultured in organoid growth medium for 24 h. The cultures are incubated with oxaliplatin at various concentrations (0.03, 0.1, 0.3, 1, 3, 10, 30, and 100 μM) for 96 h, and subjected to CellTiter-Glo (Promega) cell viability assay. The slope of the dose-response curve was calculated. IC₅₀ (half maximal inhibitory concentration) values are calculated with the four-parameter nonlinear logistic equation. E_max values were calculated as the percentage of inhibition at the maximum included concentration (100 μM). Based on IC₅₀ values, PDOs were classified into four response categories, and the OR group includes the categories of non -responders and minor responders, and OS group include moderate and strong responders.

Total RNA was extracted from 8 OR and 8 OS PDO lines, using RNeasy Plus Mini Kit (QIAGEN). RNA-sequencing libraries were constructed using TruSeq Stranded mRNA Library Prep Kit (Illumina). The purified libraries were amplified by 10 cycles of PCR, and the resulted library profile was 250-500 bp, peaking at ~300 bps. Paired-end 2 x 101-nt sequencing was conducted on the Illumina HiSeq 2500 System (NGS High Throughput Genomics Core Facility at BRACS, Academia Sinica). The short-reads in FAST-Q format, ~85 M reads per sample, were processed by a computational pipeline as described (18). Briefly, the short reads were aligned and mapped to the Homo sapiens (human) genome assembly GRCh38 (hg38) using HISAT2. The aligned transcripts were assembled using StringTie by annotated GENCODE version 29. Gene expression levels were converted to counts by a python script (prepDE.py) provided by StringTie team. RNA-seq datasets consisted of 60,714 gene tags (25,213 Ensembl gene ids and 35,501 StringTie MSTRG ids) and 225,295 transcript tags (206,696 Ensembl transcript ids and 18,599 StringTie MSTRG ids), respectively, from 16 CRC PDOs. Fragments per kilobase of transcript per million mapped reads (FPKM) expression values from individual organoids were obtained. For comparison of differentially expressed genes (DEGs) and tags (DETs) between different groups, gene expression levels were further normalized by DESeq method and analyzed using Bioconductor package DESeq2.

The classification of the 16 PDOs was visualized by principal component analysis (PCA) (R package: ropls) and unsupervised hierarchical clustering analysis (Morpheus; Metric: One minus spearman rank correlation and Linkage method: average; https://software.broadinstitute.org/morpheus) using the DETs (fold change ≥ 2 and adjusted p-value < 0.05) as the identifier input. Analysis of the hallmarks and KEGG pathway enriched in the OR and OS PDOs were conducted by Gene Set Enrichment Analysis (GSEA). False discovery rate (FDR) q value < 0.05 was considered significant.

Leave-One-Out Cross Validation (LOOCV) was conducted to select oxaliplatin response gene predictors (19). Machine learning methods, such as Random Forest, K-Nearest Neighbors, Naïve Bayesian, Decision Trees, Neural Networks, Support Vector Machines, and Logistic Regression (LR), were applied to validate the prediction power of the sets of DEGs and DETs identified by LOOCV.

Connectivity MAP (C-MAP) online tool (https://www.broadinstitute.org/cmap/) was conducted to predict inhibitors targeting oxaliplatin resistant signatures (20, 21). Perturbagens ranked with the most high average negative enrichment scores and least variations were selected for further evaluation in combination therapy.

In vivo oxaliplatin response was assessed using PDO-based xenograft tumor model. All experimental protocols were approved by the Institutional Animal Care and Utilization Committee (IACUC) at Academia Sinica, Taiwan. Single cell suspension was prepared from each PDO, and 1×10⁶ cells were mixed with organoid growth medium/matrigel (1:1) in a total volume of 100 uL and injected subcutaneously to the right flank of NSG (NOD scid gamma) mice (8-9 weeks, 23-28 g). Tumors were allowed to grow to 80 ± 100 mm³ in size and randomly assigned into the control (2.5% DMSO/PBS) and treatment (oxaliplatin at 2 mg/kg in 2.5% DMSO/PBS) groups with 5 mice per group. Treatment was given by intra-peritoneal injection for 10 times at 3-day intervals. Tumor size and body weight were measured. TGI (Tumor Growth Inhibition) was determined at the end point of treatment period for each xenograft by the following formula: %TGI = {[1 – (T_t/T₀)/(C_t/C₀)]/[1 – (C₀/C_t)]} x 100, in which Tt and T0 were median tumor volumes of treated animal at time t and time 0, and C_t and C₀ median tumor volumes of control group at time t and time 0, respectively (22). Median %TGI was calculated and reported for each group. Significant anti-tumor activity was defined as achievement of a median %TGI of at least 50%. Repeated two-way ANOVA was used to define the statistical difference between the groups using Graph Pad Prism 9.

Statistical analysis was performed in GraphPad Prism 9 (GraphPad Software, San Diego, CA). IC₅₀ for in vitro drug assay was determined by non-linear regression log(inhibitor) vs response variable slope (four parameters). To calculate p-value in Supplementary Figure 2 , Mann-Whitney test (un-paired two tailed) was used to compare the difference in IC₅₀, E_max, and AUC between the groups of minor/non-responders and moderate/responders. To calculate p-value in Figure 3A , we used Repeated Measure two-way ANOVA with the Geisser-Green house Correction and Bonferroni multiple comparison test. For all required studies, experimental mean was calculated and error was presented by standard deviation of the mean. For Supplementary Figure S3 , significant difference of relapse-free survival depending on the PDO resistance or sensitivity to oxaliplatin was determined using Gehan-Breslow-Wilcoxon test. P<0.05 was considered as significant.

Tumorous and non-tumorous tissues, plus liver metastases (n=3), were collected from 148 CRC patients, for organoid preparation. These CRC patients, categorized based on TNM staging, included 3 cases of stage 0, 24 of stage I, 42 of stage II, 57 of stage III, and 22 of stage IV, corresponding to 2%, 16.2%, 28.4%, 38.5% and 14.9%, respectively, of total samples collected. Supplementary Table S1 lists the description of samples collected, and the 3 pairs of primary and metastatic tumor samples are indicated. The success rate of establishing organoid cultures was 93% (136 of 146) and 76% (115 of 151) for the non-tumorous and tumorous tissues, respectively ( Figure 1A ), which was comparable to previous reports (9, 23, 24). The unsuccessful cases included the initial attempts for optimization of culture conditions and those with tissues full of necrotic lesions or contaminated with pathogens. Among the established patient-derived tumor organoids (PDOs), 3 quarters of them, spanning stage 0 (n = 1), stage I (n = 13), stage II (n = 26), stage III (n = 30), and stage IV (n = 17), were characterized with good growing condition, and can be expanded for continuous passages (passage number > 10). Wnt3A was required for culturing organoids derived from non-tumorous tissues, but not for most, if not all, of the tumor-derived organoids. Non-tumor organoids can only be cultured and maintained under normoxic condition, whereas tumor organoids were able to proliferate under both normoxic and hypoxic conditions. We routinely cultured tumor organoids under normoxia in the absence of Wnt3A. Immunofluorescence staining was performed to show the viability of the organoids ( Figure 1B ). Almost all the cells in both non-tumor and tumor organoids were stained positively by the live-cell-permeable green fluorescent dye Calcein-AM. The organoids were stained positive for the epithelial cell adhesion molecule (EpCAM). The non-tumor organoids were readily stained for mucin 2, cytokeratin-20 and chromagranin A, demonstrating the presence of differentiated cell types of intestinal epithelial lineage, whereas most of the tumor organoids were stained weakly or negative for the differentiation markers ( Figure 1C ). CRC-derived organoids cultured under in vitro or in vivo xenograft condition retained architectural features resembling the tumor tissues from which they were derived ( Supplementary Figure S1 ). Well and moderately differentiated tumors gave rise to cyst-like organoids with clear lumen formation, whereas the organoids derived from poorly differentiated tumors were characterized with highly compact structures ( Figure 1D ). Sequencing analyses revealed that PDOs preserved most of the mutational status of the driver genes in the primary tumors ( Figure 1E ). We noted that there were mutations either gained or lost in a small number of PDO lines during the in vitro culturing process.

Adjuvant 5-FU in combination with oxaliplatin serves as the first-line chemotherapy for patients of advanced CRC (25, 26). With the aim to address whether PDO lines can serve as a preclinical model to guide therapeutic decisions for advanced CRC patients, we performed drug sensitivity assessment of advanced CRC PDOs towards oxaliplatin. Forty-two of the 47 advanced CRC PDOs were successfully recovered from frozen stocks without showing any decline in growth rate, and subjected to the treatment of increasing concentrations of oxaliplatin. The dose response curve of individual PDOs was constructed ( Figure 2A ). The IC₅₀ ranged from 1.37 to 100 μM, varying in 2 orders of magnitude in individual PDOs, and E_max ranged from 45 to 93%. An inverse correlation was established between the potency and efficacy of oxaliplatin in the 42 PDOs ( Figure 2B ). The wide range of differences in the potency and efficacy of oxaliplatin emphasizes the need for personalized medicine. According to IC₅₀, the PDOs were grouped into four categories ( Figure 2C ), strong responders (log₁₀IC₅₀ 0.14 – 0.57), moderate responders (log₁₀IC₅₀ 0.57 – 0.98), minor responders (log₁₀IC₅₀ 0.98 – 1.24), and non-responders (log₁₀IC₅₀ 1.24 – 2.00), as previously described⁶. Strong and moderate responders were defined as oxaliplatin-sensitive (OS) group whereas minor and non-responders were defined as oxaliplatin-resistant (OR) PDOs. Individual dose response curves of oxaliplatin-sensitive and -resistant PDOs are presented in Supplementary Figure S2A . Quantification of responses to oxaliplatin by calculating the IC₅₀, E_max and the area under the dose response curve (AUC_DRC) showed a significant difference between the resistant (non-responders/minor responders) and sensitive (moderate/strong responders) groups ( Supplementary Figure S2B ). Supplementary Table S2 summarizes the clinical data of the 42 advanced CRC PDOs, along with their in vitro oxaliplatin response.

We next examined whether the oxaliplatin sensitivity determined in PDOs in vitro can be validated in vivo. As shown in Figure 3A , oxaliplatin-resistant (OR) and -sensitive (OS) PDOs were subjected to xenograft tumor model in NSG mice, followed by oxaliplatin treatment. Oxaliplatin-mediated tumor growth inhibition (TGI) was calculated as described (22). As shown, oxaliplatin treatment yielded significant inhibition (%TGI > 50%) in all of the 4 OS PDO xenografts. To the contrast, none of the OR PDO xenografts responded well to oxaliplatin-mediated inhibition. These results demonstrate a good concordance between the oxaliplatin sensitivity in PDOs and drug response in PDO xenograft tumors in vivo.

Next, we analyzed the clinical outcomes of CRC patients who had received oxaliplatin-based therapy. Among the 42 advanced CRC PDOs, 17 were derived from patients (12 of stage III and 5 of stage IV) who underwent FOLFOX adjuvant therapy. Supplementary Table S3 describes the details of the patient cohort from whom these 17 PDOs were obtained. According to oxaliplatin sensitivity assessed in vitro, the 17 PDOs were 4 non-responders, 4 minor responders, 5 moderate responders, and 4 strong responders to oxaliplatin. Notably, the patients who developed lung or liver metastasis after adjuvant therapy belonged to 3 of the 4 non-responders, 2 of the 4 minor responders, 1 of the 5 moderate responders, and 1 of the 4 strong responders ( Figure 3B ), demonstrating that FOLFOX treatment response in advanced CRC patients followed the trend of oxaliplatin sensitivity determined in PDOs (P = 0.047, using single tailed Cochran-Armitage trend test). Further, Kaplan Meier relapse-free curve of CRC patients depending on PDO resistance or sensitivity to oxaliplatin showed a significant difference (P = 0.047, using Gehan-Breslow-Wilcoxon test) ( Supplementary Figure S3 ). To summarize clinical validity of PDO-based drug screen results for predicting treatment response, we found 70% sensitivity, 71.4% specificity, 77.8% positive predivtive value, and 62.5% negative predictive value in predicting response to oxaliplatin-based regimen (Fisher’s exact test p = 0.15).

To investigate the molecular mechanisms modulating oxaliplatin response, we analyzed the expression profiles of 8 OS and 8 OR PDOs. RNA sequencing was performed. The short reads in Fast-Q format were processed through the pipeline of HISAT2, StringTie and DESeq2 algorithms, and 2555 differentially expressed transcripts (DETs, adjusted p < 0.05 and fold change ≥ 2) were obtained. By principal component analysis and unsupervised hierarchical clustering analysis ( Figures 4A, B ), these DETs yielded a clear separation between the OS and OR PDOs. Morpheus Marker selection method (https://software.broadinstitute.org/morpheus) was applied, using a permutation test with 1,000 replications, and transcripts correlated to the OR and OS phenotypes were also identified. Supplementary Figure S4 shows the heatmap of top 25 up- and 25 down-regulated transcripts associated with the OR phenotype in the form of hierarchical clustering (Spearman’s correlation), demonstrating again that OS and OR PDOs exhibited distinct expression profiles.

By Gene Set Enrichment Analysis (GSEA) (www.broad.mit.edu/gsea), hallmarks of mitotic spindle, G2/M checkpoint, UV response and TGF-β signaling were identified to be tightly associated with OR phenotype, whereas oxidative phosphorylation hallmark associated with OS phenotype ( Figure 4C ). To validate our findings, we further queried RNA-sequencing data of human CRC (459 of COAD and 170 of READ) database in The Cancer Genome Atlas (TCGA) program. Among the patients who received oxaliplatin-based chemotherapies with clinical outcomes recorded, 33 were evaluable for having RNA-Seq data ( Figure 5A ). According to treatment response, 23 patients were considered to be resistant to the treatment for developed stable disease (SD) or progressive disease (PD), whereas 10 exhibited complete response (CR) or partial response (PR) as responding or sensitive to the treatment. By GSEA analysis, the pathways and gene sets associated with CRC patients resistant and sensitive to oxaliplatin-based treatments were shown to be similar to those enriched in the OR and OS PDOs ( Figures 5B, C and Supplementary Figure S5 ).

To identify oxaliplatin response predictor signatures, Leave-One-Out Cross Validation (LOOCV) algorithm was performed to compare the expression profiles of 8 OR vs. 8 OS PDOs ( Figure 6A ). Sixteen new expression datasets were created by alternatively leaving out the RNA-Sequencing dataset of one PDO. Comparative analysis was performed in each dataset, and the top 100 DEGs and DETs between the OR and OS PDOs were identified. Among the 16 sets of top 100 DEGs and DETs, 36 DEGs and 33 DETs were found commonly present in all the 16 sets of genes and transcripts. Several machine learning programs were applied to evaluate the accuracy of using the sets of 36 DEGs and 33 DETs as predictors for oxaliplatin response. As shown, the set of 36 DEGs or 33 DETs gave a perfect or near perfect value of 1 as prediction accuracy to serve as drug response predictor signatures. We further exercised to test whether any single gene in these DEGs and DETs can serve as a drug response predictor. Using Logistic Regression (LR) model, we identified 18 models using one single gene and 20 models using one single transcript that reached the score of 1 to predict oxaliplatin response in the PDOs ( Figure 6B ). The identity of the 18 oxaliplatin response predictor signatures and their expression status in OR and OS were shown ( Figures 6C, D ). These predictor signatures also fell within the most significantly differentially expressed genes identified by GSEA ( Supplementary Table S4 ).

By GSEA, hallmarks of mitotic spindle and G2/M checkpoint were identified to be tightly associated with the OR phenotype. Supplementary Table S5 lists the gene sets associated with these two hallmarks that were enriched in the OR phenotype, and a cohort of genes were commonly enriched in these two hallmarks, including c-ABL and NOTCH2. We hypothesized that inhibition of these two hallmarks may attenuate oxaliplatin resistance. Small molecules of c-ABL inhibitor imatinib and DAPT which targets NOTCH pathway by inhibiting γ- secretase were chosen to be tested.

In a second approach, we submitted the top 150 up-regulated and 150 down-regulated genes associated with OR phenotype by GSEA to the Connectivity Map (CMAP) (20), a genomics-based drug discovery framework, and identified 24 perturbagens with enrichment scores of ≤ -95, implicating that they may target OR signatures. Similarly, the 15 predictor genes derived from LOOCV that were up-regulated in the OR PDOs were also subjected to C-MAP analysis, and 128 compounds with enrichment scores ≤ -95 were identified. Eighteen compounds were found commonly shared in these two groups of perturbagens, including HDAC inhibitors, DNA damaging agents, DNA/RNA synthesis inhibitor, and inhibitors for druggable targets, such as CDK inhibitor and HGFR inhibitor ( Figure 7A ). Six compounds with relatively high negative enrichment scores, including irinotecan, mitomycin-c, HDAC inhibitors vorinostat, scriptaid, and trichostatin A (TSA), and mycophenolate-mofetil (MPM), a selective inhibitor of inosine monophosphate dehydrogenase that may inhibit DNA/RNA synthesis, were chosen to be tested for their activity to increase oxaliplatin sensitivity. Cytotoxicity of these compounds in PDOs was assessed, and a minimum effective dose, between IC₁₀ – IC₄₀, was chosen for each compound in combination therapy with oxaliplatin. Dose response curve for oxaliplatin alone and in combination with individual compounds was constructed, and relative IC₅₀ of oxaliplatin resulted from the combination therapy and the Emax achieved in OR PDOs were shown in Figure 7B . In comparison to oxaliplatin monotherapy, most of the combination therapies significantly increased the potency and efficacy of oxaliplatin, as evidenced by the decrease of relative IC₅₀ by 1 to 2 orders of magnitude (top panel), and the increase of E_max by 10 - 20% (bottom panel). Since these candidate drugs were selected based on their potential to target OR phenotype, as expected, favorable outcomes were not expected when treating the OS PDOs ( Supplementary Figure S6 ). Notably, a much worse effect was observed when applying the combination therapies to the OS PDOs, as evidenced by the dramatically increased IC₅₀ and decreased E_max. Supplementary Table S6 summarizes the IC₅₀ values of individual compounds and in combination with oxaliplatin. These data further emphasize the importance of gene signature-based approach for personalized therapy. Most importantly, these data strongly support the usefulness of PDOs as a platform for the development of gene-based therapies.