Global Regulation of the Histone Mark H3K36me2 Underlies Epithelial Plasticity and Metastatic Progression

Cancer cells are capable of dramatic changes in phenotype and function, a phenomenon known as cellular plasticity (1). Epithelial-to-mesenchymal transition (EMT; i.e., epithelial plasticity) and its reverse process (MET) are among the most widely studied examples of cellular plasticity; these programs are utilized during embryonic development for tissue morphogenesis and are frequently reactivated during tumor progression (2). In particular, histologic and transcriptomic signatures relevant to EMT are frequently associated with metastasis, therapeutic resistance, and other factors portending poor prognosis and patient mortality (3).

The cellular hallmarks of EMT, which include repression of epithelial proteins such as E-cadherin (Ecad), Occludin, and tight junction proteins, and induction of mesenchymal proteins such as N-cadherin (Ncad), Vimentin, and Fibronectin, are classically induced by a small cohort of master EMT-related transcription factors (“EMT-TF”) including SNAI1/2, TWIST1/2, and ZEB1/2 (2). However, the relative importance of each of these transcription factors varies depending on experimental and clinical contexts, leading to seemingly conflicting reports on the roles they play in EMT and cancer progression (4–6). Such findings raise the possibility that a more generalized mechanism or common molecular feature might exist that could unify the activity of these EMT-TFs.

Epigenetic modifiers and complexes can enact broad, genome-wide changes in transcription; such global epigenetic alterations have been reported to underlie other cell fate changes, particularly during normal embryonic development (7, 8). In tumor cells, several EMT-TFs have previously been reported to functionally cooperate with various epigenetic modifiers (9); however, such studies have focused on interactions at the promoters of preselected epithelial and mesenchymal genes (10–16). And although previous profiling efforts have identified genome-scale changes in DNA methylation (17) and certain histone marks (18) during experimentally induced EMT, it remains uncertain whether these epigenetic changes are necessary and/or sufficient to induce cancer cell plasticity or whether they represent correlated (secondary) events. In this study, we took an unbiased genetic approach to identify epigenome-wide mechanisms driving epithelial plasticity and metastatic progression, and directly interrogated the functional role of specific epigenetic modifications in cancer cell identity.

To better understand epigenetic mechanisms underlying epithelial plasticity, we established an experimental system relying on stochastic changes in cell state rather than genetic or chemical induction. We began by utilizing a panel of single-cell clones derived from the KPCY (LSL-Kras^G12D; Tp53^loxP/+; Pdx1-cre; LSL-Rosa^YFP/YFP) mouse model of pancreatic ductal adenocarcinoma (PDAC; ref. 19). Clones retained both epithelial and mesenchymal subpopulations, as defined by the presence or absence of surface Ecad (20) and other classic epithelial or mesenchymal markers (Fig. 1A and B; Supplementary Fig. S1A). To confirm that epithelial plasticity is associated with global differences in chromatin landscape, we performed Assay for Transposase-Accessible Chromatin using sequencing (ATAC-seq) on FACS-sorted Ecad⁺ versus Ecad⁻ cells and identified chromatin regions exhibiting differential accessibility. Ecad⁻ cells exhibited greater accessibility in genes associated with a mesenchymal phenotype and/or EMT (Fig. 1C; Supplementary Fig. S1B). In contrast, Ecad⁺ cells exhibited greater accessibility in genes associated with an epithelial phenotype, including Ecad itself (Fig. 1C). Thus, FACS sorting based on surface Ecad expression isolates cell populations with distinct chromatin profiles, reflecting transcriptional programs related to epithelial-to-mesenchymal state.

Next, we sought to ensure that the Ecad⁺ and Ecad⁻ fractions do not represent fixed, nonplastic subpopulations. To this end, we sorted pure populations of Ecad⁺ or Ecad⁻ cells. Over a 2-to-4-week period after initial sorting, all cultures exhibited phenotypic reversion (i.e., Ecad⁺ cells became Ecad⁻, and Ecad⁻ cells became Ecad⁺; Fig. 1D). Remarkably, all cultures returned to an equilibrium that recapitulated the Ecad⁺/Ecad⁻ ratio observed in the parental culture (Fig. 1D, dotted line). Thus, these cell clones exhibit dynamic plasticity (i.e., both epithelial and mesenchymal transitions) under standard culture conditions.

To identify epigenetic modifiers important for plasticity, we carried out a focused CRISPR/Cas9 screen using a single guide RNA (sgRNA) library targeting approximately 600 epigenetic modifiers and controls in a clonal cell line (3077 c4) stably expressing Cas9 (Fig. 1E). As expected, control sgRNAs targeting Zeb1, which we had identified as the most potent EMT-TF in this cell line (Supplementary Fig. S2A–S2D), were significantly enriched in the Ecad⁺ population and depleted in the Ecad⁻ fraction. Conversely, sgRNAs targeting Ecad exhibited the opposite pattern (Fig. 1F). Furthermore, sgRNAs targeting ribosomal proteins were depleted in both fractions, whereas the representation of control nontargeting sgRNAs did not change (Supplementary Fig. S3A and S3B).

Among the sgRNAs showing significant enrichment or depletion in the Ecad⁺ or Ecad⁻ populations (Supplementary Table S1), nuclear receptor binding SET domain protein 2 (Nsd2) emerged as the top gene for which all targeting sgRNAs were significantly enriched in the Ecad⁺ subpopulation (Fig. 1F; Supplementary Fig. S3C), second only to Zeb1. NSD2 is a histone methyltransferase that dimethylates histone H3 at lysine-36 (H3K36me2), a mark associated with actively transcribed genes (21), and acts as an oncogenic driver in multiple myeloma (22) and a regulator of invasion and metastasis in prostate cancer (23–25). Conversely, sgRNAstargeting lysine-specific demethylase 2A (Kdm2a) were significantly overrepresented in the Ecad⁻ subpopulation (Fig. 1F; Supplementary Fig. S3D). KDM2A is a demethylase that preferentially targets H3K36me2 (26–28). Both Nsd2 and Kdm2a have paralogs with similar catalytic activities, but none of these were enriched or depleted in our screen (Supplementary Fig. S3C–S3F), possibly due to low levels of expression (Supplementary Fig. S3G). These results suggest that NSD2 and KDM2A act nonredundantly on H3K36 as determinants of the mesenchymal or epithelial state, respectively.

