Selective gene dependencies in MYCN-amplified neuroblastoma include the core transcriptional regulatory circuitry.

Childhood high-risk neuroblastomas with MYCN gene amplification are difficult to treat effectively1. This has focused attention on tumor-specific gene dependencies that underlie tumorigenesis and thus provide valuable targets for the development of novel therapeutics. Using unbiased genome-scale CRISPR-Cas9 approaches to detect genes involved in tumor cell growth and survival2-6, we identified 147 candidate gene dependencies selective for MYCN-amplified neuroblastoma cell lines, compared to over 300 other human cancer cell lines. We then used genome-wide chromatin-immunoprecipitation coupled to high-throughput sequencing analysis to demonstrate that a small number of essential transcription factors-MYCN, HAND2, ISL1, PHOX2B, GATA3, and TBX2-are members of the transcriptional core regulatory circuitry (CRC) that maintains cell state in MYCN-amplified neuroblastoma. To disable the CRC, we tested a combination of BRD4 and CDK7 inhibitors, which act synergistically, in vitro and in vivo, with rapid downregulation of CRC transcription factor gene expression. This study defines a set of critical dependency genes in MYCN-amplified neuroblastoma that are essential for cell state and survival in this tumor.

We conducted systematic genome-scale CRISPR-Cas9 screens in 341 human cancer cell lines, including nine neuroblastoma lines with MYCN-amplification[7]. Our cell line panel also included an additional two neuroblastoma cell lines, SKNFI and SKNAS, which do not harbor MYCN-amplification but do express high levels of MYCN or c-MYC, respectively[8]. We identified 147 candidate genes that were selectively essential to the growth and survival of MYCN-amplified neuroblastoma cell lines (Supplementary Table 1). The guide RNAs targeting these genes showed greater depletion in the nine MYCN-amplified neuroblastoma cell lines than in the other human cancer cell lines representing 26 tumor types. Classification of the 147 dependencies by gene ontology molecular function revealed a diverse set of genes encoding proteins in numerous functional categories (Fig. 1a)[9,10], with overrepresentation of ontologies for nucleic acid binding (12%, p=8.48×10−3) and transcription factor activity or binding (20%, p=6.78×10−3). Many dependency genes have not been previously associated with the sympathoadrenal cell lineage. Others have been reported to be involved in neural crest cell specification and neuroblastoma pathogenesis[1,11-20]. This set of transcriptional regulatory genes is highly connected by putative protein-protein interactions in the STRING database[21] (Fig. 1b). Representative dependency scores in MYCN-amplified neuroblastoma compared to other tumor cells are shown in Fig. 1c, d for the genes HAND2 and ISL1, and for a subset of other genes with nucleic acid binding and/or transcription factor activities in Supplementary Fig. 1. Our results show that many transcription factor genes with selective dependency in MYCN-amplified neuroblastoma cells also displayed dependency in the two MYCN non-amplified cell lines studied so far, although additional MYCN non-amplified lines will need to be studied before firm conclusions can be drawn. We next performed transient small-interfering RNA (siRNA) knockdown targeting mRNAs that encode five transcription factors dependencies (HAND2, ISL1, PHOX2B, GATA3 and TBX2) to independently verify that these genes represent dependencies in MYCN-amplified neuroblastoma cell lines. Depletion of expression of these genes resulted in suppression of colony formation and induction of apoptosis (Fig. 1e, f and Supplementary Figs. 2, 3). Thus, MYCN-amplified neuroblastoma cells rely on a discrete set of genes that are either unique to, or only minimally shared with, other tumor types.

Figure 1

Genome-scale CRISPR-Cas9 screening identifies selective neuroblastoma gene dependencies enriched for transcription factors

a. Gene ontology classification of terms associated with selective dependencies in neuroblastoma reveals enrichment for transcription factor activity or binding and nucleic acid binding proteins. n=147 genes, 252 GO-slim molecular class assignments. p=8.48×10−3 nucleic acid binding, p=6.78×10−3 transcription factor activity or binding, by 1-sided Fisher exact test (adjusted using Benjamini-Hochberg correction) comparing the ontologies of dependencies to all assayed genes in the genome.

