Genome-wide DNA methylation analysis of KRAS mutant cell lines

Oncogenic RAS mutations are associated with DNA methylation changes that alter gene expression to drive cancer. Recent studies suggest that DNA methylation changes may be stochastic in nature, while other groups propose distinct signaling pathways responsible for aberrant methylation. Better understanding of DNA methylation events associated with oncogenic KRAS expression could enhance therapeutic approaches. Here we analyzed the basal CpG methylation of 11 KRAS-mutant and dependent pancreatic cancer cell lines and observed strikingly similar methylation patterns. KRAS knockdown resulted in unique methylation changes with limited overlap between each cell line. In KRAS-mutant Pa16C pancreatic cancer cells, while KRAS knockdown resulted in over 8,000 differentially methylated (DM) CpGs, treatment with the ERK1/2-selective inhibitor SCH772984 showed less than 40 DM CpGs, suggesting that ERK is not a broadly active driver of KRAS-associated DNA methylation. KRAS G12V overexpression in an isogenic lung model reveals >50,600 DM CpGs compared to non-transformed controls. In lung and pancreatic cells, gene ontology analyses of DM promoters show an enrichment for genes involved in differentiation and development. Taken all together, KRAS-mediated DNA methylation are stochastic and independent of canonical downstream effector signaling. These epigenetically altered genes associated with KRAS expression could represent potential therapeutic targets in KRAS-driven cancer.


