Skip to main content

“Proteotranscriptomic analysis of advanced colorectal cancer patient derived organoids for drug sensitivity prediction”

Abstract

Background

Patient-derived organoids (PDOs) from advanced colorectal cancer (CRC) patients could be a key platform to predict drug response and discover new biomarkers. We aimed to integrate PDO drug response with multi-omics characterization beyond genomics.

Methods

We generated 29 PDO lines from 22 advanced CRC patients and provided a morphologic, genomic, and transcriptomic characterization. We performed drug sensitivity assays with a panel of both standard and non-standard agents in five long-term cultures, and integrated drug response with a baseline proteomic and transcriptomic characterization by SWATH-MS and RNA-seq analysis, respectively.

Results

PDOs were successfully generated from heavily pre-treated patients, including a paired model of advanced MSI high CRC deriving from pre- and post-chemotherapy liver metastasis. Our PDOs faithfully reproduced genomic and phenotypic features of original tissue. Drug panel testing identified differential response among PDOs, particularly to oxaliplatin and palbociclib. Proteotranscriptomic analyses revealed that oxaliplatin non-responder PDOs present enrichment of the t-RNA aminoacylation process and showed a shift towards oxidative phosphorylation pathway dependence, while an exceptional response to palbociclib was detected in a PDO with activation of MYC and enrichment of chaperonin T-complex protein Ring Complex (TRiC), involved in proteome integrity. Proteotranscriptomic data fusion confirmed these results within a highly integrated network of functional processes involved in differential response to drugs.

Conclusions

Our strategy of integrating PDOs drug sensitivity with SWATH-mass spectrometry and RNA-seq allowed us to identify different baseline proteins and gene expression profiles with the potential to predict treatment response/resistance and to help in the development of effective and personalized cancer therapeutics.

Background

Colorectal cancer (CRC) is the third most common and second deadliest cancer [1]. Even when diagnosed at early stage, 30–50% of patients will experience relapse and die of the disease [2]. Whether relapsed or diagnosed as stage IV disease, survival has significantly improved over the last ten years, going from 6 to more than 30 months median overall survival [3]. However, despite improvements in cancer therapy, resistance to chemotherapeutic and novel targeted therapies limits treatment success.

In light of the increasing notion of CRC molecular complexity, the role of useful preclinical models for predictive biomarkers and new drug discovery is gaining importance. Furthermore, there is an urgent need to go beyond genomics, as its utility to guide therapy has been less successful than anticipated [4]. Indeed, the presence of a molecular alteration at the genomic level does not automatically predict treatment response, as multiple factors can influence response in each individual patient to a greater or lesser extent [5].

Functional testing using human-derived cancer models offers a valuable complementary approach to help in personalizing treatments. Patient-derived organoids (PDOs) are employed to fill the gap left by traditional preclinical models. They recapitulate the biology of original tumor tissue, including therapy response [6], and can be rapidly generated for many cancer types [7,8,9,10] with a lower cost than animal models. A great number of PDO models have been generated for CRC [6, 11,12,13]. However, the evidence obtained until now with these models comes mainly from cultures deriving from primary tumors, while data obtained from metastatic samples are less consistent.

Understanding the phenotypic alterations that play a role in differential response to drugs will be pivotal in identifying predictive biomarkers and developing efficacious therapies for advanced CRC. The bulk of research to date has focused on the value of organoids as a predictive tool [6, 12, 14]; some have characterized the transcriptomic profile in relation to drug sensitivity [13, 15] in an effort to identify predictive biomarkers, but only a limited number of studies have used quantitative proteomic analysis to characterize PDOs [16, 17]. Extensive literature indicates that transcriptomic and proteomic expression profiles lack correlation due to significant post-transcriptional regulation, which could explain why RNA signatures are rarely useful as drug response predictors in the clinic [18].

In this study we generated PDO models from advanced CRC patients including heavily pre-treated patients and a model of microsatellite unstable- (MSI) high CRC, for which pre- and post-chemotherapy liver disease PDOs were established and characterized to ensure biological fidelity with the original tumors. Our main objective was to integrate functional drug assays conducted on PDOs with proteotranscriptomic study to explore the mechanisms underlying drug response. To this end, we used SWATH-MS (sequential window acquisition of all theoretical mass spectrometry), a highly accurate proteomic quantification technique and a bulk RNA-sequencing analysis and performed an integrative functional network analysis of both omics. This strategy could represent a useful tool to discover new therapeutic targets in advanced CRC and to guide personalized therapies.

Methods

Tissue processing and organoid culture

Locally advanced or metastatic colorectal samples were collected after patients had signed written informed consent at Hospital Clínico Universitario de Valencia (HCUV). The HCUV Ethics Committee (EC) approved the study (2018/063, 2021/083). Samples were collected at surgery or with image-guided tissue biopsy. Both naïve and pre-treated patients were included.

Fresh tissues were collected in PBS with antibiotics and quickly processed as follows. After serial washings and antibiotic incubation, tissue samples were minced in small fragments. Some of them were stored for DNA and RNA extraction depending on the total amount of tissue. All the remaining fragments were mechanically and enzymatically digested with a collagenase-based solution (Sigma-Aldrich, cat. No. 269395). Subsequently, the resulting digestion solution was filtered, single cells were counted, resuspended in basal media with 50% reduced-growth factor basement membrane matrix (BME Type 2, R&D, cat. No. 3533–010-002), plated in prewarmed culture plates and stored at 37 °C. After the formation of domes, complete medium was added.

Complete medium composition was as follows: 1 × N2 (cat. No. 17502048), 1 × B27 (cat. No. C21-H), 1 mM N-acetyl L-cysteine (cat. No. A9165), 10 mM nicotinamide (cat. No. N0636), 50 ng/mL hEGF (cat. No. BMS320), 100 ng/mL hNoggin (cat. No. 6057-NG-025), 0.5 mM A-83–01 (cat. No. SML0788), 10 mM SB202190 (cat. No. S7067), 10 nM gastrin (cat. No. G9145), 500 ng/mL hRSPO1 (cat. No. 4645-RS), 10 ng/mL hFGF10 (cat. No. 100–26), 10 nM PGE2 (cat. No. P6532), 10 mM Y-27632 (cat. No. Y0503).

Medium was changed three times a week. Organoid’s formation was observed in a fraction of cultures (see Results section). After reaching appropriate volume organoids were trypsinized with TrypLE™ Express (Life Technologies, cat. No. 12605–010) to expand them. Aliquots were stored in liquid nitrogen to constitute a Biobank.

All samples were named with a codified nomenclature indicating the site of biopsy (metastasis or primary) and with a serial progressive number: metastatic colon tumor organoid (mCTO), metastatic colon tumor tissue (mCTT), rectal tumor organoid (RTO) rectal tumor tissue (RTT) colon normal organoid (CNO) and colon normal tissue (CNT). Collectively, all models are referred as patient-derived organoids (PDOs).

Morphologic characterization

PDOs domes were collected in 4% neutral buffered formalin and paraffin embedded (FFPE). 4 mm slides were cut, dewaxed and hematoxylin and eosin staining (Dako, cat. No. CS700 and CS701) was performed. For IHC sodium citrate antigen retrieval was performed (Target Retrieval Solution, Citrate pH 6, S236984-2 and pH 9, S236784-2), followed by peroxidase blocking (DakoREAL™, Dako, cat. No. S2023) and by incubation with the following primary antibodies in EnVision FLEX antibody diluent (Dako, cat. No. K8006): Ki67 (Dako, cat. No. IR626, ready to use), CDX2 (Dako, cat. No. IR080, ready to use), CK20 (Dako, cat. No. IR777, ready to use). In some PDOs the following additional stainings have been perfomed: ERBB2 (Roche, cat. No. 790–2991, 4B5 clone), PTEN (Dako, cat. No. M3627, 1:50). Peroxidase-conjugated secondary antibodies were used (EnVision Flex/HRP, Dako, cat.no. GV82311-2). The slides were also counterstained with hematoxylin.

Mismatch repair proteins MLH1 (Dako, cat. No. 7R079, ready to use), PMS2 (Dako, cat. No. 7R087, ready to use), MSH2 (Dako, cat. No. IR086, ready to use), MSH6 (Dako, cat. No. IR085, ready to use), PDOs’ staining was performed when clinical data indicated that the patient had a loss of expression of these proteins.

The same procedures were performed also for FFPE tissue sections. All slides from both PDOs and tissues were independently reviewed by dedicated pathologists.

Genomic characterization

Fresh tissues and PDOs DNA was extracted with the QIAamp DNA Micro Kit (Qiagen, cat. No. 56304) to perform both targeted next-generation sequencing (NGS) and copy number variation analysis with CytoScan HD. NGS was performed using OncoSpot v1, an in-house customized 87-gene panel [19]. A minimum of 80 ng of DNA was used for library preparation, using KAPA HyperPlus Library Preparation (Roche Diagnostics). Libraries were paired-end sequenced in an Illumina MiSeq platform with a MiSeq Reagent Kit v2 (300-cycles) kit.

