RB1 loss triggers dependence on ESRRG in retinoblastoma.

Retinoblastoma (Rb) is a deadly childhood eye cancer that is classically initiated by inactivation of the RB1 tumor suppressor. Clinical management continues to rely on nonspecific chemotherapeutic agents that are associated with treatment resistance and toxicity. Here, we analyzed 103 whole exomes, 20 whole transcriptomes, 5 single-cell transcriptomes, and 4 whole genomes from primary Rb tumors to identify previously unknown Rb dependencies. Several recurrent genomic aberrations implicate estrogen-related receptor gamma (ESRRG) in Rb pathogenesis. RB1 directly interacts with and inhibits ESRRG, and RB1 loss uncouples ESRRG from negative regulation. ESRRG regulates genes involved in retinogenesis and oxygen metabolism in Rb cells. ESRRG is preferentially expressed in hypoxic Rb cells in vivo. Depletion or inhibition of ESRRG causes marked Rb cell death, which is exacerbated in hypoxia. These findings reveal a previously unidentified dependency of Rb cells on ESRRG, and they implicate ESRRG as a potential therapeutic vulnerability in Rb.

INTRODUCTION

Retinoblastoma (Rb) is the most common pediatric eye cancer and an important cause of childhood cancer death worldwide (). Despite considerable improvements in treatment and patient survival over the past century (), severe visual impairment and loss of the eye are still common, due in part to treatment resistance and toxicity associated with currently used chemotherapeutic agents. A better understanding of molecular dependencies in Rb could lead to more specific and effective targeted therapies. Rb is almost always initiated by mutational inactivation of the RB1 tumor suppressor in a susceptible retinal progenitor cell (, ), and it is generally thought that subsequent aberrations are required for full malignant transformation (, ). Several recurrent genomic events have been described in Rb, including mutations in BCOR and CREBBP, and chromosome copy number variations (CNVs) affecting 1q, 2p, 6p, 13q, and 16q (–). However, there has been limited functional validation and no clinically relevant targeted therapies associated with any proposed secondary drivers. Here, we searched for functional drivers of Rb progression that could represent potential therapeutic targets by performing a large integrative multi-omics analysis of data from whole-exome sequencing (WES), whole-genome sequencing (WGS), RNA sequencing (RNA-seq), single-cell RNA-seq (scRNA-seq), and other methods. We found that the estrogen-related receptor (ESRR) gamma (ESRRG) is an essential mediator of hypoxic adaptation and cell survival in Rb that is constitutively activated by RB1 loss and is subsequently affected by recurrent genomic aberrations in Rb. These findings suggest a selective pressure to increase ESRRG activity during Rb progression and offer a potential target for therapy.

RESULTS

Genomic landscape of Rb implicates ESRRG in tumor progression

WES data from 103 primary Rb samples were analyzed (Fig. 1, fig. S1, and data S1). Mutations were detected in RB1 in 94% of samples, including a stop-gain in 41 cases (39.8%), isodisomy of chromosome 13q in 36 cases (35.0%), loss of heterozygosity of RB1 in 19 cases (18.5%), splicing mutations in 20 cases (19.4%), homozygous deletion of RB1 in 9 cases (8.8%), and chromothripsis involving the RB1 locus in 9 cases (8.8%). Other recurrent mutations were found in BCOR (17%); FCGBP (5%); NSD1 (5%); BRCA2, CREBBP, DST, MACF1, and PI4KA (4% for each); SCN5A (3%); and CHD1, NCOA3, and PML (2% for each). The most common CNVs included 1q gain (67%), 6p gain (66%), 16q loss (56%), 2p gain (44%), 16p loss (29%), and 7q gain (26%). The MYCN locus at 2p24 was amplified in six cases, three of which had no detectable RB1 mutation. We next searched for functional patterns in the genomic landscape using Ingenuity Pathway Analysis (). Strikingly, most of the proteins encoded by recurrently mutated genes participate in an estrogen receptor/ESRR (ESR/ESRR) protein interaction network involved in development, neurogenesis, metabolism, hypoxia, and cell cycle regulation (Fig. 2A). ESRRG was the only ESR/ESRR family member that was highly expressed by RNA-seq in 16 primary Rb tumor samples (Fig. 2B). By scRNA-seq, ESRRG was the predominant ESR/ESRR family member expressed in primary Rb tumor cells (Fig. 2, C to E, and fig. S2). ESRRG is located on chromosomal region 1q41, which was the most common target of arm-level gain (Fig. 1). In addition, two of four Rb samples available for analysis by WGS and matched RNA-seq revealed more complex chromosomal rearrangements near the ESRRG locus associated with increased ESRRG mRNA expression compared to the two samples without 1q alterations (fig. S3). These findings suggest that mutations and genomic aberrations that affect ESRRG are common in Rb.

Fig. 1.

Genomic landscape of 103 primary Rbs.

Oncoprint of 103 primary Rb samples analyzed by WES. Data include status of the most commonly mutated genes, types of mutations, International Intraocular Retinoblastoma Classification (IIRC) status, family history of Rb, gender, laterality, and chromosome copy number aberrations. LOH, loss of heterozygosity.

Fig. 2.

Recurrent genomic aberrations point to ESRRG as a common target for deregulation in Rb.

(A) Ingenuity Pathway Analysis of the most common mutations in 103 Rb samples. (B) RNA-seq analysis of mRNA expression for ESR/ESRR family members in 16 primary Rb tumors. (C) Uniform manifold approximation and projection (UMAP) plots of scRNA-seq data from five primary Rb tumor samples obtained at enucleation, sorted by tumor sample (left) and Seurat cluster (right). (D) Dot plot analysis of scRNA-seq data showing relative single-cell expression of ESRRA, ESRRB, and ESRRG by tumor sample (left) and Seurat cluster (right). (E) Violin plots of scRNA-seq data showing ESRRG expression by tumor sample (left) and Seurat cluster (right). CPM, counts per million.

