Skip to main content

KrasG12D induces changes in chromatin territories that differentially impact early nuclear reprogramming in pancreatic cells

Abstract

Background

Pancreatic ductal adenocarcinoma initiation is most frequently caused by Kras mutations.

Results

Here, we apply biological, biochemical, and network biology methods to validate GEMM-derived cell models using inducible KrasG12D expression. We describe the time-dependent, chromatin remodeling program that impacts function during early oncogenic signaling. We find that the KrasG12D-induced transcriptional response is dominated by downregulated expression concordant with layers of epigenetic events. More open chromatin characterizes the ATAC-seq profile associated with a smaller group of upregulated genes and epigenetic marks. RRBS demonstrates that promoter hypermethylation does not account for the silencing of the extensive gene promoter network. Moreover, ChIP-Seq reveals that heterochromatin reorganization plays little role in this early transcriptional program. Notably, both gene activation and silencing primarily depend on the marking of genes with a combination of H3K27ac, H3K4me3, and H3K36me3. Indeed, integrated modeling of all these datasets shows that KrasG12D regulates its transcriptional program primarily through unique super-enhancers and enhancers, and marking specific gene promoters and bodies. We also report chromatin remodeling across genomic areas that, although not contributing directly to cis-gene transcription, are likely important for KrasG12D functions.

Conclusions

In summary, we report a comprehensive, time-dependent, and coordinated early epigenomic program for KrasG12D in pancreatic cells, which is mechanistically relevant to understanding chromatin remodeling events underlying transcriptional outcomes needed for the function of this oncogene.

Background

Epigenomic regulators are emerging as key critical downstream effectors of driver mutations during the process of cancer initiation, progression, metastasis, and tumor heterogeneity [1]. Pathogenic variants in KRAS act as driver mutations in many cancers, with pancreatic cancer being the most frequent. In this work, we use a mechanistically oriented, multi-omics approach based on the integrative analyses and modeling of several next-generation sequencing technologies to define epigenomic changes that underlie KrasG12D-mediated effects. Using cells derived from genetically engineered mouse models (GEMM) expressing this oncogene in an inducible manner, we systematically describe the impact of KrasG12D on both methylation changes and chromatin remodeling that account for transcriptional and non-transcriptional responses to signaling by this oncogene. First, using ChIP-seq, we demonstrate that while KrasG12D induces changes in heterochromatic and euchromatic histone marks, its effects on gene expression can be primarily explained by remodeling of enhancers and super-enhancers (H3K27ac and H3K4me1), promoters (H3K4me3), and gene body-associated functions (H3K36me3). Integration of all multi-omics data sets show that changes in heterochromatin marks (H3K9me3 and K27me3) have a minor impact on the transcriptome. Our data also provides insight into the repertoire of writer, reader, and eraser proteins, for which expression is consistent with deposition of their target marks. This data confirms and extends previous studies, by integrating more omics methodologies that are concomitantly performed in a tightly controlled and inducible GEMM-derived cell system [2]. Thus, this systematic investigation highlights the ability of KrasG12D to give rise to an epigenomic landscape that is conducive to pancreatic cell growth. Mechanistically, the finding that KrasG12D induced chromatin remodeling, primarily involving similar pathways for both gene activation and repression, helps to focus future experiments aimed at drugging these pathways to antagonize KrasG12D functions. Moreover, we describe chromatin remodeling events that do not appear to be directly associated with transcription, though may still account for KrasG12D effects. Thus, because the pathways studied here have increasingly available druggable options, this result bears biomedical relevance for chemoprevention and therapeutic studies.

Results

Cells derived from genetically engineered mouse models recapitulate early oncogenic Kras G12D signaling in vitro

First, we sought to better understand how KrasG12D coordinates over time, the reorganization of the genome. We used GEMM-derived cell models that carry a doxycycline-inducible KrasG12D transgene (iKras cell lines designated 4292, 9805, and 1012) [3, 4]. However, the ability of these cells to serve as models for studying molecular signaling and transcription coupling events is not known. Here, we show that induction of the oncogene by doxycycline increases their growth linearly after 48 hours (hrs), which was considered here as early events in KrasG12D function (Fig. 1a). We initiated our studies by using the 4292-cell line, although aspects of the data were also validated in the whole subset of cell lines. Total levels of histone marks were measured by western blot for enhancers (H3K27ac, H3K4me1), super-enhancers (H3K27ac), active promoters (H3K4me1 and H3K4me3), transcriptionally active gene bodies (H3K36me3) and silent heterochromatin (H3K27me3 and H3K9me3). We observed significant (p value < 0.05) increases in the levels of enhancers, super-enhancers, and active promoters through the H3K27ac (3.1 ± 0.7 fold) and H3K4me3 (2.7 ± 0.5 fold) histone marks, with the levels of H3K4me1 and H3K36me3 remaining relatively stable throughout the time course of KrasG12D induction (Fig. 1b, c). These changes were similarly observed in other iKras cell lines tested (9805, 1012 Additional file 1: Fig S1). Silent chromatin likewise increased in overall levels with H3K27me3 (2.9 ± 0.3 fold, p value < 0.01) and H3K9me3 (1.7 ± 0.1 fold, p value < 0.05; Fig. 1b, c and Additional file 1: Fig S1). We hypothesized that beyond the global changes of histone marks measured by western blot, there is a redistribution of histone marks in the genome that leads to epigenetic alterations with the ability to be targeted and antagonize the KrasG12D effect, a hypothesis that we explored using multi-omics methodologies.

Fig. 1
figure1

Induction of KrasG12D leads to changes in cell features. a iKras cells (4292, 9805, and 1012) were cultured without doxycycline (0 hr) or treated with doxycycline (48 hrs) to induce KrasG12D. Proliferation was measured by cell confluence in all four lines and expressed as a ratio of KrasG12D-expressing cells to control for each cell line. Cell confluency is an average of three separate experiments for each line. b Western blot analysis was performed in 4292 iKras cell line at 0, 12, 24, and 48 hrs time points. Cell lysates were probed with the indicated antibodies for KRASG12D and the histone marks. c Densitometry of western blots was performed on three separate experiments. * and ** indicates p value < 0.05 and 0.01 respectively. All data is expressed as mean ± SEM

Next, to assess signaling phosphorylation events induced by KrasG12D, we used a proteomic approach through phase-absorbed phospho-antibody (PAPA) arrays. Analysis of phospho-protein densitometric signals at 0, 12, 24, and 48 hrs after KrasG12D induction revealed a total of 66 differential protein phosphorylations (Fig. 2a, b). KrasG12D induction by western blot (Fig. 1c) shows moderate expression and levels continuing to increase out to 24 hrs. By mirroring these time points in the PAPA array, we may miss early signaling events, but we can observe phosphorylation changes that are initiated at 12 hrs and maintained by KrasG12D in these cells. There were 23 phosphorylation states in common for all three time points, whereas 2, 7, and 10 were exclusive to 12, 24, and 48 hrs, respectively (Fig. 2a). Comparing the signal ratios of protein phosphorylation across the time points demonstrates a complete reversal of the phosphorylation states between 0 and 48 hrs, with intermediate values following this trend at 12 and 24 hrs (Fig. 2b). We then categorized the protein phosphorylation status using two criteria, for each time point, based on their functional effect: (1) whether the magnitude of phosphorylation increased or decreased and (2) whether they are known to be activating or inhibiting modifications. Phosphorylation states were categorized as “activity down” if either an activating protein phosphorylation decreased in magnitude, or an inhibitory phosphorylation increased (Fig. 2c, d, Additional file 2: Table S1). Phosphorylation states were categorized as “activity up” if an activating phosphorylation increased in magnitude, or an inhibitory phosphorylation decreased. We performed pathway enrichment analysis of these two groups for each time point using the R package Rapid Integration of Term Annotation and Network resources (RITAN) [5] and the Molecular Signatures Database (MSigDB) hallmark gene set collection [6]. We confirmed that in these cell models, KrasG12D induced processes related to cell growth and survival via the PI3K-AKT-mTOR-NF-κB pathway (Fig. 2c, d). We also identified induction of signals for cell proliferation that include those which are G2/M checkpoint proteins and E2F targets (Fig. 2c, d), marking the arrival of the KrasG12D signal to the nucleus. Congruently, 24 hrs after KrasG12D induction, we noted activation of early response gene products, such as ELK-1 (pSer383, fold change = 1.6) and JUN (pSer73, fold change = 2.46), as well as activation of immediate early genes like MYC by downregulation of inhibitory phosphorylation (pSer62, fold change = 0.4, pThr58, fold change = 0.6) (Fig. 2d, bold letters and Additional file 2: Table S1) [7,8,9,10,11]. Additionally, downstream of MYC, phosphorylation of RB1 (pSer780, fold change = 1.8) (Fig. 2d, bold letters and Additional file 2: Table S1) releases E2F to mediate cell cycle progression and DNA synthesis [12, 13]. Another important pathway, recapitulated by our cell model, was that of PI3K/AKT signaling, including PDK1 (pSer241, fold change = 1.5) [14]. Cells also indicated an activated AKT1 (pSer473, fold change = 1.7; pThr208, fold change = 1.5) to promote cell survival by inhibitory phosphorylation events of BAD (pSer112, fold change = 1.6; pSer115 fold change, = 1.7) [15, 16] (Fig. 2d, bold letters and Additional file 2: Table S1). Next, evidence of AKT1 activation of MDM2 (pSer166, fold change = 2.6) was seen, which inhibits p53 [17] (Fig. 2d and Additional file 2: Table S1) as well as triggers pro-survival signals by NF-κB, through IKK-α (pThr23, fold change = 1.6) [18,19,20]. Interestingly, we found activation of several subunits of NF-κB (RELA, pThr254, fold change = 2.4; NFKB1, pSer893, fold change = 2.2; NFKB2, pSer869, fold change = 1.6) (Fig. 2d, bold letters and Additional file 2: Table S1). Lastly, we detected markers for oncogene-mediated replication stress, such as CHEK1 (pSer317, fold change = 2.1), CHEK2 (pThr68, fold change = 1.5), and BRCA1 (pSer1423, fold change = 1.5; pSer1524, fold change = 1.7) [21] (Fig. 2d, bold letters and Additional file 2: Table S1). As a control, western blots confirm the upregulation of phosphorylated AKT1, CDC25C, and NF-κB and downregulation of inhibitory MYC phosphorylation in the 4292, 9805, and 1012 KrasG12D cell lines at 48 hrs (Additional file 1: Fig S2). Thus, this data validates the usefulness of our cell models by showing that KrasG12D induces a proliferative phenotype, by triggering the canonical cytoplasmic signaling that reaches the nuclei to likely regulate, in a comprehensive and coordinated manner, a chromatin-mediated effect, which we investigate below.

Fig. 2
figure2

Analysis of cancer signaling phospho-antibody array containing 269 antibodies following oncogenic KrasG12D induction. a Venn diagram of differentially phosphorylated proteins (DPPs) at 12, 24, and 48 hrs with signal fold change ≥ |1.5|. b Signal ratios of DPPs normalized to the Z scale and plotted for 0, 12, 24, and 48 hrs. c Pathway enrichment analysis of DPPs for 12, 24, and 48 hrs using the Molecular Signatures Database (MSigDB) hallmark gene set collection. Color scale represents standardized –log10(FDR) values. d Protein-Protein interaction network obtained from DPP using the STRING database and interaction confidence score ≥ 999. Center blue circle: downregulated protein phosphorylation log2(fold change). Center yellow circle: upregulated protein phosphorylation log2(fold change). Some proteins have multiple phosphorylation sites and, hence, will have multiple colors in the center circle. Outer concentric circles represent activating (orange) or inhibiting (light blue) protein activity. The effect on activity by the combination of the mark’s effect and its level were interpreted. For instance, an inhibitory mark decreasing should lead to a more active protein. Some proteins have both inhibitory and activating phosphorylations and, hence, will have multiple colors in the outer concentric circle. Protein names highlighted in bold are mentioned in the text with their fold changes

