Human Organoids Share Structural and Genetic Features with Primary Pancreatic Adenocarcinoma Tumors

This article is featured in Highlights of This Issue, p. 1

In the United States, pancreatic cancer is the third leading cause of cancer-related death, accounting for over 50,000 cases annually, with a survival rate of less than 10% (1). The most common form of pancreatic cancer, pancreatic ductal adenocarcinoma (PDAC), is often diagnosed at a late stage, limiting effective therapeutic interventions. Only 10%–15% of patients are eligible for surgery, the only potentially curative option (2). Although adjuvant chemotherapy has shown incremental improvements in resected patients, the vast majority recurs and succumbs to metastatic disease with a 5-year survival of only 1%–2%. This poor survival underscores the need for novel tools to rapidly and precisely match patients with effective therapies.

The limited number of available preclinical PDAC models has been a major hurdle for discovery and translational research. Development of in vitro cell culture models, as well as transgenic and patient-derived xenograft (PDX) mouse models, enables investigation of biological mechanisms and new treatments. However, traditional two-dimensional (2D) immortalized monolayer cell lines are limited, in that they do not reflect the variability or the structure of PDAC tumors. Primary cells cultured as monolayers do not reflect the heterogeneity of the primary tumor due to selection in culture, lack of stromal–stromal–cell communication, tissue-specific architecture, and mechanical cues. Combined, these impact gene expression, which is key for modeling therapeutic responses (3–5). Genetically engineered mice, such as LSL-Kras^G12D/+; LSL-Trp53^R172H/+; and Pdx-1-Cre are useful in that they are models of carcinogenesis within the pancreas; however, they only represent specific genetic mutations and fail to recapitulate the variability that is characteristic of PDAC (6–8). PDX models have complex architecture and utilize heterogeneous patient tissue, but the clinical applicability is challenging due to prohibitive time requirements (up to 6 months to establish tumor growth), costs associated with in vivo systems, the large sample size required (∼100 mm³), and the influence of infiltrating murine stromal cells on the tumor (9, 10). The latter factor may influence the recent observations that the more times a PDX tumor is passaged through mice, the more transcriptionally “mouse-like” it becomes (11). Because of the challenges related to existing systems, improved models of PDAC are essential.

Recently developed three-dimensional (3D) organoids may provide dramatic benefits compared with 2D cell culture and PDX models for their use in personalized medicine and drug discovery. Organoid models are rapidly and inexpensively developed from primary or metastatic tumor tissue for expansion and molecular profiling. Organoids also recapitulate biological features of a 3D environment, as demonstrated in pancreatic, colorectal, prostate, and mammary cancers (12–15). Importantly, PDAC organoids can be grown with a high rate of success from very small amounts of tissue collected from diagnostic fine-needle biopsies (16, 17), core-needle biopsies, and from excisional intraoperative biopsies (18, 19).

Some barriers to the use of organoids in personalized therapies include lack of in-depth assessment of the histopathology, genetic stability, and molecular profiling of organoid models. Comparisons are needed between organoids and primary tumors of origin, as well as among individual organoids derived from the same patient. Before organoids can be routinely employed as model systems, more characterization of cellular and molecular properties are required. Specifically, morphologic characterization and comparisons between organoids and primary tumors of origin, as well as among individual organoids derived from the same patient, has been insufficiently reported. Therefore, we performed histopathologic profiling on organoid and PDX models for comparison with their corresponding primary tumors. In addition, we performed a thorough genomic characterization of organoid cultures by deep sequencing to determine whether the models represent the genomic constitution of the primary tumor, and if they retain their genetic characteristics over time. Furthermore, we performed ex vivo drug testing to improve our understanding of the degree to which organoids model the therapeutic response of the corresponding tumor of origin. We have defined the histopathology, genetic heterogeneity, and therapeutic sensitivity profiles of organoid models derived from patients with PDAC as an initial effort to provide thorough pathologic and genetic comparison between PDAC organoid models and the tumors from which they are derived.

Between 2014 and 2017, 10 tumor samples were collected from patients with PDAC under IRB12-1108 and IRB13-1149. Clinical information is provided (Supplementary Table S1). Clinical data were obtained from the electronic medical record (Epic). PDAC and adjacent, uninvolved pancreas were obtained from patients undergoing pancreatic resections at The University of Chicago Medicine (UCM, Chicago, IL) facilities. Samples were confirmed to be tumor or benign based on pathologic assessment. Normal pancreas tissue was obtained from patients who underwent resection for benign lesions and demonstrated no features of pancreatitis.

Research was performed under protocols IS00000556 and IS00000424 Institutional Animal Care and approved by Northwestern University (Chicago, IL). Human tumor samples were obtained from UCM pancreatic cancer patients and deidentified. In brief, freshly resected human tumor samples (0.2 g) were sent to Northwestern University (Chicago, IL) and transplanted subcutaneously into non-obese diabetic/severe combined immunodeficient (NOD-SCID) gamma (NSG) mice (The Jackson Laboratory). Tumor volumes [(length × width²)/2] were measured weekly. Default endpoints for any animal were loss of 15% body weight compared with pretumor weight, or other signs of distress. When PDX tumor reached 1.5 cm diameter, the mouse was euthanized and freshly resected 2 mm² tumor pieces were retransplanted to NSG mice for in vivo studies. A piece of subcutaneous PDX tumor was fixed in 10% formalin and processed to paraffin-embedding. Sections (5 μm) were stained with hematoxylin and eosin (H&E) for histopathologic evaluation.