b. STRING database analysis demonstrates 25 of 30 transcription factor dependencies have putative protein-protein interactions. Edge widths correspond to the level of confidence in interactions (medium confidence STRING score 0.4; high confidence STRING score 0.7; highest confidence STRING score 0.9). Red nodes indicate transcription factors previously co-reported with neuroblastoma in a literature search.

c and d. Representative scatter plots showing neuroblastoma relative dependency on HAND2 (c) and ISL1 (d) with p<1×10−9 by permutation testing with n=9 MYCN-amplified neuroblastoma cell lines compared to n=330 non-neuroblastoma cancer cell lines. Y-axis shows the gene’s dependency rank in an individual cell line. X-axis shows the gene’s dependency score in each cell line.

e and f. Colony formation assay at 10 days in BE2C (e) and Kelly (f) cells treated with two independent small-interfering RNAs targeting ISL1, PHOX2B, HAND2, GATA3 and TBX2. BE2C/Kelly siRNA#1 n=7 (control), 3 (HAND2, ISL1, PHOX2B), 4 (GATA3), 6 (TBX2); BE2C/Kelly siRNA#2 n=3 (control, HAND2, ISL1, PHOX2B), 6 (GATA3, TBX2) biologically independent samples; BE2C p<0.05 for all transcription factors relative to control siRNAs by 2-sided t-test. Center values represent mean and error bars represent SD.

The gene dependencies we identified suggest that diverse sets of complex pathways are activated to support the viability of MYCN-amplified neuroblastoma. Given the overrepresentation of transcription factor and nucleic acid binding activities among the 147 gene dependencies (Fig. 1a), we hypothesized that these genes were specifically required for transcriptional regulation in neuroblastoma cells. Previous studies have demonstrated that lineage-specific differences in cell state during development are determined by tissue-specific core regulatory circuits (CRCs) that may be subverted to drive malignant transformation[22,23]. Core regulatory circuitry is formed by transcription factors assembled into positive, feed-forward autoregulatory loops that establish and maintain cell lineage and identity through their extended regulatory networks[24,25]. CRC transcription factor genes are typically marked by extensive stretches of acetylation at histone H3-lysine 27 (H3K27ac), a genomic feature characteristic of high levels of gene expression, termed super-enhancers (SEs)[26,27] (Fig. 2a). To investigate the relevance of a CRC framework to the dependency genes identified earlier, we performed chromatin-immunoprecipitation coupled to high-throughput sequencing (ChIP-seq) for H3K27ac in five MYCN-amplified neuroblastoma cell lines, revealing a shared set of 77 genes associated with SEs across all five cell lines (Fig. 2b and Supplementary Table 2). As observed in the selective dependencies, the SE-associated genes in MYCN-amplified neuroblastoma cells were overrepresented among ontologies[9,10] including nucleic acid binding (21%, p=3.93×10−3) and transcription factor activity or binding (22%, p=4.29×10−8). Examples of the H3K27ac signals of highly ranked transcription factor genes are shown in Fig. 2c and Supplementary Fig. 4.

Figure 2

Several transcription factor dependency genes are marked by extensive H3K27ac and are enriched for dependency in human MYCN-amplified neuroblastoma

a. Normalized ChIP-seq tracks for H3K27ac demonstrate a SE at the GATA3 locus compared to a typical enhancer at the RAD23B locus in BE2C; tracks represent a combination of 2 independent experiments. ChIP-seq read densities (y-axis) were normalized to reads per million reads sequenced in each sample.

b. Gene ontology classification of terms associated with 77 shared SE-associated genes across five MYCN-amplified neuroblastoma cell lines reveals enrichment for transcription factor activity or binding (p=1.48×10−6), and nucleic acid binding proteins (p=1.40×10−2); n=77 genes, 130 GO-slim molecular class assignments, using 1-sided Fisher exact test (adjusted using Benjamini-Hochberg correction) comparing ontologies in SE-associated genes to all genes in the genome.