Kras G12D changes histone chromatin accessibility and induces mobilization of histone marks to coding and non-coding areas of the genome

Using ATAC-Seq, we identified areas of chromatin that are remodeled in KrasG12D-induced cells. These analyses revealed 6123 regions that gain and 559 regions that lost accessibility (FDR ≤ 0.05; Fig. 3a, i). Most regions that change accessibility were in introns (33%), promoters (within 10 kb, 29%), and intergenic (29%) (Fig. 3a, iii). Hence, a characteristic of the KrasG12D response is predominantly a gain in accessible chromatin, affecting 2711 genes (2564 in their introns, 1563 in promoters, from 3595 to 1830 peaks, respectively), while only 323 genes became inaccessible (312 in their introns, 79 in promoters, from 346 and 82 peaks, respectively; Fig. 3a, ii). In fact, when considering the ± 10-kb region surrounding gene transcription start sites (TSS), the global chromatin accessibility in KrasG12D-induced cells increased the most in areas adjacent to the TSS (7.7% increase in the −1 kb promoter region and 10.4% in first intron) (Fig. 3b, i–ii). Thus, we investigated whether the genes within newly opened chromatin begin to account for the oncogenic signal induced by KrasG12D. Genes mapped to NF-κB and growth regulatory expression networks (Fig. 3c, i), encoding transcription factors that work downstream of these pathways, including AP1, ETS2, SRF, and MAZ [22,23,24] (Fig. 3c, ii). Figure 3d–e shows two examples of TSS accessibility changing after KrasG12D induction at single gene loci; the upregulated Etv4 gene has increased accessibility and the Cxcl15 locus shows decreased accessibility.

Fig. 3
figure3

Global changes in accessible chromatin following oncogenic KrasG12D induction. a (i) ATAC-seq was used to measure chromatin accessibilities. Each point represents a region of accessible chromatin comparing KrasG12D off (0 hr) versus on (24 hrs). a (ii) Venn diagrams comparing the number of gene promoters gained or lost in differentially accessible regions at 24 hrs. a (iii) Pie graphs showing the genomic distribution of the significantly differentially accessible regions. b (i) Average profile plot of normalized reads ± 10 kb around gene transcription start sites (TSSs). Y-axis represents read count per million (RPM) mapped reads. Orange and green shaded areas represent the standard error of the mean. b (ii) Heatmap of normalized reads around the TSS for each gene for 0 and 24 hrs. c (i) Pathway enrichment and (ii) transcription factor enrichment analysis of genes annotated to accessible chromatin for 0 and 24 hrs. Color scale represents standardized −log10(FDR) values. The tags following the transcription factors in c (ii) represent the TRANSFAC database nomenclature, which indicates either a different motif and/or the quality of the data that was used to make the motif. d, e Normalized ATAC-seq read coverage tracks for: d Etv4 (upregulated RNA-seq gene) and e Cxcl15 (downregulated RNA-seq gene)

Due to the fact that these open areas of chromatin may correspond to different epigenomic regulatory areas, such as promoters, enhancers, and super-enhancers, we next performed ChIP-seq on the six histone marks for both KrasG12D negative (0 hr) and KrasG12D-positive (24 hrs) cells (Fig. 1b, c). ChIP-Seq for H3K27ac revealed an overall gain in enrichment of differentially marked regions, which mapped to specific genes and were located at intronic, intergenic and promoter regions after 24 hrs of KrasG12D expression in cells (Fig. 4a, i–iii). We found a similar enrichment in differentially marked regions around gene TSSs for H3K4me3, which in our data set marked most promoters responsive to KrasG12D (Fig. 4b, i–iii). In contrast, KrasG12D had little effect on mobilizing the H3K4me1 mark, which was differentially bound at fewer gene loci as compared to previous marks and was present at intergenic, intronic, and promoter regions (Fig. 4c, i–iii). In response to KrasG12D, the H3K36me3 mark became equally divided between areas that show gain or loss of enrichment mapped to promoters, 3′UTRs, exons, and 5′UTRs (Fig. 4d, i–iii). Next, we studied marks that associate with inactive chromatin. We found that the facultative heterochromatin mark, H3K27me3, locates primarily outside of promoters in intergenic regions upon KrasG12D expression (Fig. 4e, iii), which is congruent with published cell biology data that show accumulation of non-promoter regions and heterochromatin at the nuclear periphery and lamina-associated domains (LADs) [25,26,27]. The other heterochromatic mark, H3K9me3, demonstrated a modest enrichment at differentially marked regions with KrasG12D at 24 hrs (Fig. 4f, i), mapping to 90 genes enriched and to 1 gene depleted (Fig. 4f, ii). These sites mostly mapped to intergenic, intronic, and exonic regions of the genome (Fig. 4f, iii). Promoters were poorly represented within H3K9me3 rich areas, consistent with a role for this mark in heterochromatic regions to form centromeres and membrane-to-heterochromatin attachment regions [26, 28]. KrasG12D activation also increased super-enhancers from 298 at 0 hr (Fig. 5a) to 415 after 24 hrs (Fig. 5b). When super-enhancers were mapped to genes, 87 were unique at 0 hr versus 184 at 24 hrs, with an additional 184 overlapping, although depending on the time point, signal intensity varied (Fig. 5c and Additional file 2: Table S2). We found that unique super-enhancers, as predicted using the ROSE algorithm [29, 30], before KrasG12D expression (0 hr) (Fig. 5d) were near genes regulating the mesenchymal phenotype and interferon responses, while super-enhancers after KrasG12D expression (24 hrs) were near genes encoding transcription factors that mediate oncogenic functions related to proliferation and cell cycle progression [31,32,33,34,35,36,37,38] (Fig. 5d). Super-enhancer formation also encompassed regulating pathways downstream of KrasG12D, including RAS, NF-κB, and G2M checkpoint-mediated gene expression networks [6, 39] (Fig. 5d). An independent ChIP assay performed on all three KrasG12D cell lines (4292, 9805 and 1012) in the promoter regions of Btc, Etv4, Cdkn1a, and Npm1 showed an increase in the deposition of H3K27ac and H3K4me3 marks and a decrease in the deposition of those same marks at Itgb5, Pdgfrb, Wnt10b, and Uba7. This was confirmation that KrasG12D induction leads to the same remodeling of activating chromatin marks as seen with ChIP-seq analysis (Additional file 1: Fig S3). These results indicate that super-enhancers make an important contribution to the KrasG12D response and its phenotypic transitions.

Fig. 4
figure4

Global changes in histone marks following oncogenic KrasG12D induction. a H3K27ac, b H3K4me3, c H3K4me1, d H3K36me3, e H3K27me3, f H3K9me3. a (i)–f (i): Each point represents a binding region comparing KrasG12D off (0 hr) versus on (24 hrs) condition. Loss or gain of enrichment set by an FDR ≤ 0.05. a (ii)–f (ii) Venn diagrams comparing the number of genes gained or lost in differentially bound regions at 24 hrs. a (iii)–f (iii) Pie graphs showing the genomic distribution of the significantly differentially bound regions

Fig. 5
figure5

Identification of super-enhancers following KrasG12D induction. a, b Super-enhancers identified at a 0 hr and b 24 hrs. c Venn diagram of super-enhancers annotated with genes comparing unique and shared super-enhancers at 0 and 24 hrs. d Pathway enrichment analysis of annotated super-enhancers for 0 and 24 hrs using the MSigDB hallmark gene set collection. Color scale −log10 (FDR). e Gene expression heatmap of all genes associated with super-enhancers at (i) only 0 hr (87 genes) and (ii) only 24 hrs (184 genes). f Gene expression heatmap of the transcription factors associated with super-enhancers at (i) only 0 hr and (ii) only 24 hrs. RPKM values were normalized to the Z scale. g Pathway enrichment analysis of transcription factors for 0 and 24 hrs using MSigDB hallmark gene set collection. Color scale −log10 (FDR)

We next investigated how these KrasG12D induced super-enhancer changes related to gene expression changes, proximally and via transcription factors. We plotted expression levels (reads per kilobase per million mapped, RPKM) of the 87 genes proximal to super-enhancers that are present only at 0 hr (Fig. 5e, i) and the 184 genes associated to those present only at 24 hrs (Fig. 5e, ii). Generally, genes associated with super-enhancers have higher expression compared to the same genes when the super-enhancers were not formed. This consistent change in gene expression supports a clear functional role for these changes in super-enhancer architecture, rapidly after oncogene induction. Next, we investigated how many of the super-enhancer-associated genes exclusive to each time point were transcription factors and had significant changes in gene expression, to better understand how super-enhancer assembly and disassembly could signal throughout the genome. At 0 hr (no KrasG12D), super-enhancers associated with 15 different transcription factors, out of which 13 were highly expressed and then downregulated by KrasG12D induction at 24 hrs (Fig. 5f, i). Notably, these transcription factors are all regulators of cell differentiation and xenobiotic metabolism [40,41,42,43,44]. There were 6 transcription factors associated with super-enhancers at only 24 hrs in KrasG12D and half of them had corresponding transcriptional upregulation at 24 hrs (Fig. 5f, ii). Therefore, gene expression was relatively concordant with super-enhancer regulation and further implicated transcription factors could be initiating genome-wide reorganization after super-enhancer remodeling. Ongoing studies that consider the impact of three-dimensional chromatin structure and regulation of super-enhancers on more distal genes will elucidate additional levels of regulation that may explain differences from the simple model that super-enhancers always upregulate cis-acting genes. Transcription factors among the differentially expressed genes at both time points function in various signaling pathways necessary for immune and oncogenic functions (Fig. 5g). Therefore, KrasG12D expression triggers an active remodeling of super-enhancers, which act in concert with the expression of transcription factors to contribute to the effects of this oncogene. Furthermore, combined ATAC-Seq and ChIP-Seq data demonstrates that KrasG12D induces a robust global remodeling of primarily activating chromatin marks for enhancers (H3K27ac), promoters (H3K4me3), and gene bodies (H36Kme3), with less of a contribution from H3K4me1. On the other hand, heterochromatin reorganization (H3K27me3, H3K9me3) is less associated with promoters, suggesting that these marks play an additional role apart from a direct influence on transcription, as revealed below by integrative analyses of all our omics datasets. These changes are a dominant feature associated to the KrasG12D phenotype in this GEMM-derived cell model.

Kras G12D induces rapid and global changes in DNA methylation, primarily at promoter regions and distinct repetitive areas of the genome

Due to the fact that many of the histone-based pathways work in coordination with 5-cytosine methylation of CpGs, we performed reduced representation bisulfite sequencing (RRBS) [45] (Additional file 2: Table S3). We found that KrasG12D deploys a rapid and robust methylation response with thousands of differentially methylated CpGs (DMCs) identified in KrasG12D-expressing cells (Fig. 6a, b). More than half of these DMCs were in CpG islands or flanking shores (Fig. 6c). Further annotation revealed that majority of these modifications occurred within promoters, followed by intergenic, intronic, and exonic regions (Fig. 6d). Pathway enrichment analysis of genes annotated to DMCs upon KrasG12D induction revealed that overall hypermethylation with concomitant gene silencing of promoters were primary involved in epithelial-to-mesenchymal transition (EMT), congruent with transcriptional silencing of genes in this pathway identified by RNA-seq experiments (Fig. 6e, f genes with yellow circles, Fig. 7c) [46,47,48,49].

Fig. 6
figure6

