Introduction

Metabolic reprogramming is a hallmark of cancer^1-3 which allows rapidly proliferating tumor cells to acquire nutrients to meet their bioenergetic, biosynthetic, and redox demands⁴. One of the primary driving forces in reprograming cancer cell metabolism is the deregulated MYC family proto-oncogenes (C-MYC, MYCN, and MYCL)⁵, which are known to encode master transcriptional factors that regulate metabolic gene expression. MYC coordinates nutrient acquisition to produce ATP and key cellular building blocks that increase cell mass and promote DNA replication and cell division⁶. The increase in total RNA and protein synthesis by overactive MYC signaling leads to dysregulation of macromolecular processing machineries including the spliceosome⁷, and consequently pre-mRNA splicing^8-10, another hallmark of MYC-driven cancers^7,9-12, for the purpose of cellular stress adaptation. MYCN amplification is one of the most important biological features of high-risk neuroblastoma¹³. Transgenic mouse and zebrafish models have demonstrated that MYCN is a neuroblastoma driver^14,15. In tumors without MYCN amplification, C-MYC is overexpressed, further indicating that neuroblastoma is a MYC-driven cancer. The metabolic dependency of neuroblastoma has been widely studied by us and others^16-23. A larger number of splicing changes have also been noticed in high-stage neuroblastomas^24-26. Splicing alterations lead to a spliceosomal vulnerability that provides a new opportunity to develop transformative therapies by disrupting aberrant pre-mRNA splicing^7,9-12. We and others have shown that targeting the splicing factor RBM39 by indisulam, a “molecular glue” that bridges RBM39 to E3 ubiquitin ligase DCAF15 for proteasomal degradation, achieved an exceptional anti-tumor activity in neuroblastoma models^27,28. Disruption of spliceosome by Pladienolide B also resulted in significant anti-tumor effect in neuroblastoma models²⁶ However, how the dysregulated pre-mRNA splicing machinery and metabolism are orchestrated in MYC-driven neuroblastoma has not been well elucidated. Whether metabolism modulates the anti-cancer effect of splicing inhibition remains to be answered.

Next-generation sequencing studies have revealed only a few recurrent somatic mutations in neuroblastoma at the time of diagnosis^29,30. However, copy number alterations of chromosomal segments such as 17q gain, 1p36 or 11q23 loss frequently occur in high-risk neuroblastoma. While attempts to understand the functions of individual genes in these chromosomal segments have been reported (i.e., BIRC5³¹, PHB³², PPM1D³³, TRIM37³⁴ in 17q; ARID1A³⁵, CAMTA1³⁶, CASZ1³⁷, CHD5^38-40, KIF1Bβ^41-43, miR-34a⁴⁴, RUNX3⁴⁵ in 1p36), the biological consequences of these genetic events in MYC-driven tumors still remain largely unknown. Gain of 17q is the most frequent genetic event in high-risk neuroblastoma and is associated with MYCN amplification⁴⁶. In addition, in the transgenic MYCN mouse model of neuroblastoma, the chromosomal locus syntenic to human 17q is partially amplified⁴⁷, indicating that chromosome 17q is needed for MYC-mediated tumorigenesis.

JMJD6 is a JmjC domain–containing nuclear protein with iron- and 2-oxoglutarate–dependent dioxygenase activity⁴⁸, whose coding gene is located on chromosome 17q25. While the histone arginine demethylase activity of JMJD6 that catalyzes demethylation of H4R3me1/me2 is controversial⁴⁹, JMJD6 is a lysyl-5-hydroxylase that catalyzes 5-hydroxylation on specific lysine residues of target proteins⁵⁰. JMJD6 has pleiotropic functions in normal physiology and in cancer^51-54. We previously found that JMJD6 is essential for the survival of neuroblastoma cells (including MYCN-amplified and C-MYC– overexpressed cells)⁵⁵, which was further validated by an independent study⁵⁶, indicating that neuroblastoma has JMJD6 dependency. However, the exact mechanism of JMJD6 in MYC-driven cancers remains elusive. One study has shown that JMJD6 and BRD4 co-bind at antipause enhancers, regulating promoter-proximal pause release of a large subset of transcription units⁵⁷. By harnessing a similar mechanism, JMJD6 promotes cell survival of glioblastoma in vivo⁵⁸. These findings are particularly interesting because BRD4 occupies exceptionally large super-enhancers associated with genes, including C-MYC and MYCN^59-61, and the expression of those enhancers can be disrupted by BRD4 inhibitors, which have a potent antitumor effect^59-62. Here we show a new mechanism by which JMJD6 promotes tumorigenesis mediated by the MYC oncogene in that JMJD6 interacts with a subset of RNA binding proteins including RBM39 in neuroblastoma cells and regulates the alternative splicing of metabolic genes that are involved in mitochondrial metabolism. “Glutamine addiction” is one key feature of MYC-driven tumors. Glutaminase (GLS) is the enzyme responsible for conversion of glutamine to glutamate in the process of glutaminolysis to feed the tricarboxylic acid (TCA) cycle and has two splice isoforms, GAC (glutaminase C) and KGA (kidney-type glutaminase). We show that JMJD6 controls the alternative splicing of KGA and GAC, and, consequently, impacts the central carbon metabolism in neuroblastoma. Further we show that JMJD6 is correlated with the anti-cancer activity of indisulam, a “molecular glue” that degrades the splicing factor RBM39. The indisulam-mediated cancer cell killing is at least partly dependent on the glutamine-related metabolic pathway mediated by JMJD6. Our findings demonstrate a new mechanism by which JMJD6 coordinates metabolic programs and alternative pre-mRNA splicing, providing a rationale to target JMJD6 as a therapeutic target for MYC-driven cancers.

Results

The essential genes for neuroblastoma cell survival on chromosome 17q target pre-mRNA splicing and metabolism

An incomplete understanding of the biological consequences of chromosome 17q gain remains a barrier to the understanding of high-risk neuroblastoma. 1132 genes are located on 17q (Supplementary table 1). We surmised that some of the 17q genes are particularly important for neuroblastoma cell survival. Analysis of the cancer dependency genes in neuroblastoma cell lines screened with the Avana sgRNA library⁶³ revealed that 114 were essential to neuroblastoma (mean score <-0.4) (Figure 1a, Supplementary table 2). Protein interaction network analysis followed by functional annotation revealed that proteins encoded by these 114 essential genes formed distinct but interconnected modules including RNA splicing (i.e., SRSF2, DDX5, DDX42, DHX8), mitochondrial metabolism (i.e., NDUFA8, COX11, SLC25A10, SLC35B1), protein homeostasis (i.e., UBE2O, PSMB3, PSMC5), DNA repair (i.e., BRIP1, BRCA1, RAD51C, RAD51D) and transcriptional regulation (i.e., PHF12, CBX1, SMARCE1, MED1), as well as endocytosis (i.e., CHMP6, CTLC, EPN3, HGS, SNF8, VPS25) (Figure 1b). Using these 114 genes as a signature, we found that 81 of them were highly expressed in high-risk neuroblastomas, which were enriched with MYCN amplification (Figure 1c). Correspondingly, neuroblastomas with high expression levels of this gene signature were associated with a poorer event-free and overall survival of patients (Figure 1d, e). These data demonstrate that 17q genes are involved in essential biological processes and highly expressed in high-risk neuroblastomas.

Figure 1.

JMJD6 is required for neuroblastoma growth