CpG methylation in a panel of 47 cell lines shows clustering of cell lines with similar tissue of origin independent of KRAS mutation status.
Given the essential role of oncogenic KRAS in the great majority of pancreatic cancer 15,29 (see cell line information, Supplementary Fig. S1), we investigated whether the presence of an activating KRAS mutation correlates with specific patterns of global DNA methylation. We first performed genome-wide DNA methylation profiling of 11 KRAS-dependent pancreatic cancer cell lines using the Infinium HumanMethylation450 BeadChip Array 31 . We also surveyed the CpG methylation patterns in low passage, immortalized lung epithelial cells transduced with KRAS G12V (SAKRAS cells) and non-transformed empty vector controls (SALEB cells). We compared the panel of 11 KRAS-mutant pancreatic cancer cell lines to DNA methylation data collected from SALEB and SAKRAS lung epithelial cells and published Infinium methylation data from ENCODE 32 (Fig. 1). The published ENCODE data include three non-transformed human cell lines (HGPS and IMR-90 fibroblasts, and two different MCF 10 A breast epithelial cell lines) and 30 cell lines of varying cell types, genetic backgrounds, and tumorigenicity. As the pancreatic cancer cell lines were transduced with non-silencing (NS) shRNA, which could potentially affect the methylome of the transduced cells, we performed the same analysis while excluding these cells ( Supplementary Fig. S2). After unsupervised hierachial clustering of the top 1,000 most variable CpG probes across all 47 cell lines, the pancreatic cancer cell lines formed a distinct cluster with the exception of CFPAC-1_NS and PANC-1_NS cells. These data suggest that the panel of KRAS-mutant pancreatic cancer cell lines contain similar overall basal DNA methylation patterns. Other KRAS mutant lines were clustered in the same branch of the dendrogram. However, in general, the cell lines formed clusters based on cell type with a few exceptions, and this was true regardless of the exclusion of the transduced pancreatic cancer cell lines. This suggests that even as KRAS may influence some key changes to the epigenome, DNA methylation patterns observed are more influenced by cell type.
Unsupervised hierachical clustering shows cell line specific differential CpG methylation associated with KRAS suppression in pancreatic cancer cells. We have previously shown that silencing KRAS caused distinct molecular changes in pancreatic cancer cell lines 29 . Silencing of KRAS may therefore also lead to differential DNA methylation. To test this, we performed RNA-seq and genome-wide DNA methylation analysis using Illumina's Infinium arrays to determine the effect of silencing of KRAS in the 11 KRAS-mutant and -dependent pancreatic cancer cells. Briefly, cells were harvested for RNA and genomic DNA 4 to 7 days following infection with lentivirus shRNA targeting KRAS. Despite being KRAS-dependent, KRAS knockdown was not sufficient to cause dramatic cell death in pancreatic cell lines. This has been observed previously, and these cells lines were shown to be able to activate compensatory pathways in response to KRAS suppression 29 . Reduced KRAS mRNA levels were observed in KRAS-depleted cells relative to NS controls as determined by RNA sequencing (Fig. 2A). We then performed GSEA to compare the KRAS-depleted cells to the NS controls, and found a reduction in KRAS signaling ( Supplementary Fig. S3), as evident from decreased enrichment in genes which are upregulated by KRAS (HALLMARK_KRAS_UP) and increase enrichment in genes downregulated by KRAS (HALLMARK_KRAS_DN). There was also a decrease in both PI3K/AKT and mTORC signaling, which are pathways downstream of KRAS.
Unsupervised hierarchical clustering of genome-wide DNA methylation data using the top 1000 most variable CpG probes revealed co-clustering of isogenic cell line pairs in that all 11 KRAS-depleted cell lines and their isogenic controls appear more similar to each other than any other cell line (Fig. 2B). There was also no clear separation based on specific KRAS mutations (G12D vs G12V). There were two distinct branches separating the 11 isogeneic pairs (Fig. 2B), representing a group with a lower degree of methylation than the other. To identify common methylation changes between the pancreatic cell lines, we performed heirachical clustering using the union of differentially methylated probes (Δβ values ≥0.2 or ≤−0.2) appearing in at least 3 out of the 11 cell line pairs (a total of 204 CpG probes) and observed co-clustering of Pa16C, Pa01C, and Panc 10.05 ( Supplementary Fig. S4A).
www.nature.com/scientificreports www.nature.com/scientificreports/ From this list of 204 CpG probes, which represents the most frequently differentially methylated (DM) probes, we compiled the top 10 DM hypo or hypermethylated CpGs per cell line into a list ( Supplementary Fig. S4B).
Next we examined the number of DM probes per cell line as a measure of the extent of DNA methylation response due to KRAS inhibition. The DNA methylation profiles of Pa16C, Pa01C, PANC-1, and Panc 10.05 cells showed the most robust response to KRAS suppression. These four responsive cell lines showed at least 5-fold more DM CpGs compared to the seven other pancreatic cell lines tested ( Fig. 2C and Table 1). Although Pa16C cells are derived from Panc 10.05 cells 33 , Pa16C cells had more than 4-fold the number of DM CpGs ( Fig. 2C and Table 1). The four responsive cell lines showed a significant number of DM CpGs located in the promoter region (200-1500 nt upstream of the transcription start site) of dozens of functionally important genes ( Table 2). The methylation changes associated with KRAS suppression appeared to be cell line specific and were not generalizable within pancreatic cell lines. Although the methylation patterns in the NS shRNA-treated control cells were similar (Fig. 1), each cell line responded differently to the depletion of KRAS. Furthermore, two distinct groups emerged from the pancreatic cell lines, with seven lines displaying significantly less differential methylation compared to the four responsive cell lines (Pa16C, Pa01C, PANC-1, and Panc 10.05 cells). Taken together these results suggest that depleting oncogenic KRAS expression is cell line specific but also stochastic in nature. It is possible Eleven human pancreatic cancer cell lines were transduced with non-silencing (NS) shRNA (black bar above). DNA methylation patterns in these pancreatic cells were compared to the DNA methylation in lung epithelial SALEB/SAKRAS cells and Infinium methylation data obtained from ENCODE (www.encodeproject.org). The β value for each probe is represented with a color scale as shown in the key. Values closer to 1 represent highly methylated CpGs, while values closer to zero represent least methylated CpGs.
www.nature.com/scientificreports www.nature.com/scientificreports/  www.nature.com/scientificreports www.nature.com/scientificreports/ that whether a cell's CpG methylation profile is responsive or refractory to KRAS suppression likely depends on its genetic background and other factors.
Inhibitor treatment shows limited role for ERK in differential CpG methylation of Pa16C pancreatic cancer cells. Next, we investigated whether the methylation changes associated with KRAS suppression are dependent on ERK signaling, a major downstream effector of KRAS. We used Pa16C cells, the cell line with the greatest number of DM CpGs associated upon KRAS suppression (Fig. 2C), to test the effects of ERK inhibition on DNA methylation. Pa16C cells were treated with the ERK1/2-selective inhibitor, (ERKi, SCH772984) 34 (Supplementary Fig. S5) and the cells were harvested for protein and genomic DNA 3 and 7 days after treatment. The dose response of SCH772984 on Pa16C cell growth was determined ( Supplementary  Fig. S5A). Based on this, Pa16C cells were treated with 0.25 μM (3.6 on log scale), which resulted in the highest inhibition of cell growth. ERKi treatment led to growth arrest as evidenced by the lower cell confluency compared to DMSO control ( Supplementary Fig. S5B, Right) and also reduced total ERK protein and phosphorylated ERK protein as measured by western blot (Supplementary Fig. S5B, Left). Three and 7 days of ERKi treatment resulted in 29 and 37 DM probes, respectively. Only 1 CpG probe cg18988094 was hypomethylated in both 3-and 7-day ERKi-treated samples. This DM CpG is found near the gene STIP1, which has been reported to activate ERK signaling ( Supplementary Fig. S5C). We compared DM CpG profiles of the ERKi-treated Pa16C cells to the Pa16C cells transduced with shKRAS. However, there were no overlapping DNA methylation changes between the ERKi-treated and the KRAS shRNA-transduced Pa16C cells despite the similar effects on cell growth observed in both conditions ( Supplementary Fig. S5A, Right) 15 . Furthermore, KRAS shRNA induced >8000 DNA methylation changes compared to <40 DM CpGs after pharmacological ERK inhibition. These observations suggest that targeted ERK inhibition leads to Pa16C cell growth arrest similar to the growth arrest observed in KRAS shRNA transduced Pa16C cells. However, ERK does not appear to be consequential to the thousands of KRAS-associated DM CpGs present in the KRAS shRNA transduced Pa16C cells, at least after 7 days and suggests that KRAS suppression leads to sustained DM changes not affected by inhibition of downstream targets like ERK. However, it remains possible that ERK could still be responsible for KRAS-associated methylation changes that occur over a longer time frame. Fig. S4), we attempted to isolate biological processes associated with KRAS knockdown that are common between KRAS-depleted cell lines. First, we grouped the KRAS-depleted cell lines into "responsive" cells (Pa16C, Pa01C, PANC-1, and Panc 10.05 cells) and "refractory" cells referring to the other seven pancreatic cell lines in our panel, based on the number of DM CpGs identified (Fig. 2C). To focus our analysis on genes with DM CpGs most likely to produce transcriptional effects, we isolated DM CpGs found within promoter regions, 200-1500 bases upstream of the transcription start site of a gene, and within 4 kb of a CpG island, including shores and shelves. We then kept only the gene promoters that had consistent CpG differential methylation, where all of the CpGs were either hypermethylated or hypomethylated. Genes encoding transcription factors, oncogenes, kinases, and growth factors showed differential DNA methylation at their promoters in KRAS-depleted cells ( Table 2). Gene ontology analysis was performed using lists of promoters from each KRAS-depleted cell line that were hypermethylated or hypomethylated. The top ≤ 20 overlapping biological processes were compiled in a heat map (Fig. 3). Hypermethylated promoters in the responsive cells were enriched for genes involved in development and differentiation ( Fig. 3A, bold); however, the number of hypermethylated promoters was significantly reduced in the refractory cells (Table 1), which limited the number of associated biological processes (Fig. 3B). The hypomethylated promoters in both the responsive cells and the refractory cells were enriched for genes involved in development and differentiation (Fig. 3A,B). Gene ontology analysis produced a significantly lower number of biological processes for the refractory cell lines compared to responsive cells due to the paucity of DM CpGs in the these cell lines (Fig. 3C). A total of 18 biological processes were found exclusively in the responsive lines with 6 of these related to development (Fig. 3D, Top, bold). In addition, our analysis showed 7 processes that were potentially affected by KRAS suppression in both responsive and refractory cell lines (Fig. 3D, Bottom). Together these results suggest that KRAS suppression leads to differential DNA methylation affecting genes involved in development and differentiation, especially in responsive cell lines, and corroborates previous gene ontology analyses of DM genes in HRAS-transformed fibroblasts, which also showed an enrichment for genes involved in development and differentiation 28 .