To validate the role of NSD2 and KDM2A in epithelial plasticity, we performed gain- and loss-of-function studies in a panel of mouse and human PDAC cell lines. Nsd2 loss led to an increase in Ecad levels and decrease in H3K36me2 levels (Fig. 1G; Supplementary Fig. S4), whereas loss of Kdm2a resulted in reduced Ecad and increased H3K36me2 (Fig. 1H). Both of these phenotypes could be rescued by reexpression of the appropriate gene (Fig. 1I and J). Importantly, despite lower expression of the catalytically active mutant NSD2 (E1099K; ref. 29), this mutant potently reversed the effects of Nsd2 loss, suggesting that NSD2′s histone methyltransferase activity is critical for EMT. Notably, these manipulations had little effect on H3K36me3 levels, in line with previous observations (22, 30), implying that gain or loss of the H3K36me2 mark is regulated independently of effects on higher-order methylation. The opposing catalytic activities of NSD2 and KDM2A, together with their shared substrate specificity, thus implicate a mechanism wherein the writing and erasing of H3K36me2 regulates a cell's position along the epithelial-to-mesenchymal spectrum.

To understand the transcriptional changes regulated by NSD2 and KDM2A, we carried out RNA sequencing (RNA-seq) and gene set enrichment analyses (GSEA) in clonal wild-type (WT), sgNsd2, and sgKdm2a cells (Supplementary Tables S2 and S3). “EMT” was among the most significantly affected hallmark gene sets (Fig. 2A), and many classic markers of EMT were downregulated with Nsd2 loss and upregulated with Kdm2a loss (Supplementary Fig. S5A–S5D). Importantly, many other hallmark and curated gene sets enriched following Nsd2 loss were negatively enriched by Kdm2a loss, and vice versa (Supplementary Fig. S6A–S6C), suggesting that these enzymes regulate a shared transcriptomic program in a reciprocal manner.

In the clinical setting, the squamous subtype of PDAC describes a subset of tumors that are poorly differentiated and of particularly poor prognosis (31–33). Using gene- expression signatures that discriminate the squamous subtype from the well-to-moderately-differentiated progenitor subtype in patient samples (34), we found that Nsd2 loss results in negative enrichment of the squamous signature and positive enrichment of the progenitor signature, with a converse pattern observed in the setting of Kdm2a loss (Fig. 2B). These subtype changes are also reflected in the differentiation status of cultured cells and tumors arising from these cells (Fig. 2C). These results suggest that NSD2 and KDM2A exert opposing effects on cellular transcriptomes, leading to changes in tumor subtype and differentiation.

Because NSD2 appears to drive a PDAC subtype of particularly poor prognosis, we wished to evaluate its prognostic value across various cancer types. NSD2 and its paralogs are all able to dimethylate H3K36 and are each catalytically dominant in different tissue types (35–37). We thus used a gene signature representing the average expression of NSD1, NSD2, and NSD3 to more broadly capture the consequences of increased H3K36me2 deposition on clinical outcomes. Upon applying this gene signature to transcriptomic databases across various carcinomas (38, 39), we found that high expression of the signature predicted poor prognosis (Fig. 2D). In contrast, the expression of individual members of the NSD family had substantially less predictive value (data not shown). These results suggest that elevated expression of writers of H3K36 dimethylation is associated with poor clinical outcome, underscoring the importance of considering shared catalytic activities or molecular functions in relation to tumor behavior.

Tumor-initiating cells (TIC), a subpopulation of tumor cells capable of seeding and propagating new tumors, are commonly enriched for cells that are transcriptionally characterized by EMT signatures (40). Furthermore, induction of an EMT program has been shown to directly regulate the TIC phenotype (10). Accordingly, we found that forcing cells to be epithelial by loss of Nsd2 decreased their ability to form 3-D tumor spheres, whereas locking cells in the mesenchymal state by loss of Kdm2a increased their sphere-forming ability, suggesting that these enzymes reciprocally regulate tumor initiation capacity (Fig. 2E). Collectively, these findings demonstrate that regulators of the histone mark H3K36me2 drive attributes of tumor progression that are associated with poor clinical outcome.

Although epithelial-to-mesenchymal plasticity is thought to drive metastatic dissemination, recent studies have called this paradigm into question (5, 6, 41). To investigate the contributions of the epithelial and mesenchymal states to metastasis, we orthotopically transplanted tumor cells and found that Kdm2a-deficient tumors (locked in a mesenchymal state) released more circulating tumor cells (CTC) into the bloodstream compared with Cas9 control tumors (Fig. 3A), which is consistent with their increased migratory behavior in vitro (Supplementary Fig. S7A and S7B). Conversely, Nsd2-deficient tumors (locked in an epithelial state) released few CTCs into the bloodstream (Fig. 3A). These results support the notion that cells with a mesenchymal phenotype enter the bloodstream far more efficiently than cells with an epithelial phenotype.

Although Kdm2a-deficient tumor cells readily intravasated into circulation, the resulting metastatic lesions found in the lungs were significantly smaller than those arising in their WT counterparts (Fig. 3B).To further study the effect of epithelial-to-mesenchymal state on metastatic colonization, we injected Nsd2- and Kdm2a-deficient tumor cells directly into the tail veins of mice, thereby bypassing the first steps of the metastatic cascade. Nsd2-deficient cells were more efficient than control cells at forming large lung metastases (Fig. 3C), whereas Kdm2a-deficient cells remained poor colonizers (Fig. 3D). Importantly, these differences in metastatic competence were not due to differences in tumor growth or cell proliferation (Supplementary Fig. S7C and S7D). Taken together, these results demonstrate that epithelial and mesenchymal phenotypes each have their own distinct advantages during metastatic spread, supporting the concept that EMT and MET contribute at different points in the metastatic cascade and that plasticity is more important for effective metastasis than either the epithelial or mesenchymal state alone (42–46).

NSD2 and KDM2A act antagonistically to regulate H3K36me2, strongly indicating that this histone mark is a key determinant of epithelial identity. To test this hypothesis, we first performed histone mass spectrometry to look for global changes associated with epithelial plasticity. Remarkably, H3K36me2 was the only histone post-translational modification (PTM) to exhibit a significant difference (fold change >1.5, P_adj < 0.05) between Ecad⁻ and Ecad⁺ cells across multiple cell lines (Fig. 4A; Supplementary Fig. S8A and S8B). In this setting, the enrichment of H3K36me2 in mesenchymal cells was accompanied by higher expression of Nsd2 (Fig. 4B). Differences in global levels of H3K36me2 were readily detected by Western blots of sorted Ecad⁺ and Ecad⁻ cells (Supplementary Fig. S8C), and in cultured cells (Supplementary Fig. S8D) and PDAC tumors (Supplementary Fig. S8E) stained for H3K36me2. Thus, transition to a mesenchymal state is reproducibly associated with global accumulation of H3K36me2 in multiple in vitro and in vivo settings.

