Secondary resistance to anti-EGFR therapy by transcriptional reprogramming in patient-derived colorectal cancer models
Genome Medicine volume 13, Article number: 116 (2021)
The development of secondary resistance (SR) in metastatic colorectal cancer (mCRC) treated with anti-epidermal growth factor receptor (anti-EGFR) antibodies is not fully understood at the molecular level. Here we tested in vivo selection of anti-EGFR SR tumors in CRC patient-derived xenograft (PDX) models as a strategy for a molecular dissection of SR mechanisms.
We analyzed 21 KRAS, NRAS, BRAF, and PI3K wildtype CRC patient-derived xenograft (PDX) models for their anti-EGFR sensitivity. Furthermore, 31 anti-EGFR SR tumors were generated via chronic in vivo treatment with cetuximab. A multi-omics approach was employed to address molecular primary and secondary resistance mechanisms. Gene set enrichment analyses were used to uncover SR pathways. Targeted therapy of SR PDX models was applied to validate selected SR pathways.
In vivo anti-EGFR SR could be established with high efficiency. Chronic anti-EGFR treatment of CRC PDX tumors induced parallel evolution of multiple resistant lesions with independent molecular SR mechanisms. Mutations in driver genes explained SR development in a subgroup of CRC PDX models, only. Transcriptional reprogramming inducing anti-EGFR SR was discovered as a common mechanism in CRC PDX models frequently leading to RAS signaling pathway activation. We identified cAMP and STAT3 signaling activation, as well as paracrine and autocrine signaling via growth factors as novel anti-EGFR secondary resistance mechanisms. Secondary resistant xenograft tumors could successfully be treated by addressing identified transcriptional changes by tailored targeted therapies.
Our study demonstrates that SR PDX tumors provide a unique platform to study molecular SR mechanisms and allow testing of multiple treatments for efficient targeting of SR mechanisms, not possible in the patient. Importantly, it suggests that the development of anti-EGFR tolerant cells via transcriptional reprogramming as a cause of anti-EGFR SR in CRC is likely more prevalent than previously anticipated. It emphasizes the need for analyses of SR tumor tissues at a multi-omics level for a comprehensive molecular understanding of anti-EGFR SR in CRC.
The addition of epidermal growth factor receptor (EGFR) targeting monoclonal antibodies (mABs) cetuximab and panitumumab to treatment protocols has improved survival rates in patients with metastatic colorectal cancer (CRC) as has been shown in numerous clinical trials [1, 2]. An early finding in the treatment of CRC patients with anti-EGFR targeting antibodies was that only a subgroup benefitted from treatment. Clinical trials report a median overall survival benefit of approximately 8 months for CRC without KRAS and NRAS mutations for anti-EGFR antibody-containing regimens . Despite initial clinical benefit from treatment with anti-EGFR mABs, virtually all patients acquire resistance to anti-EGFR therapy within 10 to 12 months, eventually leading to disease progression.
Currently, the acquisition of tissue and subsequent genotyping of metastatic CRC lesions with suspected acquired or secondary resistance (SR) to anti-EGFR antibodies presents a challenge even in the most advanced clinical practice and is not routinely performed in clinical workup at disease progression. Therefore, molecular mechanisms of acquired SR to anti-EGFR have been primarily addressed by cell culture models, in a limited number of tissue biopsies from cetuximab (CET)-treated patients and very few patient-derived xenograft (PDX) models [4,5,6,7,8,9,10,11,12,13,14]. By far the largest fraction of the current molecular data on SR was generated indirectly by analyzing the emergence of mutated circulating tumor DNA (ctDNA) in the blood of CRC patients under anti-EGFR treatment. This led, apart from the identification of FGFR1, ERBB2, NRG, GNAS, and MET amplifications, to the discovery of emerging KRAS, NRAS, PIK3CA, BRAF, EGFR, and IRS1 mutations in resistant CRC subpopulations mediating SR towards anti-EGFR therapy [4, 6, 7, 12, 13, 15,16,17,18]. However, several contradictory reports on the prevalence of SR driver mutations derived from primary CRC tissue appeared in recent years [4,5,6,7,8, 11,12,13,14]. Furthermore, the emergence of RAS mutations in ctDNA at the time of disease progression in patients undergoing anti-EGFR therapy was reported to be neither associated with a shorter progression-free survival nor predictive for any cytoreduction [8, 19]. This suggests that as tumors undergo clonal expansion in response to targeted therapy, not every mutation identified in ctDNA may be driving SR and that other molecular mechanisms may play a role in acquiring resistance to anti-EGFR antibodies. Furthermore, mathematical modeling of the dynamic of resistance development in PDXs under vertical blockade of the EGFR pathway by Misale et al. concluded that their data supported a model which beyond point mutations in driver genes other non-genetic or genetic mechanisms are required to achieve the high number of cetuximab-resistant cells at the start of treatment .
Here we asked the question, whether in vivo selection of anti-EGFR SR tumors from multiple CRC-based PDX models is an efficient strategy to generate broad collections of SR PDX tumors. The availability of such collections would provide an abundant source of tissues from SR PDX tumors to allow in-depth molecular examinations that are currently not feasible to perform from on-treatment patient tissue biopsies. Furthermore, comprehensive pre-testing of next-generation compounds and various targeted therapy strategies aimed at overcoming SR mechanisms would greatly improve the identification of promising therapeutic combinations to be considered for clinical evaluation.
PDX biobank and treatment cohorts
Tissue samples to establish the colon cancer PDX bank were collected from 154 patients following surgical intervention for colon cancer at the Ruhr University Comprehensive Cancer Center or the Comprehensive Cancer Center of the Ludwig-Maximilians University Munich. Each PDX model derived from an individual patient tumor received a unique identification number (ID), preceded by the letters “BoC” (“Bo” for Bochum as the site where the PDX tumor was established, and “C” for colorectal carcinoma). The ID is then followed by one of the three-letter codes indicating the treatment cohort the individual mouse was randomly assigned to (K, for untreated control; 5dC for mice treated for 5 days with cetuximab; C for mice chronically treated with cetuximab) followed by another ID number for each of the two implantation sites and tumors potentially growing per site and mouse. From the PDX bank, we randomly selected 29 PDX models representing 29 individual patient’s tumors for the treatment study reported herein. A total of 21 PDX models with wildtype sequences based on targeted NGS sequencing in the genes KRAS, NRAS, PIK3CA, and BRAF (therefore called quadruple wildtype) as well as 8 PDX models harboring an activating mutation in one of the aforementioned genes were selected. Informed and written consent was obtained from all patients. The study was approved by the Ruhr University Bochum Ethics Committee, Registry no. 3841-10 & 16-5792, and the Ludwig-Maximilians University Munich Ethics Committee, project no. 131-16. All animal experiments were approved by the local authorities (84-02.04.2012.A360 & 84-02.04.2015.A135) and performed in accordance with the guidelines for Ethical Conduct in the Care and Use of Animals.
To establish treatment cohorts, tumor pieces (1–2 mm) from early passage PDXs (≤ F4 generation) were soaked in undiluted matrigel (Becton Dickinson) for 15 to 30 min and subsequently implanted subcutaneously onto 5- to 10-week-old female mice (NMRI-Foxn1nu/Foxn1nu, Janvier, St Berthevin Cedex, France) at two sites (scapular region) using as many as 4 pieces per site. Tumors were allowed to grow to a size of approx. 100–200 mm3, at which time mice were randomized in the treatment and control groups with five to six mice in each group. Tumor volumes were estimated from 2-dimensional tumor measurements by bi-weekly caliper measurements during the first 21 days and weekly measurements thereafter using the following formula: tumor volume (mm3) = [length (mm) × width (mm)2]/2. Complete response (CR) was defined as an undetectable tumor by macroscopical inspection and partial response (PR) by at least a 30% reduction in mean tumor volume compared to the mean tumor volume at the start of treatment. Disease progression was defined as a more than 20% increase in mean tumor volume determined at least at two consecutive time points compared to the tumor volume at the beginning of treatment. All other measurements were defined as stable disease. Growth curves were established by determining mean tumor volumes at different time points relative to the mean tumor volume at treatment start. For BoC106, time to progression (20% increase in mean tumor volume) was calculated relative to the mean tumor volume at week 24 (re-initiation of CET treatment). The identity of corresponding PDX tumors was secured via STR analyses using the GenomeLab Human STR Primer Set (Beckman Coulter, Krefeld, Germany) according to the manufacturer’s protocol. STR reactions were run on a CEQ8800 sequencer (Beckman Coulter). Mice were treated with CET (Merck KGaA) by intraperitoneal (i.p.) injection twice per week and dosed at 25 mg/kg, rolipram (Hycultec) daily i.p. at 1.6 mg/kg, vorinostat (Hycultec) daily i.p. at 25 mg/kg, refametinib (Selleckchem) daily by oral gavage at 25 mg/kg, ruxolitinib (Hycultec) twice daily by oral gavage at 180 mg/kg, Ly2874455 (MedKoo Biosciences) twice daily by oral gavage at 3 mg/kg, SCH772984 (MedChemTronica ) twice per day i.p. for five times per week with two consecutive days of treatment pause at 12.5 mg/kg, and trametinib (MedKoo Biosciences) five times per week with two consecutive days of treatment pause by oral gavage dosed at 0.5 mg/kg.
BoC20 was derived via the standard outgrow method. The identity of BoC20 was confirmed by comparing STR data between the BoC20 PDX founder tumor and the derived primary cell line. DiFi cells were provided by P. Bastiaens, Systemic Cell Biology, Max Planck Institute of Molecular Physiology, Dortmund, Germany. The authenticity of DiFi cells was confirmed by STR profiling and comparison to the STR data in the Expasy database  (https://web.expasy.org/cellosaurus/CVCL_6895). STR analyses were performed as described by Dirks W.G. et al.  and analyzed on a CEQ8800 sequencer (Beckman Coulter). BoC20 cells were cultured in DMEM/DMEM-F12 medium (Thermo Fisher Scientific Inc.) supplemented with Rock-Inhibitor Y-27632 (Selleckchem) and 5% fetal calf serum (FCS) and penicillin/streptomycin (Thermo Fisher Scientific Inc.). DiFi cells were cultured in standard DMEM supplemented with penicillin/streptomycin, L-glutamine, sodium pyruvate, and 10% FCS (Thermo Fisher Scientific Inc.). Cells were maintained in a humidified incubator with 5% CO2 at 37 °C. BoC20 and DiFi cells were confirmed negative for mycoplasma by routine testing.
Two hundred fifty nanograms genomic DNA was used per sample to produce sequencing libraries with the TruSeq Amplicon Cancer Panel (Illumina) according to the manufacturer’s protocol. The chosen panel covered 48 cancer-related genes (Additional file 1: Table S1), interrogating among others the known mutations in genes involved in primary and secondary anti-EGFR resistance. The 212 amplicons were simultaneously amplified in a single tube reaction. Briefly, the regions of interest were enriched by hybridizing sequence-specific oligonucleotides to the genomic DNA followed by ligation extension of the bounded oligos. The marked regions were further amplified by PCR with primers containing index barcodes for sample multiplexing. Finally, libraries were normalized by bead normalization prior to sequencing. Pooled libraries were sequenced on a MiSeq instrument (Illumina) using 2 × 150 bp paired-end reads. For data processing, fastq files were analyzed with the NextGENe V2.3.4 (SoftGenetics) software. For variant calling, raw reads were aligned to the human hg 19 assembly and primer sequences were soft-clipped prior to variant calling. Variants with a minor AF of ≥ 5% within the coding region and a minimum coverage of 10 variant reads were considered as alteration and visually confirmed with NextGENe.
DNA quality was assessed on the Agilent TapeStation System (Agilent Technologies, Inc., Santa Clara, USA). Illumina libraries were constructed with SureSelect Human All Exon V7 (Agilent Technologies, Inc.) enrichment panel following the manufacturer’s instructions. All libraries were sequenced on the Illumina NovaSeq 6000 device in paired-end mode with 2 × 101 bp reads targeting an average coverage of > 200×. Reads were preprocessed with SeqPurge  (version 2019_05) using default parameters to remove adapter contamination and low-quality bases. For read alignment, a graft genome was constructed by concatenating the human reference genome (build GRCh37) and the mouse reference genome (build GRCm38p6). Alignment to the graft genome was performed with BWA-MEM  (version 0.7.17; https://arxiv.org/abs/1303.3997). Alignments were then reduced to human chromosomes only and with mapping quality of at least q = 30, in addition alignments filtered to allow a maximum of 3 mismatches and 1 insertion or deletion per alignment. Alignment quality was analyzed using MappingQC  (version 2019_05). FreeBayes  (https://arxiv.org/abs/1207.3907) was used to call variants on the filtered alignment data, with an AF cut-off of 5%. The remaining variant calls were annotated with Ensembl Variant Effect Predictor  to assess variant impact and to add population frequencies from gnomAD, ExAC, and 1000 Genomes databases [27,28,29] and variants with a population frequency of more than 1% in any of the databases were discarded. Mutation significance was assessed with Cancer Genome Interpreter  (https://www.cancergenomeinterpreter.org/home), further annotating alterations as known or predicted cancer drivers. Only known and predicted (tier 1) drivers are reported as likely pathogenic. Tier 2 predicted drivers are reported as of unclear significance. Exome sequencing has been deposited into the NCBI BioProject database under the BioProject ID PRJNA596887 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA596887) .
RNA quality was assessed with the Agilent 2100 Bioanalyzer RNA Nano kit and Agilent Fragment Analyzer total RNA kit (Agilent Technologies, Inc., Santa Clara, USA). Samples with RNA integrity number > 7 were selected for library construction. QuantSeq 3′ mRNA-Seq Library Prep Kit FWD UMI (Lexogen, Greenland, USA) was used with 100 ng of total RNA input to construct Illumina libraries according to the manufacturer’s manual. Libraries were sequenced on Illumina NextSeq 550 platform in single-end mode with 150 bp reads targeting 10 million reads per sample. Library preparation and sequencing procedures were performed by the same individual, and a design aimed to minimize technical batch effects was chosen. Quality of raw RNA-seq data in FASTQ files was evaluated to identify potential sequencing cycles with low average quality and base distribution bias and reads were then preprocessed using fastp  (version 0.20.0). Cleaned reads were aligned using STAR  (version 2.7.3a) allowing spliced read alignment to a graft genome containing human and mouse chromosomes (see DNA). Alignment quality was analyzed using MappingQC  (version 2019_08) and visually inspected with Broad Integrative Genome Viewer  (version 2.7.2). Alignments were collapsed based on the unique molecular identifier sequence (UMI) with UMI-tools  (version 1.0.0) to obtain one alignment per original fragment. Based on the Ensembl genome annotation (GRCh37 and GRCm38, Ensembl release 97) and using only uniquely aligned reads across the whole grafted reference, read counts for all human and mouse genes were obtained using subread  (version 2.0.0). The RNA sequencing data from our study have been deposited in the NCBI’s Gene Expression Omnibus (GEO) database (accession number GSE141861; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141861) .
Copy number variation analysis
For BoC10 C3 & BoC20 C9, copy number variation (CNV) data were extracted from an independently generated set of exome sequencing data. DNA was isolated from the SR PDX tumors BoC10 C3 and BoC20 C9, respectively, as well as from one untreated control PDX tumor each (BoC10 K14 and BoC20 K11), and from corresponding patient normal tissue of both PDX models using the DNeasyBlood & Tissue Kit or QIAamp DNA MiniKit (Qiagen), and DNA concentrations were assessed by QubitTM dsDNABR AssayKit500 (Life Technologies). Tissue authentication was performed using the Investigator ESSplex Plus Kit (Qiagen), and DNA quality was verified using an Agilent2200 TapeStation. Library preparation was performed on an Agilent NGS Workstation (version F.0) using the SureSelectXT Automation Reagent Kit (Agilent) and Human All Exon V5 + UTRs reagents (Agilent). Paired-end sequencing (2 × 100 bp) was done on a HiSeq2000 applying TruSeqPEClusterKitv3 and TruSeqSBSKitv3 (Illumina). Paired-end DNA sequencing reads were mapped to a concatenated human and murine genome assembly to reduce false positive mutations introduced by murine contamination in the xenograft samples (1000 genomes  phase 2 reference assembly (hs37d5) (ftp://ftp.1000genomes.ebi.ac.uk/vol1/ftp/technical/reference/phase2_reference_assembly_sequence/) and the UCSC mm10 genome assembly  (http://hgdownload.cse.ucsc.edu/goldenPath/mm10/bigZips/chromFa.tar.gz)). Reads were aligned using BWA  (version 0.6.2-r126-tpx) allowing for soft clipping of bases of up to Phred score Q20 with a Q20 Phred score quality threshold for read trimming down to 35 bp (-q 20). Subsequent alignments were merged using samtools  (version 0.1.19-44428 cd), and PCR duplicates were marked using Picard tools  (version 1.90). Exome sequencing has been deposited into the European Genome-phenome Archive (EGA) database under the EGA study ID EGAS00001005320 (https://ega-archive.org/studies/EGAS00001005320) . VarScan2  was used to identify somatic CNAs in the exome samples (with the preceding sample used as the reference control) using the recommended workflow. Regions were filtered for unmappable genomic stretches and merged by requiring at least 70 markers per called copy number event. We selected regions with a log ratio of tumor coverage over control coverage higher than 0.55 or lower than − 0.55 as copy number gains and losses, respectively, and annotated them with RefSeq genes using BEDTools . We searched for structural variants such as translocations that might lead to gene fusions with CREST  on the DNA level.
For BoC35 C1, BoC60 C5, BoC69 C5, and BoC106 C2, the SNPs (single-nucleotide polymorphisms) represented on the Infinium® MethylationEPIC BeadChip Array were used for CNV analyses. Un-normalized probe signal intensities were extracted from the ChAMP Bioconductor package  after default filtering and removal of potential mouse contaminated probes. Therefore, mouse tail DNA was hybridized to the EPIC array and 325,107 probes that passed a detection p value cut-off of 0.01 were excluded. Preprocessed data were imported into Partek Genomic Suite (version 7.18.0130). SR sample intensities were normalized to their corresponding 5-day CET treated controls, respectively. After log2 transformation copy numbers were adjusted and genomic copy number segmentation was performed. Different from default settings, a minimum of 50 genomic markers was used, while gender information was ignored. A copy number < 1 was termed as deleted, and an amplification had to exhibit a copy number > 3.
Gene expression analyses and data processing.
An amount of 100 ng of every total RNA sample was hybridized to Agilent whole-genome expression microarrays (Human GE 4x44K, v2 G4845A, AMADID 026652, Agilent Technologies). RNA labeling, hybridization, and washing were carried out according to the manufacturer’s instructions. Images of hybridized microarrays were acquired with a DNA microarray scanner (Agilent G2505B), and features were extracted using the Agilent Feature Extraction image analysis software (AFE) version A.10.7.3.1 with default protocols and settings. The AFE algorithm generates a single intensity measure for each feature, referred to as the total gene signal (TGS), which was used for further data analyses using the GeneSpring GX software package version 14.9.1. AFE - TGS were normalized by the quantile method. Subsequently, data were filtered on normalized expression values. The gene expression data from our study have been deposited in the NCBI’s Gene Expression Omnibus (GEO) database (accession number GSE140973; https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE140973) .
Vectors and site-directed mutagenesis
The coding sequence of PPP2R5A was inserted into the pSFFV-d4EGFP vector (gift from L. Naldini, ) by replacing the d4GFP coding sequence. Site-directed mutagenesis was performed using QuikChange II XL (Agilent) and custom-designed primers (PPP2R5A-R112L-mut-s: TGGTTGAGTATGTTTCAACTAATCTTGG-TGTAATTGTTGAATCAG and PPP2R5A-R112L-mut-as: CTGATTCAACAATTACACCAAGATTAGTTGAAACAT-ACTCAACCAG) to generate the PPP2R5A R112L mutant. The full-length sequence of each vector was assessed using Sanger sequencing to confirm the presence of the wildtype or the intended mutation and that no other mutations had been inserted.
Transfection and transduction
HEK293T cells were transfected with pSFFV-d4EGFP (control), pSFFV-PPP2R5A, or pSFFV-PPP2R5A-R112L lentiviral vectors (12 μg) in combination with packaging plasmids pCMV delta R8.2 (12 μg, gift from D. Trono, Addgene #12263) and pHIT/G (6 μg, gift from M.H. Malim, School of Immunology & Microbial Sciences, King’s College, London, Great Britain) via standard calcium phosphate transfection. Thirty-six hours post transfection, lentiviral-supernatants were collected and filtered (0.45 μm pore size) for subsequent infection of DiFi or BoC20 cells in the presence of 4 μg/ml polybrene.
Purified DNA of PDX mouse model tumors was sent to Sysmex Inostics, Hamburg, Germany. DNA content was quantified using a modified version of LINE-1 Real-time PCR as previously described [48, 49]. An amount of 330 ng DNA per sample was inputted for each tissue BEAMing analysis. BEAMing digital PCR is a technique in which individual DNA molecules are attached to magnetic beads in water-in-oil emulsions and then subjected to compartmentalized (digital) PCR amplification . The mutational status of DNA bound to beads was then determined by hybridization of fluorescently labeled allele-specific probes for mutant or wildtype KRAS G12V (g35t) sequences. Finally, the bead population was analyzed by flow cytometry to count and sort wildtype and mutant beads. The result is reported as the fractional abundance of mutant DNA alleles relative to wildtype DNA alleles per sample. To generate the ratio of mutant to wildtype DNA alleles (mutant allelic fraction, MAF), an average of 3 × 106 beads were interrogated in each BEAMing analysis (approximately 90,000 beads per mutation). The cut-off for calling a mutation in xenograft derived DNA with BEAMing was set to 0.02% MAF.
Droplet digital PCR
All ddPCR assays were carried out in duplicate according to the manufacturer’s instructions (BioRad Laboratories Inc). For each reaction, a standard volume of 20 μl was prepared using 6 μl DNA containing 50–100 ng DNA and 14 μl mastermix (containing 10 μl 2× Supermix, 1 μl mutation wildtype multiplex assay mixture, 2.75 μl DNA-free water, 0.25 μl (40 U) HaeI restriction enzyme). KR2-12/13 and KR3-61 multiplex screening assays were proved by titration of standard samples purchased from Horizon discovery (H701, Wien; data not shown). The ddPCR amplifications were performed with the following protocol in the BioRad cycler: initial denaturation 95 °C 10 min, 45 cycles with 94 °C 30 s, 55 °C for 1 min, with a ramp rate of 2 °C per sec, and 10 min 98 °C final step. The fluorescence of each droplet was analyzed by QX200 droplet reader, and the data were calculated by QuantaSoft software (BioRad) with Poisson’s algorithm. The assay information of primer sequence context is given at the homepage of BioRad systems.
Real-time cell analyses (RTCA)
RTCA was performed using 16-well E-plates on the Dual Plate xCELLigence instrument (ACEA Biosciences Inc. San Diego, CA, USA, RTCA software v2.0). This system measures a dimensionless parameter called cell index (CI), which evaluates the ionic environment at an electrode/solution interface and integrates information on cell number. The assays were performed according to the manufacturer’s instructions. Briefly, 100 μl of culture media was added to each well, incubated at room temperature for 30 min and the background impedance was measured. Prior to the RTCA experiments, cells were cultivated in low-serum-medium containing 0.5% FCS for 24 h. For RTCA analyses, cells were seeded at a density of between 5000 and 30,000 cells per well in a final volume of 200 μL into E-plates (16 wells, ACEA Biosciences) and allowed to settle for 30 min at RT before the plates were inserted into the xCELLigence instrument. RTCA assays were performed in a humidified incubator at 37 °C with 5% CO2. The growth was monitored at 15-min time intervals for up to 165 h. All measurements were performed in duplicates. After 24 h, selected murine growth factors and/or CET were added. Substances used were cetuximab (Merck KGaA), Fgf9 & 10, Igf2, Pdgfc (RND Systems), Fgf18 (NOVUS Biologicals), Pdgfb, and Hbegf (Sigma-Aldrich).
Quantitative real-time PCR (qRT-PCR)
cDNA was synthesized using 2 μg of total RNA, oligo(dT)18 primers, and M-MLV reverse transcriptase (Promega) following the manufacturer’s protocol and diluted to a final volume of 50 μl with 1× first strand buffer. Intron spanning primer sets for qRT-PCR were designed using Primer Express 2.0 software (Applied Biosystems) (PPP2R5A-1141-s: AGGATGAACCCACGCTTGAG, PPP2R5A-1233-as: TAGGCTGG AAATCAGGGCTCT, FGF3-693-s: CTGGAGAACAGCGCCTACAGT, FGF3-815-as: CGAAGC ATAGAGTCGTCCCCT). qRT-PCR was performed using a SYBR Green I reaction mixture containing 75 mM Tris-HCl (pH 8.8), 20 mM ammonium sulfate, 0.01 % (v/v) Tween 20, 2 mM magnesium chloride (all Sigma-Aldrich), 1 μl of a 600-fold dilution of SYBR Green I (BioWhittaker), 2.5 U Taq polymerase (NEB), 0.2 mM dNTP (Promega), and 0.2 μM of forward and reverse primer in a final reaction volume of 20 μl. Reactions were run on a CFX Connect Real Time System (BioRad). The cycling conditions consisted of 3 min initial denaturation at 94 °C and 40 cycles of 94 °C for 30 s, 60 °C for 30 s, 72 °C for 30 s, and 80 °C for 5 s. Fluorescence was measured at the last step of each cycle. Melting curves were obtained after each PCR run and showed single PCR products. cDNAs were run in triplicate, non-RT (without reverse transcriptase), and no-template controls were run in duplicates. Expression ratios were calculated as described by M. Pfaffl  using the geometric mean expression of the housekeeping gene COX6C (COX6C-108-s: CAGGCGTCTGCGAAATCATA; COX6C-265-as: CAGCCTTCCTCATCTCCTCA GAPDH-s: TGCACCACCAACTGCTTAGC, GAPDH-as: GGCATGGACTGTGGTCATGAGA) to normalize the expression data for the gene of interest.
Protein extraction and immunoblotting
Cryoconserved xenograft tumors were sectioned at 30 μm using a Leica CM1900 cryostat. Between 20 and 30 sections or cell pellets from 1 × 107 cells were placed in 1.5-ml tubes and lysed in 0.5 ml cold RIPA lysis buffer, supplemented with protease and phosphatase inhibitors (Sigma-Aldrich), the lysates were subsequently sonicated and cellular debris was removed by centrifugation. Protein concentrations of the lysates were determined by the Bradford protein assay system (BioRad). Equal amounts of protein (10–50 μg each lane) were separated by SDS-PAGE and transferred to PVDF membranes (Carl Roth). Immunoblots were blocked with 5% BSA in tris-buffered saline and tween-20 (0.1%, v/v) for 1 h at room temperature. The membrane was incubated overnight at 4 °C with primary antibodies. Antibodies were purchased from cell signaling technology: Total-STAT3 (#4904), P-Y705-STAT3 (#4113), P-T202/Y204-ERK1/2 (#9101), P-S9-GSK-3b (#9323), P-S473-AKT (#4060), P-S62-c-Myc (#13748), P-Y1045-EGFR (#2237), P-Y1068-EGFR (#3777); Abcam: anti-PPP2R5A (ab89621). Anti GAPDH (#2118) and Anti-β-actin (#4967S; Cell Signaling Technology) were used as protein loading controls. The membrane was incubated with the corresponding secondary antibody conjugated with horseradish peroxidase (Dianova). Bands were visualized with enhanced chemiluminescence western blot detection system (Thermo Scientific) together with the ChemiDoc imaging system (BioRad). Signal intensities were quantified using Image Lab software (BioRad). Band intensity was normalized to that of GAPDH or β-actin.
Gene set enrichment analyses (GSEA)
GSEA software [52, 53] was provided by the Broad Institute of the Massachusetts Institute of Technology and Harvard University (http://www.broad.mit.edu/gsea/). The hallmark gene sets (V6.1) were used with default parameters of the GSEA software package; gene set permutation was used. Gene sets with FDR q-val ≤ 0.05 were considered appropriate.
Gene expression data were loaded as data matrix in a linear fashion into the CRIS NTP classifier  and run under standard parameters using RStudio (version 1.2.1335) and the CRIS NTP classifier program which is part of the CRISclassifier R package (available via https://www.nature.com/articles/ncomms15107#Sec29).
Unless stated otherwise, results are expressed as the means ± SEM. For the identification of differentially expressed genes, only entities where at least 3 out of the total no. of samples had values within the selected cut-off (50th–100th percentile) were further included in the data analysis process. Using the GeneSpring GX software package version 14.9.1, differentially expressed genes were identified via moderated t-test. The p values were adjusted for multiple testing according to Benjamini and Hochberg [false detection rate—FDR] and results were considered statistically significant at adjusted p values below 0.05. Lastly, only mRNAs with a fold change ≥ 1.7 in the microarray analyses were further considered. All other differences were evaluated by the nonparametric Mann-Whitney U test. Statistical computations were performed using Prism software (Graph Pad, La Jolla, CA, USA).
Characterization of the in vivo anti-EGFR secondary resistance PDX platform
Since May 2012, 158 patients with CRC undergoing surgery at our centers were prospectively included in our study to establish a representative colorectal cancer PDX cohort (Fig. 1a). We were able to reach an overall engraftment rate of 80%. Passaging of tumors from mice to mice to establish the PDX lines was highly successful (98%) excluding a strong selection bias in our xenograft cohort. Xenografted tumors retained the histopathologic characteristics of original samples (Additional file 2: Fig. S1). We randomly selected 21 PDX models with wildtype sequences based on targeted NGS sequencing in the genes KRAS, NRAS, PIK3CA, and BRAF (therefore called quadruple wildtype) as well as 8 PDX models harboring an activating mutation in one of the aforementioned genes to establish their anti-EGFR treatment response characteristics (Fig. 1b, c). In line with previous reports on genetic primary resistance mediators, PDX models with activating mutations in KRAS, NRAS, PIK3CA, or BRAF did not show any disease control with CET treatment during the primary observational period of 21 days (Fig. 2a, Additional file 2: Fig. S2). Within the quadruple wildtype group, three growth patterns could be observed: progressive disease for 6 PDX models (29%), objective response for 7 PDX models (33%), and stable disease for 8 PDX models (38%) (Figs. 1c and 2a; Additional file 2: Fig. S3-5). The observed primary response pattern frequencies are in line with previous PDX model results and clinical response data . Next, we asked if it is possible to generate in vivo selected secondary resistant tumors in the PDX model via chronic CET treatment for the 15 models (Fig. 1b) showing initial disease control (stable disease and objective response). Indeed, 12 of the 15 models developed secondary resistant tumors defined by disease progression following an initial phase of disease control of at least 3 weeks. The efficiency with which SR PDX tumors evolved in the 12 models varied (Additional file 3: Table S2). In 10 models, we were able to generate two to five independently arising SR, in two models only one SR tumor each, could be established. The main reason for the latter was either health issues of the animals forcing to terminate the experiment (BoC35) or a low number of initial tumors growing in the cohort (BoC60). Similarly, the time to progression in the individual parallel PDX tumors within each model varied by several weeks (Fig. 2b, c; Additional file 2: Fig. S6).
Genetic drivers of primary and secondary resistance
First, we re-sequenced untreated control tumors from the first five consecutively arising models with primary resistance via whole-exome sequencing (WES) (Additional file 4: Table S3). The mutational load was similar in all models, except three models with genetic instability due to inactivating mutations in mismatch repair (MMR) genes (Fig. 3a, Additional file 4: Table S3). These MMR-deficient (dMMR) tumors were from the group of primary resistant tumors, with two dMMR models harboring among others inactivating NF1 mutations paired with either a potentially activating FGFR1 mutation (K659N, BoC104) or an HRAS mutation (R161H, BoC258) and biallelic inactivation of PTEN (BoC258). From the two primary resistant and MMR proficient PDX models, BoC215 harbored a mutation in the ERBB2 gene (S310F), previously described as activating . Taken together, in three of the five primary resistant models, WES revealed potential resistance-conferring mutations, not detected in our initial targeted sequencing approach for selecting quadruple WT tumors. Furthermore, dMMR tumors did not benefit from anti-EGFR therapy. Next, we sequenced all the available 31 SR tumors from the first 10 consecutively arising models either by WES (N = 29) or by targeted NGS (BoC71 C4, BoC106 C2) to identify resistance drivers not present in the corresponding control tumors (Fig. 3b, Additional file 5: Table S4). Furthermore, selected SR tumors from five models (BoC10 C3, BoC20 C9, BoC35 C1, BoC60 C5, BoC69 C5, and BoC106 C2) were analyzed for copy number variations (CNV). In one model (BoC71) with four SR tumors, we found two acquired KRAS mutations (G12V & G12C in BoC71 C3 & C4, respectively) and one MAP2K2 (Q60P) mutation (BoC71C1), while the fourth tumor (BoC71 C6) did not reveal any alterations known to drive SR. A similar genetic heterogeneity in SR conferring driver mutations was found in model BoC106, with one SR tumor (BoC106 C8) carrying a BRAF mutation (N518I, described as low-activity or kinase-dead BRAF mutant , and therefore not considered as a driver mutation herein), a second (BoC106 C2) developing an amplification (chr12:22615401-26669741) including the KRAS gene locus (Fig. 4a, b), while the third tumor (BoC106 C6) did not reveal an obvious driver gene alteration. Additional SR driver gene alterations were identified, including MAP2K1 (K57N; BoC237; 1 out of 4 tumors), KRAS (G12V; BoC209; 1 out of 3 tumors), and EGFR (S492R; BoC60 C5). We also observed a recurrent mutation in the tumor suppressor protein phosphatase 2 regulatory subunit B'alpha (PPP2R5A; R112L) in 4 SR tumor models (BoC10, 32, 209, and 237). This missense mutation in the beta subunit of the protein phosphatase 2A has been described in COSMIC and is classified as inactivating and highly deleterious in the Cancer Genome Interpreter (CGI). However, overexpression of both the wt and the R112L mutated version of PPP2R5A (Fig. 5a–c) had a similar negative impact on colon cancer cell growth and did not alter CET treatment sensitivity (Fig. 5d). Furthermore, pAKT and pERK levels remained at lower levels following EGF stimulation (at 6 h) in cells overexpressing the PPP2R5A mutation than in control and wt cells (Fig. 5e). This together with the generally low allele frequency (AF) of mutated PPP2R5A (< 35%) found herein, does not support that this mutation conveys PPP2R5A loss of function and is unlikely relevant for SR development in our tumors. In the remaining five models including 15 SR tumors, no potential driver gene alteration was identified. Additional mutations in genes with a low predicted driver gene status according to the Cancer Genome Interpreter (CGI) (mainly with tier 2 status and an AF below 20%) were identified in some of the SR tumors. They were not considered SR drivers due to the lack of proof of their functional relevance according to the current literature (Additional file 2: Fig. S7, Additional file 5: Table S4).
Clonal stability of KRAS G12V SR driver mutation independent of treatment pressure
One of the mutated tumors (BoC71 C3) with a KRAS mutation AF of 22 and 26% with BEAMing emulsion digital PCR (Fig. 6a) and WES, respectively, was selected to search for pre-existing KRAS G12V mutation by BEAMing analysis up to the 5th generation of treatment-naive BoC71 tumors. Even when using an analytical cut-off of 0.02%, no mutation was detected (Fig. 6a). This data is compatible with both de novo KRAS mutation developing under treatment pressure or a very rare pre-existing mutated cell clone below the detection limit of BEAMing as a source of the KRAS mutated cells in this anti-EGFR-treated tumor. Previously, ctDNA analyses suggested that KRAS driver mutations selected during anti-EGFR resistance development have a fitness disadvantage, as their frequency declines in ctDNA, following a hiatus from anti-EGFR treatment [17, 58]. BoC71 C3 tumor was chosen to test if a similar treatment-dependent alteration in the mutated KRAS allele frequency could be observed. We monitored the mutated KRAS AF with high sensitive droplet digital PCR (ddPCR) under treatment and treatment hiatus in BoC71 C3 following re-implantation. First, there was a striking stability of the mutated KRAS AF of around 33% between the initially developed SR tumor C3 and all subsequent tumors generated via re-implantation (note that the AF of around 33% in the ddPCR analysis is somewhat higher than in BEAMing and WES, suggesting an overestimation of the AF in the ddPCR for technical reasons). Second, a regimen where initial treatment was followed by a treatment pause and continued tumor growth had no influence on the AF of mutated KRAS (Fig. 6b, c). These data support a high stability of the proportion of wt to mutated cells.
Transcriptional reprogramming during SR development
The unexpected low prevalence of SR conferring driver mutations in our set of tumors (23%, 7/31), let us ask whether transcriptional changes may play a role in driving SR in some of the tumors. Therefore, global gene expression profiles of all 31 tumors with acquired SR were established via standard array technology. Comparing the profiles of SR tumors with the corresponding profiles of the CET-sensitive control tumors, collected at day 5 of treatment, identified numerous transcriptional changes in the SR tumors (Additional file 3: Table S2). First, we classified sensitive tumors as well as their corresponding SR tumor according to CRIS classification , a categorization system optimized for xenograft tumors (Additional file 2: Fig. S8). As expected, none of the untreated control tumors fell into CRIS-class A, which was associated with molecular alteration such as BRAF or KRAS mutations . Class switches between the CET-sensitive short-term treated tumors (5dC) and the corresponding SR tumors were with the exception of four SR tumors not prevalent. Furthermore, AREG and EREG, but not the other three known EGFR ligands, were shown to be positive modulators of primary anti-EGFR response and treatment outcome, both in patients and CRC PDX models [59,60,61]. These EGFR ligands are considered to induce a stronger EGFR dependency in colon tumors which renders them particularly sensitive to anti-EGFR treatment. Accordingly, AREG was expressed at significantly higher levels in control tumors of CET responders versus tumors of non-responders (Fig. 7a). A similar correlation could not be demonstrated for the remaining EGFR ligands. We also searched for overexpression of genes of the ERBB receptor family, the MET receptor, the RAS-RAF-MEK-ERK or the PI3K-AKT pathway, some of them were previously shown to be amplified in anti-EGFR SR tumors. In line with the demonstrated Chr. 12p12.1 amplification, BoC106 C2 showed a strong increase in KRAS expression (Fig. 4b). However, there was no indication of an amplification of the aforementioned genes in any of the remaining SR tumors, because gene expression levels did not vary by more than 2-fold between untreated controls (K) and resistant tumors (Additional file 2: Fig. S8). We also noted three PDX models (BoC10, 69 and 209, Additional file 2: Fig. S8) with a high basal IGF2 expression. In line with previous reports, all three models belonged to the SD response group .
Gene set enrichment analyses to uncover resistance pathways
Having excluded known mechanisms of resistance for the majority of our SR tumors, we next searched for classes of genes that are over-represented in the transcriptome of SR tumors and which may, therefore, have an association with the resistance phenotype via Gene Set Enrichment Analyses (GSEA) for “hallmarks of cancer.” In agreement with previous reports, showing that reactivation of the RAS pathway is key to anti-EGFR SR development, the “KRAS_SIGNALLING_UP” set was found to be enriched in SR tumors without a genetic driver alteration in five of ten tumor models (Fig. 7b; Additional file 3: Tables S2, Additional file 6: Table S5). Dual specificity phosphatase (DUSP6) is central in the negative feedback regulation of the KRAS signaling pathway and its expression level was therefore used as a surrogate marker gene for RAS pathway activity. In line with RAS signaling attenuation, we observed in all 10 models a 2.2- to 16-fold reduced DUSP6 expression during the initial CET treatment phase (5dC) as compared to baseline (untreated) tumors (Fig. 7c; Additional file 2: Fig. S9). Upon SR development, DUSP6 was significantly higher expressed in tumors with driver mutations (Fig. 7d). Importantly, DUSP6 expression levels rose above baseline expression in tumors with near clonal to clonal selected driver mutations (here AF ≥ 0.4) or amplifications, but remained below baseline in tumors with an AF ≤ 0.3 (Fig. 7c). Furthermore, tumors with enrichment of the “KRAS_SIGNALLING_UP” set were more likely to have a higher DUSP6 expression as compared to the group of tumors not showing this enrichment, corroborating the GSEA data (Fig. 7d). Among the identified and confirmed upregulated genes that are able to activate the KRAS signaling pathway, growth factors (GFs) and their receptors initiating an autocrine growth loop via FGF3, 9, 18, 19, and 20; FGFR1 and FGFR2; IGF1, IGF2, KIT, and PDGFA are according to our data likely playing a prominent role in acquired SR (Fig. 7e; Additional file 2: Fig. S8, S10, and S11, Additional file 3: Table S2, Additional file 7: Table S6). Further, differentially regulated candidate genes in the RAS pathway were the ras guanine nucleotide exchange factor (GEF) RASGRF1, the calcium-dependent phospholipid-binding protein PLA2G4A, and the downregulation of the RAS-GTPase activating protein RASAL1 (Fig. 7e; Additional file 2: Fig. S10, Additional file 3: Table S2, Additional file 7: Table S6).
Additional, potentially cooperating pathway alterations in five of the ten resistant tumors (Boc10, 60, 69, 106, and 237) could provide TGF-beta signaling inhibition via the WNT target gene BAMBI and/or the IRX3 and IRX5 genes (Additional file 2: Fig. S10, Additional file 7: Table S6). These genes have been shown to play a role in colon cancer development and may thus provide the cell with some growth advantage under anti-EGFR therapy, supporting resistance development [62,63,64].
A subgroup of three resistant tumors (BoC20, 60, and 69) showed markedly higher expressed genes involved in cancer stem cell signaling, such as OLFM4, CXCR4, ALDH1A1, or ALDH1A3. This was accompanied by a striking upregulation of defensins (DEFA5 and DEFA6) together with an activation of the Notch pathway via the ligand DLL1 in two of the three tumors (Additional file 2: Fig. S10, Additional file 7: Table S6). This gene expression pattern shows similarities to what has been described for colonic paneth-like secretory cells. These cells are known to play a role in stem cell maintenance and may thereby support resistance towards anti-EGFR antibodies [65, 66]. Interestingly, a strong reduction of HOXA5 gene expression was also observed in BoC60, a gene recently shown to be involved in intestinal stem cell fate control .
Paracrine support of SR from the mouse stroma
It is well established that mouse PDX models maintain the architecture of the original tumor and replace human stroma with mouse stroma within the first two generations of PDX , supporting that human epithelial and mouse stroma cells establish a functional communication. Therefore, we used high coverage 3’-RNA-seq to address whether, i.e., GF expression is altered in the stroma derived from SR tumors as compared to control tumors at day 5 of treatment. Indeed, a rise in GF expression, mainly Fgf- but also Egf-, Igf-, and Pdgf-family members, was frequently observed (Fig. 8a), suggesting an additive paracrine activity of the stroma cells towards SR development, at least in a subset of tumors. We used real-time cell analyses (RTCA) to test several recombinant mouse GFs (Fgf9, 10, and 18, Pdgfb and c, Igf2, Hbegf) for their ability to reduce CET response in a quadruple wt primary cell line derived from BoC20. Indeed, mouse GFs Fgf9, 18, and Hbegf strongly interfered with the growth-inhibiting activity of CET (Fig. 8b–h), with Fgf9 and 18 showing the highest homology at the amino acid level (> 98%) with the corresponding human GF among the tested GFs.
Baseline GF expression levels are linked to anti-EGFR Response
We noted that three treatment-naive PDX models (BoC10, 69, and 209) exhibited in addition to the expression of EGFR ligands a high expression of IGF2 as well as several FGFs, and PDGFs GF family members not observed in the other models (GFhigh) (Additional file 2: Fig. S8). Furthermore, all three models responded to CET treatment with stable disease. We hypothesized that the additional GF signaling is, among others, likely increasing the RAS pathway activity, not targeted by CET. Therefore, we tested for two models (BoC69 and 209) whether the addition of a MEK inhibitor would deepen the response and furthermore delay SR development. Indeed, we observed, that both PDX models could be shifted in the PR response group via this combination (Fig. 9a, b). Importantly, while three of five BoC69 tumors under chronic CET monotherapy gave rise to SR tumors, none of the seven tumors developed SR under the combination therapy (Fig. 9c). Of note, while prolonged trametinib monotherapy was well tolerated, its combination with CET led to an increase in toxicity with two of eight mice experiencing weight loss forcing us to terminate the experiment while in the remaining 6 mice body weight remained stable throughout the therapy (Additional file 2: Fig. S12). This underpins that careful dose-finding is required to balance therapeutic efficiency with toxicity for this combination.
Functional proof of selected secondary resistance mechanisms
In BoC71 C3, we addressed the KRAS G12V mutation, a well-known secondary resistance mechanism. This model was therefore treated with the combination of an anti-EGFR antibody (CET) and a MEK inhibitor (refametinib) to achieve optimal vertical inhibition. In line with previous publications, BoC71 C3 reached partial response after 4 weeks of treatment, proving the efficacy of MEK-EGFR co-targeting (Fig. 9d) [10, 69,70,71]. In BoC237, the identified MAP2K1 C171G mutation was addressed with the highly selective and ATP-competitive ERK1/2 inhibitor SCH772984. In contrast to previous reports in a BRAF mutated cell line , ERK-I monotherapy failed to control tumor growth at a dosing of 12.5 mg/kg bid whereas the combination of CET with SCH772984, albeit it did not reach our criteria for stable disease, somewhat controlled tumor growth, but was unable to induce tumor shrinkage (Fig. 9e). Doubling of the SCH772984 dose in combination with CET was poorly tolerated, and an initiated treatment trial was discontinued.
To show that the identified transcriptional changes are indeed conferring SR, we chose to therapeutically address in a proof-of-concept experiment two novel anti-EGFR secondary resistance mechanisms, which to date have not been described for CRC in this context, involving the cAMP and FGFR signaling pathways. Secondary resistant tumor BoC10 showed a strong upregulation of two family members of cAMP-specific phosphodiesterases (PDE4B and PDE4D) of 206 and 15-fold, respectively, compared to the 5-day controls. The phosphodiesterase family is known to degrade the secondary messenger adenosine and guanosine 3′,5′-cyclic monophosphates (cAMP and cGMP, respectively). PDE4D may support proliferation in prostate and lung cancer [73, 74], and cAMP is able to augment or suppress ERK activity in a cell type-dependent manner . This secondary resistant tumor was therefore treated with a combination treatment of CET and rolipram, an unselective PDE4 inhibitor. This treatment controlled the growth of this SR tumor for 2 weeks (Fig. 9f) supporting a previously undescribed role of phosphodiesterases in SR development for a subgroup of tumors. Since this combination was unable to sustain tumor growth control over a longer period of time (Fig. 9f), we searched for alternative treatment options. In this tumor, we also observed not only upregulation of FGF ligands (FGF19 and FGF20; 88- and 68-fold, respectively) but also a very strong upregulation of the phospholipase PLA2G4A (345-fold). These genes are able to activate among others the MEK-ERK signaling pathway [76, 77]. We therefore also tested the combination of CET and refametinib. We achieved disease control until treatment was stopped after 28 days (Fig. 9f). This indicates in agreement with the GSEA data that vertical inhibition via blocking KRAS signaling in combination with controlling the feedback activation of the EGFR receptor phosphorylation via CET was critical to control this tumor. Another tumor model BoC69 developed a remarkably high expression of FGF20 and FGF13 (59- and 18-fold, respectively). As an autocrine GF loop driven by FGF appeared to be the SR mechanism most likely in this tumor, we chose to test FGFR signaling blockade via the pan-FGFR inhibitor (FGFR-I) LY2874455. FGFR-I monotherapy achieved partial response which was surprisingly good and this was not improved by adding CET, suggesting a switch in GF dependencies for tumor growth from EGFR to FGFR signaling in this tumor (Fig. 9g).
In the tumor model BoC32, we observed in SR tumors an upregulation of IL-6 family members such as Ciliary Neurotrophic Factor Receptor (CNTFR) and Interleukin 6 receptor (IL6R) (Additional file 7: Table S6), known to be able to signal via Janus kinase (JAK)/(STAT) . Therefore, we performed in addition phosphoprotein Western blot analyses for STAT3. We found a strong increase in the STAT3 phosphorylation signal in the secondary resistant tumor in three of the four SR PDX tumors generated (Fig. 9h, i; Additional file 2: Fig. S13). Since clinically usable STAT3 inhibitors are still under development, STAT3 activated resistant tumor BoC32 C2 was treated in combination with CET and ruxolitinib, a JAK1/2 inhibitor or with the histone deacetylase (HDAC)-inhibitor vorinostat both of which have been shown to cause an inhibitory effect on STAT3 signaling [79, 80]. With both substances, we observed disease control, which was more pronounced with vorinostat than with ruxolitinib (Fig. 9j). Vorinostat was therefore tested for an extended treatment period. In line with the expected effect of both drugs, we found a sustained reduction of pSTAT3 following both treatments (Fig. 9h, i). In contrast to previous reports linking high pre-treatment pSTAT3 activity as well as the primary failure of CET to reduce STAT3 transcriptional activity to poor primary anti-EGFR response [81, 82], our data show that rewiring of signaling during chronic CET therapy can induce a high STAT3 activity status likely supporting SR development.
Secondary resistance is the main reason for acquired treatment failure to targeted therapies. Here, we investigated the feasibility to generate in vivo selected CET SR PDX tumors from CRC patients as a tissue source for in-depth molecular analysis of SR mechanisms. We found that the pre-selection of CRC tumors for their KRAS, NRAS, BRAF, and PI3K wildtype status efficiently selected anti-EGFR sensitive tumors in the majority of PDX models tested (15/21). Furthermore, SR tumors were established with a high success rate (80%) from CET-sensitive PDX models. This enabled us to generate altogether 31 SR tumors and thus the largest to date available CET SR PDX collection of in vivo selected SR tumors derived from primary patient tumor tissue. WES identified a high mutational load due to MMR gene inactivation only in the CET primary resistant cases, adding a note of caution for treating MMR-deficient tumors with CET. A recent single-center retrospective analysis supports this observation . Importantly, WES detected driver gene alterations only in 23% of 31 CET SR tumors, including acquired EGFR, KRAS, MEK1, and MEK2 mutations as well as one case with a KRAS amplification. The unexpected low frequency of driver mutations in SR let us ask whether transcriptome analyses may help to explain how driver mutation-negative PDX tumors acquired resistance. We observed a widespread transcriptional reprogramming of tumors during SR development. Using the CRIS classification system, we noted only a few occasions of class switches from sensitive to SR tumors, indicating that a consistent shift in cancer cell-intrinsic transcriptomic subtypes cannot be linked to SR development. In agreement with the concept that RAS pathway activation is critical for anti-EGFR SR development , GSEA and pathway analyses supported that in half of our SR tumors the observed transcriptional reprogramming leads to reactivation of the KRAS signaling pathway also accompanied by activation of the negative feedback via DUSP6. Transcriptional reprogramming and selection of tumor cells with an increased FGF-, IGF-, and PDGF- GF-family member gene expression inducing a signaling bypass to the CET EGFR blockade was identified as likely driver of SR in a subset of SR tumors. Constitutive overexpression of a few GFs such as IGF1, IGF2, and FGF9 prior to anti-EGFR therapy has previously been described to modulate the intrinsic resistance towards anti-EGFR therapy [61, 84, 85]. In addition, overexpression of downstream tyrosine kinase receptors including FGFR1, FGFR2, KIT, and PDGFRA was identified as a genetic primary resistance mechanism [13, 86]. Very few studies have to date analyzed pairs of pre-treatment and anti-EGFR SR tumors lacking driver gene mutations at the transcriptome level. One study from Woolston et al. very recently suggested that overexpression of GFs from the FGF- and TGF-β family are likely playing a role in SR development . Importantly, they proposed cancer-associated fibroblasts (CAFs) as the main source for the observed rise in GF expression, supported by the switch from the CMS2 to the CMS4 subtype in SR tumors, a subtype known to be enriched for fibroblasts . Thus, in his model, SR is driven by stromal remodeling and paracrine activity of a GF-rich stroma. The use of PDX models enabled us to differentiate between human and mouse cells as a source of the GFs and therefore differentiate between auto- and paracrine signaling, which is difficult to assess in primary human tumors. In contrast to Woolston et al., our data show for the first time that SR development is supported via transcriptional reprogramming of the tumor cell itself establishing autocrine GF signaling, which is further augmented by GF derived from tumor stroma. We also noted a subset of tumors with a higher baseline epithelial GF expression (GFhigh)—mainly FGF-, IGF-, and PDGF-family ligands—with some tumors developing SR by further increasing the expression of GF belonging to the same GF family and thus activating the RAS-RAF-ERK pathway independent of EGFR signaling. In a first pilot experiment, combined MEK and EGFR inhibition was not only able to improve primary response in GFhigh tumors but prevented SR development, suggesting that transcriptomic data of pre-treatment primary tumors may help to identify such GFhigh tumors, which may benefit from this combination therapy. Clearly, a more systematic analysis is needed to fully evaluate the impact of baseline GF expression on CET response and SR development.
In the only additional study available probing SR patient tumors via transcriptomic analyses, one tumor was found with enrichment in epithelial-to-mesenchymal transition (EMT) signature, and in another, a drop in the stromal infiltration signature accompanied by an increase in the immune infiltration signature, while in the third analyzed tumor, no signature enrichment was observed . In line with this data, we also observed enrichment of the EMT signature in a subgroup of our SR models (3/10), only (BoC10, 20 and 106, Additional file 6: Table S5), while Woolston et al. reported an EMT expression signature enrichment in the majority of his SR tumors (5/7). Future transcriptome analyses need to include higher numbers of cases in order to define the respective ratios regarding the transcriptional landscape in anti-EGFR SR and the pathways involved.
Lastly, in a proof-of-concept experiment to evaluate the functionality of identified targets, we selected prominent resistance mechanisms including activation of the cAMP-, FGF-, and STAT3-pathway for in vivo targeting. In all instances, we were able to regain control over tumor growth following targeted inhibition of these pathways, proving that the observed transcriptional reprogramming supports SR development. Of note, in these experiments, PDX tumor growth to reach the required volume for treatment initiation was done without the presence of CET. Nevertheless, in all instances, the anti-EGFR SR phenotype was readily established upon CET treatment, suggesting that the non-genetic resistance mechanisms discovered herein are likely stable and heritable.
Moreover, our strategy to generate in vivo SR models via chronic anti-EGFR treatment was able to successfully mimic the parallel evolution of multiple resistant lesions per patient tumor, represented in our setting by the multiple SR tumors evolving in the majority of our tested PDX models. Furthermore, heterogeneous molecular SR mechanisms among different SR tumors originating from the same primary tumor were also prevalent. This supports that SR PDX tumors are a good approximation to the human situation. Besides, the low frequency of resistance driver mutations found in this study is in good agreement with several analyses in primary CET SR tumor tissues reporting a low to moderate KRAS mutation frequency (9–21%) [4, 5, 8, 12,13,14]. However, it is in marked contrast to the current circulating tumor DNA (ctDNA) literature which suggests that pre-existing or de novo generated RAS and de novo EGFR mutations are the primary cause of SR development in up to 96% of CRCs [9, 15, 17, 58, 88,89,90,91,92,93]. The discrepancy between ctDNA and primary tissue mutation data has so far not been resolved. While mutation data from primary tissues are generally derived from one representative tissue biopsy only per patient, ctDNA analyses are considered to be an integrated molecular proxy to the overall patient’s SR lesions. However, ctDNA analyses are currently technically limited to the detection of mutations and in consequence unable to detect non-mutational SR mechanisms, as found herein. Our data support that a good proportion of the resistance driver mutation-negative tissue biopsies are from truly mutation-negative SR lesions and are likely driven by transcriptional reprogramming. Therefore, optimizing treatment by targeting the mutation pattern identified via ctDNA analysis is likely to fail in patients with progressing lesions lacking a resistance driver mutation.
The percentages of mutated KRAS alleles in ctDNA have repeatedly been reported to decline when anti-EGFR treatment is withdrawn [9, 17, 58], indicating in these cases a fitness disadvantage of KRAS mutated cells. We found a striking stability of the mutated KRAS AF after re-implantation of the SR tumor in all treatment schedules tested. Our data do not argue against the existence of such a fitness disadvantage but suggest that it may be relevant in a subgroup of tumors, only. In light of the discrepancy in driver mutation frequencies detected in ctDNA and primary tissue, a more direct confirmation of the decline of driver mutation AFs in SR primary tissues is warranted, not least because anti-EGFR re-treatment strategies have been suggested based on the ctDNA data. Considering that the observed decline of mutated KRAS alleles in ctDNA upon CET treatment cessation was not consistently accompanied with tumor shrinkage , another intriguing possibility may also in part explain a proportion of the observed discrepancy between ctDNA and primary tissue mutation data. It is well conceivable that driver gene activation, as well as transcriptional reprogramming, can arise in parallel in the same lesion and compete with each other or happily coexist.
The high frequency of SR tumors lacking a resistance mutation identified in our PDX models is compatible both with the presence of proliferating cancer cell subpopulations with a reversible drug-tolerant state and slow proliferating drug-tolerant persister (DTP) cells [94,95,96,97,98,99,100]. In the first model, anti-EGFR treatment in CRC would induce stepwise transcriptional reprogramming of pre-resistant cells into a stably resistant state . In the second model, slow cycling DTP cells initially survive anti-EGFR treatment and later start to undergo cell division and acquire de novo resistance mechanisms, including driver and non-driver resistance mechanisms in a subset of persister cells, as recently shown for a number of tumor types [94, 97,98,99,100]. In this scenario, dynamic chromatin remodeling processes inducing histone modifications at promoters and distal regulatory elements and other epigenomic alterations would be the most likely molecular mechanism inducing anti-EGFR SR in our PDX models. Clearly, unidentified genetic aberrations such as noncoding regulatory mutations, and complex genetic/chromosomal changes including disruption of insulated neighborhood boundaries may also account for some of the observed transcriptional changes. Our in vivo models provide the unique platform not only to study the underlying molecular mechanisms driving transcriptional reprogramming but also allows the search for vulnerabilities of the cancer cells to interfere with its resistance evolutionary path, ultimately leading to more sustained treatment success.
A difference between patient SR tumors and our PDX model is the fact that 9 of our 12 SR models were derived from patients with UICC Stage I-III tumors, who do not receive anti-EGFR treatment in the clinical setting. Activated EGFR pathway signaling is known to support the growth of a subgroup of CRCs independent of their clinical stage and therefore more important for establishing the primary cancer growth and less critical for the development of the metastatic disease. This notion is also supported by our observation of response towards anti-EGFR therapy in earlier tumor stage CRC PDX models. Therefore, anti-EGFR SR mechanisms are likely stage independent and the transcriptional SR mechanism discovered in our study in non-stage IV models should also translate into later stage CRCs. Another limitation of our study is the fact that none of the patients included received anti-EGFR treatment throughout the duration of the study, precluding a direct comparison of the data from the PDX with the patient’s tumor progressing under treatment. However, Woolston et al., who specifically analyzed metastatic CRC lesions progressing under anti-EGFR monotherapy recently published that they were, in agreement with our data, unable to find mutations driving resistance in 64% of their biopsies . Importantly, they used ultra-deep sequencing and included a substantial number of biopsies taken from lesions, which more than doubled their volume relative to their volume at the time of best response. This should have considerably reduced the risk that sequencing sensitivity or biopsy sampling errors are responsible for the detection failure of resistance mutations. Therefore, we are confident that our data is a true reflection of what can be found in the patient’s tumor.
In summary, we demonstrate that PDX models are a suitable tool to efficiently establish SR tumors. They develop intertumoral heterogenic genetic and non-genetic anti-EGFR SR mechanisms comparable to what is observed in human CRC. We show that transcriptional reprogramming is an important SR mechanism in CRC PDX models under chronic anti-EGFR treatment. While non-genetic resistance mechanisms, including transcriptional reprogramming, are increasingly being recognized in several cancer types , they are currently largely ignored in the context of anti-EGFR SR in CRC. Our data emphasize the need for analyses of SR tumor tissues at a multi-omics level for a more differentiated molecular understanding of anti-EGFR SR in CRC. Lastly, identified non-genetic SR mechanisms could be addressed by targeted treatment, proving that SR PDX models are a valid tool enabling iterative treatment optimization targeting SR mechanisms, not possible in the patient. Ultimately, this should help to tailor more informed treatment strategies addressing anti-EGFR SR in patients.
Availability of data and materials
Exome sequencing has been deposited into the NCBI BioProject database under the BioProject ID PRJNA596887 (https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA596887) . The RNA sequencing and Agilent array data have been deposited in the NCBI’s Gene Expression Omnibus (GEO) database (accession number GSE141861(https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141861)  and GSE140973 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE140973) , respectively). Exome sequencing used for CNV analyses has been deposited into the European Genome-phenome Archive (EGA) database under the EGA study ID EGAS00001005320 (https://ega-archive.org/studies/EGAS00001005320) .
Circulating tumor DNA
Epidermal growth factor receptor
Van Cutsem E, Kohne CH, Hitre E, Zaluski J, Chang Chien CR, Makhson A, et al. Cetuximab and chemotherapy as initial treatment for metastatic colorectal cancer. N Engl J Med. 2009;360(14):1408–17. https://doi.org/10.1056/NEJMoa0805019.
Bokemeyer C, Bondarenko I, Makhson A, Hartmann JT, Aparicio J, de Braud F, et al. Fluorouracil, leucovorin, and oxaliplatin with and without cetuximab in the first-line treatment of metastatic colorectal cancer. J Clin Oncol. 2009;27(5):663–71. https://doi.org/10.1200/JCO.2008.20.8397.
Heinemann V, von Weikersthal LF, Decker T, Kiani A, Kaiser F, Al-Batran SE, et al. FOLFIRI plus cetuximab or bevacizumab for advanced colorectal cancer: final survival and per-protocol analysis of FIRE-3, a randomised clinical trial. Br J Cancer. 2021;124(3):587–94. https://doi.org/10.1038/s41416-020-01140-9.
Montagut C, Dalmases A, Bellosillo B, Crespo M, Pairet S, Iglesias M, et al. Identification of a mutation in the extracellular domain of the Epidermal Growth Factor Receptor conferring cetuximab resistance in colorectal cancer. Nat Med. 2012;18(2):221–3. https://doi.org/10.1038/nm.2609.
Dienstmann R, Patnaik A, Garcia-Carbonero R, Cervantes A, Benavent M, Rosello S, et al. Safety and activity of the first-in-class Sym004 anti-EGFR antibody mixture in patients with refractory colorectal cancer. Cancer Discov. 2015;5(6):598–609. https://doi.org/10.1158/2159-8290.CD-14-1432.
Misale S, Yaeger R, Hobor S, Scala E, Janakiraman M, Liska D, et al. Emergence of KRAS mutations and acquired resistance to anti-EGFR therapy in colorectal cancer. Nature. 2012;486(7404):532–6. https://doi.org/10.1038/nature11156.
Arena S, Bellosillo B, Siravegna G, Martinez A, Canadas I, Lazzari L, et al. Emergence of multiple EGFR extracellular mutations during cetuximab treatment in colorectal cancer. Clin Cancer Res. 2015;21(9):2157–66. https://doi.org/10.1158/1078-0432.CCR-14-2821.
Siena S, Sartore-Bianchi A, Garcia-Carbonero R, Karthaus M, Smith D, Tabernero J, et al. Dynamic molecular analysis and clinical correlates of tumor evolution within a phase II trial of panitumumab-based therapy in metastatic colorectal cancer. Ann Oncol. 2018;29(1):119–26.
Khan KH, Cunningham D, Werner B, Vlachogiannis G, Spiteri I, Heide T, et al. Longitudinal liquid biopsy and mathematical modeling of clonal evolution forecast time to treatment failure in the PROSPECT-C Phase II Colorectal Cancer Clinical Trial. Cancer Discov. 2018;8(10):1270–85.
Misale S, Bozic I, Tong J, Peraza-Penton A, Lallo A, Baldi F, et al. Vertical suppression of the EGFR pathway prevents onset of resistance in colorectal cancers. Nat Commun. 2015;6(1):8305. https://doi.org/10.1038/ncomms9305.
Pietrantonio F, Vernieri C, Siravegna G, Mennitto A, Berenato R, Perrone F, et al. Heterogeneity of acquired resistance to anti-EGFR monoclonal antibodies in patients with metastatic colorectal cancer. Clin Cancer Res. 2017;23(10):2414–22.
Bray SM, Lee J, Kim ST, Hur JY, Ebert PJ, Calley JN, et al. Genomic characterization of intrinsic and acquired resistance to cetuximab in colorectal cancer patients. Sci Rep. 2019;9(1):15365.
Bertotti A, Papp E, Jones S, Adleff V, Anagnostou V, Lupo B, et al. The genomic landscape of response to EGFR blockade in colorectal cancer. Nature. 2015;526(7572):263–7. https://doi.org/10.1038/nature14969.
Woolston A, Khan K, Spain G, Barber LJ, Griffiths B, Gonzalez-Exposito R, et al. Genomic and Transcriptomic Determinants of Therapy Resistance and Immune Landscape Evolution during Anti-EGFR Treatment in Colorectal Cancer. Cancer Cell. 2019;36(1):35–50 e9. https://doi.org/10.1016/j.ccell.2019.05.013.
Diaz LA Jr, Williams RT, Wu J, Kinde I, Hecht JR, Berlin J, et al. The molecular evolution of acquired resistance to targeted EGFR blockade in colorectal cancers. Nature. 2012;486(7404):537–40. https://doi.org/10.1038/nature11219.
Yonesaka K, Zejnullahu K, Okamoto I, Satoh T, Cappuzzo F, Souglakos J, et al. Activation of ERBB2 signaling causes resistance to the EGFR-directed therapeutic antibody cetuximab. Sci Transl Med. 2011;3(99):99ra86.
Morelli MP, Overman MJ, Dasari A, Kazmi SMA, Mazard T, Vilar E, et al. Characterizing the patterns of clonal selection in circulating tumor DNA from patients with colorectal cancer refractory to anti-EGFR treatment. Ann Oncol. 2015;26(4):731–6. https://doi.org/10.1093/annonc/mdv005.
Bardelli A, Corso S, Bertotti A, Hobor S, Valtorta E, Siravegna G, et al. Amplification of the MET receptor drives resistance to anti-EGFR therapies in colorectal cancer. Cancer Discov. 2013;3(6):658–73.
Kim TW, Peeters M, Thomas A, Gibbs P, Hool K, Zhang J, et al. Impact of emergent Circulating tumor DNA RAS mutation in panitumumab-treated chemoresistant metastatic colorectal cancer. Clin Cancer Res. 2018;24(22):5602–9.
Bairoch A. The Cellosaurus, a Cell-Line Knowledge Resource. J Biomol Tech. 2018;29(2):25–38. https://doi.org/10.7171/jbt.18-2902-002.
Dirks WG, Drexler HG. STR DNA typing of human cell lines: detection of intra- and interspecies cross-contamination. In: Helgason CD, Miller CL, editors. Basic Cell Culture Protocols. Totowa, NJ: Humana Press; 2013. p. 27–38. https://doi.org/10.1007/978-1-62703-128-8_3.
Sturm M, Schroeder C, Bauer P. SeqPurge: highly-sensitive adapter trimming for paired-end NGS data. BMC Bioinformatics. 2016;17(1):208. https://doi.org/10.1186/s12859-016-1069-7.
Li H, Durbin R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doi.org/10.1093/bioinformatics/btp324.
Verbruggen S, Menschaert G. mQC: a post-mapping data exploration tool for ribosome profiling. Comput Methods Programs Biomed. 2019;181:104806. https://doi.org/10.1016/j.cmpb.2018.10.018.
Garrison E, Marth G. Haplotype-based variant detection from short-read sequencing. 2012. https://arxiv.org/abs/1207.3907.
McLaren W, Gil L, Hunt SE, Riat HS, Ritchie GR, Thormann A, et al. The Ensembl variant effect predictor. Genome Biol. 2016;17(1):122.
Lek M, Karczewski KJ, Minikel EV, Samocha KE, Banks E, Fennell T, et al. Analysis of protein-coding genetic variation in 60,706 humans. Nature. 2016;536(7616):285–91.
Karczewski KJ, Francioli LC, Tiao G, Cummings BB, Alfoldi J, Wang Q, et al. The mutational constraint spectrum quantified from variation in 141,456 humans. Nature. 2020;581(7809):434–43. https://doi.org/10.1038/s41586-020-2308-7.
The 1000 Genomes Project Consortium, Auton A, Brooks LD, Durbin RM, Garrison EP, Kang HM, et al. A global reference for human genetic variation. Nature. 2015;526(7571):68–74.
Tamborero D, Rubio-Perez C, Deu-Pons J, Schroeder MP, Vivancos A, Rovira A, et al. Cancer Genome Interpreter annotates the biological and clinical relevance of tumor alterations. Genome Med. 2018;10(1):25. https://doi.org/10.1186/s13073-018-0531-8.
Vangala D, Ladigan S, Liffers ST, Noseir S, Maghnouj A, Götze TM, et al. Secondary resistance to anti-EGFR therapy by transcriptional reprogramming in patient-derived colorectal cancer models BioProject PRJNA596887, NCBI Sequence Read Archive. https://www.ncbi.nlm.nih.gov/bioproject/?term=PRJNA596887.
Chen S, Zhou Y, Chen Y, Gu J. fastp: an ultra-fast all-in-one FASTQ preprocessor. Bioinformatics. 2018;34(17):i884–i90.
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. https://doi.org/10.1093/bioinformatics/bts635.
Robinson JT, Thorvaldsdottir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6. https://doi.org/10.1038/nbt.1754.
Smith T, Heger A, Sudbery I. UMI-tools: modeling sequencing errors in Unique Molecular Identifiers to improve quantification accuracy. Genome Res. 2017;27(3):491–9. https://doi.org/10.1101/gr.209601.116.
Liao Y, Smyth GK, Shi W. The Subread aligner: fast, accurate and scalable read mapping by seed-and-vote. Nucleic Acids Res. 2013;41(10):e108. https://doi.org/10.1093/nar/gkt214.
Vangala D, Ladigan S, Liffers ST, Noseir S, Maghnouj A, Götze TM, et al. Secondary resistance to anti-EGFR therapy by transcriptional reprogramming in patient-derived colorectal cancer models GSE141861, NCBI's Gene Expression Omnibus database. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE141861.
Navarro Gonzalez J, Zweig AS, Speir ML, Schmelter D, Rosenbloom Kate R, Raney BJ, et al. The UCSC Genome Browser database: 2021 update. Nucleic Acids Res. 2021;49(D1):D1046–D57.
Li H, Handsaker B, Wysoker A, Fennell T, Ruan J, Homer N, et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics. 2009;25(16):2078–9.
Picard Toolkit. Broad Institute, GitHub Repository https://broadinstitute.github.io/picard/; Broad Institute. 2019.
Vangala D, Ladigan S, Liffers ST, Noseir S, Maghnouj A, Götze TM, et al. Secondary resistance to anti-EGFR therapy by transcriptional reprogramming in patient-derived colorectal cancer models EGAS00001005320, European Genome-phenome Archive database. https://ega-archive.org/studies/EGAS00001005320.
Koboldt DC, Zhang Q, Larson DE, Shen D, McLellan MD, Lin L, et al. VarScan 2: somatic mutation and copy number alteration discovery in cancer by exome sequencing. Genome Res. 2012;22(3):568–76.
Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics. 2010;26(6):841–2. https://doi.org/10.1093/bioinformatics/btq033.
Wang J, Mullighan CG, Easton J, Roberts S, Heatley SL, Ma J, et al. CREST maps somatic structural variation in cancer genomes with base-pair resolution. Nat Methods. 2011;8(8):652–4. https://doi.org/10.1038/nmeth.1628.
Morris TJ, Butcher LM, Feber A, Teschendorff AE, Chakravarthy AR, Wojdacz TK, et al. ChAMP: 450 k Chip Analysis Methylation Pipeline. Bioinformatics. 2014;30(3):428–30. https://doi.org/10.1093/bioinformatics/btt684.
Vangala D, Ladigan S, Liffers ST, Noseir S, Maghnouj A, Götze TM, et al. Secondary resistance to anti-EGFR therapy by transcriptional reprogramming in patient-derived colorectal cancer models GSE140973, NCBI’s Gene Expression Omnibus database. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE140973.
Gentner B, Schira G, Giustacchini A, Amendola M, Brown BD, Ponzoni M, et al. Stable knockdown of microRNA in vivo by lentiviral vectors. Nat Methods. 2009;6(1):63–6. https://doi.org/10.1038/nmeth.1277.
Medina Diaz I, Nocon A, Mehnert DH, Fredebohm J, Diehl F, Holtrup F. Performance of Streck cfDNA blood collection tubes for liquid biopsy testing. Plos One. 2016;11(11):e0166354. https://doi.org/10.1371/journal.pone.0166354.
Rago C, Huso DL, Diehl F, Karim B, Liu G, Papadopoulos N, et al. Serial assessment of human tumor burdens in mice by the analysis of circulating DNA. Cancer Res. 2007;67(19):9364–70. https://doi.org/10.1158/0008-5472.CAN-07-0605.
Diehl F, Li M, He Y, Kinzler KW, Vogelstein B, Dressman D. BEAMing: single-molecule PCR on microparticles in water-in-oil emulsions. Nat Methods. 2006;3(7):551–9. https://doi.org/10.1038/nmeth898.
Pfaffl MW. A new mathematical model for relative quantification in real-time RT-PCR. Nucleic Acids Res. 2001;29(9):e45–445. https://doi.org/10.1093/nar/29.9.e45.
Subramanian A, Tamayo P, Mootha VK, Mukherjee S, Ebert BL, Gillette MA, et al. Gene set enrichment analysis: a knowledge-based approach for interpreting genome-wide expression profiles. Proc Natl Acad Sci U S A. 2005;102(43):15545–50. https://doi.org/10.1073/pnas.0506580102.
Mootha VK, Lindgren CM, Eriksson KF, Subramanian A, Sihag S, Lehar J, et al. PGC-1alpha-responsive genes involved in oxidative phosphorylation are coordinately downregulated in human diabetes. Nat Genet. 2003;34(3):267–73. https://doi.org/10.1038/ng1180.
Isella C, Brundu F, Bellomo SE, Galimi F, Zanella E, Porporato R, et al. Selective analysis of cancer-cell intrinsic transcriptional traits defines novel clinically relevant subtypes of colorectal cancer. Nat Commun. 2017;8(1):15107. https://doi.org/10.1038/ncomms15107.
Bertotti A, Migliardi G, Galimi F, Sassi F, Torti D, Isella C, et al. A molecularly annotated platform of patient-derived xenografts (“xenopatients”) identifies HER2 as an effective therapeutic target in cetuximab-resistant colorectal cancer. Cancer Discov. 2011;1(6):508–23. https://doi.org/10.1158/2159-8290.CD-11-0109.
Greulich H, Kaplan B, Mertins P, Chen TH, Tanaka KE, Yun CH, et al. Functional analysis of receptor tyrosine kinase mutations in lung cancer identifies oncogenic extracellular domain mutations of ERBB2. Proc Natl Acad Sci U S A. 2012;109(36):14476–81. https://doi.org/10.1073/pnas.1203201109.
Yao Z, Yaeger R, Rodrik-Outmezguine VS, Tao A, Torres NM, Chang MT, et al. Tumours with class 3 BRAF mutants are sensitive to the inhibition of activated RAS. Nature. 2017;548(7666):234–8.
Siravegna G, Mussolin B, Buscarino M, Corti G, Cassingena A, Crisafulli G, et al. Clonal evolution and resistance to EGFR blockade in the blood of colorectal cancer patients. Nat Med. 2015;21(7):795–801.
Khambata-Ford S, Garrett CR, Meropol NJ, Basik M, Harbison CT, Wu S, et al. Expression of epiregulin and amphiregulin and K-ras mutation status predict disease control in metastatic colorectal cancer patients treated with cetuximab. J Clin Oncol. 2007;25(22):3230–7. https://doi.org/10.1200/JCO.2006.10.5437.
Jacobs B, De Roock W, Piessevaux H, Van Oirbeek R, Biesmans B, De Schutter J, et al. Amphiregulin and epiregulin mRNA expression in primary tumors predicts outcome in metastatic colorectal cancer treated with cetuximab. J Clin Oncol. 2009;27(30):5068–74. https://doi.org/10.1200/JCO.2008.21.3744.
Zanella ER, Galimi F, Sassi F, Migliardi G, Cottino F, Leto SM, et al. IGF2 is an actionable target that identifies a distinct subpopulation of colorectal cancer patients with marginal response to anti-EGFR therapies. Sci Transl Med. 2015;7(272):272ra12.
Fritzmann J, Morkel M, Besser D, Budczies J, Kosel F, Brembeck FH, et al. A colorectal cancer expression profile that includes transforming growth factor beta inhibitor BAMBI predicts metastatic potential. Gastroenterology. 2009;137(1):165–75. https://doi.org/10.1053/j.gastro.2009.03.041.
Martorell O, Barriga FM, Merlos-Suarez A, Stephan-Otto Attolini C, Casanova J, Batlle E, et al. Iro/IRX transcription factors negatively regulate Dpp/TGF-beta pathway activity during intestinal tumorigenesis. EMBO Rep. 2014;15(11):1210–8.
Barry GS, Cheang MC, Chang HL, Kennecke HF. Genomic markers of panitumumab resistance including ERBB2/ HER2 in a phase II study of KRAS wild-type (wt) metastatic colorectal cancer (mCRC). Oncotarget. 2016;7(14):18953–64. https://doi.org/10.18632/oncotarget.8006.
Rothenberg ME, Nusse Y, Kalisky T, Lee JJ, Dalerba P, Scheeren F, et al. Identification of a cKit(+) colonic crypt base secretory cell that supports Lgr5(+) stem cells in mice. Gastroenterology. 2012;142(5):1195–205 e6.
Brodrick B, Vidrich A, Porter E, Bradley L, Buzan JM, Cohn SM. Fibroblast growth factor receptor-3 (FGFR-3) regulates expression of paneth cell lineage-specific genes in intestinal epithelial cells through both TCF4/beta-catenin-dependent and -independent signaling pathways. J Biol Chem. 2011;286(21):18515–25. https://doi.org/10.1074/jbc.M111.229252.
Ordonez-Moran P, Dafflon C, Imajo M, Nishida E, Huelsken J. HOXA5 counteracts stem cell traits by inhibiting Wnt signaling in colorectal cancer. Cancer Cell. 2015;28(6):815–29.
Blomme A, Van Simaeys G, Doumont G, Costanza B, Bellier J, Otaka Y, et al. Murine stroma adopts a human-like metabolic phenotype in the PDX model of colorectal cancer and liver metastases. Oncogene. 2018;37(9):1237–50.
Misale S, Arena S, Lamba S, Siravegna G, Lallo A, Hobor S, et al. Blockade of EGFR and MEK intercepts heterogeneous mechanisms of acquired resistance to anti-EGFR therapies in colorectal cancer. Sci Transl Med. 2014;6(224):224ra26.
Troiani T, Napolitano S, Vitagliano D, Morgillo F, Capasso A, Sforza V, et al. Primary and acquired resistance of colorectal cancer cells to anti-EGFR antibodies converge on MEK/ERK pathway activation and can be overcome by combined MEK/EGFR inhibition. Clin Cancer Res. 2014;20(14):3775–86. https://doi.org/10.1158/1078-0432.CCR-13-2181.
Troiani T, Napolitano S, Martini G, Martinelli E, Cardone C, Normanno N, et al. Maintenance Treatment with Cetuximab and BAY86-9766 Increases Antitumor Efficacy of Irinotecan plus Cetuximab in Human Colorectal Cancer Xenograft Models. Clin Cancer Res. 2015;21(18):4153–64.
Morris EJ, Jha S, Restaino CR, Dayananth P, Zhu H, Cooper A, et al. Discovery of a novel ERK inhibitor with activity in models of acquired resistance to BRAF and MEK inhibitors. Cancer Discov. 2013;3(7):742–50. https://doi.org/10.1158/2159-8290.CD-13-0070.
Rahrmann EP, Collier LS, Knutson TP, Doyal ME, Kuslak SL, Green LE, et al. Identification of PDE4D as a Proliferation Promoting Factor in Prostate Cancer Using a Sleeping Beauty Transposon-Based Somatic Mutagenesis Screen. Cancer Res. 2009;69(10):4388–97. https://doi.org/10.1158/0008-5472.CAN-08-3901.
Pullamsetti SS, Banat GA, Schmall A, Szibor M, Pomagruk D, Hanze J, et al. Phosphodiesterase-4 promotes proliferation and angiogenesis of lung cancer by crosstalk with HIF. Oncogene. 2013;32(9):1121–34. https://doi.org/10.1038/onc.2012.136.
Dumaz N, Marais R. Integrating signals between cAMP and the RAS/RAF/MEK/ERK signalling pathways. Based on the anniversary prize of the Gesellschaft fur Biochemie und Molekularbiologie Lecture delivered on 5 July 2003 at the Special FEBS Meeting in Brussels. Febs J. 2005;272(14):3491–504.
Hua S, Yao M, Vignarajan S, Witting P, Hejazi L, Gong Z, et al. Cytosolic phospholipase A2alpha sustains pAKT, pERK and AR levels in PTEN-null/mutated prostate cancer cells. Biochim Biophys Acta. 2013;1831(6):1146–57. https://doi.org/10.1016/j.bbalip.2013.02.011.
Desnoyers LR, Pai R, Ferrando RE, Hotzel K, Le T, Ross J, et al. Targeting FGF19 inhibits tumor growth in colon cancer xenograft and FGF19 transgenic hepatocellular carcinoma models. Oncogene. 2008;27(1):85–97. https://doi.org/10.1038/sj.onc.1210623.
Hu X, Zhao Y, He X, Li J, Wang T, Zhou W, et al. Ciliary neurotrophic factor receptor alpha subunit-modulated multiple downstream signaling pathways in hepatic cancer cell lines and their biological implications. Hepatology. 2008;47(4):1298–308. https://doi.org/10.1002/hep.22163.
Gupta M, Han JJ, Stenson M, Wellik L, Witzig TE. Regulation of STAT3 by histone deacetylase-3 in diffuse large B-cell lymphoma: implications for therapy. Leukemia. 2012;26(6):1356–64. https://doi.org/10.1038/leu.2011.340.
Lee HJ, Zhuang G, Cao Y, Du P, Kim HJ, Settleman J. Drug resistance via feedback activation of Stat3 in oncogene-addicted cancer cells. Cancer Cell. 2014;26(2):207–21. https://doi.org/10.1016/j.ccr.2014.05.019.
Dobi E, Monnien F, Kim S, Ivanaj A, N'Guyen T, Demarchi M, et al. Impact of STAT3 phosphorylation on the clinical effectiveness of anti-EGFR-based therapy in patients with metastatic colorectal cancer. Clin Colorectal Cancer. 2013;12(1):28–36. https://doi.org/10.1016/j.clcc.2012.09.002.
Ung N, Putoczki TL, Stylli SS, Ng I, Mariadason JM, Chan TA, et al. Anti-EGFR therapeutic efficacy correlates directly with inhibition of STAT3 activity. Cancer Biol Ther. 2014;15(5):623–32. https://doi.org/10.4161/cbt.28179.
Yu Y, Ying J, Zhang W, Li W, Sun Y, Yang L, et al. Outcome of chemotherapy with or without targeted agents in metastatic colorectal cancer patients with deficient DNA mismatch repair: a single center, cohort study. Asia Pac J Clin Oncol. 2019;15(3):128–35. https://doi.org/10.1111/ajco.13130.
Scartozzi M, Mandolesi A, Giampieri R, Pierantoni C, Loupakis F, Zaniboni A, et al. Insulin-like growth factor 1 expression correlates with clinical outcome in K-RAS wild type colorectal cancer patients treated with cetuximab and irinotecan. Int J Cancer. 2010;127(8):1941–7.
Mizukami T, Togashi Y, Naruki S, Banno E, Terashima M, de Velasco MA, et al. Significance of FGF9 gene in resistance to anti-EGFR therapies targeting colorectal cancer: A subset of colorectal cancer patients with FGF9 upregulation may be resistant to anti-EGFR therapies. Mol Carcinog. 2017;56(1):106–17.
Medico E, Russo M, Picco G, Cancelliere C, Valtorta E, Corti G, et al. The molecular landscape of colorectal cancer cell lines unveils clinically actionable kinase targets. Nat Commun. 2015;6(1):7002. https://doi.org/10.1038/ncomms8002.
Becht E, de Reynies A, Giraldo NA, Pilati C, Buttard B, Lacroix L, et al. Immune and stromal classification of colorectal cancer is associated with molecular subtypes and relevant for precision immunotherapy. Clin Cancer Res. 2016;22(16):4057–66. https://doi.org/10.1158/1078-0432.CCR-15-2879.
Thierry AR, Pastor B, Jiang ZQ, Katsiampoura AD, Parseghian C, Loree JM, et al. Circulating DNA demonstrates convergent evolution and common resistance mechanisms during treatment of colorectal cancer. Clin Cancer Res. 2017;23(16):4578–91.
Strickler JH, Loree JM, Ahronian LG, Parikh AR, Niedzwiecki D, Pereira AAL, et al. Genomic landscape of cell-free DNA in patients with colorectal cancer. Cancer Discov. 2018;8(2):164–73.
Bettegowda C, Sausen M, Leary RJ, Kinde I, Wang Y, Agrawal N, et al. Detection of circulating tumor DNA in early- and late-stage human malignancies. Sci Transl Med. 2014;6(224):224ra24.
Van Emburgh BO, Arena S, Siravegna G, Lazzari L, Crisafulli G, Corti G, et al. Acquired RAS or EGFR mutations and duration of response to EGFR blockade in colorectal cancer. Nat Commun. 2016;7:13665.
Yamada T, Matsuda A, Takahashi G, Iwai T, Takeda K, Ueda K, et al. Emerging RAS, BRAF, and EGFR mutations in cell-free DNA of metastatic colorectal patients are associated with both primary and secondary resistance to first-line anti-EGFR therapy. Int J Clin Oncol. 2020;25(8):1523–32.
Esposito C, Rachiglio AM, La Porta ML, Sacco A, Roma C, Iannaccone A, et al. The S492R EGFR ectodomain mutation is never detected in KRAS wild-type colorectal carcinoma before exposure to EGFR monoclonal antibodies. Cancer Biol Ther. 2013;14(12):1143–6.
Ramirez M, Rajaram S, Steininger RJ, Osipchuk D, Roth MA, Morinishi LS, et al. Diverse drug-resistance mechanisms can emerge from drug-tolerant cancer persister cells. Nat Commun. 2016;7(1):10690. https://doi.org/10.1038/ncomms10690.
Shaffer SM, Dunagin MC, Torborg SR, Torre EA, Emert B, Krepler C, et al. Rare cell variability and drug-induced reprogramming as a mode of cancer drug resistance. Nature. 2017;546(7658):431–5. https://doi.org/10.1038/nature22794.
Sharma SV, Lee DY, Li B, Quinlan MP, Takahashi F, Maheswaran S, et al. A chromatin-mediated reversible drug-tolerant state in cancer cell subpopulations. Cell. 2010;141(1):69–80. https://doi.org/10.1016/j.cell.2010.02.027.
Knoechel B, Roderick JE, Williamson KE, Zhu J, Lohr JG, Cotton MJ, et al. An epigenetic mechanism of resistance to targeted therapy in T cell acute lymphoblastic leukemia. Nat Genet. 2014;46(4):364–70.
Liau BB, Sievers C, Donohue LK, Gillespie SM, Flavahan WA, Miller TE, et al. Adaptive chromatin remodeling drives glioblastoma stem cell plasticity and drug tolerance. Cell Stem Cell. 2017;20(2):233–46 e7. https://doi.org/10.1016/j.stem.2016.11.003.
Risom T, Langer EM, Chapman MP, Rantala J, Fields AJ, Boniface C, et al. Differentiation-state plasticity is a targetable resistance mechanism in basal-like breast cancer. Nat Commun. 2018;9(1):3815. https://doi.org/10.1038/s41467-018-05729-w.
Hata AN, Niederst MJ, Archibald HL, Gomez-Caraballo M, Siddiqui FM, Mulvey HE, et al. Tumor cells can follow distinct evolutionary paths to become resistant to epidermal growth factor receptor inhibition. Nat Med. 2016;22(3):262–9. https://doi.org/10.1038/nm.4040.
Marine JC, Dawson SJ, Dawson MA. Non-genetic mechanisms of therapeutic resistance in cancer. Nat Rev Cancer. 2020;20(12):743–56. https://doi.org/10.1038/s41568-020-00302-4.
We thank F.X. Real for support and critical discussion. We thank Merck KGaA who supported the initial anti-EGFR sensitivity tests with providing CET (Erbitux) and the DKFZ-Heidelberg Center for Personalized Oncology (DKFZ-HIPO)/NCT Molecular Diagnostics for whole-exome sequencing through HIPO-B012. We thank the Omics IT and Data Management Core Facility (ODCF), the Genomics and Proteomics Core Facility (GPCF) and the Sample Preparation Laboratory (SPL) of the German Cancer Research Center (DKFZ, Heidelberg) and the Molecular Diagnostics Program of the National Center for Tumor Diseases (NCT) for excellent technical assistance.
The work was supported by grants from the German Cancer Aid (70111971) to A. Jung, J. Werner, V. Heinemann, S. Stintzing, U. Mansmann, A. Tannapfel., C.R. Ball, S. Fröhling, H. Glimm, S.A. Hahn, and W. Schmiegel and from the Ministry of Science, North Rhine-Westphalia, Germany (PURE/PRODI) to W. Schmiegel. J.T. Siveke is supported by funding from the German Cancer Consortium (DKTK), the German Cancer Aid (grant no. 70112505 (PIPAC) and grant no. 70113834 (PREDICT-PACA)), and the German Research Foundation (DFG) through grant SI1549/3-1 (Clinical Research Unity KFO337). This work was supported by the German Cancer Research Center-Heidelberg Center for Personalized Oncology (DKFZ-HIPO) and the Molecular Diagnostics Program of the National Center for Tumor Diseases (NCT) in Heidelberg. Open Access funding enabled and organized by Projekt DEAL.
Ethics approval and consent to participate
Informed and written consent was obtained from all patients. The study was approved by the Ruhr University Bochum Ethics Committee, Registry no. 3841-10 & 16-5792 and the Ludwig-Maximilians University Munich Ethics Committee, project no. 131-16. The research conforms to the principles of the Helsinki Declaration.
All animal experiments were approved by the local authorities (84-02.04.2012.A360 & 84-02.04.2015.A135) and performed in accordance with the guidelines for Ethical Conduct in the Care and Use of Animals.
Consent for publication
D.V. received speaker’s honoraria from Roche and Falk foundation and travel grants from Celgene and Gilead. D.L.E. is an employee of Sysmex Inostics. J.K. has had a consulting or advisory role, received honoraria and/or travel/accommodation expenses funding from AstraZeneca, Novartis, Quality Initiative in Pathology (QuIP), and Roche. M.P. has received consulting fees/ honoraria and has served as a speaker or advisory board member for Amgen, Merck Serono, Roche, Lilly, Sanofi-Aventis, Baxalta, Celgene, MCI München, Novartis, Alexion, Janssen-Cilag, MSD, BMS, Abbvie, Kite, and Gilead. S.S. has received honoraria as a speaker or advisory board member for Amgen, Bayer, Lilly, Merck KGaA, Sanofi, Takeda, Taiho, Pierre-Fabre, Roche, Samsung, and MSD. SF has had a consulting or advisory role, received honoraria, research funding, and/or travel/accommodation expenses funding from the following for-profit companies: Bayer, Roche, Amgen, Eli Lilly, PharmaMar, AstraZeneca, and Pfizer. J.T.S. received research funding from BMS, Celgene, and Roche, consulting and personal fees from AstraZeneca, Baxalta, BMS, Celgene, Immunocore, Lilly, Novartis, Roche, and Shire, holds ownership in FAPI Holding (< 3%) and Pharma15 (< 3%) and is a member of the Board of Directors for Pharma15 outside the submitted work. The remaining authors declare that they have no competing interests. We thank Merck KGaA who supported the initial anti-EGFR sensitivity tests with providing CET (Erbitux).
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
TruSeq Amplicon - Cancer Panel (Illumina).
Pdf file with all supplementary figures (Fig. S1-S13) with corresponding figure legends.
Growth curves of all chronically cetuximab treated tumors, GSEA snapshots of enrichment results, and summary of gene expression data for each individual xenograft tumor analyzed.
Tables of SNVs identified in primary resistant PDX models.
Summary of CGI predicted driver mutations and tables of SNVs identified in secondary resistant PDX models.
Selected gene sets enriched in phenotype “secondary resistant” (p ≤ 0.05).
Selected genes with potential function in SR for the 10 models tested.
About this article
Cite this article
Vangala, D., Ladigan, S., Liffers, S.T. et al. Secondary resistance to anti-EGFR therapy by transcriptional reprogramming in patient-derived colorectal cancer models. Genome Med 13, 116 (2021). https://doi.org/10.1186/s13073-021-00926-7