Gene ontology analysis of DM CpGs reveals enrichment of genes involved in development and differentiation associated with changes in KRAS expression. To focus our analysis on genes with
DM CpGs most likely to produce transcriptional effects in SAKRAS cells overexpressing KRAS G12V, we isolated promoter regions with consistently DM CpGs as previously described for the pancreatic cell lines (Fig. 3). Five hundred and forty-seven genes met these conditions, including 196 genes with hypermethylated promoters. Genes encoding transcription factors, oncogenes, kinases, and growth factors showed differential DNA methylation at their promoters in SAKRAS cells overexpressing KRAS G12V (Table 3). Gene ontology analysis using the list of 196 hypermethylated gene promoters, and 351 hypomethylated gene promoters in SAKRAS cells overexpressing KRAS G12V, showed an enrichment for genes involved in development and differentiation (Fig. 5), consistent with our previous mutant KRAS loss-of-function studies performed in pancreatic cancer cells (Fig. 3 and Supplementary Fig. S6D,E). Gene ontology analysis using the list of hypermethylated and hypomethylated gene promoters from both SAKRAS KRAS G12V expressing cells and Pa16C KRAS knockdown cells, showed the common enrichment for genes involved in differentiation and development ( Supplementary Fig. S6D,E). It is noteworthy that while mutant KRAS knockdown and overexpression ultimately results in DM CpGs of genes involved in similar biological processes, the specific number of genes and location of DM affected are distinct and unique to each cell line.  www.nature.com/scientificreports www.nature.com/scientificreports/ To directly assess the role of mutant KRAS in maintaining DNA methylation patterns in the isogenic lung cells, we identified differentially methylated (DM) CpGs from SAKRAS cells in which KRAS expression had been suppressed transiently with KRAS siRNA and compared these to cells transfected with control siRNA ( Supplementary Fig. S6A,B). We observed 86 DM CpGs in SAKRAS cells following siRNA-mediated KRAS knockdown ( Supplementary Fig. S6B). Interestingly, only two of these CpGs were also DM in the SAKRAS vs SALEB cell comparison. This included LRRC7 and the pluripotency transcription factor, NANOG, which were both hypomethylated in SAKRAS cells compared to SALEB cells, and then hypermethylated following KRAS depletion via siRNA ( Supplementary Fig. S6B, Left). We identified 10 probes that were differentially methylated in opposite directions when comparing the Pa16C cells in which KRAS had been depleted with shRNA to the SAKRAS cells ( Supplementary Fig. S6C). Taken together, the lists of DM genes affected by changes in KRAS expression while distinct between cell lines, showed an enrichment for genes involved in development and differentiaion. www.nature.com/scientificreports www.nature.com/scientificreports/