To determine whether differences in H3K36 methylation also occur in the setting of experimentally induced plasticity, we treated cells with TGFβ, a potent EMT stimulus, and again found that H3K36me2 was the only histone PTM to exhibit a global increase across cell lines (Fig. 4C). Upon TGFβ stimulation, different cell lines upregulated the expression of different H3K36me2 methyltransferases, each of which can contribute to the observed global increase in H3K36me2 (Fig. 4D). Interestingly, despite overt morphologic changes within 4 days of TGFβ treatment (Supplementary Fig. S9A), a buildup of H3K36me2 was not apparent until 15 days of treatment (Supplementary Fig. S9B–S9E), suggesting that global H3K36me2 accumulation may be a gradual phenomenon that helps to maintain the mesenchymal phenotype. To explore this possibility, we withdrew TGFβ after 4 or 15 days of treatment and assessed how long it took cells to revert to an epithelial state. Although the effects of short-term TGFβ treatment were readily reversible upon TGFβ withdrawal, cells exposed to TGFβ for 15 days remained mesenchymal for much longer (Fig. 4E). Eventually, when these cells finally reverted to their original epithelial state (Fig. 4E; Supplementary Fig. S10A), their H3K36me2 levels also began returning to baseline (Supplementary Fig. S10B), demonstrating that epigenetic plasticity accompanies phenotypic plasticity. Collectively, these results suggest that TGFβ-induced EMT involves an acute phase, likely mediated by the early activity of EMT-TFs, and a stable “memory” phase, reinforced by a broader regulation of H3K36me2 levels.

Because TGFβ is such a powerful EMT agent, we wondered whether Nsd2-deficient cells would fail to respond to this stimulus. Surprisingly, Nsd2-deficient cells treated with TGFβ exhibited a typical EMT response, characterized by the repression of epithelial genes and induction of mesenchymal genes (Fig. 5A). Strikingly, EMT was still accompanied by an increase in H3K36me2 levels despite total absence of Nsd2 (Fig. 5B). This result suggests that in the presence of a strong EMT stimulus, Nsd2-deficient cells compensate by finding other ways to upregulate the mark, indicating that H3K36me2 accumulation itself is essential for mesenchymal identity.

In the absence of Nsd2, other enzymes can still dimethylate H3K36. In fact, upon TGFβ stimulation, Nsd2 paralogs were upregulated to a greater extent in the sgNsd2 clones compared with the WT (Fig. 5C), potentially explaining how Nsd2-deficient cells treated with TGFβ can still generate H3K36me2 and undergo EMT. To circumvent such compensatory mechanisms and more directly interrogate the biological role of this histone mark, we employed a mutant form of H3.3 with a Lys-to-Met substitution at position 36 (K36M), which reduces the overall abundance of H3K36me2 by dominantly inhibiting NSD2 and other H3K36 methyltransferases (Fig. 5D; refs. 36, 47, 48). In contrast to prior observations (35, 36, 47), we did not observe consistent decreases in H3K36me1 or H3K36me3.

Expression of H3.3K36M in cell lines from multiple tumor types, including PDAC (Panc02 and PANC1), lung carcinoma (A549), and prostate carcinoma (DU145), was sufficient to induce an epithelial phenotype (Fig. 5E; Supplementary Fig. S11). However, in contrast to Nsd2-deficient cells, K36M-expressing cells failed to accumulate H3K36me2 following TGFβ treatment (Fig. 5F) and were refractory to its EMT-inducing effects (Fig. 5G and H). These results indicate that increased H3K36me2 is essential for stabilizing the mesenchymal state across various cancer types, and that failure to accumulate H3K36me2 impedes a cell's ability to undergo EMT.

To better understand the epigenomic mechanisms by which H3K36me2 influences epithelial plasticity, we performed a series of chromatin immunoprecipitation sequencing (ChIP-seq) experiments. Although loss of Nsd2 led to global reductions in H3K36me2, some genomic regions maintained this mark (Supplementary Fig. S12A, boxed regions). We therefore performed separate analyses of genomic loci that exhibit dynamic changes in H3K36me2 (“H3K36me2-regulated”) versus loci that retain the mark (“H3K36me2-stable”) in Nsd2-deficient cells. H3K36me2-regulated domains were associated with loss of the activating marks H3K27ac and H3K4me1 and an increase in the repressive mark H3K27me3 (Fig. 6A; Supplementary Fig. S12B, top). In contrast, these changes were not observed in H3K36me2-stable regions (Fig. 6B; Supplementary Fig. S12B, bottom). This suggests that dynamically regulated H3K36me2 domains encompass a subset of H3K27ac-containing regulatory elements.

H3K27ac is a feature of both active promoters and enhancers. Because dynamic changes in H3K36me2 levels are most apparent in intergenic regions (refs. 22, 36, 47, 49; Supplementary Fig. S12C), we hypothesized that H3K36me2 modulates epithelial plasticity by regulating enhancer activity. To test this, we examined the spatial distribution of H3K36me2, analyzing separately genomic regions near transcriptional start sites (putative promoters) versus more distal elements (putative enhancers). Consistent with our hypothesis, almost all H3K36me2-regulated domains were in distal elements, whereas a more substantial fraction of H3K36me2-stable domains was found in promoter regions (Supplementary Fig. S12D). Similarly, H3K27ac peaks lost with Nsd2 depletion were more likely to be found in distal elements than those that were gained or unchanged (Supplementary Fig. S12E). In addition, we compared H3K36me2 and H3K27ac coverage in distal elements and promoters in the presence or absence of Nsd2. Promoter-localized H3K36me2 domains and H3K27ac were maintained with loss of Nsd2 (Supplementary Fig. 12F, right). In contrast, most H3K36me2 domains in distal elements disappeared following Nsd2 loss, a shift that was associated with large decreases in H3K27ac (Supplementary Fig. S12F, left). These results suggest that H3K36me2-regulated domains are enriched for distal regulatory elements (i.e., enhancers), whereas H3K36me2-stable domains are enriched for promoters.

Finally, we examined the effects of H3K36me2 on gene expression by analyzing the transcriptional consequences of Nsd2 loss in H3K36me2-regulated and H3K36me2-stable regions. Genes associated with H3K36me2-regulated domains exhibited reduced expression following Nsd2 loss compared with genes associated with H3K36me2-stable domains (Fig. 6C; Supplementary Fig. S13A). Remarkably, although H3K36me2-regulated domains were widely dispersed across the genome (>10,000 peaks), we found them to be associated with relatively few genes (n = 189), including many prometastatic and EMT master regulatory factors (Fig. 6C; Supplementary Fig. S13B and S13C). In contrast, the expression of genes associated with H3K36me2-stable domains was mostly unchanged or upregulated with Nsd2 loss and included markers of epithelial cell identity or differentiation (Fig. 6C; Supplementary Fig. S13B). Interestingly, H3K36me2-regulated and H3K36me2-stable genes are associated with distinct enhancer and promoter activities. Both enhancer (Fig. 6D, left) and promoter (Supplementary Fig. S13D, top) activities of H3K36me2-regulated genes are directly affected by H3K36me2 levels. And, as expected, the enhancer activities of H3K36me2-stable genes remained largely unchanged (Fig. 6D, right). However, although H3K36me2 levels remained unchanged with Nsd2 loss at the promoters of H3K36me2-stable genes, including Ecad, H3K27ac coverage still increased at these promoters (Supplementary Fig. S13D, bottom), suggesting that transcriptional activation at these sites is independent of H3K36me2 levels. These results indicate that H3K36me2 influences plasticity by altering enhancer activity, as well as corresponding promoter activity, to regulate the transcription of key regulatory factors (i.e., EMT-TFs).