For initial treatment studies, mice bearing subcutaneous pancreatic PDX tumors from patient 1 were staged to approximately 200 mm³ prior to initiation of treatments and randomized to two treatment groups: control (vehicle) and FOLFIRINOX (4 mg/kg oxaliplatin, Zydus Hospira, 61703-363-22; 50 mg/kg leucovorin, Teva, 0703514501; 25 mg/kg irinotecan, Pfizer, 0009752903; 50 mg/kg 5-FU, Thomas Scientific, F6627; n = 3, 2 tumors per mouse). Vehicle or drugs were injected intraperitoneally once a week for 3 weeks.

For the next set of studies, mice bearing subcutaneous pancreatic PDX tumors from patient 1 were staged to approximately 150 mm³ prior to initiation of treatments and randomized to 3 treatment groups: control (vehicle), gemcitabine (100 mg/kg, Sigma, G6423) and gemcitabine (100 mg/kg) + abraxane (30 mg/kg, Celgene, 6881713450); n = 3, 2 tumors per mouse. Vehicle or drugs were injected intraperitoneally once a week for 4 weeks.

To establish 2D cells and PDX-derived organoid cultures, a piece of subcutaneous PDX tumors was extracted, minced, and digested with collagenase type XI (0.125 mg/mL, Sigma, C9407) and dispase (0.125 mg/mL Gibco, 17105041), in DMEM (Gibco, 11995065) and incubated from 0.5 to 1 hour at 37°C with gentle shaking. Cells were spun and dissociated with TrypLE Express (Thermo Fisher Scientific, 12605-010) and 10 μg/mL DNase I (DN25, Sigma) for 10 minutes at 37°C, and washed with DMEM. For culturing 2D cells, an aliquot (1/10) of the dissociated suspension was seeded into a petri dish with complete 2D media (DMEM-F12 Advanced, GIBCO, 12634-010; HEPES buffer, Invitrogen, 15630-080; penicillin/streptomycin, Thermo Fisher Scientific, 15140-122; l-GlutaMax, Invitrogen, 35050-061; FBS 5% Life Technologies, 13028-014; EGF 10 ng/mL, GIBCO, PMG8043; bovine pituitary extract ∼ 57.6 μg/mL, hydrocortisone 2 μg/mL, Sigma, H0888 and insulin human recombinant 0.56 μg/mL, Life Technologies, 12585-014). For culturing organoids, dissociated cells were washed and embedded in growth factor–reduced (GFR) Matrigel (Corning, 356231) and cultured in complete media [Intesticult (Stemcell Technologies, 6005), A83-01 (0.5 μmol/L, Sigma, SML0788), fibroblast growth factor 10 (FGF10, 100 ng/mL, Gibco, PHG024)], gastrin I (10 nmol/L, Sigma, 17105-041), N-acetyl-l-cysteine (10 mmol/L, Sigma, A9165), nicotinamide (10 mmol/L, Sigma, N0636), B27 supplement (1×, Gibco, 17504-044), primocin (1 mg/mL, InvivoGen, ant-pm-1), and Y-27632 (10.5 μmol/L Tocris, 1254). To establish tumor-derived 2D cell lines and organoids, resected primary tumor samples from pancreatic resection surgeries were utilized as above. Organoids were passaged via mechanical dissociation with TrypLE Express (Thermo Fisher Scientific, 12605-010) and passage was performed weekly with a 1:2 ratio.

A gastrointestinal pathologist scored architecture of tumor and organoids from H&E-stained slides. The patterns were i. Simple (score = 1): tumor epithelium composed of a single layer of epithelial cells, but often with some loss of nuclear polarity, ii. Papillary (score = 2): tumor cell composed of focal stratification into multiple layers with occasional luminal projections, iii. Solid/cribriform (score = 3): tumor cells grow in a syncytium, within occasional small gland structures. An average aggregate score was calculated from the total number of fields and the percent of fields containing each of the three patterns, according to the following equation: Histology score = [Fields of pattern 1 + Fields of pattern 2 + Fields of pattern 3]/total fields.

Human samples, PDX and organoids were fixed with 10% formalin, paraffin embedded and sectioned (5 μm). The following antibodies were used: CK7 (Agilent/DAKO, M7018, mouse monoclonal, clone OV-TL 12/30), CK19 (Agilent/DALO, M0888, mouse monoclonal, clone RCK108), CK20 (abcam, ab76126, rabbit IgG), CEA (abcam, ab4451, mouse monoclonal, clone 26/3/13), Claudin-4 (abcam, ab53156, rabbit IgG.), CA19-9 (Thermo Fisher, MA5-13275, mouse monoclonal, clone 121 SLE) and P53 (Calbiochem, OP-43-100UG, mouse monoclonal, clone DO-1). Slides were scanned (ScanScope XT, 20x). Quantification was performed using Imagescope (Aperio eSlide Manager Sofware, Leica). Staining intensity was scored 0 to 3 using the positive pixel count algorithm.

