Cross-Cohort Analysis Identifies a TEAD4–MYCN Positive Feedback Loop as the Core Regulatory Element of High-Risk Neuroblastoma

Neuroblastoma is a malignancy arising from the developing sympathetic nervous system that typically affects very young children. About half of the cases present with widespread metastatic disease and/or extraordinarily focal amplification of the MYCN gene, often with >100 copies. These high-risk patients have poor 5-year survival probability (<40%), despite intensive multimodal treatment regimens (1). Diagnostically, neuroblastomas show relatively few recurrent somatic point mutations (2–5). In contrast, it has been known for decades that high-risk neuroblastomas harbor complex and recurrent somatic structural copy-number alterations (CNA), affecting large chromosomal regions, as well as focal amplification of the MYCN oncogene (6). Identification of the causal driver genes associated within these CNAs remains a challenge. Despite significant understanding of the genetic landscape of high-risk neuroblastoma, all newly diagnosed patients are empirically treated with intensive cytotoxic chemoradiotherapy to achieve disease remission. Thus, further elucidation of the molecular mechanisms responsible for high-risk neuroblastoma pathogenicity is required for guiding novel, more effective and less toxic precision-oncology strategies.

Following on recent results from the assembly and interrogation of regulatory network models (interactomes) in human cancers, using systems biology approaches (7), we focused on the discovery of more universal tumor dependencies in high-risk neuroblastomas. This approach has been highly successful in elucidating small regulatory modules, dubbed tumor checkpoints, comprising master regulator (MR) proteins that mechanistically regulate the transcriptional state of the tumor. MR proteins have been shown to elicit either essentiality or synthetic lethality in several cancers (8–16). Critically, aberrant activity of tumor checkpoint MRs arises from one or more genetic alterations in their upstream pathways (14, 17), thus providing a rationale for the heterogeneous mutational landscape of tumors with highly similar transcriptional states (15). MR proteins have been shown to represent critical molecular tumor vulnerabilities that can be effectively targeted pharmacologically in several human malignancies (10, 11, 16, 18, 19), thus providing novel therapeutic opportunities. Here, we sought to interrogate regulatory networks of high-risk neuroblastomas to identify and validate MR proteins representing highly conserved, tumor-specific vulnerabilities for specific subtypes of this malignancy.

Our approach represents the integration of computational inference and experimental validation, both in vitro and in vivo, to characterize the tumor checkpoint MRs that mechanistically control the transcriptional state of distinct, high-risk neuroblastoma subtypes.

We first dissected the heterogeneity of high-risk neuroblastoma gene expression profiles (GEP) to classify them into molecularly consistent tumor subtypes. Transcriptionally distinct subtypes were identified using consensus clustering (20) on GEPs of 219 high-risk tumors assembled by the NCI TARGET (Therapeutically Applicable Research to Generate Effective Treatments) initiative, excluding good-prognosis samples (stage 1; refs. 2, 21). Optimal cluster selection identified three distinct high-risk subtypes (Supplementary Table S1), as determined by the relative area under the cumulative distribution function (Supplementary Fig. S1A–S1C). This result was then validated in an independent cohort comprising 97 high-risk samples assembled by the European Neuroblastoma Research Consortium (NRC; Supplementary Fig. S1D–S1F). Cross-cohort analysis of differentially expressed genes in each subtype, compared with the respective low-risk stage 1 samples, revealed a 1:1 correspondence between TARGET and NRC-derived high-risk subtypes (Fig. 1A and B).

Further clinical and biological characterization of these subtypes identified key distinguishing features (Supplementary Tables S1 and S2, respectively). A first cluster cosegregated with MYCN amplification status and chromosome band 1p36 deletions (78% and 83% of TARGET and NRC samples, respectively; Fig. 1C; Supplementary Fig. S1G and S1H). This predominantly MYCN^Amp (MYCNA) subtype showed a highly undifferentiated phenotype and high mitosis–karyorrhexis index (MK; Supplementary Fig. S1I and S1J, respectively), potentially resulting from downregulation of differentiation programs and upregulation of cell proliferation and cell growth pathways, respectively (Fig. 1E and F). Although the majority of MYCNA-subtype samples harbored MYCN amplifications resulting in MYCN overexpression, a few samples appeared to compensate for low MYCN mRNA levels by significant c-MYC (MYC) overexpression. Consistent with previous findings, MYCN/MYC deregulation in neuroblastoma may occur via transcriptional deregulation mechanisms, independent of their genomic amplification status (22, 23). All samples in this subtype showed high MYCN/MYC-specific signature activity (ref. 24; Supplementary Fig. S1M and S1N). Although samples in other high-risk subtypes presented broad 2p gains, the MYCNA subtype is the only one showing highly focal amplification of the MYCN locus (Fig. 1C), consistent with the fact that nonfocal 2p gain is not associated with MYCN overexpression (25, 26).

A second cluster was characterized by hemizygous deletions of 11q (90% and 89% of cases in TARGET and NRC, respectively) and was thus referred to as the 11qLOH subtype (Fig. 1C; Supplementary Fig. S1G and S1H). Overall, the clinical profile (Supplementary Fig. S1I and S1J) and pathway enrichment pattern were similar to the MYCNA subtype (Spearman correlation = 0.79; Fig. 1E), including enrichment of proliferation and cell-cycle categories and negative enrichment of differentiation pathways. Yet, this subtype exhibited inverse association with MYCN expression and activity (Supplementary Figs. S1M, S1N and S2I), as well as with activation of immune-related pathways (Fig. 1E).

In contrast, the third subtype did not appear to be strongly associated with specific genomic alterations (Fig. 1A and C; Supplementary Fig. S1G and S1H) and did not present with hyperproliferative program activation (Fig. 1F and G). Similar to the 11qLOH subtype, however, it displayed activation of immune-related pathways. Interestingly, this subtype presented a strong mesenchymal signature, highly similar to the one previously reported for high-grade glioma (Supplementary Fig. S2E–S2H; ref. 27) and will thus be referred to as the mesenchymal (MES) subtype.

