The Meningioma Enhancer Landscape Delineates Novel Subgroups and Drives Druggable Dependencies

Meningiomas are the most common primary intracranial tumor, constituting more than one third of all primary central nervous system (CNS) neoplasms (1). Meningiomas are classified by grade according to histologic characteristics; although most have been considered benign (WHO grade I), a substantial percentage are higher grade (WHO grades II and III; refs. 1, 2). Furthermore, histologic grade fails to fully predict recurrence, and, upon failure of surgery and/or radiation, there are no effective systemic therapies for these patients.

Meningiomas exhibit a very strong and unexplained epidemiologic sex bias. Grade I tumors afflict female over male patients at a 3:1 ratio, a predilection that is lost in higher- grade tumors (1). Case reports and small studies have implicated progesterone agonists in meningioma growth, with reports of tumor regression following discontinuation of exogenous progesterone (3, 4). However, clinical trials have demonstrated inconsistent results, and predicting responders to hormone therapy remains a significant challenge (5–8). Identifying the basis for heterogeneity of hormonal responsiveness and the role of hormonal and sex-specific drivers in meningiomas will be critical to devising effective and targeted therapeutic strategies.

Most atypical (WHO grade II) and malignant (WHO grade III) meningiomas harbor NF2 mutations and monosomy at chromosome 22q, for which there are no clear therapeutic options (9). Prognostic methylation signatures implicate epigenetic dysregulation in meningioma, yet actionable therapeutic targets remain elusive (10, 11). The role of epigenetics in meningioma maintenance and progression is further highlighted by several reports of polycomb repressive complex (PRC) dysfunction, as well as global changes in methylation and recurrent mutations in epigenetic modifiers (12–17). Grade II meningiomas are reported to harbor PRC complex upregulation and histone 3 lysine 27 trimethylation (H3K27me3) hypermethylation compared with grade I tumors (14). However, global DNA hypomethylation (13) and reduced H3K27me3 (15) are associated with increased risk of recurrence. These data strongly implicate epigenetic dysregulation in meningioma growth and malignancy, but highlight the need for a more nuanced stratification that extends the current histologic classification scheme. Furthermore, the functional implication and therapeutic utility of these alterations remains unknown.

Enhancer profiling using histone 3 lysine 27 acetyl chromatin immunoprecipitation sequencing (H3K27ac ChIP-seq) can identify clinically relevant stratifications and targetable dependencies that may not be apparent from genetic and transcriptional methods alone (18–22). Enhancer signaling informs key mediators of transcriptional control and gene-expression programs. A subset of enhancers marked by particularly dense and long stretches of H3K27ac signal, designated as stretch enhancers or superenhancers (SE), highlight genes important to cell identity that require robust expression (20, 23, 24). Thus, in addition to informing transcriptional networks and altered chromatin regulation, enhancer profiles highlight novel dependencies that are not apparent from RNA-sequencing (RNA-seq) or genomic approaches.

We selected a cohort enriched for aggressive tumors spanning all three histologic grades, as these represent the meningiomas with a critical need for new medical therapies. By integrating the enhancer landscape of 33 meningiomas with clinical and transcriptional data, we derived clinically relevant subgroups demonstrating robust predictive value for recurrence irrespective of histologic grading schemes. Transcriptional regulatory networks reveal distinct lineage and hormonal drivers that appear to maintain subgroup-specific enhancer networks. Despite the predominance of NF2 mutations in our cohort, tumor-specific epigenetic and transcriptional signaling enrich for pathways regulated by other recurrent genetic mutations, implicating epigenetic processes in driving oncogenic programs independent of genomic aberrations. Finally, we elucidated druggable, pan-meningioma and tumor-specific dependencies by targeting SE-associated signatures to provide a foundation for developing new therapeutic approaches to meningioma.

We analyzed the enhancer landscape by H3K27ac ChIP-seq in an aggressive and diverse cohort of 33 meningiomas (Fig. 1A). Whole-exome sequencing (WES), RNA-seq, and DNA methylation analysis were performed for the majority of tumors (Supplementary Fig. S1). For a subset of samples, RNA-seq, WES, and DNA methylation data were already available (17). Tumors spanned all histologic grades and the cohort was evenly divided between male versus female patients, recurrent versus primary tumors, and history of radiation versus no adjuvant treatment. The median time to recurrence was 17 months and median overall survival was 35 months (Supplementary Table S1). To identify molecular features that stratify by prognosis and elucidate dependencies in the tumor subset with the greatest need for new therapeutic options, our cohort focused on an aggressive set of tumors. The majority of tumors harbored NF2 mutations and monosomy at chromosome 22. Fifty-five percent of tumors displayed loss of 1p36, a known risk factor for recurrence (Supplementary Fig. S1). DNA methylation analysis confirmed co-clustering of all except two tumors of our cohort with an independent cohort of reference meningiomas (Fig. 1B; ref. 25). Meningiomas are anatomically considered CNS tumors based upon location and are thought to derive from arachnoid granulation (AG) cells in the meninges. To determine whether meningiomas transcriptionally resembled CNS tumors or other cancers, we projected the meningioma samples onto a map generated from The Cancer Genome Atlas (TCGA) RNA-seq data, representative of more than 30 different tumor types, revealing transcriptional similarities between meningioma and mesothelioma (which is also driven by NF2 mutations), sarcoma, testicular germ cell tumors, and ovarian cancer (Fig. 1C). To investigate whether enhancers provide unique information about the relationship of meningiomas to normal cell types and tissues, we utilized the H3K27ac Roadmap dataset from ENCODE (26, 27) to generate normal sample clusters. Furthermore, we obtained enhancer RNA (eRNA) data extracted from more than 10,000 TCGA RNA-seq samples for cancer clustering (28). We then compared the enhancer profiles of meningiomas to the enhancers of normal cells and other cancer types (see Methods). On the basis of consensus clustering metrics, we identified 5 clusters in each dataset (Supplementary Fig. S2A). Each tumor or normal tissue type was composed of varying amounts of signal from each cluster. Brain tumors demonstrated high eRNA signal from one specific cluster, whereas another eRNA cluster was predominantly expressed in breast and prostate cancers, along with urothelial bladder carcinoma and diffuse large B-cell lymphoma (Fig. 1D). Relative to tumor models, meningiomas derived the most signal from the cluster of glioblastoma and low-grade glioma and the cluster containing breast and prostate tumors (Fig. 1D). In comparison with normal cells and tissues, meningiomas were most similar to stem cells and to ovarian and breast tissues (Supplementary Fig. S2B). Epigenetically, meningiomas resembled both hormonally driven cancers and primary malignant brain tumors, whereas the predominant transcriptional signal derived from other NF2-mutated tumors such as mesothelioma.

To broadly characterize the enhancer profile of meningiomas, we combined enhancers identified in individual tumors into a set of consensus meningioma enhancers. Enhancers were assigned to their putative target genes by coactivation analysis of enhancer and gene expression within the boundaries of transactivation domains, as described previously (29). Top SE-associated genes comprised multiple regulators of MAPK signaling, including two dual-specificity phosphatase (DUSP) family members, DUSP1 and DUSP5 (Fig. 1E). Gene set enrichment analysis (GSEA) of SE-associated genes identified enrichment for Hippo, WNT, Notch, and Sonic hedgehog signaling pathways, as well as MAPK and RHO signaling, repression of cell death, and hormonal signaling pathways (Fig. 1F). Activation of these pathways may represent an epigenetic “second-hit” in NF2-mutated meningiomas.