Discussion
The RAS small GTPase is the most commonly mutated oncoprotein in cancer 1 . RAS and its downstream effectors control key aspects of cancer development. However, until recently, attempts to directly target oncogenic KRAS have been unsuccessful. In addition to aberrant signaling, the expression of mutant KRAS is correlated with global differential DNA methylation 24,25 . Therefore, epigenetic changes associated with oncogenic KRAS expression could be an avenue where the survival of KRAS-dependent cancer cells may be vulnerable. Here we demonstrated that the that cell type was more impactful than mutant KRAS on DNA methylation. KRAS-mutant PDAC cell lines were also classified based on the responsiveness of their methylome to KRAS depletion. Furthermore, a number of studies suggest that the majority of differential DNA methylation associated with cancer may be stochastic in nature -contributing to low levels of overlap and high heterogeneity between cell lines, even when they share the same genetic background and/or origin 28,[36][37][38][39][40] . It is most likely due to this stochastic nature that we did not observe previously described methylation events driven by RAS, such as the silencing of proapoptotic FAS by HRAS in fibroblasts 24 , and the silencing of IRAK3 by mutant KRAS 21 . However, we did identify novel changes in genes related to development and differentiation after KRAS silencing, which was common to all our pancreatic cancer cell lines but was more pronounced in the KRAS-responsive lines. This suggests that while many DNA methylation changes could be stochastic in nature and simply "passenger" events, or a consequence of their cell   www.nature.com/scientificreports www.nature.com/scientificreports/ state and cell lineage, KRAS is likely still able to influence key changes to the epigenome that are ultimately crucial for the cancer phenotype. More studies are needed to determine whether stratification, such as by cancer subtype, will reveal more consistent changes in methylation patterns.
Another interesting observation is the variable number of DM CpGs associated with KRAS knockdown and/ or KRAS overexpression. KRAS remains crucially linked to cell proliferation through RAF-MEK-ERK mitogen activated protein kinase (MAPK) cascade, its main effector pathway, and inhibition of this pathway reliably leads to growth arrest 15 . However, we showed that ERK was not responsible for changes in the methylome, at least over the time frame observed. While it is possible that ERK could play a role in methylation changes over a longer period of time, the question remains, if not ERK, which KRAS effectors are leading to short term DNA methylation changes. In mouse lung adenocarcinoma cells, YAP1 was able to rescue KRAS depleted cells, suggesting a relevant mechanism to bypass loss of KRAS signaling 41 . In the same study, KRAS also induced PI3K expression, and yet, the subsequent suppression of KRAS has no effect on the upregulated AKT activation. PI3K has been shown to compensate for KRAS suppression in pancreatic cancer cells and regulate epigenetic modifiers including DNMTs 42 . Cells in which KRAS levels have been genetically reduced display sensitivity to PI3K inhibitors and dual PI3K and MEK inhibitors have been found to be more effective than blocking the individual pathways alone 43 . PI3K/AKT signaling has been shown to be an epigenetic regulator in multiple cancers by modulating the activity of DNA methyltransferase I (DNMT1) 44 . It is possible that persistent PI3K-AKT activation, even after KRAS suppression, may be able maintain the majority of methylation changes induced by mutant KRAS. This kind of sustained activity by effector pathways could maintain the methylation status of the majority of the changes initially induced by mutant KRAS expression, but were not reversed upon KRAS knockdown (Fig. 6). www.nature.com/scientificreports www.nature.com/scientificreports/ DM CpGs associated with KRAS overexpression in our study have been localized to the promoters of important tumor suppressors, oncogenes, transcription factors, and regulators of differentiation, with gene ontology analysis revealing an enrichment for differentially methylated genes involved in differentiation and development. The regulation of pluripotency and lineage-specific genes requires the integration of multiple signaling pathways, epigenetic modifiers, and transcription factors 45 . In response to KRAS suppression, KRAS-driven cells may rely on compensatory survival pathways such as the PI3K pathway. PI3K-AKT has been shown to affect the expression of differentiation and stemness genes. In our pancreatic cells, particularly KRAS-responsive cells such as Pa16C, we identified many differentially methylated genes associated with stemness following KRAS knockdown, suggesting that KRAS could be involved in inducing stemness in cancer cells through PI3K/AKT. This includes promoter hypomethylation upstream regulators of AKT signaling, such as FGF9 46,47 , and NRG3, a ligand that activates HER3, an EGFR member of receptor tyrosine kinase (RTK) signaling upstream of PI3K-AKT 48 . POU3F2 and OLIG2 were both hypomethylated -two out of the four genes, from a core set of neurodevelopmental transcription factors (POU3F2, SOX2, SALL2, and OLIG2) essential for GBM propagation. These transcription factors coordinately bind and activate regulatory elements sufficient to fully reprogram differentiated GBM cells into tumor propagating stem-like cells 49 . Another promoter hypomethylated upon KRAS knockdown is HOXA9, a major transcription factor that regulates stem cells during development. Aberrant expression of HOX genes occurs in various cancers, and HOXA9 transcriptomes are specifically associated with cancer stem cell features 50 . Hypomethylation was also found at the BMP3 promoter. BMPs are implicated in activation of signaling pathways that drive epithelial-mesenchymal transition (EMT), including WNT signaling, TGFB signaling and PI3K signaling, all important pathways in pancreatic cancer cells 41,51,52 . And finally, another promoter which appeared as hypomethylated was TWIST1, a canonical EMT transcription factor shown to promote cancer stem cell properties 53 . Overexpression of TWIST1 is reported to override Myc-induced apoptosis in tumor cells and along with the other changes, could be a compensatory response by the Pa16C KRAS-mutant pancreatic cells to survive KRAS suppression.
Together, our findings suggest that while oncogenic KRAS-associated DNA methylation changes may be stochastic in nature and superseded by cell type, the changes nevertheless converge on biological processes most notably involving pathways of development and differentiation. That ERK inhibition was not analogous to Activating KRAS mutations lead to persistant induction of effector pathways that drive the cancer phenotype including the differential DNA methylation of genes involved in development and differentiation. In some cell lines, effector pathways such as PI3K and others, are able to maintain their abberant activity independent of KRAS signaling. As a consequence of feed forward loops initiated by mutant KRAS, kinome reprogramming, or the establishment of stable epigenetic patterns, the majority of DNA methylation changes associated with mutant KRAS activity remains refractory to KRAS suppression. However, independent of the changes in DNA methylation, KRAS knockdown and ERK inhibition still both lead to growth arrest in KRAS driven cell lines. SCH772984, type I and type II ERK inhibitor.
www.nature.com/scientificreports www.nature.com/scientificreports/ KRAS suppression in Pa16C cells suggests that KRAS-mediated DNA methylation are sustained independent of ERK. Taken all together, KRAS-mediated DNA methylation are stochastic and independent of canonical downstream effector signaling. This may therefore represent a non-canonical mechanism for enhancing tumorigenic potential and possibly help explain the ineffectiveness of KRAS effector inhibition in the clinic. Exploring the KRAS-mediated methylation changes in these pathways may be a deserving direction toward identifying supplementary strategies to target KRAS-driven cancers.