To further disentangle tumor-specific signatures from those of tumor-infiltrating compartments in 11qLOH and MES subtypes, we used the ESTIMATE algorithm, which allows inferring both the stromal and immune fractions in each sample (ref. 28; Supplementary Fig. S2A–S2D). Integration of immune and stromal components by ESTIMATE suggested that the MES and a subset of the 11qLOH subtypes are characterized by lower purity, likely due to either a larger Schwannian stromal and/or immune reactive cellular infiltrate (Supplementary Fig. S2A–S2D and S2I). This is consistent with the previous reports of a higher fraction of tumor-infiltrating lymphocytes in MYCN^WT high-risk neuroblastomas (29). The stromal component enriched in these subtypes was confirmed by high fraction of stroma-rich ganglioneuroblastomas (Supplementary Fig. S1K; ref. 30). Conversely, MYCNA subtype samples showed minimal immune signature contamination, confirming known MYCN-related mechanisms of tumor immune suppression (31–33).

To assess whether the signatures emerging from cluster analysis represent cell-autonomous or microenvironment-related mechanisms, we tested whether they could be recapitulated in a panel of high-risk human neuroblastoma-derived cell lines. Interestingly, the MES signature was strongly conserved in cell lines (Supplementary Fig. S2F). Single sample gene set enrichment analysis (GSEA; ref. 34) identified SKNAS as the neuroblastoma cell line with strongest MES signature enrichment (NES = 4.54, P = 5.5E−6; Supplementary Fig. S2G and S2H). This confirms the tumor cell–autonomous mesenchymal nature of this subtype, because cell line cultures lack stromal or immune cell contamination.

In summary, MYCNA and 11qLOH subtypes displayed high activity of proliferative programs; the 11qLOH and MES subtypes showed high immune and stromal infiltration whereas the MES subtype showed cell-autonomous activation of mesenchymal programs (Supplementary Fig. S2I). Overall, all three subtypes were associated with poor survival in both TARGET and NRC datasets (Supplementary Fig. S1O and S1P), suggesting fundamentally different mechanisms leading to disease metastasis and ultimate therapy resistance.

Next, we inferred subtype-specific candidate MR proteins by independent analysis of TARGET and NRC cohort datasets. We first assembled TARGET- and NRC-specific interactomes from cohort-specific neuroblastoma GEPs, using ARACNe-AP (35), the latest version of ARACNe—an established tool for the reverse engineering of transcriptional targets of regulatory proteins (36). The TARGET and NRC interactomes comprised 205,271 and 359,846 transcriptional interactions, respectively (Supplementary Table S3), 81,035 of which were overlapping (P < 1E−16 by Fisher exact test, odds ratio = 65.02; Supplementary Fig. S3A).

Candidate MR proteins for each of the high-risk subtypes were then prioritized based on the enrichment of their transcriptional target genes in the subtype-specific signature, using the VIPER algorithm (19). Specifically, we used signatures representing the differential gene expression of each high-risk subtype compared with stage 1 samples (good prognosis), which showed significant overlap between the TARGET and NRC datasets (Supplementary Fig. S3B). Furthermore, despite completely independent MR analysis in the TARGET and NRC cohorts, top-ranking MR proteins for each molecular subtype were remarkably consistent. Indeed, overlap computed by one-sided Fisher exact test of the first 50 MRs of each subtype was highly significant (P_MES = 5.18E−11; P_11qLOH = 4.73E−35; P_MYCNA = 4.34E−33; Supplementary Fig. S3C and S3D). This confirms the robustness of the analysis, independent of cohort bias and composition. This was further confirmed by GSEA (ref. 37; |NES| > 4; P < 1E−4) for the MYCNA subtype (Fig. 2A). Given the exceedingly high congruence of the MRs identified from the two cohorts, we finalized candidate MR rank by Fisher integration of their cohort-specific P values (Fig. 1D; Supplementary Table S4).

We then proceeded to experimentally validate candidate MRs of the MYCNA subtype. We focused on this subtype for several reasons: (i) It is strongly associated with the most recurrent neuroblastoma-specific genetic alteration (38); (ii) it is associated with aggressive disease and poor prognosis (39); (iii) it shows the highest purity in terms of immune and stromal infiltration (Supplementary Fig. S2), thus resulting in high-quality GEPs for MR predictions; and (iv) it has the largest number of established cell lines (as matched by both MR and genetic analysis) for validation purposes (Supplementary Fig. S4). The 25 most significant candidate MRs of the MYCNA subtype were prioritized for experimental validation (Fig. 2B).

To identify optimal subtype representatives and appropriate negative controls for validation purposes, we analyzed a panel of 25 well-characterized neuroblastoma cell lines. Cell lines harboring MYCN amplifications as well as having high MYCN expression and signature activity (Supplementary Fig. S4A–S4D) were selected as candidate representatives of the MYCNA subtype and further refined based on the activity of the top 25 VIPER-inferred MRs selected for validation (Supplementary Fig. S4E; ref. 19). Most MYCNA cell lines showed high activity of these proteins, with BE2 identified as the model with the most statistically significant MR-activity match by enrichment analysis (P = 3.7E−24; Supplementary Fig. S4E). Negative controls were prioritized based on low MYCN/MYC expression and low VIPER-inferred activity of selected MRs, including MYCN.

To validate MYCNA-subtype MRs, we performed loss-of-function assays by silencing each candidate MR in both MYCNA and control cell lines. Multiple RNAi screenings were used to mitigate false discovery resulting from off-target effects and technology-specific biases. We tested a panel of MYCN-amplified and control cell lines for transduction efficiency, and only cells exhibiting high transduction efficiency were chosen for validation studies (Supplementary Fig. S4F).

