Generation of EpiDoKOL and validation of the performance of the sgRNA library

To identify epigenetic modifier genes indispensable for cell survival, we generated a customized epigenome-wide domain-targeted pooled sgRNA library. In Epigenetic Domain-specific Knockout Library (EpiDoKOL), sgRNAs were designed against gene regions encoding functional catalytic domains of chromatin modifiers. EpiDoKOL consists of 1750 sgRNAs in total, targeting 250 different genes, and 80 non-targeting controls (Fig. 1A). Molecular cloning of sgRNAs for pooled libraries into two different backbones (pLentiGuide and pLentiCrisprv2) was achieved, leading to 800X coverage and successful lentiviral packaging (Fig. 1B, Supplementary Fig. 1A, B).

Fig. 1

Epigenetic Domain-specific Knock Out Library (EpiDoKOL) screen identified common essentiality genes in T98G and U373 cells. A Composition of EpiDoKOL including target molecules, their functions and related epigenetic modifications. B Library generation procedure. Figure created with BioRender.com. C Cas9 activity assay of U373 and T98G cells. Microscopic images of T98G and U373 cells transduced with indicated sgRNAs 18 days post-transduction. D Flow cytometric analysis of U373 and T98G cells transduced with indicated sgRNAs (g-NT1, g-NT2: non-targeting sgRNAs; g-T1 and g-T2: sgRNAs targeting GFP), up to 12 days post-transduction. E Schematic of EpiDoKOL screening procedure. Figure created with BioRender.com. F sgRNA density plots from cells transduced with pLentiCRISPRv2 plasmid containing EpiDoKOL. Cell pellets collected before puromycin selection. G Log2 counts of sgRNAs at initial and final time points in T98G and U373 cells. H Waterfall plots for Log2fold changes of genes after screening T98G cells with EpiDoKOL for a month. Domain and gene-based analysis revealed overlapping essentiality hits. I Waterfall plots for Log2fold changes of genes after screening U373 cells with EpiDoKOL for a month. Domain and gene-based analysis revealed overlapping essentiality hits. Mageck p values were calculated with negative binomial model fitting. Other p values were determined by two-tailed Student’s t-test *P < 0.05, **P < 0.01, ***P < 0.001

As Cas9-based genetic ablation may be suboptimal in certain tumor cell lines due to chromosomal aberrations and mutational burden, we first examined the efficiency of Cas9-mediated cleavage efficiency in glioblastoma cell lines. We first transduced U373 and T98G with GFP-encoding lentiviruses and then introduced g-NT1, g-NT2 (non-targeting sgRNAs), g-T1 and g-T2 (sgRNAs targeting GFP) along with Cas9 in a single lentiviral vector. Flow cytometry indicated that GFP+ cells were reduced to 20% of the total population at 12 days post-transduction in g-T1 or g-T2 transduced cells. No change in GFP signal was observed with g-NT1 and g-NT2 vectors, indicating that Cas9-activity is optimal within 2 weeks (Fig. 1C, D). We then continued with EpiDoKOL screen in both U373 and T98G cells. To identify common regulators of glioblastoma cell survival, cells were transduced with EpiDoKOL at a low multiplicity of infection (MOI = 0.4) with 500 × coverage to ensure single sgRNA intake and proper sgRNA representation (Fig. 1E). Following antibiotic selection, transduced cells were cultured for 30 days by maintaining the sgRNA coverage at each population doubling. Using standard library preparation protocols, deep sequencing and analysis [27], we compared the sgRNA composition in initial and final timepoints. Histogram of median normalized read counts for all sgRNA plotted for U373 and T98G cells confirmed normal distribution of log2-transformed normalized counts, indicating that no bias was introduced during cloning or transduction steps (Fig. 1F). We then compared the depletion scores of gene-targeting versus non-targeting sgRNAs to evaluate the overall screen efficiency. Essential gene targeting sgRNAs were depleted significantly, whereas no change in non-targeting control sgRNA abundance was observed (Fig. 1G). Altogether, these analyses illustrated that EpiDoKOL preserves normal distribution of sgRNAs in glioblastoma cells and reveals potential essential genes for cell survival.

EpiDoKOL screen in T98G and U373 cells identified common regulators of glioblastoma cell survival