Although previous studies have found that EMT-TFs like SNAI1, SNAI2, and ZEB1 physically and functionally interact with epigenetic modifiers to drive changes in gene expression (10–16), our analyses suggest that there is a level of epigenetic regulation upstream of EMT-TF activity. To test this, we overexpressed Snai1 or Zeb1 in cells lacking Nsd2 (Fig. 6E). We found that either EMT-TF was sufficient to induce an EMT transcriptional program by downregulating epithelial markers and upregulating mesenchymal markers (Fig. 6F). In agreement with our earlier finding that Zeb1 is the major driver of EMT in our PDAC cell lines, Snai1 was less effective than Zeb1 at rescuing the EMT phenotype. Interestingly, and in contrast to Nsd2 (Supplementary Fig. S5), neither Zeb1 nor Snai1 appeared important for the expression of other EMT-TFs, suggesting that the level of H3K36 dimethylation, rather than transcription factor cross-talk, is responsible for broad EMT-TF upregulation during EMT. Finally, neither Zeb1 nor Snai1 overexpression induced the expression of Nsd1, Nsd2, or Nsd3 (Fig. 6F) or resulted in increased H3K36me2 levels in Nsd2-deficient cells (Fig. 6G). Collectively, these results indicate that these master regulatory factors act downstream of H3K36me2 activity to drive the EMT transcriptional program.

Epithelial plasticity is associated with the activity of multiple EMT-TFs, whose expression levels and functions vary across different developmental, homeostatic, and malignant settings. This context-specific variation has complicated efforts to establish unifying principles in their mode of action and downstream consequences. Here, we sought to understand how epigenetic regulators, and their molecular targets, influence epithelial differentiation state. Until recently, the biological functions of specific epigenetic modifications have been inferred from gain- and loss-of-function studies involving their putative modifying enzymes. However, many of these marks, particularly those found on histones, are regulated by multiple, redundant enzymatic writers and erasers. Furthermore, these enzymes have nonhistone substrates and/or act as scaffolds to recruit other factors (50). Consequently, it has been challenging to directly assess the biological impact of these modifications (as opposed to the enzymes catalyzing their regulation). Recently, histone lysine-to-methionine (K-to-M) point mutations have emerged as a powerful tool to study a subset of histone modifications (36, 47, 51, 52). Accordingly, we employed the K36M histone mutant to directly interrogate the function of H3K36me2 in cellular plasticity and its impact on other epigenomic features. Our study shows that gain of the histone mark H3K36me2 in carcinoma cells is a conserved epigenome-wide determinant of mesenchymal identity, thereby demonstrating that epigenomic reprogramming is not merely a consequence of EMT but instead is actively involved in enforcing and/or maintaining cellular states in cancer.

Our finding that the H3K36me2 mark is critical for reinforcing the mesenchymal state is consistent with work in human chondroblastoma, where mutant H3K36M acts as an oncogenic driver by dysregulating many of the pathways involved in normal mesenchymal differentiation (36, 47). Similarly, loss of NSD2 is associated with Wolf–Hirschhorn syndrome, a condition characterized by a collection of midline phenotypes, craniofacial abnormalities, cardiac anomalies, and neurodevelopmental delays, indicative of defects in neural crest (53, 54). EMT is essential for neural crest cells to leave the neural tube during development, and recent studies suggest that NSD2 plays a direct role in neural crest cell migration (55, 56). Therefore, our finding that NSD2 promotes a mesenchymal migratory phenotype in cancer may reflect an oncogenic adaptation of normal developmental processes involving this enzyme and its histone target, implicating the regulation of H3K36me2 in both physiologic and disease-related contexts.

We found that the enzymes involved in global writing and erasing of H3K36me2 play dramatic and distinct roles at different stages of metastatic progression. Specifically, Nsd2 loss inhibits EMT, thereby reducing invasion and bloodstream entry, whereas Kdm2a loss inhibits MET, thereby reducing metastatic colonization. Neither of these perturbations influences cell proliferation or tumor growth, suggesting that flexibility in cellular state, rather than primary tumor size, is a major determinant of metastatic progression. These results join a growing body of evidence implicating epigenetic deregulation in metastasis (57–59), where the phenotypic plasticity afforded by epigenetic reprogramming may allow cells to overcome the diverse selective pressures encountered at distinct stages of metastatic spread.

An inverse relationship between H3K36me2 and repressive polycomb activity has been proposed as a mechanism by which H3K36me2 regulates gene expression (30, 36, 47). Although we also observed an inverse relationship between H3K36me2 and H3K27me3 (the target of polycomb), our study found that H3K27ac peaks residing within broad intergenic H3K36me2 domains are lost when H3K36me2 levels decrease, indicating that H3K36me2 also mediates its effects by modulating enhancer activity. We were surprised to find that despite the global distribution of H3K36me2, relatively few genes (<200) were transcriptionally affected by dynamic changes in this histone mark. Remarkably, this list included virtually all key EMT-TFs, as well as several prometastatic factors, including members of the Zeb, Snai, Prrx, Sox, and Smad gene families. Our finding that NSD2 regulates epithelial plasticity by altering enhancer activity expands on previous work implicating NSD2 in the regulation of EMT and EMT-TF expression in prostate cancer (23–25). In contrast, genes associated with the epithelial state, including various cell adhesion molecules (and Ecad itself), appeared to be regulated at their promoters independent of dynamic changes in H3K36me2. These results support a hierarchical model of epithelial plasticity involving two tiers of control: a first tier in which transcriptional regulation of epithelial and mesenchymal programs is achieved by the acute activity of EMT-TFs, and a higher-order tier dependent on H3K36me2-mediated reprogramming of enhancers associated with these EMT-TFs, thereby achieving a more stable epigenetic memory of epithelial-to-mesenchymal state.

Finally, our results have translational implications for cancer therapy. In addition to its role in cancer metastasis, epithelial plasticity is a major mechanism of therapeutic resistance to chemotherapy and targeted agents (40). However, the signaling pathways that drive EMT and MET are highly heterogeneous and converge on diverse transcriptional regulators, making it challenging to develop strategies to inhibit these processes. Our finding that H3K36me2 regulation is a conserved and required feature of carcinoma-associated epithelial plasticity may represent a way around this barrier, as the enzymes mediating this histone PTM are all potentially targetable with small molecules. There is growing evidence that inhibitors of epigenetic modifiers can be both highly specific and potent (60), paving the way for strategies to target epithelial plasticity and other transcriptional processes. For example, other potential targets from our screen, such as KDM1A, which has been shown to interact with SNAI1 to drive EMT (13), and DOT1L, which was recently demonstrated to regulate EMT in breast cancer (61), already have inhibitors that are currently being tested in clinical trials for various cancers. Nevertheless, in terms of targeting H3K36me2, the existence of at least five histone methyltransferases capable of catalyzing its formation, NSD1, NSD2, NSD3, ASH1L, and SETMAR, suggests that redundant and/or context-dependent mechanisms may yet complicate such efforts, and further work will be needed to clarify their roles.