Pooled, in vitro short hairpin (shRNA) screens were performed in BE2 and IMR5 cells (MYCNA subtype) and NLF and SKNAS cells (negative controls). NLF was included as an informative negative control because, despite its MYCN^Amp status, it showed both low MYCN expression (Z-score = −0.74) and MYCN activity (NES = −3.31; Supplementary Fig. S4D), as well as no enrichment of MYCNA MRs (Supplementary Fig. S4E).

MRs representing critical MYCNA subtype–specific dependencies were assessed by evaluating the depletion of shRNA hairpin representation at 28 days after infection. We then further confirmed the results from these assays in vivo, to ensure a more physiologic environment for MR validation, using BE2 and SK-N-AS xenografts, based on their known tumorigenicity in immune-deficient mice (Supplementary Experimental Procedures). Overall, genes were prioritized as bona fide MRs if (i) >2-fold decrease in hairpin representation was observed in MYCNA cells (P < 0.05, z score < −1.96) and (ii) the ratio between MYCNA and control group hairpin representation was >2-fold, in both the in vitro and the in vivo screens (Supplementary Table S5). Using this very conservative criterion, we identified MYCNA-specific MR dependencies, common to both in vitro and in vivo screens, including MYCN, TEAD4, HNRNPAB, HMGB2, PRDM8, E2F3, and ECSIT (Fig. 2C and D).

To further validate key MRs, we next performed lentivirus-mediated shRNA silencing of each candidate MR in an extended panel of cell lines, using the two shRNA hairpins with the highest silencing efficiency (Supplementary Fig. S5A; Supplementary Table S5), and assessed their cell viability (Supplementary Table S5; Supplementary Experimental Procedures). These shorter-term assays (at 72–96 hours after transduction) show that MYCNA cells were more sensitive to TEAD4, TAF1D, HNRNPAB, and ECSIT silencing compared with control cell lines (Fig. 2E; Supplementary Table S5). We did not detect MYCN as a significant hit in this screen at the selected time point, likely due to the very high MYCN copy number in MYCNA cells and technical difficulty to achieve sustained silencing. However, we confirmed that MYCN silencing in BE2 cells induced differentiation and viability reduction at 7 days after transduction (Supplementary Fig. S5B and S5C). We reasoned that shorter-time assays (<96 hours) are better suited for detecting MRs producing direct effect on proliferation, whereas longer-term assays (>96 hours) are optimally suited for elucidating multifunctional dependencies. Consistently, additional MYCNA-specific MR dependencies were detected in long-term (28-day) pooled shRNA screening (Fig. 2C and D).

To reduce the effect of potential off-target shRNA effects, we used an orthogonal silencing mechanism with ON-TARGETplus siRNA pools containing a mix of four siRNAs against each of the 25 MRs (Supplementary Experimental Procedures). By comparing the average relative cell viability upon siRNA-mediated MR silencing in MYCNA cell lines and control cell lines at 96 hours after transfection (Supplementary Fig. S5D; Supplementary Table S5), we confirmed TFAP4, HNRNPAB, MYBL2, TEAD4, and ZNF219 as MYCNA-specific dependencies (Fig. 2F; Supplementary Table S5).

Taken together, comparative analysis of the three screens identified eleven potential MYCNA-specific MR dependencies (Fig. 2G), which included both novel and previously reported transcription factors (TF) associated with the MYCNA tumors, including MYCN (40), MYBL2 (41), and E2F3 (42).

We have proposed that the stability of tumor-related phenotypes is controlled by tightly autoregulated MR protein modules (tumor checkpoints) that mechanistically regulate the transcriptional state of the cancer cell (9, 15). To test this hypothesis, we assessed the ability of MYCNA-specific MRs to mechanistically regulate each other, as well as the MYCNA-subtype transcriptional signature. To elucidate the causal control architecture of the corresponding tumor checkpoint, we performed lentiviral-mediated shRNA silencing of each of the 11 MYCNA-specific MRs that emerged from the in vitro and in vivo RNAi screens (Fig. 2G), followed by qRT-PCR analysis of all other MRs (Supplementary Fig. S6A).

To select the most likely transcriptional interactions, we focused on MRs whose transcripts were strongly downregulated (≥1.5-fold) following silencing of another MR (Fig. 3A). This analysis revealed a highly modular and hierarchical 10-MR tumor checkpoint architecture (Fig. 3B). Note that the ECSIT protein was eliminated because it was neither significantly regulating nor was regulated by any other MR. TEAD4 emerged as the most critical MYCN effector, strongly regulating 5 of the remaining 8 MYCNA tumor checkpoint MRs. Taken together, MYCN and TEAD4 emerged at the top of the regulatory hierarchy and were jointly responsible for regulating all but one (HNRNPAB) of the other MRs, whereas the latter was regulated indirectly by both MYCN and TEAD4 via MYBL2. MYCN and TEAD4 thus emerged as the core control unit of the MYCNA tumor checkpoint and thus bona fide MRs of the MYCNA subtype.

Immunoblotting assays confirmed that MYCN silencing in BE2 cells downregulated TEAD4 at the protein level (Supplementary Fig. S5E). Furthermore, chromatin immunoprecipitation (ChIP) assays in BE2 cells, using an MYCN-specific antibody, confirmed MYCN binding to the promoter region of 3 of 4 predicted targets, including TFAP4, MYBL2, and TEAD4 (Supplementary Fig. S5F).

To gain further insights into the regulatory role of TEAD4, we performed ChIP with a TEAD4-specific antibody, followed by next-generation sequencing (ChIP-seq) in BE2 cells. Phantom peak quality assessment (43) of the libraries confirmed a strong signal-to-noise ratio in the experiment (Supplementary Fig. S6B and S6C). Additional quality checks after peak calling indicated a high TF-binding sequence motif affinity (Supplementary Fig. S6D), as well as cell-type specificity (Supplementary Fig. S6E and S6F). A remarkable similarity of BE2 and SKNSH (both neuroblastoma-derived cell lines) was observed, compared with other cell types, supporting the consistency of ChIP-assay results with tissue of origin.