Methods
Cell culture. PDAC cell lines were obtained from ATCC and were maintained in Dulbecco's Modified Eagle Medium supplemented with 10% fetal calf serum (FCS) (HPAC and PANC-1), in RPMI 1640 supplemented with 10% FCS (CFPAC-1, HPAF-II, Panc 10.05, and SW-1990). Low passage SALEB and SAKRAS cells were generous gifts from Dr. Scott H. Randell (UNC-Chapel Hill) and were grown as described previously 54 . The SALEB cells were generated by infecting small airway lung epithelial cells with an amphotropic retrovirus that transduces SV40 ER, which encodes both the LT and small t antigens, and a neomycin drug resistance marker. These cells were subsequently infected with a retrovirus vector that transduces the hTERT gene together with the hygromycin resistance marker. Expression of these genetic elements was sufficient to immortalize the SALEB* cells. Finally, SALEB* cells were infected with retrovirus that transduces (i) the puromycin resistance marker (SALEB) or (ii) mutant KRAS G12V oncogene together with the puromycin resistance marker (SAKRAS). All other cells were maintained in Dulbecco's Modified Eagle Medium (DMEM; Gibco) supplemented with 10% fetal bovine serum (FBS; EMD Millipore). Cell lines were used for no longer than six months before being replaced. Stable cell lines were generated by selection in 2 μg/ml puromycin.
Western Blot reagents. Cells were lysed in 1% NP-40 lysis buffer (phosphatase and protease inhibitors from Sigma-Aldrich added fresh). Protein extracts were quantified by Bradford assay (Bio-Rad Laboratories) and analyzed by SDS-PAGE. Blot analyses were done with phospho-specific antibodies to ERK1/2 (T202/Y204) and antibodies recognizing total ERK1/2 to control for total protein expression. Antibody to KRAS4B was obtained from Calbiochem. Antibody for β-actin was used to verify equivalent loading of total cellular protein. Antibodies were purchased from Cell Signaling Techonology unless otherwise stated.
Small molecule inhibitors. The ERK1/2-selective inhibitor SCH772984 was provided by A. Samatar siRNA and shRNA transfection reagents. The following human siRNA (siGenome SMARTpool) was purchased from Dharmacon as a pool of four annealed dsRNA oligonucleotides: KRAS (L-005069-00) and non-targeting control #3 (D-001210-03). Dharmafect transfection reagent 1 was used to transfect 20-40 nM siRNA according to manufacturer's instruction and cells were harvested 96 hours after transfection. The target sequence for the validated shRNA construct used to target KRAS was CAGTTGAGACCTTCTAATTGG. The lentivirus vector encoding shRNA targeting KRAS (TRCN0000010369) was provided by J. Settleman (Genentech). Target cells were transduced by combining viral particle-containing medium with complete media at a ratio of 1:4 in the presence of polybrene (8 μg/ml). Media were exchanged 8-10 h later and selection was initiated following 16 h incubation in complete media. Samples were collected 72-120 h after the initiation of selection. DNA methylation analysis. Global DNA methylation was evaluated using the Infinium HumanMethylation450 BeadChip Array using more than ~450,000 Infinium CpG probes (Illumina, San Diego, CA). 1 μg of each DNA sample underwent bisulfite conversion using the EZ DNA Methylation Kit (Zymo Research, Irvine, CA) according to the manufacturer's recommendation for the Illumina Infinium Assay. Bisulfite-treated DNA was then hybridized to arrays according to the manufacturer's protocol. We used GenomeStudio V2011.1 (Illumina) for methylation data assembly and acquisition. Methylation levels for each CpG residue are presented as β values, estimating the ratio of the methylated signal intensity over the sum of the methylated and unmethylated intensities at each locus. The average β value reports a methylation signal ranging from 0 to 1, representing completely unmethylated to completely methylated values, respectively. Methylation data was preprocessed using the DMRcate package 55 . Data preprocessing included background correction, probe scaling to balance Infinium I and II probes, quantile normalization, and logit transformation. A logit transformation converts otherwise heteroscedastic beta values (bounded by 0 and 1) to M values following a Gaussian distribution. Additionally, detection p-values>0.05 in 25% of samples, probes on X and Y chromosomes, and probes situated within 10 bp of putative SNPs were removed. Differential methylation analysis on logit-transformed values was performed to compare samples in IMA. Wilcox rank test was conducted between experimental and control samples and p-values were corrected by calculating the false discovery rate by the Benjamini-Hochberg method. Probes with adjusted p-values <0.05, and delta β values ≥0.2 or ≤ −0.2 to 4 significant figures are considered statistically significant and differentially methylated. The methylation data discussed in this publication have been deposited in NCBI's Gene Expression Omnibus and are accessible at https://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE119548. The ENCODE methylation data used in this publication were retrieved from the ENCODE Data Coordination Center and are accessible at https://www.encodeproject. org/ucsc-browser-composites/ENCSR037HRJ. RNA sequencing and analysis. RNA sequencing was performed as described in Bryant et al. 29 . Briefly, a panel of human PDAC cell lines was infected with lentiviral vectors encoding shRNA targeting KRAS or a scrambled control construct for 8-10 h, and then selected for 48-96 h (depending on cell line). Following selection, cells were washed twice in ice cold phosphate-buffered saline (PBS), scraped in ice cold PBS, collected www.nature.com/scientificreports www.nature.com/scientificreports/ by centrifugation, and flash frozen. Total RNA (50 ng) for the pancreatic cell lines was used to generate whole transcriptome libraries for RNA sequencing using Illumina's TruSeq RNA Sample Prep. Poly-A mRNA selection was performed using oligo(dT) magnetic beads, and libraries were enriched using the TruSeq PCR Master Mix and primer cocktail. Amplified products were cleaned and quantified using the Agilent Bioanalyzer and Invitrogen Qubit. The clustered flowcell was sequenced on the Illumina HiSeq. 2500 for paired 100-bp reads using Illumina's TruSeq SBS Kit V3. Lane level fastq files were appended together if they were sequenced across multiple lanes. These fastq files were then aligned with STAR 2.3.1 to GRCh37.62 using ensembl.74.genes.gtf as GTF files. Transcript abundance was quantified and normalized using Salmon in the unit of transcripts per million (TPM). Clustering was performed using R heatmap.2 package with Euclidean Distance and McQuitty clustering method. Binary sequence alignment/map (BAM) files of RNA-seq data is available from the EMBL-EBI European Nucleotide Archive (ENA) database -http://www.ebi.ac.uk/ena/ with accession number PRJEB25797. The data are accessible at http://www.ebi.ac.uk/ena/data/view/PRJEB25797. The sample accession number is ERS2363485-ERS2363504.
Gene ontology analysis. The differentially methylated (DM) CpGs (i) in a promoter region (200-1500 bases upstream of the transcription start site of a gene) and (ii) within 4 kb of a CpG island (including CpGs at shores and shelves) are referred to as Promoter CpGs in this study. If a gene contains Promoter CpGs that did not all change in the same direction (all hypermethylated or all hypomethylated), that gene was excluded from analysis. Gene sets with hypermethylated or hypomethylated Promoter CpGs are loaded into Molecular Signature Database (MSigDB) 56 (http://www.broad.mit.edu/gsea/) and members of each gene set are categorized by gene families. The gene ontology analyses were generated using IPA (QIAGEN Inc., https://www.qiagenbioinformatics. com/products/ingenuity-pathway-analysis) 57 . The gene set of interest was uploaded into IPA (Ingenuity Systems, Redwood City, CA) and the Core Analysis workflow was run with default parameters. The Core Analysis provides an assessment of significantly altered pathways, molecular networks, and biological processes represented in the samples' gene list. The relative ranking order of biological processes were determined using a negative log10 scale of their p-values. The most enriched (top 20) biological processes with p-value <0.01 were picked. The gene sets used for analysis either contained hypermethylated Promoter CpGs only or hypomethylated Promoter CpGs only. Individual promoters with both hypermethylated and hypomethylated Promoter CpGs were excluded from gene set enrichment analysis.