ESRRG regulates genes involved in retinogenesis and oxygen metabolism in Rb cells

To gain further insight into how ESRRG regulation may drive Rb progression, we performed chromatin immunoprecipitation (ChIP) using an anti-ESRRG antibody followed by DNA sequencing [ChIP sequencing (ChIP-seq)] in our RB006 and RB018 low-passage Rb cell lines (fig. S4A). A high-confidence dataset of 8163 significantly enriched peaks shared between both RB006 and RB018 was used for downstream analyses (fig. S4B and data S2). ESRRG was enriched not only at ERRE transcriptional regulatory motifs, as expected, but also at binding motifs associated with retinal and neuronal transcription factors CRX, OTX2, NEUROD, and LHX family members (fig. S4C and data S2). Many ESRRG peaks occurred in regulatory regions containing one or more of these other binding motifs but not an ERRE motif, suggesting that ESRRG can interact with chromatin independently of its canonical binding site. Using peak location analysis, ~15% of ESRRG peaks occurred in promoter regions [±3 kb from the transcription start site (TSS)], whereas almost half occurred at presumed enhancer regions between 10 and 100 kb from the TSS, with ESRRG peaks occurring at 446 (24.3%) of 1834 annotated retinal/neural super-enhancer regions (fig. S4C) ().

To correlate ESRRG chromatin localization with gene expression, we performed RNA-seq before and after depletion of ESRRG in RB006 and RB018 cells, and in two more of our recently established Rb cell lines (RB016 and wRB6). Of the 8163 significant ESRRG ChIP-seq peaks, 1226 were associated with 738 differentially expressed genes [false discovery rate (FDR) < 0.05], including 249 up-regulated genes (353 peaks) and 489 down-regulated genes (873 peaks) (data S3). Up-regulated genes were enriched for pathways related to cell cycle regulation and oxidative phosphorylation, whereas down-regulated genes were enriched for pathways related to neurogenesis, hypoxia, and estrogen signaling. CRX, OTX2, NEUROD, and LHX binding motifs were enriched in both up- and down-regulated genes, with 18.0% of peaks located in promoter regions (±3 kb from the TSS) and 52% located 10 to 100 kb from the TSS (Fig. 3, A to D, and fig. S5). Consistent with these findings, scRNA-seq analysis demonstrated a significant correlation between Rb cells expressing ESRRG and those expressing CRX and OTX2 (fig. S6).

Fig. 3.

ESRRG regulates genes involved in retinogenesis and oxygen metabolism in Rb.

(A) Integrated analysis of ESRRG ChIP-seq peaks and RNA-seq differentially expressed genes with or without shRNA-mediated knockdown of ESRRG (shESRRG). ChIP-seq was performed in RB006 and RB018 cells, and RNA-seq was performed in RB006, RB016, RB018, and wRB6 cells. The pie chart displays the percentage of peaks located within gene regions that are associated with differentially expressed genes (inner circle) or up-regulated and down-regulated genes (outer donut). The rectangular plot (top right) exhibits peak distance from the TSS. The bar plot shows the presence and co-occurrence of the most significantly enriched transcription factor binding motifs (FDR < 0.05) found within the ESRRG ChIP-seq peaks associated with up-regulated (blue) or down-regulated (red) genes after ESRRG knockdown. UTR, untranslated region. (B) ESRRG ChIP-seq track plots and corresponding RNA-seq data with (shESRRG) and without (shGFP) ESRRG knockdown at key retinal transcription factors involved in retinal development and differentiation. (C) Presence and co-occurrence of ERRE, NEUROD, LHX, and CRX/OTX2 binding motifs within ±3 kb from the TSS of all significantly up-regulated (darker colors) and down-regulated (lighter colors) genes (FDR < 0.05) following ESRRG knockdown. Individual genes are depicted as radii, and each donut represents the indicated binding motif.

RB1 interacts with and inhibits ESRRG

One of the frequently mutated genes implicated in the ESR/ESRR interaction network was RB1 itself (Fig. 2A). ESR1, ESRRA, and ESRRG contain a VXXLYD motif, similar to the IXXLFD RB1-binding motif in E2F1 and E2F2 (Fig. 4A) (–). In coimmunoprecipitation experiments, we found that endogenous RB1 and ESRRG interact and that interaction of exogenously expressed RB1 and ESRRG is disrupted by alteration of the VXXLYD motif in ESRRG (Fig. 4, B and C). Among the three ESR/ESRR family members with a VXXLYD motif, only ESRRG contains this motif within its DNA binding domain (Fig. 4A), raising the possibility that RB1 inhibits ESRRG’s transcriptional activation by blocking its interaction with ERRE regulatory motifs (). Consistent with this idea, the ectopic expression of wild-type ESRRG (V5-ESRRG-WT) activated a luciferase reporter downstream of the ERRE promoter, and this activation was enhanced by short hairpin RNA (shRNA)–mediated knockdown of RB1 in RB1-WT human embryonic kidney (HEK) 293 cells (Fig. 4D). Expression of an ESRRG mutant in which the VXXLYD motif was mutated (V5-ESRRG-MUT) resulted in greater activation of the reporter, which was only minimally enhanced by RB1 knockdown, indicating that loss of the RB1-ESRRG interaction results in constitutive ESRRG activity. Further supporting this finding, the ERRE luciferase reporter was activated by WT ESRRG in RB1-null C33A cells, and this activation was abrogated by ectopic expression of RB1. The RB1-binding mutant V5-ESRRG-MUT also activated the reporter, and it was resistant to inhibition by RB1. Moreover, we showed using ChIP–quantitative polymerase chain reaction (qPCR) that exogenous expression of RB1 in C33A cells decreased interaction of ESRRG with promoters of canonical ESRRG-regulated genes (SDHD, IDH3A, and ATP5G3) (, ), whereas knockout (KO) of RB1 in RB1-WT SH-SY5Y neuroblastoma cells resulted in increased ESRRG interaction with these promoters (Fig. 4E and fig. S7). Consistent with previous work in prostate cancer cells (), we found that ectopic expression of V5-ESRRG-WT in HEK293 cells induced p21 expression associated with hypophosphorylation (activation) of RB1, whereas knockdown of endogenous ESRRG resulted in depletion of p21 and hyperphosphorylation (inactivation) of RB1 (Fig. 4F). Similarly, in Rb cells, ESRRG was enriched at the CDKN1A gene (encoding p21), and knockdown of ESRRG decreased CDKN1A RNA expression, indicating direct transcriptional activation of CDKN1A by ESRRG (Fig. 4G). Thus, RB1 directly inhibits ESRRG, which, in turn, negatively regulates itself by activating p21, but these inhibitory mechanisms are lost in Rb. After RB1, the second most commonly mutated gene in Rb is the BCL6 corepressor (BCOR) (Fig. 1), and we found that BCOR represses ESRRG transcriptional activity independently of RB1 (fig. S8), suggesting that mutations in this gene also serve to increase ESRRG activity.