Identification of DNA methylation changes following oncogenic KrasG12D induction using RRBS. a Venn diagram of unique and overlapping differentially methylated cytosines (DMCs) for 12, 24, and 48 hrs. DMCs were classified based on p value < 0.01 and methylation difference > |10%|. b Methylation ratio of DMCs normalized to the Z scale and plotted for 0, 12, 24, and 48 hrs. Yellow: hypermethylated DMC. Blue: hypomethylated DMC. Black: no change. c Annotation of DMCs with CpG islands and shores. d Annotation of DMCs with genic features. e Pathway enrichment analysis of DMCs located within ± 3 kb of gene transcription start sites (TSSs) for 12, 24, and 48 hrs using the MSigDB hallmark gene set collection. Color scale: −log10(FDR). f Protein-protein interaction network obtained from DMCs for the epithelial to mesenchymal transition (EMT) pathway using the STRING database. Confidence score: 700. g Distribution of DMCs within various repetitive element types. h Methylation ratio of DMCs distributed within non-repeats and repetitive elements. Non-repeat n = 10,864, LINE n = 1411, LTR n = 791, Simple repeat n = 465, SINE n = 347. Wilcoxon signed rank test was used to test for differences between 0 and 24 hrs for each category. * significant at p value < 0.05

Fig. 7
figure7

RNA-seq analysis following oncogenic KrasG12D induction. a Venn diagram of differentially expressed genes (DEGs) at 12, 24, and 48 hrs. DEGs were classified based on expression fold change ≥ |2| and FDR ≤ 0.05. Total number of DEGs across time points: 2091. b RPKM expression levels of DEGs normalized to the Z scale and plotted for 0, 12, 24, and 48 hrs. Yellow: positive change. Blue: negative change. Black: no change. c Pathway enrichment analysis of DEGs for 12, 24, and 48 hrs using RITAN and MSigDB hallmark gene set collection. Color scale − log10(FDR). Gene expression networks obtained from upstream regulatory analysis at 24 hrs in Ingenuity Pathway Analysis (IPA) software for downregulated d TGFB1, and upregulated e MYC and f KRAS. The blue and green gene shapes and lines indicate predicted or observed inhibition respectively, while the yellow line indicates the predicted relationship is inconsistent with gene expression. The orange and red gene shapes and lines indicate predicted activation or observed upregulation respectively. Gray lines indicate no predicted effect

We next analyzed methylation of repetitive elements such as long interspersed nuclear elements (LINE), short interspersed nuclear element (SINE), and long terminal repeats (LTR), which can be differentially methylated in cancers and linked to genomic instability [50]. The majority of the DMCs found in KrasG12D-expressing cells were in non-repetitive elements, although a significant number of LINE, LTR, simple repeats, and SINE were also represented (Fig. 6g and Additional file 1: Fig S4a). Changes in DMCs within non-repeats and repetitive elements revealed that LINE are hypermethylated in KrasG12D cells at 12 hrs (Additional file 1: Fig S4b), but hypomethylated, along with LTR and SINE elements, after 24 hrs (Fig. 6h). At 48 hrs, LINE and SINE elements continued to be hypomethylated, whereas hypermethylation marked simple repeats (Additional file 1: Fig S4c). DMCs in non-repetitive regions became hypermethylated in KrasG12D cells at 24 and 48 hrs (Fig. 6h and Additional file 1: Fig S4c). Thus, we conclude that KrasG12D causes an overall decrease in methylation of LINE, SINE, and LTR by elements, which has yet to be recognized in response to this oncogene. Combined, these analyses demonstrate that KrasG12D triggers a rapid methylation response that varies in genomic location, although quantitatively similar through time. Simultaneously, KrasG12D causes demethylation of repetitive elements that associate to pathways such as reactivation of endogenous retroelements and genomic instability [51]. Thus, these phenomena are an important reflection of the potential of KrasG12D to participate in different aspects of carcinogenesis, such as genomic instability, which may underlie the accumulation of genetic alterations in pancreatic tissues over time, increasing the likelihood for a second hit.

Transcriptional outcome of the combined remodeling of nuclear functions induced by Kras G12D

RNA-seq, performed at 0, 12, 24, and 48 hrs after doxycycline treatment, detected 2091 differentially expressed genes (DEG) (Fig. 7a and Additional file 2: Table S4). A total of 646 DEGs were identified at all post-induction time points. However, distinct subsets only changed at specific times, with 248, 245, and 208 DEGs at 12, 24, and 48 hrs, respectively (Fig. 7a). Hierarchical clustering defined four expression patterns (Fig. 7b, color bars). Two patterns (labeled red and green) were primarily represented by transcripts that were downregulated in KrasG12D cells at 12 hrs and remained low through 24 hrs (1395 genes). At 48 hrs, 421 genes in the red pattern were upregulated, separating from the green pattern. The third pattern (labeled pink) contained 241 genes minimally expressed at 0 hr, but gradually increased to their maximum expression at 48 hrs. Finally, the 455 genes in the fourth pattern (labeled blue) increased in expression at 12 hrs and returned to near-basal levels by 48 hrs. Thus, in contrast with changes to chromatin, where we found a gain of enrichment for marks involved in transcriptional activation, nearly two-thirds of the DEGs exhibited robust downregulation. Pathway enrichment analyses using RITAN indicated that, at all three time points, downregulated genes, which included transcripts such as Snai2, Acta2, and multiple collagen genes, participate in EMT (Fig. 7c) [47]. Another set of downregulated genes were represented by the hallmark interferon alpha and gamma responses (Fig. 7c). Other algorithms such as IPA and cluster Profiler show that these pathways ranked among the top 5 in KrasG12D cells at 12 and 24 hrs (Additional file 1: Fig S5 and S6). Genes exclusively upregulated by KrasG12D at 12 and 24 hrs were enriched for proliferation and survival pathways, such as Myc targets and Mtorc1 signaling, both involved in protein synthesis and growth [52,53,54,55,56,57] (Fig. 7c and Additional file 1: Fig S5 and S6). Using upstream regulatory analysis (URA), we found Tgfb1 as a key node for regulating a mesenchymal network (Fig. 7d and Additional file 1: Fig S7). URA also shows that KrasG12D induction maintained the proliferative and growth-related MYC and KRAS gene networks at 24 hrs (Fig. 7e, f and Additional file 1: Fig S7). KrasG12D induction increased expression of EGF-like pathways [58], which reinforce its oncogenic potential [59,60,61] (Additional file 2: Table S4). To confirm the gene expression patterns and consistency of alterations in key pathways, transcriptomics were completed on the 1012 and 9805 cell lines (Additional file 2: Table S4). Filtered to match the set of 2091 genes altered by the 4292 cells, we observed that approximately 85% of the time the direction of the DEG was matched across all 3 cell lines (Additional file 1: Fig S8). Additionally, genes of the IFN-alpha, IFN-gamma, EMT, Myc, Mtorc1, and Kras up-/downregulation were similarly altered in these three KrasG12D-inducible cell lines. In summary, RNA-seq complements our cellular studies at a quantitative level by showing that KrasG12D induces the downregulation of a sizable percentage of gene networks with selective upregulation of those which promote oncogenicity.

Integrative analyses uncover single and combinatorial chromatin events underlying the early transcriptional response to Kras G12D

Our integrative analyses were initiated by filtering RNA-seq data and overlaying this data with peaks derived from ChIP-seq. We plotted H3K27ac, H4K4me3, H3K36me3, and H3K4me1 ChIP-seq signals for DEGs obtained from our RNA-seq data set in KrasG12D cells at 24 hrs (Fig. 8). We found an increase in H3K27ac and H3K4me3 at gene TSSs and H3K36me3 within gene bodies (Fig. 8a, i–iii). In fact, two-thirds of the genes expressed in KrasG12D cells displayed enrichment of these 3 marks at 24 hrs compared to 0 hr (Fig. 8a, v–vii). Concomitantly, nearly two-thirds of genes downregulated at 24 hrs in KrasG12D-expressing cells exhibited a decrease in these marks (Fig. 8b, i–iii, b, v–vii). In contrast, changes observed in other marks, such as H3K4me1 (Fig. 8a, iv, viii, b, iv, viii), H3K27me3, and H3K9me3 (Additional file 1: Fig S9), were not significantly enriched at promoter or enhancer regions and thus likely do not participate in transcriptional initiation. Moreover, we integrated information from ChIP-seq, ATAC-seq, and RRBS, using ChromHMM, to develop a model with 15 distinct chromatin states (States 1–15) as outputs, annotated as previously reported [62,63,64]. These analyses confirmed that active TSSs or flanking TSSs were marked by combinations of H3K4me3 with H3K27ac and H3K4me1 (Fig. 9a; States 1–4). Surprisingly, these regions combined only represent 1.3% of the entire genome. Of similar importance, these marks were absent in 74% of the genome (Fig. 9a; State 14). Thus, this data provides insights into the existence of an early KrasG12D-associated epigenomic landscape that accounts for the transcriptional outcome, which at a global and quantitative level facilitate the remodeling of euchromatic vs. repressive chromatin.

Fig. 8
figure8

Chromatin marks at gene transcription start sites and gene bodies correlate with up- and downregulated transcripts following KrasG12D induction. a (i–iv) Average profile plots of normalized H3K27Ac, H3K4Me3, H3K36Me3, and H3K4me1 reads around the transcription start site (TSS) and gene body for upregulated genes at 24 hrs in the RNA-seq data (446 genes). Orange and green shaded areas represent the standard error of the mean. a (v–viii) Red heatmaps show normalized reads around the TSS or gene body for each upregulated gene. b (i–iv) Average profile plots of normalized H3K27Ac, H3K4Me3, H3K36Me3, and H3K4me1 reads around gene TSS and the gene body for downregulated genes at 24 hrs in the RNA-seq data (1165 genes). Color scheme is same as a (i–iv). b (v–viii) Red heatmaps show normalized reads around the TSS or gene body for each downregulated gene

Fig. 9
figure9

Multi-omics data integration reveals single and combinatorial chromatin remodeling events that affect gene transcription and cell phenotype. a ChromHMM segmentation of the genome. TES: Transcription end site, S: States 1–15. Integration of ChromHMM with b RNA-seq, c DNA methylation, and d ATAC-seq. Chromatin segmentation of e two upregulated genes, Btc and Etv5 and f two downregulated genes, Acta2 and Pdgfrb illustrating the change in marks at 0 and 24 hrs. g Heatmap displaying log2 fold changes of genes in the RNA-seq, ATAC-seq, ChIP-seq, and methylation ratio of RRBS data set at 24 hrs. h Pathway enrichment analysis of genes represented in panel g using RITAN. Gain or loss of enrichment: RNA-seq genes enriched in (1) ATAC-seq (2), H3K27ac (3), H3K4me3 (4), H3K4me1 (5), H3K36me3 (6), H3K27me3 (7), H3K9me3, and (8) hypermethylated (gain), hypomethylated (loss) gene promoters respectively. i Gene networks generated from the genes represented in g. Center red or blue dots represent upregulated and downregulated genes in the RNA-seq data set respectively. Concentric circles around the dots represent changes in H3K27ac, H3K4me3, and H3K36me3 respectively (going inside out) with red representing gain of the mark and blue representing loss of the mark.  Enlarged figure with gene details provided as Additional file 1: Fig S10. j RPKM expression levels of differentially expressed histone readers, writers, and erasers normalized to the Z scale and plotted for 0, 12, 24, and 48 hrs