JMJD6 was among these 114 essential genes. To understand the role of JMJD6, we examined the genetic features of JMJD6 in neuroblastoma and other types of cancers. Among the genes encoding JmjC-domain containing proteins, JMJD6 was the only one that was frequently amplified in neuroblastoma (Figure 2a). High JMJD6 expression was associated with poor event-free outcome, as shown by Kaplan-Meier analysis (Figure 2b). To examine whether JMJD6 amplification is limited to specific tumor types, we explored genomic data from different cancers using the cBioportal program. JMJD6 was amplified across multiple types of adult cancers such as breast and liver cancer (Supplementary Figure 1a), and correlated with worse relapse-free survival (Supplementary Figure 1b). We further compared the RNA-seq expression of JMJD6 in 2337 samples across over 20 pediatric cancer subtypes and found that JMJD6 showed the highest expression levels in neuroblastoma (Supplementary Figure 1c), suggesting that JMJD6 might be particularly important in neuroblastoma. We validated this hypothesis using shRNA knockdown of JMJD6 in MYCN amplified cells (BE2C, SIMA, KELLY, IMR32) and non-MYCN amplified cells (SK-N-AS and CHLA20). The results showed that loss of JMJD6 greatly reduced the colony numbers in all tested cell lines (Figure 2c), demonstrating that JMJD6 is essential to neuroblastoma cells regardless of MYCN amplification. Neuroblastic tumors comprise a histologic spectrum that ranges from less-differentiated neuroblastoma to well-differentiated ganglioneuroma. The extent of differentiation in the tumor cells is correlated with prognostic significance⁶⁴. We noticed that the loss of JMJD6 led to neurite outgrowth (Supplementary Figure 1d), a unique structure of neuroblastoma cells differentiating in vitro. This morphologic change suggests that JMJD6 is required to maintain cancer cell stemness. Lastly, we validated that JMJD6 is essential to neuroblastoma growth in MYCN amplified (BE2C) and C-MYC overexpressed (SK-N-AS) xenograft models (Figure 2d, e). Taken together, these data demonstrate that loss of JMJD6 function impedes neuroblastoma cell survival and tumor growth.

Figure 2.

JMJD6 promotes MYC-mediated cellular transformation

Next, we investigated whether gain of function of JMJD6 could facilitate MYC-mediated oncogenic transformation. To test this, we used a NIH3T3 transformation assay that provides a straightforward method to assess the transforming potential of an oncogene, which may lead to morphological transformation and loss of contact inhibition, a typical feature of cellular transformation. Like the GFP control (Supplementary Figure 2a), we found that NIH3T3 cells with overexpressed JMJD6 stopped proliferation after being confluent (Supplementary Figure 2b), indicating JMJD6 alone is unable to transform NIH3T3 cells. However, overexpression of MYCN induced foci formation with enhanced cell death (Supplementary Figure 2c). Importantly, NIH3T3 cells with overexpressed MYCN and JMJD6 lost contact inhibition, accompanied with morphological change, and formed larger foci (Supplementary Figure 2d), indicating that JMJD6 enhances MYCN activity to transform NIH3T3 cells. Interestingly, MYCN alone also reprogramed metabolism of NIH3T3 cells as shown by the color change of the media, indicative of acidic pH change, which was largely rescued by co-expression of JMJD6 (Supplementary Figure 2e), suggesting that cells with enhanced lactate production by MYCN were directed to oxidative phosphorylation by JMJD6. It is believed that the cell of origin of neuroblastoma is the progeny of neural crest cells. We therefore tested the role of JMJD6 in MYC-mediated transformation using JoMa1, a cell line derived from murine neural crest, by transducing GFP, JMJD6, MYCN and JMJD6/MYCN (Figure 2f). While JMJD6 showed no difference from GFP control in regulating cell proliferation, MYCN slightly enhanced cell proliferation (Figure 2g). However, the combination of JMJD6 and MYCN remarkably increased cell proliferation, mirrored by the colony formation assay which showed JMJD6/MYCN induced rapid growth of colonies with distinct transformation morphology (Figure 2g, h). Implantation of each group into immune-deficient mice led to tumor development of MYCN and JMJD6/MYCN groups (Figure 2i). However, JMJD6/MYCN group tumors appeared to grow faster than the MYCN alone tumors. Taken together, these data indicate that JMJD6 enhances the MYC-mediated transformation, demonstrating the oncogenic role of JMJD6.

JMJD6 regulates pathways engaged in pre-mRNA splicing and mitochondrial biogenesis in neuroblastoma