Fig. 4.

RB1 interacts with and inhibits ESRRG.

(A) Protein maps showing location of RB1-binding VXXLYD motifs in ESRs (ESR1 and ESR2) and estrogen-related receptors (ESRRA, ESRRB, and ESRRG). AF, activation function domain; a.a., amino acids. (B) Western blot (WB) for endogenously expressed ESRRG or RB1 in RB1-WT SH-SY5Y neuroblastoma cells following immunoprecipitation for RB1, ESRRG or an immunoglobulin G (IgG) control. (C) Western blot for ectopically expressed hemagglutinin (HA)–tagged RB1 (HA-RB1), V5-tagged WT ESRRG (V5-ESRRG-WT), or ESRRG with mutated VXXLYD motif (V5-ESRRG-MUT) in RB1-null C33A cells following immunoprecipitation for HA (left). Western blot for ectopically expressed HA-RB1, V5-ESRRG-WT, or V5-ESRRG-MUT in RB1-null C33A cells following immunoprecipitation for V5 (right). IP, immunoprecipitation. (D) Normalized ERRE luciferase reporter activity with or without ectopic expression of V5-ESRRG-WT or V5-ESRRG-MUT as well as with (shRB1+) or without (shRB1−) doxycycline-induced RB1 knockdown in HEK293 cells (left, n = 12) or with (RB1) or without (EV2) RB1 re-expression in C33A cells (right, n = 12), respectively. Data are shown as means ± SEM. *P ≤ 0.05, **P ≤ 0.01, and ***P ≤ 0.001 (Student’s t test). (E) Chromatin immunoprecipitation for ESRRG followed by quantitative PCR in RB1-null C33A cells with and without exogenously expressed RB1 (n = 4) and in RB1-WT SH-SY5Y neuroblastoma cells with (RB1-KO) or without (RB1-WT) RB1-KO (n = 4). (F) Western blot in HEK293 cells following ectopic expression of V5-ESRRG-WT, shRNA directed against ESRRG (shESRRG), or controls (GFP and shGFP, respectively). (G) ESRRG ChIP-seq peaks at the CDKN1A locus and corresponding RNA-seq data with (shESRRG) and without (shGFP) ESRRG knockdown. DBD, DNA binding domain; EV, empty vector; LBD, ligand binding domain.

ESRRG is a hypoxic adaptation and survival factor in Rb cells

Given the known role of ESRRG in hypoxic adaptation in the retina and other tissues (, ), we investigated how RB1 loss affects regulation of ESRRG under hypoxic conditions. In RB1-WT HEK293 cells, hypoxic conditions (1% O2) induced a transient ~3-fold increase in ESRRG mRNA expression, which was associated with transient modest activation of the ERRE luciferase reporter (Fig. 5A). In HEK293 cells depleted of RB1, ESRRG expression levels were more than threefold higher at baseline and demonstrated an exaggerated increase in response to hypoxia, which was accompanied by hyperactivity of the ERRE reporter over the same time course. Conversely, RB1-null C33A cells demonstrated an exaggerated response to hypoxia in terms of both ESRRG expression and ERRE reporter activity, both of which were repressed by ectopic expression of RB1 (Fig. 5A). In Rb cell lines, shRNA-mediated depletion of ESRRG or inhibition of ESRRG using the specific inverse agonist GSK5182 () resulted in marked cell death, which was exacerbated in hypoxia (Fig. 5, B to D). Together, these findings suggest that RB1 dampens the hypoxia-induced surge in ESRRG activity and that loss of RB1 abolishes this homeostatic mechanism, resulting in nascent Rb cells becoming increasingly dependent on ESRRG in the hypoxic tumor microenvironment. Consistent with this possibility, ESRRG is normally expressed only in isolated retinal cells (fig. S9), but it is diffusely expressed in Rb tumors, particularly in hypoxic regions beyond ~100 μm of tumor blood vessels (Fig. 6A) and in hypoxic vitreous “seeds” (Fig. 6B). Intense ESRRG expression was also observed in Rb cells invading the optic nerve (Fig. 6C), which is strongly associated with metastasis and poor outcome ().

Fig. 5.

ESRRG is a hypoxic adaptation and survival factor in Rb.