As EpiDoKOL was composed of multiple sgRNAs targeting a functional domain of a single gene, it was possible that one gene was targeted by more than 10 sgRNAs in total, depending on the domain composition. Therefore, we undertook two different analysis approaches; at domain-level or gene-level and determined common depleted genes using median normalization scores (Fig. 1H, I). Accordingly, several epigenetic modifiers were found to be commonly depleted in glioblastoma cell lines. Along with epigenetic modifiers that were previously implicated in cancer cell fitness, such as CHD1, CHD4, DNMT3B, ELP3, SUPT16H, SUV39H2; novel hits ASH2L, RBX1 and SSRP1 were discovered as essential genes for both U373 and T98G cells suggesting common regulatory role (Supplementary Fig. 1C).

To validate the function of these novel hits in glioblastoma, we individually introduced selected sgRNAs and assessed their effects with several functional experiments (Fig. 2A). Accordingly, colony forming abilities of U373 and T98G cells with individual knockouts of RBX1, ASH2L, SSRP1 genes markedly decreased compared to control cells that received g-NT (Fig. 2B). To assess the effect of hit genes in a heterogenous population, we performed GFP competition assays, where GFP-negative cells were transduced with sgRNA against a hit gene and mixed in a 1:1 ratio with gNT-transduced GFP-positive cells. Monitoring both the GFP+ and GFP- fraction over 16 days, we observed gradual decreases in the cells transduced with selected sgRNA in U373 and T98G cells (Fig. 2C), and U87MG cells, as a third glioblastoma cell line (Supplementary Fig. 2A). The most significant phenotype was observed for U373 cells and we further delineated the effects of RBX1, ASH2L, SSRP1 ablation on glioblastoma cell fitness using this cell line. Depletion of ASH2L, RBX1 or SSRP1 genes slowed down the growth of cells as gauged by cell viability assays (Fig. 2D); and induced apoptosis as revealed by elevated caspase 3/7 activity (Fig. 2E) and AnnexinV positivity (Fig. 2F). Thus, apoptotic cell death is likely to contribute to the reduced cell numbers observed upon ASH2L, RBX1 and SSRP1 ablation. These proof-of-principle experiments illustrated that EpiDoKOL is a practical functional genomics tool that enables identification of epigenetic modifiers important for cancer cell fitness.

Fig. 2

Effects of candidate genes on glioblastoma cell fitness were validated with functional assays in vitro. A Scheme of validation experiments for novel EpiDoKOL essentiality hits. TD: post-transduction day. B Representative images of long-term clonogenic assay of cells infected with sgRNAs against selected hits and statistical analysis. Quantification of colonies was performed by ImageJ software. C Results of GFP competition flow cytometric assay for selected hits. Day 0 refers to the day of the cell seeding at TD 7. D-F Cell viability (D), Caspase 3/7 activity (E) and AnnexinV analysis (F) conducted on U373 cells upon depletion of ASH2L, RBX1 and SSRP1 genes. G Western blot analysis of ASH2L protein levels upon transduction of U373 cells with g-ASH2L or g-NT. H Western Blot analysis for cleaved Caspase3 and PARP in U373 cells at day 6 and 14 post-transduction with g-NT control or g-ASH2L. I Western Blot analysis for H3K4 mono and trimethylation levels in U373 cells 6 days post-transduction with g-NT control or g-ASH2L. P values determined by two-tailed Student’s t-test *P < 0.05, **P < 0.01, ***P < 0.001

ASH2L is essential for glioblastoma cell survival and regulates histone methylation and transcription