Significant TEAD4 peaks were identified in the promoters or enhancers of 2 of 5 predicted target MRs from the perturbational analysis (Fig. 3B), including HMGB2 and PRDM8, as well as of MYCN (Supplementary Table S6), further supporting the mechanistic nature of the tumor checkpoint architecture predicted by our perturbational analysis and suggesting a possible TEAD4–MYCN autoregulatory loop. Consistently, the strong effect of TEAD4 silencing on the remaining three genes (ZNF219, TFAP4, and MYBL2) can be explained as being mediated by indirect interactions via HMGB2 and PRDM8 (Fig. 3B).

Taken together, this analysis suggests a highly interconnected tumor checkpoint module, rich in autoregulatory loops, with TEAD4 representing the key effector of aberrant MYCN activity and the two proteins jointly regulating all other tumor checkpoint MRs, either physically or indirectly via another MR.

To further validate the functional role of the TEAD4–MYCN feedback loop as the core regulatory element of the MYCNA tumor checkpoint, we assessed the global effect of TEAD4 and MYCN silencing on the MYCNA signature. Lentiviral-mediated shRNA silencing of TEAD4 and MYCN was performed independently in BE2 cells, followed by RNA-sequencing (RNA-seq) profiling. To enrich for direct targets, we analyzed the results at 48 hours after transduction and confirmed that their silencing was not yet affecting the other protein levels (Supplementary Fig. S6G). MYCN and TEAD4 silencing both significantly reversed the MYCNA signature toward a stage 1 tumor signature (Fig. 3C) as confirmed by GSEA (P_shMYCN = 1.87E−5, Supplementary Fig. S6H; P_shTEAD4 = 1.81E−6; Supplementary Fig. S6I). Taken together, MYCN and TEAD4 differentially expressed genes comprise ∼70% of MYCNA subtype differentially expressed genes (Fig. 3C and D). However, individually, MYCN and TEAD4 regulated a highly complementary set of targets, with only a small fraction of common targets, thus suggesting a complementary role (Fig. 3D).

Next, to evaluate the biological programs controlled by these two TFs, we performed REACTOME pathway and gene ontology (GO) biological processes enrichment analysis of the overlap between genes differentially expressed following their silencing and the MYCNA gene expression signatures. MYCN silencing significantly reversed key activated (cell growth and proliferation) and inactivated (differentiation) programs in the MYCNA signature (Fig. 3E and F; Supplementary Table S7). In contrast, TEAD4 silencing specifically reversed key activated programs in the MYCNA signature (proliferation and DNA damage response; Fig. 3E and F; Supplementary Table S7). Consistent with previous reports, this suggests that whereas MYCN acts as both an activator and a repressor (44, 45), TEAD4 acts mainly as a transcriptional activator (46). This was confirmed by the dramatic difference in the overlap of TEAD4 ChIP-seq binding sites and TEAD4-activated targets (P = 4.05E−42) versus repressed targets (P = 1.9E−2; Supplementary Fig. S6J). Consistently, the ARACNe-AP–inferred MYCN regulon had an equivalent number of positively regulated and repressed targets, whereas the TEAD4 regulon consisted mostly of positively regulated targets (Supplementary Fig. S6K and S6L); GSEAs confirmed that their ARACNe-AP–inferred regulons were highly enriched in genes differentially expressed following their RNAi-mediated silencing (Supplementary Fig. S6M and S6N), thus confirming the overall validity of the regulatory model.

Consistent with its established role, MYCN-specific activated genes were highly enriched in cell-growth/metabolism programs, including protein biosynthesis, ribosomal biogenesis, rRNA processing, RNA processing and splicing, as well as cell proliferation programs (Supplementary Table S7). In contrast, genes repressed by MYCN were enriched in neuronal differentiation, actin cytoskeleton organization, axon guidance, and cell adhesion molecules, possibly mediated by miR-17-92 cluster activation (47) or by MIZ1-mediated repression (ref. 48; Supplementary Table S7). In line with these findings, we observed neurite outgrowth and neuronal differentiation upon MYCN silencing in these cells (Supplementary Fig. S5C). On the other hand, TEAD4 induced highly significant activation of proliferative and DNA damage response (DDR) programs but only in cells where MYCN was also active (Supplementary Table S7). Remarkably, both proliferative and DDR-related genes were strongly enriched in a subset of 125 genes that were directly activated by TEAD4 in the MYCNA subtype upregulated signature (Fig. 3G; Supplementary Tables S6 and S7). Potential involvement of TEAD4 in regulating DDR-related genes was especially intriguing as it was not previously associated with its canonical role as transcriptional effector of the Hippo pathway, suggesting an aberrant role in neuroblastoma; specifically, we observed activation of genes involved in ATR-mediated response to replication stress (Supplementary Fig. S7A), whose inhibition has been shown to reduce MYCN-driven neuroblastoma viability (49, 50). Consistently, we observed an increase in γ-H2Ax following TEAD4 silencing in the MYCNA cell lines BE2 and LAN1 (Supplementary Fig. S7B and S7C). Functional validation of regulation of proliferation by TEAD4 will be discussed in the following sections. Finally, we assessed the regulatory mechanism between TEAD4 and other TEAD proteins and between MYCN and TEAD family proteins. TEAD4 silencing induced significant TEAD1 and TEAD2 downregulation (Supplementary Table S7), suggesting that, within this subtype, TEAD4 is the dominant driver. In addition, the regulation of TEAD family members by MYCN and TEAD4 varied (Supplementary Table S7), suggesting complex regulatory mechanisms between these proteins.