Kras^LSL-G12D; Trp53^L/+; Pdx1-cre; Rosa26^YFP/YFP (KPCY) mice have been described previously (19). Female 6- to 8-week-old C57BL/6J or NOD/SCID mice were purchased from The Jackson Laboratory and Charles River Laboratory, respectively, for tumor cell injection experiments. All vertebrate animals were maintained and experiments were conducted in compliance with the NIH guidelines for animal research and approved by the University of Pennsylvania Institutional Animal Care and Use Committee.

Murine PDAC cell lines 3077c4, 483c6, and c8 were derived from primary KPCY tumors of mixed genetic background (20). The 6419 c5 and 6694 c2 cells were derived from primary tumors from KPCY mice that have been backcrossed onto the C57BL/6J mouse strain (62). The murine Panc02 cell line was provided by Dr. Robert Vonderheide, University of Pennsylvania, Philadelphia, PA. The human PDAC PANC1 cell line was provided by Dr. Anil Rustgi, Columbia University, New York, NY. The human prostate cancer cell line DU145 was provided by Dr. Irfan Asangani, and the human lung adenocarcinoma cell line A549 was provided by Dr. Kathryn Wellen, University of Pennsylvania. All human lines were originally obtained from ATCC. All murine lines and the PANC1 cell line were single-cell FACS to establish clonal cell lines.

The 3077 and 483 cell lines were cultured in PDEC media, as described previously (19). All other murine and human cell lines were cultured in DMEM supplemented with 10% heat-inactivated FBS, 10% glutamine, and 150 μg/mL gentamicin at 37°C, 5% CO₂, 21% O₂, and 100% humidity. Cell lines were maintained and passaged according to ATCC recommended procedures. Cells lines were regularly tested for Mycoplasma using MycoAlert Mycoplasma Detection Kit (Lonza).