A GI pathologist scored H&E sections of the tumor and circled areas of high tumor cellularity. DNA was extracted from the same area of an adjacent normal section using the QIAamp DNA Micro kit (56304), quantified and processed into indexed libraries using the Illumina TruSeq kit (FC-121-2001).

For 2D cell lines, DNA was extracted from approximately 3 million cells using the DNeasy Blood and Tissue kit (Qiagen, 69504) following manufacturer's protocol. DNA was extracted from a single organoid for analysis by DNA-seq. The Matrigel in one well of a 24-well tissue culture plate was depolymerized with 4°C PBS. Under 4X brightfield magnification, a single organoid was isolated in 3–4 uL of liquid Matrigel, and DNA was extracted and amplified as recommended by the manufacturer. The GenomePlex Single Cell Whole Genome Amplification Kit (Sigma-Aldrich, WGA4-10RXN) and the Repli-G Single Cell kit (Qiagen, 150343) were utilized. Following amplification, the DNA Clean and Concentrator (Zymo Research, D4013) kit was utilized and libraries were generated using KAPA LTP Library Preparation kit (KAPA Biosystems, KR0453).

The indexed libraries were pooled, and targeted capture of the exons of 1,213 known cancer genes was carried out (20). In brief, this approach utilizes a tiered assay system in which highly clinically relevant genes (tier 1, n = 316) are sequenced approximately 2.18-fold deeper than the remaining (tier 2) genes. The enriched libraries were sequenced (101 bp paired-end reads) using Illumina.