Top SEs were largely shared between different tumor grades, and DUSP1 was a top SE-associated gene in all three grades (Supplementary Fig. S3A). SEs enriched in each grade were associated with higher H3K27ac signal and with increased expression of the predicted target gene within their respective grade compared with the rest of the cohort (Supplementary Fig. S3B). Gene ontology (GO) enrichment of genes associated with the top 100 enriched SEs in each grade revealed shared signaling pathways, such as hormone-responsive programs (progesterone response, ovulation cycle, and mammary gland epithelial differentiation). Distinct pathways enriched in grade III tumors included neural programs and p38/MAPK signaling (Supplementary Fig. S3C).

To interrogate tumor-specific epigenetic and transcriptional pathways, we compared enhancer and transcriptional profiles of meningioma tissue to three normal AG cell models derived from postmortem dissections of the superior sagittal sinus. A total of 1,447 SEs were unique to tumors, while 147 were lost in comparison with normal AG cultures (Supplementary Fig. S4A). One of the most downregulated SEs in tumor samples is predicted to regulate PLCH2 (Supplementary Fig. S4A). PLCH2 resides on 1p36, a common site of copy-number loss in meningiomas. The PLCH2 SE was specifically downregulated in comparison with neighboring SEs on the same cytoband, and so may be a candidate tumor suppressor on 1p36. Gene set enrichment of top tumor-specific SE-associated genes revealed upregulation of pathways involved in epithelial-to-mesenchymal transition (EMT), Notch signaling, repression of cell differentiation, and receptor tyrosine kinase signaling. AG SEs were enriched for programs involved in mucin glycosylation, cell polarity, and vitamin metabolism (Supplementary Fig. S4B). Gene expression of predicted SE targets correlated with their enrichment in tumor versus normal samples (Supplementary Fig. S4C). Transcriptional differences recapitulated a subset of altered pathways driven by meningioma-enriched SEs including GTPase, Notch, MAPK, Hippo, and PIP signaling, although top differentially expressed genes were often distinct from the most differential SEs. Chromatin and lysine modification were specifically highlighted by RNA-seq and implicate global chromatin remodeling in meningioma pathology (Supplementary Fig. S4D). These data should be interpreted with the caveat that comparison of meningioma tissue to cultured cells may be influenced by culture-specific artifacts.

A direct comparison of motif enrichment, using normal AG cells as the background, identified progesterone receptor (PR) and glucocorticoid receptor as the most highly gained in meningiomas versus normal enhancers, whereas AP1, MED2, and YY2 were lost in tumors (Supplementary Fig. S5A). The strong tumor PR motif enrichment combined with the epidemiologic evidence implicating progesterone in meningioma pathogenesis prompted us to investigate whether PR plays a role in the tumor enhancer landscape. The PR signature was highly enriched for SE-associated genes, with nearly half of PR-regulated genes associated with SEs (hypergeometric P = 7.7e-25; Supplementary Fig. S5B and S5C). Pathways of RHO activity, phosphotidylinositol processing, and carbohydrate metabolism were also enriched, suggesting a broad role of PR-regulated programs in tumor signaling and metabolism (Supplementary Fig. S5D). To test whether progesterone stimulates growth in normal AG, AG cells were treated with progesterone, resulting in modest upregulation of cell proliferation, whereas treatment with the PR inhibitor mifepristone had opposing effects (Supplementary Fig. S5E).

We next compared irradiated versus unirradiated and recurrent versus primary tumors. Motif enrichment of differential enhancers (P < 0.05 and log₂ fold change >1) revealed that motifs targeted by the transcription factors ALX, LMX, and ARID were enriched in primary tumors, whereas EGR, TFAP, SP, and KLF family members were enriched in recurrent tumors (Supplementary Fig. S6A). KLF, SP, and AP1 motifs were enriched in irradiated samples, whereas GATA motifs were depleted (Supplementary Fig. S6B). We then analyzed pathways regulated by these differential enhancers. Recurrent tumors exhibit neuronal and glial signaling whereas primary tumors enrich for pathways involved in more mesodermal pathways (Supplementary Fig. S6C). Unirradiated tumors reflected similar signaling processes as primary tumors, potentially due to the moderate overlap of these two cohorts (Supplementary Fig. S6D). There were no significantly upregulated pathways in irradiated tumors. SEs associated with HOXB3 and SOX10 had greater signal in recurrent samples, whereas the known tumor suppressor MN1 was upregulated in unirradiated samples (Supplementary Fig. S6E and S6F).

To determine whether enhancers delineate clinically meaningful tumor subgroups, we performed non-negative matrix factorization (NMF) clustering on all 33 tumors using consensus meningioma enhancers. Clustering metrics revealed strong consensus at k = 3, with exclusion of one sample (Fig. 2A; Supplementary Fig. S7). These clusters stratified tumors into prognostic subgroups, with group 3 tumors exhibiting rapid recurrence (Fig. 2B). We integrated enhancer subgroup information with available clinical data and known predictors of recurrence, including 1p36 status, Heidelberg methylation cluster assignment, histologic grade, and history of radiation. Although some of the above factors were individual predictors of recurrence-free survival (Supplementary Table S2), in a multivariable Cox proportional hazards model, only enhancer cluster and history of radiation were significantly predictive when combined into a multivariable analysis (Supplementary Table S3). Nearly all tumors were identified as meningiomas by the Heidelberg classifier, but only approximately one third received a clear assignment into methylation subclasses. Tumors with the worst prognoses were more frequently unclassified, suggesting that expansion of the meningioma methylation classifier may be warranted, in particular with respect to rapidly recurrent tumors (Fig. 2C; ref. 11). To determine whether methylation state informed prognosis independent of existing classifiers, we performed NMF clustering on our samples using the 10% most variable probes from the 450K methylation array. K = 5 was selected due to the relatively higher cophenetic and silhouette scores (Supplementary Fig. S8A and S8B). Consensus clustering assignment was used to segregate subgroups (Supplementary Fig. S9A). The methylation subgroups were plotted against enhancer classifications in a tanglegram, with lines connecting identical samples (Supplementary Fig. S9B). In our cohort, enhancer clustering provided a cleaner segregation and a prognostic advantage over methylation clustering (Supplementary Table S2), although one small subgroup of methylation clustering with 3 samples had a poorer prognosis.

Several individual SEs were also predictive of recurrence based upon their presence or absence in a sample (Supplementary Table S4). The strongest predictor of rapid recurrence was an SE associated with CDH5 [HR: 41.9; 95% confidence interval (CI): 4.6–378.5], a biomarker of metastatic breast cancer and high-grade glioma (30, 31). A nearby SE at LRRC36 was the only positive prognostic predictor (HR: 0.1; 95% CI: 0.05–0.4; Supplementary Fig. S10). Clinical characteristics of enhancer subgroups revealed that group 3 tumors were enriched for, but not exclusive to, grade III tumors and male patients, whereas group 1 tumors had a greater proportion of grade II tumors, and group 2 tumors had a majority of grade I tumors. Both group 1 and group 2 patients had a greater proportion of female patients (Fig. 2D). History of radiation and primary versus recurrent status of the tumor did not segregate between groups (Supplementary Fig. S11).