c. Identification of enhancers by ranked H3K27ac signal across all genes in the MYCN-amplified cell lines Kelly and BE2C. Highlighted are transcription factor dependencies with shared SEs across all five MYCN-amplified cell lines.

d. The rank of CRISPR-Cas9 dependencies genome-wide demonstrates selective dependency of shared SE-associated genes. Eleven of 69 SE-associated genes were enriched for selective dependency in neuroblastoma (p=1.11×10−10 by 2-sided Fisher exact test compared to all genes assayed). Highlighted are 10 of 11 SE-associated dependency genes annotated with transcription factor activity/binding or nucleic acid binding ontologies. Closed circles reflect transcription factors evaluable by ChIP-seq; open circles represent those not evaluated. X-axis shows gene rank in enrichment analysis (modified Kolmogorov-Smirnov test with permutation testing) for MYCN-amplified neuroblastoma versus other cancer cell lines. Y-axis shows the -log10(p-value) from enrichment analysis. Enrichment p-value for MYCN, HAND2, ISL1 and LDB1 was <1×10−9.

When these 77 SE-associated genes were compared with the 147 selective dependencies in MYCN-amplified neuroblastoma identified by the CRISPR-Cas9 screen, 11 of the SE-associated genes were selective dependencies and 10 of 11 (91%) of these genes encoded transcription factors (Fig. 2d). This represents 11 of 69 SE-associated genes targeted by specific guide RNAs (15.9%, p=1.11×10−10 for enrichment), since eight of the SE-associated genes were non-coding and were not targeted in the screen. To verify that these genes are expressed in primary human neuroblastomas, we examined the expression level of each SE-associated, CRISPR-Cas9 dependency nucleic acid binding or transcription factor gene in four primary neuroblastoma tumor expression datasets. The majority of the SE-associated CRISPR-Cas9 dependency genes exhibited high levels of expression in primary neuroblastomas as compared to multiple other types of primary tumors (Supplementary Fig. 5) with LDB1 as an exception.

Recent evidence has suggested the presence of two super-enhancer-associated transcriptional networks controlling neuroblastoma cell state, one representing a multipotent neural crest cell-like or mesenchymal cell state (NCC/MES) and a second representing sympathoadrenal cells, referred to as a noradrenergic (NOR/ADRN) cell state[28,29]. The majority of cell lines and low passage primary tumors, regardless of MYCN-amplification, displayed SEs associated with this core group of transcription factors (Supplementary Fig. 6a, b). Furthermore, analysis of neuroblastoma cell lines subdivided by transcriptional network subtype demonstrated that these CRC transcription factors were also marked by SEs, although to a lesser degree in the NCC/MES cell state. Analysis of all neuroblastoma cell lines in our genome-scale CRISPR-Cas9 screen demonstrated gene expression signatures consistent with either NOR/ADRN and NCC/MES cell states (Supplementary Fig. 6c). Expression of CRC transcription factor genes was elevated in both neuroblastoma subtypes compared to cell lines derived from other tumor types (Supplementary Fig. 6d). Further, these five CRC transcription factors represented gene dependencies in nearly all neuroblastoma cell lines, regardless of MYCN-amplification status or NOR/ADRN and NCC/MES cell state, indicating a general requirement for these genes in regulating cell growth and survival in neuroblastoma (Supplementary Fig. 6e).