Among the genes validated for their essential function in glioblastoma, ASH2L was the least studied in cancer and its ablation gave a stronger phenotype; therefore we further investigated its role in glioblastoma. We knocked out ASH2L in U373 and U87MG cells and showed successful depletion of ASH2L in protein level (Fig. 2G, Supplementary Fig. 2B). Extensive cleavage of caspase-3 and PARP was observed indicating potent induction of apoptosis upon ASH2L depletion (Fig. 2H). ASH2L functions as a cofactor within MLL family of methyltransferase complexes to trigger histone H3 lysine 4 (H3K4) methylation [16,17,18,19], and its loss was previously associated with reduction of H3K4 trimethylation. Consistently, we observed reduced mono and trimethylation of H3K4 upon ASH2L depletion in U373 and U87MG cells (Fig. 2I, Supplementary Fig. 2C). Since H3K4me3 is a well-known euchromatin mark acting globally, we investigated gene expression changes upon ASH2L knock-out to gain a mechanistic understanding for its essential role. We performed RNA-seq comparing U373 cells transduced with control and ASH2L sgRNAs (Fig. 3A). ASH2L knockout resulted in 461 up-regulated (FDR < 0.05, p-value < 0.05 and log2-fold-change ≥ 1) and 1076 down-regulated (FDR < 0.05 and log2-fold-change ≤ -1) genes. Number of downregulated genes was higher than upregulated ones, consistent with the association between H3K4 trimethylation and gene activity (Fig. 3B). GSEA revealed several negatively and positively enriched pathways with significant normalized enrichment scores (Fig. 3C, Supplementary Fig. 3A). Downregulated genes upon ASH2L depletion significantly overlapped with cell cycle and mitosis-associated gene sets, whereas upregulated genes enriched for metabolism, glycolysis, and hypoxia pathways. Expression of representative downregulated hits were validated by qRT-PCR (Supplementary Fig. 3B).

Fig. 3

ASH2L transcriptionally regulates cell cycle via direct promoter interactions with G2/M-checkpoint and E2F target gene sets. A Schematic of RNAseq and greenCUT&RUN experiments. B Volcano plot of RNAseq data showing differentially expressed genes (DEGs) in ASH2L depleted U373 cells compared to controls on the 14 day post-transduction. C Normalized enrichment score results of gene set enrichment analysis (GSEA) for all gene sets available from MSigDB v7.5. Top 5 biological processes enriched in downregulated or upregulated genes upon ASH2L knockout are shown. P values were calculated by hypergeometric test. D greenCUT&RUN analysis in U373 cells revealing genomic localization of ASH2L at differentially expressed gene promoters. E Venn diagram and volcano plot of overlapped RNAseq and greenCUT&RUN data to show proportion of differentially expressed genes bound by ASH2L in U373 cells. Genes bound by ASH2L and differentially expressed upon ASH2L depletion are now called as critical genes, denoted by blue (for downregulated) and yellow (for upregulated) dots on volcano plot. F Normalized enrichment score and FDR-qval results of GSEA on critical genes for all gene sets available from MSigDB v7.5. Some of the negatively enriched pathways related to cell cycle are highlighted. G Top 10 biological processes enriched in regulatory genes upon ASH2L knockout. H RNAseq heatmaps of critical genes. I greenCUT&RUN representative igv plots of critical genes. J Flow cytometric cell cycle analysis of ASH2L-depleted U373 by PI staining on post-transduction day 14 and its statistical analysis. P values determined by two-tailed Student’s t-test *P < 0.05, **P < 0.01, ***P < 0.001

Regulation of the transcriptome by ASH2L prompted us to analyze its localization on chromatin using CUT&RUN in glioblastoma cells. Using the default criteria of the HOMER (FDR < 0.0001 for calling the peaks and fold changes ≥ 4 and p value < 0.001 against control for filtering the peaks), 13,360 and 1964 chromatin peaks were called in U373 and U87MG cells, respectively (Fig. 3D). To find the commonalities between these two cell lines, initially all the peaks were merged, and coverage was compared. In total, 13,782 merged peaks were identified, of which, 117 (0.85% of total merged peaks, 5.96% of the U87MG peaks) and 4493 (32.6% of total merged peaks and 33.6% of the U373 peaks) were found unique for the respective cell lines. Conversely, 94.1% of U87MG peaks and 66.4% of U373 peaks were common indicating direct targets of ASH2L in glioblastoma cells (Supplementary Fig. 3C).

ASH2L is critical for cell cycle regulation