This pattern of histone-based chromatin remodeling correlated with expression of upregulated genes in RNA-seq, hypomethylation of DNA promoters by RRBS, and increased accessible chromatin regions by ATAC-seq (Fig. 9b–d; State 2). States 5 and 6, which represent transcription of genes at 5′ and 3′ ends, were marked with H3K4me3, H3K27ac, H3K4me1, and H3K36me3, constituting 0.26% of the genome (Fig. 9a; States 5–6). Nine percent of the genome was also enriched in H3K36me3 reflecting the ability of KrasG12D to induce an active program of post-initiation, transcriptional processing within gene bodies (Fig. 9a, b; States 7–8). Enhancer regions were the primary remodeling target of KrasG12D and comprised about 5% of the genome marked by both H3K4me1 and H3K27ac (Fig. 9a; States 9-13). Another important finding is that 10% of the genome was marked by H3K9me3 or H3K27me3 (Fig. 9a; States 13 and 15 respectively), which as suggested by ChromHMM, relocated to LADs (Fig. 9a; States 13–15) in response to KrasG12D induction. Mapping the ChromHMM states along a 2D linear gene representation depicts how chromatin states alternate along segments of DNA, although they likely represent contacts that these genes make with chromatin proteins which aggregate in distinct 3D domains [65]. For instance, two upregulated genes, Btc and Etv5, were part of the proliferative network observed in RNA-seq (Fig. 7f), which gained transcriptionally active chromatin and lost the quiescent state within the gene body in KrasG12D cells (Fig. 9e). Under the same conditions, Etv5 gained marks of an active TSS. Conversely, close examination of genes downregulated at 0 hr, representing networks regulated by Tgfb1 (Acta2 and Pdgfrb) (Fig. 7d), demonstrated that their gene bodies lost marks and acquired a quiescent chromatin state upon KrasG12D induction (Fig. 9f). We also integrated RNA-seq expression patterns with all epigenomic data sets by a hierarchical clustering-based method which illustrated that KrasG12D cells displayed a total of 859 DEGs with upregulated genes mapping to accessible chromatin regions (ATAC-seq) and enrichment for H3K27ac, H3K4me3, and H3K36me3 marks (Fig. 9g). Conversely, downregulated genes primarily located at accessible chromatin regions lost enrichment of H3K27ac, H3K4me3, and H3K36me3. However, both up- and downregulated genes demonstrated an overall increase in H3K27me3. Pathway enrichment analyses on the integrated epigenomics dataset found that genes within accessible ATAC-seq regions and triple-marked by H3K27ac, H3K4me3, and H3K36me3 participate in cell survival, growth, and proliferative pathways (Fig. 9h). Conversely, genes that lost H3K27ac, H3K4me3, and H3K36me3 but gained H3K27me3 and underwent promoter hypermethylation primarily belong to pathways that can aid cells to acquire a pro-oncogenic capacitation by KrasG12D (Fig. 9h). We also built interacting gene expression networks of loci modified by different epigenomic pathways, overlaying fold changes in RNA expression (center circle) with the changes in H3K27ac, H3K4me3, and H3K36me3 (Fig. 9i and Additional file 1: Fig S10). Most downregulated genes that control EMT as well as IFN-α and IFN-γ responses exhibited depletion of at least one of the activating histone marks (Fig. 9i and Additional file 1: S10). Conversely, upregulated genes that control ribosomal biogenesis through Myc were enriched in at least one of these marks (Fig. 9i and Additional file 1: Fig S10). Interestingly, analyses of the effects of KrasG12D on key epigenomic regulators that are known to impact the marks studied here demonstrated the downregulation of histone deacetylases (HDACs), enzymes which remove the H3K27ac mark, in order to dissemble enhancers and super-enhancers (Fig. 9j). Hence, these changes are congruent with a hierarchical and critical role of enhancers and super-enhancers and, potentially, some regulators (e.g., HDAC, histone acetyltransferases (HATs)) to initiate transcription in response to KrasG12D. In summary, this study integrates knowledge derived from multiple datasets, generated in a controlled and time-dependent manner, to account for many features shown to characterize the KrasG12D phenotype in GEMM-derived pancreatic cells.

Discussion

The current study extends our understanding of how the epigenome functions as an effector of early KrasG12D signaling. Using an extensive battery of state-of-the-art multi-omics methodologies, we have used a step-wise design to follow the impact of this oncogene on levels of chromatin marks and their genome localization for 6 different histone marks, as well as its impact on chromatin accessibility, DNA methylation, and the transcriptome. The analyses of results from each of these datasets were mapped initially to the entire genome (global remodeling), without taking into consideration their impact on transcription. Subsequently, the integration of all the data together provided information on how distinct changes in the epigenome associate or not with the transcriptional response to KrasG12D. The new knowledge derived from these experiments represents a robust characterization of these oncogene-driven events in a pancreatic cell model. Consequently, it is important to summarize and discuss these findings in light of the role that this oncogene has in preparing pancreatic cells to begin the transition toward abnormal cell growth.

The choice of cell model for the study was carefully considered since previous studies have demonstrated that pancreatic tumor cells are heterogenous [1]. Thus, we first validated the usefulness of GEMM-derived pancreatic cells carrying an inducible KrasG12D allele as an appropriate model for our study. These modifications were accompanied by changes in a phenotypic transition from a slow proliferative cell to a more epithelioid proliferative one. Using PAPA arrays, we defined that the extensive signaling cascade triggered by KrasG12D was congruent with its biological effect and complemented the analyses that support the use of these GEMM-derived cells as a model for studying the early epigenomic landscape of KrasG12D.

While several studies seeking to understand the role of KrasG12D-mediated effects in pancreatic cells using multi-omic approaches have emerged and provide useful data, our current investigation has the advantage of describing a controlled, time-dependent, coordinated, and comprehensive design, following the signal from the oncogene to the moment it enters into the nucleus to remodel the epigenome. This allows us to build a “landscape” for the early epigenomic response to KrasG12D, which should be a key resource for future mechanistic studies to determine the role of a myriad of writers, readers, and erasers of both DNA methylation and the histone code. First, we analyzed the data in a general manner, seeking to understand the total remodeling of the epigenome, independent of considering the impact on transcription. Subsequently, we mapped the key events that can account for the early oncogenic transcriptional response of KrasG12D. These results are important, since most work in this area has been directed to underscore transcriptionally linked chromatin-remodeling events. However, epigenomic changes without an apparent link to direct transcriptional control may still influence gene expression through changes in nuclear structure and dynamics. Further support for this idea is given by the fact that the transcriptional response to KrasG12D is dominated by gene silencing events, including for transcription factors which likely have trans-activity, while ATAC-Seq shows the largest effects on chromatin opening. However, heterochromatin marks (H3K9Me3 and H3K27Me3) do not appear as major direct contributors to the transcriptional outcome. Instead, we suggest, and is predicted by ChromHMM and published cell biology data, that heterochromatin marks primarily undergo relocation to the nuclear periphery and LADs [25,26,27], rather than with mRNA synthesis and processing. Therefore, KrasG12D mounts a transcriptional response primarily characterized by gene expression silencing, for which the underlying mechanism cannot be attributed to reduced chromatin accessibility, DNA methylation, nor heterochromatin formation on transcriptionally active genes.

Notably, we find that most of the transcriptional output in response to KrasG12D signaling can be accounted for by a hierarchical cascade of changes that begin with the formation of enhancers and super-enhancers to influence the function of gene promoters and bodies. Hence, our study both confirms and extends a previous report performed in GEMM using an elegant approach, which concluded that the metastatic potential of PDAC cells correlates with changes in enhancers reflected by a global H3K27ac enrichment [2]. Interestingly, although human and mice PDAC might differ in many aspects, work from our laboratory has also demonstrated a key role for enhancers and super-enhancers in enabling the acquisition of specific pancreatic cancer subtypes [64]. Hence, these studies expand our understanding of chromatin remodeling in response to KrasG12D signaling and may help to delineate potential evolutionarily conserved mechanisms for oncogenesis from mice-to-human. Interestingly, examination of the expression of writers, readers, and erasers of both DNA and histone-based pathways demonstrate that these changes primarily correlate with the downregulation of HDACs. This observation is relevant since these molecules can deacetylate the H3K27ac mark for enhancers and super-enhancers, which are necessary to trigger and maintain other epigenomic events.

We should also consider potential impact of the current study for future biomedical experiments in mice, as well as the search for novel therapeutic strategies. Indeed, many of the writers, readers, and erasers of both DNA methylation and the histone code are emerging as important mechanistic nodes in cancer development as well as promising new therapeutic targets [66]. Emerging in vivo studies using GEMM are finding that some nuclear regulators known to either deposit, bind, or reverse both DNA and/or histone marks impact KRASG12D initiation [2]. However, the field is just at an early stage and a deeper, exhaustive understanding of the role of the epigenome as an effector of oncogenes will remain a matter of intensive investigations [65].

Conclusion

In conclusion, outlining how all the pathways studied here work in concert and in a time-regulated manner to mediate the effects of the most common initiating mutation for PDAC is important for helping to interpret mechanistic experiments and inspire the exploration of future therapeutic directions. We demonstrate that KrasG12D has a quantitatively larger effect on transcriptional activation via remodeling euchromatin-associated histone pathways while its role in repression is rather passive, by decreasing their level and reorganization. We also describe, however, chromatin remodeling events which, although not directly associated with transcriptional output, are part of the KrasG12D epigenomic landscape likely necessary for achieving its full oncogenic potential. KrasG12D-induced chromatin remodeling in GEMM pancreatic cells, as captured using multi-omics, shows the early nuclear response and describes how changes in gene expression are brought about by this oncogene, bearing both mechanistic and potential medical relevance due to promising and emerging experimental therapeutics targeting these pathways. Thus, this new knowledge should be taken into consideration when planning future mechanistic studies directed toward antagonizing early effects of KrasG12D, a critical step for inhibiting cancer initiation.

Experimental procedures

Tissue culture and reagents

iKras cell lines 4292, 9805, and 1012 were maintained in complete media (RPMI, Gibco 11875 with 10% fetal bovine serum and 1 μg/mL doxycycline) to continually express constitutively active KrasG12D as described [4, 67]. Cell lines were provided by Dr. Marina Pasca di Magliano and have regularly tested negative for mycoplasma with the last test occurring in July 2020. For RNA, DNA, ATAC, and ChIP experiments, cells were grown for 48–72 hrs without doxycycline to inactivate expression of KrasG12D. To initiate the time course, cells were counted and split to 5E5 cells per 100 mm dish and allowed to attach overnight. The following morning, 0 hr (no KrasG12D expression) was collected per RNA/DNA/ChIP protocols and doxycycline (1 μg/mL) was added to media and allowed to incubate for 12, 24, and 48 hrs prior to harvest. The inducible KrasG12D expression system utilized here does not allow for the elimination of any confounding effects that doxycycline treatment alone may initiate on the transcriptome and epigenome of these cells. For proliferation assays, cells were plated at 1000 cells per well in a 96-well plate, doxycycline added to experimental groups and then incubated in the IncuCyte S3 with images captured every 6 hrs. Cell confluence was calculated using the IncuCyte software, with each condition normalized to 0 hr and fold changes calculated between doxycycline (KrasG12D expressing) and no doxycycline condition at each time point.

Western blot

iKras proteins were collected from cells at all time points after KrasG12D induction by lysis in RIPA buffer (20 mM Tris-HCl (pH 7.5) 150 mM NaCl, 1 mM Na2EDTA, 1 mM EGTA, 1% NP-40, 1% sodium deoxycholate), supplemented with phosphatase and protease inhibitors (Thermo Fisher Scientific). Equal amounts of protein were electrophoresed on 12% SDS-PAGE gels and transferred to nitrocellulose membranes (GE Healthcare). Membranes were blocked in either 5% milk or 3% BSA for 1 hr, and primary antibody incubations were performed overnight at 4 °C with rocking. Additional file 2: Table S5 contains primary antibodies and dilutions used. Anti-mouse or anti-rabbit secondary antibodies (1:5000, Millipore) were incubated on the membranes for 1 hr at room temperature with shaking, followed by detection of bands with chemiluminescence (Thermo Fisher Scientific). Uncropped versions of all western blots are available (Additional file 1: Fig S11, Fig S12, Fig S13). Quantification of bands in three independent biological experiments was completed using ImageJ and statistical significance determined by Student’s t tests in GraphPad Prism.