Given these findings, we hypothesized that the CRC responsible for specifying neuroblastoma cell fate is formed by these SE-marked transcription factor dependencies. To test this hypothesis, we performed ChIP-seq in BE2C and Kelly cells with antibodies recognizing MYCN and five SE-regulated transcription factor dependencies. These data demonstrated that HAND2, ISL1, PHOX2B, GATA3, TBX2, and MYCN bind consensus DNA sequences adjacent to each other to form dense clusters of transcription factors comprising distinct “epicentres,”[30,31] within their own enhancer regions (Fig. 3a and Supplementary Fig. 7). These epicentres occupied discrete regions of open chromatin delineated by ATAC-seq. Binding of each transcription factor was preferentially found in gene enhancers containing enriched consensus binding motifs (Supplementary Table 3, Supplementary Fig. 8) suggesting that these transcription factors collectively regulate the expression of their own genes, and each of the other five transcription factors evaluable by ChIP, and form a feed-forward autoregulatory CRC. When overexpressed, MYCN is thought to bind active promoter and enhancer elements throughout the genome, similar to the activities of MYC[32,33], resulting in amplified transcription that overcomes rate-limiting constraints on cell proliferation and survival. Genome-wide ChIP-seq analysis demonstrated that MYCN binds coordinately with other members of the CRC to regions of open chromatin marked by enhancers across the genome (Fig. 3b), indicating that it may function as a CRC member or, alternatively, as a general transcriptional amplifier as has recently been reported[34]. Amplification of the MYCN gene raised the background signal of input DNA across the amplified region, precluding formal evaluation of CRC transcription factor binding to any potential MYCN enhancer (Supplementary Fig. 9). Analysis demonstrated a nearly two-fold statistical enrichment of binding of all 5 CRC transcription factors, with or without MYCN, to dependency genes (p<0.01 for binding to dependency genes versus remainder of expressed genes in the CRISPR-Cas9 screen by two-sided Fisher’s exact test). However, MYCN and these transcription factors are also present at genes that are not selective neuroblastoma dependencies.

Figure 3

Dependency transcription factors form the core regulatory circuitry in MYCN-amplified neuroblastoma

a. Normalized ChIP-seq tracks for H3K27ac, MYCN, and five transcription factors at the HAND2, ISL1 and PHOX2B loci in BE2C cells; H3K27ac track represents a combination of 2 independent experiments in BE2C and other tracks are representative of an independent experiment performed in BE2C and Kelly cells. ChIP-seq read densities (y-axis) were normalized to reads per million reads sequenced in each sample. SEs are noted as red bars and arrows indicate epicentres.

b. Genome-wide co-occupancy for CRC transcription factors as determined by ChIP-seq. Regions (rows) were defined as those enriched in ChIP-seq reads for at least one transcription factor and are ranked by the MYCN signal therein. Color keys for reads-per-million-normalized signal are displayed below each heatmap.

