An embryonic diapause-like adaptation with suppressed Myc activity enables tumor treatment persistence

Lead contact

Further information and requests for resources and reagents should be directed to and will be fulfilled by the Lead Contact, Constantine S. Mitsiades (constantine_mitsiades@dfci.harvard.edu).

Materials availability

The organoid models and plasmids generated in this study are available upon request to the Lead Contact.

Data and code availability

The data retrieval, generation and analysis are described in the Methods section. Microarray and RNA-seq foldchange and counts tables are available as Table S1 (log₂FC in drug-refractory organoids and PDXs, deseq without beta correction), Table S2 (read counts data from our organoids and PDX samples), and Table S3 (processed count data from the Boroviak et al. dataset through our pipeline). The raw sequencing data (RNA-seq) have been deposited to the Gene Expression Omnibus (GSE162285). RNA-seq from Boroviak et al., Scogniamiglio et al., and Bulut-Karslioglu et al., are available under the accession numbers E-MTAB-2958, GSE74337 and GSE81285, respectively. Custom scripts and pipelines used for data pre-processing are described in the method section.

Cell Lines

MDAMB-231, MCF-7, and ZR75–1 cell lines were purchased from ATTC. Lenti-X-293T cells were purchased from Clontech. The cells were routinely maintained in DMEM High Glucose, 4.5g/l D-Glucose, 110mg/l Sodium Pyruvate, 10% FBS, 10 i.u./ml penicillin and 10μg/ml streptomycin and passaged using trypsin. Authentication of the cell lines was performed by short tandem microsatellite repeat analysis. Before their use in experiments, the cells were tested for mycoplasma using a combined PCR-ELISA kit (Roche Diagnostics).

Patient-derived models

The clinical characteristics of the patients are described in Suppl. Figure S1D. The PrCa patient-derived cells for the 3D-PDO models MSK-PCa1, MSK-PCa2, MSK-PCa3, MSK-PCa7 and PM154 were obtained at the Memorial Sloan Kettering Cancer Center and were expanded similarly to previous studies (Gao et al., 2014). The BrCa patient-derived samples HCI002, HCI003, HCI005, HCI007, HCI009 and HCI011 were initially engrafted as PDXs at the Huntsman Cancer Institute, as previously described (DeRose et al., 2011). Small PDX tissue fragments (1–2mm) were enzymatically processed for 1h at 37°C and the released floating malignant cells were collected by centrifugation, resuspended in matrigel, plated in 6-well plates and expanded as 3-D organoids, as previously described (Gao et al., 2014). The 3D-PDO cultures were maintained in Advanced DMEM/F12 supplemented with 3,151 ml/L D-glucose, 110 mg/L sodium pyruvate, EGF (50 ng/mL), FGF-10 (10 ng/mL), FGF-2 (1 ng/mL), Nicotinamide (10 mM), A83–01 (0.5 μM), SB202190 (10 μM), Y-27632 (10 μM), B27 (1X), N-Acetyl-L-Cysteine (1.25 mM), GlutaMAX (2 mM), HEPES (10 mM), Primocin (0.5μg/ml), R-Spondin (10% v/v), Noggin (10% v/v), 0.1nM estradiol (for the BrCa 3D-PDOs only), 0.1nM DHT (for the PrCa 3D-PDOs only).

Reagents

PBS (Corning), DMEM High Glucose (Life Technologies), Advanced DMEM/F12 (Life Technologies), Opti-MEM (Life Technologies), trypsin 0.05% (Mediatech), penicillin/streptomycin 1X (Mediatech), matrigel (Corning), rat tail collagen (Corning), Lipofectamine 2000 (Life Technologies), Trypan blue (Sigma Aldrich), polybrene (Santa Cruz Biotechnology), estradiol (Sigma Aldrich), dehydrotestosterone (Sigma Aldrich), fulvestrant (Sigma Aldrich), vinorelbine (Sigma Aldrich), docetaxel (LC Laboratories), epirubicin (SelleckChem), carfilzomib (SelleckChem), EGF (Miltenyi Biotec), FGF-10 (R & D Systems), FGF-2 (Peprotech), nicotinamide (Acros Organics), A83–01 (R & D Systems), SB202190 (Sigma Aldrich), Y-27632 (Sigma Aldrich), B27 (Life Technologies), N-Acetyl-L-Cysteine (Sigma Aldrich), Glutamax (Life Technologies), HEPES (Life Technologies), Primocin (Invivogen), monoclonal mouse anti-human Ki-67 antibody clone MIB1 (Dako), monoclonal mouse anti human ER antibody clone 1D5 (Dako), monoclonal mouse anti-cMyc antibody clone N-262 (Santa Cruz Bio.). R-spondin and Noggin were generated by collection of conditioned media from HEK293T cells transfected with the respective expression plasmids, as previously described (Gao et al., 2014). BET inhibitor JQ1 and CDK9 degrader ZZ1–33 (previously referred to as THAL-SNS-032) and inhibitors NVP2 and SNS032 were provided by the laboratories of Drs. Nathanael Gray and James Bradner. The chemical library panel of FDA-approved oncology drugs (100 compounds) was a kind gift from the Developmental Therapeutics Program of the NCI.