(A) Normalized gene expression of ESRRG (top left, n = 3) and normalized ERRE luciferase activity (top right, n = 6) in HEK293 cells with (shRB1+) or without (shRB1−) doxycycline treatment to induce shRNA-mediated knockdown of RB1 at the indicated time points in hypoxia (1% O2). Normalized gene expression of ESRRG (bottom left, n = 4) and normalized ERRE luciferase activity (bottom right, n = 3) in C33A cells with ectopically expressed RB1 (RB1) or EV control at the indicated time points in hypoxia (1% O2). (B) Representative images of LIVE (green)/DEAD (red) cell viability assays of three recently established Rb cell lines (RB006, RB018, and RB020) stably expressing shRNA directed against ESRRG (shESRRG) or GFP control (shGFP) cultured in normoxia versus hypoxia (1% O2) for 21 days. (C) Normalized LIVE/DEAD cell ratios averaged for all three Rb cell lines illustrated in (B) in normoxia versus hypoxia O2 (1%). *P ≤ 0.05, **P ≤ 0.01, and ***P ≤ 0.001 (Student’s t test). (D) Normalized LIVE/DEAD cell ratios in wRB6 cells (left, n = 3) and RB020 cells (right, n = 3) treated with the indicated micromolar concentrations of the ESRRG inverse agonist GSK5182 for 7 days in normoxia versus hypoxia (1% O2). *P ≤ 0.05, **P ≤ 0.01, and ***P ≤ 0.001 (Student’s t test). All data are shown as means ± SEM.

Fig. 6.

ESRRG expression in human Rbs.

(A) Hematoxylin and eosin staining (left) and multiplex fluorescence immunohistochemistry (middle and right) of Rb sample #2587-18. ESRRG (white), TUBB3 (red), 4′,6-diamidino-2-phenylindole (DAPI) (blue), and HIF1A (green). V, viable tumor cells concentrated in cuffs around blood vessels; N, necrotic tumor cells beyond ~100 μm of tumor blood vessels. Arrows demarcate regions enriched for ESRRG- and HIF1A-expressing cells in hypoxic zones between viable and necrotic tumor cells. (B and C) Similar multiplex fluorescence immunohistochemistry findings in Rb samples #2647-18 (left panels) and #1109189582 (right panels). ESRRG (white), TUBB3 (red), DAPI (blue), and HIF1A (green). (D) Vitreous tumor cell cluster or “seed” in Rb sample #2647-18 demonstrating strong diffuse staining for ESRRG and HIF1A. ESRRG (white), TUBB3 (red), DAPI (blue), and HIF1A (green). (E) Hematoxylin and eosin staining (left) and multiplex fluorescence immunohistochemistry (right) in Rb sample #29-15 showing invasion of the optic nerve by tumor cells (arrows). ESRRG (white), TUBB3 (red), DAPI (blue), and PECAM1 (green).

DISCUSSION

In this study, we show that (i) recurrent mutations in Rb tumors may have in common the effect of increasing ESRRG activity, (ii) RB1 interacts with and inhibits ESRRG such that mutational inactivation of RB1 may release ESRRG from negative regulation and confer upon it oncogenic potential, and (iii) Rb cells are dependent on ESRRG for survival and hypoxic adaptation, which may create a selective pressure for acquiring genomic aberrations that increase ESRRG activity. The most commonly mutated gene in Rb, after RB1 itself, is the transcriptional corepressor BCOR, a component of the noncanonical polycomb repressive complex (PRC) 1.1 that is also mutated in other developmental disorders affecting the eye (, ). Inactivating mutations in BCOR are found in about 20% of Rb tumors and are associated with poor prognosis (, ). We show here that BCOR represses ERRE-mediated transcription independently of RB1, suggesting that BCOR loss derepresses ESRRG beyond that associated with RB1 loss. Because ESRRG is located on chromosome 1q, the most common chromosome copy number aberration in Rb, and as we found examples of more complex genomic rearrangements affecting the ESRRG locus, it is interesting to speculate that chromosome-level alterations may also contribute to increased ESRRG activity in Rb. Several other genes located on 1q have also been implicated in Rb progression (, ), suggesting that 1q gain may promote tumor progression through its effect on multiple genes. Together, recurrent mutations in Rb may have in common the effect of increasing ESRRG activity, which could make ESRRG an attractive therapeutic target.

ESRRG is an estrogen-related orphan nuclear receptor transcription factor that is normally expressed in the developing retina and central nervous system, where it regulates genes involved in development, proliferation, and oxygen metabolism (, , –). In Rb cells, we found that ESRRG localizes to and regulates the transcription of numerous genes containing CRX, OTX2, NEUROD, and LHX binding motifs that regulate normal retinal development (–). CRX and OTX2 in particular are key regulators of oxidative metabolism and mitochondrial function during early neurogenesis and retinogenesis (–), and they can become oncogenic when aberrantly activated (, ). The developing retina is a hypoxic environment owing to the proliferation of retinal precursor cells, oxygen demand of newly formed neurons, and limited blood supply (, ). This oxygen demand is markedly increased in Rb due to uncontrolled proliferation of tumor cells (, ). We showed that RB1 loss increases promoter occupation of ESRRG at genes involved in oxidative metabolism and that Rb cells require ESRRG for survival, proliferation, and hypoxic adaptation. Rb cells demonstrate constitutive ESRRG expression that is enriched in oxygen-poor regions beyond ~100 μm from tumor blood vessels () and in vitreous seeds, which consist of clusters of Rb cells that proliferate in the hypoxic vitreous cavity and portend chemoresistance and poor outcome in Rb (, ). Further, ESRRG expression clearly demarcates aggressive Rb cells invading the optic nerve, which is a major risk factor for metastatic death ().

In conclusion, Rb tumors undergo ongoing hypoxic adaptation to allow continued proliferation and invasion in an increasingly oxygen-depleted environment (, ), and our findings implicate ESRRG and its decoupling from normal homeostatic regulation by RB1 and BCOR as a key mediator of this hypoxic adaptation. Further studies are warranted to investigate the potential role for inhibiting ESRRG in the management of Rb.

MATERIALS AND METHODS

Specimens and clinical data

Human Rb samples were collected from enucleated eyes at the time of surgery, snap-frozen, and stored at −80°C until used for analysis with the approval from the University of Miami Institutional Review Board (no. 20130636). Written informed consent was obtained from each patient. Clinical and histopathologic information were obtained, and samples were deidentified for further analysis. All data were included in the analysis.