FastQ files were trimmed using Cutadapt v1.9.1 (21) to remove low quality bases and adapter sequences. The reads were aligned against human genome reference (hg19) using BWA-MEM v0.7.8 (http://bio-bwa.sourceforge.net). For PDX samples, reads arising from mouse stroma were filtered out by mapping against a custom-built genome of the mouse strain used for PDX generation, and removing all perfectly aligned reads, prior to human reference alignment. PCR duplicates were removed while sorting on-the-fly using novosort v1.03.9 (Novacraft Technologies Sdn Bhd, http://www.novocraft.com/products/novosort/). Bedtools v2.22.1 (22) was used to ascertain coverage at tier 1/tier 2 loci. A threshold of 300X mean coverage at tier1 genes was used for tumor, normal, PDX and 2D cell line samples, while organoid samples were allowed more leniency in tier 1 coverage (median coverage 226x).

Single nucleotide variants (SNVs) were called using MuTect v1.1.7 (23) while insertions and deletions (indels) were called using scalpel-discovery 0.5.3 (24). For both SNVs and indels, only calls within genomic regions targeted by the capture panel were retained for subsequent analyses. All variants were annotated using ENSEMBL's Variant Effect Predictor version 84 (25) and only variants within coding regions or disrupted splice sites were included in analyses. Calls with a variant allele frequency (VAF) < 5% for the patient tumor and normal samples, position coverage < 30, or an allele frequency > = 0.01 in ExAC (26) were removed. Exceptions were made in cases where either a low frequency variant was called in the models and present at less than 5% VAF in the tumor, and when a variant was called in the tumor but missed the cutoff in the models. To further improve the quality of indel calls, putative calls in low complexity genomic regions identified using Dustmasker (27) prone to false positives were removed.

A GI pathologist circled areas of high tumor cellularity on H&E sections of the tumor and PDX. RNA was extracted from the same area of an adjacent section using RNeasy FFPE kit (Qiagen, 73504) after treating with Deparaffinization Solution (Qiagen, 19093).

The RNeasy Mini kit (Qiagen, 74104) was used for extracting RNA from both 2D cell lines and organoids. Approximately 3 million cells and 1 million cells were used as starting material for 2D and organoids, respectively.

Exome capture was carried out using the RNA access kit (Illumina, RS-301-2001) for the majority of samples, with the exception of 2D cell lines from patients 1–3 which had total RNAseq performed with ribosome depletion (Illumina, MRZH11124). Libraries were pooled and sequenced on HiSeq or NextSeq 500 instruments, to generate paired-end 100 bp reads.

All samples were quality checked and aligned using the following general outline: FastQ files were adapter and quality-trimmed using cutadapt v 1.15 (21) and read quality was evaluated further using FastQC (28). Reads were aligned using STAR v 2.4.2a (29). Alignment quality was assessed using Picard v2.16 CollectRnaSeqMetrics and quantification was performed using eXpress (30). A minimum of 20 million reads per sample was required for quantification and subsequent differential expression (DE) analysis.

The R princomp and DESeq2 (31) libraries were used for DE analysis using the effective counts output from eXpress. PCA combined with linear regression was used to determine which covariates should be included in the final regression model for differential expression. In the set of samples used for multi-patient comparison, no significant covariates were identified besides the independent variable of interest. However, for the in-depth analysis of patient 1, there was a detectable batch effect influencing principal component 1 (adjusted r² 0.758, p-value 0.00028), and was included as a covariate in the general linear model using DESeq2. Genes were considered differentially expressed if they had an adjusted p-value < 0.05 and a log2 fold change greater than 1 or less than −1. A representative subset of the top DE genes by adjusted p-value from each biotype was obtained and the union (518 unique genes) was used in downstream analysis. Data from normal patient tissues in the study were combined and used as controls for all DE analyses.

The single cell RNA-seq dataset is comprised of 15 separate drop-seq runs. For each run, a total of 2000 mature organoids were pooled, typsinized with TrypLE in presence of 10 μmol/L Y-27632, strained into a single cell suspension and counted. The final concentration of the cell suspension varied between 100–120 cells/uL depending on the flow rate at which monodisperse droplets were formed. An aliquot from the cell suspension was stained with Trypan Blue to ensure high cell viability (>90%) at the time of establishing the cell preparation protocol. Organoids were harvested for drop-seq analysis at same growth size to ensure uniformity from run to run.

Nanoliter-sized droplets containing single cells and barcoded beads were prepared as described in Macosko and colleagues (32). Briefly, droplets generated by co-flowing 1 mL each of cell and bead solutions through a microfluidic device were collected. The droplets were broken and the beads were collected and reverse-transcribed. The cDNA obtained was then PCR amplified and quantified. Finally, the cDNA was fragmented and amplified using primers that allow amplification of only the 3′ ends, processed into RNA-seq libraries and sequenced on Illumina NextSeq 500. Sequencing data was analyzed using the pipeline created by Macosko and colleagues (32) to generate digital gene expression (DGE) matrices. Only cells with more than 250 genes were retained for subsequent clustering.

We developed a novel clustering approach based on a class of probabilistic generative models called topic models. In brief, the DGE is first normalized, scaled and log transformed. The contribution of undesirable sources of variation such as cell cycle phase, batch effect and mitochondrial gene expression is assessed by calculating scores for each of these factors and regressinged out of the data (33). Next, a training set comprised of cells that express at least 900 genes is created and highly variable genes are identified (34). A topic model is then generated from the training set using the highly variable genes as the vocabulary (35). Then, a topic distribution is calculated for the cells not included in the training set by comparing their gene expression to that of the training set. Statistical significance of topics is calculated based on KL-divergence and insignificant topics are eliminated. Finally, cells are assigned to the topic with the highest probability, and clusters with fewer than 50 cells are removed. Pairwise marker genes for each cluster are identified using genes that are expressed in least 10% of the cells in both clusters and have average log fold change >1 between clusters. Clusters that do not have at least 10 unique marker genes are merged (36).

Relative viability of organoids was determined by measuring Annexin V fluorescence in real-time with Incucyte (Essen Bioscience). Organoids were seeded (3,000 cells/well) into a U-bottom ULA 96 well-plate (Costar, 7007). Cells were incubated with Annexin V Green Reagent (Essen Bioscience, 4642) and treated with gemcitabine hydrochloride (G6423, Sigma) at 3 nmol/L, 10 nmol/L, and 30 nmol/L. Real-time fluorescence was measured with Incucyte during 72 hours.

Relative number of viable cancer cells was determined by measuring the optical density using CellTiter 96 Aqueous One Solution Cell Proliferation Assay kit [3-(4,5-dimethylthiazol-2-yl)-5-(3-carboxymethoxyphenyl)-2-(4-sulfophenyl)-2H-tetrazolium, inner salt; MTS] according to the manufacturer's instructions (Promega, G3580). To resemble clinical pharmacokinetics of SOC drugs, tumor organoids were treated with vehicle, gemcitabine, gemcitabine and paclitaxel, and FOLFIRINOX (oxaliplatin, leucovorin, irinotecan, 5-fluorouracil) for 4 hours at indicated concentrations. Post-treatment, test compounds were replaced with fresh compound-free complete culture media. Organoids were grown another 72 hours and relative cell growth and death (as compared with optical density signal before the start of the treatment) was measured by MTS cell viability assay after 72 hours of growth.

For differentially expressed genes, R PCA package princomp was used to check for batch effects. In the set of samples used for multipatient comparison, batch effects were negligible and therefore were not regressed out. However, for the in-depth analysis of patient 1, there was a detectable batch effect influencing principal component 1 (adjusted r² = 0.758), so it was included in the batch model for DESeq2. Differential expression with P values and P_adj values were calculated using default methods of negative binomial global linear model and Wald statistics. For the multipatient comparison, as sequencing was done over two lanes each for each patient, each lane was treated as a replicate to increase statistical significance; for the in-depth comparison, for all but PDX there were multiple runs (2–4). In all cases, normal tissues from all patients in the study were used as the control for differential expression. Differentially expressed genes were analyzed by Ingenuity Pathway Analysis to predict cancer related pathways. The Benjamini–Hochberg method was utilized and an FDR cutoff of q < 0.05 was applied. Two-way ANOVA Multiple comparisons were performed in all in vivo/in vitro studies.

The genomic data from this publication has been deposited at www.ncbi.nlm.nih.gov. The submission ID is SUB4032998. The BioProject ID is PRJNA471134. Private login credentials are pending with the NIH (Bethesda, MD).

To generate clinically relevant PDAC models, benign and tumor tissue were collected at the time of surgical resection (Supplementary Table S1). PDAC organoids were grown from both PDX and primary tissue (15). PDAC tumors from 5 patients were used to establish PDX tumors, followed by the growth of PDX-derived organoids from those xenografts. In addition, PDAC tumors from another 5 patients were utilized to grow primary tumor–derived organoids directly (without PDX; Supplementary Fig. S1).

To validate tumor pathology in PDX and organoid models, H&E histologic analysis of uninvolved pancreas, primary tumor, PDX, and organoid models was performed (Fig. 1A and B). A gastrointestinal pathologist compared the architecture and cell morphology in both PDX and organoid models to the primary tumors (Fig. 1A). Benign ducts revealed simple columnar epithelium without significant cytologic atypia, and areas of PDAC epithelium showed various degrees of architectural disorder and cytologic atypia. Specifically, tumors from patients 3 and 5 (Fig. 1A) showed well-differentiated tumor glands with small basally oriented nuclei and little or no cellular stratification. Tumor glands from patient 4 were slightly more atypical and demonstrated regions of stratification. Glands from patients 1 and 2 showed the greatest extent of cell atypia and cell stratification. We divided the architectural patterns of the organoids into three morphologic classifications: simple, papillary, and solid/cribriform (Fig. 1C). These architectures were quantified, and a high degree of correlation was observed between the morphologic structure of the primary tumors and the corresponding organoids (R² = 0.64; Fig. 1D, left). The architecture of the PDX models was more complex and atypical compared with the organoid models, in part, due to the presence of mouse stroma within the tumor. Nevertheless, organoids derived from the PDX models recapitulated the microscopic features of the primary tumors.

To determine whether organoids derived directly from the primary tumor also maintain tissue architecture similar to PDX-derived organoids, we grew organoids directly from the primary tumor (patients 6–10) and performed histologic analysis of tissue derived from these 5 patients. We observed a similar glandular architecture between primary tumor and organoid models (Fig. 1B). Like those derived from PDX-passaged tumors, organoids derived directly from primary tumors maintained the morphologic structure of their corresponding primary tumor. Specifically, the well-differentiated tumor glands seen in the organoids from patients 8 and 10 resembled the relatively simple glandular architecture observed in the primary tumors from these patients. Likewise, for patients 6, 7, and 9, the organoids demonstrated more cellular stratification akin to the primary tumors of these same patients. Quantitative analysis of glandular architecture demonstrated a high degree of morphologic correlation between primary PDAC and organoids (Fig. 1D, right; R² = 0.95).

To define the relevant protein expression profiles of the PDX and organoid models and to compare these protein expression patterns with both benign pancreatic and tumor tissue, we used a panel of IHC markers. Because the majority of benign tissue obtained from the tumor resections demonstrated evidence of chronic pancreatitis, we also stained normal pancreatic samples selected from patients without PDAC and with no evidence of pancreatitis. Cytokeratin markers 7 and 20 (CK7 and CK20) were selected to assess for pancreatic epithelial differentiation, and cellular tumor antigen p53 (p53), Claudin-4 and Cancer Antigen 19-9 (CA19-9) were selected due to their frequent expression in neoplastic pancreatic ductal epithelium (refs. 37–40; Fig. 2A and B).

To quantify the extent that protein expression in the organoid and PDX models was recapitulated by the staining observed in the primary tumor, pixel count and intensity were measured for all cases. Cytokeratin proteins, normally found in the intracytoplasmic cytoskeleton of epithelial tissue, change their protein profile expression in PDAC. CK7 and CK19 are expressed in 90%–100% of PDAC (41, 42). In contrast, CK20 is found in less than 20% of PDAC cases (43). We found that the protein expression profiles of these cytokeratins are maintained in the PDX and organoid models (Fig. 2A–C; Supplementary Fig. S2). We also examined nuclear labeling of the mutated tumor suppressor p53, as aberrant expression occurs in at least 70% of PDAC (39, 44). Consistent with these previous observations, we found that normal pancreatic epithelium showed low levels of p53 expression while p53 is increased in primary tumors and corresponding PDX and organoid models (Fig. 2A–C). We also stained for Claudin-4, which is often overexpressed in PDAC (45, 46), and observed that its pattern of overexpression is elevated in both PDX and organoids models (Fig. 2A–C). Next, we examined tumor epithelial cells for the expression of the sialyl Lewis motif that corresponds to CA19-9 (47). We found that all of the primary tumors studied expressed CA19-9, and that both the PDX and organoid models stained positive for CA19-9 as well. In normal pancreatic epithelium, we saw that CA19-9 was restricted to the luminal surface. However, in tumor cells, the antigen was localized in the cytoplasm, basolateral membrane, and occasionally at the apical cell membrane. In addition, the staining intensity was increased in tumor cells compared with normal epithelium (ref. 48; Fig. 2A–C). Finally, we investigated the expression of CEA, which has been reported to have positive staining in PDAC (49). We observed positive staining in the tumor group, as well as the PDX and organoid models (Supplementary Fig. S2). In summary, we observed a strong concordance of expression between the primary tumor, organoid, and PDX models for all the protein biomarkers examined.

We next investigated whether PDX and organoid models preserve the genomic and transcriptomic characteristics of their corresponding tumors. We extracted DNA from the primary tumor, peripheral blood, PDX, organoids, and 2D cell lines for 7 of the patients in our study and performed targeted sequencing of a 1,213 genes panel comprised of known cancer-related genes (20). We achieved an average coverage of 930X for tier 1 genes and 491X for tier 2 genes (20). Single organoids (5–9) were sequenced from each patient to assess the extent of intra-patient model heterogeneity. Consistent with previous reports (50–52), primary tumors from all patients had mutations in KRAS and TP53, which were reproduced in PDX, organoid, and 2D cell lines. Furthermore, mutations specific to individual patients such as CDKN2A, NALCN, ZBTB16, and PARP1 were also present in corresponding tumor models (Fig. 3A). Remarkably, these same patterns were observed in primary tumor–derived organoids (Fig. 3B). Genetic mutations were stable over time and passage number in PDX-derived organoids, as evidenced by sequencing (patient 2) at passages 4 and 10 (approximately 8 weeks apart; Fig. 3C). Interestingly, the variant allele frequency for mutations in the tumor models was higher than that of the primary tumor (median value: 57.69 and 12.44, respectively).

To understand differences in transcriptional profiles between PDAC from individual patients and their corresponding PDX and organoid models, we performed RNA sequencing (RNA-seq) on the primary tumor, adjacent normal tissue, PDX models, and 2D cell lines from patients 1–4. Principal Component Analysis (PCA) of the RNA-seq data indicated that the primary tumor and related models (PDX and 2D cell lines) are more similar in their gene expression profiles than to adjacent normal tissue (Supplementary Fig. S3). Results showed a high correlation between primary tumor and PDX (R²: 0.91–0.96) and 2D cell lines (R²: 0.71–0.96), in contrast to primary tumor and matched normal tissue (R²: 0.09–0.63; Fig. 4A). To verify that the normal pancreatic tissue collected is representative of normal pancreas, we compared transcriptional profiles of normal pancreas to brain, muscle, and pancreatic tissue from the Genotype-Tissue Expression project (GTEx). These analyses verified that the normal pancreas samples collected in this study were indeed representative of larger normal pancreatic tissue cohorts (Supplementary Fig. S4).

To further investigate RNA expression profiles, pathway analysis was performed on the tumor, PDX, and 2D cell line RNA-seq data. We discovered that cancer-related gene expression pathways, including p53 signaling and cell-cycle regulation are similarly up- or downregulated across all sample types (Supplementary Fig. S5A). Notably, upregulated p53 signaling is consistent with an increase in observed TP53 protein levels (Fig. 2C). Differences in immune signaling dependent on model type were apparent (Supplementary Fig. S5B), which is not surprising because PDXs are grown in immunodeficient mice and the 2D cell lines lack any immune component.

Finally, we compared organoid RNA-seq analysis from organoids (patient 1) to the models and primary tumor. We observed a strong correlation between organoids and primary tumor (R²: 0.66), and with both PDX (R²: 0.84) and 2D cell lines (R²: 0.85; Supplementary Fig. S6). The Venn diagram representation comparing tumor, 2D cell line, organoids, and PDX shows that the majority of the differentially expressed genes (1,815 genes) overlap between the three models and the tumor sample (Fig. 4B). Furthermore, a majority of differentially expressed genes are similarly up- or down-regulated across all models in comparison with the primary tumor (Fig. 4C). Functional analysis of differentially expressed genes indicated that proliferation and cancer-related genes are upregulated across primary tumor and models when compared with normal, whereas some models show decreased expression of invasion-related genes (Fig. 4D). Overall, we observed that 2D and 3D models maintain much of the gene expression that contributes to tumor phenotypes, although differences are apparent.

To understand the cellular composition of the organoid models, we performed high-throughput single-cell RNA-seq using the Drop-seq platform (32). We sequenced the transcriptomes of 7,675 single cells derived from organoids (patient 1). We detected a median of 417 genes per cell (mean 715 genes per cell) at an average sequencing depth of approximately 20,000 reads per cell (Fig. 5A).

Drop-seq measurements were compared with the bulk RNA-seq data discussed above, and the data were correlated (R²: 0.65; Fig. 5B). The single-cell transcriptional results were clustered to identify distinct subgroups of cells, using a computational approach based on a class of probabilistic generative topic models. Clustering of the Drop-seq data resulted in three transcriptionally distinct groups of cells (Fig. 5C). The largest cluster comprised 97.2% of the total population of cells, while the remaining 2.8% of the cells were divided between cluster 2 (0.9%) and cluster 3 (1.9%). Interestingly, cluster 3 expressed Prominin-1 (CD133), a marker for pancreatic cancer stem cells (Table 1; ref. 53). Cluster 2 expressed markers amphiregulin (AREG) and epiregulin (EREG), which are EGF ligands that have been implicated in tumorigenesis and poor outcomes (Table 1; refs. 54, 55). A list of all marker genes for each cluster is provided in the Supplementary Materials (Supplementary List S1). These data indicate that the pancreatic cancer–derived organoids are largely clonal populations derived from the primary tumor.

To determine whether PDX- and primary tumor–derived organoid models can be used for drug response assays, we tested organoids for dosage-dependent and drug-specific responses and investigated whether organoid drug responses could recapitulate PDX results.

Organoids were treated with gemcitabine, a chemotherapeutic agent used in PDAC, and resulting levels of apoptosis were measured. PDX-derived organoids (patient 1) were treated with three different concentrations of gemcitabine (3 nmol/L, 10 nmol/L, 30 nmol/L) for 72 hours and apoptosis was measured in real time (Supplementary Movie S1). Treatments caused a dose-dependent apoptotic response to gemcitabine (P = 0.003 for 10 nmol/L, P = 0.009 for 30 nmol/L; Fig. 6A). To determine whether primary tumor–derived organoid samples would yield similar results, the assay was performed on the sample from patient 6. Again, we observed a dose-dependent response, and interestingly, there was a full response to gemcitabine at 10 nmol/L indicating that the sensitivity threshold was higher in this sample compared with patient 1 (Fig. 6B).

To investigate the in vivo drug response to FOLFIRINOX, we utilized a PDX from a patient that became resistant to FOLFIRINOX postsurgery (patient 1, Supplementary Table S1). Following PDX engraftment, tumor volume was measured over 4 weeks. We did not find a significant benefit of FOLFORINOX treatment (Fig. 6C). To mimic patient treatments after FOLFIRINOX resistance, we performed a second PDX experiment employing the PDX tumors that become resistant to FOLFORINOX in vivo. Mice were treated with a standard-of-care combination therapy for PDAC, gemcitabine with or without abraxane (protein-bound paclitaxel; Fig. 6D). Tumor volume was measured over 5 weeks and showed that combination therapy reduced the growth by approximately 85% compared with controls (gemcitabine vs. controls, P = 0.0017), and reduced tumor growth almost twice as effectively as gemcitabine alone (gemcitabine/abraxane vs. gemcitabine, P = 0.0025). To determine whether organoid models recapitulated this response, PDX-derived organoids (patient 1) were treated with the same agents followed by analysis of cell viability. A combination of gemcitabine and paclitaxel dramatically decreased cell viability, whereas gemcitabine alone only reduced approximately 20% cell growth, which was not significantly different from untreated organoids (Fig. 6E). Overall, these data demonstrated that treatment of PDX-derived organoids recapitulates the effects of in vivo treatments. Therefore, organoids derived from either PDX or primary tumors are not only suitable for drug treatment studies, but results may also reflect patient-specific sensitivities to drugs.

Here we have defined the histopathology, genetic heterogeneity, and therapeutic sensitivity profiles of organoid models derived from patients with PDAC. The results indicate strong concordance at the morphologic and molecular levels, and are promising for the development of personalized organoid models that could be used to help guide therapeutic decisions. We have shown that organoids can recapitulate the morphologic architecture of the primary tumor of origin. Specifically, we identified three major epithelial growth patterns that are found in the primary tumor and are well reflected in organoid models. These structures include simple, papillary, and cribriform morphologies that can be indicators of degree of tumor differentiation. The use of protein expression profiles in histopathologic classification is used along with structural information to determine tumor type and degree of differentiation. We demonstrated that commonly used protein markers for PDAC analysis are maintained between both organoid and PDX models and primary tumors. Importantly, because organoids can be grown from a fine-needle biopsy, information such as structure and protein expression patterns may be obtainable even prior to surgical resection, opening the possibility to facilitate rapid interventions based on personalized organoid molecular, cellular, and therapeutic response characterizations.

Genetic sequencing of PDAC is challenging due to the large proportion of stromal cells in the tumor and heterogeneity of epithelial cancer cells. Because organoids have only an epithelial component, they have an amplified signal compared with primary tumor, allowing deeper sequencing of tumor cells and increased detection of mutations. As observed previously (12–15), we have shown genetic concordance between PDAC organoids with the primary tumor. Most mutations were replicated in the organoid models and the primary tumor, but we observed some genetic variation that could be due to the heterogeneity of the primary tumor and the outgrowth of different clones. Prior studies have not investigated the heterogeneity of the organoid population or the stability of genetic characteristics over time. This is the first report of single organoid DNA sequencing over time using multiple passages of organoids from the same patient tumor. Our results indicate that the genetic composition of the organoids was consistent over time, for at least ten passages.

Understanding the heterogeneity of cellular states within individual patient derived organoids is a key step toward determining how they might be used as personalized models. The similarity of transcriptomes of individual cells within organoids from patient 1 are overall very similar and indicative of clonal origins, with >97% of cells occupying a single large cluster. However, we detected two small populations of cells that may offer insights into the biology of the disease. For instance, the group of cells distinguished by CD133 expression may be a stem-like subpopulation. Similarly, the group of cells distinguished by AREG expression may also be biologically relevant. Multiple studies have highlighted the role of AREG in tumorigenesis of epithelial malignancies, including pancreatic cancer. AREG is often overexpressed in PDAC, is predictive of poor outcome and encourages growth of pancreatic cancer cell lines (56). Therefore, the sequencing of individual cells from organoids provides valuable information about clonal populations.

Caveats for the organoid system include limitations of immune signaling and inability to test interactions with immune cells and response to immunotherapies. The development of 3D coculture systems with an immune component is an important future step. Another weakness of the organoid system is the lack of stromal cells that exist in patient tumors, limiting our understanding of how the stroma impacts drug response. The genetic context of the stroma has been shown to be important in clinical outcome. Similar to immune cell coculture, stromal cell coculture is an important next step. Nonetheless, the observation of patient-specific gemcitabine responses in organoid models, as well as the differential response of both PDX and organoids to common combination therapies such as gemcitabine plus abraxane or FOLFIRINOX, offer great promise for the use of PDAC organoids to help determine patient-specific therapies. Increasing the numbers of mutational and germline genetic backgrounds represented by patient organoids, and developing methods to study immune and stroma interactions with organoids, will undoubtedly be important to fulfill the full potential of organoid models. In addition, little is known about the differences that may be present between primary tumor-derived organoids and organoids derived from metastatic sites. Similar investigations to this study are required to determine to what extent organoids grown from metastatic PDAC tumors recapitulate the histologic, proteomic, and genetic characteristics from the primary tumor and how they have evolved.

In summary, our study has demonstrated that organoids are potentially valuable for drug treatment studies because they maintain distinct phenotypes and respond differently to drug combinations and dosage. These models may allow for ex vivo drug testing of patient samples to steer treatment decisions. Furthermore, these models have potential utility for high-throughput drug discovery assays.

K.P. White is a president and has ownership interest (including stock, patents, etc.) at Tempus Labs. No potential conflicts of interest were disclosed by the other authors.

Conception and design: I. Romero-Calvo, C.R. Weber, A.P. Mazar, M.M. Buschmann, K. Roggin

Development of methodology: I. Romero-Calvo, C.R. Weber, M. Brown, A. Ugolkov, W. Qiang, J.P. Segal, A. Rzhetsky, A.P. Mazar, M.M. Buschmann

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): I. Romero-Calvo, C.R. Weber, M. Ray, K. Kirby, S.M. Sparrow, A. Ugolkov, W. Qiang, Y. Zhang, J.P. Segal, A.P. Mazar, M.M. Buschmann

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): I. Romero-Calvo, C.R. Weber, M. Ray, M. Brown, R.K. Nandi, T.M. Long, A. Ugolkov, W. Qiang, Y. Zhang, T. Brunetti, J.P. Segal, A.P. Mazar, M.M. Buschmann, R. Weichselbaum