3D cultures

For simplicity, in this manuscript we use the term “organoids” to refer to multicellular 3-D structures formed in gels by cancer cells originating from either established cell lines or patient-derived tumors. We recognize that in the literature, 3D structures formed by cancer cell lines are often referred to as “spheroids”, while the term “organoid” is frequently used, especially in recent years, for material that was recently derived from patients. The use of the term “organoid” for 3D cultures of both cell line and patient-derived samples also reflects that in this study and others (i) many cell lines do not form typical “spheroid” structures in 3D, and therefore the term “organoid” is more reflective and inclusive; (ii) other cell lines form in 3D cultures structures that are indistinguishable with those from patient-derived “organoids”.

The cells were suspended in 80% matrigel as previously described (Dhimolea et al., 2010). The gels were casted in either 12-well (for histological and molecular analyses) or multi-well plates (for cell viability assays; described below). The gels were incubated for 1h at 37°C to allow for solidification, supplemented with medium, and used for further downstream applications (see other sections of Methods). To extract the cells from 3-D organoids, the gels were first disrupted mechanically in cold PBS; next, the organoids were isolated by centrifugation and incubated in trypsin for 20–30 min to generate single-cell solution.

Whole exome sequence analysis of patient-derived cancer models

A total of 14×2 paired-end sequencing samples (Illumina HiSeq 2000) were processed for prostate cancer primary tumor (Primary tumor, Xenograft and Organoid samples for MSK-PCa1, MSK-PCa2, MSK-PCa3, MSK-PCa7). A total of 31×2 paired-end sequencing samples (Illumina HiSeq 2000) were processed for breast cancer primary tumors, xenograft and organoid models (HCI-002, HCI-003, HCI-005, HCI-007, HCI-009, and HCI-011). The exome-seq data was generated using Agilent Sureselect per manufacturer’s instructions. We aligned reads to the human genome (hg19 from the UCSC) using the BWA (Burrows-Wheeler alignment) application (Li and Durbin, 2009). We applied the Piccard and the GATK (GenomeAnalysisTK-3.6) tools (McKenna et al., 2010) for processing BAM files of aligned reads generated by the bwa application for each sample. For the analysis we followed the GATK best practices recommendations. In detail, reads in the BAM files were sorted and duplicated reads removed. The alignments were realigned around indels and quality scores were recalibrated based on GATK. We processed BAM files using the picard (picard-1.138) tools ReorderSam, SortSam, AddOrReplaceReadGroups and MarkDuplicates. The alignment recalibration was performed using the GATK tools RealignerTargetCreator, IndelRealigner, BaseRecalibrator and PrintReads. Known SNP for the base recalibration were defined based on the NCBI dbSNP database (dbsnp_138.hg19.vcf). We used samtools for the indexing of BAM files. The identification of single nucleotide polymorphisms (SNPs) was performed based on GATK tools UnifiedGenotyper, VariantFiltration and Mutect2. For Mutect2 we filter out known human variants dbSNP138 and retain SNPs of known cancer associated mutations from the COSMIC database. Unreported SNPs were annotated using SnpEff and variantAnnotator from GATK.

Copy number variation was estimated based on GATK4 tools PadTargets, CalculateTargetCoverage, CombineReadCounts, CreatePanelOfNormals, NormalizeSomaticReadCounts, PerformSegmentation and CallSegments. The panel of normal controls (PON) were generated from blood samples. For the prostate cancer dataset we used 3 paired-end sequencing samples (MSK-PCa1_BM1, MSK-PCa2_BM5 and MSK-PCa3_ST1). For the breast cancer samples normal controls we used 3 paired-end sequencing samples (HCI-002_0900570-B, HCI-003_0903293-B, HCI-005_1007496-B). The sequence processing procedures were performed on the Orchestra High Performance Compute Cluster (HPC) at Harvard Medical School. The Orchestra HPC NIH supported shared facility is partially provided through grant NCRR 1S10RR028832–01.

Persister vs Vehicle exome sequencing analysis