DNA sequencing, quality control, and alignment

WES was performed on 31 of our Rb enucleated tumor samples, 20 of which had matched blood DNA available for sequencing. DNA was extracted using the Wizard Genomic DNA Purification Kit (Promega, Madison, WI) and the QuickGene DNA Whole Blood Kit S (Fujifilm, Tokyo, Japan), respectively. Exome fragments were captured using the SureSelect Human Exon V5 (Agilent) and sequenced on the Illumina HiSeq 2500. WES data from an additional 72 Rb samples with matched blood were obtained from the European Genome-Phenome Archive (EGAS00001001690) with permission from the Data Access Committee (EGAC00001000431) (). WGS from four Rb samples with matched blood was obtained with permission from dbGaP (phs000352.v1.p1). FASTQ files from all sources underwent the following bioinformatics pipeline. Sequence data were quality-controlled using FastQC (v0.11.3). WES and WGS reads were trimmed (if required) and aligned to the human genome (hg19/GRCH37) using NovoAlign (v3.06.04), marked for duplicates using Picard (v1.128), and realigned around small and large indels using Abra (v0.94c) (), and read mate was fixed and analyzed for coverage statistics using Picard.

Mutation analysis

WES datasets underwent variant calling for single-nucleotide polymorphisms (SNPs) and indels using MuTect2 (GATK 2016-01-25 nightly build) (). For MuTect2, mutation calls for a panel of normals (n = 117 germline blood samples) were generated and pooled together to further filter out mutations that were present in at least two normals. For tumors without matched blood samples, MuTect2 was used for variant calling with a panel of normals and a high coverage blood sample to filter out likely germline mutations. Mutations were further filtered out if alternate tumor read counts were <3, if alternate tumor read counts were <8 and blood alternate read counts were >1, and if the minor allele frequency (MAF) was <5%. For all sequencing samples, the BAM files were investigated manually for the regions of interest (RB1 and BCOR) using the Integrative Genomics Viewer (v2.3.80, Broad Institute, Cambridge, MA).

For all called mutations, ANNOVAR was used for annotation (). Following annotation, mutations were filtered out if the MAF was 1% or greater in the population according to 1000 Genomes Project (August 2015) or Exome Sequencing Project (March 2015), and mutations listed in Database of Single Nucleotide Polymorphisms (dbSNP) (v138) were filtered out, with the exception of SNPs with a MAF < 1% (or unknown) in the population, a single mapping to reference assembly, or with a “clinically associated” tag. If the exact same insertion or deletion or non-frameshift multinucleotide substitution was present in at least two samples without matched germ line and were not present in any matched tumor sample, it was removed as a likely germline variant. Functional consequences of SNPs were assessed by three mutational predictor tools: PolyPhen-2 (probably damaging, possibly damaging, benign) (), Functional Analysis through Hidden Markov Models (FATHMM) (damaging, tolerated) (), and MetaLR (damaging, tolerated) (). In nonexonic regions and for synonymous mutations, SNPs were considered deleterious if two of three of the above prediction algorithms predicted a damaging mutation. For nonsynonymous exonic SNPs, mutations were not considered deleterious if two of three algorithms predicted a benign or tolerated mutation. Probably damaging, possibly damaging, or damaging calls were considered deleterious. All insertions and deletions in exonic regions and alterations in splicing junctions were considered deleterious. For identification of additional driver mutations, the unmatched tumor samples were first excluded before being evaluated if mutations were found in matched samples and were also excluded if known to be one of the benign mutations commonly seen in cancer exomes (, ). Lollipop plots displaying the distribution of driver mutations along the protein domains of each gene were plotted using trackViewer () in R (version 4.0.3).

Copy number variations

CNVs were ascertained from WES data using CNVkit (v0.8.5) () and by cgpBattenberg using default settings (). Isodisomy was determined by plotting B allele frequencies using CNVkit, Allele-Specific Copy Number Analysis of Tumors (ASCAT), and cgpBattenberg. CNVkit output was analyzed using Genomic Identification of Significant Targets in Cancer (GISTIC) () to determine significant CNVs using cutoffs of FDR < 0.05, z score ≥ 6, and chromosomal arm score > |0.1|. MYCN amplification was determined using GISTIC with a z score cutoff of 3. Structural variants from the WGS samples were determined using DELLY (version 0.7.6) ().

Pathway analysis

Genes with mutations in at least two Rb samples were loaded into the “Path Designer” tool of the Ingenuity Pathway Analysis product (QIAGEN, Valencia, CA, USA), and an unsupervised interaction network identifying additional interacting members was generated using the “path explorer” option in the “build” workbox, restricting the software’s analysis output to paths of documented interaction with a maximum of one intermediary node between the selected proteins.

Cell culture

Newly established, patient-derived low-passage Rb cell lines (wRB6, RB006, RB010, RB015, RB016, RB018, RB020, and RB021) were derived from enucleated eyes of patients of the senior author (J.W.H.) and cultured at 5% O2 in Dulbecco’s modified essential medium (DMEM)/F12 with B-27 minus vitamin A (Life Technologies), 1% penicillin-streptomycin, 2 mM GlutaMAX (Gibco), basic fibroblast growth factor (10 ng/ml; PeproTech, catalog no. 100-18B), recombinant human stem cell factor (10 ng/ml; PeproTech, catalog no. 300-07), and epidermal growth factor (20 ng/ml; PeproTech, catalog no. AF-100-15). All low-passage Rb cell lines had an average doubling time of ~72 hours and a viability of ~80%. Human C33A cells were obtained from the American Type Culture Collection (ATCC) and maintained in Eagle’s minimum essential medium supplemented with 10% heat-inactivated fetal bovine serum (HI-FBS) and 1% penicillin-streptomycin. Human HEK293 cells were obtained from ATCC and maintained in DMEM with 10% HI-FBS, 2 mM GlutaMAX, and 1% penicillin-streptomycin. Of note, although RB1-WT HEK293 cells express adenoviral E1A that partially inhibits RB1, it does not affect the interaction of RB1 with E2F1 () or with ESRRG (fig. S10). RB1-WT neuroblastoma cell line (SH-SY5Y) was generously provided by S. Dhar’s laboratory at the University of Miami. SH-SY5Y cells were cultured in DMEM/F12 (1:1) supplemented with 10% HI-FBS. Cell proliferation and viability were measured using Trypan Blue Stain (Gibco) and/or LIVE/DEAD Cell Viability Assay (Thermo Fisher Scientific). All cells were maintained in ambient O2 unless otherwise stated. Rb cells were treated with various concentrations of the ESRRG selective inverse agonist GSK5182 (Aobious Inc., Gloucester, MA) for 7 days.