Writing, review, and/or revision of the manuscript: I. Romero-Calvo, C.R. Weber, M. Ray, M. Brown, K. Kirby, T.M. Long, S.M. Sparrow, A. Ugolkov, H. Kindler, A.P. Mazar, M.M. Buschmann, R. Weichselbaum, K. Roggin

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): I. Romero-Calvo, C.R. Weber, M. Brown, K. Kirby, S.M. Sparrow, M.M. Buschmann

Study supervision: I. Romero-Calvo, M.M. Buschmann, R. Weichselbaum, K. Roggin

We would like to thank Drs. Mitchell Posner and Jeffrey Matthews for allowing us to consent patients in their clinics and for providing study guidance; Dr. William Dale for providing direction in human subjects protocol development and clinical research operations; Margaret Eber, Megan Flanagan, and Teresa Barry for consenting patients and coordinating biospecimens; Marlin Amy Halder for assistance with DNA sequencing studies; Shuang Qin Zhang for consulting on therapeutic treatment of PDX models; Bradley Long and Sabah Kadri for assistance with DNA sequencing and data interpretation; and to the Friends of Jack Karp and Barbara Turf for their generous financial support of this work. I. Romero-Calvo, C.R. Weber, M. Ray, M.M. Buschmann, K. Kirby, S.M. Sparrow, T. Brunetti, and M.M. Buschmann's work, as well as all reagents and clinical research expenses, was provided by the Friends of Jack Karp and Barbara Turf. H. Kindler is partially supported by NIH U10CA180836. The work of A. Ugolkov and W. Qiang was partially supported by NCI U54CA193419 and CCSG P30 CA060553.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.