Next, following up on the ChIP assays showing binding of MYCN in the TEAD4 proximal region and vice-versa (Supplementary Fig. S5F; Supplementary Table S6; Fig. 3B), we proceeded to further elucidate the regulatory interaction between these two proteins. Specifically, we performed time-dependent assessment of MYCN expression and protein abundance following TEAD4 silencing in BE2 cells. This study showed small but significant MYCN mRNA downregulation following TEAD4 silencing (Fig. 4B), supporting TEAD4-mediated MYCN transcriptional regulation and confirming the functional nature of the ChIP-seq findings. Yet, downregulation of MYCN at the protein level was far more prominent (Fig. 4A), suggesting a strong posttranslational regulatory interaction.

No evidence of direct TEAD4–MYCN protein–protein interaction could be detected by coimmunoprecipitation assay. Thus, to identify additional proteins that may mediate the interaction, we assessed gene expression changes of established modulators of MYCN protein turnover, such as AURKA, FBXW7, HUWE1, TRPC4AP, and CDK1 complex (51–55), following TEAD4 silencing. Both AURKA and CDK1 emerged as significantly downregulated (AURKA: 2-fold; P = 2.48E−07; CDK1: 2-fold; P = 1.38E−07; Supplementary Table S7). We further confirmed that TEAD4 regulates AURKA and CDK1 both at the protein and at the gene expression levels (Fig. 4A and B), with TEAD4 bound to CDK1's proximal region by ChIP-seq (Supplementary Table S6), thus confirming TEAD4-mediated regulation of established MYCN modulators. These interactions were confirmed in an additional MYCNA cell line, LAN1, albeit to a lesser degree (Supplementary Fig. S7D and S7E), in line with its lower enrichment in MYCNA MR signature proteins (Supplementary Fig. S4E). Consistent with these findings, treatment with cycloheximide increased MYCN turnover by 2-fold following TEAD4 silencing, compared with control cells (Fig. 4D and E), which was rescued by the proteasome inhibitor MG-132 (Fig. 4F and G).

Stabilization and degradation of both MYCN and MYC proteins require sequential phosphorylation at serine 62 and threonine 58, and the sequence around this region is conserved in both (53, 56). It has been reported that the expression of MYCN and MYC occurs in a mutually exclusive fashion by repressing each other in neuroblastoma cells (57). Therefore, we examined regulation of MYC by TEAD4 in MYCN^WT cells (SY5Y). Our data confirmed that TEAD4 also regulates MYC (Supplementary Fig. S7F and S7G). Although there are conflicting data in the literature regarding regulation of MYC by AURKA (52, 58, 59), CDK1 has been shown to regulate both MYCN and MYC (60). Furthermore, it has been reported that TEAD4 binds to the enhancer region of MYC (61), indicating transcriptional regulation as well.

To gain further insight into the interplay between TEAD4 and MYCN/MYC, we assessed their correlation in neuroblastoma patient GEPs. We observed that both TEAD4 expression and activity were increased only in MYCNA tumors or in tumors presenting overexpression of either MYCN or MYC (Fig. 5A; Supplementary Fig. S7H). This pattern is consistent across tumors from multiple histologies where TEAD4 positively correlates with MYC target hallmark activity (Supplementary Fig. S7I). TEAD4 expression correlates positively with MYCN in tumors of neuronal origin (Supplementary Fig. S7J) whereas in most other histologies TEAD4 positively correlates with MYC (Supplementary Fig. S7K).

Taken together, these data show the existence of a strong TEAD4–MYCN positive feedback loop in MYCNA subtype samples. This loop is mediated by both transcriptional and posttranslational interactions and decouples these proteins from their physiologic regulatory controls, thus inducing aberrant activity of a 10-protein tumor checkpoint identified by our analysis (Fig. 4C). Thus, aberrant TEAD4 activity is necessary to stabilize and complement the tumorigenic role of MYCN and MYC in neuroblastoma.

The observation that TEAD4 regulates MYCN/MYC and cell cycle–related programs implicates an aberrant, context-specific role of TEAD4 in high-risk neuroblastoma. To evaluate the phenotypic consequences of TEAD4 silencing in MYCN^Amp and MYCN^WT cell lines, we chose cell lines with varying degree of expression of MYCN or MYC, yet expressing comparable levels of TEAD4 protein (Supplementary Fig. S8A). TEAD4 silencing led to decrease in tumorigenic potential of MYCNA cell lines BE2 and LAN1, as shown by dramatic colony count reduction in long-term colony formation assays (Fig. 5B and C). In contrast, there was no change in colony formation in the low-MYCN-activity cell line SKNFI and only a modest decrease in SKNAS cell lines (Fig. 5B and C), which despite low MYCN levels present with higher MYC levels than SKNFI (Supplementary Fig. S8A).

To further elucidate the role of TEAD4 in regulating cell proliferation, we performed GSEA of differentially expressed genes following TEAD4 silencing in BE2 cells and observed significant enrichment of cell cycle–specific genes (P = 6.7E−05; Fig. 5D). The most downregulated genes (i.e., in the GSEA leading edge) contained almost half (49/128) of the KEGG cell-cycle pathway genes. In particular, we observed repression of several genes involved in cell-cycle progression and DNA replication. Several of these genes have a reported role in high-risk neuroblastoma pathology, often associated with MYCNA tumors. These include cyclin-dependent kinases (CDK2, CDK1, and CDC25B; refs. 53, 62, 63), cyclins (CCND1; ref. 64), E2Fs (E2F1, E2F2, and E2F3; ref. 42), DNA replication factors (PCNA, MCM7, and CDC6; refs. 65–67), checkpoint kinases (CHEK1, CHEK2, and WEE1; refs. 2, 49, 68) and components of the ubiquitin-proteasome system (SKP2; ref. 69). In addition, 19 of the 49 cell-cycle genes were shown to bind TEAD4 in their proximal region (Fig. 5D; Supplementary Table S6), suggesting direct transcriptional regulation. Further investigation of the phenotypic influence of TEAD4 on the cell cycle and proliferation of MYCNA cells showed that TEAD4 silencing induced significant accumulation of cells in G₀–G₁ with concomitant decrease of cells in S phase (Fig. 5E). Consistent with this, we observed a decrease in proliferating cells by EdUrd staining (Fig. 5F). We did not observe induction of apoptosis by Annexin V staining and only a modest increase in PARP cleavage upon TEAD4 silencing (Supplementary Fig. S8B and S8C). Collectively, these findings suggest TEAD4 as a critical component driving cellular proliferation of MYCNA cells and a novel MYC/MYCN-mediated tumor dependency in neuroblastoma.