For studies interrogating individual genes, sgRNA-encoding oligonucleotides (Supplementary Table S4) were inserted into the pSpCas9(BB)-2A-Puro (PX459) V2.0 vector, a gift from Feng Zhang, MIT and the Broad Institute, Boston, MA (Addgene plasmid #62988), using a BsmBI restriction site. shRNA oligonucleotides (Supplementary Table S5) were designed using the splashRNA algorithm (63), and cloned into the LT3REPIR vector (a gift from Scott Lowe, Memorial Sloan Kettering Cancer Center, New York, NY), and modified from the LT3GEPIR vector, Addgene, plasmid #111177, such that GFP is replaced by dsRed], as described previously (64). Briefly, XhoI and EcoRI restriction sites were PCR cloned onto the shRNA oligonucleotides using the 5′ miR-XhoI primer TACAATACTCGAGAAGGTATATTGCTGTTGACAGTGAGCG, the 3′ miR-EcoRI primer TTAGATGAATTCTAGCCCCTTGAAGTCCGAGGCAGTAGGCA, and the NEBNext High-Fidelity 2x PCR Master Mix (New England Biolabs), according to manufacturer's protocol. Amplified shRNA oligos and LT3REPIR vector were digested with XhoI and EcoRI and ligated together.

H3.3 K36M mutation was introduced to the pcDNA4/TO-Flag-H3.3 plasmid (a gift from Bing Zhu, National Institute of Biological Sciences, Beijing, China; Addgene plasmid #47980) by site-directed mutagenesis. The Q5 site-directed mutagenesis kit (New England Biolabs) was used according to the manufacturer's protocol with the following mutagenic primers: forward, GGAGGGGTGATGAAACCTCATC; reverse, AGTAGAGGGCGCACTCTT. Following confirmation of the point mutation by Sanger sequencing, the epitope-tagged H3.3 K36M was cloned into the pCDH-EF1-FHC lentiviral vector (a gift from Richard Wood, The University of Texas MD Anderson Cancer Center; Addgene plasmid #64874) for constitutive expression using NotI and BamHI restriction sites.

Cultured tumor cells were dissociated into single cells with Hank's based enzyme-free cell dissociation solution (EMD Millipore) and washed in Hank's Balanced Salt Solution with 5% FBS and DNase I (Sigma). Cells were then stained with rat anti–E-cadherin (1:250; akara, M108), followed by APC donkey anti-rat (1:100; Jackson Immunoresearch). For FACS, samples were filtered through a 70-μm strainer to form single-cell suspensions. Flow cytometric analysis was performed on a LSR II flow cytometer (BD Biosciences) and analyzed using FlowJo software (Treestar). FACS was performed on a FACS Aria II sorter (BD Biosciences).

A list of sgRNA sequences directed against a comprehensive murine epigenetic modifier gene list was a gift from Scott Lowe, and a sgRNA list targeting the functional and catalytic domains of murine epigenetic modifiers was taken from a previous study (65). In addition, 108 sgRNAstargeting EMT regulatory factors, 500 nontargeting sgRNAs, and 50 sgRNAs targeting essential genes were added. In total, the library contains 4,803 sgRNAs. All sgRNAs were designed using design principles previously reported and were filtered to minimize predicted off-target effects (66). sgRNAs were synthesized in a pooled format on an array platform (CustomArray, Inc.) and then PCR cloned into the LRG2.1 vector, a gift from Christopher Vakoc, Cold Spring Harbor Laboratory, Cold Spring Harbor, NY (Addgene plasmid #108098), using previously described methods (67). To ensure proper representation of sgRNAs in the pooled lentiviral plasmids, the library was analyzed by deep sequencing on a NextSeq500/550 (Illumina, 75 cycles High Output Kit v2.0) and MAGeCK (68), which confirmed that 100% of the designed sgRNAs were cloned in the LRG2.1 vector with relatively even read-count distribution.

Transfection of 293T cells for lentivirus production was performed using Opti-MEM I (Gibco), DNA (expression plasmid, psPAX2, and pVSVg mixed in a 4:2:1 ratio), and polyethylenimine (PEI; Polysciences) at a 3:1 ratio with total DNA. Recipient cells were transduced with filtered, unconcentrated viral supernatant in the presence of 8 μg/mL polybrene. Appropriate antibiotic selection (2–8 μg/mL puromycin, 10 μg/mL blasticidin) was applied for at least 72 hours. Viral titering for sgRNA library was done by transduction at limiting dilution and determining the GFP⁺ percentage 72 hours later using the GUAVA easyCyte HT BG flow cytometer (Millipore).

Full-length human NSD2 and E1099K NSD2 mutant pMSCV-puro expression plasmids were kindly provided by I. Asangani. Full-length KDM2A (GFP-FBXL11) was a gift from Michele Pagano, New York University, New York, NY (Addgene plasmid #126543). Full-length murine Snai1 and Zeb1 were gifts from Thomas Brabletz, University Erlangen, Germany, and cloned into the pCW57.1 vector using AvrII/BamHI and MluI/BamHI restriction sites, respectively. 293T cells were transfected using Opti-MEM I (Gibco), DNA (expression plasmid, HIV gag-pol, and pVSVg mixed in a 4:2:2 ratio), and PEI (Polysciences) at a 3:1 ratio with total DNA. Virus collection and transduction were performed as previously detailed. Puromycin selection (5 μg/mL) was begun 48 hours post-transduction and continued for 1 week.

Cas9 was introduced to the 3077 c4 cell line by transduction of the lentiCas9-Blast construct (a gift from Feng Zhang; Addgene plasmid #52962), and selected with 10 μg/mL blasticidin. A clonal 3077 c4-Cas9 line was generated by single-cell FACS and was later transduced with the epigenetic modifier library at a multiplicity of infection of 0.2 such that approximately 20% of cells were GFP+. 210 million cells were transduced to yield 1,000× coverage of the sgRNA library in both the Ecad⁺ and Ecad⁻ subpopulations. T0 GFP⁺ Ecad⁺ and GFP⁺ Ecad⁻ samples were collected by FACS 3 days post-transduction. GFP⁺ Ecad⁺ and GFP⁺ Ecad⁻ samples were collected at 2 weeks (T2weeks) and 4 weeks (T4weeks) post-transduction. Enough cells were collected to maintain at least 500x coverage of the sgRNA library in each sample.

Genomic DNA was harvested using the QIAamp DNA Mini Kit (Qiagen), and libraries were constructed according to previously described protocols (69). Briefly, sgRNAs were amplified over 18 cycles with Herculase II fusion DNA polymerase (Agilent) per manufacturer specifications with PCR#1 forward and reverse primers (Supplementary Table S6). Each reaction included 3 μg of genomic DNA, and multiple PCR reactions were run in parallel such that all extracted genomic DNA was used to maintain library coverage. PCR reactions were then pooled for each sample, and 5 μL of each pooled PCR#1 sample was used as a template for PCR#2, which added Illumina P5/P7 adapters, barcodes, and staggers for nucleotide complexity. For PCR#2, the template was amplified over 25 cycles with PCR#2 forward and reverse primers (Supplementary Table S6), and the resulting reactants were column purified with the QIAquick PCR Purification Kit (Qiagen) and gel extracted with the QIAquick Gel Extraction Kit (Qiagen). The barcoded libraries were pooled at an equal molar ratio and sequenced on a NextSeq500/550 (Illumina, 75 cycles High Output Kit v2.0) to generate 75 bp single-end reads.

MAGeCK was used for subsequent screen analysis (68). Briefly, the sequencing data were de-barcoded and the 20 bp sgRNA sequence was mapped to the reference sgRNA library without allowing for any mismatches. The read counts were calculated for each individual sgRNA and normalized with the nontargeting sgRNAs. Normalized read counts of sgRNAs in Ecad⁺ and Ecad⁻ subpopulations were log₂ transformed, and graphical representation was done using the R-package ggplot2. To exclude false positives due to off-target effects or inconsistencies in sample preparation, we focused only on hits for which >80% sgRNAs targeting a gene were enriched.

Pancreatosphere assays were performed as described by Rovira and colleagues (70). Briefly, cell lines were trypsinized and dissociated into single cells. Cells (2,500) were suspended in tumorsphere media: DMEM/F-12 with 1x B-27 supplement (Gibco), 3% FBS, 100 mmol/L B-mercaptoethanol (Gibco), 1x nonessential amino acids (Sigma), 1x N2 supplement (Gibco), 20 ng/mL mouse EGF (Gibco), 20 ng/mL mouse FGF-basic (Gibco), and 10 ng/mL mouse LIF (Gibco). Cells were plated in 6-well ultra-low attachment culture plates (Costar) and spheres were counted after 5 days. Five independent wells were used per cell line.

Mice were anaesthetized using continuous-flow isoflurane, and their abdomens were sterilized. An incision was then made over the left upper quadrant of the abdomen, and the pancreas was exteriorized onto a sterile field.Tumor cells (1.0 × 10⁴) in 50 μL sterile DMEM were injected into the tail of the pancreas via a 27 5/8″ gauge needle. Successful injection was confirmed by the formation of a liquid bleb at the site of injection with minimal fluid leakage. The pancreas was then gently placed back into the peritoneal cavity, and the peritoneum and overlying skin were closed with 4-0 coated Vicryl violet FS-2 sutures (Ethicon). Tumors, lungs, and livers were harvested, weighed, and measured 4 to 6 weeks following implantation. Two hundred microliters of blood was drawn from tumor-bearing animals via cardiac puncture with a 1 mL insulin syringe coated with 0.5 mol/L EDTA pH 8.0 (Gibco) to prevent coagulation and was immediately placed in a 150 mm gridded plate containing RBC lysis buffer (BD Biosciences). After 10 minutes of lysis at room temperature, PBS was added to the plate and CTCs were directly visualized on a fluorescent microscope and counted.

Tumor cells (1.0 × 10⁵) in 100 μL sterile DMEM were injected into the tail veins of mice via a 27 5/8-gauge needle. After 2 weeks, lungs were harvested and weighed.

Cells were plated to 70% confluency in 6-well plates and grown in serum-free media for 24 to 48 hours. A horizontal scratch was then made in each well with a 200 μL pipette tip. Cells were photographed at identical points at 0, 24, and 48 hours post-scratch. Cell migration was assessed by measuring the area of the scratch at each time point using ImageJ.

Two thousand cells were seeded into each well of a 96-well plate. EdU was spiked into the media at 10 μmol/L, and 2 hours later cells were fixed and permeabilized with the fixation/permeabilization buffer set (eBioscience) for 30 minutes at room temperature. Cells were then washed, and EdU was labeled with Alexa Fluor 647 azide (1:2,000; Thermo Fisher Scientific) by a click chemical reaction (100 mmol/L Tris-HCl pH 8.5, 4 mmol/L CuSO₄, 100 mmol/L ascorbic acid) for 30 minutes at room temperature, washed, and analyzed by flow cytometry.

Cells (1 × 10⁵) were plated into 6-well plates and 24 hours later were treated with either citric acid (control) or 10 ng/mL TGFβ (Cell Signaling Technology) over a given time course. Cells were continuously dosed every 3 days and were passaged in the given agents over longer time points.

RNA was prepared from cultured tumor cells or sorted cells using RNeasy Mini Kit or RNeasy Micro Kit (Qiagen). For qRT-PCR, cDNA was generated using the High-Capacity cDNA Reverse Transcription Kit (Life Technologies). qRT-PCR analysis was performed with SsoAdvanced SYBR (Bio-Rad) using a CFX384 Real-Time System (Bio-Rad). Transcript quantities were determined using the difference of C_t method, and values were normalized to the expression of Tbp. Primer sequences are listed in Supplementary Table S7.

Sequencing libraries were generated using the NEBNext Ultra RNA Library Prep Kit (New England Biolabs) and sequenced in technical duplicates on an Illumina HiSeq platform to generate 150 bp paired-end reads (GSE137382) by Novogene. Reads were aligned to the mm10 reference genome using STAR (v2.5.1; ref. 71), and raw counts of gene transcripts were obtained using featureCounts (72), both with default settings. The raw count matrix was subsequently imported into R-studio and used as input for DESeq2 (73) for normalization and differential gene expression analysis. Salmon was used in parallel to normalize and quantify gene expression in transcripts per million through quasi-alignment. To generate ranked gene lists for preranked GSEA (http://www.broad.mit.edu/gsea/), genes were ranked by their DESeq2-derived Wald statistic. The gene signatures for Squamous and Progenitor subtypes were derived from (ref. 34). Briefly, the top 500 and bottom 500 protein-coding genes of each subtype were identified from RNA-seq data derived from patient samples (31), and used to define each subtype's gene signature.

ATAC-seq samples were prepared by sorting 50,000 Ecad⁺ and 50,000 Ecad⁻ subpopulations from three clonal cell lines. The libraries were prepared as previously described with minor modifications (74). Briefly, nuclei were isolated from sorted cells using a solution of 10 mmol/L Tris-HC pH 7.4, 10 mmol/L NaCl, 3 mmol/L MgCl2, and 0.1% IGEPAL CA-630. Immediately following nuclei isolation, the transposition reaction was conducted using Tn5 transposase and TD buffer (Illumina) for 45 minutes at 37°C. Transposed DNA fragments were purified using a Qiagen MinElute Kit, barcoded with primers based on Illumina TruSeq indices (Supplementary Table S8), and PCR amplified for 11 cycles using NEBNext High Fidelity 2x PCR master mix (New England Biolabs). Libraries were purified by extraction from a 6% TBE gel, followed by column purification with the Qiagen PCR Cleanup Kit. Sequencing was performed on biological triplicates using a HiSeq2500 (Illumina) on rapid run mode to generate 50 bp paired end reads (GSE137520).

After adapter trimming with cutadapt, reads were aligned to the mouse reference genome (GRCm38; May 23, 2014), which was downloaded from the UCSC repository, by Bowtie2 (75). Picard (Broad Institute) was used to identify and discard reads that aligned to the mitochondrial genome as well as reads mapping to multiple genomic loci. MACS2 with ‘-p 1e-7–nolambda–nomodel' was applied on each ATAC-seq replicate separately to identify accessible chromatin. Peaks were subsequently merged using BEDTools (76), and ATAC-seq read counts were calculated in the merged peaks for every replicate. The resulting count table was used to identify differentially accessible loci with DESeq2 in R-studio, which were then associated with their putative target genes by GREAT (77). The GREAT algorithm associates genomic regions with genes by defining a ‘regulatory domain’ for each gene in the genome, then computes ontology term enrichments using a binomial test. Visualization of these loci was achieved by uploading the reads-per-million-aligned-read-normalized bigwig files to Integrative Genomics Viewer (78).

For ChIP, approximately 2 × 10⁷ cells were cross-linked with 1% formaldehyde for 10 minutes at room temperature, followed by quenching with 125 mmol/L glycine for 5 minutes at room temperature. They were then washed with PBS and lysed with nuclei isolation buffer. Samples were sonicated with a Diagenode Bioruptor at high intensity setting for 15 five-minute cycles of 30 seconds on/off, to yield approximately 300-bp fragments. After high-speed centrifugation, soluble chromatin from Drosophila S2 cells was spiked in at 2.5% of the mouse chromatin. The mixed soluble chromatin was diluted 10x and precleared with magnetic protein A Dynabeads (Invitrogen). Immunoprecipitation with anti-H3K36me2 (Cell Signaling Technology, 2901S), anti-H3K27me3 (Cell Signaling Technology, 9733S), anti-H3K27ac (Abcam, ab4729), anti-H3K4me1 (Abcam, ab8895), or rabbit IgG Isotype (Abcam, ab171870) was performed overnight at 4°C with rotation, with 10% kept as input DNA. Immunocomplexes were recovered with magnetic Protein A dynabeads and washed sequentially using a low-salt buffer, a high-salt buffer, lithium chloride buffer, and TE buffer. Samples were eluted twice with sodium bicarbonate. Eluted ChIP DNA was treated with RNase A (0.2 mg/mL) overnight at 65°C, then incubated with Proteinase K (0.2 mg/mL) for 2 hours at 45°C. DNA samples were purified using QIAquick PCR Purification Kit (Qiagen). For ChIP sequencing, libraries were prepared according to the NEBNext Ultra II DNA Library Prep Kit (New England Biolabs) protocol, using barcoded primers from NEBNext Multiplex Oligos Index Primers Sets 1 and 2. Libraries were pooled and sequenced in technical duplicates on the NextSeq500/550 (Illumina, 75 cycles High Output Kit v2.0) to generate 75 bp single-end reads (GSE137521).

Sequencing reads were aligned to the mouse (GRCm38; May 23, 2014) and the Drosophila melanogaster reference genomes (dm6; August 2014) using STAR (v2.5.1) with default settings. Duplicate and low-quality (PHRED score <30) reads were removed using SAMtools (79). The normalization factor was determined from the Drosophila reads as described previously (80).

The bam files were converted to ChIPRX-normalized bigwig tracks using deepTools (81). Single-end reads were extended up to the fragment length (200 bp) along the read direction. Enrichment peaks were called with MACS2 with default settings (H3K27ac and H3K4me1) and broad settings (H3K36me2 and H3K27me3). Blacklisted peaks (https://sites.google.com/site/anshulkundaje/projects/blacklists) were excluded in all subsequent analysis. To identify peaks differentially bound between WT and sgNsd2 samples, we first identified the consensus peak set that includes the unique set of peaks enriched across both samples. Differentially bound peaks were then identified using DESeq2 with the ChIPRX normalized tag counts at the consensus peak regions as input, and only peaks with P_adj <0.001 were used for further analysis. Overlap analysis of peaks was performed using an in-house python script, and peaks were annotated using HOMER (http://homer.ucsd.edu/homer/ngs/annotation.html).

Tissues were fixed in zinc-formalin and embedded in paraffin for histologic analysis and immunofluorescence staining. Sections were deparaffinized, rehydrated, and prepared by antigen retrieval. For cell lines, cells were seeded into 8-well Nunc Lab-Tek II chamber slides (Thermo Fisher Scientific) and fixed in 4% paraformaldehyde for 15 minutes. For immunofluorescence staining, sections or fixed cells were blocked in PBS with 0.3% Triton-X and 5% donkey serum for 1 hour, stained with primary and secondary antibodies, and mounted with Aqua Polymount (Polysciences, Inc). Primary antibodies used include goat anti-GFP (Abcam, ab6673), rat anti–E-cadherin (Takara Bio, M108), and rabbit anti-H3K36me2 (Abcam, ab9049). Slides were visualized using an Olympus IX71 inverted multicolor fluorescent microscope equipped with a DP71 camera. Images were quantified using ImageJ software.

For hematoxylin and eosin staining, sections were deparaffinized, rehydrated, stained with hematoxylin, differentiated with acidic ethanol, stained for eosin, dehydrated, and mounted with Permount. Slides were visualized using the Keyence BZ-X710 all-in-one fluorescence microscope.

Histones were acid-extracted according to standard protocols. Briefly, nuclei were isolated from cells lysed with a hypotonic lysis buffer. Histones were extracted from the nuclei with 0.2 N H₂SO₄, precipitated in 33% trichloroacetic acid, washed with acetone, and resuspended in diH2O.

For whole-cell lysates, cells were washed with cold PBS and lysed in RIPA lysis buffer. Extracted histones or whole-cell lysates were separated by SDS-PAGE, transferred to polyvinylidene difluoride membrane, blocked in 5% nonfat milk in PBS plus 0.1% Tween-20, probed with primary antibodies, and detected with horseradish peroxidase–conjugated secondary antibodies (Jackson Immunoresearch). Primary antibodies used include: anti-H3 (Cell Signaling Technology, 9715S); anti-H3K36me2 (Cell Signaling Technology, 2901S); anti-H3K36me3 (Abcam, ab9050); anti-H3K36M (Millipore, ABE1447); anti-NSD2 (Millipore, MABE191); anti-KDM2A (Abcam, ab191387); anti-alpha-tubulin (Cell Signaling Technology, 3873S); anti–E-cadherin (Takara, M108).

Histones were derivatized according to standard protocol (82). Briefly, 5 to 10 ug of acid-extracted histones were resuspended in 100 mmol/L ammonium bicarbonate (pH 8.0) and mixed with freshly prepared propionic anhydride with acetonitrile for 15 minutes at 37°C. Histones were then digested with trypsin (enzyme: sample ratio 1:20) overnight at 37°C. The peptides were then desalted and stored dried. They were resuspended in 0.1% formic acid just before mass spectrometry.

Samples were placed in a TriVersa NanoMate (Advion) and acquired either manually or by using a sequence coordinated with MS acquisition by a contact closure. The NanoMate was set up with a spray voltage of 1.7 kV and a gas pressure of 0.5 psi. Samples were acquired in the Orbitrap Fusion Tribrid (Thermo Fisher Scientific). All scans were acquired in the orbitrap, at 240,000 resolution for the full MS and at 120,000 resolution for MS-MS. The AGC target for the tSIM-MSX scans was set to 10E6. The full description of the DI-MS acquisition method has been previously described (83).

Raw files were searched with a modified version of the software EpiProfile 2.0 (84). Histone peptides were collected in MS scans, and isobaric peptides were collected in targeted preset MS-MS scans. The software reads the intensities from MS scans to calculate the percentage of all peptides with the same amino acid sequence. The unique fragment ions in the MS-MS scans were extracted to discriminate isobaric peptide intensities from the MS scans. The software EpiProfileLite is available on GitHub at https://github.com/zfyuan/EpiProfileLite, including the user guide.

Comparisons between two groups were performed using Student unpaired t test. The Holm–Sidak correction for multiple testing was applied where indicated. All statistical analyses were performed with Graphpad Prism 7 and 8 (GraphPad). Error bars show SEM, as indicated in the legends, and P < 0.05 was considered statistically significant. (* indicates P < 0.05; **, P < 0.01; ***, P < 0.001; ****, P < 0.0001 unless otherwise indicated; ns, not significant).

For sequencing experiments, DESeq2 (73) used to generate P_adj (FDR). For The Cancer Genome Atlas survival analysis, KM plotter (40) was used to generate log-rank P values.

All sequencing data has been deposited in the Gene Expression Omnibus under the series GSE137523, which includes the following subseries: GSE137382 (RNA-seq), GSE137520 (ATAC-seq), and GSE137521 (ChIP-seq).

PRISM software was used for the statistical analysis and data visualization (http://www.graphpad.com). The R language and environment for statistical computing and graphics (https://www.r-project.org) was utilized in this study for the statistical and bioinformatics analysis of RNA-seq, ATAC-seq, and genome-seq data. The R packages used for the analysis described in the method section were obtained from the Bioconductor (https://www.bioconductor.org) and CRAN (https://cran.r-project.org/web/packages/).

No potential conflicts of interest were disclosed.

Conception and design: S. Yuan, B.Z. Stanger

Development of methodology: S. Yuan, F.J. Sanchez-Rivera, J. Li, N.V. Bhanu, J. Shi, S.W. Lowe, I.A. Asangani, B.Z. Stanger

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): S. Yuan, J. Li, N.V. Bhanu, T. Yamazoe, A.J. Merrell, Y. Sela, S.K. Thomas, Y. Jiang, J.B. Plesset, E.M. Miller, B.A. Garcia, I.A. Asangani

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): S. Yuan, R. Natesan, J. Li, N.V. Bhanu, J.H. Lin, I.A. Asangani, B.Z. Stanger

Writing, review, and/or revision of the manuscript: S. Yuan, R. Natesan, S.W. Lowe, I.A. Asangani, B.Z. Stanger

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): S. Yuan, B.Z. Stanger

Study supervision: B.Z. Stanger

Other (student research assistant): E.M. Miller

Other (input into the functional genomics methodology and reviewing the manuscript): S.W. Lowe

This work was supported by NIH grants CA229803 and DK083355 (to B.Z. Stanger), CA224970 (to S. Yuan), CA196539 and AI118891 (to B.A. Garcia), CA249210 (to I.A. Asangani), a V Foundation Translational award (to B.Z. Stanger), a Leukemia and Lymphoma Robert Arceci Award (to B.A. Garcia), and DOD W81XWH-17-1-0404 (to I.A. Asangani), the Abramson Family Cancer Research Institute, the Abramson Cancer Center, and the NIH/Penn Center for Molecular Studies in Digestive and Liver Diseases. We thank John Tobias in the Penn Molecular Profiling Facility, Jonathan Schug in the Penn Next Generation Sequencing Core, Kathryn Sun in the Penn Pancreatic Cancer Research Center, Ken Zaret, Andres Blanco, and Liling Wan for helpful discussion of the study and manuscript, Grace Anderson and Donita Brady for technical guidance of the screen, and all members of the Stanger laboratory for helpful discussions.