We explored genes bound by ASH2L on chromatin and differentially expressed upon ASH2L depletion as major targets of ASH2L. By overlapping RNAseq and greenCUT&RUN results, we generated a list of genes that may be directly linked with the essential function of ASH2L. When ASH2L peaks on promoter-TSSs were overlapped with DEGs from RNA-seq data, a total 150 of up-regulated and 350 down-regulated genes were identified, which together constitute a “critical gene list” (Fig. 3E). GSEA conducted on this list led to the identification of G2/M checkpoint, E2F targets and mitotic spindle hallmark gene sets as enriched (Fig. 3F-G). Heatmap for differential expression of cell cycle and mitosis related genes as well as greenCUT&RUN IGV tracks for selected genes are shown as examples (Fig. 3H-I, Supplementary Fig. 3D). Based on our finding of cell cycle regulatory genes as ASH2L targets, we conducted cell cycle analysis on ASH2L-depleted U373 and U87MG cells. Propidium iodide staining revealed decreased number of cells in the G0/G1 and accumulation at S and G2/M phases of cell cycle upon ASH2L knockout, resembling mitotic arrest (Fig. 3J, Supplementary Fig. 2D). Collectively, these results suggest that ASH2L plays a critical role in maintenance of glioblastoma cell fitness through regulating cell cycle progression.

Glioblastoma cells have differential dependency on SET1/MLL family of transcription factors interacting with ASH2L

To address the essential role of ASH2L for glioblastoma cell survival further, we focused our attention to chromatin modifying complexes, specifically, SET1/MLL family of histone methyltransferases SET1 (SETD1A), SET1B (SETD1B), MLL1 (KMT2A), MLL2 (KMT2B), MLL3 (KMT2C) and MLL4 (KMT2D) [28]. Each of these methyltransferases can form distinct multi-subunit complexes with diverse functions. All SET1/MLL families contain four core subunits, WDR5, RBBP5, ASH2L and DPY30, which are collectively named as the WRAD module. Other subunits of SET1/MLL family complexes are specific to one or a few SET1/MLL families and involved in chromatin recruitment [29]. Being part of WRAD module, ASH2L acts as a cofactor to support trimethylation [16,17,18,19]. To assess which SET1/MLL family complexes ASH2L is mainly part of in glioblastoma cells, we performed quantitative mass spectrometry (MS) of the ASH2L interactome for U373 and U87MG cell lines. Cells were transduced with inducible N-terminal GFP fusions of ASH2L gene encoding N-GFP-ASH2L fusion protein. Nuclear and cytoplasmic proteins were extracted, GFP-affinity purified and subjected to tandem MS (Fig. 4A). Accordingly, we identified common members of MLL-family and SETD1A/B complexes (DPY30, WDR5, ASH2L, RBBP5, HCFC1); as well as with MLL-specific (KMT2A, KMT2B, KMT2C, KMT2D, NCOA6, KDM6A, MEN1, PSIP1, PAGR1, PAXIP) and SETD1A/B-specific (SETD1A, SETD1B, WDR82, BOD1, CXXC1) complex members in nuclear fraction. Based on LFQ values, we determined the stoichiometry of interacting proteins (Fig. 4B-D). While in cytoplasmic fraction, interaction was evident only with SETD1A/B complex (Supplementary Fig. 4A, B); in nuclear fraction ASH2L interacted with MLL family and SETD1A/B members with similar stoichiometry (Fig. 4B-D). To assess the functional importance of SET1/MLL family complex members for glioblastoma cell viability, we performed siRNA mediated knock-down of WDR5 and KMT2A (MLL1) in U373 and U87MG cells and checked their colony forming ability (Fig. 4E). Downregulation of MLL1 significantly reduced colony forming potential of both U373 and U87MG cells, whereas WDR5 silencing affected only U87MG cells (Fig. 4F). Taken together, these results implicate ASH2L-MLL1 complex as a specific dependency in GBM cell proliferation.

Fig. 4

Glioblastoma cells have differential dependency to SET1/MLL family of transcription factors interacting with ASH2L. A Scheme of tandem mass spectrometry experiments performed on glioblastoma cells. B, C Proteomic analyses of GFP-tagged ASH2L in U373 (B) and U87MG (C) cells. Volcano plots of significant interactors of GFP-ASH2L isolated from nuclear extracts are shown. Stoichiometry plots of bound SET1/MLL family members (MLL family specific: yellow, SETD1A/B specific: red, MLL & SETD1A/B common members: blue) are depicted. All interactors are normalized to the GFP-ASH2L bait. Results shown represent Intensity Based Absolute Quantification with standard deviations. D Distribution of ASH2L common interactors in terms of subcellular localization in U373 and U87MG cells. E qRT-PCR analysis for mRNA levels upon siRNA mediated knockdown of WDR5 and KMT2A genes in U373 and U87MG cells. F Representative images of long-term clonogenic assay upon siRNA mediated knockdown of WDR5 and KMT2A genes and statistical analysis. Quantification of colonies was performed by ImageJ software. P values determined by two-tailed Student’s t-test *P < 0.05, **P < 0.01, ***P < 0.001