To validate our prognostic subgroups, we performed H3K27ac ChIP-seq on an additional 14 meningiomas (Supplementary Table S5). Tumors were subgrouped with single-sample GSEA (ssGSEA) using the top 100 differential enhancers from each subgroup from the discovery cohort as a signature. Unsupervised hierarchical clustering based on the scaled ssGSEA score, in combination with the maximal scaled subgroup score, was used to classify the validation cohort tumors. This method assigned the correct subgroup to the original cohort in 31 of 32 tumors. In every case, the subgroup assignment using the maximal ssGSEA scores was concordant between tumors on a given leaf (Fig. 2E). We then assessed whether this classification was prognostic in the validation cohort (Fig. 2F-I; Supplementary Table S6). Tumors in group 3 were more likely to recur irrespective of grade (P = 0.1; Fig. 2F and G). Within grade II, group 3 tumors were significantly more likely to recur (P = 0.0025; Fig. 2H and I).

We also confirmed that our subgroup classifier was relevant in a published independent cohort of tumors using gene expression data (17). We used the top 100 or 250 group 3 SE-associated genes to form a gene signature and stratified meningiomas by high versus low expression, cut at the median (Supplementary Fig. S12A and S12B). Tumors with high expression of these signatures tended to have a higher recurrence rate (top 250 SEs: P = 0.138; top 100 SEs: P = 0.161).

To determine whether enhancer subgroups harbor distinct molecular drivers, we compared typical enhancer, SE, and transcriptional profiles. Although the top SEs were generally shared across subgroups (Fig. 3A), subgroup-enriched SEs highlighted specific programs and differential SEs within in each group (Fig. 3B; Supplementary Fig. S13A). Gene expression changes recapitulated a subset of SE signaling, although each data type presented additional unique pathway enrichments. Genes involved in cholesterol metabolism and repression of stem cell differentiation were notable within group 1. We therefore designated this group to be “adipogenesis/cholesterol.” Mesodermal and muscle development pathways were upregulated in group 2, whereas neuronal programs were enriched in group 3 tumors (Fig. 3B and C). Two differentially expressed genes, CRYGN and SPOCK1, effectively stratified subgroups 1 and 2 versus 3. CRYGN was upregulated in subgroups 1 and 2 versus 3 (ANOVA P = 0.00015, Student t test: 1 vs. 3 P = 1.7e-4, 2 vs. 3 P = 7e-4; Supplementary Fig. S13B and S13C). SPOCK1 marked subgroup 3 versus 1 and 2 (ANOVA P = 1.15e-5, Student t test: 1 vs. 3 P = 2e-5, 2 vs. 3 P = 4.2e-4; Supplementary Fig. S13D). Combinatorial expression of the PR network gene CRYGN versus SPOCK1 effectively stratified the poorly prognostic subgroup 3 from subgroups 1 and 2 and may serve as a potential future clinically implementable assay to identify tumors with poor prognosis (Supplementary Fig. S13E).

There was a notable skewing of H3K27ac signal toward group 1 tumors, which was particularly evident in SEs (Supplementary Fig. S13A). Given the reported alterations in PRC complex members, we investigated whether gene expression changes or mutations could explain this bias. There was no enrichment of PRC complex mutations in the group, nor were there clear alterations of gene expression in these pathways. However, subgroup SEs enriched for targets of SUZ12, a member of the PRC complex (Supplementary Fig. S13F). In addition, 3 of 4 tumors with the highest number of SEs had mutations in the deacetylase SIRT2, predicted by 3 different algorithms to be deleterious (Supplementary Fig. S13G).

As enhancers are associated with transcription factor (TF) binding, we reasoned that enhancers with similar transcriptional drivers could be segregated into distinct modules based upon correlation of H3K27ac signal. Enhancer modules empower reconstruction of transcriptional circuitry by integrating TF expression data and enhancer motif enrichment for different modules. Degree of connectivity in both protein and transcriptional networks predicts critical genes in cell survival and disease prognosis (32–34). Thus, enhancer hubs within a given module may highlight enhancers and genes of particular importance. Module eigen-enhancers compared across clinical or biological characteristics identify modules and predict transcriptional programs that are differentially enriched within a given trait. We therefore generated a weighted correlation network of enhancers from the meningioma enhancer peakset (Fig. 4A; refs. 35, 36). Enhancer networks exhibited clear scale-free topology and segregated cleanly into modules (Supplementary Fig. S14A and S14B). SEs were more highly connected than typical enhancers, both within and between modules (Supplementary Fig. S14C). SEs may therefore be of particular importance in maintaining the structure of epigenetic networks. Using the corresponding eigen-enhancer signal, we identified modules that were significantly altered between enhancer subgroup, sex, recurrence-free survival, or Heidelberg methylation classification (Supplementary Fig. S14D). Enhancer subgroup and methylation classification were associated with the most differential module signal throughout the network relative to other features. There was a significant overlap between modules that were differential by enhancer subgroup and those that were distinct by methylation group, highlighting the interplay between these epigenetic regulatory mechanisms in network structure (hypergeometric P = 0.005; Supplementary Fig. S14E).

We next explored individual modules that were differentially enriched by tumor recurrence, 1p36 status, or patient sex for further interrogation (Supplementary Fig. S14D). The subnetwork that was associated with rapid recurrence was highly enriched for regulators of neural crest differentiation and stem cell maintenance. HOXD family members and FOXM1, which were previously associated with malignant behavior (17), were also enriched in this module (Supplementary Fig. S15A). The 1p36-associated module was upregulated in tumors harboring 1p36 deletion. Motifs of two TFs on the 1p36 cytoband were highly enriched in this enhancer module. These TFs, HES5, a NOTCH-responsive repressive factor that regulates brain and cardiac development, and PAX7, which regulates myogenesis, may therefore be important transcriptional regulators associated with loss of 1p36 in meningioma (Supplementary Fig. S15B). Three sex-differential modules were enriched in female versus male patients (Supplementary Fig. S15C). TFs predicted to regulate these modules were highly enriched for adipogenesis and circadian rhythm regulation. Furthermore, these TFs were upregulated in Genotype–Tissue Expression (GTEx) adipose and breast tissue and downregulated in GTEx brain tissue. These subnetworks revealed distinct transcriptional programs in male versus female tumors.

To predict which TFs were driving subgroup differences, we generated subgroup-specific networks by selecting all modules that demonstrated subgroup-specific changes (ANOVA P < 0.1) and assigning them to the subgroup with the highest eigen-enhancer enrichment (Fig. 4B; Supplementary Fig. S14D). Correlating TF expression with the module eigen-enhancer for each subgroup-specific module restratified enhancer subgroups despite being unsupervised, supporting a functional interaction between TF enrichments and enhancer subnetworks (Supplementary Fig. S16). We intersected this list with TF motifs enriched in each subgroup-specific module. The top ten TFs ranked by motif enrichment were used to build subgroup TF networks (Fig. 4C). Group 1 TFs were enriched for programs involved in cholesterol and adipogenesis as well as hormonal nuclear receptors (Fig. 4D). Group 2 TFs were enriched in mesodermal development, consistent with differential SE (Fig. 3B) and RNA-seq (Fig. 3C) programs (Fig. 4D). Group 3 TFs regulated neural crest development (Fig. 4D).