Phospho-antibody array

Cancer Signaling Phospho-Antibody Array from Full Moon BioSystems were used per manufacturer recommendations and sent back to the company for scanning and quantification. The protein phosphorylation state ratio was defined as the signal intensity of phospho site-specific antibody/signal intensity of site specify antibody. The ratio change between samples was defined as the treatment sample/control sample. RITAN [5] and the Molecular Signatures Database (MSigDB) hallmark gene set collection [6] were used to perform pathway enrichments.

Assay for transposase-accessible chromatin sequencing (ATAC-seq)

iKras 4292 pellets of 5E3 cells were collected at 0 and 24 hrs and resuspended in ATAC buffer and processed as described by Volk et al. [68]. Briefly, the genome was fragmented by tagmentation, DNA fragments purified, and sequencing and indexing primers added by PCR. Additional cycles of amplification, ranging from 10 to 14 cycles, was completed to minimally increase library content. Libraries were recovered with the MinElute PCR Cleanup kit (Qiagen, 28004) and size selected for fragment sizes of 100–500 bp with AMPure XP beads (Beckman Coulter, A63881). Sequencing was completed by the GSPMC at the Medical College of Wisconsin on the Illumina HiSeq-2500 with 126 bp paired end read length and an average of 125 million paired end reads per sample, completed in biological triplicate. The ATAC-seq pipeline from Encode consortium was used for adapter trimming, read alignment, and peak calling of libraries [69]. Briefly, paired end reads were mapped by Bowtie using the mm10 reference genome. ATAC-seq peaks were called using MACS2 software at p value ≤ 0.01. The R package DiffBind was used to identify differentially accessible regions with FDR ≤ 0.05. As defined in DiffBind, the coordinates and the size of these regions were defined by obtaining the narrowest region of overlapping peaks. The identified regions were annotated with genes if they were located within ± 10 kb of gene transcription start sites (TSSs) using the R package ChIPseeker [70]. R package RITAN and the Molecular Signatures Databases [6, 71] was used to perform pathway and transcription factor enrichment analysis. To identify global differences in accessible chromatin at each time point, ngsplot software was used to generate normalized tag density profiles around ± 10 kb of gene transcription start sites. ATAC-seq reads around specific genes were visualized using the integrative genomics viewer (IGV) [72]. Prior to visualization, reads were normalized (RPKM) using the deepTools software [73].

ChIP-seq

At time points of 0 and 24 hrs, DNA-protein interactions were crosslinked in 4292 iKras cells using 10% formaldehyde. The reaction was quenched with glycine and cells harvested in PBS containing protease inhibitors. Chromatin was sheared, immunoprecipitated, and processed as previously described in Lomberk et al. [64]. Supernatant was incubated with various histone mark antibodies (Additional file 2: Table S6) and protein G-agarose beads (Roche 11719416001) for 3 days with mixing at 4 °C. Following DNA purification, samples underwent real-time PCR at positive and negative loci and confirmed histone mark enrichment over 1% input. Sequencing was completed at the Mayo Clinic Medical Genomics Core on the Illumina HiSeq 2000 with 51 bp paired end and 25 million reads per sample, completed in biological duplicate. Data was analyzed using the HiChIP pipeline [74]. Briefly, paired end reads were mapped by BWA [75] using the mm10 reference genome and pairs with one or both ends uniquely mapped were retained. H3K4me3, H3K4me1, and H3K27ac peaks were called using the MACS2 software package [76] at FDR ≤ 1%. SICER [35] was used to identify enriched domains for H3K36me3, H3K27me3, and H3K9me3. The R package DiffBind [77, 78] was used to identify differentially bound regions (DBRs) with FDR ≤ 0.05. The DBRs were annotated with genes using the R package ChIPseeker [70]. The following cut-offs from the transcription start site was used for gene annotation: H3K27ac (± 10 kb), H3K4me3 (± 3 kb), H3K36me3 (− 3 kb,10 kb), H3K4me1 (± 10 kb), H3K27me3 (any distance from TSS), H3K9me3 (any distance from TSS). The ROSE software [29, 30] was used to identify stitched enhancers and to separate super-enhancers from typical enhancers based on K27ac bam files and peak files. The stitching distance used was 12,500 kb. Regions around TSS 2500 kb were not considered while calling super-enhancers. R package RITAN [79] and the MSigDB hallmark gene set collection [6] was used to perform pathway enrichment analysis of genes annotated to super-enhancers.

qPCR of histone marks at specific ChIP locations

At time points of 0 and 24 hrs, DNA-protein interactions were crosslinked in 4292, 9805, and 1012 iKras cells in triplicate and prepared for ChIP assay as described above. The immunoprecipitation was performed against H3K4me3 and H3K27ac marks (Additional file 2: Table S6). BED files generated by H3K4me3 and H3K27ac ChIP-seq performed on 4292 iKras were imported to IGV and enriched bound regions were selected for primer design (Additional file 2: Table S7). As sample reference, we utilized 1% of input for all ChIP performed, and to normalize the samples, we employed the following standard equation: 1% input = 2 ^ [(mean Ct input − 6.64) − mean Ct ChIP] × 100. The results are shown as fold change (sample/control).

Reduced representation bisulfite sequencing (RRBS)

At each time point (0, 12, 24, and 48 hrs), cells were washed once in PBS, detached with a scraper in 1 mL PBS, pelleted, and preserved at −20 °C. DNA was isolated from cell pellets using the QIAamp DNA Mini and Blood Mini Kit (Qiagen, 51304). A single biological condition was collected and provided to the Mayo Clinic Medical Genomics Core for bisulfite conversion, library preparation, and sequencing. Briefly, DNA (250 ng) was digested with Msp1 (New England Biolabs (NEB), R0106M) and purified using Qiaquick Nucleotide Removal Kit (Qiagen, 28004). End-repair A tailing was performed (NEB, M0212L) and TruSeq methylated indexed adaptors (Illumina, 15025064) were ligated with T4 DNA ligase (NEB, M0202L). Size selection was performed with Agencourt AMPure XP beads (Beckman Coulter, A63882). Bisulfite conversion was performed using EZ-DNA Methylation Kit (Zymo Research, D5001) as recommended by the manufacturer with the exception that incubation was performed using 55 cycles of 95 °C for 30 s and 50 °C for 15 min. Following bisulfite treatment, the DNA was purified as directed and amplified using Pfu Turbo C Hotstart DNA Polymerase (Agilent Technologies, 600414). Library quantification was performed using Qubit dsDNA HS Assay Kit (Life Technologies, Q32854) and the Bioanalyzer DNA 1000 Kit (Agilent Technologies, 5067-1504). The final libraries from RRBS were prepared for sequencing as per the manufacturer’s instructions in the Illumina cBot and HiSeq Paired end cluster kit v3. Samples were sequenced at 51 bp paired end reads using Illumina HiSeq 2000 with TruSeq SBS sequencing kit v3. Data was collected using HiSeq data collection v1.5.15.1 software, and bases were called using Illumina’s RTA v1.13.48. Raw data was further analyzed using SAAP-RRBS [80], a streamlined analysis and annotation pipeline for RRBS. Briefly, FASTQ files were trimmed to remove adaptor sequences and reads less than 15 bp were discarded. Trimmed FASTQ files were then aligned against the reference genome mm10 using BSMAP [81]. Methylation was reported along with custom CpG annotation. A minimum of five reads was required for inclusion of a cytosine in subsequent high-level analyses. Differential methylation of individual CpG loci was detected by performing a Fisher exact test with Methylkit [82] between a pair of samples at different time points (0, 12, 24, and 48 hrs) after selecting only the CpGs which have a data available across all the samples. The differentially methylated CpGs (DMCs) were selected according to the combination of p value (p value < 0.01) and absolute methylation difference > 10%. DMCs within ± 3 kb of gene TSS were used for pathway enrichment analysis using RITAN and the Molecular Signatures Database (MSigDB) hallmark gene set collection [6] for each time point. Gene networks were generated using RITAN and Cytoscape [83]. R package genomation [84] was used for DMC annotation with CpG islands and genic elements. Annotation of DMCs with repetitive elements was performed using repeat element coordinates (mm10) obtained from RepeatMasker [85] and BEDTools [86].

RNA-seq

At each time point (0, 12, 24, and 48 hrs for 4292 and 0, 24 hrs for 1012 and 9085), cells were washed once in PBS and harvested with RLT + βME as per RNeasy Mini Kit (Qiagen, 74106). RNA isolation included an on-column DNA digestion step and yielded > 1 μg per sample. For 4292 cells, three independent biological replicates were collected and provided to the Mayo Clinic Medical Genomics Core for sequencing. RNA libraries were prepared with the Illumina TruSeq RNA v2 kit and sequencing completed on the Illumina HiSeq-2000 with 101 bp paired end reads. For 1012 and 9085 cells, two independent biological replicates were collected and provided to the Genomic Science and Precision Medicine Center (GSPMC) at the Medical College of Wisconsin. Libraries were prepared with the Illumina TruSeq RNA stranded kit and sequencing completed on the Illumina NovaSeq6000 with 100 bp pair-end reads. Raw sequencing reads were processed through the Mayo RNA-seq bioinformatics workflow, MAPR-Seq v1.2.1.3 [87]. Raw and normalized (RPKM) counts for 23,398 genes and corresponding exons, expressed single nucleotide variants (SNVs) as well as gene fusions were obtained per sample. The R package from Bioconductor, edgeR v3.8.6 [88], was used for differential analysis comparing gene expression at all time points (12, 24, and 48 hrs) to 0 hr. Prior to the analysis, lowly expressed genes (raw counts less than 30 reads) were removed. Genes with false discovery rate (FDR) less than 5% and an absolute fold change ≥ 2 were considered to be significantly differentially expressed. The Ingenuity Pathway Analysis software (IPA, Qiagen) [89], RITAN and clusterProfiler [90] software was used to analyze canonical pathways, upstream regulators, diseases, and toxicity functions using the subset of genes that were significantly differentially expressed.

Data integration

RNA-ChIP integration

To identify local differences in occupancy of H3K27ac, H3K4me3, and H3K36me3 around differentially expressed genes (DEGs), we merged bam files from both replicates at 0 and 24 hrs and used ngsplot to generate normalized tag density profiles around gene TSSs or the gene body. The upregulated and downregulated DEGs were identified from 24 hrs in RNA-seq analysis.

Integration of ChromHMM with RNA-seq, ATAC-seq, and RRBS