YAP and TAZ are established TEAD4 cotranscriptional activators in the Hippo pathway. They are known to bind TEAD family proteins, including TEAD4, to promote cell proliferation, growth, and survival (70). We thus explored the potential role of YAP/TAZ in TEAD4-mediated regulation of MYCNA signature genes. First, we assessed whether differential activity of YAP/TAZ could be detected in neuroblastoma samples. Surprisingly, YAP and TAZ activity, as defined by the expression of their target genes—based on REACTOME (71) and perturbation assays (72)—was not correlated with the TEAD4 expression or MYCNA status in neuroblastoma samples (Supplementary Fig. S9A). Similarly, we observed that the active form of YAP/TAZ, as indicated by their nuclear localization, was not correlated with TEAD4 protein expression or MYCNA status in a panel of neuroblastoma cell lines (Supplementary Fig. S9B).

To further validate these findings, we performed lentiviral-mediated shRNA silencing of TAZ in YAP-null cells, BE2 (Supplementary Fig. S9C), and compared the gene expression signature of TAZ and TEAD4. First, we observed only minimal overlap in genes differentially expressed following their silencing, supporting independent transcriptional regulatory activity in MYCNA cells (Supplementary Fig. S9D–S9F). Furthermore, ARACNe-inferred TEAD4 transcriptional targets, which displayed strong enrichment in the TEAD4 knockdown signature (Supplementary Fig. S6N), showed no significant enrichment in the TAZ silencing gene expression signature (Supplementary Fig. S9G). We also confirmed that bona fide targets of YAP/TAZ such as CTGF and CYR61 (71) were significantly downregulated following TAZ silencing (Limma test, P_CTGF = 5.28E−5; P_CYR61 = 2.1E−4; Supplementary Table S7), and yet the same genes were significantly upregulated following TEAD4 silencing (Limma test, P_CTGF = 4.5E−5; P_CYR61 = 9.5E−6, Supplementary Table S7), suggesting orthogonal TAZ and TEAD4 regulatory programs. Finally, TAZ silencing in BE2 and LAN1 cells showed minimal effect on cell viability compared with TEAD4 (Supplementary Fig. S9H and S9I). Taken together, our results show that aberrant TEAD4 activity in MYCNA neuroblastoma is largely independent of YAP/TAZ expression and nuclear localization.