Given the overrepresentation of hormone receptor TFs in multiple subgroups and the reported intertumoral heterogeneity of PR, we investigated whether hormone receptor expression or activity distinguished meningioma subgroups. PR signature scores were significantly higher in the adipogenesis/cholesterol and mesodermal group (Fig. 4E). Conversely, androgen receptor (AR) signature was upregulated in the neural crest group (Fig. 4E). However, individual hormone receptors including AR (Supplementary Fig. S17A), estrogen receptor (ER; Supplementary Fig. S17B), and PR (Supplementary Fig. S17C) demonstrated trends in expression, but did not differ consistently and significantly between groups. The difference between AR and PR expression separated subgroup 2 versus 3 (P < 0.05) (Supplementary Fig. S17D). PR and the PR-regulated gene expression were positively and significantly correlated (r = 0.64, P = 0.0013; Supplementary Fig. S17E). Hormone receptor activity, rather than expression alone, delineated subgroups. Thus, subgroup-enriched TF networks implicate distinct transcriptional programs and highlight hormonal players in subgroup-specific tumor maintenance.

Given the lack of effective medical therapies for meningioma and our data demonstrating that SEs mark important genes for tumor cell survival, we sought to identify novel, druggable targets using SE-associated genes. First, we assessed the top consensus SEs to identify individual druggable genes. Among the top 10 SE-associated genes, only DUSP1 had an available inhibitor (Fig. 5A; ref. 37). We also derived signatures from the top SE-associated genes in each subgroup and predicted drugs to antagonize these signatures (see Methods; Fig. 5A). We tested this panel of drugs against 6 meningioma cell lines and 2 normal AG cell lines to identify drugs with antitumor activity (Supplementary Table S7). Two existing meningioma models (CH157-MN and IOMM-Lee) were used along with 4 additional models derived in this study. Cells were maintained in DMEM with 7% FBS. Of the 6 tumor cell lines tested, 3 were aggressive models to ensure drug efficacy in the tumors in greatest need of new therapeutics: DI-134, derived in-house from the poor prognosis/neural crest subgroup; CH157-MN, an existing grade III model; and IOMM-Lee, a model of unknown grade harboring genomic instability. Drug responses were benchmarked against the FAK inhibitor GSK2256098, which is currently in clinical trials to target the NF2-mutated subset of meningiomas based on its efficacy in vitro in inhibiting growth of merlin-deficient mesothelioma (https://clinicaltrials.gov/ct2/show/NCT02523014).

Of the 17-drug panel, 7 drugs had an average maximum effect (E_max) of >50% reduction in cell viability versus DMSO control (Fig. 5B). The DUSP1/6 inhibitor BCI had the strongest effect on cell viability with an E_max < 0.25 in all 6 cell lines. Additional effective drugs included inhibitors of bromodomain (JQ1), MDM (RITA), aurora kinase (AT2983), FGFR (orantinib), HSPs (KRIBB11), HIV protease (ritonavir), and RHO kinase (GSK269962A). FAK inhibition was minimally effective in any cell model with an average E_max > 0.5 (Fig. 5C).

The SE assigned to DUSP1 was strong and present throughout the cohort, making it a rational target for pan-meningioma therapy (Supplementary Fig. S18A). We validated the screen results by treating 5 meningioma cultures and a normal AG culture with the same inhibitor (Fig. 6A). The BET inhibitor JQ1 reduced expression of DUSP1, confirming DUSP1 regulation by an SE (Fig. 6B). To confirm DUSP1 as a dependency in meningioma, we targeted DUSP1 by shRNA in the meningioma cell line IOMM-Lee (Fig. 6C and D) and the newly derived meningioma model DI-134 (Fig. 6E and F). Depletion of DUSP1 with two nonoverlapping shRNAs reduced cell viability versus a nontargeting control (Fig. 6C–F). Knockout of DUSP1 via CRISPR/Cas9 in CH157-MN (Supplementary Fig. S18B and S18C) or IOMM-Lee (Supplementary Fig. S18D and S18E) impaired cell viability in vitro. Consistent with the reported phosphatase function of DUSP1, treatment with BCI increased phosphorylation of ERK, JNK, and p38, and induced expression of cleaved caspase-3 and cleaved PARP, consistent with induction of apoptosis via known targets of DUSP1 (Fig. 6G). Depletion of DUSP1 in an orthotopic xenograft model using CH157-MN cells prolonged mouse survival versus nontargeting control (Fig. 6H and I). To test the efficacy of DUSP1/6 inhibition in vivo in the absence of restricted delivery, CH157-MN cells were implanted subcutaneously into flanks of immunodeficient mice. After 1 week of tumor growth, mice were treated twice daily with 10 mg/kg of BCI or an equivalent volume of DMSO via intraperitoneal injection. No toxicity was observed over the course of treatment (Supplementary Fig. S18F). Treatment with BCI decreased tumor weight (P = 0.043; Fig. 6J) and tumor volume (P = 0.029; Fig. 6K and I) compared with DMSO. Using enhancer profiling, we identified a set of novel drug candidates in meningioma and demonstrated the therapeutic potential of a DUSP1/6 inhibitor in vivo. DUSP1/6 inhibition is a promising therapeutic candidate that can be targeted for further compound development.

Although the majority of meningiomas are histologically benign, these tumors cause significant morbidity and mortality in many patients. Meningiomas that fail standard-of-care therapy of surgical resection and radiation lack effective medical therapies. Furthermore, clinicians face challenges selecting which patients should undergo radiotherapy after surgery, potentially exposing patients to unnecessary and permanent toxicity, based upon histologic grading alone. Methylation clustering provides additional prognostic stratification, but such studies do not provide actionable targets (11). Mutational profiling is one option for stratifying tumors, largely low-grade ones (12, 38). However, although mutations may serve as initiating events, the loss of a tumor suppressor, such as NF2, does not provide an actionable target, nor sufficiently account for oncogenic programs. We now propose an additional strategy to delineate meningiomas by risk of recurrence based upon novel prognostic subgroups of meningioma driven by distinct transcriptional drivers (Fig. 7). This stratification effectively discriminated high- versus low-risk subgroups in a discovery and validation ChIP-seq cohort. A limitation of this study was the modest sample size, and thus future work will be necessary to further validate our findings regarding the prognostic potential of the identified subgroups. High-risk signatures applied to gene expression data trended toward significance for predicting rapid recurrence, suggesting that future work could investigate improved methods to extend the enhancer classification to apply to additional data types. Enhancer subgroups exhibit signatures of hormonal drivers and may be amenable to hormonal therapies using progesterone antagonists. Across tumor grade and subgroup, SE signatures reveal novel druggable targets, which can inform future therapeutic development for tumors that have not responded to resection and radiation.

Meningiomas are anatomically classified as CNS tumors, but contextualizing them in the larger cancer landscape allows for improved molecular understanding and novel therapeutic strategies. Transcriptional and epigenetic data provide complementary but unique insights into tumor biology. Although transcriptional analyses reveal that meningiomas most closely resemble other NF2-mutated tumors, the enhancer profile of meningiomas demonstrates characteristics of both CNS malignancies and hormonally driven tumors. Interrogation of the SE landscape implicates PR as a core transcriptional driver of a subset of tumors. These data provide a putative explanation for the strong epidemiologic sex bias and case reports of meningiomas receding following discontinuation of progesterone agonists (3–7). The uneven responses reported in clinical trials may thus be driven by distinct epigenetic subgroups of tumors. Further investigation into the specific biology of arachnoid cap cells that promote their responsiveness to progesterone is warranted.

Epigenetically dysregulated programs recapitulated pathways regulated by known meningioma driver mutations. Our cohort, which consisted largely of NF2-deficient tumors, demonstrated strong upregulation of not only Hippo, but also Notch, WNT, and Sonic hedgehog signaling pathways. Dysregulation at the epigenetic level may reveal additional oncogenic drivers following NF2 loss in merlin-deficient tumors. SE pathways also highlight altered MAPK signaling and coordinated changes in the PIP pathway enzymes relative to normal AG cells.

By deriving enhancer networks, we identified epigenetic and transcriptional modules that correlate with a wide range of important clinical and genetic features of our cohort. Rapid recurrence was highly associated with neural crest progenitor and stem cell programs, as well as neuronal signaling pathways. The meninges are thought to derive from both mesenchymal and neural crest lineages, and thus malignant tumors may reactivate developmental programs distinct from less aggressive, more mesenchymal tumors, with this distinction reflected in our reported enhancer-derived subgroups.

Finally, targeted therapies for meningiomas are lacking. There is a critical lack of drug options for patients with aggressive tumors for whom resection and radiation are insufficient. We demonstrate the utility of enhancer profiling to identify potent drug candidates for meningioma treatment. SE signatures effectively elucidate meningioma dependencies and predict drugs that potently and selectively kill meningioma cells in vitro and in vivo, paving the way for further investigation into novel therapeutic options in treatment of meningioma.

Meningioma samples were acquired from excess surgical resection tissues from patients at the Cleveland Clinic. All specimens were reviewed by a neuropathologist upon resection prior to banking and again prior to sample processing. Appropriate written informed consent was obtained from patients in accordance with a Cleveland Clinic Institutional Review Board (IRB)–approved protocol (090401), and the study was performed in accordance with the Declaration of Helsinki. Clinical data were collected under IRB 16956 from Cleveland Clinic. All studies involving human patients were conducted in accordance with the Declaration of Helsinki. Clinical information was not available for deidentified samples. Tissue was snap-frozen and stored at −80°C. For cell models, tissues were dissociated using the Papain dissociation system (Worthington Biomedical Corp; LK003150) and cultured in DMEM supplemented with 1% penicillin/streptomycin and 7% FBS. Short tandem repeat analyses were performed annually to authenticate the identity of each tumor model used in this article. Cells were stored at −160°C when not being actively cultured. All cells were thawed within 1 month of these experimental procedures. All experiments conformed to the relevant regulatory standards. CH157-MN was provided by the Gillespie lab at the University of Alabama-Birmingham and IOMM-Lee was provided by the Jensen lab at the University of Utah.

Arachnoid cells were derived from postmortem tissue via the body donation program at the Cleveland Clinic. Arachnoid granulations from the superior sagittal sinus were dissected and dissociated using the Papain dissociation system (Worthington Biomedical Corp.; LK003150) and cultured in DMEM with 1% penicillin/streptinomcyin and 10% FBS. Cells were immortalized at passage 6 by transfection with SV40-Large T antigen (Addgene plasmid #21826).

TRIzol reagent (Sigma-Aldrich) was used to isolate total cellular RNA from cell pellets. The qScript cDNA Synthesis Kit (Quanta BioSciences)was used for reverse transcription into cDNA. Quantitative real-time PCR was performed by using Applied Biosystems 7900HT cycler using SYBR-Green PCR Master Mix (Thermo Fisher Scientific). Quantitative PCR (qPCR) primers used in this study were human DUSP1-fwd: 5′-ACCACCACCGTGTTCAACTTC-3′, DUSP1-rev: 5′-TGGGAGAGGTCGTAATGGGG-3′, beta-2 microglobulin-fwd: 5′-GAGGCTATCCAGCGTACTCCA-3′ and beta-2 microglobulin-rev: 5′-CGGCAGGCATACTCATCTTTT-3′, GAPDH-fwd: 5′-TGACAACTTTGGTATCGTGGAAGG-3′, GAPDH-rev: 5′-AGGCAGGGATGATGTTCTGGAGAG-3′

Cells were collected and lysed in RIPA buffer (50 mmol/L Tris-HCl, pH 7.5; 150 mmol/L NaCl; 0.5% NP-40; 50 mmol/L NaF with protease inhibitors) and incubated on ice for 30 minutes. Lysates were centrifuged at 4°C for 15 minutes at 14,000 rpm, and supernatant was collected. The Pierce BCA Protein Assay Kit (Thermo Fisher Scientific, catalog no. 23225) was utilized to determine protein concentration. Equal amounts of protein samples were mixed with SDS Laemmli loading buffer, boiled for 10 minutes, and electrophoresed using NuPAGE Bis-Tris Gels, then transferred onto polyvinylidene difluoride membranes. TBS-T supplemented with 5% nonfat dry milk was used for blocking for a period of 1 hour followed by blotting with primary antibodies at 4°C for 16 hours. Blots were washed 3 times for 5 minutes each with TBS-T and then incubated with appropriate secondary antibodies in 5% nonfat milk in TBS-T for 1 hour. For all Western immunoblot experiments, blots were imaged using Bio-Rad Image Lab software and subsequently processed using Adobe Illustrator to create the figures. The following antibodies were used for Western blot analysis: p44/42 MAPK (Erk1/2), 1:1,000 (Cell Signaling Technology, catalog no. 4695, 137F5); phospho-p44/42 MAPK (Erk1/2; Thr202/Tyr204), 1:1,000 (Cell Signaling Technology, catalog no. 4370, D13.14.4E); p38 MAPK, 1:1,000 (ProteinTech, 14064-1-AP); phospho-p38 (Thr180/Tyr182), 1:2,000 (Cell Signaling Technology, catalog no. 4511, D3F9); SAPK/JNK, 1:1,000 (Cell Signaling Technology, catalog no. 9252); phospho-SAPK/JNK (Thr183/Tyr185), 1:1,000 (Cell Signaling Technology, catalog no. 4668T, 81E11); cleaved caspase-3 (Asp175), 1:1,000 (Cell Signaling Technology, catalog no. 9664S, 5A1E); HRP-conjugated GAPDH, 1:40,000 (ProteinTech, catalog no. HRP-60004); and DUSP1, 1:1,000 (LifeSpan BioSciences, catalog no. LS-C332288).

For in vitro cell viability assays, 1,000 cells/well were plated in a 96-well plate in DMEM with 7% FBS on day 0, then treated 24 hours later with the corresponding drug or DMSO at an equivalent percent volume to the highest drug concentration. Progesterone (Sigma catalog no. P8783-5G, CAS #57-83-0) was resuspended in 100% ethanol. Mifepristone/RU486 (Thermo Fisher Scientific, catalog no. 501011574, CAS #84371-65-3) was resuspended in DMSO. Cell viability was assayed 48 hours later using CellTiter-Glo reagent (Promega; G7572). Plates were shaken at room temperature for 15 minutes and then read for luminescence output.

For the compound drug screen assay, compounds and DMSO vehicle controls were transferred in a volume of 10 nL to black, clear-bottom 1,536-well plates (Greiner; 789092) using an acoustic transfer system Gen4+ instrument (EDC Biosystems). Cells were dispensed in a volume of 10 μL and a density of 500 cells per well using a Multidrop Combi liquid handler (Thermo Fisher Scientific; 5840300). Plates were incubated at 37°C and 5% CO₂ for 72 hours, at which point CellTiter-Glo reagent (Promega; G7572) was added to each well using the Multidrop Combi instrument. Plates were shaken at room temperature for 15 minutes and then read for luminescence output using an EnVision plate reader (Perkin Elmer). Maximum effect was then calculated as the greatest reduction in cell viability versus the average of the DMSO controls for the corresponding plate for each cell line. The following drugs were tested: GNE-0877 (catalog no. S7367, CAS #1374828- 69-9, SelleckChem), orantinib (TSU-68, SU6668; catalog no. S1470, CAS #252916- 29-3, SelleckChem), (+)-JQ1 (catalog no. S7110, CAS #1268524-70-4, SelleckChem), SB525334 (catalog no. S1476, CAS #356559-20-1, SelleckChem), ritonavir (catalog no. S1185, CAS #155213-67-5, SelleckChem), ropivacaine Hcl (catalog no. S4058, CAS #132112-35-7, SelleckChem), (+)-bicuculline (catalog no. S7071, CAS #485-49-4, SelleckChem), fenofibrate (catalog no. S1794, CAS #49562-28-9, SelleckChem), gatifloxacin (catalog no. S1340, CAS #112811-59-3, SelleckChem), piperacillin sodium (catalog no. S4222, CAS #59703- 84-3, SelleckChem), dual specificity protein phosphatase 1/6 Inhibitor, BCI (catalog no. 317496, CAS #15982-84-0, Millipore Sigma), RITA (catalog no. S2781, CAS #213261-59-7, Selleck-Chem), KRIBB11 (catalog no. S8402, CAS #342639-96-7, SelleckChem), AT9283 (catalog no. S1134, CAS #896466-04-9, SelleckChem), GSK2699 [catalog no. S7687, CAS #850664-21-0 (free base), SelleckChem], reboxetine mesylate (catalog no. S3199, CAS #98769-84-7, SelleckChem), and GSK2256098 (catalog no. S8523, CAS #1224887-10-8, SelleckChem).

Lentiviral constructs expressing nonoverlapping shRNAs directed against DUSP1 (TRCN0000355637, TRCN0000367631) or a nontargeting control shRNA (TRCN0000231489) with no targets in the human genome were obtained from Sigma-Aldrich. DUSP1-targeting shRNAs were assayed for knockdown efficiency by qPCR and were then used for all following experiments. For CRISPR/Cas9 experiments, sgRNAs were cloned into the lentiCRISPR v2 plasmid (Addgene #52961). DUSP1 sgRNA: 5′-CGTCCAGCAACACCACGGCG-3′. 293T cells were used to generate lentiviral particles by cotransfecting the packaging vectors psPAX2 and pMD2.G using a standard calcium phosphate transfection method in DMEM containing 1% penicillin/streptomycin. For cell viability assays, meningioma models were transduced with the corresponding lentiviral constructs and selected 48 hours later using 1 μg/mL of puromycin. Cells were selected for 72 hours, then plated at a density of 1,000 cells/well in a 96-well format with 12 wells per condition. Cell viability was assayed by incubating with CellTiter-Glo while shaking for 20 minutes and then read for luminescence.

In vivo drug studies were performed by implanting 1 million cells into the flank of NSG (NOD.Cg-Prkdcscid Il2rgtm1Wjl/SzJ; The Jackson Laboratory) mice. Once tumors had formed mice were randomly assigned into drug versus treatment group by a blinded investigator. Mice were treated with BCI (catalog no. 317496, CAS #15982-84-0, Millipore Sigma) or an equivalent percentage of DMSO resuspended in 0.9% saline by intraperitoneal injection daily. All murine experiments were performed under an animal protocol approved by the University of California, San Diego Institutional Animal Care and Use Committee. Healthy, wild-type male mice of NSG background, 4 to 6 weeks old, were randomly selected and used in this study for flank injection. For intracranial xenograft studies, mice were orthotopically implanted with 10,000 cells and were sacrificed upon observation of neurologic signs by a blinded investigator. None of the mice had experienced any treatment or procedures before the experiments described. Mice were housed together in a controlled environment with 14 hours of light and 10 hours of dark per day. Animal husbandry staff at the University of California, San Diego regularly observed all animals, and no more than 5 mice were housed in each cage. Animals were monitored until tumor size exceeded 2.0 cm in the longest dimension.

RNA-seq was performed as described previously (22). Total RNA was extracted from flash-frozen, pulverized tissue using the miRNeasy Mini Kit (Qiagen, catalog no. 217004) in accordance with the manufacturer's protocols. Libraries were prepared and sequenced by Genewiz. Stranded RNA library preparation was performed with ribosomal RNA depletion according to instructions from the manufacturer (Epicentre). Paired-end sequencing was performed on the Illumina HiSeq 2500 with 2 × 150 bp paired-end read configuration.

WES was performed as described previously (21). DNA from flash-frozen tumor samples was extracted using the Quick-DNA 96 Plus kit (Zymo). The Agilent SureSelect Human All Exon 50-Mb Target Enrichment Kit v4 was used to capture all human exons for deep sequencing using the manufacturer's protocol v2.0.1. The SureSelect Human All Exon Kit targets regions of 50 Mb in total size, which is approximately 1.7% of the human genome. In brief, 3 μg genomic DNA was sheared with a Covaris S2 to a mean size of 150 bp. Five hundred nanograms of library DNA was hybridized for 24 hours at 65 °C with the SureSelect baits. Fifteen nanograms final exome-enriched library (without barcode) was used as a template in a 50-μL PCR reaction. The Herculase II Fusion enzyme (Agilent) was used together with the NEBNext Universal PCR primer for Illumina and NEBNext Index primer (NEB #E7335S) under the following conditions: the initial denaturation step for 2 minutes at 98°C was followed by four cycles of 30 seconds at 98°C, 30 seconds at 57°C, 1 minute at 72°C, and a final step of 10 minutes at 72°C. Barcoded samples were then sequenced on the HiSeq2000 in 2 × 100-bp paired-end mode.

ChIP was performed as described previously (22). Briefly, 5 to 10 mg of flash-frozen primary meningioma tumors was pulverized, cross-linked using 1.5% formaldehyde, and sonicated to fragment sizes of 200 to 800 bp. Samples were incubated overnight with 5 μg H3K27ac antibody per ChIP experiment (Active Motif; 39133). Enriched DNA was quantified using PicoGreen (Invitrogen), and ChIP libraries were amplified and barcoded using the Thruplex DNA-seq Library Preparation Kit (Rubicon Genomics) according to the manufacturer's recommendations. Following library amplification, DNA fragments were agarose gel (1.0%) size selected (<1 kb), analyzed using a Bioanalyzer (Agilent Technologies), and sequenced at Genewiz using Illumina HiSeq 2500 2 × 150-bp paired-end reads.

For the subset of samples for which previously published RNA-seq data were available, data were downloaded from the Gene Expression Omnibus (GEO; GSE101638). To match the data format of this series, only single-end reads were used and were trimmed to 75 bp using BBDuk from the BBMap toolset. Trim Galore was used to trim adaptors and remove low-quality reads (39). Reads were quantified against Gencode v29 using Salmon with correction for fragment-level GC bias, positional bias, and sequence-specific biases (40). Transcripts were summarized to gene level and processed to transcripts per million (TPM) using the R/Bioconductor (http://www.R-project.org/) package DESeq2 (41) and batch-corrected using ComBat from the R sva package (42). Comparisons were performed using contrasts in DESeq2 followed by Benjamini–Hochberg adjustment to correct for false discovery rate (FDR). For GSEA comparing tumor versus normal, a preranked list was generated and weight by the inverse of the FDR multipled by the sign of the log₂ fold change. This list was used as input to the desktop GSEA software (43). Gene sets tested were GO biological processes, Reactome, and KEGG. Enrichment maps for all gene sets significant at an FDR < 0.2 with an edge cutoff (gene set similarity) of 0.375 were visualized using the Bader Lab Enrichment Map plugin (44) in Cytoscape v3.6 (45).

For the subset of samples for which previously published whole-exome sequencing data were available, data were downloaded from GEO (GSE101638). Low-quality reads and adaptors were trimmed using Trim Galore (39). BWA-MEM version 0.7.17 was used to align paired-end exome sequencing reads to the hg19 reference genome (46). SAMtools (47) was used to sort and index BAMS. PCR duplicates were removed with PicardTools (http://broadinstitute.github.io/picard/). Single-nucleotide variants (SNV) and insertions/deletions (indels) were identified with the Genome Analysis Toolkit v3.8 (48) in accordance with the Genome Analysis Toolkit best practices with the principal steps of base quality score recalibration, variant genotyping for single-nucleotide variants and indels, and variant hard-filtering with standard recommendations (49, 50). Variants with predicted deleterious effects were annotated using ANNOVAR (51). SNVs or indels predicted to be deleterious by 2/3 of PolyPhen2, whole-exome SIFT, or MutationTaster with a frequency of <1% in the 1000 Genomes Project, NHLBI-ESP, and CompleteGenomes were annotated as mutations.

Genomic DNA was extracted as described for WES and samples were sent for Infinium MethylationEPIC BeadChip Kit (Illumina). Probes from the Infinium MethylationEPIC BeadChip were downsampled to only those intersecting with the 450K chip. These probes were processed with the R packages minfi (52) and CopywriteR (53), which were used to identify copy-number variations for NF2 and the corresponding regions of chromosome 22 and for 1p36.

Methylation data was uploaded to the Heidelberg meningioma classification tool (25) at https://www.molecularneuropathology.org/mnp. Meningioma methylation group was assigned to the classification with the highest clustering score. Meningiomas not matching any classification (score < 0.3), were considered unclassified.

Normalized methylation data for CNS tumors was downloaded from GSE109379, and meningioma cohort data were processed using the same pipeline as described previously (25). Data from GSE109379 were clustered with the UMAP algorithm using the R package uwot (https://github.com/jlmelville/uwot) using the parameters: n_neighbors = 100, learning_rate = 0.5, init = “random,” pca = 50, min_dist = 0.4. New samples were embedded onto the map.

H3K27ac ChIP-seq was processed following the ENCODE guidelines. FASTQ reads were trimmed to remove low-quality reads and adaptors with Trim Galore (39) and uniquely mapped reads were aligned to the human reference genome hg19/GRCh37 with the Burrows-Wheeler Aligner to generate BAMs (46). BAMs were sorted and indexed with SAMtools (47). PCR duplicates were removed using PicardTools (http://broadinstitute.github.io/picard/). Peaks were called using the MACS2 (54) callpeak function on the ChIP-seq BAM file using the following parameters: BAMPE for paired-end reads, scaling to the larger dataset, the default log₂ fold change enrichment of 2 versus input and a P value cutoff of 1e-5. Consensus peaksets and normalized H3K27ac densities were generating using the R/Bioconductor package DiffBind (https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf) using the following parameters: score = score = DBA_SCORE_TMM_MINUS_FULL, bUseSummarizeOverlaps = T. Peaks were required to be present in at least 2 of 33 tumor samples or 1 of 3 normal samples. Overlapping peaks were merged. Peaks intersecting ENCODE blacklisted regions v1 (55) on haplotype chromsomes or on chromosomes X or Y were excluded from analysis. Peaks within 2.5 kb of a known transcription start site were also excluded. Bigwig tracks were generated using the DeepTools (v3.1.2; ref. 56) bamCoverage command with RPKM normalization. Genomic coverage heat maps were generated using the DeepTools (56) computeMatrix followed by plotHeatmap functions. Peaks were visualized using Integrative Genomics Viewer software (57). SEs were called using ROSE (24) and default parameters. Gene set enrichment plots for tumor SEs were generated using ClueGO (58) for GO BP, KEGG, and Reactome gene sets and visualized in Cytoscape (45). Tumor-specific and normal-specific SEs were called by intersecting consensus tumor SEs with consensus normal SEs. Waterfall plots of log₂ fold change were derived using the normalized values from DiffBind. Tumor versus normal SE enrichments were generated by inputting the genes associated with SEs with a log₂ fold change > 0.5 into ClueGO (58).

Gene–enhancer pairs were called using the R/Bioconductor (http://www.R-project.org/) package InTad (29). Genes with a significant (P < 0.1), positively correlated enhancer were assigned to the enhancer within the transactivation domain that had the most significant correlation. Enhancers with no significantly correlated genes were assigned to the nearest expressed gene.

BAM files for tumor or normal H3k27ac ChIP-seq or the corresponding input samples were merged using sambamba (59) and centered around nucleosome-free regions using the Homer (60) findpeaks function with the “histone” and “nfr” flags to center peaks around the nucleosome-free regions. Enriched motifs in tumor peaks were identified using findMotifsGenome with the normal peaks as background and vice versa to identify normal-enriched motifs.

Nonnegative matrix factorization was performed using the R (http://www.R-project.org/) package NMF (61) on normalized enhancer and SE peak densities. Ranks from 2 to 10 were analyzed with “nrun = 100.” For enhancers, one sample was excluded from the final clustering due to a silhouette width < 0.5. For NMF clustering of methylation data, the top 10% most variable probes were used.

The PR network was derived by identifying genes highly correlated with PR (r > 0.4). A genome-wide list of peaks containing a PR motif was derived using HOMER (60) scanMotifGenomeWide function with default settings. This list was intersected with the enhancer peakset and genes that were correlated with PR and had a PR motif in their assigned enhancer. This final gene list was used to infer the PR-regulated gene network. This gene list was intersected with the identified SE list to generate the overrepresentation analysis for PR-regulated genes and SE-associated genes. The gene set enrichment plot for PR-regulated genes was generated using ClueGO for GO BP, KEGG, and Reactome gene sets and visualized in Cytoscape.

Drug predictions were performed using the genes associated with the top 100 SEs with the greatest mean H3K27ac signal density for each subgroup. This list was input to the LINCS Clue database (62) using the “Query” tool. The top 10 compound perturbations that were inversely correlated with the gene signature were selected for each subgroup. The most negatively correlated drug representing that category was selected to represent the compound in the drug panel.

TCGA pan-cancer data was downloaded from freeze 1.3 of the TCGA PanCan Atlas at the Synapse website (https://www.synapse.org/#!Synapse:syn4557014). Data were logged and plotted using the R (http://www.R-project.org/) package uwot (https://github.com/jlmelville/uwot) with the following paramters: n_neighbors = 100, learning_rate = 0.5, init = “random,” pca = 50, min_dist = 0.4, ret_model = T. ComBat batch-corrected, log₂-transformed meningioma RNA-seq was projected onto this map.

Roadmap H3K27ac ChIP-seq bigwig files were downloaded from https://www.encodeproject.org/matrix/?type=Experiment&award.project=Roadmap&searchTerm=H3K27ac&assembly=hg19. FASTQ files from this study were trimmed to 36 bp using BBDuk from BBMap (https://sourceforge.net/projects/bbmap/) to match the Encode Roadmap data and were processed as described in the ENCODE workflow(https://github.com/ENCODE-DCC/chip-seq-pipeline2). ENCODE blacklisted regions (55) were subtracted and bigwigs were merged across all 213 samples using DeepTools (56) multiBigwigSummary with a bin size of 10,000. The top 10% most variable enhancers were input into consensusClusterPlus (ref. 63; maxK = 25, reps = 50). Optimal cluster size was determined by the greatest change in area under the cumulative density function curve at K > 3.

Enhancer RNA (eRNA) data processed from TCGA pan-cancer RNA-seq data (28) were kindly provided by Chen and colleagues (28). Single-end FASTQ files were trimmed with Trim Galore (39) and aligned using the workflow provided by NCI Genomic DNA Commons (https://docs.gdc.cancer.gov/Data/Bioinformatics_Pipelines/Expression_mRNA_Pipeline/). The STAR (64) 2-pass tool was used to generate BAM files. Gene expression values were derived using HT-Seq (65) count and summarized to gene level according to Gencode v22 annotation. The data were imported into R using DiffBind (https://bioconductor.org/packages/release/bioc/vignettes/DiffBind/inst/doc/DiffBind.pdf) and signal was called restricted to a bed file of the ∼10,000 enhancers used in the original report (28). Data were normalized using log₂ RPKM to match the original report. The datasets were merged, quantile normalized across sample, and then mean-centered across each eRNA. The optimal cluster size was determined using the R/Bioconductor (http://www.R-project.org/) package ConsensusClusterPlus (ref. 63; maxK = 25, reps = 50). Optimal cluster size was determined by the greatest change in area under the cumulative density function curve at K > 3.

Grade of membership models were derived using the R/Bioconductor (http://www.R-project.org/) package CountClust (66) at the default tolerance of 0.1. Cluster membership for each sample was averaged across tissue type, and the resulting tissue cluster scores were mean-centered for each cluster across tissue type. Z-scores were plotted in a heat map.

PR signature scores were derived using ssGSEA from GenePattern (67) using the derived PR signature. AR signature scores were derived using ssGSEA using the PID_AR_PATHWAY from C2-curated gene sets from the Molecular Signatures Database (MsigDB; ref. 43).

Enhancer networks were derived using the normalized values for the tumor consensus enhancer peakset and the R (http://www.R-project.org/) package WGCNA (36). The pickSoftThreshold function was used to select the lowest soft thresholding power β that demonstrated a scale-free topology model fitted with an R² > 0.9. An adjacency matrix accounting for only positive correlations was generated with β = 8. The dynamicTreeCut method was used with a minimum cluster size of 40 and height of 0.998 to create a dendrogram and modules from the dissimilarity matrix, generating 53 modules. The moduleEigengenes function was used to calculate the eigen-enhancer for each module. The eigen-enhancers were used in clinical correlation analyses with relevant statistics described in figure legends.

Ternary plots for group- or grade-enriched SEs were derived by summing the log₂ fold change for each SE to 1 and creating a ratio by dividing the log₂ fold change for each subgroup or grade over the summed log₂ fold change. These ratios were then squared and plotted with each subgroup ratio indicating one dimension along the ternary axes. Gene set enrichment plots for tumor SEs were generated using ClueGO²⁷ for GO BP, KEGG, and Reactome gene sets.

TF networks were derived using the Regulatory Genomics Toolbox (RGT; ref. 68). All tumor BAMs were merged using sambamba (59) merge and the RGT HINT tool was used to identify TF footprints in the merge tumor BAM across the tumor enhancer peakset. Motif enrichment was calculated for all tumor enhancers and for the enhancers of each individual module using the motifanalysis matching tool. Motifs from HOCOMOCO (69). Homer (60) and JASPAR (70) were used. Motif enrichment for each module versus the background of all tumor enhancers was then calculated using Fisher exact test. To account for multiple testing across many enhancers, P values were adjusted using the Benjamini–Hochberg method. If multiple motifs for the same TF were identified, the lowest motif enrichment FDR was used.

For each module, enriched motifs were filtered using RNA-seq data to include expressed TFs that were correlated with the module eigen-enhancer r > 0.1. To determine TF networks for a subgroup, TFs filtered by the above criteria were ranked by motif enrichment within each subgroup-enriched module. TF ranks were then summed across subgroup-enriched modules and the top 10 TFs with the lowest rank value (i.e., most enriched) that were present in at least 50% of the subgroup-enriched modules were selected for visualization. TF networks were visualized using Cytoscape (45). GSEA for TFs in a module or subgroup was performed using Enrichr (71).

All data has been deposited and can be downloaded using the SRA accession PRJNA579990 (WES) or GEO accession GSE139652 (ChIP-seq and RNA-seq).

A.R. Morton reports grants from NIH (CA236313) during the conduct of the study. M.A. Vogelbaum reports personal fees from Tocagen, other from Infuseon Therapeutics (indirect equity, royalty rights), and other from Celgene (funding to institution for a clinical trial) outside the submitted work. J.N. Rich reports grants and personal fees from Synchronicity Pharma outside the submitted work. No potential conflicts of interest were disclosed by the other authors.

B.C. Prager: Conceptualization, data curation, validation, investigation, methodology, writing-original draft, writing-review and editing. H.N. Vasudevan: Resources, data curation, investigation. D. Dixit: Investigation. J.A. Bernatchez: Investigation. Q. Wu: Investigation. L.C. Wallace: Investigation. S. Bhargava: Investigation, visualization. D. Lee: Investigation. B.H. King: Investigation. A.R. Morton: Formal analysis. R.C. Gimple: Investigation. M. Pekmezci: Investigation. Z. Zhu: Investigation. J.L. Siqueira-Neto: Supervision. X. Wang: Investigation. Q. Xie: Investigation. C. Chen: Supervision, writing-review and editing. G. Barnett: Resources.M.A. Vogelbaum: Resources. S.C. Mack: Conceptualization, supervision. L. Chavez: S, formal analysis, supervision, writing-original draft, writing-review and editing. A. Perry: Supervision, writing-originaldraft, writing-review and editing. D.R. Raleigh: Resources, supervision, writing-original draft, writing-review and editing. J.N. Rich: Conceptualization, resources, supervision, funding acquisition, writing-original draft, project administration, writing-review and editing.

This work was supported by grants provided by NIH (CA217066, to B.C. Prager; CA217065, to R.C. Gimple; CA197718, CA238662, NS103434, to J.N. Rich), CPRIT award, ALSF Young Investigator Award, Rally Research Grant, BEAR Necessities Pediatric Cancer Foundation Grant, Children's Cancer Research Fund award, and Baylor College of Medicine Junior Faculty Award (to S.C. Mack). We appreciate critical input and feedback from the members of the Rich laboratory. We also thank the Gillespie lab at the University of Alabama-Birmingham and the Jensen lab at the University of Utah for providing CH157-MN and IOMM-Lee cells, respectively. We appreciate the assistance of Mary McGraw and Dr. Christopher Hubert in obtaining samples. We also thank the patients and donors who provided the tissue used in this study and acknowledge Dr. Richard Drake and Dr. Jennifer McBride from the Cleveland Clinic for their assistance in deriving arachnoid granulation models.

The costs of publication of this article were defrayed in part by the payment of page charges. This article must therefore be hereby marked advertisement in accordance with 18 U.S.C. Section 1734 solely to indicate this fact.