Abstract
Neuroblastoma is an embryonic cancer that contributes disproportionately to death in young children. Sequencing data have uncovered few recurrently mutated genes in this cancer, although epigenetic pathways have been implicated in disease pathogenesis. We used an expression-based computational screen that examined the impact of deubiquitinating enzymes on patient survival to identify potential new targets. We identified the histone H2B deubiquitinating enzyme USP44 as the enzyme with the greatest impact on survival in patients with neuroblastoma. High levels of USP44 significantly correlate with metastatic disease, unfavorable histology, advanced patient age, and MYCN amplification. The subset of patients with tumors expressing high levels of USP44 had significantly worse survival, including those with tumors lacking MYCN amplification. We showed experimentally that USP44 regulates neuroblastoma cell proliferation, migration, invasion, and neuronal development. Depletion of the histone H2B ubiquitin ligase subunit RNF20 resulted in similar findings, strongly implicating this histone mark as the target of USP44 activity in this disease. Integration of transcriptome and epigenome in analyses demonstrates a distinct set of genes that are regulated by USP44, including those in Hallmark MYC target genes in both murine embryonic fibroblasts and the SH-SY5Y neuroblastoma cell line. We conclude that USP44 is a novel epigenetic regulator that promotes aggressive features and may be a novel target in neuroblastoma.
Implications: This study identifies a new genetic marker of aggressive neuroblastoma and identifies the mechanisms by which its overactivity contributes to the pathophysiology of this disease.
Introduction
Neuroblastoma is the most common solid tumor outside the brain that affects children. The disease occurs almost exclusively in young children, with more than 90% of cases occurring within the first decade of life. Neuroblastoma most commonly develops in the adrenal medulla, the largest sympathetic ganglion in humans. Based on studies in humans and animal models, neuroblastoma is known to result from defective development of the peripheral nervous system with cells failing to undergo apoptosis and encountering a differentiation block after ganglion formation—leading to nascent tumors (1). In most cases, tumors found in the neonatal period resolve spontaneously or require comparatively low-intensity chemotherapy or surgical resection to result in cure (2). In a distinct subset, however, the disease will metastasize and/or harbor aggressive biological features such as MYCN amplification, or segmental chromosome alterations that lead to an aggressive cancer that is very difficult to cure. Therefore, there is a tremendous need for more insight into the pathogenesis of this disease in the hope of developing novel and less toxic therapies. High-throughput sequencing of exomes, or whole genomes, has uncovered relatively few recurrent single-gene variants (3–7). In addition to amplification of the MYCN oncogene in approximately 20% of tumors, which is one of the most potent negative prognostic findings in this disease, the genomic alterations most associated with aggressive disease are segmental chromosome alterations (8, 9). Alterations in portions of several chromosomes such as losses of chromosome 1p, or 11q, or gain of chromosome 17q are common and associated with poor outcomes (8, 9). Rare cases of germline predisposition to neuroblastoma have uncovered genes such as ALK and PHOX2B, which, when altered, drive changes in peripheral nervous system development and lead to tumors (10–15). A better understanding of the mediators of the early steps in neuronal development may lead to an improved understanding of tumor pathogenesis and novel therapies.
Over the last several years, several studies have implicated epigenetic pathways in the pathogenesis of neuroblastoma (16–19). Based on the observation of two distinct morphologies of the neuroblastoma cell line SK-N-SH, a super-enhancer-defined gene expression network was uncovered—defining two differentiation states that resemble adrenergic or mesenchymal cells (17, 20). Most tumors are thought to be comprised of cells in both states—though this is challenging to demonstrate using histologic techniques. Recent work has shown that the relative chemoresistance of cells in the mesenchymal state leads to the enrichment of this cell type in recurrent or refractory tumors (17, 21–23). Recent evidence indicates that cells may transit between the adrenergic state and the mesenchymal in vivo through epigenetic reprograming (24). Epigenetic regulators have also been implicated in the differentiation blockade in neuroblastoma—with histone deacetylase inhibitors leading to enhanced retinoic acid–induced differentiation (25, 26).
Deubiquitinating enzymes (DUB) are encoded by a family of nearly 100 genes in mammals that reverse selective protein ubiquitination, thereby regulating myriad processes from protein degradation to DNA damage repair. There are several examples (21–23) of DUBs acting as oncogenes (UCHL1, USP7, and others), tumor suppressors (A20 (TNFAIP3), BAP1, USP24, and others), and an extensive list of DUBs that under certain contexts may both promote or protect from cancer (27). The increasing body of evidence linking these enzymes with cancer has led to growing interest in this family as a new therapeutic target in cancer. To identify novel potential drug targets in neuroblastoma, we utilized the PREdiction of Clinical Outcomes from Genomic profiles tool (PRECOG; https://precog.stanford.edu/), which correlates gene expression with patient survival across multiple human cancers. Using this tool, we profiled the impact of 95 known or predicted DUBs on outcomes in neuroblastoma (28, 29). The results indicated that high expression of the DUB USP44 had the greatest impact on patient survival. Here, we show that USP44 is a powerful predictor of poor outcomes in neuroblastoma independently of other established risk factors and provide mechanistic data implicating it as an epigenetic regulator and potential drug target in this aggressive childhood disease.
Materials and Methods
Data analysis
The analysis of the impact of DUBs on cancer outcomes was performed using the PRECOG web tool (https://precog.stanford.edu/) as previously published (19). The datasets used in the study are listed in Supplementary Table S1. All survival analyses were initially performed using the R2: Genomics Analysis and Visualization platform (http://r2.amc.nl). Gene expression and clinical data from the datasets were downloaded from the R2 platform using the datagrabber tool and analyzed using embedded tools in GraphPad Prism or using the BlueSky statistics program. Where indicated, gene expression was stratified by the median expression, or using a percentile ranking derived from the R2 platform’s KaplanScan tool to find the cutoff with the most significant P value. The CIN25 signature score (21) was generated by determining the percentile rank expression for each of the 25 genes within the signature score using the sample Ranked GeneSetScores tool within R2. A high CIN25 score was defined as a greater than median CIN25 score. The chr. 4p score was generated by identifying the top 25 genes from chromosome 4p that were upregulated in MYCN amplified versus nonamplified tumors (see Supplementary Table S3). Tumors from the SEQC dataset were then scored using the sample Ranked GeneSetScores as described above. A high chr. 4p score was defined as a greater than median score. Unless otherwise indicated, graphs were reproduced using GraphPad Prism by exporting data from R2. Single-cell RNA-seq data were analyzed using the R2 webtool.
Cell culture, proliferation, migration, invasion, immunoblotting
Murine embryonic fibroblasts (MEF) were generated in our laboratory from day 13.5 embryos using standard techniques as reported and characterized (10). The genotype of MEFs was assessed by PCR. All MEFs were cultured as described previously (30). All experiments using primary MEFs were performed between passages 2 and 5. Mycoplasma testing was not routinely performed on MEF cultures. Immortalization of MEFs was performed by lentiviral transduction with the SV40 large T antigen as described (31). The passage number of immortalized MEFs is not monitored. Neuroblastoma cell lines were obtained from ATCC and identity verified every 3 years using the ATCC verification service. Proliferation was monitored for three independent cell cultures for each genotype/condition in replicate (minimum = 3; range 3–9) using an IncuCyte (Essen Bioscience) by monitoring five regions of each well every 2 hours. The device was maintained in a water-jacketed incubator with 5% carbon dioxide. Retinoic acid (RA)-induced differentiation was performed as described (32). Briefly, RA was added to SH-SY5Y cells, and proliferation, confluence, nuclei numbers, and neurite length were monitored using ImageJ or the IncuCyte with built-in tools where indicated. Cell migration and invasion were measured using the transwell system with the migrated/invaded cells quantitated using Cultrex cell migration and invasion kits according to the manufacturer’s instructions (R&D systems). Immunoblotting was performed using standard sodium dodecyl sulfate–polyacrylamide gel electrophoresis, semidry transfer, and the following Cell Signaling Technology antibodies: histone H2B (8135), mono-Ub-H2B (5546), RNF20 (9425), histone H3 (4499), and MYCN (9405). The USP44 antibody is from Santa Cruz (#sc-377203). All blots were developed using film and blots were scanned.
CRISPR editing and RNAi
IMR-32 cells with USP44 gene deletion were generated using the pCas-Guide system (Origene) using predesigned guide RNAs targeting USP44. Following drug selection, cells were plated using limiting dilution, and single-cell clones were selected. The deletion of the USP44 allele was verified by PCR. Although determining whether USP44 was essential was not part of the experimental design, only single-allele deletions were identified. Where indicated, the expression of USP44 was further reduced using previously characterized lentiviral shRNA (33). In experiments targeting RNF20, IMR-32 cells were stably transduced to express Cas9 using lentiCas9-Blast (Addgene, #52962). IMR-32Cas9 cells were transduced with RNF20 targeting guide RNAs (gRNA) cloned into the LentiArray vector (Invitrogen). Pooled puromycin-selected transductions were used for experiments. Lentiviral transductions of MEFs or neuroblastoma cell lines were conducted as described (11).
ChIP-seq and qPCR
Chromatin immunoprecipitation was performed by crosslinking cells with 1% formaldehyde for 20 minutes and quenched by glycine for 5 minutes. Cells were lysed with lysis buffer (150-mmol/L NaCl, 20-mmol/L EDTA, 50-mmol/L Tris-HCl, 0.5% v/v NP-40, 1% v/v Triton X-100, and 10-mmol/L NaF). Nuclear pellets were resuspended and sonicated in sonication buffer (150-mmol/L NaCl, 20-mmol/L EDTA, 50-mmol/L Tris-HCl, 1% v/v NP-40, 0.5% v/v sodium deoxycholate, 10-mmol/L NaF, and 0.1% SDS) using a Bioruptor Pico (Diagenode) and a cycle setting of 30 seconds on/off until fragment sizes around 300 to 600 bp was reached. Protease inhibitor (Roche, 11836153001) and DUB inhibitors iodoacetamide and N-ethylmaleimide were added to both lysis and sonication buffers. Samples were precleared by 50% slurry of Sepharose 4B (GE Healthcare) resuspended in sonication buffer. Antibodies were added and incubated rotating overnight at 4°C. Antibodies include 1-µg H3K4me3 (Diagenode, C15410003), 1-µg H3K27ac (Diagenode, C15410196), 2-µg H3K27me3 (Diagenode, C15410195), and 20-µL Ubiquityl-Histone H2B (Lys120; H2Bub1; 5546, Cell Signaling). Protein A-sepharose beads (Cytiva) were added to samples and incubated for 2 hours, washed, decrosslinked, and RNAse A and proteinase K treated, and DNA was extracted. Studies were performed in triplicate for each condition. Extracted DNA was sent to the Genome Analysis Core at Mayo Clinic to be library prepped and sequenced paired-end 50 bp on a HiSeq4000 (Illumina). Library preparation was carried out with ThruPLEX DNA-seq Kit V2 (Rubicon Genomics).
ChIP-seq bioinformatic analysis
Paired-end 50-bp sequencing reads were mapped to the reference genome assembly mm10 using bowtie2/2.3.3.1 (34). Triplicates were merged using samtools merge (samtools/1.9; ref. 35). Bigwig files for individual and merged bam files were generated ignoring duplicates and 1× depth (reads per genome coverage, RPGC) using bamCoverage (deeptools/3.3.2; ref. 36). Localization profiles were viewed using the Integrated Genomics Viewer (2.12.2; ref. 37). MACS2 (38) callpeak (python/3.9.2) was used to call the significant peaks with—broad cutoff 0.01 for H3K4me3 and H3K27ac using input files from respective conditions as background. ChIP occupancy was evaluated by the computeMatrix (deeptools/3.3.2) tool, and the average profiles and heatmaps were generated based on computeMatrix values with plotHeatmap (deeptools/3.3.2). mm10 transcriptional start site (TSS) coordinates were retrieved from the UCSC table browser [group: Genes and Gene Predications; track: NCBI RefSeq; table: RefSeq All (ncbiRefSeq)]. To ascertain TSS-specific regions in our system, we intersected TSS coordinates with our H3K4me3 bed file (MACS2 output) to retrieve active TSS regions using bedtools/2.27.1 (39). H3K27ac differential binding analysis was performed using DiffBind (Galaxy Version 2.10.0 + galaxy0). Default settings were used for Genomic Regions Enrichment of Annotations Tool (GREAT) and ChIP-Atlas.
RNA extraction coupled with NGS sequencing (RNA-seq)
RNA was extracted using the RNeasy Plus Mini Kit (Qiagen, 74134) and cDNA conversion using SuperScript III First-Strand Synthesis System (ThermoFisher, 18080051) per manufacturer’s protocol. Before RNA-sequencing, RNA integrity was validated by gel electrophoresis, and 1 µg of purified RNA was sent to the Genome Analysis Core at Mayo Clinic to be library prepped and sequenced stranded paired-end 100 bp on a HiSeq4000 (Illumina; MEF data) or sent to the Robert Bosch Center for Tumor Diseases Sequencing Core and sequenced stranded paired-end 59 bp on a NextSeq 2000 (SH-SY5Y data). Library preparation was carried out with TruSeq Stranded mRNA Sample Prep Kit (Illumina). Studies were performed in triplicate.
RNA-seq bioinformatic analysis
Paired-end sequencing reads were mapped to the reference genome assembly hg38 or mm10 using STAR (star/2.7.3a; ref. 40). Reverse (sense strand) reads were quantified using htseq-count (htseq/0.9.1; ref. 41) and used for differential gene expression analysis via DESeq2 (Galaxy Version 2.11.40.7 + galaxy2; ref. 42). Gene set enrichment analysis (GSEA 4.2.2; ref. 43) was performed with default settings using normalized counts from DESeq2 for expressed genes. ChIP platform Mouse_Gene_Symbol_Remapping_Human_Orthologs_MSigDB.v2022.1. Hs.chip was used to convert mouse gene symbols to humans. Z-scores were calculated with the following formula: [(samplenc − row mediannc)/row StDevnc] (nc = normalized count). Subsequent Z-scores were plotted using GraphPad Prism software.
Data availability
Data are available in a public, open-access repository: GSE262506 (SH-SY5Y RNA-seq), GSE228130 (MEF RNA-seq), GSE228131 (MEF ChIP-seq).
Statistics
Quantitative in vitro experiments involving MEFs were performed with a minimum of three independently derived MEF lines that were nonimmortalized or immortalized as indicated. Graphs of these data represent the mean ± SEM for the results obtained for each line—not the individual values. The test used is listed in each respective figure legend. Unless otherwise noted in each figure, a P value of less than 0.05 is used to denote significance. Unless otherwise noted in the figure legend, statistics were performed using the functions embedded in the software GraphPad Prism. Where indicated, calculations were performed using the statistical functions embedded in web tools including the z-scores in the PRECOG tool (https://precog.stanford.edu/) and the genomic correlations on the R2 website (http://r2.amc.nl) which were corrected for multiple comparisons using the FDR. The scan function to find the ideal cutoff values for survival based on gene expression uses P values corrected by the Bonferroni method.
Results
High levels of USP44 identify an aggressive subset of neuroblastoma tumors in children
We recently published the results of a computational analysis utilizing the PRECOG profiles tool (https://precog.stanford.edu/) in which we profiled the impact of 95 known or predicted DUBs on outcomes in neuroblastoma (28, 29). The results revealed that the enzyme USP44 has the highest positive z-score value of 8.33, which is similar to that of the canonical neuroblastoma oncogene MYCN (z-score = 8.18; Fig. 1A). To further investigate the relationship between USP44 and survival in neuroblastoma, we used public gene expression datasets that have associated clinical outcomes (Supplementary Table S1). We examined the expression of 94 DUBs in tumors classified as high-risk compared with non-high-risk tumors in the SEQC dataset (n = 498; Supplementary Table S1). USP44 was the most overexpressed DUB in high-risk tumors (Fig. 1B). Next, we stratified the cohort based on USP44 expression to examine its impact on survival using the R2: Genomics Analysis and Visualization Platform (http://r2.amc.nl). Patients whose tumors had high USP44 expression (>median) had significantly worse event-free and overall survival (EFS and OS, respectively) compared with those with lower expression (Supplementary Fig. S1A and S1B). We observed a similar adverse survival impact when we analyzed two distinct datasets [Primary NRC n = 283 (Supplementary Fig. S1C and S1D), TARGET n = 249 (Supplementary Fig. S1E and S1F)] despite substantial differences in the rates of MYCN amplification and metastatic disease in these cohorts (Supplementary Table S1). We next used the KaplanScan function of the R2 platform and identified a cutoff of 76 percentile USP44 expression that produced the most significant separation of both EFS and OS (76th and 77th percentile, respectively; Fig. 1C and D). Unless otherwise stated, USP44 normal (NL) includes tumors with transcript between 0 and 76th percentiles, whereas USP44 high (HI) includes tumors between the 77th and 100th percentiles.
We examined the relationship between USP44 expression and other established prognostic factors. There was a significant and stepwise increase in the expression of USP44 with increasing clinical stage and was highest in those tumors from patients with metastases (International Neuroblastoma Staging System stage 4; Fig. 2A). There was also a significant correlation between USP44 expression and older patient age at diagnosis, MYCN amplification, and unfavorable histology (Fig. 2B–D). We confirmed the impact of USP44 as an independent prognostic factor through univariable and multivariable analyses in which MYCN status, histology, age, and stage were included (Table 1). Despite the association of high USP44 expression with MYCN amplification, a significant impact of USP44 on survival was observed in patients with MYCN nonamplified tumors and those with non-high-risk tumors (Fig. 2E and F). These data lead us to conclude that increased USP44 expression is an independent prognostic factor for poor survival in children with neuroblastoma.
Covariate . | HR . | 95% CI . | P value . |
---|---|---|---|
Overall survival | |||
USP44 HI | 7 | 4.58–10.69 | <0.0001 |
Age at diagnosis (>18 months) | 10.04 | 6.03–16.73 | <0.0001 |
INSS stage 4 | 22.82 | 13.83–37.65 | <0.0001 |
MYCN amplified | 21.14 | 13.08–34.15 | <0.0001 |
Event-Free survival | |||
USP44 HI | 7.05 | 4.62–10.75 | <0.0001 |
Age at diagnosis (>18 months) | 10.47 | 6.28–17.46 | <0.0001 |
INSS stage 4 | 25.03 | 14.95–41.93 | <0.0001 |
MYCN amplified | 17.09 | 10.76–27.15 | <0.0001 |
Covariate . | HR . | 95% CI . | P value . |
---|---|---|---|
Overall survival | |||
USP44 HI | 7 | 4.58–10.69 | <0.0001 |
Age at diagnosis (>18 months) | 10.04 | 6.03–16.73 | <0.0001 |
INSS stage 4 | 22.82 | 13.83–37.65 | <0.0001 |
MYCN amplified | 21.14 | 13.08–34.15 | <0.0001 |
Event-Free survival | |||
USP44 HI | 7.05 | 4.62–10.75 | <0.0001 |
Age at diagnosis (>18 months) | 10.47 | 6.28–17.46 | <0.0001 |
INSS stage 4 | 25.03 | 14.95–41.93 | <0.0001 |
MYCN amplified | 17.09 | 10.76–27.15 | <0.0001 |
Data from SEQC (GSE62564) were analyzed using the Cox proportional hazard, multiple models. Cases that lacked all variables were omitted (final n = 270 patients).
Abbreviations: HR, hazard ratio; CI, confidence interval; INSS, international neuroblastoma staging system.
Because USP44 regulates H2BK120ub1 (44, 45), we speculated that it may play a role in regulating the transition between super-enhancer-driven adrenergic and mesenchymal cell states (17). However, we found no correlation between USP44 expression and the mesenchymal or adrenal scores in the neuroblastoma datasets used in our study (Supplementary Fig. S2). We previously demonstrated that USP44 is important to ensure accurate mitosis, as both USP44 loss and its overexpression lead to chromosome instability (CIN; refs. 30, 46). We have also shown that neuroblastoma tumors with a higher CIN gene expression signature (CIN25) have a worse survival (28, 47). Looking at USP44, there is a significant correlation between increasing USP44 expression and the CIN25 score (Supplementary Fig. S3A and S3B). Next, we plotted the EFS in the SEQC dataset stratified on both CIN and USP44. There was an additive effect of both parameters on survival in this analysis, with the worst outcome seen in tumors with a positive CIN25 signature and high USP44 (Supplementary Fig. S3C). The absence of either of these features resulted in an improved outcome, and patients whose tumors lacked both features did generally well. This suggests that there are likely CIN-dependent and CIN-independent effects of USP44 on neuroblastoma biology.
USP44 expression promotes features of aggressive neuroblastoma
To establish a mechanistic connection between USP44 and disease aggressiveness in neuroblastoma, we identified cell lines with low or high levels of USP44 in gene expression datasets for further experiments. The dichotomous expression of the USP44 protein was confirmed in these cells by Western blotting (Supplementary Fig. S4A and S4B). We then generated cell-line pairs by introducing USP44 into cells where it was low and depleting USP44 in cells where it was expressed (Supplementary Fig. S4C). We depleted USP44 using a two-step process. First, cells stably expressing Cas9 were transduced with a lentivirus encoding USP44 targeting gRNA and drug resistance. We isolated several single-cell clones from these transductions—but did not observe a biallelic deletion of USP44. As cells lacking one allele expressed intermediate amounts of USP44 protein, in some experiments, we additionally transduced the cells with USP44 shRNA to reduce expression further. We monitored the impact of these manipulations on cell proliferation using the IncuCyte live cell imaging system. We observed that depletion of USP44 slows proliferation in a dose-dependent fashion whereas overexpression promotes it (Fig. 3A and D). As the transcript levels are increased in metastatic tumors, we examined the impact of USP44 overexpression or depletion on in vitro surrogates of metastasis, including assays for cell migration and invasion. Again, we observed that overexpression of USP44 increased both features whereas USP44 depletion suppressed them (Fig. 3B, C, E, and F). As USP44 is both a regulator of histone ubiquitination with likely impacts on global gene expression (44, 45) and as it plays an important role in ensuring accurate mitosis (30, 46), we manipulated H2BK120ub1 independently by depleting an essential component of the ubiquitin E3 ligase complex responsible for this modification. Using stable Cas9-expressing IMR-32 cells, we introduced a control gRNA or one that targeted RNF20—an essential subunit of the RNF20/40 complex (Supplementary Fig. S4D; ref. 48). In results very similar to those seen when USP44 is overexpressed and consistent with previous work (48), we observed that RNF20 depletion enhanced the proliferation, migration, and invasion of IMR-32 cells (Fig. 3G–I). These data strongly indicate that USP44 drives malignant behavior in neuroblastoma through the suppression of H2BK120ub1.
Consistent with the reported role of USP44 and nuclear co-repressor (N-CoR) in neuronal differentiation (49), we found that USP44 expression decreased in embryonic stem cells differentiating along the neural lineage in both bulk and single-cell RNA-seq datasets (Fig. 4A and B). The drop of USP44 expression in these experiments occurred concurrently with other pluripotency markers (POU5F1 and LIN28A) and the acquisition of neural markers (NCAM1 and NTRK2/3). We also examined the expression of USP44 in single cells derived from developing mouse embryos focusing on the developing neural crest (50). The cells with the highest levels of USP44 were the premigratory neural tube cells, with expression low in most other cell types (Supplementary Fig. S5A). Interestingly, USP44 is not expressed, or is expressed at very low levels, in all cell types analyzed in a recently published analysis of the developing adrenal medulla in humans (30, 51; Supplementary Fig. S5B). A similar pattern of USP44 expression was also observed in bulk transcript data from the brain across the developmental spectrum (31; Supplementary Fig. S6A) as well as early human embryogenesis (52; Supplementary Fig. S6B). We conclude from these data that USP44 is expressed at very early stages of development within the sympathetic nervous system and declines rapidly prior to the differentiation of the adrenal medulla.
The reduction in USP44 observed during development led us to hypothesize that its continued expression in neuroblastoma tumors may contribute to the differentiation block seen in this cancer. To test this hypothesis, we introduced USP44 into SH-SY5Y cells and examined its effect on RA-induced differentiation. When exposed to RA, SH-SY5Y cells elongated and developed neuritic projections (Fig. 4C). Concurrent with this there is an increase in H2BK120ub1, as described during neural differentiation from pluripotent cells (44). We observed a significant USP44-induced reduction in RA-induced neuritic development and H2BK120ub1 (Fig. 4D and E; Supplementary Fig. S7). Taken together, these data lead us to conclude that USP44 expression declines during embryonic development of the peripheral nervous system and that aberrant continued expression, such as that observed in neuroblastoma tumors and cell lines, contributes to the differentiation block seen in this disease.
USP44 expression is correlated with copy number of chromosome 4p
Although classified by COSMIC as a tier 2 cancer gene and a likely cancer-causing gene in mice, systematic genomic sequencing projects have not demonstrated recurrent sequence variants in USP44 in neuroblastoma. Therefore, we hypothesized that USP44 expression may be driven by another genomic event. As mentioned previously, high USP44 levels were associated with MYCN-amplification. This led us to investigate whether USP44 expression was induced by MYCN. At the transcript level, there was a slight but significant negative correlation between USP44 and MYCN in tumors (Supplementary Fig. S8A and S8B). Further evidence was observed in cells with manipulated MYCN expression, in which no effect on USP44 is seen (Supplementary Fig. S8C and S8D). These data strongly indicate that the high levels of USP44 seen in MYCN-amplified tumors are not due to a direct impact of MYCN on USP44.
The long arm of chromosome 12, which includes the USP44 locus, is gained in up to 20% of tumors, with the highest proportion with gains at 12q-ter (53, 54). In contrast, copy number loss at any portion of chromosome 12 is rare (53, 54). In a distinct dataset in which copy number data were accompanied by gene expression (4), there was no correlation between copy number at the USP44 locus and its transcript level (Supplementary Fig. S9). When we broadened the analysis and looked for copy number variations that correlated with USP44, there was a significant positive correlation (FDR < 0.01) between USP44 transcript levels and copy number alterations of genes on chromosomes 4p and 14p and a significant negative correlation with genes on chromosomes 1p and 17p (Fig. 5A). When USP44 expression was correlated with conventionally defined segmental chromosomal alterations, the correlation with chr. 4p remained, with loss correlated with reduced expression (Fig. 5B), whereas there was no significant correlation between USP44 and chr. 1p status (Fig. 5C). Chromosome 4p is commonly lost in MYCN nonamplified tumors observed in MYCN-amplified tumors, whereas this region is frequently retained in MYCN amplified tumors (53, 54). In MYCN-nonamplified tumors, there was a significant reduction in USP44 in those with chr. 4p loss (Fig. 5D). Although the number of MYCN amplified tumors with chr. 4p loss is small, there remains a trend (P = 0.08) toward lower USP44 expression in MYCN amplified cases in which chr. 4p was lost (Fig. 5E). We conclude that the expression of USP44 is not a target of MYCN activity. Instead, USP44 expression is best accounted for based on the copy number of chr. 4p and that the correlation of MYCN amplification with retention of chr. 4p is the cause of the apparent correlation between USP44 and MYCN status.
As there is currently no widely available marker for cases with high USP44 expression, we next wondered whether the disomy of 4p itself may be of useful prognostic value. Using a list of the top 25 genes from chr. 4p, which are overexpressed in MYCN amplified tumors compared with nonamplified tumors, we generated a chr. 4p gene set score using the R2 embedded tool. When we applied this score across all cases within the SEQC dataset, there was a highly significant survival separation between those cases with high (>median) or low chr. 4p scores (Supplementary Fig. S10A). The negative prognostic impact of the chr. 4p score was substantially reduced in tumors where USP44 was low (<median; Supplementary Fig. S10A). As the negative prognostic impact of this score may influenced by the correlation with MYCN amplification, we examined the survival in MYCN nonamplified cases and again found those with a high chr. 4p score to have significantly worse survival (Supplementary Fig. S10B). The score is itself not an accurate test for USP44 status, as it would have a high false positive and negative rate (Supplementary Fig. S10C). These data indicate that retention of both copies of chr. 4p may be a poor prognostic biomarker for MYCN nonamplified neuroblastoma.
Histone H2B monoubiquitination correlates with gene expression patterns after loss of USP44
To directly examine the impact of USP44 on gene expression, we performed mRNA sequencing (RNA-seq) in SH-SY5Y cells overexpressing (OE) USP44. SH-SY5Y cells have low USP44 expression and high MYC expression. Heatmap analysis shows there are 303 genes (Log2FC > 0.5; base mean > 10; P-Adj < 0.05) and 476 genes (Log2FC < −0.5; base mean > 10; P-Adj < 0.05) that are upregulated and downregulated, respectively, in USP44 OE SH-SY5Y cells when compared with an empty vector (EV) control (Supplementary Fig. S11A). Interestingly, GSEA demonstrated significant enrichment of MYC target genes (Fig. 6A) in USP44 OE SH-SY5Y cells. To understand if this is a conserved function of USP44 across species, we next sought to validate this data in an orthogonal approach in MEFs. We compared wild-type (WT) MEFs to those lacking the Usp44 gene (Usp44−/−), as well as Usp44−/− MEFs into which USP44 was reintroduced (Usp44res). There were 4,694 genes that changed at least 1-fold (±0.5 Log2FC; base mean > 10; P-Adj < 0.05) between WT and Usp44−/− MEFs (2,587 up and 2,107 down in Usp44−/−). With the exception of nine genes, all genes were restored to within 0.5 log2FC of WT levels by the reintroduction of USP44 (Supplementary Fig. S11B). GSEA in Usp44−/− MEFs revealed MYC target genes are downregulated (Fig. 6B) upon loss of USP44 in MEFs, thus showing the same directionality as the USP44 OE SH-SY5Y cells in which MYC target genes are enriched. In patient neuroblastoma tumors, we also found these gene sets to be significantly enriched in USP44 HI comparted to USP44 NL tumors, highlighting the patient relevance of these findings (Fig. 6C). To further cross-validate these datasets, we performed GSEA using up- and downregulated genes in SH-SY5Y USP44 OE cells (gene set) and observed how Usp44−/− MEFs regulated these genes. We found genes downregulated in USP44 OE SH-SY5Y cells were slightly enriched (FDR = 0.082) in Usp44−/− MEFs showing similar directionality of gene regulation. Genes upregulated in SH-SY5Y USP44 OE cells were also significantly upregulated in Usp44−/− MEFs showing incorrect directionality (i.e., opposite trend of enrichment) of gene regulation (Supplementary Fig. S11C).
To understand the epigenetic impact on USP44-dependent gene regulation, we performed chromatin immunoprecipitation followed by high-throughput sequencing (ChIP-seq) for H3K4me3, H2BK120ub1, H3K27me3, and H3K27ac in WT, Usp44−/−, and Usp44res MEFs. Because we observed similar regulation of two independent MYC target gene sets, we intersected core enrichment genes from the GSEA in both SH-SY5Y USP44 OE and WT MEF (i.e., downregulated upon USP44 loss in MEFs) conditions, which resulted in 47 commonly regulated MYC target genes (Fig. 6D). Interestingly, we did not observe a change in H3K4me3 levels at these loci, however, there was a decrease in H2BK120ub1 occupancy at the TSS and along the gene bodies (GB) downstream of the TSS (Fig. 6E) suggesting ubiquitin-mediated gene regulation. This is in contrast to the tight crosstalk between H3K4me3 and H2BK120ub1 as previously described (55), although we observed a change in H3K4me3 occupancy at downregulated and upregulated genes (Supplementary Fig. S12A–S12B). Integrated Genomics Viewer tracks at the top two commonly regulated MYC target genes with the highest H2BK120ub1 occupancy, Pprc1 and Nop56, are shown in Fig. 6F.
Upon further global analysis, there was a significant alteration in H2BK120ub1 occupancy of genes that were altered in Usp44−/− cells (Supplementary Fig. S12C). Consistent with prior reports about H2BK120ub1, there was a substantial increase along GB downstream of the TSS in genes upregulated by USP44 loss. This change in H2BK120ub1 was accompanied by an increased occupancy of H3K4me3 near the TSS of these genes (Supplementary Fig. S12B). These changes were again restored by the reintroduction of USP44 into Usp44−/− cells. For unregulated genes, the loss of USP44 did not alter H2BK120ub1 or H3K4me3 in the gene body (Supplementary Fig. S12B and S12C). Interestingly, downregulated, upregulated, and unregulated genes had lower occupancy of H2BK120ub1 at the TSS in Usp44−/− MEFs compared with WT MEFs with only changes of H2BK120ub1 in the gene body correlating with gene expression change. We hypothesize that H2BK120ub1 spreading into the gene body promotes transcriptional elongation to control gene expression mediated by USP44. Lastly, given the finding that H2BK120ub1 facilitates the removal of the repressive mark H3K27me3 during differentiation (56), a histone mark deposited by the Polycomb Repressive Complex-2 that negatively regulates gene expression, we examined the effects of Usp44 expression on this mark as well. As expected, we observed an increase in genes downregulated and a decrease in genes upregulated in Usp44−/− cells, which were restored upon USP44 reintroduction (Supplementary Fig. S12D).
USP44 alters H3K27ac levels at distal elements that are bound by differentially regulated transcription factors after USP44 loss
As enhancers have been shown to drive cell states and mediate key pathways in cancer, we sought to explore differential occupancy of H3K27ac at intergenic enhancer regions as a result of USP44 loss. We performed differential H3K27ac binding analysis in WT and Usp44−/− MEFs and found 4,200 (fold change < −2; FDR < 0.01) and 1,428 (fold change > 2; FDR < 0.01) regions that were reduced or enriched in Usp44−/− MEFs, respectively. These regions were also restored upon USP44 rescue (Supplementary Fig. S13A), indicating that they are USP44-dependent. As enhancers are distal regulatory elements, we sought to determine where the differentially bound H3K27ac regions were localized using the Genomic Regions Enrichment of Annotations Tool (57). Interestingly, we observed regions that gain and lose H3K27ac localize distal to the TSS, suggesting that USP44 may influence enhancer activity to alter target gene expression (Supplementary Fig. S13B). We next utilized the ChIP-Atlas tool (58) to identify transcription factors that may localize to these regions. This tool utilizes publicly available datasets for transcription factor occupancy to identify statistically significant enrichment within specific genomic regions of interest. In regions that lost H3K27ac, we observed enrichment of transcription factors that are upregulated (Prdm16, Fos, Rnf2, Smarca2, and Fosl2) and downregulated (Fosl1, Myc, Rela, and Smarca4) in Usp44−/− MEFs. In contrast, in regions that gained H3K27ac, we observed enrichment of transcription factors that are upregulated (Vdr, Prdm16, Pparg, Cebpa, Rnf2, Smarca2, and Cebpb) and downregulated (Smarca4) in Usp44−/− MEFs (Supplementary Fig. S13C; Supplementary Table S2). Interestingly, we observed an enrichment of Myc at regions that lost H3K27ac and is downregulated in Usp44−/− MEFs, suggesting a possible feedback loop. Future studies depleting USP44 with subsequent ChIP-seq for MYC and H2BK120ub1 and HiChIP studies will help dissect the USP44-MYC enhancer–promoter interaction axis to mediate H2BK120ub1 at target genes.
Taken together, these results show that USP44 upregulation is a feature of aggressive neuroblastoma that regulates key oncogenic pathways such as MYC target genes. USP44 may regulate these pathways by influencing chromatin interactions by controlling H3K27ac levels at intergenic enhancer regions that interact with target genes to promote H2BK120ub1-mediated gene regulation. The transcription factors bound at these regions from ChIP-Atlas analyses may be interesting for future follow-up studies.
Discussion
USP44 is a DUB that was first identified in a screen designed to uncover novel regulators of the metaphase-to-anaphase transition (59). We and others have shown that this enzyme plays a complex role in genomic instability, with its activity required for accurate mitosis and DNA repair (33, 46, 59–61). These data, combined with the observed increased incidence of spontaneous and carcinogen-induced tumors in mice lacking Usp44, lead to it being considered as a tumor suppressor gene (46, 61). USP44 has also been identified as a DUB that acts on histone H2B monoubiquitinated at lysine 120 (H2BK120ub1), as part of the N-CoR complex (44, 45). H2BK120ub1 occurs predominantly in the body of genes undergoing active transcription and is catalyzed by the ubiquitin (Ub) ligase complex formed by RNF20/RNF40 (62). H2BK120ub1, USP44, and N-CoR have been implicated in differentiation (44, 49, 56, 63, 64). The role of histone ubiquitination and its regulation in the pathogenesis of neuroblastoma have not been investigated. Here, we report a novel role of USP44 as a prognostic biomarker and oncogene in neuroblastoma through its ability to regulate histone ubiquitination—and promote a MYC-like gene expression program.
Neuroblastoma is an aggressive form of childhood cancer for which novel therapies are urgently needed to improve outcomes and reduce the toxic burden of therapy. In this report, we describe a previously undescribed role of the DUB USP44 in neuroblastoma biology. Risk stratification in this disease helps to guide patient therapy, reserving more intense therapy for those with higher-risk diseases. Unfortunately, the stratification is imperfect, with many patients in the high-risk group who likely do not need intense therapy and conversely those in the low- or intermediate-risk groups who need more therapy (65). Our data show that tumors with elevated USP44 levels are at high risk for relapse and that USP44 is an independent prognostic factor in this disease (Figs. 1 and 2). This is despite the significant coexistence of established prognostic factors, including age >18 months, unfavorable histology, and MYCN amplification (Table 1). We found that USP44 mRNA and protein were correlated in immunoblotting (Supplementary Fig. S4), but there is no clear way to identify USP44 HI cases using standard pathologic techniques. However, our analysis of genomic alterations revealed that increased USP44 expression correlated with disomy at chr. 4p (Fig. 5). An expression signature score derived from genes in this region was used as a surrogate to identify tumors with both copies of chr. 4p, we found that chr. 4 disomy is a powerful predictor of poor survival in MYCN nonamplified cases (Supplementary Fig. S10). However, only 70% of USP44 HI cases can be identified with this method. Further work is therefore required to better identify cases within this new high-risk group.
Our previous work on this enzyme showed that USP44 is a nonessential enzyme, with mice lacking both copies of the gene born at Mendelian frequencies. With age, these mice develop tumors at an accelerated rate indicating that USP44 has tumor suppressor activity. Similar tumor suppressor and tumor-promoting effects have been reported for the H2B ubiquitin ligases RNF20 and RNF40 (66). The long latency, however, indicates that the short-term inhibition of USP44 that would be encountered during anticancer therapy may be well tolerated. Additionally, it is not clear whether the loss of USP44 after birth, or during the juvenile period, would have a similar tumor suppressive effect. Embryos and cells lacking USP44 have an increased rate of CIN and aneuploidy due to the mal-positioning of the centrosomes in the early phases of mitosis. This molecular phenotype is rescued by the replacement of WT USP44, but not a mutant enzyme that lacks a critical motif needed for binding to the centriole protein, centrin (46). Proteomics work has shown that this USP44-centrin complex is distinct from the complex that includes USP44 and the N-CoR subunits TBL1X and TBL1XR1, which regulates the ubiquitination status of histone H2B. This separation of function may help to explain—at least in part—how the same enzyme acts as a tumor suppressor and oncogene. Although appealing to think of these activities of mitosis regulator and epigenetic regulator as separate, in vivo, it is likely that both high and low levels of USP44 impact both processes. In support of this, we have also observed that elevated levels of USP44 lead to chromosome segregation errors and aneuploidy (34) and that overexpression of USP44 has opposite effects on the transcription of genes affected upon Usp44 knockout (Fig. 6; Supplementary Fig. S10). Our analysis of clinical data here shows that many tumors with high USP44 also have a high CIN25 signature. The additive effect of high USP44 and CIN producing worse outcomes than that seen with a high CIN25 signature alone (Supplementary Fig. S3), however, indicates that USP44 has CIN-independent effects on the biology of neuroblastoma.
In addition to its reported tumor suppressor activity, high levels of USP44 have also been observed in certain cancers. We and others have reported overexpression in T-cell leukemia (33, 67). Others have found USP44 to promote malignant glioma, as well as gastric and prostate cancers (68–71). USP44 activity toward H2BK120ub1 is required for full transcription repression mediated by N-CoR, as depletion of USP44 reduced the repressive activity of N-CoR toward known targets (45). Histone H2B monoubiquitination regulates processes including gene transcription and DNA damage repair (72). The addition of Ub to histone H2B at K120 is primarily catalyzed by the RNF20/40 E3 ligase that is recruited to chromatin through an interaction between the RNF20/40 adaptor protein WAC (WW domain-containing adapter with coiled-coil) and RNA polymerase II and coordinates transcription-coupled H2B ubiquitination (56, 62). Deregulation of H2BK120ub1 has been implicated in the regulation of early steps in organism development and in the pathogenesis of cancer. Alterations in the coding or copy numbers of the genes encoding the RNF20/40 complex, as well as H2B DUBs, have been commonly observed in several human cancers (72–74). Unlike many histone marks that cluster near the transcription start site of target genes, H2BK120ub1 is characteristically enriched downstream of the TSS in the transcribed regions of GB—a pattern that we also observed (48, 75). Increased H2BK120ub1 occupancy is associated with increased H3K4me3 and increased gene transcription (76, 77). Our RNA-seq and ChIP-seq data reinforce this, as transcripts that were upregulated in response to USP44 had increased H2BK120ub1 and H3K4me3 occupancy in Usp44 null MEFs (Supplementary Fig. S12). Both the height and width of H3K4me3 peaks were increased for upregulated genes. Increased breadth of H3K4me3 occupancy has been linked with enhanced transcriptional elongation (78) downstream of H2BK120ub1 (55, 64, 79). Genes with particularly broad H3K4me3 peaks are enriched for tumor suppressor genes and additionally may mark genes involved in establishing cell identity (78, 80). In tumors with high USP44, the loss of H2BK120ub1 and the corresponding decrease in the peak and breadth of H3K4me3 may therefore lead to reduced expression of such genes.
In back-to-back publications, the Johnsen and Oren groups also reported that, during the differentiation of stem cells, levels of H2BK120ub1 increase with time and that deregulation of H2BK120ub1 homeostasis during this process disrupts differentiation (44, 56). The association of high USP44 with unfavorable histology in neuroblastoma (Fig. 2D) is consistent with this notion, as these tumors are poorly differentiated. Our finding that USP44 overexpression suppresses H2BK120ub1 and neuritic differentiation response to RA (Fig. 4) further bolsters this conclusion. Further supporting a critical role for reduced H2BK120ub1 in those tumors with high USP44 is our data showing that RNF20 depletion also promotes proliferation, migration, and invasion (Fig. 3). As recently summarized, H2BK120ub1 protein (81) or genomic alterations in RNF20/40 (72) have been observed in several cancers suggesting that similar mechanisms may occur broadly in malignancy. USP22 is also frequently altered in human cancer with deletions and overexpression seen (72). Given the distinct gene expression patterns observed between various H2BK120ub1 DUBs, it is not surprising that they may have different impacts on cancer pathogenesis.
We provide additional compelling data that it is the suppression of histone H2B ubiquitination that drives malignant behavior, as depletion of the ubiquitin ligase component RNF20 has a very similar phenotype as that seen upon USP44 overexpression (Fig. 3). Furthermore, we observed additional factors that are enriched at differentially bound H3K27ac regions in Usp44 null MEFs (Supplementary Fig. S13C) that USP44 may directly recruit or set the epigenetic landscape for additional factors to bind. For example, MYC-associated feedback loops have been described in neuroblastoma (82–84). In Usp44 null MEFs, Myc expression is downregulated and was significantly enriched at H3K27ac regions that are lost upon USP44 loss. Thus, there may be a negative feedback loop in which downregulation of Myc leads to the loss of enhancer–promoter interactions, which normally drives a USP44-MYC-driven transcriptional program. Although USP44 OE in SH-SY5Y cells did not increase MYC expression, possibly due to the already high expression in these cells, USP44 may help promote MYC binding to enhancer regions to drive its respective transcriptional program and promote aggressiveness in neuroblastoma. Future chromatin interaction and association studies with USP44 and MYC will be needed to address these questions as well as ChIP-seq studies for USP44 to relate transcriptional changes to genome occupancy will be critical to complement our novel finding of USP44 mediating neuroblastoma aggressiveness.
We report that the histone H2B DUB USP44 drives malignant behavior and is an independent predictor of poor outcomes in childhood neuroblastoma. High levels of USP44 drive neuroblastoma cell proliferation, migration, and invasion indicating that the poor outcomes are a direct effect of USP44 activity in this disease. The driving factor(s) leading to high levels of USP44 are unclear though we have identified that tumors with two intact copies of chromosome 4p have higher expression of USP44 than those where this region is lost. USP44-driven regulation of histone H2B ubiquitination is an important mediator of these effects. Our data support the need for further investigation of the effects of USP44 inhibition as a strategy in neuroblastoma and into the safety of USP44 inhibition in animals. The data validate the need for the development of small-molecule inhibitors of this enzyme to facilitate these studies.
Authors’ Disclosures
P.J. Galardy reports grants from Hyundai Hope on Wheels, Fraternal Order of Eagles, and Mayo-Sanford Collaborative Research Fund during the conduct of the study. No disclosures were reported by the other authors.
Authors’ Contributions
T.L. Ekstrom: Investigation, writing–review and editing. S. Hussain: Data curation, investigation, writing–review and editing. T. Bedekovics: Data curation, investigation, writing–review and editing. A. Ali: Data curation, investigation, writing–review and editing. L. Paolini: Investigation, writing–review and editing. H. Mahmood: Investigation, writing–review and editing. R.M. Rosok: Resources, data curation, formal analysis, investigation. J. Koster: Conceptualization, resources, data curation, formal analysis, writing–review and editing. S.A. Johnsen: Conceptualization, data curation, supervision, methodology, writing–review and editing. P.J. Galardy: Conceptualization, resources, data curation, supervision, funding acquisition, writing–original draft, project administration, writing–review and editing.
Acknowledgments
This work was supported by grants from the Hyundai Hope on Wheels Foundation, the Mayo Clinic-Sanford Health Collaborative Research Fund, The Fraternal Order of Eagles, and the Mayo Clinic.
Note: Supplementary data for this article are available at Molecular Cancer Research Online (http://mcr.aacrjournals.org/).