ASH2L is essential for tumor progression in vivo

To assess the role of ASH2L in more clinically-relevant settings, we tested the effect of ASH2L KO in patient derived primary glioblastoma cells, GBM4. GBM4 cells, upon transduction with g-ASH2L, exhibited efficient reduction of ASH2L expression (Fig. 5A), and displayed significantly lower proliferation and sphere formation ability (Fig. 5B, C). Encouraging results obtained in primary GBM cells led us to investigate the role of ASH2L in glioblastoma growth in vivo. We utilized orthotopic xenograft models using Luciferase (Fluc)-expressing U87MG cells, stably transduced with control sgRNA or ASH2L sgRNA. Following transduction, cells were intracranially injected and tumor growth was monitored until 33 days post-implantation (Fig. 5D). Accordingly, U87MG cells with ASH2L knockout could not form tumors as compared to the control cells (Fig. 5E). Efficient block of tumor forming capacity by lack of ASH2L encouraged us to check available patient data. Indeed, ASH2L expression was higher in glioblastoma in comparison to low grade gliomas (LGG-Grade 3 and Grade 3); ASH2L was expressed in all glioblastoma subtypes, but enriched for proneural subtype, based on TCGA (Fig. 5F-G). We next examined ASH2L protein in patient samples by immunohistochemistry utilizing a glioblastoma tissue microarray with 80 cores representing 40 different glioblastoma cases. 70% of glioblastoma specimens had strong ASH2L expression, whereas nonmalignant tissue specimens had none/low expression (Fig. 5H, I, Supplementary Fig. 5). Collectively, these clinically relevant results illustrate that ASH2L plays a critical role in glioblastoma progression.

Fig. 5

ASH2L is essential for glioblastoma tumor growth in vivo. A qRT-PCR analysis of ASH2L expression levels upon transduction of GBM4 primary glioblastoma cells with g-ASH2L or g-NT. B Representative images of GBM4 primary neurospheres 16 days post-transduction. C Cell viability analysis of GBM4 cells 1, 3 and 5 days after seeding. D Scheme for in vivo validation of ASH2L essentiality for GBM via intracranial injection of g-ASH2L or g-NT transduced U87MG cells (n = 5 per group) and bioluminescence detection. E Representative images of 3 mice taken during bioluminescence measurements at days 0-12-33 were illustrated. Bioluminescence signal of tumors formed by g-ASH2L or g-NT transduced U87MG cells were compared 33 days post-injection. F The boxplots displaying the ASH2L gene expression ratio (RSEM normalization values) in low grade gliomas (LGG, grade 2 and grade 3) and GBM based on TCGA. G The boxplots displaying the ASH2L gene expression ratio in GBM subtypes based on TCGA. H Representative core images from Brain Glioblastoma tissue microarray (TMA) stained with anti-ASH2L antibody. Scale bar 500 μm. I Percentage of ASH2L-positive cores were shown (n = 80 cores). P values determined by two-tailed Student’s t-test *P < 0.05, **P < 0.01, ***P < 0.001. J Model of ASH2L essentiality for glioblastoma cell survival. ASH2L together with DPY30, RBBP5 and WDR5 forms WRAD module, which acts as cofactor of SET1/MLL family transcription factors (MLL1/2, MLL3/4, SETD1A, SETD1B) to bind promoters of target genes and induce methylation of H3K4, a mark for euchromatin state. Active RNA polymerase binds to open chromatin to initiate transcription of genes involved in cell cycle progression, regulation of mitotic spindles and survival of cells (e.g. TRA2B, BARD1, KIF20B, ARID4A, SMARCC1). Deregulation of ASH2L levels interferes with cell cycle and leads to cell cycle arrest and apoptosis