Plasmids and lentivirus expression vectors

shRNA plasmids were created by integrating validated small interfering RNA (siRNA) sequences 5′-CCTGTCAGGAAACTGTATGAT-3′ () and 5′-AATGGCCATCAGAACGGACTT-3′ targeting human ESRRG (shESRRG) into pLKO.1-Puro vector (Addgene, #8453). Nonspecific siRNA sequence 5′-AACAGCCACAACGTCTATATC-3′ or siRNA sequence 5′-CAACAGCCACAACGTCTATAT-3′ targeting green fluorescent protein (GFP) was used as controls for shRNA. Tetracycline-inducible TET-shRB1 and constitutive shRB1 constructs were created by integrating validated siRNA sequences 5′-GAAAGGACATGTGAACTTA-3′ () and 5′-GAACGATTATCCATTCAAA-3′ () targeting human RB1 (shRB1) into TET-pLKO.1-Puro vector (Addgene #21915) and pLKO.1-Puro vector, respectively. The plasmids were packaged into lentiviral particles by transient cotransfection into HEK293T cells with pMD2G and psPAX2 packaging plasmids using jetPRIME reagent (Polyplus). The pCMV–ESRRG-WT-V5 vector was created by PCR amplification of ESRRG full-length cDNA fragment (Horizon, #MHS6278-202806209) and subsequent recombination into a pCMV plasmid encoding C-terminal V5. An RB1-binding site mutant of ESRRG (V5-ESRRG-MUT) was generated using site-directed mutagenesis to replace 5′-GTCAGGAAACTGTATGATGAC-3′ with 5′-GACAAAGTCGATGCGAGGTAT-3′. A construct encoding RB1 was created by PCR-amplifying cDNA for human RB1 (Horizon, #MHS6278-202758572) and cloning into an N-terminal hemagglutinin (HA)–tagged expression vector. SH-SY5Y cells with RB1-KO were created by transducing viruses encoding the CRISPR-Cas9 with guide RNAs targeting RB1 (5′-CCGAAAAACGGCCGCCACCGC-3′ and 5′-CCGCGGAACCCCCGGCACCGC-3′). The transduced cells were selected with hygromycin and puromycin and clonally selected for RB1-KO. The pLL-FLAG-BCOR plasmid was a gift from V. Bardwell. The ERRE promoter luciferase reporter plasmid 3xERRE-luciferase (pERRE-LUC; Addgene, # 37851) was a gift from R. Riggins. CMV promoter–driven beta-galactosidase plasmid (pCMV-beta-gal) was obtained from ATCC (#77177). Lentiviral particles were concentrated by precipitation from culture supernatants using PEG-it (System Biosciences) and quantified with the QuickTiter Lentivirus Titer Kit (Cell Bio Labs). Gene knockdown and overexpression were verified with Western blot or qPCR using appropriate controls (Fig. 4 and figs. S10 and S11).

Luciferase assays

RB1-null C33A cells were cotransfected with indicated quantities of pCMV-HA-RB1, pLL-FLAG-BCOR, pCMV-V5-ESRRG-WT or pCMV-V5-ESRRG-MUT, and pERRE-LUC using jetPRIME reagent (Polyplus). Stoichiometric quantities of empty vectors pCMV-V5-EV, pCMV-HA-EV, and pLL-FLAG-EV were used for C33A cell transfection normalization. HEK293 cells expressing TET-shRB1 (HEK293-TET-shRB1) were induced with doxycycline for 48 hours and then cotransfected with indicated quantities of pCMV-HA-RB1, pCMV-V5-ESRRG-WT or pCMV-V5-ESRRG-MUT, and pERRE-LUC. Cells not treated with doxycycline and stoichiometric quantities of EVs pCMV-V5-EV and pLL-FLAG-EV were used for HEK293 cell transfection normalization. Luciferase activity was measured 24 hours after plasmid transduction by lysing the cells in the passive buffer (Promega) and combining equal volumes of cell lysates with the Bright-Glo reagent (Promega). Cotransfection with a CMV-beta-gal plasmid and bicinchoninic acid protein quantification assays were used for luciferase signal normalization.

Coimmunoprecipitations

RB1-null C33A and RB1-WT HEK293 cells (107 cells per condition) and SH-SY5Y neuroblastoma cells (207 cells per condition) were cultured in 15-cm dishes, transfected with the indicated expression plasmids, and lysed with ice-cold NP-40 buffer supplemented with protease and phosphatase inhibitors (Roche). All cells were treated with dithiobis succinimidyl propionate (Thermo Fisher Scientific) directly on plastic. Cell lysates were passed 20 times through a 25-gauge needle, incubated on ice for 1 hour, pelleted by centrifugation at 12,000 relative centrifugal force (RCF) for 10 min, precleared by 30 min coincubation with magnetic beads (Life Technologies) at 4°C, and incubated overnight at 4°C with 5 μg of the indicated antibodies per 1 mg of quantified protein lysates. Antibody-protein conjugates were captured by incubation with protein G Dynabeads for 1 hour, washed with ice-cold coimmunoprecipitation buffer, and eluted by boiling in SDS running buffer. Total precleared lysates were used as positive control. Antibodies used for these experiments are listed in data S4.

Western blot

Cell pellets were resuspended in lysis buffer, which consisted of 50 mM Hepes (pH 7.2), 400 mM NaCl, 0.1% NP-40, 0.5 mM EDTA (pH 8), 2.5 mM dithiothreitol plus protease, and phosphatase inhibitors. Resulting cell lysates were incubated on ice for 10 min, sonicated at medium power for 10 s, centrifuged for 10 min at 10,000 RCF to remove cellular debris, and quantified, and equal quantities were denatured by boiling with SDS loading buffer for 5 min at 90°C. Proteins (30 μg per condition) were resolved on SDS–polyacrylamide gel electrophoresis gradient gels and transferred to a polyvinylidene difluoride membrane (Thermo Fisher Scientific). Blots were blocked in 5% bovine serum albumin (BSA) in tris-buffered saline (TBS; pH 7.6), probed overnight with indicated primary antibodies in 5% BSA and 0.15% Tween 20 in TBS (TBS-T), washed in TBS-T, and incubated for 1 hour with a secondary antibody in 5% BSA/TBS-T solution. Proteins were visualized using SuperSignal chemiluminescence (Thermo Fisher Scientific) or LI-COR imaging technology (LI-COR Biosciences). Antibodies used for these experiments are listed in data S4.

Hypoxia experiments

Recently established Rb cell lines (RB006, RB018, RB020, and RB021) were transduced with either shESRRG or shCTRL lentiviral particles and treated with puromycin (2 μg/ml) 48 hours after transduction to select for stable lentiviral integration. At 7 days after transduction, cells were maintained in normoxia or hypoxia (1% O2). Cell viability at day 21 was measured with the LIVE/DEAD Cell Viability Assay (Thermo Fisher Scientific). Whole wells were scanned using the EVOS FL Auto Imaging System (Thermo Fisher Scientific), and image fluorescence was quantified using particle analysis with ImageJ software. HEK293 cells stably integrated with TET-shRB1 were induced with doxycycline hyclate (1 μg/ml) for 48 hours, transfected with pERRE-LUC plasmid, and maintained in either normoxia or hypoxia (1% O2) for an additional 24 hours. C33A cells were transfected with pCMV-HA-RB1 and pERRE-LUC plasmids for 24 hours and exposed to either normoxia or hypoxia (1% O2) for 24 hours. Uninduced HEK293-TET-shRB1 and C33A cells transfected with pCMV-HA-EV were used as experimental controls.

Quantitative real-time PCR

Total RNA was extracted from cells with a combination of TRIzol and an RNeasy Mini RNA isolation kit (QIAGEN, Valencia, CA, USA) and reverse-transcribed using a high-capacity cDNA transcription kit (Applied Biosystems, Foster City, CA, USA). qPCR was performed in triplicate using the 7300 Real-time reverse transcription polymerase chain reaction (RT-PCR) system (Applied Biosystems) according to the manufacturer’s description using the following thermocycler parameters: 1 min of initial activation at 98°C, 40 cycles of 5-s denaturation at 98°C, and 55-s annealing and extension at 60°C. The relative gene expression data were analyzed by the comparative CT (ΔΔCT) method. The results were normalized to GAPDH, RNA18S, and ACTB as internal controls. Oligonucleotide primer sequences used are listed in data S4.

ChIP sequencing and DNA sequencing

Recently established low-passage RB006 and RB018 cell lines with 1q gain, C33A cells with and without exogenous RB1 expression, and SH-SY5Y neuroblastoma cells with and without RB1-KO were grown to approximately 5 × 107 cells, collected, and pelleted. Cell pellets were fixed in 1% formaldehyde for 10 min to cross-link protein-DNA interactions. Cell nuclei were extracted and subjected to sonication in a Bioruptor Pico (Diagenode). Shearing was optimized to yield chromatin fragments of 200 to 500 base pairs. Chromatin was incubated in ChIP-qualified antibody against ESRRG (data S4). Antibodies were captured using protein G–coupled magnetic beads, and protein-DNA was eluted. DNA was purified for downstream analysis. Successful ChIP was validated with RT-PCR for known binding genes and known negative controls. Technical controls included nonspecific immunoglobulin G incubation of cell lysates and no-antibody (input) controls for nonspecific DNA binding to protein G–coupled beads. DNA was sent for library preparation and sequencing in the Oncogenomics Shared Resource of the University of Miami. Bioinformatic analysis was performed by established ChIP-seq pipelines including quality control (FastQC), adapter trimming and alignment (NovoAlign), and normalization to an input control followed by peak calling (MACS2). Heatmap and band plots were generated using deepTools. Overlapping high-confidence peaks that were called in both the RB006 and RB018 ESRRG pulldown ChIP-seq datasets as well as a pooled dataset that combined the RB006 and RB018 data were determined using HOMER (P < 0.001) and used for downstream analyses. Euler diagram of peaks called for each dataset and overlapping datasets was plotted using eulerAPE (version 3.0) (). Motif and gene set enrichment analyses (GSEAs) were conducted using HOMER. Super-enhancer regions identified in retinal and neural tissue were downloaded from the comprehensive human super-enhancer database (version 1.03) (). Peaks with pooled peak calls, super-enhancer locations, exon and intron locations, and direction of transcription were plotted for visualization using our custom program SparK (version 2.6.2) (). Peak location, motif analysis, and co-occurrence visualization were conducted using ChIPseeker, UpSetR, and circlize in R (–). ChIP-qPCR primer sequences are listed in data S4.

RNA-seq and analysis

RNA was extracted from tumor samples or Rb cell lines (RB006, RB016, RB018, and wRB6) with or without shESRRG knockdown using a combination of TRIzol and an RNeasy Mini RNA isolation kit (QIAGEN, Valencia, CA, USA). RNA-seq libraries were prepared using the TruSeq Stranded Total RNA Prep Kit with Ribo-Zero Gold to remove cytoplasmic and mitochondrial ribosomal RNA according to the manufacturer’s recommendation (Illumina, San Diego, CA). Total RNA-seq libraries were run on an Illumina NextSeq 500 sequencing instrument according to the manufacturer’s protocol. RNA-seq from the four Rb samples that underwent WGS was obtained with permission from dbGaP (phs000352.v1.p1). Reads were aligned and gene counts were made using STAR, data quality was assessed using FastQC and RSeQC, and gene expression was normalized, batch-corrected, and determined using EdgeR. Significant differences (FDR < 0.05) between cells expressing shESRRG versus shGFP control were calculated using EdgeR. For tumor samples that underwent RNA-seq without a detectable biallelic loss of RB1 seen in exome analysis, RNA tracks were manually curated to look for any RB1 mutations.

Integration of ChIP-seq and RNA-seq data

Significant ESRRG ChIP-seq peaks overlapping among RB006, RB018, and pooled samples were correlated with significantly up-regulated or down-regulated genes after shESRRG knockdown in four primary cultured cell lines (RB006, RB016, RB018, and wRB6). Genes that were significantly up-regulated or down-regulated in association with ESRRG ChIP-seq peaks were analyzed with GSEA and MSigDB. Differentially expressed genes after ESRRG knockdown were evaluated for ESRRG-associated transcription factor motifs.

scRNA-seq analysis

Single-cell suspensions were counted using the Cellometer K2 Fluorescent Viability Cell Counter (Nexcelom), verified manually using a hemocytometer, and adjusted to 1000 cells/μl. For RNA-seq, samples were prepared using the Chromium Single Cell 5′ Library & Gel Bead Kit v2 (10X Genomics) according to the manufacturer’s protocol with a capture target of 10,000 cells. Each sample was processed on an independent Chromium Single Cell A Chip (10X Genomics) and run on a thermocycler (Eppendorf). 5′ gene expression libraries were sequenced using NextSeq 500 150-cycle high-output flow cells or NovaSeq 6000 200-cycle S4 flow cells. Raw base call (BCL) files were analyzed using CellRanger (version 3.0.2). The “mkfastq” command was used to generate FASTQ files, and the “count” command was used to generate raw gene barcode matrices aligned to the GRCh38 Ensembl build 93 genome. The data from all five samples were combined in R (3.5.2) using the Seurat package (3.0.0), and an aggregate Seurat object was generated (, ). To ensure our analysis was on high-quality cells, filtering was conducted by retaining cells that had unique molecular identifiers greater than 400, expressed 100 to 6000 genes inclusive, and had mitochondrial content less than 10%. This resulted in 6371 cells for RB025, 12,808 cells for RB026, 10,675 cells for RB027, 1354 cells for RB028, and 1445 cells for RB029. Doublets were assessed using the DoubletFinder (2.0.3) algorithm () using the multiplet rates in the 10X Genomics Chromium Single Cell 3′ Reagent Kits User Guide (v3 Chemistry). Doublets were assessed by individual samples based on the number of captured cells. Data for all five samples were combined using the Standard Integration Workflow (https://satijalab.org/seurat/v3.0/integration.html). Data from each sample were normalized using the NormalizeData() function, and variable features were identified using FindVariableFeatures() with 5000 genes and the selection method set to “vst,” a variance stabilizing transformation. To identify integration anchor genes among the five samples, the FindIntegrationAnchors() function was used with 30 principal components (PCs) and 5000 genes. Using Seurat’s IntegrateData(), the samples were combined into one object. The data were scaled using the ScaleData() function to reduce the dimensionality of this dataset. Principal components analysis (PCA) was performed, and the first 30 PCs were summarized using uniform manifold approximation and projection (UMAP) dimensionality reduction (). Thirty PCs were chosen on the basis of results from analysis using JackStraw() and elbow plots. The DimPlot() function was used to generate the UMAP plots (Fig. 2C). Clustering was conducted using the FindNeighbors() and FindClusters() functions using 30 PCA components and a resolution parameter set to 1.8. The original Louvain algorithm was used for modularity optimization (). The resulting 36 louvain clusters were visualized in a two-dimensional UMAP representation. Differentially expressed genes were identified for each cluster using FindAllMarkers() and filtered for |log fold change| greater than 0.5 and adjusted P value less than 0.05. Dot plots were generated using the DotPlot() function, violin plots were generated with the vlnplot() function, and scatterplots were generated with FeatureScatter() in Seurat with “assay” set to “RNA.”

Immunohistochemistry

For immunohistochemistry, 4-μm paraffin sections of Rb tumor specimens were deparaffinized in xylene, rehydrated in an ethanol gradient, permeabilized for 30 min in 0.5% Triton X-100 (v/v in phosphate-buffered saline), and subjected to antigen retrieval using sodium citrate buffer (pH 6.0) at 95°C for 20 min. Sections were then washed, equilibrated in water for 1 hour, and blocked with 1% BSA for 1 hour at room temperature. Samples were then incubated in primary antibody overnight at 4°C, washed 3× in distilled water, and then incubated in secondary antibody for 2 hours at room temperature. Samples were then washed 3× in distilled water, mounted with SlowFade Diamond Antifade mounting medium with 4′,6-diamidino-2-phenylindole (Thermo Fisher Scientific), and imaged using an SP8 Leica laser scanning confocal microscope (Leica Microsystems Inc., Buffalo Grove, IL).

Statistical analyses

Data statistics and bioinformatics analyses were performed using R (version 4.0.3) and Bioconductor packages (www.bioconductor.org/) unless otherwise indicated. Multiple comparisons were adjusted by the Benjamini-Hochberg correction.