We used the software ChromHMM which uses hidden Markov model to segment the genome into distinct states [91]. ChIP-sequencing replicates at each time point were merged before binarization step to learn models from the data. Fifteen states of segmentation were visualized for both 0 and 24 hrs time points. The 15-state segmentation was used to identify fold enrichment across various genomic features (note y-axis of fold enrichment plots such as CpG islands and refseq Exon). The coordinates for LADs were obtained from Peric-Hupkes et al. [92]. For integration of ChromHMM with RNA-seq, ATAC-seq, and RRBS, bed files containing coordinates for each region of interest were generated using the following methods: (1) RNA-seq: DEGs identified at 24 hrs were annotated with chromosome number, start and end position using the R Bioconductor package -biomaRt [93], (2) ATAC-seq: peak files for 0 and 24 hrs were sorted using BEDTools [86] and the three replicates for each time point merged using BEDOPS [94] to generate bed files, (3) RRBS: coordinates of the DMCs identified at 24 hrs (p value < 0.01 and methylation difference > │10%│) that were within ± 3 kb of gene TSSs were extracted into hyper- and hypomethylated bed files. ChromHMM output files containing the segmented genome for 0 and 24 hrs were used to view chromatin states of gene loci on the UCSC genome browser [95]. Heatmap of RNA-seq and combined epigenomic markers was generated using log2fold changes for RNA-seq, ATAC-seq, and histone marks, while methylation ratio was used for RRBS. Pathway enrichment analysis and gene networks were generated using RITAN [79], MSigDB hallmark gene set collection [6], and Cytoscape [83].

Availability of data and materials

All sequencing datasets generated for this manuscript are available in the ArrayExpress public database repository [96,97,98,99,100] under the following accession codes: E-MTAB-10909, ATAC-seq of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10901, ChIP-seq of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10900, RRBS of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10896, RNA-seq of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10897, RNA-seq of pancreatic cancer cells (1012 and 9805) derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele.

Abbreviations

GEMMs:

Genetically engineered cells and mouse models

RPM:

Read count per million

TES:

Transcription End Site

PAPA:

phase-absorbed phospho-antibody

RE:

Regulatory elements

ATAC-seq:

Transposase-accessible chromatin sequencing

ChIP-seq:

Chromatin immunoprecipitation followed by sequencing

RRBS:

Reduced representation bisulfite sequencing

HMM:

Hidden Markov modeling

LADs:

Lamina-associated domains

DB:

Differentially bound

TSS:

Transcription start site

RPKM:

Reads per kilobase million

DMCs:

Differentially methylated CpGs

EMT:

Epithelioid to mesenchymal transition

LINE:

Long interspersed nuclear elements

SINE:

Short interspersed nuclear element

LTR:

Long terminal repeats

DEGs:

Differentially expressed genes

RITAN:

Rapid integration of term annotation and network resources

MSigDB:

Molecular signatures database

IGV:

Integrative genomics viewer

IPA:

Ingenuity Pathway Analysis

URA:

Upstream regulatory analysis

HATs:

histone acetyltransferases

HDACs:

Histone deacetylases