c. Quantitative RT-PCR of BE2C cells treated with transient siRNA to each CRC gene - HAND2, ISL1, PHOX2B, GATA3 and TBX2 - results in decreased expression of mRNA transcripts for all of the CRC members, which was not observed for non-targeting control siRNA transfection (n=3 independent biological experiments, all siRNA-treated transcription factor gene expression are significantly different from both control siRNAs with p<0.05 by 2-sided t-test. Horizontal lines demonstrate the median with upper and lower box boundaries demonstrating the 25-75th centiles. Upper and lower bounds represent the 10th-90th centiles.

d. HAND2, ISL1, PHOX2B, GATA3, and TBX2 form a positive feedback, interconnected co-regulatory loop. MYCN regulates each of these genes as a part of the CRC. SEs and gene loci are represented by rectangles, and proteins are represented by oval symbols.

To establish co-regulation of CRC component genes, we determined the gene expression levels by qRT-PCR after transient siRNA-mediated knockdown of each gene (Fig. 3c and Supplementary Fig. 10). We observed decreased expression levels of all CRC transcription factors after knocking down each individual CRC gene. The impact of CRC transcription factor depletion on the super-enhancer landscape of the cell remains to be studied; however, these data implicate these genes and their encoded transcription factors in a feed-forward autoregulatory mechanism that drives MYCN-amplified neuroblastoma (Fig. 3d).

Efforts to target aberrant transcriptional activity in tumors such as MYCN-amplified neuroblastoma have centered on agents that inhibit cyclin-dependent kinases (CDKs) and BET-bromodomain-dependent pathways that act as mediators of MYCN-driven initiation and elongation of transcription[35,36]. Based on the feed-forward, autoregulatory CRC we have identified, we postulated that MYCN-amplified neuroblastoma cells could be targeted synergistically by combining JQ1 (BRD4 inhibitor) and THZ1 (CDK7 inhibitor). While treatment with JQ1 or THZ1 alone reduced the growth of Kelly and BE2C cells (Fig. 4a), combination treatment dramatically reduced cell numbers and induced apoptosis (Fig. 4a, Supplementary Fig. 11). Comparison of results after treatment with multiple concentrations of JQ1 and THZ1 alone and together demonstrated synergistic inhibition of viability of a panel of both MYCN-amplified and non-amplified neuroblastoma cell lines by the Chou-Talalay combination index (Fig. 4b and Supplementary Fig. 12)[37-39]. Next, we subcutaneously xenografted BE2C cells in nude mice and treated tumor-bearing mice with either vehicle control, JQ1 and/or THZ1. In contrast to control and single agent-treated mice, combination inhibition with JQ1 and THZ1 led to reduced tumor progression and increased survival compared to either single agent alone without a substantial effect on animal weight (Fig. 4c, Supplementary Fig. 13). These data indicate that compounds disrupting transcriptional initiation and elongation can serve as a synergistically potent and effective method of tumor inhibition in MYCN-amplified neuroblastoma in vitro and in vivo.

Figure 4

Pharmacological disruption of SE-mediated transcription selectively suppresses the core regulatory circuitry and neuroblastoma cell growth

a. BE2C and Kelly cells demonstrate decreased growth with JQ1 and/or THZ1 treatment (JQ1 3 μM (all); THZ1 125 nM (BE2C) and 78 nM (Kelly); alone or in combination or DMSO control). n=4 biologically independent experiments. p<0.001 for combination treatment at all timepoints relative to control, JQ1 or THZ1 alone; p<0.001 for JQ1 or THZ1 alone relative to DMSO at 72h by 2-way ANOVA with Tukey post hoc correction. Mean relative growth is plotted, error bars represent SD.

b. Chou-Talalay normalized isobolograms depicting combination index (CI) scores over a range of THZ1 and JQ1 concentrations in BE2C and Kelly (CI scores <1=synergy, >1=antagonism; r ed line represents additivity, CI=1).

c. BE2C xenografts demonstrate reduced growth with combinations of JQ1 and THZ1 compared to JQ1, THZ1 or DMSO vehicle controls. n=8 mice per treatment group. Mean tumor volume is plotted with error bars representing SEM. N.S.=not significant, *** vehicle vs combination p=0.0001, JQ1 vs combination p=0.021, THZ1 vs combination p=0.0002 by 2-way ANOVA with Tukey post hoc correction.

d. Gene expression analysis in BE2C cells treated with JQ1 and THZ1 for 1, 4 and 12h. Data is normalized to untreated cells at 0h and levels of ERCC spike-in RNAs.

e. Log2-fold change of CRC gene transcripts normalized to the top 1% highest expressed genes at 0h in BE2C cells; n=3 independent experiments. CRC genes display reductions in expression relative to the top 1% highest expressed genes in the genome at all timepoints (1h p=0.0059, 4h p=0.0021, 12h p=0.0032 by 2-sided t-test).

To evaluate the effects of combination treatment on both CRC and global gene expression, we treated BE2C cells with a synergistic combination of JQ1 for 0, 1, 4 and 12h, and performed gene expression analysis with spike-in controls. In comparison to the top 1% of highest expressed transcripts, CRC mRNA levels rapidly decreased in expression by 1h of treatment (Fig 4d, e). We observed that the majority of downregulated genes at 4h of combined JQ1 and THZ1 treatment were bound by all 6 CRC members, suggesting they represent a network of genes directly regulated by the CRC (Supplementary Fig. 14). By 12h of treatment, the majority of all mRNA transcripts were downregulated, indicative of general transcriptional collapse. These results were confirmed by qRT-PCR analysis of expression levels of CRC and non-CRC member genes (ACTB, HPRT) in BE2C and Kelly cells (Supplementary Fig. 15). The expression level of each of the six transcription factor genes was dramatically reduced by the combination of JQ1 and THZ1, with more limited effects noted with either drug alone. These results underscore the effects of JQ1 and THZ1 combination treatment in MYCN-amplified neuroblastoma, though the complete effects of combination transcriptional disruption remain to be elucidated.

Our findings identify a collection of genes representing diverse vulnerabilities in MYCN-amplified neuroblastoma, including the transcription factors that form an autoregulatory CRC specific for these cells. At least six of the transcription factors encoded by members of the dependency group belong to this CRC, each of which directly regulates the expression of its own gene as well as those encoding the other CRC transcription factors. This CRC also controls an extended regulatory network of genes, which we suggest contributes to initiation and maintenance of the transformed phenotype in MYCN-amplified neuroblastoma. We hypothesize that a similar set of transcription factors form a CRC in MYCN non-amplified neuroblastoma, due to reliance on high level expression of MYC or MYCN in these tumors. Our results underscore the interconnected autoregulatory mechanisms underlying the expression and activity of CRC transcription factors. This dependence on transcriptional activity leads to the enhanced vulnerability of these cells to combinatorial targeting of key nuclear proteins that mediate transcriptional initiation and elongation. A subset of the selective dependencies that we have identified in this tumor likely reflect the consequences of transformation by the CRC, and it will be important to investigate these genes and the pathways they control as a source of additional “druggable” targets in this tumor.

Online Methods

Cell Lines

Kelly and BE2C neuroblastoma cell lines were obtained from ATCC (BE2C) and DSMZ (Kelly). All cell lines used for the genome-scale CRISPR-Cas9 screen are detailed in Meyers et al[7]. All cell lines were STR tested for identity. Neuroblastoma cell lines were cultured in RPMI media containing 10% FBS and 1% penicillin-streptomycin and validated to be free of Mycoplasma species.

In vivo tumor models

We adhered to animal protocols approved by the Dana-Farber Cancer Institute Animal Care and Use Committee. Animals were maintained according to institutional guidelines. Eight-week old female nude mice were subcutaneously implanted with 1×106 BE2C cells in 30% matrigel/PBS. Mice were assigned to four groups: vehicle, JQ1, THZ1, or the combination of JQ1 and THZ1. Treatment with small molecule inhibitors was initiated once tumors engrafted and reached 100-150 mm3. Compounds were solubilized in 10% DMSO and 90% 5% dextrose in water. JQ1 was delivered at 50 mg/kg daily by intraperitoneal injection, and THZ1 was delivered at 10 mg/kg twice daily by intraperitoneal injection. Mice were treated for 14 days and followed for survival. Tumors were measured by calipers twice weekly. Mice were weighed twice a week. Animals were sacrificed according to institutional guidelines when tumors reached 2000 mm in length or width or if animals became moribund.

CRISPR-Cas9 Screening and Analysis

The CRISPR-Cas9 screen was performed with the Avana library containing approximately 70,000 guides and an average of 4 guides per gene[4]. Details of library construction, cell transduction and analysis are found in the Supplementary Note. This version of the Avana dataset contains 341 cell lines, including nine MYCN-amplified neuroblastoma lines (CHP212, IMR32, Kelly, KPNYN, MHHNB11, NB1, SKNBE2, SKNDZ, and SIMA). This dataset includes two MYCN non-amplified neuroblastoma cell lines, SKNAS and SKNFI, which were excluded from the enrichment analysis for a total of 339 cell lines analyzed. Other cell lines in the screen represent 26 different types of tumors (breast 26, chondrosarcoma 1, colorectal 25, endometrial 12, esophageal 10, Ewing 1, gastric 10, glioma 31, head and neck 8, leukemia 9, liver 13, lung non-small cell 45, lung small cell 2, lymphoma 4, medulloblastoma 4, melanoma 29, osteosarcoma 4, ovarian 31, pancreatic 22, renal 16, rhabdoid 1, soft tissue sarcoma 5, thyroid 1, urinary tract 19) and were combined to create the comparison group, as detailed in Meyers et al[7].

Dependency scores were generated for each gene in each cell line using the CERES algorithm[7]. Briefly, CERES is a computational method to estimate gene dependency levels from CRISPR-Cas9 essentiality screens while accounting for the copy-number-specific effect, as well as variable sgRNA activity. The gene dependency scores generated by CERES are absolute dependency scores with each cell line scaled to have an average dependency score of 0 for negative control sgRNAs and an average dependency score of -1 for a set of positive control core essential genes.

To generate a list of selective dependencies in MYCN-amplified neuroblastoma, the CERES data were mean centered for each gene. The genes were then selected as dependencies if they had more than 2 MYCN-amplified neuroblastoma lines that scored below 4 standard deviations from the mean, 3 lines that scored below 3 standard deviations from the mean, or were enriched for dependency below a threshold of 2 standard deviations from the mean using the Kolmogorov–Smirnov test with a false discovery rate (FDR) less than 0.05 using multiple testing correction of the empiric p-values. Candidate dependencies were filtered to only include those that also had an absolute dependency effect prior to the mean centering described above. Genes that had overall negative scores across all cancer lines in the screen were considered to be likely pan-essential and were excluded. Candidate dependency genes were further filtered to include those with expression in MYCN-amplified neuroblastoma in the Cancer Cell Line Encyclopedia (CCLE)[40]. Thus, final scores reflect selective dependencies relative to findings with all other tumor types in the screen.

Gene Ontology Classification

Gene ontologies were assigned using the GO-slim molecular function ontologies. These 42 GO-slim molecular functions were combined to create broader classifications of genes such as transcription factor activity or binding (Supplementary Table 4). In each class of ontologies, a gene was counted only once for membership to any of the defining ontologies. A Fisher’s exact test was performed to look for enrichment in the set of 147 dependency genes or 77 super-enhancer associated genes versus the genome as a whole. P-values were adjusted using the Benjamini-Hochberg multiple testing correction.

Quantitative PCR

Total RNA was harvested using Trizol (Life Technologies) according to the manufacturer’s protocol. First-strand synthesis was performed with Superscript III (Invitrogen). Quantitative PCR analysis was conducted on the Vii7 system (Life Technologies) with SYBR Green PCR Master Mix (Roche) using validated primers specific to each target each gene. Primer sequences are displayed in Supplementary Table 5.

Microarray Analysis

BE2C cells grown in 6 well plates were treated with JQ1 (3 μM) and THZ1 (125 nM) and collected at 0, 1, 4, and 12 h with Trizol (Life Technologies). ERCC spike-in RNA was added directly to Trizol prior to RNA extraction as previously described[41]. Samples were prepared using the 3′ IVT assay (Affymetrix) and labeled by biotinylation. Following labeling, biotinylated cRNA was hybridized to the Primeview human 3′ GeneChips for 16 h at 45°C in an Affymetrix Hybridization Oven 645. Arrays were washed and stained in the Affymetrix Fluidics Station 450 using the GeneChip Hybridization, Wash, and Stain Kit. The arrays were scanned on GeneChip Scanner 3000 7G. For further details of data analysis, refer to the Supplementary Note.

Western Blotting

Whole cell lysates were prepared in Buffer C as previously described[23]. Equivalent amounts of protein were separated by SDS-PAGE and blocked by standard methodologies. Immunoblotting was carried out with the following antibodies: PHOX2B-B11X (Santa-Cruz Biotechnology), ISL1 clone 39.4D5 (Developmental Hybridoma Studies Bank), HAND2-A12X (Santa-Cruz Biotechnology), GATA3 (Pierce Technologies), TBX2 (Abcam Biotechnology), cleaved caspase-3 and PARP1 (Cell Signaling Technology). β-actin (Cell Signaling Technology) was used to confirm equal loading. Secondary horseradish peroxidase (HRP)-linked IgG antibodies to mouse or rabbit were from Cell Signaling Technology. Enhanced chemiluminescence agents were Pico and Femto ECL (Thermo-Fisher Scientific) and used according to the manufacturer’s instructions. All data displayed are representative of at least 3 independent protein lysates and blots.

siRNA Knockdown

Control and target siRNAs were purchased from Dharmacon (GE Biosciences). Cells were subjected to 100 nM concentrations of siRNA in OPTIMEM-I (Gibco) for 6h using lipofectamine 2000 (Invitrogen) as previously reported[42]. Knockdown was confirmed by western blotting at 24-48h post-transfection, and RNA was extracted at the time of maximal protein knockdown. Two independent single siRNA oligonucleotides were utilized per gene assayed, and all experiments were completed at least in triplicate.

Colony Formation Assay

Cells were transfected with siRNAs as described. After treatment, cells were replated at 500 cells/well in 6 well dishes, grown in regular growth media for 10 days prior to 100% methanol fixation, 0.05% crystal violet staining and subsequent quantitation. Experiments were completed in triplicate; data shown is the average of three independent experiments.

ATAC-sequencing

Assay for Transposase Accessible Chromatin using sequencing (ATAC-seq) was performed as previously described[43]. Briefly, Kelly and BE2C cells were cultured under standard conditions and 50,000 cells were trypsinized and washed twice with PBS. The Tn5 transposition reaction was performed with the Nextera kit (Illumina) according to the manufacturer’s protocol. Transposed DNA was then purified using the MinElute kit (QIAgen). Library amplification was performed using NEBNext High-Fidelity 2X PCR Master Mix (NEB), with SYBR Green I dye (Invitrogen) and Nextera primers Ad1 and Ad2.[x] at 25 μM. Amplified libraries were purified with the PCR purification kit (QIAgen), and quality control was performed using the 2200 Tapestation (Agilent). Sequencing was performed with Illumina NextSeq 500 paired-end 75bp reads.

Genome-Wide Occupancy Analysis

Chromatin immunoprecipitation (ChIP) coupled with high-throughput DNA sequencing (ChIP-seq) was performed as previously described[22,23]. The following antibodies were used for ChIP: GATA3 (Thermo Fisher MA1028), HAND2 (Santa Cruz sc9409X), ISL1 (Developmental Hybridoma Studies Bank 40.3A4), MYCN (Invitrogen Biotechnology, MA-1-16638), PHOX2B (Santa Cruz sc376997X), TBX2 (Abcam ab33298), and H3K27ac (Abcam ab4729). Details of ChIP protocols and analysis, as well as SE identification and assignment, are found in the Supplementary Note.

Synergy Studies and Pharmacologic Inhibitors

Synergy screening was performed with THZ1 synthesized by the Gray Laboratory (Dana-Farber Cancer Institute) and JQ1 synthesized by the Bradner Laboratory (Dana-Farber Cancer Institute). Experiments were validated with compounds from commercially available sources: JQ1 (Sigma-Aldrich) and THZ1 (EMD Millipore). For assessment of synergy, neuroblastoma cell lines were plated in 384-well tissue culture treated plates at a density of 25,000 cells/ml. Cells were then treated with single or combinations of compounds and subsequently analyzed for cell viability on days 0, 1 and 3 using the Cell-TiterGlo luminescent assay kit (Promega) per the manufacturer’s instructions. Luminescence was read on a Fluostar Omega Reader (BMG Labtech). For cell growth assays and gene expression studies, JQ1 was used at 3 μM and THZ1 at 78 nM (Kelly) or 125 nM (BE2C), with DMSO as a vehicle control. Synergy was assessed by the Chou-Talalay Combination index, isobologram and excess over Bliss methods[38].

Statistical Analysis

Data from the ChIP-seq and CRISPR-Cas9 screens were analyzed as described above. Animal experiments were analyzed by two-sided ANOVA for tumor volume and weight means, and by the log-rank test for survival. Other data were analyzed with one- or two-sided ANOVA with post hoc Tukey tests, two-sided t-tests, or one- or two-sided Fisher exact tests as appropriate for multiple or pairwise comparisons. Statistical significance was defined as a p-value <0.05. Data were analyzed with GraphPad Prism 7.01, and all error bars represent SD unless otherwise noted.

Data and Code Availability

Gene dependency scores from the CRISPR-Cas9 screen are publically available at https://doi.org/10.6084/m9.figshare.5520160.v1[7,44]. ChIP-seq (GSE94822), ATAC-seq (GSE94823) and spike-in microarray data (GSE108914) are available through GEO; for further details see Supplementary Table 6. Custom code written in R/python to perform analyses of CRISPR-Cas9 and ChIP-seq are available upon request.