Read quality control was performed with FASTQC v0.11.8 (available online at: http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). Sequencing adapters and low-quality reads were filtered out with fastp v0.11.8 [20]. Filtered reads were aligned against the human reference genome GRCh38 with BWA-mem v0.7.17 (Li, H. (2013). Aligning sequence reads, clone sequences and assembly contigs with BWA-MEM. arXiv: Genomics). PicardTools v2.18.6 (http://broadinstitute.github.io/picard) was applied to eliminate duplicated reads. A combination of LoFreq v2.1.3.1 [21] and Mutect2 v 4.0.5.0 [22] was used to call variants for individual patients. Variants were merged with vcftools [23] and normalized to avoid multi-allelic positions with BCFtools [24]. Then, variants were functionally annotated using Variant Effect Predictor v96 [25]. Additional information was added from the Cancer Hotspot database [26] and a proprietary cancer mutation database containing information from different resources such as PCT MDAnderson database (https://pct.mdanderson.org), commercial cancer panels (Oncomine® and Sequenom), relevant literature and our own expertise. Variants with a depth of coverage greater than 100 × and an allelic fraction of 5% or higher were reported. Pathological variants were selected if present in our proprietary database or located in a cancer hotspot with a population frequency below 1%. Discordant results and low VAFs were manually curated with IGV, v. 2.9.4. Pathogenic variants were summarized using the ComplexHeatmap (v2.4.3) R package [27]. Data available after acceptance.

For copy-number analysis, hybridization-based genomic profiling was performed using CytoScan HD Array (Affymetrix, CA, USA) according to manufacturer’s protocol. CytoScan HD CEL files were processed through Chromosome Analysis Suite (ChAS) software version 4.1 with a single sample algorithm. All samples were manually reviewed, and unbalanced samples were reprocessed with normal diploid algorithm. Obtained Chyp files were analyzed with ChAS and IGV software v. 3.0. Copy number and loss of heterozygosity (LOH) data were retrieved for further analysis.

RNA-sequencing analysis (RNA-seq)

Total RNA was extracted from fresh tissues and PDOs with RNeasy Micro Kit (cat. no. 74004) and integrity verified with TapeStation RNA Analysis ScreenTape (Agilent Technologies). Sequencing libraries were prepared with the NEBNext® Ultra (TM) II Directional RNA Library Prep Kit for Illumina® Module (New England BioLabs) and NEBNext Poly(A) mRNA Magnetic Isolation Module for mRNA enrichment following the manufacturer's instructions. Libraries were paired end sequenced in an Illumina NextSeq 550 platform with a NextSeq 500/550 High Output 300 cycles kit.

Sequence quality control was performed as in gene panels. Filtered reads were mapped against the human reference GRCh38 genome using STAR v2.7.3a [28]. Isoform quantification was performed with RSEM v1.3.3 [29] and further processed with Tximport v1.16.1 [30] to summarize counts per gene. Differential expression analysis was conducted with the DESeq2 v1.28.1 package [31] using an adjusted p-value cutoff of 0.05.

Gene-set enrichment analysis (GSEA) analysis [18] was run against the hallmark gene sets from the Human Molecular Signatures Database (MSigDB). The significance threshold was set at P-value below 0.05 and nominal false discovery rate (FDR) below 0.05.

To determine the agreement between gene expression in tumor tissues and their derived organoids, first stromal genes described in Isella et al. [32] and non-protein coding genes were removed. To minimize the effect of the artificial environment in which organoids are grown, differentially expressed genes between tumor tissues and organoids were identified with DESeq2 and filtered out. The remaining genes were used to calculate Pearson correlation coefficients between organoids and tissues and were plotted in a heatmap along with an unsupervised hierarchical clustering using the heatmap v1.0.12 R package (https://CRAN.R-project.org/package=pheatmap).

Droplet-digital PCR

DNA from PDOs was digested with the restriction enzymes Mse I (TakaRa, cat. No. 1247A) and Hae III (TakaRa, ca. no. 1051A). PCR droplets were generated using the QX200 droplet generator (Bio-Rad), and the PCR reaction was run in a C1000 Touch thermal cycler (Bio-Rad) according to the manufacturer’s protocol. Results were analyzed with the QuantaSoft software (Bio-Rad). The following probes have been purchased from Bio-Rad: ERBB2 (cat. No. dHsaCP1000116), TP53 (cat. No. dHsaCP1000586), AURKA (cat. No. dHsaCNS193384431), CDKN2A (cat. No. dHsaCP1000581), FGFR1 (cat. No. dHsaCP2500319), MET (cat. No. dHsaCP2500321), SMAD4 (cat. No. dHsaCP2500468). AGO2 was used as a reference assay (cat. No. dHsaCP2500349). DNA from fresh normal mucosa tissue was used as diploid control sample.

Proteomics analysis by LC–MS/MS-SWATH

Sample preparation

PDOs pellets were lysed in UTC buffer (7 M urea, 2 M thiourea and 4% CHAPS) and the protein concentration was measured with the RC DC Protein Assay (Bio-Rad). A pooled sample containing 5 µg of protein from each sample was used to prepare the peptide library. For the SWATH analysis, 20 µg of protein from each sample were used. Samples were denatured at 95 °C for 5 min in sample buffer and subjected to sodium dodecyl sulfate–polyacrylamide gel electrophoresis (SDS-PAGE). Pooled samples were resolved in gel and lanes were cut into five pieces. Gels containing individual samples were not resolved and whole samples were sliced into a single piece. Protein digestion and subsequent analysis were performed as published elsewhere [33] on Proteomics Service (SCSIE) of the University of Valencia.

LC–MS/MS

Spectral peptide library was obtained by liquid chromatography and tandem mass spectrometry (LC–MS/MS) analysis, operating the instrument in a data-dependent acquisition mode. Peptide mixtures were loaded onto a trap column (3 µ C18-CL, 350 μm × 0.5 mm; Eksigent) and desalted with 0.1% TFA at 5 µl/min for 5 min. Peptides were then loaded onto an analytical column (3 µ C18-CL 120 Ᾰ, 0.075 × 150 mm; Eksigent) equilibrated in 5% acetonitrile 0.1% FA (formic acid). Elution was carried out with a linear gradient of 7 a 40% buffer B in A (A: 0.1% FA; B: 5% acetonitrile (ACN), 0.1% FA) at a flow rate of 300 nL/min over 60 min. Eluted peptides were analyzed in a mass spectrometer nanoESI-qQTOF (TripleTOF 6600; SCIEX). Up to 100 ions were selected for fragmentation after each survey scan. Data files were processed using the ProteinPilot search engine (version 5.0, SCIEX) to search the Swiss-Prot database with the following parameters: Trypsin specificity, IAM cys-alkylation, Species Homo sapiens, and the search effort set to through and FDR correction.

Quantitative analysis of individual samples was performed by sequential window acquisition of all theoretical spectra-mass spectrometry (SWATH-MS). Peptides were analyzed by LC using a trap column and analytical column as previously described but operating the TripleTOF 6600 mass spectrometer instrument in SWATH mode. We used 100 variable windows from 400 to 1250 m/z with a total cycle time of 2.79 s. Quantitative data was extracted from.wiff files with Peak View 2.2 (SCIEX) using the peptide library generated by ProteinPilot as indicated above. For every protein in the spectral library, a maximum of 20 peptides were quantified among those with a 95% confidence threshold and an FDR lower than 1%. Shared peptides were excluded. For every peptide, a maximum of 6 transitions (fragment ions) were quantified. All proteins contained in the library were monitored in all samples, producing complete quantitative matrices.

Data analysis

SWATH quantitative data (protein areas) were median normalized and Log2 transformed. Principal component analysis (PCA) was performed with the Orange data mining toolkit to inspect sample grouping and similarity. Differentially expressed proteins were selected by Student t-test using a p-value of 0.05. Functional analysis of differential proteins was performed using STRING database [34] and Cytoscape StringApp [35] within Cystoscape software [36].

Proteotranscriptomic data integration

Differential proteins identified by RNA-seq and proteomics were pooled in a single protein list with unique UniProt accession codes. Functional analysis of differential proteins was performed using STRING database [34] and Cytoscape StringApp [35] within Cystoscape software [36]. To select the most relevant interactions, only high (score 0.7) and highest (score 0.9) confidence interactions were used for oxaliplatin and palbociclib datasets, respectively. A large network of interacting proteins was observed in each dataset containing about 30 (oxaliplatin) and 45 (palbociclib) % of proteins. The networks were further characterized by performing a functional enrichment analysis of the main modules observed. In the case of palbociclib dataset, a part of the network, containing mainly metabolic proteins, has a more diffuse appearance. In order to gain some information regarding the main metabolic pathways included in this part of the network, their proteins were clustered using the MCL algorithm within StringApp using a Granularity parameter of 2. A functional enrichment analysis was performed on the resulting clusters containing 4 or more proteins. In addition, the metabolic proteins were mapped into KEGG Metabolic pathways (hsa01100).

Western blot analysis

To verify some of the differentially expressed proteins only detected by proteomics, we performed a Western Blot analysis. We followed a previously published protocol [37]. Briefly, whole-cell protein extracts were prepared using RIPA buffer (50 mmol/L Tris–HCl pH 7.5, 150 mmol/L NaCl, 0.1% SDS, 1% Triton x-100, 0.5% deoxycholic acid sodium salt (w/v)) supplemented with 2 μL/mL protease inhibitor cocktail (Sigma) and 10 μL/mL phosphatase inhibitor cocktail (Sigma). Samples were sonicated and centrifuged at 14,000 × g speed for 20 min at 4 °C. The protein concentration was quantified by BCA protein assay kit (Thermo Scientific). 20ug of total protein were loaded per well on 12% SDS-PAGE, transferred onto a nitrocellulose membrane and incubated with the following primary antibodies from Santa Cruz: anti-TCP1 (Cat. No. sc-374088), anti-CCT2 (Cat. No. sc-374152) and anti-CCT8 (Cat. No. sc-377261). Immunoreactive bands were recognized using peroxidase-conjugated secondary antibody (DAKO). Immunoblots were visualized using the “ECL Western Blotting detection kit” reagent (GE Healthcare) and the ImageQuant LAAS 400 (Healthcare Bio-Sciences) system. Quantification of detected bands was carried out using ImageJ.

Drug assays

PDOs were trypsinized till single cells, plated in 50% BME domes with complete medium in 96-well plates and 48 h later, after having observed the formation of organoids, treated with increasing doses of different drugs. 5-Fluorouracil and oxaliplatin were kindly provided by HCUV Pharmacy. The following drugs were purchased from Selleck: SN-38 (cat. No. S4908), erlotinib (cat. No. S7786), trametinib (cat. No. S2673), tepotinib (cat. No. S7067), erdafitinib (cat. No. S8401), TAS-102 (cat. No. S8539), alpelisib (cat. No. S2814), palbociclib (cat. No. S1579), olaparib (cat. No. S1060), GSK458 (cat. No. S2658), birabresib (cat. No. S7360), AZD6738 (cat. No. S7693), regorafenib (cat. No. S1178). After 120 h of treatments, CellTiterGlo 3D (Promega, cat. No. G9681) viability assay was performed following the manufacturer’s instructions. All experiments were performed with three technical replicates and with at least two independent biological experiments. Area under the curve (AUC) was calculated with trapezoid rule. Statistical analysis and graphical representation were performed with GraphPad Prism software v. 8.2.1. Experimental replicates of drug screening were compared with one-way ANOVA assuming a p < 0.05 as statistically significant and represented as linear regression. For each drug a pattern was assembled containing the corresponding Ln-AUCs across the treated organoids. Next, Z-scores were calculated for each drug, considering median Ln-AUC and standard deviation for each point (LnAUC-mean/standard deviation). Finally, the patterns were aggregated column-wise into a matrix. The obtained matrix was used to assess the relative sensitivity/resistance of each organoid line.

Results

PDOs can be established from pre-treated advanced colorectal cancer patient samples

A total of 34 patients with advanced CRC were included in the present study between November 2018 and November 2020, and 50 samples were collected, mainly from liver metastasis. The primary tumor was biopsied only when metastases were not accessible. In the presence of more than one liver metastasis, a sample was taken from each one. We included both naïve and heavily pre-treated patients. Organoids’ growth was observed in 29 lines from 22 patients, with an overall success rate of 59%. All sample characteristics are shown in Supplementary Table S1. Most biopsies were performed at surgery. PDOs establishment was successful in 90% of cases without chemotherapy during the previous 6 months, dropping to 51% if preoperative chemotherapy had been administered, although the difference was not statistically significant. The PDOs’ establishment rate was not affected by primary tumor location (PTL, left vs. right vs. rectum), tissue biopsy site (primary vs. metastasis), or RAS/RAF mutational status. In 19 samples from 12 patients no organoids could be grown, this culture failure owing mainly to lack of initial growth, arrested growth or initial contamination. We observed heterogeneous growth behavior among our PDOs. Some could be established as long-term cultures (more than 3 months in culture and over 10 passages) while some others could only grow as short-term cultures (1–3 months in culture and 1–9 passages) [38] (Supplementary Fig. S1). Long-term cultures displayed an exceptional capacity for multiple freeze thaw cycles. No specific molecular features were identified as predicting long-term versus short-term culture establishment.

PDOs recapitulate morphology and immunohistochemistry and could be derived from tissues with low cellularity

PDOs were stained with hematoxylin and eosin (H&E) to show cellular architecture, which faithfully reproduced the morphology seen in culture (Fig. 1A, top and middle panel). Cellular architecture presents similarities with corresponding tissues (Fig. 1A, bottom panel).

Fig. 1
figure 1

PDOs recapitulate morphologic features of original tissue and can be obtained from low cellularity biopsies. (A) Culture images of some PDOs lines (top panel; scale bar 50 μm, 20X), compared with corresponding H&E staining (hematoxylin and eosin) (middle panel; scale bar 50 μm, 20X) and with matched tissue stained with H&E (bottom panel; scale bar 100 μm, 10X). (B) Corresponding percentage of tumor cells assessed as percentage of neoplastic cells with respect to the total amount of viable cells in the tissue sample from which the culture was derived, assessed by the pathologist. (C) IHC for CDX2 and CK20 antibodies in a selection of PDOs. CNO75 is a normal mucosa organoid from a CRC patient, used as control for both markers. (Scale bar 50 μm, 20X). (D) Lack of expression of MLH1 and PMS2 in a core biopsy from a liver metastasis of patient 50 and its corresponding derived PDOs. RTO2 is a pMMR model used as positive control. (Scale bar 50 μm, 20X). T: tissue

To elucidate whether tumor cell percentage could impact on the capacity to obtain organoids, each tissue section destined for organoid culture was revised by a dedicated pathologist. A cut-off of 30% tumor cells was applied to classify our samples. Our data indicate that PDOs can be grown independently of tumor cell proportion found in the original tissue (two-tailed Fisher’s exact test p = 0.4568, Fig. 1B).

IHC staining shows that PDOs are positive for CDX2 and CK20, intestinal markers that are both employed in standard diagnosis of CRC (Fig. 1C). The only exception was seen in mCTO66S3, which tested CK20 negative. The original tissue displayed weak and parceled expression, with most cells clearly negative (Supplementary Fig. S2). Indeed, loss of CK20 expression may be considered a negative prognostic marker in some settings [39]. CDX2 and CK20 expression of corresponding tissues is showed in Supplementary Fig. S2, where a weak expression of CK20 and no CDX2 staining in mCTT47 was observed, probably due to low cellularity, high level of fibrosis and intratumoral heterogeneity. On the other hand, our PDO mCTO50 and mCTO50B, derived from a patient with MSI high showed lack of expression of MLH1 and PMS2 proteins, as was seen in the original tissue (Fig. 1D, Supplementary Fig. S2).

PDOs recapitulate the genomic and transcriptomic profile of the original tumor

To genomically compare organoids and their original tissues, we employed NGS analysis using a customized panel [19]. A high correlation was found regarding germline variants, single nucleotide variants (SNVs), short insertions and deletions (r > 0.8, p < 0.0001, Supplementary Fig. S3A-B). The percentage of tumor cells, necrosis, mucinous proportion, and tumor stroma ratio were calculated by a dedicated pathologist. These indices did not interfere with the ability of PDOs to reproduce the original tissue genomic profile (Supplementary Fig. S3C), except when low cellularity was present (Supplementary Fig. S3B).

As expected, the most prevalent mutations affected TP53, KRAS and APC (Fig. 2A), and besides these, PIK3CA, PTEN and SMAD4 mutations were also detected, encompassing all relevant mutations in CRC. The main driver oncogenic mutations and hotspot variants are represented in Table 1. In most cases, allelic fraction (AF) was higher in PDOs (paired t test, p < 0.0001), reflecting their enrichment in epithelial cells. Moreover, CytoScan HD SNP-array was conducted on PDOs and matched fresh tissue to detect copy number alterations and LOH. Likewise, PDOs recapitulate the overall copy number variation profile, and are significantly enriched (Fig. 2B). CytoScan results were validated by ddPCR (Supplementary Fig. S3D).

Fig. 2
figure 2

PDOs recapitulate genomic and transcriptomic profile of original tissues. (A) Concordance of SNVs and InDels between PDOs and corresponding tissues. On the left percentage of genomic alterations detected across the study subjects. (B) Whole Genome View representation of long-term PDOs and corresponding tissues according to Cytoscan HD. Data is expressed as the weighted log2 ratio of the copy number on the left Y-axis, and the chromosome number on the X-axis. 0 corresponds to diploid, upper and lower spikes indicate gain and loss regions, respectively. mCTO50 tissue was not available. (C) Heatmap showing the Pearson correlation coefficient (color key) between tumors (rows) and organoids (columns) based on the normalized counts as described in materials and methods section. Dendrograms show the hierarchical clustering based on the complete method and Euclidean pair-wise distance. Different passages from the same organoid culture are included

Table 1 Allelic fraction (AF) distribution in PDOs and matched tissues of main driver, hotspot and not hotspot mutations with potential relevance. Each mutation has been manually reviewed in Integrative Genome Viewer (IGV). Only mutations with an AF of at least 5% in PDOs have been included

RNA-seq was also performed to compare gene expression between PDOs and corresponding tissue and to analyze whether it was stable over time in culture. This is a key issue, considering PDOs as in vitro “avatars” of patient tumors. Three biological replicates deriving from distinct culture passages were employed for each organoid line, except for RTO2 and RTO7 for which two replicates were employed, due to technical issues. Overall, PDOs and tissues cluster separately because of the co-existence of more cell types in tissues and stromal genes (Supplementary Fig. S3E). However, narrowing down the analysis to non-stromal genes contribution (supplementary methods), significant expression correlation can be observed between each PDOs line and its corresponding tissue (Fig. 2C). The two exceptions were mCTO50B and mCTO47, which were derived from a tissue with low cellularity (less than 10%). Interestingly, PDOs from different passages exhibit a similar gene expression profile, indicating that the influence of culture condition is minimal and that this profile is stable over time.

PDOs show differential anti-tumor response to standard and experimental treatments

To study whether PDOs could be an appropriate model to assess differential response to standard therapies, we exposed our long-term models to several approved drugs. The reproducibility of drug treatment results in PDOs was confirmed across several experiments (Supplementary Fig. S4A). Table 2 summarizes the clinical and molecular characteristics of patients whose models were used.

Table 2 Clinical and molecular characteristics of patients from which long-term cultures have been generated

All the established PDOs exhibited a heterogeneous response to chemotherapy (Fig. 3A). As an example, mCTO50B (a PDO obtained after chemotherapy) had a worse response than mCTO50, which was generated from the initial tumor before any therapy was given, so the former could represent a potential model of chemo-resistant disease.

Fig. 3
figure 3

PDOs response to drug screening. (A) Log transformed dose–response curves in selected standard drugs and non-standard (C) drugs. (B) Z-score Ln-AUCs heatmap (red: no response; blue: good response) for standard treatment. (D) Z-score Ln-AUCs heatmap (red: no response; blue: good response) for non-standard treatment (lower panel), matched with genomic data. In the upper panel are depicted main gene mutations, in the following loss of heterozygosity (LOH) and in the last copy number variations

SN-38 (the active metabolite of irinotecan) was highly active in vitro, and only RTO2 displayed lower sensitivity compared to the other lines in our panel. Adding oxaliplatin to 5-FU resulted in a general decrease in PDOs viability (Supplementary Fig. S4B). Along these lines, combination with SN-38 seemed to be more effective than single agent exposure (Supplementary Fig. S4B, not statistically significant). Anti-EGFR treatment with erlotinib showed modest activity in all lines, including the RAS wild type RTO7 PDO. Regorafenib also showed modest in vitro activity, whereas trifluridine/tipiracil (TAS-102) was more active in all PDOs except mCTO66S3 (Supplementary Fig. S4C). A matrix representation of normalized Z-score Ln-AUCs (Fig. 3B) identifies relative resistance/sensitivity to a single drug among different PDOs.

To further explore sensitivity to other anti-cancer agents, our PDOs were exposed to different drugs not currently included in the standard of care for advanced CRC. Agents targeting pathways altered in our models were selected. Response was matched with genomic alterations in all models (Fig. 3C-D). Differential responses were seen for palbociclib, trametinib, alpelisib, the BET inhibitor birabresib (Fig. 3C). For example, RTO7 showed a dramatic response to palbociclib, a selective inhibitor of CDK4/6 kinases, in comparison with the other PDOs, particularly RTO2 and mCTO66S3. Considering genomic data, the lower responder mCTO66S3 has the higher copy number level of both CDK4 and CDK6, reported as a mechanism of resistance [40]. Additionally, RTO7 has the highest copy number of c-MYC, potentially indicating cell cycle activation [41,42,43]. Moreover, this line has a SMAD4 p.P356S missense variant, in a hotspot region with LOH, thus explaining the nearly 100% mutant allelic fraction (Table 1). Despite the lack of functional validation, it is predicted to be oncogenic by several databases, leading to a loss-of-function of the protein [44] and is associated with increased c-MYC activity [45]. Indeed, the high response to palbociclib observed in RTO7 certainly warrants further research, as all these data point to a “MYC-addicted” phenotype of RTO7.

The BET inhibitor birabresib (Fig. 3C, middle) showed good activity in all lines, particularly in SMAD4 loss RTO7 and mCTO66S3 (Fig. 3D, Supplementary Fig. S4D, S5A), confirming synthetic lethality via restoration of MYC inhibition [45]. MEK-inhibitor trametinib showed modest activity in all lines except for KRAS mutant RTO2 (Fig. 3C), while PI3Kα-inhibitor alpelisib showed modest activity in all PDOs except for a slightly more in PIK3CA mutant line RTO2 (Fig. 3C, upper right). Intriguingly, the PI3K-mTORC1/2 inhibitor GSK458 was significantly active, not only in PIK3CA mutant RTO2, but also in mCTO66S3 (Fig. 3C, lower right). In contrast, mCTO50B had a significantly worse response. Indeed, this PDO has a frameshift mutation of the PTEN gene (p.K267Rfs*9) reported as likely oncogenic in OncoKB™, and we showed that it leads to reduced gene expression and a loss of PTEN expression per IHC (Supplementary Fig. S5B). We also performed additional drug testing without observing relevant activity (Supplementary Fig. S4D).

PDOs differentially replicate patient response depending on the site of disease

To determine whether our PDOs would also reproduce patient response to treatment, we retrieved clinical information and matched drug response observed for mCTO50 and RTO7 (Fig. 4), in which follow up data events were available. mCTO50 was generated from liver metastasis of a patient with MSI-high CRC before initiating any anti-tumor therapy. This patient received cytoreductive chemotherapy with FOLFOXIRI (a combination of 5-FU, leucovorin, oxaliplatin and irinotecan), presenting appropriate tumor shrinkage and allowing liver metastases resection. When tumor response was tested in vitro, mCTO50 was significantly more sensitive to 5-FU and oxaliplatin than mCTO50B, which was generated from the remaining residual disease tissue obtained at surgery. The drug treatment effect in this PDO was comparable to that observed in the patient as mCTO50. The patient had disease progression shortly after surgery, indicating that chemotherapy failed to exert a long-term effect. A similar degree of sensitivity to irinotecan was observed in both PDOs, however, prompting us to hypothesize that resistance mechanisms to this drug were not yet established.

Fig. 4
figure 4

Patient’s response matched with PDOs drug assays, RTO7 in the upper panel and mCTO50 in the lower panel. Each CT scan response evaluation is compared with corresponding treatment administered in the corresponding PDOs. ≠ means discordant response, = means concordant response. Mo(s): months in terms of progression-free survival; St: stage; FOLFIRI: 5-fluorouracil and irinotecan; FOLFOX: 5-fluorouracil and oxaliplatin; FOLFOXIRI: 5-fluorouracil, oxaliplatin and irinotecan; BSC: best supportive care; SD: stable disease; PR: partial response; PD: progressive disease; R1: microscopic residual tumor

RTO7 exhibited a good response to irinotecan-based chemotherapy, received as part of the patient’s first line treatment regimen (Table 2). Conversely, patient response was poor (showing progression at first CT scan). The only concordant response observed was for the anti-EGFR agent and for TAS-102 (i.e. lack of response). However, it should be taken into consideration that the organoid line derived from the primary location rather than liver metastases, which could not be biopsied. This highlights the intrinsic biological heterogeneity of CRC and the need to biopsy metastatic sites in the context of translational studies.

Proteotranscriptomic analysis in PDOs in relation to oxaliplatin or palbociclib sensitivity

Our aim was to uncover functional correlates between baseline proteotranscriptomic expression and drug sensitivity in PDOs. To assess the steady-state protein abundance, a SWATH-MS based label-free quantitative proteomics analysis was performed as a more reproducible strategy in a translational context. The relative level of the 1157 proteins quantified in all samples (FDR < 1%) helped us to distinguish each PDO by principal component analysis (PCA) (Supplementary Fig. S5C) and these were consistently reproduced at different culture passages. We focused on understanding the molecular mechanisms involved in sensitivity to palbociclib and the lack of response to oxaliplatin, where the most striking differences were observed among PDOs. Three biological replicates proceeding from different culture passages were used for each PDOs line.

RTO2, RTO7 and mCTO50 showed better response to oxaliplatin than mCTO50B and mCTO66S3 (unpaired t-test p = 0.0038) (Fig. 3A). A proteomics analysis (Supplementary file S1) identified 95 differentially expressed proteins (51 up-regulated, 44 down-regulated) in non-responder compared to responder organoids. Among the proteins identified with the greatest positive and negative fold change (Fig. 5A), we found several previously associated with oxaliplatin resistance. As an example, resistant PDOs showed higher levels of PRDX6, which is a negative regulator of ferroptotic cell death [46], a process that enhances CRC sensitivity to oxaliplatin [47]. Another up-regulated protein was ALDH9A1, a member of the ALDH family of proteins which are involved in aldehyde detoxification which in turn is associated with acquired chemoresistance in colorectal cancer cells [48]. These same models exhibited lower levels of NDRG1, a replication stress response protein which can inhibit epithelial-mesenchymal transition (EMT), a process that has been associated with phenotypes of chemoresistance [49], and of CDH17, which has been reported as a marker of good response to 5-fluorouracil and oxaliplatin-based chemotherapy [50].

Fig. 5
figure 5

Proteotranscriptomic characterization of oxaliplatin responder lines in comparison with no-responder ones. (A) Network mapping of 95 proteins showing the 4 sub-clusters composed of more than 3 proteins. Relevant proteins with no associations to others are represented as isolated nodes. Colors are depicted according to the protein abundance (log2ratio) compared to responder PDOs (left panel). Bar-chart of GO terms represented as percentage of annotated proteins using the same color coding (right panel). (B) Hierarchical clustering heatmap of differentially expressed genes per transcriptomic. (C) GSEA hallmarks analysis of RNA-seq data. NES: normalized enrichment score. FDR: false discovery rate. R: responder; NR: non responder

To gain insight into alternative mechanisms of oxaliplatin resistance, a functional protein–protein interaction network of differentially expressed proteins was evaluated by STRING database (http://string-db.org, version 11). Using a high confidence score, some interesting interactions were identified, with some subclusters remaining at the highest confidence score (Fig. 5A). Gene ontology (GO) analysis revealed an increase of proteins related with translation in non-responder models, with an interesting subcluster enriched in tRNA-aminoacylation process that included different aminoacyl-tRNA synthetase (ARSs) and an auxiliary protein AIMP2, joined to form the cytoplasmic multi-tRNA synthetase complex [MSC] [51]. Beyond their role in protein synthesis, ARSs were related to the induction of unfolded protein response, which in turn was associated with escape from chemotherapy induced senescence [52].

Another enriched process in non-responder models was mitochondrial import, including among other proteins two subunits of the ATP synthase (ATP5A and ATP5B), a mitochondrial membrane protein complex mediating ATP synthesis, and citrate synthase (CS), one of the key enzymes in the tricarboxylic acid cycle (TCA), indicating a switch from a glycolytic based to a mitochondrial metabolism with oxidative phosphorylation (OXPHOS) as main source of energy. This could represent a way for some cancer cells to repair platinum-induced DNA damage more efficiently [53].

In the RNA-seq analysis, 597 differentially expressed genes were detected (Fig. 5B). GSEA study showed that non-responding PDOs were enriched in G2M checkpoint, TGFbeta and DNA repair hallmarks (Fig. 5C), somewhat expected as regards oxaliplatin response. Indeed, it acts by inducing DNA adducts, therefore an increase in the capacity to repair damaged DNA could help the cancer cell to survive. A key moment of DNA repair occurs during G2 to M phase transition. In addition, TGFbeta contributes to oxaliplatin resistance in CRC [54]. We noted that non-responder PDOs showed enriched unfolded protein response and PI3K-Akt-mTOR and mTORC1 signaling hallmarks (Fig. 5C). As previously indicated, the former could be involved in the escape from chemotherapy induced senescence [58].

After observing high sensitivity to palbociclib in RTO7 (unpaired t-test p = 0.0002), and to understand the molecular mechanisms involved in this, a comparison of protein expression profiles was made, revealing 245 significantly changing proteins in RTO7 versus non-responder models (RTO2 and mCTO66S3) (Supplementary file S1). Functional protein–protein interaction networks were evaluated using the STRING database with a highest confident score (Fig. 6A). Gene ontology (GO) functional enrichment analysis matched the largest identified cluster (55 of 111 proteins) to biological processes related to protein synthesis, folding and degradation, and with mRNA splicing process via spliceosome, with most of the proteins upregulated in our RTO7 PDO as a palbociclib- responder model (Supplementary Fig. S5). We identified a subcluster comprising all the eight proteins that form the T-complex protein Ring Complex (TRiC) (Supplementary Fig. S5D), an essential eukaryotic molecular chaperonin that aids in the folding of ~ 10% of the proteome including oncoproteins [55].

Fig. 6
figure 6

Proteotranscriptomic characterization of Palbociclib responder line in comparison with no-responder ones. (A) Network mapping of 111 proteins showing the 4 clusters composed of more than 3 proteins. Proteins with no associations to others were removed and nodes are colored according to the protein abundance (log2ratio) compared to responder PDO (left panel). Bar-chart of GO terms represented as percentage of annotated proteins using the same color coding (right panel). (B) Hierarchical clustering heatmap of differentially expressed genes per transcriptomic. (C) GSEA hallmarks analysis of RNA-seq data. NES: normalized enrichment score. FDR: false discovery rate. R: responder; NR: non responder

Other top processes included in the remaining clusters were related to regulation of cellular component organization, cell adhesion, and metabolic processes, among others (Fig. 6A). All these processes point towards a high proliferative phenotype supported by the up regulation of different mechanisms that preserve proteome integrity.

Transcriptional-wide changes were also observed in palbociclib sensitive model (Fig. 6B). GSEA analysis showed that RTO7 is enriched in MYC targets hallmarks, fatty acids and lipid biosynthesis (which together are a consequence of MYC activation [56]) and unfolded protein response (Fig. 6C), a process that indicates an alteration of protein homeostasis. In the latter hallmark, NPM1, a c-MYC activator [57], is one of the most significantly enriched genes. We could not quantify Myc protein, by SWATH-MS analysis, maybe due to relative scarcity being it a transcription factor. Nevertheless, we were able to detect higher levels of NPM1 among non-clustering proteins.

We observed a weak correlation between protein abundance and transcript quantification among differentially expressed proteins in both comparisons (palbociclib: r = 0.1477, p < 0.01; oxaliplatin: r = 0.1446, p = 0.02; Supplementary Fig. S6A), indicating the occurrence of complex post-transcriptional regulation. For instance, the relevant proteomic data regarding TRiC and ARSs roles in palbociclib responder and oxaliplatin non-responders’ models respectively, were not confirmed by RNA-seq (Supplementary Figure S6B). Nevertheless, some relevant processes related to protein folding, biosynthesis and proliferative features were captured by both RNA and protein analysis.

Integrative functional network (IFNA) analysis of proteotranscriptomic data

To gain deeper insight into the processes involved in the differential drug response in PDOs we performed an integrative analytical approach [58]. Thus, we fused the transcriptomics and proteomics differentially expressed datasets and extracted the common functional context through an integrative network analysis by STRING application via Cytoscape.

Data fusion and integrative functional network analysis (IFNA) were conducted for both oxaliplatin and palbociclib drug response comparisons. The functional clusters shown in the integrative network, representing a high level of interaction and integration of the proteomic and transcriptomic data, are referred as modules in the rest of the manuscript. High confidence interacting proteins are shown in Fig. 7 and Fig. 8, where relevant modules are marked with a shading. In both cases, large networks containing approximately 30–45% of the differential proteins were observed. Interestingly, most of the proteomic clusters and RNA-seq hallmarks (Fig. 5B and Fig. 6B) are now contained in the networks, including some previously isolated clusters. The functional enrichment analysis of the different modules observed in IFNA are listed in Additional File 2 for oxaliplatin and palbociclib comparisons.

Fig. 7
figure 7

Integrative functional network analysis of oxaliplatin RNA/protein fused dataset. High confidence interactions (score 0.7) showing 7 modules. Nodes in red correspond to proteins identified as differentially expressed by RNA-seq only, in blue those identified by proteomics only and in green those identified by both omics. Key for modules: 1: Glutathione metabolism & antioxidant activity; 2: Nitrogen compound metabolism; 3: Translation & tRNA aminoacylation; 4: TCA & Oxidative phosphorylation; 5: Cell cycle regulation; 6: Antigen presentation; 7: Cell adhesion

Fig. 8
figure 8

Integrative functional network analysis of palbociclib RNA/protein fused dataset. Highest confidence interactions (score 0.9) showing modules. Nodes in red correspond to proteins identified as differentially expressed by RNA-seq only, in blue those identified by proteomics only and in green those identified by both omics. Grey line marks on the left 7 metabolism related modules with the highest integration. Key for metabolic modules: 1: Oxidative phosphorylation; 2: Pentose phosphate pathway, and Glycolysis; 3: Metabolism of nucleotides, and purinergic nucleotide receptor signaling pathway; 4: PPARA signaling pathway; 7: Lipid biosynthetic process; 9: Oxidoreductase activity; 10: Organonitrogen compound metabolic process. Key for modules on the right: 1: Regulation of signal transduction; 2: Regulation of cell communication and motility; 3: Translation; 4: RNA splicing; 5: Transcription and DNA repair; 6: Cell cycle; 7: Membrane trafficking; 8: Proteasome; 9: Protein ubiquitination; 10: Interferon signaling; 11: ER-Golgi vesicle-mediated transport; 12: TRiC/CCT complex. Details of the modules are indicated in Additional file 2

In the context of oxaliplatin response, the whole dataset contained 12 proteins commonly represented by both omics, 83 only by proteomics and the remaining by RNA-seq. Only high confidence interactions are shown in Fig. 7 containing the main network and two isolated clusters. Globally it contained 191 proteins, 50 of them were detected only by protein, 136 only by RNA while 5 were commonly detected.

Within the main network, 5 modules were selected to perform a functional enrichment analysis that confirmed the processes highlighted in individual omics profiles and now appear connected into this new integrative network. For instance, the relevant translation and tRNA aminoacylation processes appear in module 3 which is now complemented with RNA data. Another interesting finding is module 4, where the proteomic cluster mitochondria protein import is now complemented with TCA cycle and Oxidative phosphorylation processes by RNA data. Proteasome proteomic cluster and DNA repair RNA-seq hallmark are now enriched processes in module 5.

Some previously unconnected differential proteins now appear integrated into the network. Such is the case of CDH17 in module 7, and, interestingly, PRDX6 that is now part of module 1 enriched in redox metabolism, which may be involved in oxaliplatin resistance.

Regarding palbociclib sensitivity, the fused dataset contained 172 proteins detected by proteomics, 73 detected by both proteomic and RNA-seq, and the remaining detected only by RNA-seq. Figure 8 shows the large network of the highest confidence interacting proteins in palbociclib comparison. Certain modules, related to splicing, translation, proteostasis and metabolism showed a high complementarity between RNA-seq and proteomics. This complementarity reinforces the functional enrichment already defined. Interestingly, TRiC complex, that was only detected as differential at protein level (some of their subunits verified by western blot analysis, Supplementary Fig. S5E), is now complemented by BAG5, involved in endoplasmic reticulum stress regulation and unfolded protein response activation [59]. In fact, it has been described that the protein BAG5 is a co-chaperone involved in protein folding [60].

A new central module related to cell cycle and mainly represented by RNA-seq data, appeared highly connected with the proteasome, TRiC, and (to a lesser extent) with splicing processes. Moreover, the genes that define the GSEA Myc targets signature, a hallmark enriched by RNA-seq analysis, map into this network encompassing all these processes.

The large diffuse metabolism module constitutes a good example of differential RNA and protein integration. Applying a clustering algorithm as described in Methods section, we found a clear enrichment in the pentose phosphates and glycolysis pathways, and also in lipid and nucleotide metabolism in relation to palbociclib sensitivity. When we mapped differentially expressed metabolic genes into KEGG Metabolic pathway map, we observed a nice complementation of both omics (Supplementary Figure S7). These cellular pathways can satisfy a high energy and anabolic precursors demand to sustain enhanced growth and proliferation in this putative Myc-addictive phenotype.

Discussion

Advanced CRC is a complex and heterogeneous disease in which resistance development limits the success of drug treatments. PDOs have been shown to maintain the same genomic and phenotypic features of original tumor tissue and can be used for high-throughput drug screening, enabling us to anticipate clinical response. Although the number of patients included in this study is not high, we demonstrate the potential of integrating PDO drug response with a deep proteotranscriptomic analysis, providing novel insights into sensitivity to oxaliplatin and palbociclib.

We report the generation and characterization of a library of PDOs from advanced CRC patients, including both naïve and pretreated subjects, showing that it is possible to obtain organoids even from heavily pretreated samples. We were able to generate a basal and progressive MSI-high PDO, to our knowledge the first study to achieve these results. MSI lines were generally derived for primary CRC tumors [11].

Our organoids recapitulate their matched tumor tissue at genomic and transcriptomic level. We observed an increase in PDO variant allele frequencies (VAFs) and copy number variations compared to corresponding tissues, due to the enriched epithelial nature of organoids, avoiding the dilution effect of the DNA of other cell types such as fibroblast, endothelial and immune cells, etc. that are present in the tumor microenvironment (TME). Therefore, we could detect rare genomic events, such as the previously characterized SMAD4 and PTEN mutations (Table 1), uncovering their potential oncogenic relevance (PTEN mutation, for instance, leads to protein expression loss). A few reports have been published on the comparison of transcriptomic similarity between PDOs and original tissues [61] and direct comparison showed a clear difference between organoids and tissues, as expected. Nevertheless, after adjusting for stromal gene contribution, organoids displayed a highly similar expression profile compared with that of tissue, which is stable over time. This proves an only minimal effect of culture condition, and that expression findings are robust and result from cancer cells.

Regarding tissue amount required for organoid generation, we demonstrated the feasibility of deriving reliable cultures even from tissues with a very low cellularity, although the genomic and phenotypic features of these PDOs may show lower correlation with the original tissue. Demonstrating that these models could reproduce and therefore predict patient outcomes is crucial for their use in precision medicine. Indeed, our mCTO50B PDO derived from extremely low cellularity tissue could effectively replicate patient response. Biopsy site could be viewed as a more important issue. Given that not all lesions are easily accessible for biopsy, the use of a primary site as a substitute should be considered with caution as it might not replicate ex vivo the response observed in the clinic. This was the case with RTO7, which could only be used a primary site biopsy, while the response observed in the PDOs failed to perfectly match what it was observed in liver metastasis.

We captured heterogeneous drug response to several agents and performed a proteotranscriptomic analysis of responder and non-responder models to identify baseline molecular features associated with response. This is particularly important for backbone chemotherapeutic agents such as oxaliplatin, for which no biomarkers have been brought into clinical use yet [62]. SWATH-MS was selected for proteomic study as an emerging strategy for biomarker discovery due to the coverage achieved, the reliable quantitation obtained and the possibility to re-interrogate the data [63]. Few studies have characterized organoids at the quantitative proteomic level. We report on the use of SWATH in PDOs in this regard, underlining its potential utility to identify a bona fide biomarker of response to anticancer agents, in contrast to most current efforts focused on RNA signatures. Indeed, when individually analyzed, we observed cases where both techniques go in the same direction as well as discrepancies between protein abundance and their corresponding RNAs, probably due to mRNA post-transcriptional regulation. However, transcriptomic experiments have a much wider genome coverage than proteomics. So, we consider that both methodologies are complementary and that a multi-omics approach can provide richer biological information for a better comprehension of the complex pathological molecular mechanisms involved in drug response. In addition, after fusing RNA and protein data into a unique dataset and thanks to IFNA, we could confirm those processes which were detected at RNA or protein level and found new modules, enriched in processes which could have a putative role in differential drug response and that were highlighted after dataset integration.

Some of the differential proteins identified in relation to oxaliplatin response have been previously associated with resistance, but surprisingly, functional enrichment analysis highlights the t-RNA aminoacylation process, including some aminoacyl-tRNA synthetases (ARSs) and an auxiliary protein AIMP2 which forms the MSC complex [51]. ARSs are traditionally considered housekeeping molecules since they catalyze the aminoacylation of tRNAs in protein synthesis, an essential process for maintaining cell homeostasis. Nevertheless, ARSs and AIMPs are closely associated with tumor biology [52]. Moreover, according to recent findings in the literature [64], specific ARSs are involved in the induction of unfolded protein response, which has been associated with escape from chemotherapy-induced senescence. This could therefore represent a novel mechanism of resistance to chemotherapy in our models [65]. As a proof of concept, the unfolded protein response process was also captured in the transcriptomic study as one of the most significant hallmarks, and fused data network analysis allowed us to confirm the putative involvement of these processes into a differential oxaliplatin response. In addition, IFNA showed that TCA and OXPHOS genes cluster together with the mitochondrial protein import, a process which was already identified by proteomic approach. Indeed, increased mitochondrial import of proteins involved in ATP synthesis such as ATP5A and ATP5B, together with augmented citrate synthase, a key enzyme of TCA cycle, could mark a shift in tumor metabolism from glycolysis to oxidative phosphorylation pathway (OXPHOS) that could represent a mechanism of oxaliplatin resistance, such as has been unraveled in other tumor models [53]. We could hypothesize that up-regulation of these energy supplier processes guarantees efficient repair of oxaliplatin-induced DNA-damage, since the enzymes involved in this repair consume large amounts of ATP. In fact, targeting ATP synthases has been proposed as a therapeutic strategy to defuel cancer growth. ATP synthase subunit ATP5B has recently been identified as the specific target of apoptolidin A, a glycomacrolide that selectively addressed OXPHOS-dependent cancer [66]. Moreover, data fusion and IFNA evidenced that proteasome proteomic cluster is integrated into a module that contains also genes involved in DNA repair, and highlighted a functional enrichment related to oxidative metabolism and detoxification that connects a previously isolated protein such as PRDX6, a ferroptosis negative regulator [46].

Not all these processes or expression of specific genes were detected by both proteomic and transcriptomic analysis. Two hallmarks, unfolded protein response and PI3K-AKT-mTOR and mTORC1 signaling, were detected by transcriptomic data, but not at the proteomic level. Both were shown to play a significant role in chemotherapy-induced senescence escape [65] and concur with the process highlighted by proteomic results. All this data showed that although transcriptomic and proteomic data were not completely correlated, integrating both analyses could be effective to identify a bona fide mechanism of drug response.

Although CDK4/6 inhibitor palbociclib is not currently used in the clinic to treat CRC, we used it as a proof-of principle of targeted therapy to be assessed in our PDOs. RAS wild type RTO7 showed resistance to anti-EGFR drugs, but was extremely sensitive to palbociclib, and differential gene expression analysis showed that it harbors significant MYC activation (GSEA analysis). Analyzing proteomic data, RTO7 has higher levels of Nucleophosmin (NPM), a c-MYC activator [57]. Interestingly, palbociclib has been reported to suppress NPM phosphorylation at threonine 199, thereby reducing cell proliferation in preclinical models of endometrial cancer [67]. MYC expression has been shown to correlate with poor prognosis in RAS wild type CRC [68], but probably this is not sufficient to justify response to palbociclib. SWATH studies indicated that this model has a strong activation of several processes involved in proteostasis, among which the whole protein folding TCP-1 ring complex (TRiC) was up-regulated in palbociclib-sensitive PDO. TRiC altered expression was not detected at the mRNA level. However, the unfolded protein response process was one of the most significantly enriched hallmarks found in transcriptomics. Moreover, both omics indicate the same process, while in this case proteomic data also pointed more precisely at the presence of this molecular mechanism. TRiC aids in folding approximately 10% of the proteome, with cytoskeletal proteins actin and tubulin among its best-known substrates, which also include several oncoproteins such as P53 and MYC itself [69]. TRiC proteins were shown to promote cell cycle progression and an invasive phenotype [55] and their overexpression was related with poor prognosis in CRC [70]. Taking all of this into account, altered TRiC complex expression may have prognostic relevance and could be a potential new therapeutic target. In fact, highly proliferative cancer cells increase their dependence on several stress response pathways, including the heat-shock response to alleviate competition among proteins for access to chaperones [71]. These results are reinforced and highlighted after data fusion and IFNA. Indeed, this integrative approach showed that BAG5, detected by RNA-seq, appears to strictly interact with TRiC complex. BAG5 is involved in endoplasmic reticulum stress regulation and unfolded protein response activation [60]. In particular, BAG5 interacts with heat-shock proteins acting as a cochaperone in protein folding [72]. Moreover, some TRiC subunits (TCP1 and CCT5) are among Myc Targets v.1 signature genes (GSEA) and, indeed, TRiC complex overexpression has been correlated with Myc activation [73].

IFNA also highlighted two interesting findings that are in line with a putative Myc addicted phenotype. On one hand, a new module appeared corresponding to cell cycle, mainly represented by RNA-seq data, which was strictly connected to modules involved in proteostasis, mainly represented by proteomic data. On the other hand, the concordance of RNA and protein data indicated an enrichment of metabolic pathways in RTO7, such as pentose phosphate, glycolysis, nucleotides and lipid metabolism, indicating that this PDO is characterized by a high energy demand.

All these results obtained by proteotranscriptomic integration allowed us not only to better identify those functional modules with a direct overlap between omics, but also to define new modules enriched in processes which could have a role in drug response or resistance.

Besides the limitation represented by the number of patients included, which could have underrepresented the heterogeneous mechanisms involved in drug response, we should consider the need to work on functional validation of proteo-transcriptomic data and to validate findings in additional models. In addition, the lack of TME in our models does not allow to study the interplay between cancer cells and their natural microenvironment, thus the possible effect of TME on drug response was not captured in our PDO cultures.

To our knowledge, this is among the few studies employing mass-spectrometry quantitative proteomics by SWATH to characterize PDOs and is the first to integrate proteome study with gene expression profile in the search for potential biomarkers of response/resistance to drug treatment on PDOs. Our findings are confirmed by data fusion and integrative functional network analysis, supporting the inclusion of proteomics in multi-omic characterization and its feasibility in PDOs.

Conclusions

In conclusion, we have demonstrated the feasibility of building a PDO library from advanced and pretreated CRC patients. We provided a proof-of-concept validation that PDOs, combined with a genomic and proteotranscriptomic approach, can be used to better screen and identify relevant mechanisms involved in sensitivity to anticancer agents in a functional precision medicine context.

Availability of data and material

We are in process of submitting all genomics, transcriptomics, and proteomics raw data in public repositories.

Abbreviations

AF:

Allelic fraction

ARSs:

Aminoacyl-tRNA synthetase

AUC:

Area under the curve

CRC:

Colorectal cancer

CS:

Citrate synthase

ddPCR:

Droplet-digital PCR

EMT:

Epithelial-mesenchymal transition

FDR:

False discovery rates

FFPE:

Formalin-fixed paraffin embedded

GSEA:

Gene set enrichment analysis

H&E:

Hematoxylin and eosin

IHC:

Immunohistochemistry

IFNA:

Integrative functional network analysis

LOH:

Loss of heterozygosity

MSC:

Multi-tRNA synthetase complex

MSI:

Microsatellite instability

NGS:

Next-generation sequencing

NPM1:

Nucleophosmin

OXPHOS:

Oxidative phosphorylation

PCA:

Principal component analysis

PDOs:

Patient-derived organoids

SNVs:

Single nucleotide variants

SWATH-MS:

Sequential window acquisition of all theoretical mass spectrometry

TCA:

Tricarboxylic acid cycle

TME:

Tumor microenvironment

TRiC:

T-complex protein Ring Complex

VAFs:

Variant allele frequencies

References

  1. Sung H, Ferlay J, Siegel RL, Laversanne M, Soerjomataram I, Jemal A, et al. Global Cancer Statistics 2020: GLOBOCAN Estimates of Incidence and Mortality Worldwide for 36 Cancers in 185 Countries. CA Cancer J Clin. 2021;71(3):209–49. https://doi.org/10.3322/caac.21660.

    Article  Google Scholar 

  2. Argilés G, Tabernero J, Labianca R, Hochhauser D, Salazar R, Iveson T, et al. Localised colon cancer: ESMO Clinical Practice Guidelines for diagnosis, treatment and follow-up. Ann Oncol Off J Eur Soc Med Oncol. 2020;31(10):1291–305. Available from: http://www.ncbi.nlm.nih.gov/pubmed/32702383.

    Article  Google Scholar 

  3. Van Cutsem E, Cervantes A, Adam R, Sobrero A, Van Krieken JH, Aderka D, et al. ESMO consensus guidelines for the management of patients with metastatic colorectal cancer. Ann Oncol. 2016;27(8):1386–422. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0923753419347544.

    Article  Google Scholar 

  4. Marquart J, Chen EY, Prasad V. Estimation of the Percentage of US Patients With Cancer Who Benefit From Genome-Driven Oncology. JAMA Oncol. 2018;4(8):1093. https://doi.org/10.1001/jamaoncol.2018.1660.

    Article  Google Scholar 

  5. Gambardella V, Tarazona N, Cejalvo JM, Lombardi P, Huerta M, Roselló S, et al. Personalized medicine: Recent progress in cancer therapy. Vol. 12, Cancers. MDPI AG; 2020.

  6. Vlachogiannis G, Hedayat S, Vatsiou A, Jamin Y, Fernández-Mateos J, Khan K, et al. Patient-derived organoids model treatment response of metastatic gastrointestinal cancers. Science. 2018;359(6378):920–6. https://doi.org/10.1126/science.aao2774.

    Article  Google Scholar 

  7. Driehuis E, van Hoeck A, Moore K, Kolders S, Francies HE, Gulersonmez MC, et al. Pancreatic cancer organoids recapitulate disease and allow personalized drug screening. Proc Natl Acad Sci. 2019;116(52):26580–90. https://doi.org/10.1073/pnas.1911273116.

    Article  Google Scholar 

  8. Gao D, Vela I, Sboner A, Iaquinta PJ, Karthaus WR, Gopalan A, et al. Organoid cultures derived from patients with advanced prostate cancer. Cell. 2014;159(1):176–87. https://doi.org/10.1016/j.cell.2014.08.016.

    Article  Google Scholar 

  9. Yan HHN, Siu HC, Law S, Ho SL, Yue SSK, Tsui WY, et al. A Comprehensive Human Gastric Cancer Organoid Biobank Captures Tumor Subtype Heterogeneity and Enables Therapeutic Screening. Cell Stem Cell. 2018;23(6):882-897.e11. Available from: https://linkinghub.elsevier.com/retrieve/pii/S1934590918304806.

    Article  Google Scholar 

  10. Sachs N, de Ligt J, Kopper O, Gogola E, Bounova G, Weeber F, et al. A Living Biobank of Breast Cancer Organoids Captures Disease Heterogeneity. Cell. 2018;172(1–2):373-386.e10. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0092867417313193.

    Article  Google Scholar 

  11. Van De Wetering M, Francies HE, Francis JM, Bounova G, Iorio F, Pronk A, et al. Prospective derivation of a living organoid biobank of colorectal cancer patients. Cell. 2015;161(4):933–45. https://doi.org/10.1016/j.cell.2015.03.053.

    Article  Google Scholar 

  12. Ooft SN, Weeber F, Dijkstra KK, McLean CM, Kaing S, van Werkhoven E, et al. Patient-derived organoids can predict response to chemotherapy in metastatic colorectal cancer patients. Sci Transl Med. 2019;11(513):eaay2574. https://doi.org/10.1126/scitranslmed.aay2574.

    Article  Google Scholar 

  13. Bruun J, Kryeziu K, Eide PW, Moosavi SH, Eilertsen IA, Langerud J, et al. Patient-Derived Organoids from Multiple Colorectal Cancer Liver Metastases Reveal Moderate Intra-patient Pharmacotranscriptomic Heterogeneity. Clin Cancer Res. 2020;26(15):4107–19. https://doi.org/10.1158/1078-0432.CCR-19-3637.

    Article  Google Scholar 

  14. Ooft SN, Weeber F, Schipper L, Dijkstra KK, McLean CM, Kaing S, et al. Prospective experimental treatment of colorectal cancer patients based on organoid drug responses. ESMO Open. 2021;6(3):100103. Available from: https://linkinghub.elsevier.com/retrieve/pii/S2059702921000600

    Article  Google Scholar 

  15. Tiriac H, Belleau P, Engle DD, Plenker D, Deschênes A, Somerville TDD, et al. Organoid profiling identifies common responders to chemotherapy in pancreatic cancer. Cancer Discov. 2018;8(9):1112–29.

    Article  Google Scholar 

  16. Cristobal A, van den Toorn HWP, van de Wetering M, Clevers H, Heck AJR, Mohammed S. Personalized Proteome Profiles of Healthy and Tumor Human Colon Organoids Reveal Both Individual Diversity and Basic Features of Colorectal Cancer. Cell Rep. 2017;18(1):263–74. https://doi.org/10.1016/j.celrep.2016.12.016.

    Article  Google Scholar 

  17. Boj SF, Hwang C-I, Baker LA, Chio IIC, Engle DD, Corbo V, et al. Organoid Models of Human and Mouse Ductal Pancreatic Cancer. Cell. 2015;160(1–2):324–38. Available from: https://linkinghub.elsevier.com/retrieve/pii/S009286741401592X.

    Article  Google Scholar 

  18. Fabbri L, Chakraborty A, Robert C, Vagner S. The plasticity of mRNA translation during cancer progression and therapy resistance. Nat Rev Cancer. 2021;21(9):558–77. Available from: https://www.nature.com/articles/s41568-021-00380-y.

    Article  Google Scholar 

  19. Gambardella V, Lombardi P, Carbonell-Asins JA, Tarazona N, Cejalvo JM, González-Barrallo I, et al. Molecular profiling of advanced solid tumours. The impact of experimental molecular-matched therapies on cancer patient outcomes in early-phase trials: the MAST study. Br J Cancer [Internet]. 2021 Oct 26;125(9):1261–9. Available from: https://www.nature.com/articles/s41416-021-01502-x

  20. Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884-90. Available from: http://www.ncbi.nlm.nih.gov/pubmed/30423086.

    Article  Google Scholar 

  21. Wilm A, Aw PPK, Bertrand D, Yeo GHT, Ong SH, Wong CH, et al. LoFreq: a sequence-quality aware, ultra-sensitive variant caller for uncovering cell-population heterogeneity from high-throughput sequencing datasets. Nucleic Acids Res. 2012;40(22):11189–201. Available from: http://www.ncbi.nlm.nih.gov/pubmed/23066108.

    Article  Google Scholar 

  22. McKenna A, Hanna M, Banks E, Sivachenko A, Cibulskis K, Kernytsky A, et al. The Genome Analysis Toolkit: A MapReduce framework for analyzing next-generation DNA sequencing data. Genome Res. 2010;20(9):1297–303. https://doi.org/10.1101/gr.107524.110.

    Article  Google Scholar 

  23. Danecek P, Auton A, Abecasis G, Albers CA, Banks E, DePristo MA, et al. The variant call format and VCFtools. Bioinformatics. 2011;27(15):2156–8. Available from: http://www.ncbi.nlm.nih.gov/pubmed/21653522.

    Article  Google Scholar 

  24. Li H. A statistical framework for SNP calling, mutation discovery, association mapping and population genetical parameter estimation from sequencing data. Bioinformatics. 2011;27(21):2987–93. Available from: http://www.ncbi.nlm.nih.gov/pubmed/21903627.

    Article  Google Scholar 

  25. McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GRS, Thormann A, et al. The Ensembl Variant Effect Predictor. Genome Biol. 2016;17(1):122. Available from: http://www.ncbi.nlm.nih.gov/pubmed/27268795.

    Article  Google Scholar 

  26. Chang MT, Bhattarai TS, Schram AM, Bielski CM, Donoghue MTA, Jonsson P, et al. Accelerating Discovery of Functional Mutant Alleles in Cancer. Cancer Discov. 2018;8(2):174–83. Available from: http://www.ncbi.nlm.nih.gov/pubmed/29247016.

    Article  Google Scholar 

  27. Gu Z, Eils R, Schlesner M. Complex heatmaps reveal patterns and correlations in multidimensional genomic data. Bioinformatics. 2016;32(18):2847–9. Available from: http://www.ncbi.nlm.nih.gov/pubmed/27207943.

    Article  Google Scholar 

  28. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29(1):15–21. Available from: http://www.ncbi.nlm.nih.gov/pubmed/23104886.

    Article  Google Scholar 

  29. Li B, Dewey CN. RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics. 2011;12:323. Available from: http://www.ncbi.nlm.nih.gov/pubmed/21816040.

    Article  Google Scholar 

  30. Soneson C, Love MI, Robinson MD. Differential analyses for RNA-seq: transcript-level estimates improve gene-level inferences. F1000Research [Internet]. 2015;4:1521. Available from: http://www.ncbi.nlm.nih.gov/pubmed/26925227

  31. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15(12):550. Available from: http://www.ncbi.nlm.nih.gov/pubmed/25516281.

    Article  Google Scholar 

  32. Isella C, Terrasi A, Bellomo SE, Petti C, Galatola G, Muratore A, et al. Stromal contribution to the colorectal cancer transcriptome. Nat Genet. 2015;47(4):312–9. Available from: http://www.nature.com/articles/ng.3224.

    Article  Google Scholar 

  33. Pinheiro T, Lip KYF, García-Ríos E, Querol A, Teixeira J, van Gulik W, et al. Differential proteomic analysis by SWATH-MS unravels the most dominant mechanisms underlying yeast adaptation to non-optimal temperatures under anaerobic conditions. Sci Rep. 2020;10(1):22329. Available from: http://www.nature.com/articles/s41598-020-77846-w.

    Article  Google Scholar 

  34. Szklarczyk D, Gable AL, Lyon D, Junge A, Wyder S, Huerta-Cepas J, et al. STRING v11: protein-protein association networks with increased coverage, supporting functional discovery in genome-wide experimental datasets. Nucleic Acids Res. 2019;47(D1):D607-13. Available from: http://www.ncbi.nlm.nih.gov/pubmed/30476243.

    Article  Google Scholar 

  35. Doncheva NT, Morris JH, Gorodkin J, Jensen LJ. Cytoscape StringApp: Network Analysis and Visualization of Proteomics Data. J Proteome Res. 2019;18(2):623–32. Available from: http://www.ncbi.nlm.nih.gov/pubmed/30450911.

    Article  Google Scholar 

  36. Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. Available from: http://www.ncbi.nlm.nih.gov/pubmed/14597658.

    Article  Google Scholar 

  37. Gambardella V, Gimeno-Valiente F, Tarazona N, Ciarpaglini CM, Roda D, Fleitas T, et al. NRF2 through RPS6 Activation Is Related to Anti-HER2 Drug Resistance in HER2 -Amplified Gastric Cancer. Clin Cancer Res. 2019;25(5):1639–49. Available from: https://aacrjournals.org/clincancerres/article/25/5/1639/82123/.

    Article  Google Scholar 

  38. Shi R, Radulovich N, Ng C, Liu N, Notsuda H, Cabanero M, et al. Organoid Cultures as Preclinical Models of Non-Small Cell Lung Cancer. Clin Cancer Res. 2020;26(5):1162–74. https://doi.org/10.1158/1078-0432.CCR-19-1376.

    Article  Google Scholar 

  39. Lugli A, Tzankov A, Zlobec I, Terracciano LM. Differential diagnostic and functional role of the multi-marker phenotype CDX2/CK20/CK7 in colorectal cancer stratified by mismatch repair status. Mod Pathol. 2008;21(11):1403–12. Available from: http://www.nature.com/articles/modpathol2008117.

    Article  Google Scholar 

  40. Pandey K, An H, Kim SK, Lee SA, Kim S, Lim SM, et al. Molecular mechanisms of resistance to CDK4/6 inhibitors in breast cancer: A review. Int J Cancer. 2019;145(5):1179–88. https://doi.org/10.1002/ijc.32020.

    Article  Google Scholar 

  41. Yi J, Liu C, Tao Z, Wang M, Jia Y, Sang X, et al. MYC status as a determinant of synergistic response to Olaparib and Palbociclib in ovarian cancer. EBioMedicine. 2019;43:225–37. Available from: https://linkinghub.elsevier.com/retrieve/pii/S2352396419301677.

    Article  Google Scholar 

  42. Ji W, Zhang W, Wang X, Shi Y, Yang F, Xie H, et al. c-myc regulates the sensitivity of breast cancer cells to palbociclib via c-myc/miR-29b-3p/CDK6 axis. Cell Death Dis. 2020;11(9):760. Available from: https://www.nature.com/articles/s41419-020-02980-2.

    Article  Google Scholar 

  43. Wang T-H, Chen C-C, Leu Y-L, Lee Y-S, Lian J-H, Hsieh H-L, et al. Palbociclib induces DNA damage and inhibits DNA repair to induce cellular senescence and apoptosis in oral squamous cell carcinoma. J Formos Med Assoc. 2021;120(9):1695–705. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0929664620306094.

    Article  Google Scholar 

  44. Kuang C, Chen Y. Tumor-derived C-terminal mutations of Smad4 with decreased DNA binding activity and enhanced intramolecular interaction. Oncogene. 2004;23(5):1021–9. Available from: http://www.ncbi.nlm.nih.gov/pubmed/14647410.

    Article  Google Scholar 

  45. Shi C, Yang EJ, Liu Y, Mou PK, Ren G, Shim JS. Bromodomain and extra-terminal motif (BET) inhibition is synthetic lethal with loss of SMAD4 in colorectal cancer cells via restoring the loss of MYC repression. Oncogene. 2021;40(5):937–50. Available from: http://www.nature.com/articles/s41388-020-01580-w.

    Article  Google Scholar 

  46. Lu B, Chen X, Hong Y, Zhu H, He Q, Yang B, et al. Identification of PRDX6 as a regulator of ferroptosis. Acta Pharmacol Sin. 2019;40(10):1334–42. Available from: http://www.nature.com/articles/s41401-019-0233-9.

    Article  Google Scholar 

  47. Yang C, Zhang Y, Lin S, Liu Y, Li W. Suppressing the KIF20A/NUAK1/Nrf2/GPX4 signaling pathway induces ferroptosis and enhances the sensitivity of colorectal cancer to oxaliplatin. Aging (Albany NY). 2021;13(10):13515–34. https://doi.org/10.18632/aging.202774.

    Article  Google Scholar 

  48. Durinikova E, Kozovska Z, Poturnajova M, Plava J, Cierna Z, Babelova A, et al. ALDH1A3 upregulation and spontaneous metastasis formation is associated with acquired chemoresistance in colorectal cancer cells. BMC Cancer. 2018;18(1):848. https://doi.org/10.1186/s12885-018-4758-y.

    Article  Google Scholar 

  49. Moudry P, Lukas C, Macurek L, Hanzlikova H, Hodny Z, Lukas J, et al. Ubiquitin-activating enzyme UBA1 is required for cellular response to DNA damage. Cell Cycle. 2012;11(8):1573–82. https://doi.org/10.4161/cc.19978.

    Article  Google Scholar 

  50. Tan IB, Ivanova T, Lim KH, Ong CW, Deng N, Lee J, et al. Intrinsic Subtypes of Gastric Cancer, Based on Gene Expression Pattern, Predict Survival and Respond Differently to Chemotherapy. Gastroenterology. 2011;141(2):476-485.e11. Available from: https://linkinghub.elsevier.com/retrieve/pii/S001650851100597X.

    Article  Google Scholar 

  51. Hyeon DY, Kim JH, Ahn TJ, Cho Y, Hwang D, Kim S. Evolution of the multi-tRNA synthetase complex and its role in cancer. J Biol Chem. 2019;294(14):5340–51. Available from: https://linkinghub.elsevier.com/retrieve/pii/S002192582035496X.

    Article  Google Scholar 

  52. Zhou Z, Sun B, Huang S, Yu D, Zhang X. Roles of aminoacyl-tRNA synthetase-interacting multi-functional proteins in physiology and cancer. Cell Death Dis. 2020;11(7):579. Available from: http://www.nature.com/articles/s41419-020-02794-2.

    Article  Google Scholar 

  53. Huang D, Chowdhury S, Wang H, Savage SR, Ivey RG, Kennedy JJ, et al. Multiomic analysis identifies CPT1A as a potential therapeutic target in platinum-refractory, high-grade serous ovarian cancer. Cell Reports Med. 2021;2(12):100471. Available from: https://linkinghub.elsevier.com/retrieve/pii/S2666379121003438

    Article  Google Scholar 

  54. Mao L, Li Y, Zhao J, Li Q, Yang B, Wang Y, et al. Transforming growth factor-β1 contributes to oxaliplatin resistance in colorectal cancer via epithelial to mesenchymal transition. Oncol Lett. 2017;14(1):647–54. https://doi.org/10.3892/ol.2017.6209.

    Article  Google Scholar 

  55. Liao Q, Ren Y, Yang Y, Zhu X, Zhi Y, Zhang Y, et al. CCT8 recovers WTp53-suppressed cell cycle evolution and EMT to promote colorectal cancer progression. Oncogenesis. 2021;10(12):84. Available from: https://www.nature.com/articles/s41389-021-00374-3.

    Article  Google Scholar 

  56. Goetzman ES, Prochownik E V. The Role for Myc in Coordinating Glycolysis, Oxidative Phosphorylation, Glutaminolysis, and Fatty Acid Metabolism in Normal and Neoplastic Tissues. Front Endocrinol (Lausanne) [Internet]. 2018 Apr 12;9. Available from: http://journal.frontiersin.org/article/https://doi.org/10.3389/fendo.2018.00129/full

  57. Li Z, Boone D, Hann SR. Nucleophosmin interacts directly with c-Myc and controls c-Myc-induced hyperproliferation and transformation. Proc Natl Acad Sci. 2008;105(48):18794–9. https://doi.org/10.1073/pnas.0806879105.

    Article  Google Scholar 

  58. Haider S, Pal R. Integrated Analysis of Transcriptomic and Proteomic Data. Curr Genomics. 2013;14(2):91–110. Available from: http://www.eurekaselect.com/openurl/content.php?genre=article&issn=1389-2029&volume=14&issue=2&spage=91.

    Article  Google Scholar 

  59. Bruchmann A, Roller C, Walther TV, Schäfer G, Lehmusvaara S, Visakorpi T, et al. Bcl-2 associated athanogene 5 (Bag5) is overexpressed in prostate cancer and inhibits ER-stress induced apoptosis. BMC Cancer. 2013;13(1):96. https://doi.org/10.1186/1471-2407-13-96.

    Article  Google Scholar 

  60. Gupta MK, Randhawa PK, Masternak MM. Role of BAG5 in Protein Quality Control: Double-Edged Sword? Front Aging [Internet]. 2022 Mar 3;3. Available from: https://www.frontiersin.org/articles/https://doi.org/10.3389/fragi.2022.844168/full

  61. Li X, Francies HE, Secrier M, Perner J, Miremadi A, Galeano-Dalmau N, et al. Organoid cultures recapitulate esophageal adenocarcinoma heterogeneity providing a model for clonality studies and precision therapeutics. Nat Commun. 2018;9(1):2983. Available from: http://www.nature.com/articles/s41467-018-05190-9.

    Article  Google Scholar 

  62. Huang D, Savage SR, Calinawan AP, Lin C, Zhang B, Wang P, et al. A highly annotated database of genes associated with platinum resistance in cancer. Oncogene. 2021;40(46):6395–405. Available from: https://www.nature.com/articles/s41388-021-02055-2.

    Article  Google Scholar 

  63. Narasimhan M, Kannan S, Chawade A, Bhattacharjee A, Govekar R. Clinical biomarker discovery by SWATH-MS based label-free quantitative proteomics: impact of criteria for identification of differentiators and data normalization method. J Transl Med. 2019;17(1):184. https://doi.org/10.1186/s12967-019-1937-9.

    Article  Google Scholar 

  64. Zhou Z, Sun B, Nie A, Yu D, Bian M. Roles of Aminoacyl-tRNA Synthetases in Cancer. Front Cell Dev Biol [Internet]. 2020 Nov 27;8. Available from: https://www.frontiersin.org/articles/https://doi.org/10.3389/fcell.2020.599765/full

  65. Guillon J, Coquelet H, Leman G, Toutain B, Petit C, Henry C, et al. tRNA biogenesis and specific aminoacyl-tRNA synthetases regulate senescence stability under the control of mTOR. Holz M, editor. PLOS Genet. 2021;17(12):e1009953. https://doi.org/10.1371/journal.pgen.1009953.

    Article  Google Scholar 

  66. Reisman BJ, Guo H, Ramsey HE, Wright MT, Reinfeld BI, Ferrell PB, et al. Apoptolidin family glycomacrolides target leukemia through inhibition of ATP synthase. Nat Chem Biol. 2022;18(4):360–7. Available from: https://www.nature.com/articles/s41589-021-00900-9.

    Article  Google Scholar 

  67. Lin C-Y, Lee L-Y, Wang T-H, Hsu C-L, Tsai C-L, Chao A, et al. Palbociclib Promotes Dephosphorylation of NPM/B23 at Threonine 199 and Inhibits Endometrial Cancer Cell Growth. Cancers (Basel). 2019;11(7):1025. Available from: https://www.mdpi.com/2072-6694/11/7/1025.

    Article  Google Scholar 

  68. Strippoli A, Cocomazzi A, Basso M, Cenci T, Ricci R, Pierconti F, et al. c-MYC Expression Is a Possible Keystone in the Colorectal Cancer Resistance to EGFR Inhibitors. Cancers (Basel) [Internet]. 2020 Mar 10;12(3). Available from: http://www.ncbi.nlm.nih.gov/pubmed/32164324

  69. Spiniello M, Steinbrink MI, Cesnik AJ, Miller RM, Scalf M, Shortreed MR, et al. Comprehensive in vivo identification of the c-Myc mRNA protein interactome using HyPR-MS. RNA. 2019;25(10):1337–52. https://doi.org/10.1261/rna.072157.119.

    Article  Google Scholar 

  70. Sun H, Wang Y, Jing H-Y, Yang X-Y, Shi X-X, Zhang J-H, et al. Chaperonin-Containing TCP1 Subunit 6A Is a Prognostic Potential Biomarker That Correlates With the Presence of Immune Infiltrates in Colorectal Cancer. Front Genet [Internet]. 2021 May 4;12. Available from: https://www.frontiersin.org/articles/https://doi.org/10.3389/fgene.2021.629856/full

  71. Solimini NL, Luo J, Elledge SJ. Non-Oncogene Addiction and the Stress Phenotype of Cancer Cells. Cell. 2007;130(6):986–8. Available from: https://linkinghub.elsevier.com/retrieve/pii/S0092867407011518.

    Article  Google Scholar 

  72. Stricher F, Macri C, Ruff M, Muller S. HSPA8/HSC70 chaperone protein. Autophagy. 2013;9(12):1937–54. https://doi.org/10.4161/auto.26448.

    Article  Google Scholar 

  73. Yao L, Zou X, Liu L. The TCP1 ring complex is associated with malignancy and poor prognosis in hepatocellular carcinoma. Int J Clin Exp Pathol. 2019;12(9):3329–43. Available from: http://www.ncbi.nlm.nih.gov/pubmed/31934176.

    Google Scholar 

Download references

Acknowledgements

We thank all patients participating in this study and their families. We wish to particularly acknowledge the collaboration of the INCLIVA Biobank (PT20/00029; B.000768 ISCIII) integrated in the Biobanks and Biomodels ISCIII Platform.

We thank the Proteomic Unit of the Central support service for experimental research (SCSCIE) of the University of Valencia for assistance with proteomic analysis, Dr. J. Portero-Trigo from the Central Unit of Research in Medicine, University of Valencia for assistance with the CytoScan HD analysis and the Pathology Department of HCUV for technical assistance with IHC analysis.

We would also like to thank Zahara Garzón-Lloria and Enrique Seda-Garcia for technical assistance, and the Precision Medicine Unit and the Bioinformatic and Biostatistics Units from INCLIVA for their support of this project.

Funding

This work was supported by grants from the Carlos III Health Institute [grant number PI18/01508 and PI21/0693] to TF and [grant number PI18/01909 and PI21/00689] to AC and NT. MCS was supported with a predoctoral grant by the Spanish Association against Cancer (AECC, Valencia-2019). FGV was supported by Generalitat Valenciana postdoctoral grant APOSTD/2021 (APOSTD/2021/168). TF was supported by Joan Rodés contract JR17/00026 from the Carlos III Health Institute. VG was funded by a Joan Rodés contract JR21/00042 from the Carlos III Health Institute; NT was funded by Joan Rodés contract JR20/00005 from the Carlos III Health Institute. DR was funded by Joan Rodés contract JR16/00040 from the Carlos III Health Institute. JC and MMSP were awarded a VLC Bioclinic grant 2021/257 from the University of Valencia/INCLIVA. SZT was supported by a Carlos III grant in 2018. FP was supported by ESMO Translational Research Project with the aid of a grant from Amgen; any views, opinions, findings, conclusions, or recommendations expressed in this material are those solely of the authors and do not necessarily reflect those of ESMO or Amgen. MFGB was supported with a predoctoral grant by the Conselleria d'Educació, Investigació, Cultura i Esport (Grisolía contract GRISOLIAP/2017/161). MMSP was supported by grants PID2020-119111 GB-I00 from the Spanish Ministry of Science and Innovation and PROMETEU/2019/065 from Generalitat Valenciana. Part of the equipment used in this study has been funded by Generalitat Valenciana and co-financed with ERDF funds (OP ERDF of Comunitat Valenciana 2014–2020). INCLIVA BioBank PT17/0015/0049 B.000768 ISCIII Valencian Biobanking Network Spanish National Biobanks Network. The Precision Medicine Unit of INCLIVA was supported by a donation from Fundación FERO and the Fundación para la Promoción de Acciones Solidarias. NanoString and Illumina equipment are co-financed by the European Union through the Operational Program of the European Regional Development Fund (ERDF) of the Valencian Community 2014–2020.

Author information

Authors and Affiliations

Authors

Contributions

FP conceived and designed the study, developed the methodology, acquired the data, carried out data analysis and interpretation, wrote, revised, and edited the manuscript. BGM developed the methodology and provided technical support. FGV acquired the data, carried out data analysis and interpretation and provided technical support. MCS developed the methodology and provided technical support. VG, MFGB, CMC, CAC, LS, TF, MH, SR, DR, NT acquired the data. PRG and SZT carried out data interpretation and analysis and provided technical support. JACA carried out statistical analysis. MMSP acquired the data, carried out data analysis and interpretation and provided technical support. AC conceived, designed, and supervised the study, wrote and reviewed the manuscript. JC conceived, designed, and supervised the study, developed the methodology, carried out data analysis and interpretation, wrote and reviewed the manuscript. All authors read and approved the final manuscript.

Corresponding authors

Correspondence to Federica Papaccio, Andrés Cervantes or Josefa Castillo.

Ethics declarations

Ethics approval and consent to participate

The study was approved by the ethic committee of Hospital Clínico Universitario de Valencia (2018/063, 2021/083), and all patients signed an informed consent.

Consent for publication

Not applicable.

Competing interests

AC declares institutional research funding from Genentech, Merck Serono, BMS, MSD, Roche, Beigene, Bayer, Servier, Eli-Lilly, Novartis, Takeda, Astellas, Natera and Fibrogen and advisory board or speaker fees from Amgen, Merck Serono, Roche, Bayer, Servier and Pierre Fabre in the last five years. DR declares institutional research funding from Tesaro, AbbVie, Novartis, Merus, Roche, Relay Therapeutics, HiFiBiO inc. VG declares institutional research funding from Cellcentric, Roche, Black Diamond Therapeutics, Seagen, Kinnate Biopharma, Genentech, Bayer, Eli-Lilly, Boehringher. TF declares institutional research funding from Zymeworks, AstraZeneca, Amgen, Beigene, Adaptimmune, Daiichi Sankyo. SR declares institutional research funding from Novartis, Nouscom, Pfizer, Scandion Oncology, Ability Pharmaceuticals, Pierre Fabre, Mirati Therapeutics. FP declares institutional research funding from Merck Serono. All remaining authors have declared no conflict of interest.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1 

.

Additional file 2

.

Additional file 3: Supplementary Fig. S1

. PDOs show a heterogeneous growth behavior in culture. Maximum number of days in culture of all models. Each vertical line represents a single passage. Pink: long-term cultures (more than 3 months in culture and more than 10 passages); light pink: do not meet criteria for long-term cultures; blue: short-term cultures (1-3 months in culture and 1-9 passages); light blue: models that do not meet criteria for short-term cultures; ❄ : cryopreservation.

Additional file 4: Supplementary Fig. S2

. (A) IHC for CDX2 and CK20 antibodies in corresponding tissues. (Scale bar 50 mm, 20X). (B) Lack of expression of MLH1 and PMS2 in tissue corresponding to mCTO50B and RTO2 (Scale bar 50 mm, 20X).

Additional file 5: Supplementary Fig. S3

. (A) Pearson correlation of SNVs, insertions and deletions across all our cohort. Mutated genes have been filtered per frequency (at least 5% in PDOs). (B) forest plot representing the distribution of r Pearson correlation of mutated genes in each PDO line and corresponding tissue. p<0.0001. (C) Linear regression between pathological features and NGS concordance. Tumor cells (%) assessed as the percentage of neoplastic cells with respect to the total amount of viable cells in the tissue sample from which the culture was derived. TSR (tumor-stroma ratio) assessed as the percentage of stroma in a 10x hotspot field, when tumor cells are present in four fields (TSH hotspot) or as average TSR evaluated in up to ten 10X fields. (D) Copy number variation (CNV) profile of ddPCR assays of main driver genes across selected PDOs lines. (Red line: diploid status). CNT33 is a normal genomic DNA sample and is located at the right of the figure in green; samples are depicted in CNV ascending order. The CNV is assumed to follow a Poisson distribution and values represent the estimated number of copies with a 95% confidence interval. Copy number above two means amplification in that region and copy number below two means deletion in that region. (E) Individual plot of first two dimensions using principal componentanalysis of normalized VST showing the distribution of organoid lines and tissues.

Additional file 6:

Supplementary Fig. S4. (A) Scatterplot for technical replicates of drug screening data. Correlation between the three different technical replicates. Each data point represents the normalized value for an individual organoid line. (B) Log transformed dose-response curves for 5Fluorouracil, oxaliplatin and SN38 and combination respectively and ΔAUC calculation for each line (a negative value or positive indicates presence or absence of additive or synergistic effect, respectively). (C) Log transformed dose-response curves in selected standard drugs. (D) Log transformed dose-response curves in selected non-standard drugs.

Additional file 7: Supplementary Fig. S5

. (A) Log2 TPM+1 of CDK4 and CDK6 expression in RTO7, BRD2 and BRD4 expression in RTO2, c-MYC and SMAD4 in RTO7 compared with the other cultures. (B) PTEN IHC in MSI PDOs and RTO7 as positive control. (C) Principal component analysis of protein expression showing the distribution of organoid lines and tissues. Variance absorption from PC3 and PC4: 21.98%. (D) Log2 protein abundance (PA) of TriC complex proteins in RTO7 versus RTO2/mCTO66S3. (E) Western blot analysis of TCP1, CCT2 and CCT8 proteins. GADPH is included as control. TPM: transcripts per million; PA: protein abundance.

Additional file 8: Supplementary Fig. S6

. (A) Correlation between gene (Log2 TPM) and protein (Log2 PA) expression. Differentially expressed proteins have been matched with corresponding genes. (B) Log2 TPM+1 of TRiC complex and ARSs gene expression in palbociclib (upper panel) andoxaliplatin (lower panel) comparisons. TPM: transcripts per million; PA: protein abundance.

Additional file 9: Supplementary Fig. S7

. Differentially expressed proteins mapped in KEGG metabolic pathways (hsa01100) in the palbociclib comparison related to the metabolic enrichment group. In red those pathway reactions catalyzed by proteins identified by RNA-seq only, in blue those by proteomics only, in green those identified by both omics.

Additional file 10: Supplementary Table S1.

All sample characteristics are represented among our PDOs cohort. P value refers to two-sided Fisher’s exact test for pre-treated and RAS/RAF status variables and to Chi-Square test for PTL and site of biopsy, being Chi-square results 0,4759, 2 and 2,671, 2 respectively.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Papaccio, F., García-Mico, B., Gimeno-Valiente, F. et al. “Proteotranscriptomic analysis of advanced colorectal cancer patient derived organoids for drug sensitivity prediction”. J Exp Clin Cancer Res 42, 8 (2023). https://doi.org/10.1186/s13046-022-02591-z

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s13046-022-02591-z

Keywords

  • Colorectal cancer
  • Organoids
  • Quantitative proteomics
  • Precision medicine
  • Drug resistance
  • Transcriptomics
  • Proteotranscriptomics integrative functional network analysis