References

  1. 1.

    Lomberk G, Dusetti N, Iovanna J, Urrutia R. Emerging epigenomic landscapes of pancreatic cancer in the era of precision medicine. Nat Commun. 2019;10(1):3875. https://doi.org/10.1038/s41467-019-11812-7.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  2. 2.

    Roe JS, Hwang CI, Somerville TDD, Milazzo JP, Lee EJ, Da Silva B, et al. Enhancer reprogramming promotes pancreatic cancer metastasis. Cell. 2017;170:875–888.e820.

    CAS  Article  Google Scholar 

  3. 3.

    Collins MA, Bednar F, Zhang Y, Brisset JC, Galbán S, Galbán CJ, et al. Oncogenic Kras is required for both the initiation and maintenance of pancreatic cancer in mice. J Clin Invest. 2012;122(2):639–53. https://doi.org/10.1172/JCI59227.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  4. 4.

    Collins MA, Brisset JC, Zhang Y, Bednar F, Pierre J, Heist KA, et al. Metastatic pancreatic cancer is dependent on oncogenic Kras in mice. PloS one. 2012;7(12):e49707. https://doi.org/10.1371/journal.pone.0049707.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  5. 5.

    Zimmermann M. RITAN: Rapid Integration of Term Annotation and Network resources. In: R package version 1.8.0, vol. 7; 2019. https://doi.org/10.7717/peerj.6994.

    Chapter  Google Scholar 

  6. 6.

    Liberzon A, Birger C, Thorvaldsdottir H, Ghandi M, Mesirov JP, Tamayo P. The Molecular Signatures Database (MSigDB) hallmark gene set collection. Cell systems. 2015;1(6):417–25. https://doi.org/10.1016/j.cels.2015.12.004.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  7. 7.

    Sommerlandt FMJ, Brockmann A, Rössler W, Spaethe J. Immediate early genes in social insects: a tool to identify brain regions involved in complex behaviors and molecular processes underlying neuroplasticity. Cell Mol Life Sci. 2019;76(4):637–51. https://doi.org/10.1007/s00018-018-2948-z.

    CAS  Article  PubMed  Google Scholar 

  8. 8.

    Sakamoto KM, Frank DA. CREB in the pathophysiology of cancer: implications for targeting transcription factors for cancer therapy. Clin Cancer Res. 2009;15(8):2583–7. https://doi.org/10.1158/1078-0432.CCR-08-1137.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  9. 9.

    Besnard A, Galan-Rodriguez B, Vanhoutte P, Caboche J. Elk-1 a transcription factor with multiple facets in the brain. Front Neurosci. 2011;5:35. https://doi.org/10.3389/fnins.2011.00035.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  10. 10.

    Babu GJ, Lalli JM, Sussman MA, Sadoshima JI, Periasamy M. Phosphorylation of elk-1 by MEK/ERK pathway is necessary for c-fos gene activation during cardiac myocyte hypertrophy. J Mol Cell Cardiol. 2000;32(8):1447–57. https://doi.org/10.1006/jmcc.2000.1185.

    CAS  Article  PubMed  Google Scholar 

  11. 11.

    Sears R, Nuckolls F, Haura E, Taya Y, Tamai K, Nevins JR. Multiple Ras-dependent phosphorylation pathways regulate Myc protein stability. Genes Dev. 2000;14(19):2501–14. https://doi.org/10.1101/gad.836800.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  12. 12.

    Johnson DG, Schneider-Broussard R. Role of E2F in cell cycle control and cancer. Front Biosci. 1998;3(4):d447–8. https://doi.org/10.2741/A291.

    CAS  Article  PubMed  Google Scholar 

  13. 13.

    Lee T, Yao G, Nevins J, You L. Sensing and Integration of Erk and PI3K Signals by Myc. PLOS Comput Biol. 2008;4(2):e1000013. https://doi.org/10.1371/journal.pcbi.1000013.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  14. 14.

    Casamayor A, Morrice NA, Alessi DR. Phosphorylation of Ser-241 is essential for the activity of 3-phosphoinositide-dependent protein kinase-1: identification of five sites of phosphorylation in vivo. Biochem J. 1999;342(Pt 2):287–92. https://doi.org/10.1042/bj3420287.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  15. 15.

    Persad S, Attwell S, Gray V, Mawji N, Deng JT, Leung D, et al. Regulation of protein kinase B/Akt-serine 473 phosphorylation by integrin-linked kinase: critical roles for kinase activity and amino acids arginine 211 and serine 343. J Biol Chem. 2001;276(29):27462–9. https://doi.org/10.1074/jbc.M102940200.

    CAS  Article  PubMed  Google Scholar 

  16. 16.

    Datta SR, Dudek H, Tao X, Masters S, Fu H, Gotoh Y, et al. Akt phosphorylation of BAD couples survival signals to the cell-intrinsic death machinery. Cell. 1997;91(2):231–41. https://doi.org/10.1016/S0092-8674(00)80405-5.

    CAS  Article  PubMed  Google Scholar 

  17. 17.

    Mayo LD, Donner DB. A phosphatidylinositol 3-kinase/Akt pathway promotes translocation of Mdm2 from the cytoplasm to the nucleus. Proc Natl Acad Sci U S A. 2001;98(20):11598–603. https://doi.org/10.1073/pnas.181181198.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  18. 18.

    Bai D, Ueno L, Vogt PK. Akt-mediated regulation of NFkappaB and the essentialness of NFkappaB for the oncogenicity of PI3K and Akt. Int J Cancer. 2009;125(12):2863–70. https://doi.org/10.1002/ijc.24748.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  19. 19.

    Israël A. The IKK complex, a central regulator of NF-kappaB activation. Cold Spring Harb Perspect Biol. 2010;2:a000158.

    Article  Google Scholar 

  20. 20.

    Adli M, Merkhofer E, Cogswell P, Baldwin AS. IKKalpha and IKKbeta each function to regulate NF-kappaB activation in the TNF-induced/canonical pathway. PloS one. 2010;5(2):e9428. https://doi.org/10.1371/journal.pone.0009428.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  21. 21.

    Kotsantis P, Petermann E, Boulton SJ. Mechanisms of oncogene-induced replication stress: jigsaw falling into place. Cancer Discov. 2018;8(5):537–55. https://doi.org/10.1158/2159-8290.CD-17-1461.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  22. 22.

    Liu F, Yang X, Geng M, Huang M. Targeting ERK, an Achilles’ heel of the MAPK pathway, in cancer therapy. Acta Pharm Sin B. 2018;8(4):552–62. https://doi.org/10.1016/j.apsb.2018.01.008.

    Article  PubMed  PubMed Central  Google Scholar 

  23. 23.

    Ray A, Ray BK. Induction of Ras by SAF-1/MAZ through a feed-forward loop promotes angiogenesis in breast cancer. Cancer Med. 2015;4(2):224–34. https://doi.org/10.1002/cam4.362.

    CAS  Article  PubMed  Google Scholar 

  24. 24.

    Cargnello M, Roux PP. Activation and function of the MAPKs and their substrates, the MAPK-activated protein kinases. Microbiol Mol Biol Rev. 2011;75(1):50–83. https://doi.org/10.1128/MMBR.00031-10.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  25. 25.

    Buchwalter A, Kaneshiro JM, Hetzer MW. Coaching from the sidelines: the nuclear periphery in genome regulation. Nat Rev Genet. 2019;20(1):39–50. https://doi.org/10.1038/s41576-018-0063-5.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  26. 26.

    van Steensel B, Belmont AS. Lamina-associated domains: links with chromosome architecture, heterochromatin, and gene repression. Cell. 2017;169(5):780–91. https://doi.org/10.1016/j.cell.2017.04.022.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  27. 27.

    Wang J, Kumar RM, Biggs VJ, Lee H, Chen Y, Kagey MH, et al. The Msx1 homeoprotein recruits polycomb to the nuclear periphery during development. Dev Cell. 2011;21(3):575–88. https://doi.org/10.1016/j.devcel.2011.07.003.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  28. 28.

    Allshire RC, Ekwall K. Epigenetic regulation of chromatin states in Schizosaccharomyces pombe. Cold Spring Harb Perspect Biol. 2015;7(7):a018770. https://doi.org/10.1101/cshperspect.a018770.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  29. 29.

    Whyte WA, Orlando DA, Hnisz D, Abraham BJ, Lin CY, Kagey MH, et al. Master transcription factors and mediator establish super-enhancers at key cell identity genes. Cell. 2013;153(2):307–19. https://doi.org/10.1016/j.cell.2013.03.035.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  30. 30.

    Loven J, Hoke HA, Lin CY, Lau A, Orlando DA, Vakoc CR, et al. Selective inhibition of tumor oncogenes by disruption of super-enhancers. Cell. 2013;153(2):320–34. https://doi.org/10.1016/j.cell.2013.03.036.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  31. 31.

    Yu KR, Park SB, Jung JW, Seo MS, Hong IS, Kim HS, et al. HMGA2 regulates the in vitro aging and proliferation of human umbilical cord blood-derived stromal cells through the mTOR/p70S6K signaling pathway. Stem Cell Res. 2013;10(2):156–65. https://doi.org/10.1016/j.scr.2012.11.002.

    CAS  Article  PubMed  Google Scholar 

  32. 32.

    Vleugel MM, Greijer AE, Bos R, van der Wall E, van Diest PJ. c-Jun activation is associated with proliferation and angiogenesis in invasive breast cancer. Hum Pathol. 2006;37(6):668–74. https://doi.org/10.1016/j.humpath.2006.01.022.

    CAS  Article  Google Scholar 

  33. 33.

    Ahlin C, Lundgren C, Embretsen-Varro E, Jirstrom K, Blomqvist C, Fjallskog M. High expression of cyclin D1 is associated to high proliferation rate and increased risk of mortality in women with ER-positive but not in ER-negative breast cancers. Breast Cancer Res Treat. 2017;164(3):667–78. https://doi.org/10.1007/s10549-017-4294-5.

    CAS  Article  PubMed Central  Google Scholar 

  34. 34.

    Hurlin PJ, Dezfouli S. Functions of myc:max in the control of cell proliferation and tumorigenesis. Int Rev Cytol. 2004;238:183–226. https://doi.org/10.1016/S0074-7696(04)38004-6.

    CAS  Article  PubMed  Google Scholar 

  35. 35.

    Zang C, Schones DE, Zeng C, Cui K, Zhao K, Peng W. A clustering approach for identification of enriched domains from histone modification ChIP-Seq data. Bioinformatics. 2009;25(15):1952–8. https://doi.org/10.1093/bioinformatics/btp340.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  36. 36.

    Wu S, Du Y, Beckford J, Alachkar H. Upregulation of the EMT marker vimentin is associated with poor clinical outcome in acute myeloid leukemia. J Transl Med. 2018;16(1):170. https://doi.org/10.1186/s12967-018-1539-y.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  37. 37.

    Dash S, Sarashetti PM, Rajashekar B, Chowdhury R, Mukherjee S. TGF-beta2-induced EMT is dampened by inhibition of autophagy and TNF-alpha treatment. Oncotarget. 2018;9(5):6433–49. https://doi.org/10.18632/oncotarget.23942.

    Article  PubMed  PubMed Central  Google Scholar 

  38. 38.

    Panda D, Gjinaj E, Bachu M, Squire E, Novatt H, Ozato K, et al. IRF1 maintains optimal constitutive expression of antiviral genes and regulates the early antiviral response. Front Immunol. 2019;10:1019. https://doi.org/10.3389/fimmu.2019.01019.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  39. 39.

    Cui X, Shen D, Kong C, Zhang Z, Zeng Y, Lin X, et al. NF-κB suppresses apoptosis and promotes bladder cancer cell proliferation by upregulating survivin expression in vitro and in vivo. Sci Rep. 2017;7(1):40723. https://doi.org/10.1038/srep40723.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  40. 40.

    Neavin DR, Liu D, Ray B, Weinshilboum RM. The role of the aryl hydrocarbon receptor (AHR) in immune and inflammatory diseases. Int J Mol Sci. 2018;19(12):3851. https://doi.org/10.3390/ijms19123851.

    CAS  Article  PubMed Central  Google Scholar 

  41. 41.

    Roy S, Guler R, Parihar SP, Schmeier S, Kaczkowski B, Nishimura H, et al. Batf2/Irf1 induces inflammatory responses in classically activated macrophages, lipopolysaccharides, and mycobacterial infection. J Immunol. 2015;194(12):6035–44. https://doi.org/10.4049/jimmunol.1402521.

    CAS  Article  PubMed  Google Scholar 

  42. 42.

    Date D, Das R, Narla G, Simon DI, Jain MK, Mahabeleshwar GH. Kruppel-like transcription factor 6 regulates inflammatory macrophage polarization. J Biol Chem. 2014;289(15):10318–29. https://doi.org/10.1074/jbc.M113.526749.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  43. 43.

    Doi A, Ishikawa K, Shibata N, Ito E, Fujimoto J, Yamamoto M, et al. Enhanced expression of retinoic acid receptor alpha (RARA) induces epithelial-to-mesenchymal transition and disruption of mammary acinar structures. Mol Oncol. 2015;9(2):355–64. https://doi.org/10.1016/j.molonc.2014.09.005.

    CAS  Article  PubMed  Google Scholar 

  44. 44.

    Fenizia C, Bottino C, Corbetta S, Fittipaldi R, Floris P, Gaudenzi G, et al. SMYD3 promotes the epithelial-mesenchymal transition in breast cancer. Nucleic Acids Res. 2019;47(3):1278–93. https://doi.org/10.1093/nar/gky1221.

  45. 45.

    Deaton AM, Bird A. CpG islands and the regulation of transcription. Genes Dev. 2011;25(10):1010–22. https://doi.org/10.1101/gad.2037511.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  46. 46.

    Kasashima H, Yashiro M, Kinoshita H, Fukuoka T, Morisaki T, Masuda G, et al. Lysyl oxidase is associated with the epithelial-mesenchymal transition of gastric cancer cells in hypoxia. Gastric Cancer. 2016;19(2):431–42. https://doi.org/10.1007/s10120-015-0510-3.

    CAS  Article  PubMed  Google Scholar 

  47. 47.

    Anastassiou D, Rumjantseva V, Cheng W, Huang J, Canoll PD, Yamashiro DJ, et al. Human cancer cells express Slug-based epithelial-mesenchymal transition gene expression signature obtained in vivo. BMC Cancer. 2011;11(1):529. https://doi.org/10.1186/1471-2407-11-529.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  48. 48.

    Murillo-Garzón V, Gorroño-Etxebarria I, Åkerfelt M, Puustinen MC, Sistonen L, Nees M, et al. Frizzled-8 integrates Wnt-11 and transforming growth factor-β signaling in prostate cancer. Nat Commun. 2018;9(1):1747. https://doi.org/10.1038/s41467-018-04042-w.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  49. 49.

    Morel A-P, Lièvre M, Thomas C, Hinkal G, Ansieau S, Puisieux A. Generation of breast cancer stem cells through epithelial-mesenchymal transition. PloS one. 2008;3(8):e2888. https://doi.org/10.1371/journal.pone.0002888.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  50. 50.

    Anwar SL, Wulaningsih W, Lehmann U. Transposable elements in human cancer: causes and consequences of deregulation. Int J Mol Sci. 2017;18(5):974. https://doi.org/10.3390/ijms18050974.

    CAS  Article  PubMed Central  Google Scholar 

  51. 51.

    Ishak CA, De Carvalho DD. Reactivation of endogenous retroelements in cancer development and therapy. Annu Rev Canc Biol. 2020;4(1):159–76. https://doi.org/10.1146/annurev-cancerbio-030419-033525.

    Article  Google Scholar 

  52. 52.

    van Riggelen J, Yetil A, Felsher DW. MYC as a regulator of ribosome biogenesis and protein synthesis. Nat Rev Cancer. 2010;10(4):301–9. https://doi.org/10.1038/nrc2819.

    CAS  Article  PubMed  Google Scholar 

  53. 53.

    Dai M-S, Lu H. Crosstalk between c-Myc and ribosome in ribosomal biogenesis and cancer. J Cell Biochem. 2008;105(3):670–7. https://doi.org/10.1002/jcb.21895.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  54. 54.

    Schlosser I, Hölzel M, Mürnseer M, Burtscher H, Weidle UH, Eick D. A role for c-Myc in the regulation of ribosomal RNA processing. Nucleic Acids Res. 2003;31(21):6148–56. https://doi.org/10.1093/nar/gkg794.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  55. 55.

    Ben-Sahra I, Manning BD. mTORC1 signaling and the metabolic control of cell growth. Curr Opin Cell Biol. 2017;45:72–82. https://doi.org/10.1016/j.ceb.2017.02.012.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  56. 56.

    Saxton RA, Sabatini DM. mTOR signaling in growth, metabolism, and disease. Cell. 2017;168(6):960–76. https://doi.org/10.1016/j.cell.2017.02.004.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  57. 57.

    Kim J, Guan K-L. mTOR as a central hub of nutrient signalling and cell growth. Nat Cell Biol. 2019;21(1):63–71. https://doi.org/10.1038/s41556-018-0205-1.

    CAS  Article  PubMed  Google Scholar 

  58. 58.

    Wieduwilt MJ, Moasser MM. The epidermal growth factor receptor family: biology driving targeted therapeutics. Cell Mol Life Sci. 2008;65(10):1566–84. https://doi.org/10.1007/s00018-008-7440-8.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  59. 59.

    Navas C, Hernández-Porras I, Schuhmacher AJ, Sibilia M, Guerra C, Barbacid M. EGF receptor signaling is essential for k-ras oncogene-driven pancreatic ductal adenocarcinoma. Cancer Cell. 2012;22(3):318–30. https://doi.org/10.1016/j.ccr.2012.08.001.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  60. 60.

    Schneeweis C, Wirth M, Saur D, Reichert M, Schneider G. Oncogenic KRAS and the EGFR loop in pancreatic carcinogenesis-a connection to licensing nodes. Small GTPases. 2017;9(6):457–64. https://doi.org/10.1080/21541248.2016.1262935.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  61. 61.

    Logsdon CD, Lu W. The significance of Ras activity in pancreatic cancer initiation. Int J Biol Sci. 2016;12(3):338–46. https://doi.org/10.7150/ijbs.15020.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  62. 62.

    Kundaje A, Meuleman W, Ernst J, Bilenky M, Yen A, Heravi-Moussavi A, et al. Integrative analysis of 111 reference human epigenomes. Nature. 2015;518(7539):317–30. https://doi.org/10.1038/nature14248.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  63. 63.

    Ernst J, Kellis M. Chromatin-state discovery and genome annotation with ChromHMM. Nat Protoc. 2017;12(12):2478–92. https://doi.org/10.1038/nprot.2017.124.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  64. 64.

    Lomberk G, Blum Y, Nicolle R, Nair A, Gaonkar KS, Marisa L, et al. Distinct epigenetic landscapes underlie the pathobiology of pancreatic cancer subtypes. Nat Commun. 2018;9(1):1978. https://doi.org/10.1038/s41467-018-04383-6.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  65. 65.

    Lomberk GA, Urrutia R. The triple-code model for pancreatic cancer: cross talk among genetics, epigenetics, and nuclear structure. Surg Clin North Am. 2015;95(5):935–52. https://doi.org/10.1016/j.suc.2015.05.011.

    Article  PubMed  PubMed Central  Google Scholar 

  66. 66.

    Lomberk GA, Iovanna J, Urrutia R. The promise of epigenomic therapeutics in pancreatic cancer. Epigenomics. 2016;8(6):831–42. https://doi.org/10.2217/epi-2015-0016.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  67. 67.

    Zhang Y, Morris JPT, Yan W, Schofield HK, Gurney A, Simeone DM, et al. Canonical wnt signaling is required for pancreatic carcinogenesis. Cancer Res. 2013;73(15):4909–22. https://doi.org/10.1158/0008-5472.CAN-12-4384.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  68. 68.

    Volk A, Liang K, Suraneni P, Li X, Zhao J, Bulic M, et al. A CHAF1B-dependent molecular switch in hematopoiesis and leukemia pathogenesis. Cancer Cell. 2018;34:707–723.e707.

    CAS  Article  Google Scholar 

  69. 69.

    Consortium, E. P. An integrated encyclopedia of DNA elements in the human genome. Nature. 2012;489(7414):57–74. https://doi.org/10.1038/nature11247.

    CAS  Article  Google Scholar 

  70. 70.

    Yu G, Wang LG, He QY. ChIPseeker: an R/Bioconductor package for ChIP peak annotation, comparison and visualization. Bioinformatics (Oxford, England). 2015;31:2382–3.

    CAS  Article  Google Scholar 

  71. 71.

    Liberzon A, Subramanian A, Pinchback R, Thorvaldsdóttir H, Tamayo P, Mesirov JP. Molecular signatures database (MSigDB) 3.0. Bioinformatics (Oxford, England). 2011;27:1739–40.

    CAS  Article  Google Scholar 

  72. 72.

    Robinson JT, Thorvaldsdóttir H, Winckler W, Guttman M, Lander ES, Getz G, et al. Integrative genomics viewer. Nat Biotechnol. 2011;29(1):24–6. https://doi.org/10.1038/nbt.1754.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  73. 73.

    Ramírez F, Dündar F, Diehl S, Grüning BA, Manke T. deepTools: a flexible platform for exploring deep-sequencing data. Nucleic Acids Res. 2014;42(W1):W187–91. https://doi.org/10.1093/nar/gku365.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  74. 74.

    Yan H, Evans J, Kalmbach M, Moore R, Middha S, Luban S, et al. HiChIP: a high-throughput pipeline for integrative analysis of ChIP-Seq data. BMC Bioinformatics. 2014;15(1):280. https://doi.org/10.1186/1471-2105-15-280.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  75. 75.

    Li H, Durbin R. Fast and accurate short read alignment with Burrows–Wheeler transform. Bioinformatics. 2009;25(14):1754–60. https://doi.org/10.1093/bioinformatics/btp324.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  76. 76.

    Zhang Y, Liu T, Meyer CA, Eeckhoute J, Johnson DS, Bernstein BE, et al. Model-based Analysis of ChIP-Seq (MACS). Genome Biol. 2008;9(9):R137. https://doi.org/10.1186/gb-2008-9-9-r137.

    CAS  Article  PubMed Central  Google Scholar 

  77. 77.

    Ross-Innes CS, Stark R, Teschendorff AE, Holmes KA, Ali HR, Dunning MJ, et al. Differential oestrogen receptor binding is associated with clinical outcome in breast cancer. Nature. 2012;481(7381):389–93. https://doi.org/10.1038/nature10730.

    CAS  Article  PubMed  Google Scholar 

  78. 78.

    Stark R, Brown G. (2011) DiffBind: differential binding analysis of ChIP-Seq peak data. http://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf.

    Google Scholar 

  79. 79.

    Zimmermann MT, Kabat B, Grill DE, Kennedy RB, Poland GA. RITAN: rapid integration of term annotation and network resources. Peer J. 2019;7:e6994. https://doi.org/10.7717/peerj.6994.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  80. 80.

    Sun Z, Baheti S, Middha S, Kanwar R, Zhang Y, Li X, et al. SAAP-RRBS: streamlined analysis and annotation pipeline for reduced representation bisulfite sequencing. Bioinformatics. 2012;28(16):2180–1. https://doi.org/10.1093/bioinformatics/bts337.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  81. 81.

    Xi Y, Li W. BSMAP: whole genome bisulfite sequence MAPping program. BMC bioinformatics. 2009;10(1):232. https://doi.org/10.1186/1471-2105-10-232.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  82. 82.

    Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13:R87.

    Article  Google Scholar 

  83. 83.

    Shannon P, Markiel A, Ozier O, Baliga NS, Wang JT, Ramage D, et al. Cytoscape: a software environment for integrated models of biomolecular interaction networks. Genome Res. 2003;13(11):2498–504. https://doi.org/10.1101/gr.1239303.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  84. 84.

    Akalin A, Franke V, Vlahoviček K, Mason CE, Schübeler D. genomation: a toolkit to summarize, annotate and visualize genomic intervals. Bioinformatics (Oxford, England). 2014;31:1127–9.

    Article  Google Scholar 

  85. 85.

    Tarailo-Graovac M, Chen N. Using RepeatMasker to identify repetitive elements in genomic sequences. Curr Protoc Bioinformatics. 2009;4:4.10.

    Google Scholar 

  86. 86.

    Quinlan AR, Hall IM. BEDTools: a flexible suite of utilities for comparing genomic features. Bioinformatics (Oxford, England). 2010;26:841–2.

    CAS  Article  Google Scholar 

  87. 87.

    Kalari KR, Nair AA, Bhavsar JD, O'Brien DR, Davila JI, Bockol MA, et al. MAP-RSeq: Mayo Analysis Pipeline for RNA sequencing. BMC Bioinformatics. 2014;15(1):224. https://doi.org/10.1186/1471-2105-15-224.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  88. 88.

    Robinson MD, McCarthy DJ, Smyth GK. edgeR: a Bioconductor package for differential expression analysis of digital gene expression data. Bioinformatics. 2010;26(1):139–40. https://doi.org/10.1093/bioinformatics/btp616.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  89. 89.

    Kramer A, Green J, Pollard J Jr, Tugendreich S. Causal analysis approaches in Ingenuity Pathway Analysis. Bioinformatics (Oxford, England). 2014;30:523–30.

    Article  Google Scholar 

  90. 90.

    Yu G, Wang L-G, Han Y, He Q-Y. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS. 2012;16(5):284–7. https://doi.org/10.1089/omi.2011.0118.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  91. 91.

    Ernst J, Kellis M. ChromHMM: automating chromatin-state discovery and characterization. Nat Methods. 2012;9(3):215–6. https://doi.org/10.1038/nmeth.1906.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  92. 92.

    Peric-Hupkes D, Meuleman W, Pagie L, Bruggeman SW, Solovei I, Brugman W, et al. Molecular maps of the reorganization of genome-nuclear lamina interactions during differentiation. Mol Cell. 2010;38(4):603–13. https://doi.org/10.1016/j.molcel.2010.03.016.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  93. 93.

    Durinck S, Spellman PT, Birney E, Huber W. Mapping identifiers for the integration of genomic datasets with the R/Bioconductor package biomaRt. Nat Protoc. 2009;4(8):1184–91. https://doi.org/10.1038/nprot.2009.97.

    CAS  Article  PubMed  PubMed Central  Google Scholar 

  94. 94.

    Neph S, Kuehn MS, Reynolds AP, Haugen E, Thurman RE, Johnson AK, et al. BEDOPS: high-performance genomic feature operations. Bioinformatics (Oxford, England). 2012;28:1919–20.

    CAS  Article  Google Scholar 

  95. 95.

    Karolchik D, Hinrichs AS, Kent WJ. The UCSC Genome Browser. Curr Protocols Bioinformatics. 2009;1:1.4.

    Google Scholar 

  96. 96.

    Mathison, A. J., Kerketta R., De Assuncao, T. M., Leverence, E., Zeighami, A., Urrutia, G., Stodola, T., Di Magliano, M. P., Iovanna, J. L., Zimmermann, M. T., Lomberk, G., and Urrutia, R. (2021) ATAC-seq of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10909. ArrayExpress. http://www.ebi.ac.uk/arrayexpress.

    Google Scholar 

  97. 97.

    Mathison, A. J., Kerketta R., De Assuncao, T. M., Leverence, E., Zeighami, A., Urrutia, G., Stodola, T., Di Magliano, M. P., Iovanna, J. L., Zimmermann, M. T., Lomberk, G., and Urrutia, R. (2021) RNA-seq of pancreatic cancer cells (1012 and 9805) derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10897. ArrayExpress. http://www.ebi.ac.uk/arrayexpress.

    Google Scholar 

  98. 98.

    Mathison, A. J., Kerketta R., De Assuncao, T. M., Leverence, E., Zeighami, A., Urrutia, G., Stodola, T., Di Magliano, M. P., Iovanna, J. L., Zimmermann, M. T., Lomberk, G., and Urrutia, R. (2021) RNA-seq of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10896. ArrayExpress. http://www.ebi.ac.uk/arrayexpress.

    Google Scholar 

  99. 99.

    Mathison, A. J., Kerketta R., De Assuncao, T. M., Leverence, E., Zeighami, A., Urrutia, G., Stodola, T., Di Magliano, M. P., Iovanna, J. L., Zimmermann, M. T., Lomberk, G., and Urrutia, R. (2021) RRBS of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10900. ArrayExpress. http://www.ebi.ac.uk/arrayexpress.

    Google Scholar 

  100. 100.

    Mathison, A. J., Kerketta R., De Assuncao, T. M., Leverence, E., Zeighami, A., Urrutia, G., Stodola, T., Di Magliano, M. P., Iovanna, J. L., Zimmermann, M. T., Lomberk, G., and Urrutia, R. (2021) ChIP-seq of pancreatic cancer cells derived from a genetically engineered mouse model, which harbor the inducible KrasG12D allele. E-MTAB-10901. ArrayExpress. http://www.ebi.ac.uk/arrayexpress.

    Google Scholar 

