Abstract
Pancreatic ductal adenocarcinoma (PDAC) is characterized by extensive desmoplasia, which challenges the molecular analyses of bulk tumor samples. Here we FACS-purified epithelial cells from human PDAC and normal pancreas and derived their genome-wide transcriptome and DNA methylome landscapes. Clustering based on DNA methylation revealed two distinct PDAC groups displaying different methylation patterns at regions encoding repeat elements. Methylationlow tumors are characterized by higher expression of endogenous retroviral transcripts and double-stranded RNA sensors, which lead to a cell-intrinsic activation of an interferon signature (IFNsign). This results in a protumorigenic microenvironment and poor patient outcome. Methylationlow/IFNsignhigh and Methylationhigh/IFNsignlow PDAC cells preserve lineage traits, respective of normal ductal or acinar pancreatic cells. Moreover, ductal-derived KrasG12D/Trp53−/− mouse PDACs show higher expression of IFNsign compared with acinar-derived counterparts. Collectively, our data point to two different origins and etiologies of human PDACs, with the aggressive Methylationlow/IFNsignhigh subtype potentially targetable by agents blocking intrinsic IFN signaling.
The mutational landscapes of PDAC alone cannot explain the observed interpatient heterogeneity. We identified two PDAC subtypes characterized by differential DNA methylation, preserving traits from normal ductal/acinar cells associated with IFN signaling. Our work suggests that epigenetic traits and the cell of origin contribute to PDAC heterogeneity.
This article is highlighted in the In This Issue feature, p. 521
Introduction
With an overall five-year survival rate of about 9%, pancreatic ductal adenocarcinoma (PDAC) still belongs to the cancer entities with the poorest prognosis (Cancer Research UK and American Cancer Society; ref. 1). A hallmark of PDAC is the immense contribution of desmoplastic stroma to the total tumor mass. Typically, epithelial PDAC cells constitute only a minority of cells (as low as 10%) within the tumor mass. The extensive stromal presence represents a major technical challenge when aiming to interrogate the specific molecular signals contributed by the neoplastic epithelial cell compartment from data generated using bulk tumor samples. Some studies have addressed this question in the field of gene expression using sophisticated bioinformatics, labor-intensive laser microdissection (2, 3), and, more recently, single-cell sequencing approaches (4–7). PDAC methylation studies rely, so far, mostly on bulk tumor analyses (8) or comparisons with cell lines (9) and xenografts (10–12), and have focused on predefined regions of the genome. Here, we report whole-genome bisulfite sequencing combined with transcriptome and functional analyses of FACS-purified human PDAC and normal pancreas epithelial cells. We provide insights on how the altered DNA methylome landscape of PDAC cells affects tumor aggressiveness and suggest the existence of different cells of origin for two distinct PDAC subsets.
Results
Transcriptome and Methylome of Purified Human PDAC and Normal Pancreatic Epithelial Cells
We aimed to determine the genome-wide DNA methylation and transcriptomic landscapes of purified epithelial cancer cells present within human PDAC tumor masses. To avoid confounding effects of stromal cells in PDAC tumors, we established a protocol to FACS-isolate epithelial cells (as EpCAM+/CD45−) from primary untreated PDAC tumors (PDAC-Ep; n = 7). For comparison, we FACS-isolated epithelial cells (as EpCAM+/CD45−, mainly composed of ductal and acinar cells) from tumor-free adjacent (normal) pancreatic tissues (NP-Ep; n = 6; Fig. 1A; Supplementary Table S1; Extended Data Fig. 1). Using low-input technologies, we derived the transcriptomes [RNA sequencing (RNA-seq); ref. 13] and the DNA methylomes [tagmentation-based whole-genome bisulfite sequencing (TBWGBS); ref. 14] of these epithelial cells (Fig. 1A; Supplementary Tables S2 and S3).
As expected, isolated epithelial cells showed high expression of EPCAM and very low or undetectable expression of stromal markers (CD45, PDGFRB, and VWF; Supplementary Fig. S1A). Compared with NP-Ep, PDAC-Ep cells showed robust expression of KRT19 and decreased expression of SMAD4, as expected for PDAC cells (refs. 15, 16; Supplementary Fig. S1B). 48.7% of the RNA-sequenced reads from PDAC-Ep cells harbored KRAS mutations, whereas no KRAS mutations were found in NP-Ep RNA reads (Supplementary Fig. S1C). Altogether, these data demonstrate the purity of the isolated epithelial cell populations. Additionally, although we cannot exclude that in some PDACs EpCAM− epithelial tumor cells coexist, only 0% to 8% mutant KRAS reads were detected in RNA-seq data from paired EpCAM−/CD45− samples (which comprise fibroblastic cells, vascular endothelial cells, nerve tissue as well as potentially epithelial cells that have lost EpCAM expression), indicating that our sorting strategy captured the majority of the PDAC epithelial cells (Supplementary Fig. S1D).
The methylome analysis identified 56,177 regions showing differential DNA methylation [differentially methylated regions (DMR)] between PDAC-Ep and NP-Ep (cutoff 40% differences in methylation; Supplementary Fig. S1E; Supplementary Table S4). These DMRs covered 1.78 × 107 bp of the genome. Of these, a larger proportion was hypomethylated in PDAC-Ep compared with NP-Ep cells than hypermethylated (66% versus 34%, respectively; Supplementary Fig. S1E and S1F). This reflects the different global distribution of DNA methylation levels of single CpGs (CG) along the genomes of these samples (Supplementary Fig. S1G). The majority of DMRs (79%) were localized outside CpG islands (CGI; 10%) or CGI shores (11%; Supplementary Fig. S1H). The vast majority of CGIs showed a gain of methylation in PDAC-Ep, whereas the regions outside of CGIs showed a marked loss of methylation (Supplementary Fig. S1I). As these non-CGI sites make up the majority in the genome, a global hypomethylation pattern is observed in PDAC-Ep.
Next, we correlated the methylomes and transcriptomes by closest gene analysis (13, 17–19). Only 13% of the DMRs (7,397) correlated with upregulation or downregulation of gene expression (FDR < 0.05; Supplementary Fig. S2A–S2C). Enrichment analysis of the DMRs correlating to gene-expression changes revealed deregulated pathways including axon guidance, WNT, TGFβ, NOTCH signaling, and apoptosis, confirming previous reports (refs. 8–10; Supplementary Table S5). DMRs correlating to differential gene expression were also mainly hypomethylated in PDAC-Ep compared with NP-Ep (69%; Supplementary Fig. S2B). The 7,397 DMRs correlating with changes in gene expression could be classified into four groups according to the direction of the changes (Supplementary Fig. S2B). Inverse correlation between DNA methylation and gene expression (hyper/loss and hypo/gain) comprised the majority of DMRs (63.5%). These DMRs were located closer to transcription start sites (TSS) and had a larger overlap with promoter regions (Supplementary Fig. S2B and S2C, top plots). Hypermethylated regions showed a bigger overlap to CGIs and enhancers compared with hypomethylated regions (Supplementary Fig. S2B and S2C, bottom plots). The hypermethylation/loss of expression group included genes of the pancreatic progenitor lineage (NR5A2, PDX1, and PROX1; Supplementary Fig. S2D and S2E and data not shown). Pancreatic progenitor genes were also included in the hypomethylation/loss of expression group (e.g., ONECUT1), which also comprised genes involved in tight junctions (e.g., CDH5). Gain of expression correlating to loss of methylation (hypo/gain) occurred in genes previously reported to be upregulated in PDAC, such as SERPINB5 (20), genes involved in cytoskeletal organization (ITGB1 and ACTB), as well as genes related to inflammatory response (CCL20, CCL22, and IL1A; Supplementary Fig. S2F–S2H and data not shown). Lastly, hypermethylation/gain of expression occurred in gene body regions distant to TSS enriched in CGI (Supplementary Fig. S2B and S2C), and it was found in genes related to embryonic development (HOXA2, HOXA3, and HOXB9) and matrix metalloproteases (MMP) including MMP9 (Supplementary Fig. S2I and S2J). In summary, we provide a comprehensive data set comprising whole-genome methylome and transcriptome data of freshly isolated pure human PDAC and normal pancreatic epithelial cells. Our data confirm, but also massively expand, the number of known DMRs and associated transcriptional changes in purified PDAC epithelial cells to the whole-genome level (Supplementary Table S4; Supplementary Fig. S2K; refs. 8–10, 21).
Methylation Profiling Separates Two PDAC Clusters Associated with IFN Signaling and Survival
The large amount of DMRs identified between PDAC-Ep and NP-Ep cells illustrates the profound changes that occur at the level of DNA methylation in epithelial cells during PDAC carcinogenesis. Interestingly, by performing principal component analysis (PCA) using PDAC-Ep versus NP-Ep DMRs, the tumor samples separated into two clearly distinguishable groups: methylation cluster 1 (MC1) and methylation cluster 2 (MC2; Fig. 1B). A separation of the PDAC-Ep samples was also revealed by unsupervised PCA (Supplementary Fig. S3A). The existence of two PDAC methylation groups was further validated by consensus clustering (Fig. 1C and Extended Data Fig. 2). Comparison of the transcriptomes of MC1 and MC2 PDAC cells identified 320 differentially expressed genes (FDR < 0.01; Supplementary Fig. S3B; Supplementary Table S6). No changes in expression of EPCAM or KRT19 were observed, indicating that the presence of the two groups was not due to contamination with other cell types (Supplementary Fig. S3C). Gene set enrichment analysis (GSEA; refs. 22, 23) using the hallmark gene sets revealed a marked enrichment for IFNα and IFNγ response signatures in MC2 samples (Fig. 1D; Supplementary Fig. S3D). Among other genes, the transcription factor STAT1, a master regulator of IFN signaling, as well as IFI44L and IFITM1 were highly upregulated in MC2 compared with MC1 tumor cells (Fig. 1E and F; Supplementary Fig. S3E). Interestingly, several immune-checkpoint modulators were also higher expressed in MC2 as compared with MC1 (Supplementary Fig. S3F). To further analyze the status of IFN signaling in public PDAC data sets, we built a 47-gene signature (IFNsign) by combining genes enriched in MC2 from both IFNα and IFNg response gene sets (see Methods for details; Supplementary Fig. S3G). As expected, the IFNsign was robustly overexpressed in MC2 compared with MC1 (Fig. 1G; Supplementary Fig. S3H). Analysis of expression data sets from human PDAC tumors revealed a strong correlation between the expression of different individual IFN-related genes and our IFNsign (Fig. 1H; Supplementary Fig. S4A). These data demonstrate that an entire IFN program is highly expressed in a subset of patients with PDAC. Importantly, patients with an IFNsignhigh or STAT1high status showed a significantly overall worse survival compared with patients with low expression of this IFN program or STAT1 (Fig. 1I; Supplementary Fig. S4B). Of note, there were no significant differences on tumor epithelial content between IFNsign high and low patients (Supplementary Fig. S4C). Increased immune infiltration was associated with IFNsignhigh tumors, but this was not linked to changes in survival (Supplementary Fig. S4D and S4E). Altogether, these data suggest that the differences observed in survival are linked to the differential activation of an IFN program in PDAC epithelial cells.
Previous studies, mostly performed by bulk analysis, have investigated the possibility to stratify PDACs by differential gene-expression signatures (2, 12, 24–26). Although some tumors might have a blended correlation to subtypes (27, 28), associations between the expression of the progenitor-like/classical subtype signature with more favorable patient outcome, and the squamous-/basal-/quasi-mesenchymal-like subtype signature with poor patient outcome, have been intensely investigated (2, 12, 24–27, 29, 30). MC1 samples showed enrichment for progenitor-like/classical signatures, whereas MC2 samples were enriched for squamous-/basal-like gene signatures (Supplementary Fig. S4F). Accordingly, STAT1 and IFNsign expression were higher in tumors classified as squamous-/basal-like compared with those classified as progenitor-like/classical (Supplementary Fig. S4G and S4H). Together, these results demonstrate an association of an activated IFN program specifically within the PDAC epithelial cell compartment and the previously described squamous-/basal-like subtype. Nevertheless, STAT1 or IFNsign expression levels were able to further stratify patients with tumors classified as progenitor-like or classical subtype in a more favorable or poor outcome group (Supplementary Fig. S4I and S4J), suggesting both overlapping and complementary stratification power of the IFN signature.
In summary, whole-genome methylome and transcriptome analysis of purified tumor cells revealed two groups of PDACs with distinct methylation landscapes connected to a differential activation of an IFN program and worse overall survival.
IFNsignhigh Status Correlates with Low DNA Methylation at Non-CGI Sites and Upregulation of Repeat-Derived Transcripts
Next, we explored whether the observed differences in the expression of IFN-related genes in PDAC epithelial cells could be explained by differential methylation of the loci encoding IFN target genes themselves. Toward this aim, we first determined DMRs between MC2/IFNsignhigh and MC1/IFNsignlow PDAC-Ep (intertumor DMRs; Supplementary Table S7). We observed no association between changes in methylation and changes in expression of IFN-related genes. Moreover, from the 3,070 DMRs identified between MC1 and MC2, only 114 (3.7%) were located close to differentially expressed genes, and the correlation did not reach significance. These data suggest that the major differences in DNA methylation between tumor groups (>3,000 DMRs between MC1 and MC2) do not directly influence gene expression. Accordingly, DMRs between MC1 and MC2 were located at larger distances to TSS compared with DMRs between normal and PDAC epithelial cells (Supplementary Fig. S5A).
We next performed a genome-wide distribution analysis of DNA methylation, which revealed an overall lower methylation level of MC2 PDAC cells compared with MC1 ones (Fig. 2A; Supplementary Fig. S5B and S5C). Global hypomethylation is a known epigenetic alteration occurring in some cancers (31, 32), and it is typically associated with large genomic regions that coincide with lamina-associated, late-replicating regions termed partially methylated domains (PMD; ref. 33). Indeed, PMD estimation in our cohort revealed that around 47% of the genome of MC2 PDAC cells was occupied by PMDs, in contrast to only 15% in MC1 samples (Supplementary Fig. S5C and S5D). Although MC1-PMDs seemed to have a rather random genomic distribution, as shown by the poor overlap of the PMDs of each individual MC1 sample, MC2-PMDs were highly conserved as shown by a strong overlap between the individual MC2-PMDs (Supplementary Fig. S5E and S5F). These data highlight that hypomethylation in PDAC epithelial cells does not equally occur in all patients.
Because we noticed that hypomethylation in MC2 did also occur outside PMDs (Supplementary Fig. S5G), we next explored methylation levels by genomic location. Methylation at CGI sites showed no major difference between MC2 and MC1 (Fig. 2B). In fact, although clustering based on CGI methylation did separate normal from PDAC samples, this was not sufficient to discriminate tumor groups (Supplementary Fig. S5H; Extended Data Fig. 2). In striking contrast, clustering using regions outside CGIs effectively discriminated MC1 and MC2 tumor groups (Supplementary Fig. S5I–S5K; Extended Data Fig. 2). More specifically, MC2 and MC1 differed remarkably in their methylation levels at sites harboring repetitive elements including long interspersed nuclear elements (LINE), short interspersed nuclear elements (SINE), and long terminal repeats/endogenous retroviruses (LTR/ERV), all of which were consistently lower methylated in MC2/IFNsignhigh PDAC cells (Fig. 2B). These results highlight the value of using whole-genome analyses to include areas beyond regions with high CGI content.
Global demethylation imposed by pharmacologic inhibition of DNA methyltransferases or histone deacetylases has been reported to induce transcriptional derepression of repeats. Especially LTR/ERV transcripts can form double-stranded RNA (dsRNA) structures that activate intracellular pattern recognition receptors leading to increased IFN signaling (34–37). Analysis of the transcriptome of MC1 and MC2 samples revealed that transcripts derived from LTR/ERV and LINE but not SINE sequences showed a higher enrichment in MC2/IFNsignhigh compared with MC1/IFNsignlow PDAC cells (Fig. 2C; Supplementary Fig. S6A). Moreover, MC2 samples were enriched in gene signatures related to dsRNA response (Fig. 2D; Supplementary Fig. S6B). The level of the IFNsign correlated with the expression of dsRNA recognition receptors and with the expression of a subset of ERVs (Fig. 2E and F; Supplementary Fig. S6C and S6D). In summary, based on methylation of purified epithelial cells, we have identified a subset of PDAC tumors showing lower methylation at repeat sites, increased expression of LTR/ERV- and LINE-derived transcripts, an engaged dsRNA response and high IFN signaling.
Validation of MC1 and MC2 Subgroups in an Independent Cohort
To validate our data in an independent cohort, we analyzed samples from additional 11 patients with PDAC (Supplementary Table S1). From these, we determined the transcriptomes of their purified EpCAM+ epithelial cells as for the initial seven samples described above. In addition, PDAC samples were orthotopically transplanted into NSG mice and the PDAC methylomes of the 1° PDXs were determined using Infinium MethylationEPIC arrays (validation set; Supplementary Fig. S7A), which, compared with previous array platforms, additionally cover mainly non-CGI sites (38). 1° PDX tumors derived from the initial MC1/IFNsignlow and MC2/IFNsignhigh samples were also profiled as controls (discovery set; Supplementary Fig. S7A). As previously observed using WGBS data, consensus clustering of the discovery PDX samples confirmed the separation between MC1 and MC2 tumors based on methylation at non-CGIs but not at CGIs (Supplementary Fig. S7B and S7C). Strikingly, three of the new tumors clustered with the three MC2/IFNsignhigh PDX samples, whereas eight clustered with the four MC1/IFNsignlow samples (Supplementary Fig. S7B and S7C). As observed with TBWGBS, the PDX MC2 PDACs also showed (i) lower methylation at LINE, SINE, and LTR/ERV sites (Supplementary Fig. S7D) and (ii) enrichment in expression of IFN signatures and increased STAT1 and ERV expression compared with MC1 tumors (Supplementary Fig. S7E–S7I). Association of the MC1 and MC2 methylation clusters to the previously described classical/progenitor-like and basal-like/squamous subtypes existed, although it was less notable than the association with the STAT1 status (Supplementary Fig. S7B).
Collectively, these data validate the MC1 and MC2 PDAC groups and provide further evidence for the link between hypomethylation-mediated derepression of repeat elements and the activation of the IFN pathway in a total of 18 human PDACs (12 MC1/IFNsignlow and 6 MC2/IFNsignhigh).
Cell-Autonomous Activation of IFN Signaling in Primary PDAC Cells
To expand our findings, we investigated a panel of PDAC patient-derived primary cell lines (PDC; ref. 39) by expression and MethylationEPIC array analyses. We could identify IFNsignhigh and IFNsignlow PDCs according to RNA and protein expression (Fig. 3A; Supplementary Fig. S8A). IFNsignhigh cells were enriched in IFNα and IFNγ response gene sets (Fig. 3B) as well as in previously reported squamous-like/basal/quasi-mesenchymal-like subtype signatures (Supplementary Fig. S8B). In agreement with the poor overall survival associated with IFNsignhigh patients (Fig. 1I), mice injected orthotopically with IFNsignhigh PDCs succumbed earlier (i.e., reached endpoint criteria; Fig. 3C), formed higher-grade tumors (Supplementary Fig. S8C), and displayed a higher incidence of lung metastasis formation compared with mice injected with IFNsignlow PDCs (Fig. 3D; Supplementary Fig. S8D and S8E). These results provide functional in vivo evidence that IFNsignhigh PDAC cells are more aggressive compared with IFNsignlow tumor cells.
DNA methylation analyses showed that IFNsignhigh PDCs had lower genome-wide methylation compared with IFNsignlow PDCs with the strongest differences at sites containing LTR/ERV and LINE repeats (Fig. 3E; Supplementary Fig. S8F). Treatment of IFNsignlow PDCs with the DNA demethylating agent decitabine (5-Aza-2′-Deoxycytidine; 5-Aza-CdR) was sufficient to induce the expression of IFN-related genes (Supplementary Fig. S8G and S8H). Importantly, IFNsign expression in PDCs also correlated with the expression of dsRNA sensors and signatures of response to dsRNA (Fig. 3F; Supplementary Fig. S8I). Indeed, higher levels of cytoplasmic dsRNA were detected in IFNsignhigh PDCs (Fig. 3G and H). Finally, knockdown of MAVS, a master regulator of dsRNA response, using two differentially efficient shRNAs, caused a proportional downregulation of IFN-response genes (Fig. 3I). Together with the data obtained from the 1° PDXs and primary tumors, these data show that IFNsignhigh status in PDAC epithelial cells correlates with lower DNA methylation and higher expression of repetitive elements/dsRNA and indicate that the IFNhigh status in MC2-type PDACs is sustained in a tumor cell–autonomous manner. Interestingly, type I and II IFNs were below detection limit in concentrated conditioned medium (CM) of PDCs in both IFNsignhigh and IFNsignlow PDCs. Blocking of type III IFN receptor (IFNLR1) did not reduce the expression of IFN-related genes in PDCs (Supplementary Fig. S9A and S9B) and neither affected the growth of PDCs in cocultures with stellate cells (StC; see below; Supplementary Fig. S9C). These data exclude that secreted IFNs themselves are responsible for the robust activation of the IFN pathway in PDCs and are consistent with a cell-intrinsic and ERV/dsRNA-mediated mechanism of IFN activation.
Inhibition of IFN Signaling Impairs PDAC Tumor Growth
We next interrogated if the disruption of IFN signaling would affect tumor growth in vivo. Knockdown of STAT1 reduced the expression of IFN-related genes only in IFNsignhigh, but not in IFNsignlow, PDCs (Supplementary Fig. S10A). Strikingly, STAT1 downregulation impaired orthotopic tumor growth of IFNsignhigh but not of IFNsignlow cells in mice (Fig. 3J and K). Similarly, downregulation of MAVS negatively affected the in vivo growth of IFNsignhigh (Fig. 3I and L) but not of IFNsignlow tumors (Fig. 3L; Supplementary Fig. S10B). These results suggest that patients with PDAC bearing IFNsignhigh tumors could benefit from inhibition of the STAT1 pathway. To test this hypothesis, we established subcutaneous tumor models using both PDC types. When tumors were palpable, mice were treated with the FDA-approved JAK–STAT inhibitor ruxolitinib. Ruxolitinib treatment reduced the expression of IFN-related genes in vivo (Supplementary Fig. S10C) and significantly impaired tumor growth of IFNsignhigh cells, whereas it did not interfere with the growth of IFNsignlow PDCs (Fig. 3M). These data establish a preclinical proof of concept that tailored inhibition of IFN signaling might be beneficial for patients with MC2-type PDAC.
IFNsignhigh PDAC Epithelial Cells Reprogram the Stroma toward Tumor Promotion
We next sought a possible explanation contributing to the observed higher aggressiveness of MC2 tumors and the poorer survival of these patients. Cancer-associated fibroblasts (CAF) are typically the most abundant stromal cell population in PDAC (40). The existence of different types of CAFs with diverse or even opposing roles has been recognized in different types of cancer (41), including, more recently, PDAC (42–45). The existence of two main subsets of PDAC-CAFs, the inflammatory (iCAF) and the myofibroblastic (myCAF), has been shown in mouse models and human samples (5, 46, 47, 48). We aimed to investigate if IFNsignhigh/low PDAC epithelial cells could differentially educate or reprogram the fibroblastic stroma, which in turn may have an impact on tumor aggressiveness. Toward this end, we treated normal human pancreatic StCs with tumor cell CM derived from IFNsignhigh or IFNsignlow PDCs, or CM from StCs as control, and analyzed changes in gene expression (Fig. 4A). Indeed, transcriptomes of StCs exposed to IFNsignhigh CM clustered apart from the ones treated with IFNsignlow CM or control medium (Fig. 4A). The top enriched signatures in IFNsignhigh CM-treated StCs were IFNa and IFNγ response programs that included classic IFN-related genes (such as MX1 or IFI44L; Fig. 4B; Supplementary Fig. S11A). Along with classic IFN-related genes, CM of IFNsignhigh tumor cells also induced the expression of inflammatory programs (Supplementary Fig. S11B) and proinflammatory cytokines (e.g., CCL2, IL8, IL1A, and IL1B), some of them previously shown to be expressed by iCAFs (Supplementary Fig. S11A and S11C). In line with these data, patients with MC2 PDAC showed a trend toward higher systemic inflammation as measured by C-reactive protein levels in serum (Supplementary Fig. S11D).
To functionally test the impact of this reprogramming of StCs mediated by IFNsignhigh PDAC cells, we first assessed the in vitro growth of tumor cells alone (monoculture) versus in indirect coculture with StCs. IFNsignhigh PDCs increased the expression of both IFN and proinflammatory cytokines in cocultured StCs as previously shown in the CM experiments (Supplementary Fig. S11E). IFNsignhigh cells grown in the presence of StC expanded more than IFNsignlow cells (Fig. 4C; Supplementary Fig. S11F), suggesting a stronger positive feed-forward loop provided by IFNsignhigh-activated StCs. Treatment with ruxolitinib partially reduced this growth-promoting effect in the coculture setting, whereas it had no effect on PDCs or StCs alone (Supplementary Fig. S11G and S11H). Finally, to test the functional consequences of this reprogramming in vivo, we coinjected IFNsignlow PDC cells together with StCs that had been pretreated with either IFNsignhigh, IFNsignlow, or StC control CM into recipient mice and monitored PDAC tumor growth (Fig. 4D; Supplementary Fig. S11I and S11J). In line with our in vitro results, StCs reprogrammed by IFNsignhigh CM supported PDC-mediated tumor growth to a greater extent compared with StCs reprogrammed by IFNsignlow CM or control CM (Fig. 4D; Supplementary Fig. S11J).
In summary, our data demonstrate that IFNsignhigh cells elicit a distinct activation program in normal StCs toward a proinflammatory and tumor-promoting phenotype, which includes the upregulation of several molecules previously assigned to protumorigenic microenvironments (Supplementary Fig. S11K; refs. 48, 49). These data provide an explanation for the higher aggressiveness and overall poorer survival of patients with MC2 PDAC tumors.
Methylation Traits and Transcriptomic Data Suggest Two Different Types of Cells of PDAC Origin
Despite a markedly different DNA methylome between tumor and normal cells, specific methylation signatures retained from normal cells may exist in tumors, a characteristic previously used to trace back the potential cell of cancer origin (50–55). The healthy adult human exocrine pancreas is composed of two main epithelial cell types, with acinar cells comprising the majority of the pancreatic tissue (70%–60%) and ductal cells representing 20% to 30% (56). Both cell types arise from a common progenitor during pancreas development (57, 58). During transformation, acinar cells can dedifferentiate and/or transdifferentiate to progressively acquire ductal-cell characteristics [so-called acinar to ductal metaplasia (ADM); refs. 57, 59, 60]. Although genetic mouse models have shown that both acinar (through ADM) and ductal cells can give rise to PDAC lesions (61–66), the cell of origin for human PDAC remains highly controversial (63, 67, 68).
We investigated if the DNA methylome of MC1/IFNsignlow and MC2/IFNsignhigh tumors contained remnant footprints of healthy ductal or acinar cells. For this, and since our initial FACS sorted EpCAM+ NP-Ep cells (Fig. 1) are composed of an unknown representation of ductal and acinar cells, we compiled a total of 20 new samples derived from 10 healthy human pancreas donors. This collection comprised (i) primary acinar and ductal cells isolated by FACS (ref. 69; n = 4 samples); (ii) FACS-purified ductal (n = 3 samples) and dedifferentiated acinar cells (n = 3 samples) after in vitro culture of an exocrine mixed population; (iii) the original mixed exocrine population of these samples (refs. 59, 70; n = 3 samples); and (iv) normal bulk pancreas tissue (which is composed of a mix of ductal and acinar cells with possibly few endocrine cells; n = 3 samples; Supplementary Fig. 12A). We extracted the DNA from these samples, performed Infinium MethylationEPIC arrays, and supplemented these data with a previously published methylation data set of human acinar and ductal cells (n = 4; ref. 71). Integration of the data revealed that acinar and ductal cells have very distinct methylation patterns and that these epigenetic traits are maintained during the early stages of acinar dedifferentiation (Supplementary Fig. S12B and S12C). Normal bulk pancreas showed an intermediate profile with slightly higher similarity to acinar/dedifferentiated acinar cells, consistent with the observation that normal pancreas consists of a majority of acinar cells (Supplementary Fig. S12B and S12C; ref. 56).
Next, we joined this new methylome data set of normal pancreas cell types with our collection of 18 PDAC tumors (Supplementary Fig. S7A). Unsupervised clustering of the complete EPIC array data set separated MC2/IFNsignhigh from MC1/IFNsignlow tumors in PC3, as previously observed with the WGBS data (Fig. 5A; Fig. 1B). At the global level, the methylomes of PDAC-Ep cells largely diverged from those of normal epithelial cells regardless of the tumor group (MC1 or MC2) or of the normal epithelial type (acinar or ductal). This was seen in PC1, which explained the largest proportion of the variance between samples (38%; Fig. 5A). This large tumor to normal epithelial cell methylome divergence was further corroborated when using the top differentially methylated sites between acinar and ductal or MC1 and MC2 samples (Supplementary Fig. S13A and S13B, respectively). Although PC1 could clearly be associated with the differences between normal and PDAC-Ep cells, both PC2 and PC3 separated MC1 from MC2 as well as the normal acinar from ductal cells, suggesting that, although globally different, both types of PDAC-Ep cells share some similarities with both acinar and ductal cells at the level of DNA methylation (Fig. 5A; Supplementary Fig. S13C and S13D). Indeed, when checking the methylation status of mature acinar and ductal markers in the tumor samples, all tumors showed a profile similar to that of acinar cells as they showed increased methylation of ductal genes, and a profile similar to ductal cells as they showed increased methylation of acinar genes (Supplementary Fig. S13E). These data indicate that, during transformation, tumor cells massively rearrange their DNA methylome, resulting in a pattern that strongly diverges from the methylation profiles of mature normal acinar or ductal cells. Conserved lineage-related traits associated with cell of origin might therefore be obscured in the tumors by the profound changes happening at the global level.
To further decipher this, we reasoned that methylation marks linked to a ductal cell of origin and maintained during transformation should also be present in human pancreatic ductal epithelial (HPDE) cells, an immortal cell line derived from normal pancreatic ductal cells (72, 73). The methylome of HPDE cells was indeed similar to that of normal ductal cells at CGs associated with classic acinar and ductal markers (Supplementary Fig. S13F). Strikingly, HPDE cells cluster with normal ductal cells when using the top 1,000 CGs driving PC3 but not with those of PC2, suggesting that the information on ductal identity is contained in PC3 (Fig. 5B; Supplementary Fig. S13G). Methylation at these CGs also associated MC2/IFNsignhigh tumors to ductal and HPDE cells, whereas MC1/IFNsignlow tumors intermingled with normal acinar and dedifferentiated acinar cells (Fig. 5B). Additionally, annotation of these CGs to chromatin states identified in purified normal human acinar and ductal cells (69) revealed that PC3-CGs were significantly more associated with acinar, and especially ductal, regulatory regions (i.e., enhancers) than CGs driving other components of the PCA (Supplementary Fig. S13H). Collectively, these data strengthen the hypothesis that PC3 can distinguish between acinar and ductal identity independent of transformation and uncover commonalities at the level of DNA methylation between MC2/IFNsignhigh PDAC epithelial cells and healthy as well as immortalized cells of ductal origin. Conversely, MC1/IFNsignlow PDAC cells show higher similarities to healthy acinar cells in PC3.
To further explore a possible relationship between ductal cells and MC2/IFNsignhigh tumors, we interrogated the status of IFN signaling in healthy human acinar and ductal cells at the transcriptomic level by analyzing publicly available gene-expression data of acinar and ductal cells at the bulk and single-cell level (Fig. 5C; refs. 74–77). GSEA revealed an enrichment of IFN-related signatures in ductal compared with acinar cells in all data sets (Fig. 5D; Supplementary Fig. S14A–S14C). Accordingly, expression of the IFNsign and individual IFN-related genes was consistently higher in healthy ductal compared with acinar cells (Fig. 5E; Supplementary Fig. S14D). Ductal cells also showed enrichment for dsRNA response signatures (Fig. 5F; Supplementary Fig. S14B and S14C) and indeed showed higher expression of repeats, including LTR/ERVs (Fig. 5G; Supplementary Fig. S14E). Interestingly, normal ductal cells also expressed higher levels of the gene signature associated with basal-like PDAC subtype (Supplementary Fig. S14F), which is also enriched in MC2/IFNsignhigh tumor epithelial cells (Supplementary Figs. S4F and S7B). Genes upregulated in MC2/IFNsignhigh tumors but not part of the IFNsign were also enriched in normal ductal cells (Supplementary Fig. S14G). Collectively, these data indicate a broader similarity of MC2/IFNsignhigh tumors to ductal cells at the transcriptomic level.
Taken together, transcriptomic and DNA methylome data link MC2/IFNsignhigh tumors to ductal cells, whereas MC1/IFNsignlow tumors are more similar to acinar cells. These data support the possibility that MC1-like tumors may preferentially arise from acinar cells, whereas MC2-like PDACs originate from ductal cells in the exocrine pancreas.
To explore this possibility further, we analyzed the expression of IFN-response genes in a panel of mouse PDAC cell lines derived from two genetically engineered mouse models both driven by KrasG12D and Trp53 loss, but targeted to either ductal or acinar cells (ref. 63; Fig. 5H; Supplementary Fig. S15A). Strikingly, PDAC cells from ductal-derived tumors (Sox9CreER; KrasLSL-G12D; Trp53f/f mice) expressed significantly higher levels of IFN-related genes and dsRNA sensors compared with the ones derived from tumors of acinar origin (Ptf1aCreER; KrasLSL-G12D; Trp53f/f mice; Fig. 5H; Supplementary Fig. S15A).
Collectively, our data derived from human normal and cancerous pancreatic epithelial cells, as well as from mouse transgenic models, point to a ductal origin for MC2 Methylationlow/IFNsignhigh PDACs, whereas MC1 Methylationhigh/IFNsignlow tumors appear rather derived from acinar cells (Fig. 5I).
Discussion
Based on whole-genome DNA methylome analysis of purified epithelial tumor cells isolated from fresh human PDAC tissue, and analyses of PDAC-derived primary PDX and PDC models, our data uncovered an unexpected link between hypomethylation at repeat elements, the activation of an IFN program, the generation of a proinflammatory tumor-promoting microenvironment, poor prognosis, and a ductal cell of origin in a subset of human PDAC tumors (Fig. 5I). Expression of IFN-related genes (78, 79) or repeats (80) has previously been noted in PDAC, supporting our observations here. However, these descriptive observations had not yet been linked to differential DNA methylation, cell of origin, proinflammatory microenvironment, or PDAC aggressiveness. Activation of IFN signaling through intracellular recognition of derepressed endogenous repeats has been described upon treatment with epigenetic drugs such as decitabine or HDAC inhibitors, a phenomenon termed viral mimicry (34, 35, 81, 82). This viral mimicry is typically linked to growth arrest mainly due to the cytotoxic effects exerted by IFN signaling (83). The extent and makeup of the repeat expression and IFN activation acutely induced by epigenetic drugs may, however, differ from the chronic situation we identified here in untreated naïve MC2-type PDACs. Interestingly, p53 function has been shown to suppress mobilization and expression of retrotransposons, whereas loss of p53 can overcome the IFN-mediated apoptosis induced by deregulated repeat expression in mouse embryonic fibroblasts (37, 84). Indeed, TP53 mutations have a higher incidence in squamous/basal-like PDAC tumors (26, 29, 30), which associate with the MC2 type. This phenomenon seems to be also related to a ductal origin of PDAC, as only PDACs derived from KrasG12D/Trp53−/− ductal cells, but not acinar cells, show high expression of the IFN network (Fig. 5G; Supplementary Fig. S15). Additional mechanisms such as histone modifications or the joint action of certain transcription factors might control repeat expression and mobilization beyond DNA methylation. Additionally, this control might differ between normal and tumor cells, as we observed that ductal cells have higher methylation levels than acinar cells despite bearing higher expression of repeats. Further studies in the field of genome regulation will shed light on the molecular details of these complex mechanisms.
IFNsignhigh/MC2 tumor cells activate and reprogram StCs much more robustly compared with IFNsignlow/MC1 tumor cells, generating a proinflammatory microenvironment that ultimately promotes tumor growth. This process initiated by hypomethylated repeat expression links to the known important roles of the PDAC stroma for tumor cell survival and growth (46, 48, 49, 85). The IFN-mediated stromal activation in MC2 PDACs may be conceptually related to the situation in colorectal cancers, where epithelial cells escape TGFβ-mediated cytostatic effects by mutating SMAD4 or TGFβ receptor, while benefiting from the tumor-promoting effect of TGFβ on the stromal microenvironment (86).
Activation of IFN signaling has been recently proposed as a mechanism to increase immune-checkpoint blockade (ICB) response (87–89). Indeed, IFNsignhigh tumor cells expressed higher levels of several immune-checkpoint molecules (Supplementary Fig. S3F). Although we could not evaluate the functional contribution of the adaptive immune system to the IFNsignhigh/low scenario in the immunodeficient PDX mouse models, this observation suggests potentially targetable differences between the tumor subgroups. However, a sustained chronic IFN activation status, similar to the one we observe in IFNsignhigh tumors, has also been reported to cause ICB resistance (90). Thus, approaches to use ICB to potentially treat IFNsignhigh/low PDACs deserve further functional investigation. Possibly, efficient inhibition of intrinsic IFN signaling may be an approach to target these tumors with few side effects on normal cells (91).
The cell of origin in human PDAC has been debated for decades. PDAC mouse models originated from acinar cells (61, 92) and data from tracing experiments (60, 93) show that mutant acinar cells transdifferentiate to a ductal-like phenotype (with mixed characteristics of embryonic progenitor cells; ref. 67) to form precancerous lesions called pancreatic intraepithelial neoplasia (PanIN), which ultimately progress into PDAC. Mouse models where PDAC originates from ductal cells have been challenging as these cells seem more refractory to transformation, requiring additional insults, such as TP53 mutations, to progress toward PDAC (65, 66). However, as we and others showed, ductal-derived PDAC models associate with worse outcome and show significantly fewer and, if at all, only late-stage PanINs (63, 64). Starting from ductal cells, an acinar–ductal metaplasia would not be required, and classic PanIN formation may not occur in such PDAC subtypes. Indeed, little or no occurrence of PanIN lesions has also been observed in some patients with PDAC and was associated with poor survival (94–96). Future studies need to extend such observations to serial samples and association with specific molecular features of MC1- and MC2-type tumors. In addition to our analysis on resectable PDAC tumors, it will also be interesting to investigate these patterns in the spectrum of unresectable tumors as well as metastatic lesions. This would be of special interest in light of the poor effects achieved by the combination of ruxolitinib with capecitabine to treat patients with metastatic PDAC (97, 98).
Our understanding in pancreatic cancer biology has increased exponentially in recent years. Numerous studies have contributed to our knowledge of PDAC genetics and transcriptomics, and new studies start to shed light on the epigenomics of PDAC (11, 12, 99). Our data show that differential DNA methylation of repeat-derived genetic elements results in the execution of an epithelial cell–autonomous IFN expression program, which promotes a proinflammatory and tumor-promoting stromal feed-forward loop leading to increased aggressiveness of pancreatic cancer. Moreover, we suggest that the basis for Methylationlow/IFNsignhigh (MC2) and Methylationhigh/IFNsignlow (MC1) PDAC phenotypes is already laid down in the lineage identity of the ductal and acinar cells of the normal pancreas. Stable epigenetic traits imprinted in the cells of origin could be used to stratify different types of PDAC and may improve risk stratification and provide novel targeting opportunities.
Methods
PDAC and Normal Pancreas Samples for Epithelial Isolation
Tissue samples were obtained from patients who received partial pancreatoduodenectomy at the Department of General, Visceral and Transplantation Surgery, University of Heidelberg. Patients were part of the HIPO project. The study was approved by the ethical committee of the University of Heidelberg (case number S-206/2011 and EPZ-Biobank Ethic Vote #301/2001) and conducted in accordance with the Helsinki Declaration; written informed consent was obtained from all patients. Fixed sections were evaluated by two independent pathologists (W. Weichert and A. Muckenhuber), and only samples that did not contain normal tissue contamination were included in the study.
Normal samples were resected from sites distal to tumor lesions. A slide section from the processed sample was fixed for pathologic evaluation. Hematoxylin and eosin staining was performed in fixed sections, and slides were evaluated by two independent pathologists (W. Weichert and A. Muckenhuber). Only histologically normal samples were included in the study. Patient and tumor characteristics are summarized in Supplementary Table S1.
Donor pancreata were isolated and processed by the Beta Cell Bank of the JDRF Centre for Beta Cell Therapy in Diabetes (Brussels, Belgium) affiliated to the Eurotransplant Foundation, and written informed consent for use of donor material for research was obtained according to Belgian laws, or were obtained from organ procurement centers in the United States (including IIDP, NDPRI, IIAM, and nPOD) from deidentified postmortem samples and were excluded from the Institutional Review Board review per 45CFR56 and NIH policy.
Human Tissue Dissociation and Flow Cytometry Sorting
Tumor Tissue.
Tumor tissues were minced using a razor blade and dissociated using the human Tumor Dissociating Kit (Miltenyi Biotec) following the manufacturer's instructions. Dissociated samples were sequentially filtered through cell strainers of 100 μm → 70 μm → 40 μm (BD Falcon). To lyse erythrocytes, the cell pellet was resuspended in ACK lysing buffer (Lonza) and incubated 2 minutes at room temperature. After one wash with Hank's Balanced Salt Solution (HBSS; Thermo Fisher), cells were pelleted and resuspended in tumor staining buffer (HBSS + 1% BSA + 2 mmol/L EDTA) and blocked for 10 minutes at room temperature with FcR blocking reagent (Miltenyi Biotec). Next, the conjugated antibodies (see below) were added and the sample was stained for 15 minutes at 4°C.
Normal Tissue.
Normal samples were first dislodged by injecting fetal bovine serum (FBS) in the tissue. Next, tissues were minced using a razor blade. Pieces were placed in a 50 mL falcon with 7 mL of G Solution (ref. 100; 1 l HBSS + 0.9 g glucose + 47.6 μmol/L CaCl2) and incubated on ice for 1 minute. G solution was removed, and the process was repeated two to three times until the G Solution was clear. Pieces were added to the Digestion Mix Solution [Advanced DMEM/F12 + 25 mmol/L HEPES + 1% BSA + trypsin inhibitor (1:10; Sigma-Aldrich) + enzymes from the Tumor Dissociating Kit (Miltenyi Biotec)] and incubated at 37°C for 15 minutes. Digested tissue was homogenized by pipetting and sequentially filtered using cell strainers of 100 μm → 70 μm → 40 μm (BD Falcon). Digestion was stopped by adding one volume of FBS. If some tissue fragments were still left, the digestion was repeated and the filtered fractions were pooled. Cells were pelleted, resuspended in ACK lysing buffer (Lonza), and incubated 2 minutes at room temperature. Cells were washed using Wash Buffer (HBSS + 1% BSA + 1% FCS + 1 mmol/L EDTA + trypsin inhibitor + 2 mmol/L MgCl2) and pelleted. Pellet was resuspended in Normal Staining Buffer (Wash Buffer + 2 mmol/L CaCl2) + 2.5 μg/mL of DNase (Sigma-Aldrich) and FcR blocking reagent (Miltenyi Biotec) and incubated for 15 minutes at room temperature. Next, the conjugated antibodies (see below) were added, and the sample was stained for 15 minutes at 4°C.
Staining antibodies used were EPCAM-FITC (1:11 dilution; clone HEA-125; Miltenyi Biotec) and CD45-VioBlue (1:11 dilution; clone 5B1; Miltenyi Biotec). Propidium iodide was used to exclude dead cells. Epithelial cells were defined as EPCAM-FITC+/CD45-VioBlue−. CAFs were defined as EPCAM-FITC−/CD45-VioBlue−. An example of the gating strategy can be found in Supplementary Information.
For DNA extraction, sorted cells were collected into ice-cold CO2 independent medium (Thermo Fisher), pelleted, and stored at −80°C. For RNA extraction, cells were sorted into RNA lysis buffer (Arcturus PicoPure RNA Isolation Kit, Life Technologies) and stored at −80°C. Single-cell suspensions were sorted using a FACSFusion system (BD Biosciences).
DNA Isolation and TBWGBS
Genomic DNA from sorted cells was isolated using the QIAamp DNA Micro Kit (Qiagen) according to the manufacturer's instructions. Tagmentation-based WGBS was performed as described previously (101) using 30 ng genomic DNA as input. Four sequencing libraries with different barcodes were generated per sample and pooled in equimolar amounts to a 10 nmol/L pool. Each pool was sequenced paired-end, 125 bp, on two lanes of a HiSeq2000 sequencer (Illumina).
Mapping of Whole-Genome Bisulfite Sequencing Data and Methylation Calling
The TWGBS data were processed as described (14): The human reference genome (hg19) was in silico transformed for both the top strand (C to T) and bottom strand (G to A) using MethylCtools (55). Before alignment, trimming of adaptor sequences was performed using SeqPrep (https://github.com/jstjohn/SeqPrep). The first read in each read pair was then C to T converted, and the second read in the pair was G to A converted. The converted reads were aligned to a combined reference of the transformed top strands (C to T) and bottom strands (G to A) using BWA (bwa-0.6.2-tpx; ref. 102) with default parameters but disabling the quality threshold parameter for read trimming (-q) of 20 and the Smith-Waterman for the unmapped mate (-s). After alignment, reads were converted back to the original states, and reads mapping to the antisense strand of the respective reference were removed. Duplicate reads were marked, and the complexity was determined using Picard MarkDuplicates (http://broadinstitute.github.io/picard/). Reads with alignment scores of less than 1 were filtered before subsequent analysis. Total genome coverage was calculated using the total number of bases aligned from uniquely mapped reads over the total number of mappable bases in the genome. At each cytosine position, reads that maintain the cytosine status were considered methylated, and reads in which cytosine was converted to thymine were considered unmethylated. Only bases with Phred-scaled quality score of ≥ 20 were considered. In addition, the five base pairs at the two ends of the reads were excluded from methylation calling according to M-bias plot quality control. For the TWGBS libraries, the first nine base pairs of the second read and the final nine base pairs before the adaptor of the first read were excluded from methylation calling.
Calling of DMRs.
The raw counts of methylated and unmethylated reads for each CpG site from different libraries were merged for each sample. R (version 3.4.2) and BSmooth (ref. 103; version 1.14.0) were used with default parameters for smoothing and visualization purposes. Calling of DMRs was performed using DSS (ref. 104; version 2.26.0) with default parameters for comparison between normal and tumor samples as well as comparing MC1 and MC2.
Annotation.
Gene-related annotations were obtained from Gencode v19. Genomic annotations for DNase, CGIs, and repeats (SINEs, LINEs, and ERVs) were obtained from UCSC table browser. CGI shores were defined as 2 kb upstream and downstream flanking CGI. Chromatin states (15 states model) from normal pancreas from Roadmap data set (105) were used to annotate DMRs to enhancers. EnhG, EnhBiv, and Enh were merged and defined as “Enhancer.” The value for the overlap between one DMR and one genomic feature is the fraction of base pairs of the DMR that overlap to the genomic feature.
Overlap Differentially Methylated Genes from Different Studies
Calculation of Global Beta Values and Beta Values at CpG Islands and Repeats
Based on CpG island, LINE, SINE, and ERV annotations, mean beta values were calculated for each single feature in each sample. Beta values for each annotation category per sample were derived by calculating a mean beta value of all features of the same annotation category.
For the global beta value distribution plots, the distribution of genome-wide smoothed beta values per sample was plotted using frequency polygons. ggplot(…) + geom_freqpoly(binwidth = 0.01). In case of group-wise distribution, the mean of the group was calculated.
3-D Methylation-Based PCA
Based on all DMRs between normal versus tumor, a beta value matrix was created and a PCA using R (version 3.4.2) and prcomp(…) was performed (Fig. 1B). For an unbiased approach (Supplementary Fig. S3A), the genomes were organized into windows of 5 kb and filtered for a mean CpG coverage of 10 reads in each sample and window. The 10,000 most variable windows were used to perform PCA.
Methylation Density Plots and 2-D PCA
Density plots were generated by the densityHeatmap() function in the ComplexHeatmap package (107). For each column that represents one sample in the plot, colors were mapped to the density values in the corresponding distribution. The black dashed lines correspond to the five quantiles of the distributions and the red dashed lines correspond to the mean value of the distributions. For 2-D PCA plotting, CpG sites were randomly sampled by a probability of 0.001. For this subset of CpG sites, the top 2,000 CpG sites with the highest variance were selected to perform the PCA.
Estimation of PMDs
Genome-wide methylation profiles are represented as long vectors of numeric values, where the order of data points corresponds to chromosomal CpG positions. The Bayesian framework and dynamic programming approach of fastseg (version 1.20.0; ref. 108) were used to segment the genome based on the distribution of methylations. As fastseg does not take into account the distance between CpG sites but rather only considers the order of CpG sites, each chromosome was split into blocks if the gap between two adjacent CpG sites was larger than 100 kb. Blocks containing fewer than 50 CpG sites were removed.
Although fastseg can predict segment borders, each segment must still be classified based on its methylation. To do so, the distribution of both the mean and the standard deviation of methylation in the segments was examined (Supplementary Methods Fig. 1A). On each plot, there are at least two clear clusters, which correspond to unmethylated and fully methylated segments. We marked segments with mean methylation ≤ 0.1 and SD ≤ 0.15 as unmethylated and segments with mean methylation ≥ 0.7 and SD ≥ 0.2 as fully methylated (indicated by red rectangles in the above plots). The remaining segments were defined as intermediately methylated segments, and PMDs were detected among them.
Because the intermediately methylated segments are more heterogeneous (mean methylation ranging from 0.1 to 0.7), it is expected that there are clusters of segments, which are close to each other. The gap ratio was defined to measure the relative distance to neighboring segments in relation to the width of the segment itself. For a certain segment i, the gap ratio denoted as ri is defined as ri = di/wi, where di is the minimal distance to the two neighboring segments and wi is the width of the segment i. Segments with a smaller gap ratio tend to be close to the neighboring segments and need to be merged into a larger segment.
To select a proper cutoff of the gap ratio to merge neighboring segments, the distribution of gap ratios for all samples was computed (Supplementary Methods Figure).
On the left side of the distributions, the high peaks correspond to small gap ratios, where segments are highly close to each other. Thus, we set a cutoff of gap ratio ≤ 0.4, and if two segments were close enough, they were merged into a long segment keeping the gaps in between. Rainfall plotting (Supplementary Methods Figure) demonstrates that the process successfully merges the cluster.
Finally, the width distribution of intermediately methylated merged segments was computed, and a threshold for the width of a PMD was set to ≥ 10 kb.
Hilbert Plots
Hilbert curve visualization was performed by the HilbertCurve R package (109). The Hilbert curves for both PMDs and methylation profiles were constructed with level 9. The averaging mode for the methylation was “absolute.”
UpSet Plot
UpSet plots were generated by the UpSet() function of the ComplexHeatmap package. The intersection size of two sets of genomic regions was defined as the total number of base pairs of the intersected regions and the union size of two sets of genomic regions was defined as the total number of base pairs of the merged regions. Mode “distinct” was used to calculate the size of each combination of sets.
Isolation of Human Normal Acinar, Ductal, Dedifferentiated Acinar, and Mixed Exocrine Samples
Ductal and acinar cells from fresh tissue were isolated as previously described using HPx1 (acinar) and CD133 (ductal) markers (74).
To obtain ductal and dedifferentiated acinar cells from in vitro cultures, exocrine preparations were cultured as previously described (59, 70) in Advanced RPMI-1640 medium containing 5% heat-inactivated FBS and 1% penicillin–streptomycin solution (Life Technologies) under 5% CO2 atmosphere at 37°C. FITC-conjugated UEA-1 (Ulex Europeaus Agglutinin-1, Sigma) lectin labeling of acinar cells was performed according to Houbracken and colleagues and Baldan and colleagues (59, 70). Exocrine cell fraction labeled with UEA-1 was kept in suspension culture. At day 4 of culture, cell clusters were dissociated following the protocol of Baldan and colleagues (59). Cells were stained with carbohydrate antigen 19.9 (mouse monoclonal anti-human CA19.9, Dako). Alexa Fluor 647 anti-mouse (Jackson Laboratory) was used as secondary antibody. Analysis and cell sorting were performed on a BD FACSAria (BD Biosciences). Viable, single cells were gated based on forward and side scatter.
Infinium MethylationEPIC DNA Methylation Arrays
Genomic DNA was extracted from two independent replicates of HPDE cells using the DNeasy Blood and Tissue Kit (Qiagen) according to the manufacturer's instructions. PDX genomic DNA was extracted using Allprep DNA/RNA/Protein Mini Kit (Qiagen). Genomic DNA from acinar, ductal, dedifferentiated acinar, or mixed exocrine samples was extracted from frozen pellets using the QIAamp DNA Micro Kit (Qiagen) according to the manufacturer's instructions. Additionally, three commercial genomic DNA preps were used [#CD564011 (Origene), #CD564404 (Origene), and #A712202 (BioChain)].
DNA methylation data were generated using Infinium MethylationEPIC BeadChip arrays at the Genomics and Proteomics Core Facility of the DKFZ or were previously available (GSE122126; ref. 71). Processing of the data was done with the RnBeads Package (version 2.4.0) for R (110) using default settings. Beta values were normalized using the BMIQ (Beta MIxture Quantile dilation; ref. 111) method. Hierarchical clustering was performed by the pheatmap() function using default parameters. Annotation to acinar and ductal enhancers was done using previously published ChromHMM data from isolated human acinar and ductal cells (69). Calculations were performed with R version 3.6.0 in R-Studio.
RNA Isolation and RNA-seq
Total RNA isolation and RNA-seq were performed as previously described (13): RNA was extracted using the Arcturus PicoPure RNA Isolation Kit (Life Technologies) including a DNase treatment step (Qiagen) following the manufacturer's instructions. cDNA was generated with 1 ng of total RNA and 13 cycles of amplification using the SMARTer Ultra Low RNA Kit for Illumina sequencing (Clontech) according to the manufacturer's instructions. The sequencing libraries were generated using the NEB Next ChIP-Seq Kit (New England Biolabs) according to the manufacturer's instructions. The quality of total RNA, cDNA, and library was monitored throughout the process using a 2100 Bioanalyzer (Agilent). Sequencing was done in a HiSeq2000 sequencer (Illumina) with three samples per lane, paired-end 100 bp.
RNA-seq Data Processing
RNA-seq reads were aligned to the hg19 reference genome using the STAR alignment software (version 2.3.0e; ref. 112), with the merged transcriptome of Gencode v19 and lincRNA catalog annotations (113). Counts of reads mapped to exons were estimated by htseq-count (114). Expression values normalized using DEseq2 (115) were used for heat-map visualization. On the heat map, genes were scaled by z-score scaling. Heat maps were generated by the ComplexHeatmap R package (107). TPM values were calculated by normalizing to the sum of the length of nonoverlapping exons for each gene and to the library size for each sample.
For expression of repeats, featureCount (version 1.5.3; ref. 116) was applied on the merged annotation of Gencode v19 and RepeatMasker (UCSC) considering three classes or repetitive elements: SINEs, LINEs, and LTRs. The –fracOverlap parameter was set to 0.5 to enforce that each given single read is counted as in a repeat only if more than 50% base pairs of the read overlap to it. TPM normalization was applied to repeats expression using a library size estimated by the total reads mapped to genes. The expression difference between groups 1 and 2 was defined as (μ1 − μ2)/σ, where μ1 and μ2 are the mean in groups 1 and 2, respectively, and σ is the standard deviation by pooling values in the two groups.
The association of samples to previously described PDAC subtypes was based on published signatures (Supplementary Table S8; refs. 2, 24, 117) using the following method: Assume a signature matrix has n signature genes and K subgroups. The signature scores in subgroup k are presented as a vector denoted as sk. For sample i in another expression matrix where the expression vector is denoted as xi with the same set of genes as the signatures (also in the same order), we calculate the Euclidean distance denoted as dik between xi and sk. Over all signature vectors, dik forms a vector denoted as di. We order the vector di decreasingly and take the first two highest values di(1) and di(2) and define the difference score as
where μ is the mean of vector di. To get a significantly large ai, data points in each sk are randomly permutated nr times (e.g., nr = 1,000), following the process above. The random difference score in permutation j is calculated denoted as aijrandom. The P value for the difference score to be significantly large is pi = 1/1,000 j∑1,000I(rijrandom > ri), where I(expr) = 1 if expr is evaluated as true, or else I(expr) = 0. If pi < 0.05, the subgroup with the smallest distance to sample i is assigned, or else sample i is labeled as no subgroup is assigned.
Correlation Methylation-Expression
For each DMR, the Spearman correlation test was applied to the mean methylation of the DMR and the expression value of the gene with the nearest TSS to that DMR. The significant correlations were filtered by FDR ≤ 0.05.
Analyses of Public Data Sets
PDAC Data.
Normalized expression data, subtype annotation, tumor cellularity (QPURE) estimation, and survival data were obtained from the Supplementary Information accompanying the Bailey and colleagues publication (25). Normalized expression and survival data were downloaded from the NCBI Gene Expression Omnibus (GEO) with accession number GSE21051. The Cancer Genome Atlas (TCGA) pancan normalized values from PDAC samples were downloaded from the TCGA database, and non-PDAC samples were excluded prior to the analyses. TCGA subtyping was used as defined in Connor and colleagues (27). Kaplan–Meier curves were calculated with GraphPad Prism 8 software (Bailey and colleagues, GSE21051; ref. 25) or in KM Plotter (TCGA; ref. 118). For correlation between expression of genes and of genes with IFNsign, the average IFN signature expression for each patient was calculated, and correlations were computed via Spearman correlation (R).
Normal Acinar and Ductal RNA Expression Data.
Normalized RPKM values for bulk RNA-seq data of acinar and ductal cells were obtained from accession number GSE57973. Single-cell RNA-seq preprocessed data as raw count values for the studies with the accession numbers GSE81547, E-MTAB-5061, and GSE85241 were downloaded and further processed per study using R (version 3.4.2) and Seurat (ref. 119; version 2.3.0). The analysis used healthy control samples only. Log normalization was applied using NormalizeData() and highly variable genes defined using FindVariableGenes() with default parameters. Data were further normalized using linear regression via the function ScaleData(), regressing for the number of unique molecular identifiers in all data sets as well as patient origin for GSE81547, patient origin for E-MTAB-5061, and patient and library origin for GSE85241. For downstream analysis, the first 10 principal components were used. Cells were defined as ductal or acinar based on authors' annotation, for study GSE81547 and E-MTAB-5061. In case of unavailable cell annotation (GSE85241), cells were clustered into acinar or ductal using FindClusters(dims.use = 1:10, resolution = 0.8), and the expression of genes was defined in the original study as markers of acinar or ductal cell. A mean expression value for each cell was calculated based on an IFN signature and individual genes were visualized as violin plots, additionally applying a Mann–Whitney test.
Acinar and Ductal Methylation Data.
EPIC methylation array data from isolated acinar and ductal cells were obtained from accession number GSE122126, and normalized and processed together with the data generated in house.
Human Cells
Tumor PDCs have been previously reported (PACO2, PACO17, PACO14, and PACO18; ref. 39) or have been similarly derived (PACO3, PACO20, and PACO28). PDCs were maintained in cancer stem cell (CSC) medium, which contains advanced DMEM–nutrient mixture F-12 (DMEM/F12; Life Technologies) with N2 supplement (Life Technologies), 50 ng/mL basic fibroblast growth factor (bFGF; PeproTech), 20 ng/mL epidermal growth factor (EGF; PeproTech), 10 ng/mL LONG R3 insulin-like growth factor–1 (IGF1; Sigma), 100 μmol/L β-mercaptoethanol (Life Technologies), 2 μg/mL heparin (Sigma) and grown in Primaria plates (BD) and used at maximum passage 15. PDCs were categorized according to the expression of IFN signature genes as: IFNsignhigh: PACO2, PACO17, and PACO20. IFNsignlow: PACO3, PACO14, PACO18, and PACO28.
Human pancreatic StCs were purchased from ScienCell Research Laboratories and cultured using poly-L-lysine–coated (Sigma-Aldrich) tissue culture flasks and SteCM medium (ScienCell) with 2% serum. StC were used for a maximum of six passages.
HPDE cells (73) were obtained from the ATCC, cultured according to the supplier's information, and used at passage 14.
All cell lines used were authenticated monthly by single-nucleotide polymorphism profiling and tested for Mycoplasma contamination (both by Multiplexion).
Conditioned Medium.
PDCs were seeded in CSC complete medium. After 24 hours, cells were washed once with PBS (Sigma-Aldrich) and incubated in CSC-reduced medium (CSC minus bFGF, EGF, and IGF) for 1 hour. Subsequently, fresh CSC-reduced medium was added for 48 hours. CM was harvested and spinned at 300 × g for 10 minutes to remove cellular debris. PDCs were harvested and counted to ensure equal cell numbers. StCs were preincubated for 1 hour with CSC-reduced medium and treated with PDC CM or StC CM for 12 hours prior to RNA extraction. For long-term conditioning, StCs were treated with PDC CM or StC CM every 3 days during 10 days.
Cocultures.
PDCs were seeded in 0.4-μm pore transwells (Becton Dickinson). Cell numbers were adapted experimentally to ensure equal values during the experiment. The same number of cells was seeded in a transparent 96-well plate to monitor cell confluence during the experiment. Twenty-four hours later, pancreatic StCs were seeded onto companion plates (Becton Dickinson) coated with poly-L-lysine (Sigma-Aldrich). The next day, StCs and PDCs were preincubated in CSC-reduced medium for 1 hour prior to setting the coculture. After 2 days of coculture, 50% of the medium was replaced with fresh CSC-reduced medium. After 4 days of coculture, transwells were placed on a new plate, and cell amount was estimated using CellTiter-Blue stock solution (Promega) following the manufacturer's instructions.
Ruxolitinib (Focus Biomolecules) was added daily at a final concentration of 100 nmol/L.
Human IFN lambda receptor 1 (IFNLR1) extracellular domain (clone MMHLR-1; PBL assay science) or Mouse IgG1 Iso Control (clone P3.6.2.8.1; eBioscience) was used at a final concentration of 0.5 μg/mL. IFNLR1 or IgG control was added on days 1 and 3 of coculture. Recombinant IL28A or IL29 was used at 10 ng/mL.
5-Aza-CdR Treatment.
PDCs were treated daily with 500 nmol/L of 5-Aza-2′-deoxycytidine (5-Aza-CdR; Santa Cruz) or DMSO as control for 5 days. On day 6, proteins and RNA were extracted.
Cytokine Treatment.
PDCs were seeded in Primaria 96 wells (Corning) and treated every 2 days for a total of 5 days with recombinant human IL8, IL1A, or IL1B (PeproTech) or 0.1% BSA PBS. Cell amount was estimated using CellTiter-Blue stock solution (Promega) following the manufacturer's instructions.
IFN ELISAs.
CM from 5 × 106 cells was collected after 48 hours. Cellular debris was excluded by centrifugation at 1.500 rpm. Twelve milliliters of tumor media or only media as control was concentrated to 1.2 mL using Amicon columns (3,000 Nominal Molecular Weight Limit; Merck Millipore) following the manufacturer's instructions. The following ELISA kits were used as indicated by the manufacturer: Human IFN-alpha Platinum ELISA (BMS216; Affymetrix Bioscience), VeriKine Human IFN Beta ELISA Kit (41410, pbl assay science), Human INFγ ELISA Kit (EHIFNG; Invitrogen), Human IL29/IFN-λ1 (DY7246; R&D Systems), and Human IL28A/IFN-λ2 (DY1587; R&D Systems).
Generation of Knockdown Cells
STAT1 and MAVS knockdowns were generated with a miR-E (120) StagBFPEP lentiviral vector, a modified version of the original SGEP vector kindly provided by Johannes Zuber (IMP—Research Institute of Molecular Pathology GmbH, Vienna) in which the constitutively expressed GFP protein is replaced by the tagBFP protein. miR-E sh oligonucleotides were designed using the shERWOOD algorithm tool (121): Nonsilencing control: 5′-TGCTGTTGACAGTGAGCGCTCTCGCTTGGGCGAGAGTAA GTAGTGAAGCCACAGATGTACTTACTCTCGCCCAAGCGAGATTGCCTACTGCCTCGGA-3′; shSTAT1 #1: 5′-TGCTGTTGACAGTGAGCGCCAGAAAGAGCTTGACAGTAA ATAGTGAAGCCACAGATGTATTTACTGTCAAGCTCTTTCTGTTGCCTACTGCCTCGGA-3′; shSTAT1 #2: 5′-TGCTGTTGACAGTGAGCGCTCAGAGCACAGTGATGTTAGA TAGTGAAGCCACAGATGTATCTAACATCACTGTGCTCTGAATGCCTACTGCCTCGGA-3′; shMAVS #1: 5′-TGCTGTTGACAGTGAGCGCCCAAGTTGCCAACTAGCTCAAT AGTGAAGCCACAGATGTATTGAGCTAGTTGGCAACTTGGATGCCTACTGCCTCGGA-3′; shMAVS #2: 5′- TGCTGTTGACAGTGAGCGAACCAATCCAGCACCATCCAAAT AGTGAAGCCACAGATGTATTTGGATGGTGCTGGATTGGTGTGCCTACTGCCTCGGA-3′. Oligonucleotides were amplified by PCR using the Q5 High-Fidelity DNA Polymerase (New England Biolabs). The following primers were used for PCR amplification: miRE-Xho-fw: 5′-TGAACTCGAGAAGGTATATTGCTGTTGACAGTGAGCG-3′; miRE-EcoOligo-rev: 5′-TCTCGAATTCTAGCCCCTTGAAGTCCGAGGCAGTAGGC-3′. PCR products containing shGene and nonsilencing miR-Es were subcloned into the StagBFPEP recipient vector via EcoRI-HF and XhoI restriction sites. Lentiviral particles were produced by transfecting HEK293T cells with StagBFPEP-miR-E vectors together with plasmids encoding for the packaging proteins pMD2G and psPAX2 using polyethylenimine (Polyscience). Cells were transduced at a multiplicity of infection 5.
Gene-Expression Microarrays
Total RNA from PDCs or treated StCs was isolated using the RNeasy Kit (Qiagen) according to the manufacturer's instructions. Gene-expression analyses were performed using Illumina HumanHT12v4 BeadChips at the Genomics and Proteomics Core Facility of the DKFZ.
GSEA
GSEA (23) was conducted using the GSEA desktop application, and the gene sets were downloaded from the Broad Institute (Hallmarks) or indicated in Supplementary Information using 2,000 permutations. For microarray data, quantile-normalized data were used as input. For RNA-seq data, ranked log2 fold-change values calculated by DESeq2 (115) were used. Genes with NA FDR values were excluded from the analyses. IFNsign comprises genes from Interferon-alpha response and Interferon-gamma response Hallmark gene sets with a rank metric score > 1 in MC2 versus MC1 samples.
Gene sets used in the study are provided in Supplementary Table S8.
Immunofluorescence
PDCs cells were seeded on eight-chamber culture slides (Falcon) coated with poly-L-lysine (Sigma-Aldrich) and grown to 50% to 70% confluency. Cells were fixed in 4% freshly depolymerized formaldehyde for 10 minutes, permeabilized with 0.25% (vol/vol) Triton X-100 (Sigma) for 20 minutes, and blocked with 1% BSA for 1 hour. J2 antibody for dsRNA (ref. 122; Scicons J2 from English and Scientific Consulting Kft; 10010200) or mouse IgG2a k isotype control (eBM2a; eBioscience) was incubated overnight at 4°C (20 μg/mL) and detected by fluorescence using a goat anti-mouse IgG, F(ab')2 dylight 649 (1:300 dilution; Jackson ImmunoResearch Laboratories). Slides were mounted using ProLong Antifade GOLD with DAPI (Life Technologies), as described by the manufacturer. Pictures were taken using a Nikon Eclipse Ti microscope.
FACS dsRNA
PDCs (100,000) were fixed with Cytofix/Cytoperm solution (BD) for 10 minutes on ice. After washing with Permwash (BD), cells were stained with J2 antibody (Scicons J2 from English and Scientific Consulting Kft; 10010200) or mouse IgG2a k isotype control (eBM2a; eBioscience) overnight at 4°C (2 μg/mL). The day after, cells were washed with Permwash (BD) and stained using a goat anti-mouse conjugated to APC (1:200 dilution; R&D Systems, F0101B). Quantification of the difference between the median fluorescent intensity of J2 staining and isotype staining was calculated on singlets for each cell line (ΔMFI). An example of the gating strategy is shown in Supplementary Information. Two independent experiments were performed per cell line except for PACO20, for which only one replicate was available. FACS analyses were performed on an LSR2 (BD Biosciences), and data analysis was done with FlowJo software version 10.0.3 (FlowJo, LCC).
Western Blot
Whole-cell lysates of PACO cells were prepared using RIPA lysis buffer (Cell Signaling Technology) with 1 mmol/L PMSF (Sigma), 1 mmol/L EDTA, and 1 mmol/L Halt Protease-Phosphatase Inhibitor Cocktail (Pierce). Protein concentration was determined by the Pierce BCA assay (Thermo Scientific). Proteins were resolved on 4% to 12% TGX gels (Criterion, Bio-Rad) with TGS (Tris-Glycine-SDS) running buffer (Bio-Rad) and blotted on PVDF membranes (Trans-Blot TURBO, Bio-Rad). Membranes were blocked for 1 hour in TBS containing 0.1% (v/v) Tween-20 and 5% (w/v) BSA. Primary antibodies against MDA5 (D74E4; 1/1,000; Cell Signaling Technology, #5321), STAT1 (42H3; 1/1,000; Cell Signaling Technology, #9175), phospho-STAT1 (Tyr701; 58D6; 1/1,000; Cell Signaling Technology, #9167), phospho-STAT1 (Ser727; 1/1,000; Cell Signaling Technology, #9177), and MX1 (E-5; 1/1,000; Santa Cruz, sc-271024) were incubated overnight at 4°C in blocking solution. Secondary antibodies [anti-rabbit-IgG horseradish peroxidase–linked antibody (1/10,000 in TBS containing 0.3% Tween-20, Cell Signaling Technology, #7074) or mouse IgG binding protein horseradish peroxidase (HRP; 1/10,000 in TBS containing 0.1% Tween-20, Santa Cruz, sc-516102)] were incubated for 1 hour at room temperature. Membranes were washed in TBS-Tween 0.1%, and immunocomplexes were detected using the ECL kit (Amersham International) and acquired with a ChemiDoc XRS+ System (Bio-Rad). Uncropped images are shown in Extended Data Fig. 3.
Real-Time Quantitative PCR
Total RNA was extracted using the miRNeasy mini kit (Qiagen) and reverse transcribed using the high-capacity cDNA reverse transcription kit (Applied Biosystems) combined with RNase-Free DNase (Qiagen) treatment according to the manufacturer's protocol. RT-qPCR was performed and analyzed on ViiA7 (Applied Biosystems) using Power SYBR Green PCR Master Mix (Applied Biosystems; Supplementary Fig. S10A; Fig. 3I; Supplementary Figs. S10B and S8H) or TaqMan gene-expression assay (Applied Biosystems; Supplementary Figs. S11A, S11E, S11I, and S10C). The following primers and TaqMan probes (Applied Biosystems) were used to acquire expression data with the Viia 7 Real-Time PCR System (Applied Biosystems): Primers: HPRT1 (Fw: 5′-CCTGGCGTCGTGATTAGTGAT-3′; Rev: 5′-AGACGTTCAGTCCTGTCCATA-3′); OAS2 (Fw: 5′-GGAGCTTCCTGATTGGCAGA-3′; Rev: 5′-ATGTAGGGTGGCAAGCACTG-3′); MX1 (Fw: 5′-CTGCACAGGTTGTTCTCAGC-3′; Rev: 5′-GTTTCCGAAGTGGACATCGCA-3′); STAT1 (Fw: 5′-CAGCTTGACTCAAAATTCCTGGA-3′; Rev: 5′-TGAAGATTACGCTTGCTTTTCCT-3′); MAVS (Fw: 5′-AGGAGACAGATGGAGACACA-3′; Rev: 5′-CAGAACTGGGCAGTACCC-3′); TaqMan Probes: PPIA (Hs04194521_s1), GAPDH (Hs99999905_m1), HPRT1 (Hs02800695_m1), MX1 (Hs00895608_m1), IFI44L (Hs00915292_m1), CCL2 (Hs00234140 m1), OAS2 (Hs00942643_m1), STAT1 (Hs01013996_m1), IL8 (Hs00174103_m1), IL1A (Hs00174092_m1), and IL1B (Hs01555410_m1).
Mice
NOD.Prkdcscid.Il2rgnull (NSG) mice were bred and housed under specific pathogen-free conditions at the central animal facility of the German Cancer Research Center (DKFZ). Female mice were used for the studies. All animal experiments were approved by the governmental committee for animal experimentation (Regierungspräsidium Karlsruhe; G39/13, G115/19, and G213/18).
Orthotopic Experiment.
For the orthotopic tumor growth of PDC experiments (Fig. 3C, J, and L), 200,000 cells were mixed with Matrigel (2 mg/mL; BD) and injected into the mice's pancreas. Engraftment of tumors and subsequent growth were monitored by regular palpation of the implantation site. For the experiment shown in Fig. 3C, mice were killed as soon as they reached abortion criteria defined in the procedure (human endpoint).
StC/PDC Coinjection.
50,000 IFNsignlow cells (PACO3) were mixed with 100,000 StCs, which had been previously conditioned with StC or PDC medium and injected subcutaneously in HBSS. Mice were killed 90 days after injection.
Ruxolitinib Treatment.
100,000 IFNsignhigh (PACO2) or IFNsignlow (PACO3) PDCs were subcutaneously injected in mice. When tumors reached a size of approximately 50 mm3 (26 or 60 days, respectively), mice were randomly assigned to each treatment group (n = six per group; one tumor from the PACO2 ruxolitinib group was excluded from the analysis because no data on initial size of the tumor could be taken). Vehicle (1% DMSO in 0.1% NaCl) or ruxolitinib (8.7 mg/Kg in 0.1% NaCl; TargetMol) was administered daily by peritoneal injection during 13 days. Mice were killed 3 days after last administration. Tumors were measured with a caliper, and volume was calculated according to the formula (length × height × width) × (π/6). Fold increase was normalized to size at starting of treatment.
1° PDX.
Primary xenografts were established by implanting 1 to 2 mm3 pieces of patient tumors onto the pancreas of NSG mice.
Mouse PDAC Cells
To generate mouse PDAC cell lines, tumors identified at necropsy in Sox9-CreER;KrasLSL-G12D;Trp53f/f;R26RYFP and Ptf1aCreER;KrasLSL-G12D;Trp53f/f;R26RYFP mice were dissociated and plated using previously described culture conditions (100). RNA was isolated from low-passage cells grown on a collagen gel, and RNA-seq libraries were generated using the TruSeq kit (Illumina) on a Neoprep system (Illumina) and then sequenced using a NextSeq device (UBC Biomedical Research Centre Next-Generation Sequencing Service) at 30 to 50 × coverage. Paired-end reads were aligned to the mouse genome (UCSC mm9) using STAR alignment (112), and fragments per kilobase of transcript, per million mapped reads (FPKM) data were generated with Cufflinks on the Basespace (Illumina) platform.
IHC
Tumor specimens were fixed in 10% formalin overnight and embedded in paraffin. For IHC, slides were deparaffinized and rehydrated. Antigen retrieval was enhanced by boiling in a steam pot at pH 6 in Dako Target Retrieval Solution (Dako) for 15 minutes, followed by cooling for 30 minutes and washing in distilled water. Nonspecific binding was blocked using the Linaris Avidin/Biotin Blocking Kit (Vector Labs) according to the manufacturer's instructions. Slides were incubated with primary antibodies for 30 minutes, rinsed in PBS-T (PBS with 0.5% Tween-20), incubated for 20 minutes with the appropriate secondary antibody using the Dako REAL Detection System (Dako), and rinsed in PBS-T. After blocking of endogenous peroxidase and incubation with streptavidin–HRP (20 minutes at room temperature), slides were developed with AEC (Dako) and counterstained with hematoxylin. All antibodies were diluted in Dako antibody diluent. Anti-STAT1 (dilution 1:100; HPA000982, Sigma) and anti-SERPINB5 (dilution 1:75; HPA020136, Sigma).
Statistical Analysis
For the in vitro experiments, indicated biological replicates were used or grouped analyses were carried out. For the in vivo treatment experiments, six mice per group were used. Hence, for the reported differences, the sample size used gave sufficient power to call them reliable. For the in vivo coinjection, the number of injected mice was chosen according to the availability of cells and, therefore, grouped analysis was carried out. Quantitative results were analyzed by one-way analysis of variance (one-way ANOVA), two-sided unpaired Mann–Whitney U test, or two-sided unpaired Student or nested t test as indicated in the figure legends, using GraphPad Prism 8 software. For RNA-seq expression data of sorted epithelia cells, two-tailed Wald test was calculated by DEseq2 package (115). Survival analysis was performed by using the Mantel–Cox log-rank test and the hazard ratio (log-rank) test using the GraphPad Prism 8 software. Estimation of variation within each group was determined by SEM or 95% CI as indicated in the figure legends. Levels of significance were determined as follows: *, P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001.
Data Availability
Gene-expression microarray data of PDCs, StCs treated with CM, normalized Infinium MethylationEPIC data, and normalized WGBS data have been deposited in the GEO database under accession numbers GSE134092, GSE134093, GSE134217, and GSE161956, respectively. Processed high-throughput data (i.e., full DMR calling and differential gene expression from RNA-seq data) are provided as Supplementary Tables. Sequence data and methylation raw files have been deposited at the European Genome–Phenome Archive (EGA), which is hosted by the European Bioinformatics Institute and the Centre for Genomic Regulation, under accession number EGAS00001004660. Other raw data are available upon reasonable request.
Authors' Disclosures
E. Espinet reports grants from EMBO (post-doctoral fellowship) during the conduct of the study. N.A. Giese reports grants from Heidelberger Stiftung Chirurgie and BMBH during the conduct of the study. M. Safavi reports personal fees from F. Hofmann-La Roche (M. Safavi is employed by F. Hoffmann-La Roche) outside the submitted work. E. Backx reports grants from Fund for Scientific Research Flanders (FWO; FWO Odysseus I) outside the submitted work. M. Schlesner reports grants from German Ministry of Science and Health during the conduct of the study. N. Pfarr reports personal fees from BMS, Novartis, Illumina, Thermo Fisher Scientific, and Roche outside the submitted work. B. Brors reports grants from German Federal Ministry of Research and Education (BMBF) “de.NBI-epi (German Network for Bioinformatics Infrastructure: partner project for epigenetics); FKZ 031L0101A” during the conduct of the study. I. Rooman reports grants from Fund for Scientific Research Flanders (FWO; FWO Odysseus I) during the conduct of the study. W. Weichert reports grants from Roche, MSD, BMS, AstraZeneca, and Bruker and personal fees from Roche, MSD, BMS, AstraZeneca, Pfizer, Merck, Lilly, Boehringer, Novartis, Takeda, Bayer, Amgen, Astellas, Illumina, NewOncology, and Agilent outside the submitted work. M.R. Sprick reports grants from BMBF during the conduct of the study; in addition, Dr Sprick has patents for US2017260593 (A1) and US2019317097 (A1) issued. A. Trumpp reports grants from Dietmar Hopp Foundation, German Cancer Consortium (DKTK), and Bundesministerum für Bildung und Forschung (BMBF) during the conduct of the study. No disclosures were reported by the other authors.
Authors' Contributions
E. Espinet: Conceptualization, data curation, formal analysis, supervision, funding acquisition, investigation, visualization, methodology, writing–original draft, project administration, writing–review and editing. Z. Gu: Resources, data curation, software, formal analysis, investigation, visualization, methodology, writing–review and editing. C.D. Imbusch: Resources, data curation, software, formal analysis, investigation, visualization, writing–review and editing. N.A. Giese: Resources, funding acquisition, writing–review and editing, screened, selected, and included patients, provided human tissue specimens and advice regarding human samples. M. Büscher: Formal analysis, investigation, writing–review and editing. M. Safavi: Formal analysis and investigation. S. Weisenburger: Formal analysis and investigation. C. Klein: Formal analysis and investigation. V. Vogel: Investigation. M. Falcone: Software, formal analysis, investigation, writing–review and editing. J. Insua-Rodriguez: Formal analysis, investigation, writing–review and editing. M. Reitberger: Software, investigation, writing–review and editing. V. Thiel: Investigation. S.O. Kossi: Investigation. A. Muckenhuber: Investigation and pathological assessment. K. Sarai: Resources, investigation, and generated mouse cell line data. A.Y.L. Lee: Resources and generated mouse cell line data. E. Backx: Resources and provided normal acinar and ductal samples. S. Zarei: Resources and generated mouse cell line data. M.M. Gaida: Investigation, pathological assessment. M. Rodriguez-Paredes: Formal analysis, investigation, writing–review and editing. E. Donato: Formal analysis, investigation, writing–review and editing. H.-Y. Yen: Investigation, pathological assessment. R. Eils: Resources, funding acquisition, and investigation. M. Schlesner: Resources, funding acquisition, writing–review and editing. N. Pfarr: Investigation, pathological assessment. T. Hackert: Resources, investigation, screened, selected, and included patients, provided human tissue specimens and advice regarding human samples. C. Plass: Resources, writing–review and editing. B. Brors: Resources. K. Steiger: Resources, investigation, pathological assessment. D. Weichenhan: Investigation, methodology, and prepared the tagmentation libraries. H.E. Arda: Resources, methodology, provided normal acinar and ductal samples. I. Rooman: Resources, writing–review and editing, provided normal acinar and ductal samples. J.L. Kopp: Resources, writing–review and editing, and generated mouse cell line data. O. Strobel: Resources, funding acquisition, writing–review and editing, screened, selected, and included patients, provided human tissue specimens and advice regarding human samples. W. Weichert: Resources, funding acquisition, writing–review and editing. M.R. Sprick: Resources, funding acquisition, project administration, writing–review and editing, pathological assessment. A. Trumpp: Conceptualization, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.
Acknowledgments
We thank the EPZ-Biobank and surgery teams, especially M. Schenk and M. Fischer, for excellent technical assistance regarding tissue specimens and M. Bähr for excellent technical assistance on TBWGBS library preparation. We thank the next-generation sequencing (NGS) unit of the Genomics Core Facility at DKFZ for providing NGS and related services, the microarray unit and all members of the Flow Cytometry Core Facility at DKFZ for excellent support, and the members of the Central Animal Laboratory at DKFZ for animal husbandry. We thank DKFZ-HIPO for technical support and funding. We thank ODCF, especially J. Kerssemakers, for support on data deposition. We thank A. Descot, M. Saini, F. M. Zickgraf, and C. Eisen for discussions and U. Rochlitzer for help with translating the mouse experimental protocols. We thank L. Wang and T. TruongVo for help deriving the acinar and ductal human samples. This work was supported by HIPO-015A (A. Trumpp., M.R. Sprick, E. Espinet, O. Strobel, T. Hackert, R. Eils, and W. Weichert); the Dietmar Hopp Foundation (A. Trumpp); the German Ministry of Science and Education (BMBF) e:Med program for systems biology (PANC-STRAT consortium, grant no. 01ZX1305C and 01ZX1605C; A. Trumpp, M.R. Sprick, E. Espinet, R. Eils, O. Strobel, and W. Weichert); the German Ministry of Science and Education (BMBF; grant no. 031L0076A; B. Brors); the DKFZ-NCT program NCT3.0 (A. Trumpp, M.R. Sprick, M. Schlesner, R. Eils, N.A. Giese, and O. Strobel—NCT3.0_2015.17 PrecO-Panc); the DKTK Joint Funding Program (A. Trumpp, M.R. Sprick, and W. Weichert); the Pancreatic Cancer Canada Foundation Innovation Grant, Pancreas Centre BC, CIHR Open Operating and New Investigator grants to J.L. Kopp, and a University of British Columbia Faculty of Medicine Graduate Award (A.Y.L. Lee); the Deutsche Forschungsgemeinschaft (DFG, German Research Foundation, Project ID 329628492—SFB 1321, S02P; K. Steiger, and W. Weichert); the German Research Foundation (M.M. Gaida; GA 1818/2-1). The collection and processing of the specimens via PancoBank was supported by Heidelberger Stiftung Chirurgie. I. Rooman is a recipient of an Odysseus I fellowship of the Fund for Scientific Research Flanders (FWO). E. Espinet was a recipient of an EMBO long-term fellowship (ALTF 344-2013).