The fastq files were processed as described in “Whole exome sequence analysis” and alignments were recalibrated based on GATK best practices. The identification of single nucleotide polymorphisms (SNPs) was performed based on GATK tools Mutect2 and FilterMutectCalls using a normal female sample as comparator. SNPs were annotated using Oncotator (Ramos et al., 2015) and we considered only non-synonymous SNPs causing a change in protein sequence. We used bam-readcount to extract the nucleotide frequencies for each called mutation form the “PASS” filtered mutect2 output. The Tumor allele was defined from Oncotator annotated SNP calls and we extracted the total readcount and the tumor allele readcount from the respective BAM alignment files (Persister and Vehicle) using bam-readcount (https://github.com/genome/bam-readcount). The tumor allele frequency was estimated by dividing the number of reads encoding the tumor allele by the total number of reads for the respective loci. In order to exclude low read coverage loci we considered only loci for the HCI002 with a 50x coverage in each of both the Persister and untreated vehicle sample. For the MDAMB-231 Persistor and for the untreated vehicle sample we considered a 200x read coverage in the loci in each of both samples.

Cell viability assays

Cell line MCF-7, ZR-75–1, MDAMB-231, and 3D-PDO models HCI002, HCI003, HCI007, HCI009, HCI011, MSK-PCa1, MSK-PCa2, MSK-PCa3, MSK-PCa7 and PMI154 were stably transfected with lentiviral vectors expressing the luciferase gene, which allows the use of the emitted bioluminescence to longitudinally and non-disruptively measure cancer cell viability, as previously described(McMillin et al., 2012). The effects of chemotherapeutics, JQ1 and CDK9 inhibitors in cell lines or patient-derived organoids were also confirmed using CellTiter-Glo (CTG) assay (Promega) cell viability assays following manufacturer’s instructions. The assay results (quadruplicates) were analyzed in Excel and plotted in GraphPad Prism (error bars represent standard error of the mean). Two in vitro culture configurations were applied in these studies:

2-D culture cell viability assays: cells were seeded into 384-well plates (500–2×10³ cells/well), in supplemented media (50uL-100uL/well) and incubated for 24h prior to addition of compound/s at indicated concentrations.

3-D organoids cell viability assays: the cells were plated in matrigel similarly to single-cell format. After plating, the cells were allowed to grow for 7–10 days into 3-D organoids. Once the organoid structures were formed (confirmed by microscope observation), the medium was replenished, and the organoids were incubated for 24 more hours before adding the compound(s) at indicated concentrations.

During prolonged drug-response assays the medium and tested compounds were replenished at regular intervals. At each indicated time point (1–21 days) either luciferin (for the luciferase-positive cells; concentration per manufacturer’s instructions) or CTG reagent (per manufacturer’s instructions) was added to each well, and the plates were read using a microplate reader (BioTek Synergy 2). The linearity of bioluminescence assay in 3-D cultures was estimated by measuring the signal of known cell numbers in 3-D cultures (data not shown).

EDL organoids regrowth suppression experiment

MDAMB-231 3-D organoids were exposed to the indicated chemotherapeutic agents for 96h to generate persister organoids (phenotypically and morphologically similar to Figures 1B, 1D, respectively); next, the gels were washed in cold PBS and the persister organoids were extracted and re-embedded in new gels an incubated in the presence of JQ1 (200nM), INK128 (100nM) or DMSO-control. The MDAMB-231 cells were extracted for the newly-formed gels at indicated time points using cold PBS and, when necessary, brief incubation with trypsin at 37°C to generate single-cell solution. The growth rates of persister organoids after drug washout was estimated by cell counting (Countess automated cell counter) at indicated time points. Cell counting measurements were also confirmed by bioluminescence measurements.

PRISM-based phenotypic studies in pools of barcoded cancer cell lines

PRISM is a method that allows pooled screening of mixtures of cancer cell lines by labeling each cell line with 24-nucleotide barcodes and has been previously described in detail (Yu et al., 2016). Briefly, we used 253 adherent cancer cells lines that stably expressed DNA barcode sequences with limited sequence homology to the human genome. For the 2-D cultures, the cells were seeded in 25cm flasks (100×10³/flask in experimental triplicates) and incubated in supplemented medium for 24h before adding docetaxel (100nM) or DMSO-control. For the 3-D cultures, the cells were suspended in matrigel and dispensed in 10cm dishes (100×10³ in 1ml matrigel/dish); let at 37°C for 37h to solidify and supplemented with culture medium. The gels were incubated at 37°C for 8 days to allow for 3-D organoid structures formation (confirmed by microscope observation), followed by medium replenishment and addition of docetaxel (100nM) or DMSO-control the next day. Both 2-D and 3-D cultures were harvested 5, 10 and 15 days after the start of treatment. The cellular DNA was isolated using standard DNA extraction kit (Qiagen) and sequenced to estimate the copies of the specific barcodes corresponding to distinct cell lines, as previously described(Yu et al., 2016). The comparison of the normalized counts for each barcode between treated conditions and the respective time-point DMSO-controls, was used to estimate the docetaxel-induced reduction of viability for the corresponding cell line.

Barcode-mediated clonal tracking experiment

The Crop-seq-opti vector (available from Addgene 106280) was modified by replacing the Puro resistance marker with Vex-GFP using restriction cloning. We cloned the barcoded gRNA library on the modified Crop-seq-opti_Vex vector using used the BsmBI restriction site. The double stranded 58 bp barcoded gRNA insert library was generated by performing a PCR extension reaction with the following primer oligos:

N20_gRNA GAGCCTCGTCTCCCACCGNNNNNNNNNNNNNNNNNNNNGTTTTGAGACGCATGCTGCA

and gRNA_RevExt TGCAGCATGCGTCTCAAAAC

as follows: 98°C for 2 min, 10x (65°C for 30 s, 72°C for 10s), 72°C for 2 min, hold at 4°C. The double stranded barcode-sgRNA oligo was purified using a QIAquick PCR Purification Kit (Qiagen, #28104). The double-stranded product contains two BsmBI sites that, upon digestion, generate complimentary overhangs for ligation into the Crop-seq-opti_Vex. Assembling of the double-stranded barcode-sgRNA insert into the Crop-seq-opti_Vex vector was done using a Golden Gate assembly reaction with at a molar ratio of 1:5 and cycled 100x (42°C for 2 min and 16°C for 5 min). The assembled plasmid was then purified using the DNA Clean & Concentrator™ kit (Zymo, #D4033) and used to transform Endura electrocompetent cells (Lucigen, #60242–2). Transformed bacteria were incubated overnight at 37°C plasmid DNA was extracted using a Qiagen Plasmid Plus Midi Kit (Qiagen #12943).

Lentiviral preps were generated by transfecting Lenti-X 293T cells with the Crop-seq barcode library and packaging plasmids psPAX2 and MD2.G. The MDAMB-231 cells were transfected with the generated viral preps at low MOI (~10% transfection efficiency; assessed by flow cytometry) to ensure that most cells receive one single barcode; flow cytometry analysis indictaed 6% transfection rate. Approximately 125 thousand EGFP+ cells (roughly representing equivalent number of barcodes) were collected by flow cytometry and expanded in culture. Prior to inoculation in animals, a pellet of 2 million cells was frozen to assess the baseline distribution of barcodes (T0). The rest of the cells were inoculated in 15 SCID mice (2 million cells per animal). When tumors reached ~50–100mm³ the animals were randomized into 3 groups and treated with docetaxel (15mg/kg), epirubicin (3mg/kg) and DMSO control. After 4 weeks of treatment all tumors were harvested and genomic DNA was extracted. The gDNA was amplified by PCR (2ug gDNA per 25uL reaction; 5 reactions per harvested tumor) using Q5 PCR master mix and thermal cycles: 98°C for 2min, 22x (98°C 10s, 65°C 30s, 72°C15s), 72°C 2min, 4°C indefinitely. The PCR product was purified using the Ampure cleanup kit and sequenced using the Illumina NextSeq platform.

Analysis of clonal composition using sequenced barcodes

Counts of expressed barcode tags were extracted from FASTQ files using Pycashier (https://github.com/DaylinMorgan/pycashier), an extension of Cashier (developed by Russell Durrett). Pycashier is a wrapper that uses cutadapt to identify and remove flanking adapter sequences and starcode (v1.3) to performs fast minimum levenshtein clustering. FASTQ Files were uncompressed and renamed in order to run Pycashier. Default parameters were used including Pred quality filtering of 30, levenshtein clustering with a ratio 3, distance 1 and filtering reads that have counts less than 10. A cut-off of 10 was confirmed as an appropriate threshold in a cumulative plot of barcode frequencies as described in Bystrykh L et al., On a cumulative frequency scale, the curve noticeably bent at the turning point, indicating the change from high frequencies (true barcodes) to very low frequencies of barcodes (false) close to 10. The experiment was run in 2 batches, the second batch (D3–5, E1–5) were grown for longer (2–3 days) in culture. There was no difference in the median count between batches, but there was a difference in library depth (colSums) and upper quantiles. Quantile normalization (R package preprocessCore::normalize.quantiles) was applied to adjust the batch effect and read counts were adjusted by the median ratio method (Li et al., Genome Biology volume 15, Article number: 554 (2014)). Specifically, the adjusted read count x is calculated as the rounded value of x/s, where s is the size factor in experiment j and computed as the median of all size factors calculated from individual sgRNA read counts. In total 168,852 barcodes were observed in one or more samples. The greatest barcode diversity (n=134,488) was in T0 whereas a mean of 30,889, 22,806 and 30,752 barcodes were detected in vehicle (V), Dox (D) and Epi (E) treated animals respectively. In differential barcode analysis, barcode with 2 or fewer observations in vehicle, or in either treatment animal (min(D,E)) were excluded. Non-parametric Rank Product (geometric mean rank) analysis (Breitling et al., 2004) was used to identify barcodes that were differentially enriched in treatment compared to vehicle. P values were corrected for multiple testing using the False Discovery Rate (FDR, (Benjamini and Hochberg, 1995)). Additionally differential barcode analysis using a negative bionominal model was applied as described by Li el al (Li et al., 2014). Linear regression analysis of barcodes was performed using Bioconductor packages voom, limma and p value were adjusted for multiple testing using the FDR. Wilcox signed rank test was performed on barcode alpha (Shannon) and Beta (Bray Curtis) diversity indices, calculated using the R package functions vegan::vegdist, vegan::betadisper. Shannon index assessed the barcode diversity within an observation, whereas Bray Curtis indices measure the pairwise dissimilarity between observations (samples). The Bray Curtis index is 0 if the two observations share all the same barcodes and 1 they don’t share any common barcodes. Tukey Honest significant differences’ test was also used to examine the overlap in confidence intervals on the differences between the means of Bray Curtis dissimilarity indices. Bioconductor R version 4.0.2 was used in all statistical analyses and code is available on github (http://github.com/aedin/barcodes).

NADPH/NADP measurement assay

The cells were seeded as 2-D cultures in 10 cm dishes (1×10⁶ cells/dish). JQ1 was added the next day, followed by addition of epirubicin or DMSO-control 4h later. The cells were harvested after 12h by scrapping in cold cell recovery solution, followed by centrifugation and pellet storage at −80°C. The cellular pellets from 2-D and 3-D cultures were processed for NADPH/NADP measurement using the NADP/NADPH Quantitation Colorimetric Kit (BioVision), according to manufacturer’s instructions.

Apoptosis assay

The rate of cell apoptosis in 2D and 3D cultures was estimated by luminescence measurement, using the RealTime-Glo™ Annexin V Apoptosis assay (Promega, Madison, WI) following the manufacturer’s instructions. Briefly, the cells were incubated with two annexin V fusion proteins (Annexin V-LgBiT and Annexin V-SmBiT) which contain complementary subunits of NanoBiT® Luciferase; and a time-released luciferase. The cell surface levels of membrane phosphatidylserine in apoptotic cells brings the Annexin V-LgBiT and Annexin V-SmBiT luciferase subunits into complementing proximity, which is reflected by the strength of bioluminescence signal emitted in culture (quadruplicate measurements).

BH3 profiling

Flow cytometry-based BH3 profiling was performed as previously described (Ryan et al., 2016). Cells were incubated for 90 minutes with the peptides BIM and BID and stained with Zombie Aqua Dye (Biolegend, #423101) to assess viability and cytochrome c-Alexa Fluor 488 (Biolegend, #612308). Flow cytometry analysis was performed on a LSRFortessa X-20 machine (BD Biosciences) and the results analyzed on the FlowJo software. The percentage of cytochrome C release in persister organoids vs. control was calculated based on average of the Mean, Geometric Mean, or Median of Fluorescence Intensity (MFI) of the Cytochrome C -Alexa Fluor 488 antibody normalized to the respective values of the positive and negative controls in the same sample (respectively cells exposed to buffer or 25μM alamethicin). For the MYC-KD and JQ1-treated cells the percentage of cytochrome C release (duplicates) was calculated based on the Median of Fluorescence Intensity (MFI) of the Cytochrome C-Alexa Fluor 488 antibody normalized to the respective values of the positive and negative controls in the same sample (respectively cells exposed to buffer or 25μM alamethicin). Priming was visualized in graphs using the average percentage of cytochrome C release ± SEM of duplicates. The experiments were repeated 2 times to ensure consistency of the data.

Cell cycle analyses

The gels were initially disrupted mechanically using cold PBS. The 3-D organoids isolated by centrifugation and incubated with trypsin for 20 min at 37°C to generate single-cell suspensions. A similar incubation time was applied for the 2D cultures. The cells isolated from 2D or 3D cultures were stained using the FxCycle™ PI/RNase Staining Solution (Thermo Fisher) according to manufacturer’s instructions and analyzed by flow cytometry.

Histological analyses

Cells suspended in matrigel were casted in 12-well plates (1.5 ml/well) and the gels were let to solidify for 1h at 37°C. After gel solidification, the 3-D cultures were supplemented with medium and incubated for 7–8 days to allow for organoid formation. Gels were then extracted from the culture vessel and fixed in 10% formalin. Subsequently, the gels were either stained with carmine (whole mount preparation) or processed for paraffin embedding and immunohistochemistry (slide sections), as previously described (Dhimolea et al., 2010). Tumor tissues were fixed in 10% formalin and processed for immunohistochemistry following standard protocols.

Western Blot

Cell lysates were analyzed using 7% NuPage Tris-Acetate gels (Invitrogen) according to manufacturers’ instructions and transferred to a nitrocellulose membrane. The membrane was incubated overnight with the primary antibody [anti-Myc antibody (N-262; Santa Cruz Biotechnologies) 1:200 followed by incubation with the secondary anti-rabbit antibody (7074S; Cell Signaling Technology) 1:2000 at room temperature for 1h. The protein bands were visualized using chemiluminescence solution (RPN2232; GE Healthcare).

CRISPR-Cas9 mediated MYC gene knockout

We cloned two target sgRNA sequences for MYC, as well one target sequence for each of the olfactory receptor genes OR10G2 or OR10G3 (used as controls, given the biological inertness of these genes in epithelial cancers), into the Lenti X CRISPR-Cas9 (Clontech) according to manufacturer’s instructions. The sgRNA target sequences are as follows:

MYC_sgRNA1: ACAACGTCTTGGAGCGCCAG

MYC_sgRNA2: GCCGTATTTCTACTGCGACG

OR10G2: GGAGGCTTCTTAGATTTGGG

OR10G3: GCTTAGCAGTCATGAGCACA

Lentiviral particles were generated as described above. MDAMB-231/Cas9+ cells for 16h in cell medium, and viral prep at 2:1 ratio (polybrene was added at 8μg/ml concentration). After the end of the incubation with the viral preps, cells were washed and incubated for additional 24h. Next, the transduced cells were transferred to 6-well plates and selected with hygromycin (1mg/ml) for 72h. The effect of chemotherapeutics on the transduced cells was assessed by two cell viability assays: a) Transduced cells were plated in 384-well plate and treated the next day with chemotherapeutics at the indicated doses. Cell viability was assessed by CTG as described above. b) The transduced cells were exposed chemotherapeutics in concentrations and durations indicated. Next, the cells were harvested using trypsin and incubated with 0.4% trypan blue for 5 min.; the number of dead and live cells was estimated using the Countess II FL automated cell counter.

CRISPR-interference mediated MYC gene knockdown

We performed the CRISPR-interference experiments similarly to previously-described protocols (Kearns et al., 2014). The plasmid construct pHAGE TRE dCas9-KRAB expressing the tet-regulatable dCas9-KRAB was obtained from Addgene (Plasmid #50917) and used to generate lentiviral particles using Lenti-X 293T cells and packaging plasmids psPAX2 and MD2.G. The plasmid was introduced in MDAMB-231 cells via lentiviral transfections and the cells expressing the tet-inducible dCas9-KRAB vector were selected using G418 (1 mg/mg). The following sense (and respective antisense) sgRNA sequences against MYC or olfactory receptor genes were generated:

MYC_1: CACCgAGGCAGAGGGAGCGAGCGGG

MYC_2: CACCgCCCGGCTCTTCCACCCTAGC

MYC_3: CACCgGCAGCGCAGCTCTGCTCGCC

MYC_4: CACCgGCTGTAGTAATTCCAGCGAG

MYC_5: CACCgGCGCTGCGGGCGTCCTGGGAA

OR10G3: CACCgGTGGATACGGAATTCCTGTC

OR6S1: CACCgCAACAGAGTTCGTCCTGGCA

and subsequently cloned in the by ligation into the BsmB1I site of plasmid pXPR_050. Next, the MDAMB-231 cells expression the doxycycline-inducible dCas9-KRAB vector were lentivirally transfected with the sgRNA-expressing plasmids and selected using puromycin (1ug/ml). Doxycycline-induced dCas9 expression was titrated by western blot. Gene kockdown was also assesed by western blot. Exposure to 2ug/ml docycycline induced Myc downregulation within 3 days (e.g. Figure S6A).

Myc overexpression

Doxycycline-mediated induction of Myc was enabled by first transducing MDAMB-231 or MCF7 cells with the lentiviral vectors pLV[Exp]-EGFP/Neo-EFS>tTS/rtTA, and next with pLV[Exp]-mCherry:T2A:Bsd-TRE>hMyc or pLV[Exp]-mCherry:T2A:Bsd-TRE>mCherry (as control); followed by selection with the respective antibiotics. The lentiviral vectors were purchased from VectorBuilder.

Animal studies

The in vivo experiments were conducted in accordance with the guidelines of the DFCI Institutional Animal Care and Use Committee. MSK-PCa1 and HCI003 patient-derived organoid cells were suspended in medium containing 20% matrigel and injected subcutaneously in the flank of NOD.Cg-Prkdc^scid mice (2×10⁶ cells/injection) implanted with dehydrotestosterone or estradiol pellets, respectively. When the average tumor diameter reached approximately 1cm, the mice were separated into treatment vs. control groups as indicated and treated with either viblastine/control (2mg/kg bw, weekly i.v. injections, for MSK-PCa1 PDX) or afatinib/control (25mg/kg bw by daily oral gavage, for the HCI003 PDX). Tumor growth or response to treatment was monitored by caliper measurements or by in vivo bioluminescence measurements. Residual drug-persistent tumors were collected after 2–4 weeks, at the plateau phase of viability curves.

Proteomic analysis (Reverse phase protein analysis, RPPA)

RPPA analyses were conducted by the Functional Proteomics RPPA Core facility (Houston, TX), which is supported by MD Anderson Cancer Center Support Grant # 5 P30 CA016672–40. We examined the expression levels of 303 proteins in HCI009 and HCI011 3D-PDO samples (isolated as described above) treated with either DMSO-control or docetaxel 100nM (experimental duplicates). Lysates of these protein samples were processed by the core based on protocols and procedures outlined on the website of the facility (https://www.mdanderson.org/research/research-resources/core-facilities/functional-proteomics-rppa-core/rppa-process.html). Antibody array signal reads from replicates were averaged for each protein. The resulting matrix of RPPA data was further processed to generate log-transformed linear normalized data, linear normalized, and normalized median centered data.

Sample collection and gene expression analysis

The cells in 2-D cultures were collected by scrapping in cold cell recovery solution, pelleted by centrifugation and stored at −80°C. The 3-D gels were disrupted mechanically using cold cell recovery solution and 3-D organoids were collected by centrifugation and stored at −80°C. The 3-D gels of cell lines were collected 8–10 days after seeding to allow for organoid formation. The gels of 3-D PDO models HCI002, HCI003, MSK-PCa1, MSK-PCa3 were collected 10–12 days after drug exposure. The gels of 3-D PDO models HCI009 and HCI011 were collected 5 days after drug exposure and 5 days post drug washout (after 12 days total drug exposure). The gels of 3-D PDO MSK-PCa7 were collected 12 days after drug exposure and 5 days after drug washout. The RNA was isolated using Qiagen RNA isolation kit. The RNA isolated from MDAMB-231 cell lines 2-D cultures and 3-D organoids (experimental duplicates) was supplemented with ERCC RNA spike-in mix to adjust for the cell number, and subsequently analyzed by hybridization in Affymetrix Prime View 3’ microarrays at the Molecular Biology Core Facility (Dana-Farber Cancer Institute). We performed differential gene expression analysis for microarray data of each cell line in 2D and 3D organoid models using the bioconductor LIMMA packageA (Ritchie et al., 2015). The RNA isolated from 3D-PDOs and MDAMB-231 organoids (experimental triplicates) treated with vehicle, docetaxel, or afatinib was analyzed by RNA sequencing using Illumina NextSeq 500 Next Gen at the Molecular Biology Core Facility (Dana-Farber Cancer Institute). We performed differential gene expression analysis for the RNA-seq data using the bioconductor DEseq2 package (Love et al., 2014).

Public datasets used

The molecular profiles of bulk post-treatment residual patient tumors were obtained from datasets GSE87455 (Kimbung et al., 2018), GSE32603 (Magbanua et al., 2015), GSE43816 (Gruosso et al., 2016) and GSE32072 (Gonzalez-Angulo et al., 2012). The PROMIX trial dataset GSE87455 (which is the largest of the three clinical datasets) was analyzed at individual patient level using the residual tumor samples collected after 2 cycles of chemotherapy (69 patients; e.g. Figure 3G) or after 6 cycles of chemotherapy, at surgery (58 patients). The analysis of I-SPY trial dataset GSE32603 (39 patients) at patient level is depicted in Figure S5D.

The molecular profiles of single cell post-treatment residual patient tumors were obtained from dataset SRP114962 (Kim et al., 2018). For the single-cell RNAseq profile, sample annotations were imported from the public dataset. Raw fastq files were downloaded from ENA. Adapter and quality trimming was performed using BBMap (v.38.73), and single end reads were aligned against the GENCODE v.32 hg38 transcriptome using salmon (v.1.0.0). In cases where read pairs were present, only read 2 was used. Further processing was performed in R/Bioconductor. Salmon quantification results were imported with tximport, using the transcripts to gene symbol mappings generated from the GENCODE gtf file. Analysis was performed using Seurat after normalization with SCTransform. Signature scores were calculated by the SingleCellSignatureScorer module of the SingleCellSignatureExplorer using counts obtained from the Seurat SCTransform function. For plotting the single-cell copy numbers (Figure 3C) we used the processed data matrix of long integer copy numbers from the whole genome sequencing data provided by the authors of the paper (Kim et al., 2018).

The transcriptional profiles of embryonic stages were obtained from the Boroviak et al.(Boroviak et al., 2015) dataset with EBI ArrayExpress (Kolesnikov et al., 2015) accession E-MTAB-2958. We generated a read count matrix from the provided BAM files using feature counts from the subread package (version 1.5.2) (Liao et al., 2014). For the gene level feature counts, we used the GTF annotation file from cell ranger 10X genomics (version cellranger-mm10–1.2.0; http://www.10xgenomics.com). Mouse ensemble gene identifiers were mapped to MGI Gene symbols and retrieved using the bioconductor biomaRt package (Durinck et al., 2009).

Human to mouse orthologous gene mapping

We retrieved orthology mapping of mouse gene symbols to human gene symbols from the MGI database (Blake et al., 2017). The mapping was extracted from the report HOM_MouseHumanSequence.rpt and is available at ftp://ftp.informatics.jax.org (accessed August 2017). When multiple human genes mapped to a single mouse gene symbol, we assigned the expression value of corresponding mouse gene symbol for each human gene symbol. Mouse gene symbols without a human mapped gene symbol were excluded from the analysis.

Gene set enrichment analysis (GSEA)

The GSEA analysis was performed using the GSEA software (Subramanian et al., 2005) version 3.0 (gsea-3.0.jar) based on the preranked option. For each differential gene expression analysis we use the per gene foldchange for gene ranking. For the analysis of microarray data sets if a gene was represented by multiple probes we use the foldchange of the probe with the smallest nominal p-value. For the visualizations we use the GSEA enrichment scores (ES) and the gene set size normalized enrichment scores (NES). For generating the gene set enrichment plots in Figures 1E, 6G and 7H we used we used the top 1000 genes with FDR<0.05 in persister HCI002 model and the mouse embryonic diapause (Boroviak dataset (Boroviak et al., 2015)).

Embryonic diapause-like (EDL) signature score and survival analysis

The EDL score was estimated by the Spearman correlation coefficient for the GSEA enrichment scores of all terms with nominal significance in Embryonic Diapause (DIA vs E4.5 EPI) compared to the respective GSEA enrichment scores estimated each for each of the organoid, PDX and patient derived (e.g. Treatment vs Baseline, Surgery vs Baseline) comparisons. We performed the gene sets using the Molecular Signature Database (MSigDB v6.2) (Subramanian et al., 2005), the Gene Ontology GO database (Ashburner et al., 2000) and the ConsensusPathDB-human database (CPDB) (Kamburov et al., 2013). The GO gene set collection was extracted from the bioconductor org.Hs.eg.db package (Carlson et al., 2016) for the biological process, cellular component and molecular function domain.

Based on the intensity of EDL signature in residual tumors, patients of the PROMIX or ISPY datasets were classified into EDL-high and EDL-low groups, to compare the transcriptional profiles of baseline tumors between these two groups in each dataset and derive a transcriptional signature of the proneness of baselined tumors eventually develop EDL during treatment. An EDL-proneness signature based on the top 50 upregulated and top 50 genes downregulated in baseline tumors of EDL-high vs. EDL-low groups from the PROMIX or ISPY datasets were then used to stratify BrCa patients of the METABRIC dataset and examine potential differences clinical outome between the patients with high vs. low EDL-proneness signature. Concordant results were obtained with EDL-proneness signatures generated based on different numbers of up- or down-regulated genes in baseline tumors of EDL-high vs. EDL-low groups.

Pairwise 2D GO analysis

We defined the gene sets from the Gene Ontology database(Ashburner et al., 2000) that were extracted from the bioconductor org.Hs.eg.db package and from the comprehensive pathway and gene set collection provided from the ConsensusPathDB-human database (CPDB)(Kamburov et al., 2013). We implemented the 2D annotation enrichment analysis in the script language R as previously described(Cox and Mann, 2012). We performed Hotelling’s T² test for each gene set using the Manova function in R. We considered a total of about 8816 GO-terms and 3110 CPDB-pathways that are associated with at least 10 genes represented in a respective gene expression dataset. For each term we tested the null-hypothesis whether the composite mean value of the gene ranks of a given gene set from the two datasets is larger or smaller than the global composite mean value of the genes that are not members of the respective gene sets for the two datasets. In order to consider multiple hypothesis testing, we applied the FDR multiple testing correction procedure as previously defined(Benjamini and Hochberg, 1995). As described in Cox et al(Cox and Mann, 2012), we employed a position score for each gene set for the two input datasets. The position score is the rank average in each dataset as defined by a value ranging between −1 and 1. For the analysis we ranked the genes for each dataset based on the log₂ fold change values.

In order to estimate a similarity score between two datasets, we estimated the spearman correlation coefficient between the position scores for the significant gene sets with FDR<=0.05. The analyses were conducted both separately for GO-terms and CPDB-pathways, as well as combined. Figure 3H and Figures S4A, S5A–S4C indicate the scores of commonly significant GO-terms and CPDB-pathways between drug-persistent 3-D organoids vs. vehicle-control counterparts (plotted in the Y axis) and embryonic stages vs. E4.5 epiblast (Boroviak et al., 2015) plotted on the X axis. The Spearman correlation coefficient is indicated in each Figure.

Gene essentiality data

CERES scores (Meyers et al., 2017) as metrics of gene essentiality in conventional 2D cultures, were examined for various genes commonly upregulated/downregulated across different drug-persistent PDO and PDX models based on results from genome-scale in vitro CRISPR/Cas9 gene knockout screens of the Dependency Map program (www.depmap.org) performed with the AVANA sgRNA library. The patterns of essentiality for different genes of interest were consistent across several different releases of data from the Dependency Map studies.