To further investigate the clinical relevance of TEAD4 in neuroblastoma tumors, we assessed whether its expression was associated with neuroblastoma tumor progression and outcome. Indeed, TEAD4 expression and VIPER-inferred activity were both higher in the MYCNA subtype compared with other high-risk neuroblastoma subtypes, whereas stage 1 samples exhibited the lowest levels (Kruskal–Wallis test, P = 3.05E−11; Fig. 6A and B). Similarly, TEAD4 protein staining in tumor microarrays (TMA) from 116 neuroblastoma cases showed that the high-risk group expressed significantly higher levels of the protein compared with low-risk and normal tissues. This was assessed by computing the final score [intensity of staining (0–3) by percentage of cells stained (0–100; Supplementary Table S8; Fig. 6C]. Among the high-risk group, MYCNA tumors had higher expression of TEAD4 protein compared with non-MYCNA tumors (Fig. 6D).

We then performed Cox proportional hazards analysis on NRC-cohort samples and the SEQC 498NB independent cohort (73). Both unbiased cohorts comprise all tumor risk groups and stages. We observed that TEAD4 expression alone was a strong predictor of patient survival (P_NRC = 5.64E−11; P_SEQC = 8.99E−15; Supplementary Fig. S10A and S10B). Multivariate Cox regression analysis concluded that TEAD4 is a predictor of survival independent of currently used clinical and biological variables for risk stratification (74), stage (P_NRC = 8.97E−06, HR_NRC = 5.36; P_SEQC = 1.14E−06, HR_SEQC = 2.11), MYCN amplification (P_NRC = 9.05E−04, HR_NRC = 4.54; P_SEQC = 1.35E−06, HR_SEQC = 1.65), age (P_NRC = 1.87e−08, HR_NRC = 9.73; P_SEQC = 1.67E−06, HR_SEQC = 2.43) and a combination of all three covariates in NRC (P_NRC = 3.69E−03, HR_NRC = 3.69) but not in SEQC (Supplementary Fig. S10C). We further confirmed TEAD4 as a predictor of survival independently of a meta-PCNA proliferation signature (ref. 75; P_NRC = 6.26E−3, HR_NRC = 3.3; P_SEQC = 1.10E−06, HR_SEQC = 1.54), and at least in the NRC cohort TEAD4 expression predicts survival independently of an MYCN functional signature capturing the combined activity of MYCN and MYC in neuroblastomas (ref. 24; P_NRC = 6.26E−03, HR_NRC = 3.3; Supplementary Fig. S10D).

Greater understanding of the molecular mechanisms downstream of the genetic alterations that drive high-risk neuroblastoma subtypes is critically required to facilitate development of novel therapeutic strategies and to improve overall patient survival. Consistent with the recently proposed tumor checkpoint model (15), we elucidated a modular and hierarchical 10-protein architecture, responsible for the implementation of an aggressive neuroblastoma subtype (MYCNA) associated with aberrant activity of MYCN/MYC proteins. Experimental validation of these MR proteins confirmed their enrichment in MYCNA subtype–specific essential genes and their ability to regulate each other, as well as MYCNA-signature genes, through multiple autoregulatory loops, thus establishing their role as bona fide master regulators. A novel TEAD4–MYCN positive feedback loop, mediated by both transcriptional and posttranslational interactions, emerged as the dominant regulatory structure at the top of the tumor checkpoint hierarchy and was shown to represent a critical tumor dependency.

Critically, these MR proteins were conserved in independent neuroblastoma cohorts, and aberrant activation of the tumor checkpoint they comprise was found to be consistent across the entire MYCNA subtype, independent of MYCN amplification status. This is significant in light of recent findings suggesting that progression of malignant neuroblastoma with the most aggressive phenotype is driven by aberrant MYCN or MYC activity (22, 23). We show that aberrant interaction between TEAD4 and MYCN/MYC allows cells to undergo rapid proliferation and replication stress, while simultaneously activating complementary DDR pathways, thus providing mechanistic support for this observation. Consistently, inhibition of CHK1 and WEE1, both of which emerged as TEAD4 regulated, was previously shown to exhibit strong synthetic lethality in MYC-driven tumors—including neuroblastomas with high MYCN/MYC activity—by inducing potent cytotoxic response (49, 50, 68). As MYC proteins are capable of mediating oncogene-induced replication stress and genomic instability (76), the positive feedback loop between TEAD4 and MYCN, and possibly MYC, further supports this tumor initiation and maintenance mechanism.

Surprisingly, our data show that TEAD4 activity in MYCNA-subtype neuroblastoma is not mediated by its canonical cotranscriptional effectors in the Hippo pathway, YAP and TAZ (71, 77). This is likely because the positive feedback loop with MYCN decouples it from its normal physiologic control mechanisms by inducing saturation of its expression. Although prior studies have indicated TEAD4 driving previously identified YAP/TEAD targets in neuroblastoma cells (61, 78), these studies were performed in MYCN^WT cells, further confirming that different regulatory mechanisms are at play in MYCN^Amp versus MYCN^WT cells. Interestingly, studies in mammary tumors have shown that MYC represses YAP/TAZ activity, while also showing that MYC induction decreases binding of YAP/TAZ to its bona fide targets, CTGF and CYR61, but not of TEAD4 (79); additional studies have shown YAP/TAZ-independent regulatory activity of TEAD4 (80). In contrast to the canonical view that TEAD4 lacks independent transcriptional activation (81), TEAD4 was also recently reported to have a transcriptional activation domain, supporting transcriptional regulation independent of the YAP/TAZ DNA binding domain (82). This is consistent with the regulatory activity detected by our studies when aberrant TEAD4 protein expression is achieved due to its interaction with MYCN and MYC.

Although validation was restricted to the MYCNA subtype, our analysis identified identically conserved MR protein architectures for two additional subtypes, which largely overlapped with previous classification based on genetic alterations in neuroblastomas (83), including a subtype presenting a strong mesenchymal signature but lacking hallmark genomic alterations. MR proteins for this subtype were found to be highly overlapping with those previously reported for the mesenchymal subtype of high-grade glioma (8). Additionally, conservation of these proteins in cell lines implies their role as subtype MRs independent of stromal/immune infiltration. However, high rates of stromal and immune cell infiltration in a subset of these tumors suggests that single-cell analysis may be required to further elucidate the interaction between tumor cells and stroma, including the presence of mesenchymal cells within neuroblastoma tumors highlighted by recent studies (84). Perhaps most importantly, identification of master regulator proteins responsible for the implementation and stability of high-risk neuroblastoma subtypes, which are conserved in cell-line models, provides the opportunity for the systematic identification of subtype-specific therapeutic vulnerabilities using methodologies such as OncoTreat, which were recently validated in neuroendocrine tumors (16). Because OncoTreat is NYS CLIA certified, this may further support the design of clinical trials enriched for patients representative of specific molecular subtypes.

Taking these results together, TEAD4 emerged as a highly conserved, mutation-independent tumor vulnerability for the MYCNA subtype of high-risk neuroblastoma, as well as a highly significant prognostic factor. As such, it may represent an ideal novel target for therapeutic intervention in this high-risk subtype. The results of this study are also highly consistent with our proposed model of a recurrent tumor architecture responsible for canalizing the effect of multiple genomic alterations to implement critical programs necessary for tumor initiation and stability (15). Although MR proteins were discovered in multiple prior studies (8–10, 13), their modular regulatory architecture and the existence of core autoregulatory loops had not been previously reported. These findings, combined with previous MR protein studies, suggest that similar regulatory architectures may be responsible for the implementation and stability of other cancers.

Expression profile datasets on neuroblastoma from TARGET and NRC cohorts are described in Supplementary Table S1. Further characterization of the subtypes by pathway enrichment analyses (Supplementary Experimental Procedures) are provided in Supplementary Table S2. Detailed description of the approach, clustering analyses (20), and master regulator analyses (19) is provided in Supplementary Experimental Procedures. Details of the subtypes are provided in Supplementary Table S1, neuroblastoma interactomes in Supplementary Table S3, and subtype-specific MRs in Supplementary Table S4.

All neuroblastoma cell lines were obtained from the Children's Hospital of Philadelphia (CHOP) cell line bank, the Children's Oncology Group, or ATCC. They were maintained in DMEM or RPMI-1640 supplemented with 10% FBS, 2 mmol/L l-glutamine, and antibiotics. The 293FT cells were maintained in DMEM supplemented with 10% FBS and antibiotics. We received the cell lines in 2012 and the experiments were performed between 2012 and 2017. The genomic identity of each line was routinely tested and last confirmed in 2015 using the AmpFISTR Identifier Kit (Applied Biosystems). In addition, lines were routinely tested to confirm lack of Mycoplasma contamination.

Experimental validation of top-ranked MRs of MYCNA subtype was performed in neuroblastoma cell lines by lentivirus-mediated pooled shRNA and individual shRNA screening (Sigma) and siRNA screening (Dharmacon), using cell viability as a readout (Supplementary Experimental Procedures). For the respective shRNA and siRNA sequences, see Supplementary Table S5. RNA-based analyses were performed by RT-PCR (see Supplementary Table S9 for oligonucleotide sequences) or RNA-seq analyses; ChIP assays using kits from Millipore and ChIP-seq analyses were performed as described previously (Supplementary Experimental Procedures). Details of ChIP-seq and RNA-seq analyses are provided in Supplementary Tables S6 and S7, respectively. Details of other molecular, biochemical, and cellular assays are provided in Supplementary Experimental Procedures.

Analyses of TEAD4 protein expression in neuroblastoma tumors were performed using TMAs from CHOP (Supplementary Experimental Procedures). Available histopathologic features are summarized in Supplementary Table S8. Kaplan–Meier curve analyses and Cox proportional hazards regression analyses were performed using the R “survival” package (http://cran.r-project.org/web/packages/survival/index.html; Supplementary Experimental Procedures).

The tumor genomics data from the TARGET cohort are available through the data matrix portal (https://ocg.cancer.gov/programs/target/data-matrix). The NRC expression data are available in Gene Expression Omnibus with accession codes GSE85047. Data generated through the ENCODE project, including TEAD4 ChIP-seq data, were obtained from http://genome.ucsc.edu/ENCODE/downloads.html. Additional data generated in this study, including RNA-seq profiles from BE2 cells and ChIP-seq data in BE2 cells using TEAD4 antibody, are available in Gene Expression Omnibus with accession code GSE84389.

Mice were housed in a pathogen-free animal facility. All animal studies were approved by the Institutional Animal Care and Use Committee at Columbia University (#AAAQ2459). Mice used for subcutaneous xenograft experiments were 4- to 6-week-old male and female athymic nude (Nu/Nu, Charles River Laboratories).

B.J. Abraham has ownership interest (including patents) in Syros Pharmaceuticals. M.J. Alvarez is Chief Scientific Officer at DarwinHealth. R.A. Young has ownership interest (including patents) in Syros Pharmaceuticals and is a consultant/advisory board member for the same. A. Califano is the founder and equity holder of DarwinHealth Inc., a company that has licensed some of the algorithms used in this manuscript from Columbia University. Columbia University is also an equity holder in DarwinHealth Inc. He is also a consultant/advisory board member for Tempus. No potential conflicts of interest were disclosed by the other authors.

Conception and design: P. Rajbhandari, G. Lopez, C. Capdevila, A. Iavarone, J.M. Silva, J.M. Maris, A. Califano

Development of methodology: G. Lopez, C. Capdevila, B. Salvatori, J. Yu, R. Rodriguez-Barrueco, M. Yarmarkovich, M.J. Alvarez, A. Iyer, J.S. Wei, J.M. Silva, J.M. Maris, A. Califano

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): P. Rajbhandari, G. Lopez, R. Rodriguez-Barrueco, D. Martinez, N. Weichert-Leahey, A. Iyer, K. De Preter, J. Koster, S. Asgharzadeh, R.C. Seeger, J.S. Wei, J. Khan, J. Vandesompele, P. Mestdagh, R. Versteeg, A. Iavarone, A. Lasorella, J.M. Silva, J.M. Maris

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): P. Rajbhandari, G. Lopez, J. Yu, M. Yarmarkovich, B.J. Abraham, M.J. Alvarez, J.L. Harenza, K. De Preter, J. Vandesompele, P. Mestdagh, A. Iavarone, A. Lasorella, J.M. Silva, J.M. Maris, A. Califano

Writing, review, and/or revision of the manuscript: P. Rajbhandari, G. Lopez, B.J. Abraham, J.L. Harenza, D. Oldridge, J. Vandesompele, A.T. Look, J.M. Silva, J.M. Maris, A. Califano

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): G. Lopez, B.J. Abraham, D. Oldridge, J. Koster, R. Versteeg, J.M. Maris

Study supervision: R.A. Young, A. Iavarone, A. Lasorella, J.M. Silva, J.M. Maris, A. Califano

Other (VIPER analysis for cell lines): M.J. Alvarez

This work was supported by the Cancer Target Discovery and Development Centers National Cancer Institute (NCI) U01 CA217858 (A. Califano), U01 CA217858-01 (A. Califano), and the NIH grants RC1MD004418 to the TARGET consortium and CA98543 and U10 CA180899 to the Children's Oncology Group (J.M. Maris). The computational analyses were supported by NIH grants S10OD021764 (A. Califano) and S10OD012351 (A. Califano). Additional funding sources include T32GM008798 training grant (P. Rajbhandari), A Collaborative Pediatric Cancer Research Awards Program (Rally Foundation) Independent Investigator Award (G. Lopez), NCI Outstanding Investigator Award R35CA197745-02 (A. Califano), NCI Research Centers for Cancer Systems Biology Consortium 1U54CA209997 (A. Califano), NCI-NIH Contract No. HHSN261200800001E (J.M. Maris), R35 CA220500 (J.M. Maris), R01 CA180692 (J.M. Maris), and P01 CA217959 (J.M. Maris). The authors wish to acknowledge the Neuroblastoma Research Consortium (NRC), represented by Rogier Versteeg, Jan Koster, Alexander Schramm, Johannes H. Schulte, Angelika Eggert, Luigi Varesio, Angela Rita Sementa, Jo Vandesompele, Frank Speleman, Raymond Stallings, and Ingrid Ora, for providing gene expression and DNA copy-number data from 283 primary neuroblastoma tumors. The content of this publication does not necessarily reflect the views or policies of the Department of Health and Human Services, nor does mention of trade names, commercial products, or organizations imply endorsement by the U.S. government.

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.