Download references

Review history

The review history is available as Additional file 3.

Peer review information

Tim Sands was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Funding

This work was supported by National Institutes of Health Grants R01 CA178627 (to G.L.) and R01 DK52913 (to R.U. and G.L.), the Advancing a Healthier Wisconsin Endowment (to each G.L. and R.U.), the We Care Fund for Medical Innovation and Research (to each G.L. and R.U.), and The Linda T. and John A. Mellowes Endowed Innovation and Discovery Fund (to R.U.).

Author information

Affiliations

Authors

Contributions

R.U., G.L., and A.M. designed the study. A.M., R.K., T.M., E.L., G.U., A.Z., M.Z., and T.S. acquired and analyzed the data. A.M., R.K., M.Z., T.M., G.L., and R.U. interpreted the data. A.M., R.K., T.M., M.M, J.I., M.Z., G.L., and R.U. contributed to the draft of the paper. A.M., T.M., M.Z., G.L., and R.U. revised the paper. The authors read and approved the final manuscript.

Authors’ information

Twitter handles: @MathisonAngie (Angela J. Mathison), @theysee (Romica Kerketta), @TimStodola (Timothy Stodola), @MTZimmermann (Michael T. Zimmermann), @GLomberk (Gwen Lomberk), @precision_ru (Raul Urrutia).

Corresponding authors

Correspondence to Gwen Lomberk or Raul Urrutia.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

Additional file 1: Fig S1.

Induction of KrasG12D leads to changes in global histone mark levels. Fig S2. Induction of KrasG12D leads to changes in phosphorylation events. Fig S3. Induction of KrasG12D leads to changes in the remodeling of activating chromatin. Fig S4. Annotation of DMCs with non-repeats and repetitive elements. Fig S5. Ingenuity Pathway Analysis (IPA) of RNA-seq genes. Fig S6. Gene enrichment analysis of RNA-seq genes. Fig S7. Upstream regulatory analysis of RNA-seq genes conducted in IPA. Fig S8. RNA-seq analysis following oncogenic KrasG12D induction in 1012 and 9085 cell lines. Fig S9. Chromatin marks at gene bodies for up and downregulated transcripts following KrasG12D induction. Fig S10. Networks generated from genes associated with H3K27ac, H3K4me3 and H3K36me3 marks, enlarged figure with gene details of Fig. 9i from the main manuscript. Fig S11. Uncropped version of western blots present in Fig. 1b of the main manuscript. Fig S12. Uncropped version of western blots present in Supplementary Figure 1. Fig S13. Uncropped version of western blots present in Supplementary Figure 2.

Additional file 2: Table S1.

Signaling Phospho Antibody Array. Table S2. Super-enhancer Signal Intensities. Table S3. Statistics based on per sample analysis for RRBS. Table S4. RNA-seq fold changes. Table S5. Antibodies used for western blot. Table S6. Antibodies used for ChIP-seq. Table S7. Primers used for ChIP-seq.

Additional file 3.

Review history.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and Permissions

About this article

Verify currency and authenticity via CrossMark

Cite this article

Mathison, A.J., Kerketta, R., de Assuncao, T.M. et al. KrasG12D induces changes in chromatin territories that differentially impact early nuclear reprogramming in pancreatic cells. Genome Biol 22, 289 (2021). https://doi.org/10.1186/s13059-021-02498-6

Download citation

Keywords

  • Pancreatic cancer
  • KRAS
  • Epigenomics