We surmised that loss of function of genes in the same functional module/pathway may induce similar effects across different cell lineages, which in turn corroborates the hypothesis that JMJD6 is a player in that signaling pathway. To test this hypothesis, we analyzed the dependency correlation of JMJD6 knockout and other genes (defined as co-dependency genes if they are positively correlated), by using the DepMap data (https://depmap.org) that includes genome-wide knockout in 1027 cell lines of more than 20 cancer types (Supplementary Table 3), followed by pathway enrichment. The data showed that JMJD6 co-dependency genes were significantly and positively correlated with spliceosome/mRNA splicing (i.e., RBM39, SF3B1), ubiquitin-mediated proteolysis and endocytosis and a number of 17q25 genes, (Figure 3a, c), which mirrored the pathway network of 17q essential genes in neuroblastoma (Figure 1b). JMJD6 knockout was negatively correlated with the knockout of genes housed at chromosome 1p, which is frequently deleted in high-risk neuroblastoma, and oxidative phosphorylation as well as protein translation (Figure 3b, d).

Figure 3.

BRD4 is known to regulate MYC expression^61,65. Previous studies have shown that JMJD6 and BRD4 interact to regulate gene transcription^56,58,66, suggesting that JMJD6 might directly modulate MYC expression. To assess this possibility, we knocked down JMJD6 in neuroblastoma cell lines BE2(C) and SK-N-AS, which express MYCN and C-MYC, respectively, for RNA-seq analysis. The sequencing data were analyzed for differential gene expression (Supplementary Table 4). Interestingly, loss of JMJD6 showed minimal impact on expression of MYCN in BE2C cells or C-MYC in SK-N-AS cells (Supplementary Figure 3a), and western blot analysis did not show alteration of MYCN expression although C-MYC protein was slightly downregulated by loss of JMJD6 (Supplementary Figure 3b). However, BRD4 inhibitors drastically inhibited both MYCN and C-MYC expression in neuroblastoma cells^59-61, suggesting that JMJD6 inhibition has a distinct effect from the BRD4 inhibition. Gene set enrichment analysis for pathway engagement for the genes commonly downregulated or upregulated in both cell lines revealed that loss of JMJD6 most significantly repressed the expression of genes involved in pre-mRNA splicing, histones, and cell cycle G1/S checkpoint (Figure 3e), and enhanced the pathways involved in mitochondrial functions and heat shock response (Figure 3e). Interestingly, the genes transcribed from the mitochondrial genome were elevated in both cell lines (Supplementary Figure 3c), suggesting that JMJD6 directly or indirectly regulates the transcription of mitochondrial genome. These data are consistent with the co-dependency pathways of JMJD6 (Figure 3a-d). Depletion of JMJD6 in both cell lines led to the downregulation of MYC signaling pathways (Supplementary Figure 4), although ranked behind the pathways of splicing and metabolism, suggesting that the MYC pathways are not primarily regulated by JMJD6. These data indicate that JMJD6 does not regulate the gene expression of the MYC family of transcription factors but might indirectly regulate the MYC pathway.

To verify if JMJD6 regulates pre-mRNA splicing, we analyzed the RNA-seq using two algorithms. The first one is event-based analysis to identify the altered exon splicing (Figure 3f). We found that knockdown of JMJD6 dominantly affects the exon skipping (SE) although other splicing events were also altered, albeit with a much smaller number (Figure 3f). Pathway analysis of common events in both BE2C and SK-N-AS cells demonstrated that genes involved in metabolism and splicing are the most significantly affected by loss of function of JMJD6 (Figure 3g). Using the second algorithm of RNA splicing analysis previously developed (Figure 3h) that allows discovery of new isoforms of genes generated through alternative splicing⁶⁷, we identified 580 genes in BE2C cells and 1018 genes in SK-N-AS cells undergoing alternative splicing after JMJD6 knockdown, 133 of which were shared by both (Figure 3h, Supplementary Fig. 5, Supplementary Table 5). The alternatively spliced genes were involved in a variety of pathways (Supplementary Table 6). KEGG pathway enrichment analysis of the 133 commonly alternatively spliced genes showed that only metabolic pathway genes were significantly enriched (Figure 3f), most of which are involved in mitochondrial bioenergetics and folate metabolism (Figure 3h). Collectively, these data demonstrate that JMJD6 regulates RNA splicing of genes engaged in mitochondrial metabolism, being one of the key mediators of the 17q locus activity in neuroblastoma.

JMJD6 regulates alternative splicing of glutaminolysis gene GLS

Overactive MYC signaling leads to altered macromolecular processing machineries in response to an increase in total RNA and protein synthesis⁷. MYC is also a master regulator of cancer metabolism involved in ribosomal and mitochondrial biogenesis, glucose and glutamine metabolism and lipid synthesis, leading to the acquisition of bioenergetic substrates enabling the cancer cell to grow and proliferate^68-70. “Glutamine addiction” is one feature of MYC-driven cancer⁷¹. The pre-mRNA splicing altered by JMJD6 knockdown included GLS, the key enzyme of glutaminolysis, prompting us to investigate the GLS splicing mediated by JMJD6 knockdown. GLS is known to catalyze the conversion of glutamine to glutamate, and is alternatively spliced to form two isoforms, GAC and KGA⁷², with different cellular localization and catalytic capacities⁷³. The GAC isoform is more frequently upregulated in cancer cells than KGA⁷⁴, and has been shown to be regulated by MYC^75-77, leading to a “glutamine addiction” phenotype in MYC-driven tumors⁷¹. We found that loss of JMJD6 led to a splicing switch from the GAC isoform (with exons 1-15) to the KGA isoform (with exons 1-14 and 16-19) (Figure 4a), which was further confirmed by isoform-specific RT PCR (Figure 4b). We then investigated the expression of GAC/KGA at the protein levels after JMJD6 knockdown. Western blot showed that the KGA isoform was increased after JMJD6 knockdown in all three tested cell lines, SK-N-AS, BE2C and SIMA (Figure 4c). Then we further validated the JMJD6 effect on GLS isoform expression by using a luciferase reporter that indicates the isoforms of GAC and KGA. Indeed, JMJD6 knockdown significantly increased the expression of the KGA reporter (Figure 4d). RNA-immunoprecipitation showed that JMJD6 bound to GLS RNA (Figure 4e), suggesting that JMJD6 may directly regulate the splicing of GLS. We reasoned that if the regulation of GLS splicing by JMJD6 was a bone fide mechanism, the expression levels of JMJD6 would correlate with the levels of GAC/KGA in tumors. Indeed, JMJD6 was positively correlated with GAC and negatively correlated with KGA in two independent neuroblastoma cohorts (Figure 4f), supporting the hypothesis that JMJD6 is required to maintain the high ratio of GAC/KGA in cancer cells by controlling their alternative splicing. Clinical relevance of GAC and KGA in neuroblastoma was evidenced by the findings that high GAC was associated with a worse event-free survival and high KGA was associated with a better event-free survival (Figure 4g), suggesting that the GAC/KGA ratio may play a role in cancer progression. However, introduction of either GAC or KGA in BE2C cells or SKNAS cells promoted colony formation (Figure 4h, i), indicating that enhanced glutaminolysis by either GAC or KGA overexpression is pro-proliferative by gaining more ATPs.

Figure 4.

JMJD6 forms an interaction network with proteins involved in splicing and protein synthesis

To understand the mechanism of JMJD6 in regulating splicing in neuroblastoma, we performed an unbiased identification of JMJD6 interacting partners by introducing a flag-tagged JMJD6 into SK-N-AS and BE2C cells, followed by immunoprecipitation to pull down the JMJD6 associated complex and protein identification with mass spectrometry (Figure 5a, Supplementary Table 7). We found that JMJD6 mainly interacted with two classes of proteins which are involved in RNA splicing and protein synthesis in both cell lines, respectively (Figure 5a, Supplementary Figure 6). We then validated the interactions of JMJD6 with splicing factors using immunoprecipitation and western blot, and demonstrated JMJD6 formed a complex with these RNA binding proteins, including RBM39 (Figure 5b), a therapeutic target of high-risk neuroblastoma²⁷. Since JMJD6 also interacted with several molecules involved in protein translation, we investigated if JMJD6 also regulates protein synthesis by using an approach, Click-IT AHA to label the newly synthesized proteins, followed by western blot assessment (Figure 5c). Interestingly, overexpression of JMJD6 greatly reduced total protein synthesis (Figure 5c), suggesting that JMJD6 may antagonize protein production.

Figure 5.

Then we determined if the splicing factors interacting with JMJD6 also regulate GLS isoform expression. Among the splicing factors with which JMJD6 interacted, CPSF6 has been previously shown to regulate the alternative splicing of GLS⁷⁸. We validated the function of CPSF6 in neuroblastoma cells and found that, indeed, loss of CPSF6 led to a dramatic switch from GAC to KGA in BE2C and SK-N-AS cells (Figure 5d, 5e). Previous studies showed that JMJD6 and U2AF2 (U2AF65) interact to regulate splicing^50,79. We also found that knockdown of U2AF2 resulted in a similar phenotype to JMJD6 knockdown in that the expression of KGA isoform was greatly increased in both cell lines (Figure 5d, 5e). These data collectively support the functions of JMJD6 in regulating the splicing of metabolic genes and protein homeostasis in MYC-driven neuroblastoma.

JMJD6 regulates production of TCA intermediates and nucleoside triphosphate

To further dissect the biological consequences of loss of function of JMJD6, we created stable JMJD6 knockout clones (Supplementary Figure 7) and defined the metabolite spectrum affected by loss of function by using liquid chromatography with tandem mass spectrometry (LC-MS/MS). JMJD6 knockout greatly reduced the production of tricarboxylic acid cycle (TCA) intermediates (i.e., oxoglutarate, fumarate) and nucleoside triphosphate (ATP, CTP, GTP) (Figure 6a), indicating that JMJD6 is a key bioenergetics regulator in cancer cells. Pathway analysis revealed that the reduced metabolites were involved in the Warburg effect, TCA, pentose phosphate pathway, and mitochondrial electron transport chain (Figure 6b), all of which are critical for providing cancer cell bioenergetics for proliferation and survival. Oxoglutarate (α-ketoglutarate) and fumarate are downstream products of glutaminolysis (Figure 6c). We reasoned that cellular metabolites such as glutamate and oxoglutarate may predict the cytotoxic effect of loss-of-function of JMJD6. If cells have higher levels of glutamate and oxoglutarate, they might be less dependent on JMJD6 due to their higher capacity of buffering against reduced glutaminolysis. To test this hypothesis, we used DepMap data that included 225 metabolites in 928 cancer cell lines from over 20 lineages⁸⁰, and analyzed the correlation of each metabolite with JMJD6 gene dependency. The data showed that cells with high levels of AMP, glutamate, alanine, 2-hydroxyglutarate and 2-oxoglutarate were more resistant to JMJD6 knockout, while cells with high levels of lactose and sucrose were more sensitive to JMJD6 knockout (Figure 6d, Supplementary Table 8). High levels of AMP activate AMP kinase (AMPK), consequently leading to enhanced fatty acid oxidation to stimulate ATP production while alanine can be converted to pyruvate to provide acetyl-CoA to fuel the TCA cycle (Figure 6c). Therefore, high levels of AMP and alanine may provide cells alternative bioenergetics sources for survival. These data further indicate that JMDJ6 function is wired into the regulation of mitochondrial metabolism.

Figure 6.

JMJD6 determines the efficacy of indisulam, a molecular glue degrading splicing factor RBM39

Dysregulated splicing as a vulnerability of MYC-driven cancers provides a rationale to target neuroblastoma by using splicing inhibitors as a therapeutic approach. We and others have recently reported that indisulam, a “molecular glue” that selectively degrades the splicing factor RBM39, is exceptionally effective at causing tumor regression in multiple high-risk neuroblastoma models without overt toxicity^27,28, suggesting indisulam has translational potential. Understanding the factors determining the efficacy of indisulam or any other drug is critical for developing precision therapy, combination therapy or preventing therapy resistance. In addition to complexing together (Figure 5a, b), JMJD6 and RBM39 exhibit significant correlation of co-dependency in cancer cells, namely, cancer cells have similar dependency on JMJD6 and RBM39 for survival (Supplementary Table 3, Figure 7a). These data indicate that JMJD6 may play a role in modulating the effect of indisulam. To test this hypothesis, we performed GSEA to identify the differential pathways between indisulam-sensitive and indisulam-less sensitive neuroblastoma cells. It turned out that histone lysine demethylase (HDM) genes including JMJD6 are the in the only gene signature that is significantly enriched in indisulam-sensitive cells (Figure 7b, c, d). Indeed, knockout of JMJD6 led to partial but significant resistance to indisulam treatment of SK-N-AS cells (Figure 7e-g) and BE2C cells (Figure 7h-j), supporting that cells with high JMJD6 expression are more dependent on RBM39. Since we found that JMJD6 plays a key role in modulating glutaminolysis, we tested if expression of GAC or KGA could affect the activity of indisulam. Interestingly, overexpression of either GAC and KGA renders BE2C and SK-N-AS cells more resistant to indisulam treatment (Figure 7k-n), suggesting that enhanced glutaminolysis may confer therapeutic resistance to spliceosome inhibition.

Figure 7.

Discussion

MYC is an oncogenic driver of many types of cancer and plays a pivotal role in regulating glycolysis, glutaminolysis, nucleotide and lipid synthesis, and ribosome and mitochondrial biogenesis⁶. Recent studies have revealed that there is also an interplay between MYC and pre-mRNA splicing machinery^7-10,81-84. Pre-mRNA splicing is an essential biological process catalyzed by the spliceosome to produce mature mRNAs^85,86. Over 90% of multiexon genes in the human genome undergo alternative splicing^87,88, which significantly expands the diversity of the proteome and consequently impacts various biological functions^89,90. In this study, we found that the chromosome 17q locus, which is frequently gained in MYC-driven cancers, houses numerous genes essential to cancer cell survival that are implicated in pre-mRNA splicing, ribosome and mitochondrial biogenesis and other biological functions. Particularly, JMJD6, which is located on 17q25 and is amplified in a number of cancer types, physically interacts with a subset of splicing factors such as RBM39 and regulates the alternative splicing of metabolic genes. Depletion of JMJD6 inhibits cancer cell proliferation and impedes tumor growth while overexpression of JMJD6 promotes MYC-mediated tumorigenesis, suggesting that JMJD6 and potentially other 17q genes have oncogenic functions in cellular transformation. Previous studies suggest that JMJD6 and BRD4 interact to regulate gene transcription^56,66. However, our unbiased identification of the JMJD6 interactome only identified a subset of proteins involved in mRNA splicing and protein translation in neuroblastoma cells, suggesting that JMD6 may predominantly regulate protein homeostasis to facilitate MYC-mediated transformation. Nevertheless, we cannot exclude the possibility that our immunoprecipitation conditions were too harsh, leading to dissociation of proteins that loosely or dynamically bind to JMJD6.

MYC-induced metabolic reprogramming triggers cellular dependency on exogenous glutamine as a source of carbon for mitochondrial membrane potential maintenance and macromolecular synthesis², leading to “glutamine addiction”². Glutaminolysis is a process by which GLS and GLS2 convert glutamine to glutamate, which is, in turn, converted by glutamate dehydrogenase or transaminase to 2-oxoglutarate that is further catabolized in the TCA cycle. Additionally, glutamate is a substrate for production of glutathione, an important antioxidant. We previously showed that neuroblastoma relies on MYCN-induced glutaminolysis for survival²⁰. In this study, our RNA-seq barely detected the expression of GLS2 in the neuroblastoma cell models we used, indicating that GLS is the major enzyme that catalyzes glutamine in these model systems. GLS has two isoforms, GAC and KGA, resulting from alternative splicing. KGA is mainly localized in the cytoplasm while GAC is localized in mitochondria and has a higher basal activity⁷³. GAC mRNA levels strongly correlate with the conversion of glutamine to glutamate, as a proxy for GAC activity⁹¹. The positive correlation of JMJD6 and GAC suggests that the JMJD6-high tumors have enhanced glutaminolysis activity. CSPF6 and the noncoding RNA CCAT2 have been reported to regulate the splicing of GLS isoforms^78,92. Interestingly, we found that depletion of JMJD6 leads to a GLS isoform switch from GAC to KGA, indicating that JMJD6 is involved in alternative splicing of GLS. We further found that JMJD6 physically interacts with CPSF6 in a splicing network and validated that loss of CPSF6 results in remarkable induction of the KGA isoform. We further validated that other splicing factors such as U2AF2 also interact with JMJD6 to regulate the GLS isoform switch. These data indicate that cancer cells can adjust metabolism through alternative splicing to produce enzymes with distinct subcellular localization and activity that promote cellular transformation or progression of an oncogenic phenotype. The cooperation of JMJD6 and MYC in cellular transformation further supports the hypothesis that JMJD6 is needed for metabolic reprogramming triggered by MYC. However, overexpression of either GAC or KGA promotes cell proliferation, suggesting that the switching of KGA/GAC is a cellular fitness mechanism in response to interruption of the spliceosome by adjusting the metabolic rate. Within the tumor microenvironment (i.e., replete and deplete oxygen and nutrient supply) GLS activity is possibly finely tuned through splicing mechanism for adaption.

Additionally, we found that JMJD6 physically interacts with a subset of ribosomal proteins that are responsible for protein translation. Interestingly, overexpression of JMJD6 reduces global protein synthesis. A recent study showed that MYC overactivation leads to proteotoxic stress in cells by enhancing global protein synthesis, consequently causing cell death⁹³. The increased global protein synthesis by MYC needs to be buffered through loss of DDX3X, a regulator of ribosome biogenesis and global protein synthesis, for lymphomagenesis⁹³. Our findings suggest that, besides the functions in regulating alternative splicing for metabolism, JMJD6 is involved in MYC-mediated cell transformation by buffering unwanted proteotoxic stress due to high rate of protein synthesis induced by MYC (Figure 8).

Figure 8.

Neuroblastoma is responsible for as much as 15% of childhood cancer mortality⁹⁴. With current intensive multimodal therapies, 5-year survival rates for high-risk patients remain less than 50%^95-98. In addition, survivors of high-risk disease have a significant risk of developing long-term side effects including subsequent malignant neoplasms due to cytotoxic chemotherapy and radiotherapy^99,100. Unfortunately, developing effective precision therapies against high-risk neuroblastoma has been challenging due to the lack of targetable recurrent mutations in neuroblastoma^29,30,101. Our previous study showed that indisulam, the splicing inhibitor that targets RBM39, a JMJD6-interacting partner, induced a durable complete response in multiple high-risk neuroblastoma models, supporting its potential use in future clinical trials. Our current study showed that JMJD6 expression is positively correlated with the effect of indisulam, and knockout of JMJD6 confers resistance to indisulam treatment. In line with the biological functions of JMDJ6 in regulating GLS isoform expression and mitochondrial metabolism, overexpression of GAC or KGA also caused resistance to indisulam treatment. These data indicate that JMJD6 could serve as a biomarker that predicts response to indisulam or other splicing inhibitors.

Materials and Methods

Cell lines

KELLY, SIMA, BE2C, IMR32, SK-N-AS, CHLA20 were cultured in 1X RPMI1640 (Corning, 15-040-CV) supplemented with 10% Fetal Bovine Serum (Sigma-Aldrich, F2442), 1% L-Glutamine (Corning, A2916801). NIH3T3 and 293T cells were cultured in 1X DMEM supplemented with 10% Fetal Bovine Serum (Sigma-Aldrich, F2442), 1% L-Glutamine (Corning, A2916801). All cells were maintained at 37 °C in an atmosphere of 5% CO². JoMa1 cells were kindly provided by Dr. Schulte (Department of Pediatric Oncology and Hematology, University Children’s Hospital Essen, Essen, Germany) were cultured in NCC Medium: DMEM (4.5 mg/ml Glucose, L-Glutamine, Pyruvate): Ham’s F12 (1:1) was supplemented with: 1% N2-Supplement (Invitrogen, no. 17502-048), 2% B27-Supplement (Invitrogen, no. 17504-044), 10ng/ml EGF (Invitrogen), 1ng/ml FGF (Invitrogen), 100U/ml Penicillin–Streptomycin (Invitrogen) and 10% Chick-Embryo-Extract (Gemini Bio-Products, CA). Neural crest culture medium was supplemented with 200 nM 4-OH-tamoxifen (Sigma no. H7904) in routine culture to ensure nuclear localization of c-MycERT and JoMa1 cell proliferation. JoMa1 cells were grown on cell culture flask/dish coated with fibronectin, NCC-medium supplemented with 200nm 4-OHT was changed daily. Cells were passaged after 3–4 days in culture when 70% confluence was reached (4 X 10⁶ cells/10 cm dish).

All human-derived cell lines were validated by short tandem repeat (STR) profiling using PowerPlex® 16 HS System (Promega) once a month. Additionally, a polymerase chain reaction (PCR)-based method was used to screen for mycoplasma once a month employing the LookOut® Mycoplasma PCR Detection Kit (MP0035, Sigma-Aldrich) and JumpStart™ Taq DNA Polymerase (D9307, Sigma-Aldrich) to ensure cells were free of mycoplasma contamination.

Antibodies

GAPDH (Cell Signaling Technology, 5174s, Rabbit antibody), MYCN (Santa Cruz Biotechnology, 53993, Mouse antibody), FLAG (Sigma, F1804, Mouse antibody), Biotin (Bethyl Laboratories, A150-109A Rabbit), ACTIN (Sigma, A3854, mouse antibody), PUF60 (Thermo Fisher, PA5-21411, Rabbit antibody), U2AF2 (Novus Biologicals, NBP2-04140, Rabbit antibody), CPSF6 (Bethyl Laboratories, 357A,Rabbit antibody), DHX40 (Novus Biologicals, NBP1-91834, Rabbit antibody), DHX8 (Abcam, AB181074, Rabbit antibody), LUC7L1 (Novus Biologicals, NBP2-56401, Rabbit antibody), LUC7L2 (Novus Biologicals, NBP2-33621, Rabbit antibody), LUC7L3 (Novus Biologicals, NBP1-88053, Rabbit antibody), RBM39 (ATLAS, HPA001519, Rabbit antibody), GLS (KGA-specific), (Proteintech,20170-1-AP, Rabbit antibody), GLS(GAC-Specific) (Proteintech, 19958-1-AP, Rabbit antibody), JMJD6 (ATLAS, HAP059156, Rabbit antibody), JMJD6 (Santa Cruz biotechnology, sc-28348, Mouse antibody).

Retroviral plasmids and retrovirus packaging

MSCV-IRES-GFP and MSCV-IRES-mCherry were obtained from St Jude Vector Core. Human JMJD6 and murine MYCN were subcloned into MSCV-IRES-GFP and MSCV-IRES-mCherry, respectively. The MSCV-CMV-CMV-Flag-HA-JMJD6 was purchased from Addgene (Addgene # plasmid 31358). The retrovirus packaging was done as described in the following procedure. Briefly, HEK93T cells were transfected with viral vectors by combining 5μg of target vector, 4.4μg of pMD-old-gag-pol, and 0.6μg of VSV-G plasmids in 400uL of DMEM without serum or L-glutamine. PEIpro transfection reagent (Polyplus 115-010) was added at 2:1 (PEIpro μL: μg of plasmid) per 100mm dish of cells and mixed well, and incubated at RT for at least 20 minutes, prior to adding cells. The following day, fresh medium was added to cells. For 3-4 days, viral media was harvested and replaced twice per day. Viral media was centrifuged at 1500RPM for 10 minutes and filtered through a 0.45um vacuum filter. Virus was concentrated by ultracentrifugation at 28.5kRPM for 2 hours at 4C, aspirated, and resuspended in either OptiMEM or PBS, aliquoted, and frozen at -80C until use. Wasie, add information for GAC/KGA here.

siRNA transfection

25uM of each siRNA oligo was resuspended in 500uL of prewarmed Opti-MEM, reduced serum medium (Gibco Life technologies #31985-070) in 6 well plates. To each well, 7μL of RNAiMax (Invitrogen Lipofectamine RNAiMAX transfection reagent 13778100) was added, mixed, and left at room temperature for 10 minutes. After incubation, 100,000 cells of each indicated cell line were added to each well in a total of 2mL volume with RPMI medium supplemented with 10% FBS. JMJD6 siRNA#43, 5-CCAAAGUUAUCAAGGAAA-3; JMJD6 siRNA#45, 5-CAGUGAAGAUGAAGAUGAA-3. U2AF2 siRNA#1 AGAAGAAGAAGGUCCGU; U2AF2 siRNA#2 GUGGCAGUUUCAUAUUUG. CPSF6 siRNA#1 GGAUCACCUUCCAAGACA. CPSF6 siRNA#2 AGAACCGUCAUGACGAUU.

SDS-PAGE and Western blot

Cells were washed twice with ice-cold phosphate-buffered saline (PBS) and directly lysed on ice with 2X sample loading buffer (0.1 M Tris HCl [pH 6.8], 200 mM dithiothreitol [DTT], 0.01% bromophenol blue, 4% sodium dodecyl sulfate [SDS] and 20% glycerol). On ice, cell lysates were sonicated once with a 5 second bursts at 40% amplitude output (Sonics, VIBRA CELL) followed by 25 minutes heating at 95 °C. After the cell lysates were centrifuged at 13,000 × g at room temperature for 2 mins, 10-20 µl of the cell lysates were separated on 4-15% Mini-PROTEAN® TGX™ Stain-Free^TM Protein Gels from Bio-Rad and transferred to methanol-soaked polyvinylidene difluoride (PVDF) membranes (Millipore). Lysates for RBM39 G268V mutant cell lines and DCAF15 genetically modified cells were generated as previously described¹⁰². Membranes were blocked in PBS buffer supplemented with 0.1% TWEEN 20 and 5% skim milk (PBS-T) and incubated for 1 hour at room temperature under gentle horizontal shaking. Membranes were incubated overnight at 4 °C with the primary antibodies. The next day, membranes were washed 3 times (for 5 minutes) with PBS-T at room temperature. Protected from light, membranes were then incubated with goat anti-mouse or goat anti-rabbit HRP-conjugated secondary antibodies (1:5,000) for 1 hour at room temperature, followed by three 5-minite washes with PBS-T at room temperature. Lastly, membranes were incubated for 1 minute at room temperature with SuperSignal West Pico PLUS Chemiluminescent Substrate (34580, Thermo Fisher Scientific) and the bound antigen-antibody complexes were visualized using Odyssey Fc Imaging System (LI-COR Corp., Lincoln, NE).

RNA extraction and RT-PCR isoforms of GLS

RNA was extracted using RNeasy® Plus Mini Kit (Qiagen, reference # 74136) following the manufacturer’s protocol. cDNA was prepared in 20ul reaction from 500ng of total RNA using Superscript™ IV First Strand Synthesis System (Invitrogen, reference # 1809105) kit. Real-time PCR reactions were run in triplicates (n=3) in the 7500 Real-time PCR system by Applied Biosystems (Thermo Fisher Scientific) using power SYBR Green PCR master mix (Applied Biosystems, reference # 4367660). ΔΔCT methods were applied to analyze the results. The following primers were used to perform the quantitative Real-time PCR-GAPDH (Forward: AACGGGAAGCTTGTCATCAATGGAAA, Reverse: GCATCAGCAGAGGGGGCAGAG), GAC (Forward: GAGGTGCTGGCCAAAAAGCCT, Reverse: AGGCATTCGGTTGCCCAAACT), KGA (Forward: CTGCAGAGGGTCATGTTGAA, Reverse: ATCCATGGGAGTGTTATTCCA).

Lentiviral packaging of pLenti and shRNA

The GAC and KGA cDNAs were synthesized by Genscript company and cloned into pLenti vector. The TRC lentiviral-based shRNA knockdown plasmids for JMJD6 were purchased from Horizon Discovery (sh#46: RHS3979-201781036, TTAAACCAGGTAATAGCTTCG; sh#47: RHS3979-201781037, ATCTTCACTGAGTAGCCATCG) The lentiviral shJMJD6 and shControl (pLKO.1) particles were packaged by Vector Lab at St Jude. Briefly, HEK293T cells were transfected with shRNA constructs and helper plasmids (pCAG-kGP1-1R, pCAG4-RTR2, and pHDM-G). The 48- and 72-hr post-transfection replication-incompetent lentiviral particles were harvested and transduced into cells with 8μg/ml of polybrene. 48 hours later, 1 μg/ml of puromycin was added for selection for additional 48 hours before injection into mice or immunoblotting.

JMJD6 CRISPR KO method

Genetically modified neuroblastoma cells were generated by using CRISPR-Cas9 technology. Briefly, 400,000 NB cells were transiently co-transfected with 100pmol of chemically modified gRNA (GGACTCTGGAGCGCCTAAAA) (Synthego), 33pmol of Cas9 protein (St. Jude Protein Production Core), 200ng of pMaxGFP (Lonza), and, using solution P3 and program DS-150 in small cuvettes according to the manufacturer’s recommended protocol. Five days post-nucleofection, cells were sorted for GFP+ (transfected) cells and plated as single cells into 96-well plates. Cells were clonally expanded and screened for the desired modification using targeted next generation sequencing followed by analysis with CRIS.py (https://pubmed.ncbi.nlm.nih.gov/30862905/).

Colony formation for JoMa1 cells

Matrigel was kept at 4°C to being liquified for 6 hours. 50μL of Matrigel per 1 square centimeter area was added to 24-well plate without air bubble. The 24-well plate was kept at 37°C in cell culture incubator till it was solidified. 200 of JoMa1 cells transduced with GFP, JMJD6, MYCN, MYCN+JMJD6 in DMEM:F12 enriched media without tamoxifen were seeded onto the 24-well coated with Matrigel. This was done in triplicate. Cells were checked daily, and media were changed every 3 days without disturbing the Matrigel by removing and adding media gently. To stain the colonies, cells were fixed by formaldehyde (3.7% in PBS) for 2 min at room temperature, followed by permeabilization with 100% methanol (not ice-cold) for 20min at room temperature. The colonies were stained by 0.4% crystal violet.

Crystal Violet Staining

After removing media, cells were washed with Dulbecco’s phosphate buffered saline without calcium or magnesium (DPBS, Lonza) and treated with 4% formaldehyde in PBS (PFA) for 20 minutes. Once PFA was removed, cells were stained with 0.1% crystal violet stain for 1 hour. KGA/GAC overexpression colony formation: 5000 cells were plated of BE2C control, KGA and GAC overexpressing cells and were culture for 7 days; 10,000 cells were plated of SKNAS control, KGA and GAC overexpressing cells and were culture for 7 days (n=3). After 7 days, medium was removed and cells were washed with Dulbecco’s PBS (DPBS) (DPBS, Lonza) and treated with 4% formaldehyde in PBS [paraformaldehyde (PFA)] for 30 min. PFA was later removed and cells were stained with 0.1% crystal violet stain for 1 hour. Experiments were repeated twice. Indisulam treatment on WT and JMJD6-KO cell lines: JMJD6-WT and KO cells were plated in 12-well plate (50,000 cells/well for SKNAS) and 6-well plate (5,000 cells/well for BE2C) (n=3). Next day, cells were treated with Indisulam with indicated concentration for 7 days (SKNAS cells) and 5 days (BE2C cells). Crystal violet staining was performed to visualize and quantify the colony formation. Experiments were repeated twice. Indisulam treatment on KGA/GAC overexpressing cell lines: 10,000 BE2C and 100,000 SKNAS control, KGA and GAC overexpressing cells were plated in 6-well plate (n=3). Next day, cells were treated with Indisulam for 7 days. Crystal violet staining was performed to visualize and quantify the colony formation. Experiments were repeated twice.

Click-iT AHA labeling assay for metabolic labeling of newly synthesized proteins

Click-iT was performed as previously described. Briefly, cells were plated at 5 million cells per 100mm dish in RPMI supplemented with 10% FBS. Cells were washed with warm PBS and replaced with methionine-free medium (Thermo 21013024 supplement with glutamine and sodium pyruvate) for 1 hour at 37° C in 5% CO². Following, fresh methionine-free media containing 50μM of Click-iT AHA (L-azidohomoalanine) (Thermo C10102) was added to the cells for 2 hours at 37°C. After AHA-labeling, cells were washed with warm PBS and lysed with 1% SDS, 50mM Tris-HCl, (pH 8.0) supplemented with phosphatase inhibitors (PhosSTOP, Sigma) and protease inhibitors (cOmplete mini, Roche) by applying the buffer directly to the plate, incubating the cells on ice for 30 minutes, tilting the plates, and collecting the lysate. Lysates were briefly sonicated, vortexed for 5 minutes, and centrifuged at 18,000xg for 5 minutes at 4°C. Total protein quantification was assayed using the EZQ Protein Quantification Kit (Thermo R33200) according to the manufacturer’s protocol and results were read on a fluorescence-based microplate reader (BioTek Synergy 2). Click chemistry of the biotin-alkyne (PEG4 carboxamide-propargyl biotin) (Thermo B10185) to the AHA-labeled lysates was performed using the Click-iT Protein Reaction Buffer Kit (Thermo C10276) using a concentration of 40μM biotin-alkyne per click reaction (and no biotin-alkyne added for controls). Following the click reaction, samples were either assayed for total biotinylated protein by following the manufacturer’s protocol. For total biotinylated protein, briefly, 600μL of methanol, 150μL of chloroform, and 400μL of megaOhm water was sequentially added and vortexed, followed by centrifugation at 18,000xg for 5 minutes. The upper aqueous phase was discarded, and 450μL of methanol was added, vortexed, and centrifuged again at 18,000xg for 5 minutes. This methanol step was performed in duplicate to remove residual reaction components. Protein pellets were allowed to air dry and resuspended in a suitable volume of sample buffer and heated prior to western blot analysis.

Immunoprecipitation

5X10⁶ BE2C and SK-N-AS cells expressing Flag-JMJD6 were cultured in 150cm dish with complete RPMI media. Cells were washed twice with cold PBS after reaching 95% confluency, then lysed in 1ml lysis buffer (50mM Tris-HCl, pH 7.4, 150mM NaCl, 1mM EDTA, 1% Trion-X100 with complete protease inhibitors (Sigma 11836170001, added fresh) and PhosSTOP (Sigma 4906845001). Cells were scrapped into a 1.5mL Eppendorf tube and incubated on ice for 15min, which were mixed by vortex every 5 min. Cell lysates were spun by 135000rpm for 10min at 4C. The supernatant was transferred to a new tube. The cell lysates were subject to immunoprecipitation using M2 anti-Flag beads (Sigma, M8823) overnight by rocking at 4°C. The following day, beads were washed 3X with buffer and eluted with 5 packed gel volumes of FLAG peptide in TBS buffer (3uL of stock FLAG peptide at 5ug/uL per 100uL of TBS buffer) while rotating at 4°C for 30 minutes. Beads were briefly spun and the supernatant was removed from the beads (eluate). This elution step was repeated one more time and pooled with the first eluate. Prior to western blot, input, flow throughs, and elution samples were processed by adding 4X sample buffer supplemented with 50mM DTT and heated at 75°C for 10 minutes prior to running on a gel.

RNA-immunoprecipitation

SK-N-AS cells expressing Flag-JMJD6 were grown in a 10-cm dish in RPMI complete media. After 70% confluency, cells were washed with cold PBS twice and then were subject to lysis with Polysome Lysis Buffer (100mM KCl, 5mM MgCl², 10mM HEPES, pH7.0, 0.5% NP-40, 1mM DTT, 100 U/ml RAasin RNase inhibitor (Promega, N2511), 2mM vanadyl ribonucleoside complexes solution (Sigma, 94742), 25μL protease inhibitor cocktail for mammalian cells (Sigma, P8340)). Cell lysates were precleared with magnetic IgG beads for 1 hour. The cell lysates were subject to immunoprecipitation using M2 anti-Flag beads (Sigma, M8823) overnight by rocking at 4°C. The same amounts of lysates were saved at -80°C for input RNA extraction. The beads were washed with 250 μL Polysome Lysis Buffer for 4 times, flowed by washing with Polysome Lysis Buffer containing 1M urea. RNA was released by adding 150μL of Polysome Lysis Buffer containing 0.1% SDS and 45μg protease K (Ambion, AM2548) and incubated at 50°C for 30min. RNA was extracted with phenol-chloroform-isoamyl alcohol mixture (Sigma, 77618). RNA was recovered by adding 2μL of GlycoBlue (15mg/ml, Ambion, AM9516), 36μL of 3M sodium acetate and 750μL ethanol followed by incubation at -20°C for overnight. RNA was precipitated with 70% ethanol and air dried, followed by resuspension with RNase-free water followed by DNaseI (Promega, M6101) treatment to remove genomic DNA. The resultant RNAs were subjected to RT-qPCR analysis using 3 sets of GAC and KGA primers and 18S rRNA as control. 18S rRNA F: GCTTAATTTGACTCAACACGGGA; 18S rRNA R: AGCTATCAATCTGTCAATCCTGTC. GLS-GACiso_F: GAGGTGCTGGCCAAAAAGCCT; GLS-GACiso_R: AGGCATTCGGTTGCCCAAACT. GLS-KGAiso_F: CTGCAGAGGGTCATGTTGAA; GLS-KGAiso_R: ATCCATGGGAGTGTTATTCCA. KGA_set2_F: GCAGCCTCCAGGTGCTTTCA; KGA_set2_R: GTAATGGGAGGGCAGTGGCA. KGA_set3_F: TGCCCGACACTGCCCTTTAG; KGA_set3_R: CCTGCCAGACAGACAACAGCA. GAC_set2_F: TGCTTCTCAAGGCCTTACTGC; GAC_set2_R: AGGCATTCGGTTGCCCAAACT. GAC_set3_F: CCTTCTAGAGGTGCTGGCCAAA; GAC_set3_R: TGCAACACAAATATGCAGTAAGGC. For validation of protein immunoprecipitation, 20% of beads after overnight incubation were removed and processed as follows: Beads were washed 3X with buffer and eluted with 5 packed gel volumes of FLAG peptide in TBS buffer (3μL of stock FLAG peptide at 5μg/μL per 100uL of TBS buffer) while rotating at 4°C for 30 minutes. Beads were briefly spun and the supernatant was removed from the beads (eluate). This elution step was repeated one more time and pooled with the first eluate.

Identification of JMJD6 interacting partners by LC-MS/MS

Protein samples were run on a short gel as described in a previously published protocol¹⁰³. Proteins in the gel bands were reduced with dithiothreitol (DTT) (Sigma) and alkylated by iodoacetamide (IAA) (Sigma). The gel bands were then washed, dried, and rehydrated with a buffer containing trypsin (Promega). Samples were digested overnight, acidified and the resulting peptides were extracted. The extracts were dried and reconstituted in 5% formic acid. The peptide samples were loaded on a nanoscale capillary reverse phase C18 column by a HPLC system (Thermo EASY-nLC 1000) and eluted by a gradient. The eluted peptides were ionized and detected by a mass spectrometer (Thermo LTQ Orbitrap Elite). The MS and MS/MS spectra were collected over a 90-min liquid chromatography gradient. Database searches were performed using Sequest (version 28 revision 13) search engine against a composite target / decoy Uniprot human protein database. All matched MS/MS spectra were filtered by mass accuracy and matching scores to reduce protein false discovery rate to <1%. Spectral counts, matching to individual proteins reflect their relative abundance in one sample after the protein size is normalized. The spectral counts between samples for a given protein was used to calculate the p-value based on G-test ¹⁰⁴.

Metabolome profiling by LC-MS/MS

MJD6 knockout or parental SK-N-AS cells were cultured in 6-well plates to ∼85% confluence and washed with 2 mL ice cold 1X Phosphate-Buffered Saline (PBS). The cells were then harvested in 300 µL freezing 80% acetonitrile (v/v) into 1.5 mL tubes and lysed in the presence of 0.5mm Zirconia/silica beads by Bullet Blender (Next Advance) at 4 °C until the sample were homogenized. The resulting lysate was then centrifuged at 21,000 x g for 5 min and the supernatant was dried by speedvac. The samples were resuspended in 50 µL of 1% acetonitrile plus 0.1% trifluoroacetic acid, and separated by Ultra-C18 Micro spin columns (Harvard apparatus) into hydrophilic metabolites (flow through) and hydrophobic metabolites (eluent of 125 µL of 80% acetonitrile plus 0.1% trifluoroacetic acid). Ten µL of hydrophilic metabolites were dried, reconstituted in 3 µL of 66% acetonitrile and analyzed by a ZIC-HILIC column (150 × 2.1 mm, EMD Millipore) coupled with a Q Exactive HF Orbitrap MS (Thermo Fisher) in negative mode and metabolites were eluted within a 45 min gradient (buffer A: 10mM ammonium acetate in 90% acetonitrile (pH=8); buffer B: 10mM ammonium acetate in 100% H²O (pH=8)). Twenty µL of hydrophobic metabolites were dried and resuspend in 3 µL of 5% formic acid followed by separation with a self-packed nanoC18 column (75 μm × 15 cm with 1.9 µm C18 resin from Dr. Maisch GmbH) and detection with a Q Exactive HF Orbitrap MS (Thermo Fisher) in positive mode. Metabolites were eluted within a 50 min gradient (buffer A: 0.2% formic acid in H²O; buffer B: 0.2% formic acid in acetonitrile). MS settings for both types of samples included MS1 scans (120,000 resolution, 100-1000 m/z, 3 x 10⁶ AGC and 50 ms maximal ion time) and 20 data-dependent MS2 scans (30,000 resolution, 2 x 10⁵ AGC, ∼45 ms maximal ion time, HCD, Stepped NCE (50, 100, 150), and 20 s dynamic exclusion). A mix of all samples served as quality control was injected in the beginning, middle and the end of the samples to monitor the signal stability of the instrument. The data analysis was performed by a recently developed software suite JUMPm. Raw files were converted to mzXML format followed by peak feature detection for individual sample and feature alignment across samples. Metabolite identification was supported by matching the retention time, accurate mass/charge (m/z) ratio, and MS/MS fragmentation data to our in-house authentic compound library and the matching of m/z and MS/MS fragmentation data to, downloaded experimental MS/MS library (MoNA, https://mona.fiehnlab.ucdavis.edu/), in-silico database generated from Human Metabolome Database (HMDB), and mzCloud (https://mzcloud.org). Peak intensities were used for metabolite quantification. The data was normalized by both cell numbers (before data collection) and trimmed median intensity of all features across samples (post data collection).

Differential gene expression and gene set enrichment analysis (GSEA) for RNA-seq experiments

Total RNA from cells and tumor tissues were performed using the RNeasy Mini Kit (Qiagen) according to the manufacturer’s instructions. Paired-end sequencing was performed using the High-Seq platform with 100bp read length. Total stranded RNA sequencing data were processed by the internal AutoMapper pipeline. Briefly the raw reads were firs trimmed (Trim-Galore version 0.60), mapped to human genome assembly (GRCh38) (STAR v2.7) and then the gene level values were quantified (RSEM v1.31) based on GENCODE annotation (v31). Low count genes were removed from analysis using a CPM cutoff corresponding to a count of 10 reads and only confidently annotated (level 1 and 2 gene annotation) and protein-coding genes are used for differential expression analysis. Normalization factors were generated using the TMM method, counts were then transformed using voom and transformed counts were analyzed using the lmFit and eBayes functions (R limma package version 3.42.2). The significantly up- and down-regulated genes were defined by at least 2-fold changes and adjusted p-value < 0.05. Then gene set enrichment analysis (GSEA) was carried out using gene-level log2 fold changes from differential expression results against gene sets in the Molecular Signatures Database (MSigDB 6.2) (gsea2 version 2.2.3).

RNA splicing analysis

After mapping RNA-seq data, rMATS v4.1.0 was used for RNA alternative splicing analysis by using the mapped BAM files as input. Specifically, five different kinds of alternative splicing events were identified, i.e., skipped exon (SE), alternative 5’-splicing site (A5SS), alternative 3’-splicing site (A3SS), mutually exclusive exon (MXE) and intron retention (RI). To keep consistent, the same GTF annotation reference file for mapping was used for rMATS. For stranded RNA-seq data, the argument “--libType fr-firststrand” was applied. To process reads with variable lengths, the argument “--variable-read-length” was also used for rMATS. To select statistically significantly differential splicing events, the following thresholds were used: FDR <0.05 and the absolute value of IncLevelDifference > 0.1. For visualization, the IGV Genome Browser was used to show the sashimi plots of splicing events. To investigate the genome-wide correlations of differential splicing between two genotypes (e.g., shRNA knockdown of JMJD6 and non-target shRNA in cells), we extracted splice junctions for all samples of both genotypes of interest from the STAR¹⁰⁵ output files suffixed with “SJ.out.tab”, which contain high confidence collapsed splice junctions. Only those unique mapped reads crossing the junctions were considered. By extracting the union of the unique junction positions, we constructed a unified junction-read feature vector for each sample. Then, we normalized the junction-read vectors of each sample with TMM method in “voom” and “limma” and R package, assuming a negative binomial distribution. Next, we averaged the junction-read vectors for samples of the same genotype. The gene level expression was estimated based on the canonical junctions from the most abundant isoforms estimated for each gene. The fold changes of exon junctions significantly deviated from gene level changes were regarded as differentially spliced junctions for between cell-line comparisons.

Data mining

JMJD6, GAC and KGA expression in tumor tissues were downloaded from R2 (https://portals.broadinstitute.org/ccle), Kocak dataset GSE45547 (649 samples) and Fischer dataset GSE120572 (394 samples). In both datasets, the probe UKv4_A_23_P311616 represented JMJD6, the probe UKv4_A_23_P308800 represented GAC and UKv4_A_23_P39766 represented KGA. JMJD6 expression data from the RNA-seq data of various pediatric cancer tissues were downloaded from St Jude cloud (https://pecan.stjude.cloud/). The copy number alterations of JMJD6 and the related Kaplan-Meier analysis were downloaded from cBioportal (cbioportal.org). The data for correlation of metabolite abundance and JMJD6 knockout effect were downloaded from DepMap (https://depmap.org/portal/).

Pathway network analysis

The 114 essential fitness genes to neuroblastoma cell survival identified through genome-wide CRISPR/Cas9 library screen were uploaded into STRING program (https://string-db.org) for network interaction analysis with confidence threshold 0.15. The resulting network was then uploaded into Cytoscape program for presentation¹⁰⁶. The clusters were grouped based on the biological functions of each gene.

Copy number analysis of JMJD6 and other genes encoding JmjC domain histone demethylases from St Jude neuroblastoma cohort.

Somatic copy number alternations (SCNA) were determined by CONSERTING (PMID: 25938371) for each pair of tumor and normal samples. The normalize read depth ratio (log2 ratio) for the CNV segments with JmjC-domain containing proteins were extracted and used for CNV heatmap generation (https://CRAN. R-project. org/package= pheatmap) and hierarchical clustering of samples.

Xenograft studies

• (1) shRNA-mediated JMJD6 knockdown. Neuroblastoma cells were transduced with shRNA lentiviral particles targeting JMJD6. 48 hours later, 1 μg/ml of puromycin was added for selection for additional 48 hours. Cancer cells (5x10⁶) were mixed with Matrigel (1:1 ratio in volume) and subcutaneously injected into the flank sites of NSG mice. (2) JMJD6 and MYC-mediated transformation. After JoMa1 cells were transduced with GFP, JMJD6, MYCN and JMJD6/MYCN, 10⁴ cells per group were mixed with Matrigel (1:1 ratio in volume) and subcutaneously injected into the flank sites of NSG mice. Mice were sacrificed when they reached the humane endpoint. Tumors were measured by using electronic calipers, and volumes calculated as width ρχ/6 xd3 where d is the mean of two diameters taken at right angles.

Statistical analysis

All quantitative data are presented as mean ± SD. Unpaired Student’s t test was performed for comparison of two groups. Spearman correlation was used to assess the relationship between two variables. Kaplan-Meier method was used to estimate the survival rate. Mann-Whitney rank test (two-sided) was used to compare the tumor volume between two groups at every time point. P-values across multiple time points were adjusted for multiple comparison using the Holm-SidaK method. p<0.05 was considered as statistically significant. All the statistical analyses, except where otherwise noted, were performed using GraphPad Prism (v9).

Data accessibility

GEO accession number: GSE185867 To review GEO accession GSE185867: Go to https://nam11.safelinks.protection.outlook.com/?url=https%3A%2F%2Fwww.ncbi.nlm.nih.gov%2Fgeo%2Fquery%2Facc.cgi%3Facc%3DGSE185867&data=04%7C01%7Cjun.yang2%40stjude.org%7C56d8ac619531468173d608d98ea927d8%7C22340fa892264871b677d3b3e377af72%7C0%7C0%7C637697679446209542%7CUnknown%7CTWFpbGZsb3d8eyJWIjoiMC4wLjAwMDAiLCJQIjoiV2luMzIiLCJBTiI6Ik1haWwiLCJXVCI6Mn0%3D%7C1000&sdata=XtZl1ttoOPBvfhwmX9uCcwVN1Yg02qDsvcsHdJv58xw%3D&reserved=0 Enter token qzixwsoohbsxjij into the box

Acknowledgements

We thank the staff of the St. Jude Animal Resource Center and Hartwell Center for their dedication and expertise. The work was supported by American Cancer Society-Research Scholar (130421-RSG-17-071-01-TBG, J.Y.), National Cancer Institute (1R01CA229739-01, J.Y., R01CA266600, J.Y.). The work was also supported by the American Lebanese Syrian Associated Charities (ALSAC). The content is solely the responsibility of the authors and does not necessarily represent the official views of the National Institutes of Health.

Author contributions

C.J., W.Q., S.S., H.H., D.H.B. D.H performed experiments. G.W., H.J., T-C.C., D.F., J-H.C performed proteomic, genomic and RNA-seq analysis. S.M.S and S.M.P-M provided CRISPR knockout. R.W., K.F. provided key reagents. A.M.D, A.M., J.P., and J.Y. supervised the studies. J.Y conceived the project and wrote the manuscript with help from co-authors.