Publicly available Gene Expression Omnibus (GEO) database GSE59983 (https://www.ncbi.nlm.nih.gov/geo/) [19] from 76 retinoblastoma tissue samples, collected from patients who underwent primary enucleation without receiving previous treatment and profiled with Affymetrix human genome u133 plus 2.0 PM microarray, were analyzed using R package GEOquery [20]. Differentially expressed genes (DEGs) among primary tumors (N = 72) vs normal retinal (N = 4) were matched with previously obtained data of DAC-driven co-regulated genes from the WERI-Rb-1 cell line profiled using PIQORTM Cell Death Human Sense Microarrays [18].

GeneMania algorithm (http://www.genemania.org) was used to generate network maps of DEGs showing co-expression among connected genes. Functional enrichment analysis was obtained by the DAVID tool [21]. The resulting networks were exported to the Cytoscape platform (https://cytoscape.org/) for gene mapping and proper visualization of highlighted DEGs and their relationships. Gene ontology analysis for functional annotations of the differentially expressed genes was carried out using the Biological Networks Gene Ontology (BiNGO) and DAVID tools [21, 22]. Cytoscape was also used to properly visualize the connections among genes and their biological functions.

Public available single-cell transcriptome data from GEO database GSE142526 [23] of normal retina cells (human retinal organoids ORG_D104; ORG_D110 and retinospheres derived from human fetal retina RS_D134_pl_26FV) and from GEO database GSE196420 [24] of patient-derived non-familial retinoblastoma cells classified as E and D (Supplementary Table 1) according to the intraocular retinoblastoma classification (IIRC) (wRB6, RB006, RB010, RB015, RB016, RB018, RB020, and RB021) were also analyzed using the Seurat package [25] with a state-of-the-art pipeline previously utilized in [26, 27]. Single-cell transcriptomes data were filtered based on defined criteria when genes expressed in cells $>$200 and number of RNA read counts were within 300 and 100,000 with mitochondrial percentage in cells 10. We also regressed cell cycle scores [28]. Clustering was performed using the Louvain algorithm, and Uniform Manifold Approximation and Projection (UMAP) visualization was generated [29]. The analyzed single-cell data (https://ankushs.shinyapps.io/Retinoblastoma_GSE196420/) is made available through the ShinyCell framework [27]. Differential expression of genes was computed using the R package GEOquery [20].

The Weighted Gene Co-expression Network Analysis package (WGCNA) [30, 31] was used to reconstruct weighted gene co-expression networks for the differentially expressed genes in primary tumor and normal retinal cells. Edge weights, computed based on topology overlap measures, assigned co-expression correlation values between 0 and 1 to two connected genes (GEOquery). Subnetworks of 15 retinoblastoma hub genes were extracted from the co-expression networks of primary tumors using a threshold cutoff of 0.05 on edge weights among co-expressing genes. Cytoscape version 3.3 was used to visualize, and topological parameters were computed using Centiscape [27].

WERI-Rb-1 cells (ATCC, Manassas, VA, USA) were maintained in RPMI1640 medium (EuroClone, Milan, Italy) supplemented with 10% fetal bovine serum (EuroClone, Milan, Italy) and 2 mM L-glutamine, plus penicillin (100 U/mL) and streptomycin (100 mg/mL). The manufacture guarantee the cells are Mycoplasma Free and authenticated. Cells were split twice a week by resuspension in fresh media at a concentration of 3 $\times$ 10⁵ cells/mL. For treatments, cells were seeded at a density of 2 $\times$ 10⁵ cells/mL in six-well plates. After 24 h 2.5 µM 5-Aza-2^′-deoxycytidine (Merck KGaA, Darmstadt, Germany) was added to the culture medium of the treated cells and maintained for up to 96 h. Control cells were treated in parallel with the vehicle.

Total RNA was extracted from WERI-Rb-1 cells using the NucleoSpin RNA isolation kit (Macherey-Nagel, Duren, Germany) and from subcutaneous xenograft tumor samples using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA), according to the manufacturer’s instructions. RNA concentration and purity were determined using a NanoDrop spectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA). For each sample, 1 µg of total RNA was reverse transcribed using the Maxima H Minus First Strand cDNA Synthesis Kit (Thermo Fisher Scientific, Waltham, MA, USA). For the selected set of genes, qPCR validation was performed using the DyNAmo Flash SYBR Green qPCR Kit (Thermo Fisher Scientific, Waltham, MA, USA) with the PikoReal Real-Time PCR System (Thermo Fisher Scientific, Waltham, MA, USA). Primer specificity was checked using primer-BLAST (https://www.ncbi.nlm.nih.gov/tools/primer-blast) and confirmed by melting curve analysis. Primer pairs sequences are reported in Supplementary Table 2a. Amplification conditions were as follows: 7 minutes at 95 °C, followed by 40 cycles of 10 seconds at 95 °C, 20 seconds at 60 °C and 20 seconds at 72 °C. All samples were analyzed in triplicate. The relative expression of target genes was evaluated using the comparative cycle threshold method ($\mathrm{\Delta }$$\mathrm{\Delta }$CT), with $\beta$-actin (ACTB) used for normalization.

DNA methylation in the promoter regions of selected genes (CASP8, FAS, BIK, TP73, DAP3, and RRAD) was analyzed by MSP, based on the differences between methylated and unmethylated DNA sequences, following sodium bisulfite treatment. Genomic DNA was extracted from WERI-Rb-1 cells and subjected to bisulfite modification using the EpiJET Bisulfite Conversion Kit (Thermo Fisher Scientific, Waltham, MA, USA). Modified DNA was then used for PCR reactions. Primer pairs for methylated (M) and un-methylated (U) sequences were designed to target CpG islands in promoter regions of the selected genes, using sequence data obtained from UCSC GENOME Browser (http://genome.ucsc.edu/cgi-bin/hgGateway). Primer specificity was validated through agarose gel electrophoresis of the PCR products, confirming the presence of a single band with the expected molecular weight. Detailed primer sequences are reported in Supplementary Table 2b.

Experiments were conducted on opportunistic pathogen-free the Naval Medical Research Institute (NMRI) male athymic BALB/c Nude mice, aged six to seven weeks (Harlan Laboratories, Udine, Italy), in accordance with EU Directive 2010/63/EU and the regulation of the Italian Ministry of Health. The mice were maintained on standard laboratory food and water that were available at all times, under a 12 h artificial light/dark cycle. After induction of deep anesthesia by inhalation of isoflurane in an induction chamber, the mice were euthanized using CO₂, as recommended by attachment IV Table 3 of EU Directive 2010/63/EU; 14G00036, GU Serie Generale n.61 del 14-03-2014, and Italian Ministry of Health. All the procedures were verified by the Ethics Committees of the Toscana Life Sciences and the Istituto Superiore di Sanità (ISS) on behalf of the Italian Minister of Health (Permit Number: # CNR-030314 and # CNR-101013) following ethical ICLAS and ARRIVE guidelines.

For the subcutaneous implants, animals were anesthetized using 2.5% isoflurane during the manipulation. WERI-Rb-1 cells, at a concentration of 3.6 $\times$ 10⁷ in 100 µL 1$\times$ PBS, were injected subcutaneously in a 1:1 ratio with Matrigel TM basement membrane matrix (BD Biosciences, Franklin Lakes, NJ, USA) into the left flank of each mouse (total volume 200 µL). Once the grafts became palpable, their volume was measured using a digital caliper (length $\times$ width²/2), and animals were assigned to experimental groups through minimization following the ARRIVE guidelines [32] to start the treatments. Tumor volume (mm³) was measured biweekly and confirmed by ultrasound imaging (VisualSonics VEVO 2100, Fujifilm Sonosite Inc., Bothell, WA, USA) in the last session of measurement.

The orthotopic retinoblastoma model was established unilaterally in one eye by injecting 1 $\times$ 10⁴ WERI-Rb-1 cells in 10 µL of PBS. Briefly, under an operating microscope and using a Hamilton syringe 32 gauge, (cod: 12301828 Fisher Scientific Italia, Segrate (MI), Italy) the right eye globe of anesthetized animals was punctured laterally, through the conjunctiva and sclera, to access the vitreous cavity. Ultrasound Imaging (VisualSonics VEVO 2100, Fujifilm Sonosite Inc., Bothell, WA, USA) was used to establish the tumor appearance. Treatment began before visible leukocoria (white reflex in the eye pupil) appeared in the affected eye. Volume measurements were performed for group assignment (baseline) and then once a week to minimize animal distress due to repeated anesthesia.

In both experimental settings (xenograft and orthotopic retinoblastoma models), the treated groups received biweekly I.V. injections of 300 µL of 75 µg DAC in PBS buffer (meaning a therapeutic dose of 2.5 mg/kg), while the control groups received the vehicle. At the end of each measurement and treatment session, the animals were monitored for signs of distress and allowed to recover in the original cages. After 3 weeks of treatment, the mice were euthanized via CO₂ inhalation. No signs of distress were observed [32] and no animals died during the treatments. Resected tumor masses from subcutaneous xenografts were processed using TRIzol reagent (Thermo Fisher Scientific, Waltham, MA, USA) and stored at –20 °C for further qPCR analysis.

The VisualSonics Vevo 2100 imaging system (VisualSonics VEVO 2100, Fujifilm Sonosite Inc., Bothell, WA, USA) was utilized to measure retinoblastoma growth within the eye. Animals were anesthetized with 5% isoflurane at an oxygen flow rate of 2 L/min (maintained at 2.5% isoflurane at an oxygen flow rate of 2 L/min) and placed on a warming pad in a prone position to facilitate signal acquisition while monitoring temperature, respiratory and heart rate. B-mode images were acquired using the MS-550 Blue transducer (central frequency, 40 MHz) connected to a 3-dimensional motor collecting frames 0.5 mm apart in the eye region. In the orthotopic model, for off-line analysis, a field of interest (FOI) outlining the tumor boundaries in the eye was drawn for reconstruction in 3D B-mode. These results were normalized to baseline (volume at tumor appearance) and expressed as fold increase (V = Vtx/Vt0) $\pm$ sem.

For the in vitro experiments the results represent the mean $\pm$ SEM of at least three independent experiments. Two-Way ANOVA was used to evaluate the impact of decitabine on cell cycle phases at the different time points (24, 48, and 72 h). Similarly, Two-Way ANOVA was used to describe the significance of the longitudinal effect of the treatment on the number of apoptotic cells (mean $\pm$ SEM). For qPCR data, $\mathrm{\Delta }$Ct values were subjected to a one-sample t-test and are shown as mean $\pm$ SD. To analyze the effect of time and treatment (DAC) on the tumor growth respectively in subcutaneous and orthotopic Rb xenograft models we applied two-way ANOVA repeated measure analysis for both time and treatment parameters (recommended Sidak analysis; GraphPad-Prism 7.0e, GraphPad Software, Inc., Boston, MA, USA). Results are expressed as mean $\pm$ SEM.

To assess whether the co-regulated clusters of genes (Supplementary Fig. 1), previously identified in DAC-treated WERI-Rb-1 cells [18], could represent candidate markers and potential therapeutic targets of DAC for primary retinoblastoma tumors, we performed a comparative analysis of differentially expressed genes (DEGs) in primary tumors versus WERI-Rb-1 cells. Publicly available transcriptome datasets from patient-derived samples and low-passage patient-derived single-cells were analyzed [19, 23, 24].

Table 1 reports the most significant time-related DEGs in WERI-Rb-1 cells after exposure to DAC within the co-regulated genes clusters, previously identified through the analysis of PIQORTM Cell Death Human Sense Microarrays data [18].

A strong down-regulation of pro-survival signals is observed at a short latency (48 h) following DAC treatment. These signals include several pro-mitogenic mediator surface molecules (RASL10B, TOLLIP, and TRAF2), intracellular MAPK members (MAPK8IP1, MAPK8IP2, MAPK12, IRAK1, IKKB, and NFKB3), and intracellular anti-apoptotic mediators (BCL2L1, MKL1, JUNB, BIRC5, and FOS). Conversely, genes playing a role in cell cycle arrest, like PPP1R15A, STAT1, and CDKN1A, are overexpressed starting 48 h after DAC treatment. Furthermore, significant activation of pathways associated with pro-apoptotic signals becomes evident 72 h after DAC treatment. In particular, a significant over-expression of several apoptotic mediator surface antigens, such as FAS, DAP3, and TRAMP receptors as well as intracellular pro-apoptotic members like PSMD2, ALG2, and GAPD, along with members of the caspases cascade (e.g., CASP8, CASP3, and CASP6) is observed. Concurrently, the p53-dependent pathway, which includes CDKN1A, CDC10, PIG8, and BAX, is activated, indicating a potential triggering of cytochrome C release in combination with other mechanisms (e.g., p73-dependent BCL2L1 repression via BIK). The simultaneous down- and up-regulation of all these pathways may explain the effective anti-cancer activity of DAC in WERI-Rb-1 cells.

Interestingly, a notable overexpression of genes associated with DNA repair signaling, such as PARP1 and LIG4 (DNA LIGASE IV), along with tissue-specific genes like PPM1D (a photoreceptor-related gene) and AKAP12 (a gene related to hemato-retinal barrier), both involved in visual signals [19, 33], was observed at a later time point (96 h), suggesting a possible reprogramming action of DAC in the surviving cells. To investigate the potential of these DAC-driven genes as markers of primary retinoblastoma tumors, a gene expression profile analysis was performed comparing whole transcriptome microarray datasets of tumor cone photoreceptor lineages vs normal retinal cells. Both cell types were derived from frozen tissue samples collected after primary enucleation in patients with familial retinoblastoma [19]. The analysis identified 463 up-regulated genes, with a Log fold change greater than 2, and 282 down-regulated genes, with a Log FC $\le$–2 (see Supplementary Table 3).

Among the 463 up-regulated genes, 136 are mainly linked to biological processes involved in cell proliferation, including RASL10B, TRAF2, RELA, FOS, BCL2L1, MAPKs, and NFKBs. These genes were found to be highly down-regulated in WERI-Rb-1 cells after DAC treatment. Additionally, 24 genes that act as positive regulators of DNA damage and repair signaling, including PARP1 and LIG4, were up-regulated in the surviving WERI-Rb-1 cells after 96 h of DAC treatment. Interestingly, in the cone photoreceptor lineages of primary tumors, the top up-regulated genes (Log FC $\ge$4) included those playing a key role in tumor progression, such as BIRC5 oncogene. We found these genes highly down-regulated at early time points in WERI-Rb-1 cells following DAC treatment. Among the 282 down-regulated genes in primary tumors, more than 200 are associated with cell cycle arrest and pro-apoptotic signals, such as FAS, BIK, CASP8, and CASP6. Additionally, 27 of these genes are functionally related to retina development, visual perception, and rhodopsin-mediated signaling pathways. Examples include PPM1D, a photoreceptor-related gene, and AKAP12, which plays a role in the hemato-retinal barrier [20, 21]. These genes were significantly up-regulated in WERI-Rb-1 cells after DAC treatment, suggesting that DAC, as an epigenetic drug, may help restore retinal function by reversing the epigenetic silencing of key factors involved in retina development. A network map of DEGs common to WERI-Rb-1 cells and primary tumors was generated to identify potential hub genes and their crosstalk among pathways. The analysis revealed 15 hubs/driver genes that are tightly interconnected, sharing numerous pathways primarily associated with the regulation of TNF-, FAS-, p53-dependent apoptotic signaling, and NF-$\kappa$B pro-survival activity (Fig. 1).

In addition, we developed a protein-protein interaction network map that highlights the central role of these 15 hub genes in regulating numerous proteins/factors (Supplementary Fig. 2). Notably, the oncogenes TRAF2, RELA, and BIRC5, which play pivotal pro-survival roles and exhibit numerous interactions with their neighboring proteins, are highly up-regulated in primary tumors but down-regulated early after DAC treatment. Conversely, FAS and CASP8, which are involved in regulating extrinsic apoptosis pathways, along with BAX, BCL2L1, HSF, and CASP6, key players in intrinsic apoptosis, are highly up-regulated in DAC-treated WERI-Rb-1 cells and are expressed at lower levels in primary tumors. These genes exhibit numerous interactions with proteins/factors that regulate cell fate. Since bulk RNA-seq studies are primarily designed to provide general transcriptomic data for entire tissue samples, they inherently do not capture the tumor heterogeneity of Rb [34]. Therefore, single-cell RNA sequencing (scRNA-seq) analysis was applied to verify the central role of the 15 hub/DAC-driven genes in the genesis and malignant progression of retinoblastoma. Two publicly available single-cell transcriptome datasets [23, 24], comprising human normal retina cells (organoids and retinospheres) and patient-derived sporadic retinoblastoma cells of different intraocular retinoblastoma classification (IIRC) were analyzed (Supplementary Table 1). Comparative computational analysis at the single-cell level identified 25 subclusters annotated with highly expressed genes (see https://ankushs.shinyapps.io/Retinoblastoma_GSE196420/). For each set of normal cells (ORG_D104; ORG_D1110; and RS_D134_pl_26FV) and primary retinoblastoma cells (RB025; RB026; RB027; RB028; RB029) (Fig. 2a) we describe the expression profiles of DEGs, including the 15 hub genes. UMAP images in Fig. 2b,c illustrate the expression of two representative genes, namely BIRC5 and CASP8.

Most of the key genes found to be up- or down-regulated in WERI-Rb-1 after DAC treatment, as listed in Table 1, show an inverse correlation in expression compared to those identified in subclusters of patient-derived single-cells (Fig. 2d). Among the 15 hub genes shown in Fig. 1, the BIRC5 and TRAF2 genes, are mostly overexpressed in primary tumor cells (RB025, RB026, RB027, RB028) and are highly down-regulated in DAC-treated WERI-Rb-1 cells. Similarly, DAXX and HSF1, which are down-regulated after DAC treatment, are primarily overexpressed in primary tumors (e.g., RB028) and in fetal retinas (RS-D134_lp_26FV). Meanwhile, the regulation of PPP1R15A, RRAD, FAS, CASP6, CASP3, BAX, CASP8, and BIK genes, mainly involved in cell growth arrest and apoptotic signaling and significantly overexpressed in DAC-treated WERI-Rb-1 cells, are expressed at low levels in organoids (ORG_D104 and ORG_D110) and in all primary tumor samples (RB025, RB026, RB027, RB028, and RB029) (Fig. 2d). The overall data suggest that the epigenetic alterations may be involved in tumor progression and that the identified 15 hub/driver genes could serve as key targets for DAC treatment in primary retinoblastoma tumors.

To test the robustness of the computational data analysis, quantitative PCR was performed on the selected set of 15 hub/driver genes, which showed high differential expression in WERI-Rb-1 cells after DAC epigenetic treatment (Fig. 3).

We noted a significant up-regulation of several pro-apoptotic genes, namely FAS, CASP8, BIK, and the tumor suppressor gene RRAD, at all time-points analyzed (p 0.05 at 48, 72, and 96 h). Gene expression changes were observed for other pro-apoptotic genes only at certain time points, such as HSF1 (p 0.05 at 72 h), and DAP3, TP73, PSMD2, and CASP6 (p 0.05 at 48 and 72 h). It is also noteworthy that other genes, such as BAX and BCL2L1, which encode apoptotic activators that regulate mitochondria membrane potential, were significantly overexpressed only at the latest time point (p 0.05 at 96 h). On the other hand, the BIRC5 oncogene underwent early significant down-regulation that was maintained over a long period (p 0.05 at 48 and 72 h), suggesting an apoptotic response. Our qPCR analysis of DAC-treated WERI-Rb-1 cells confirmed the overall trend of mRNA expression regulation observed in the previous microarray analysis (Table 1).

To assess the potential effect of DAC on modulating the expression of key genes, we employed Methylation-Specific PCR (MSP) to examine the epigenetic status of selected up-regulated genes, that are commonly epigenetically silenced in tumors. As shown in Fig. 4, changes in DNA methylation patterns in CpG-rich regions of promoters were observed for the CASP8 and BIK genes. For both genes, the promoter methylation levels (M) in untreated samples (CTRL) shifted to unmethylated patterns (U) following DAC treatment. Interestingly, for BIK, two bands of similar intensity appeared in the control sample: one corresponding to the methylated pattern and the other to the unmethylated one, indicating the possibility of allele-specific methylation, a common occurrence in tumors. In contrast, the CpG islands of FAS, TP73, DAP3, and RRAD gene promoters appeared unmethylated in both control and DAC-treated cells. Overall, these results suggest that DAC modulates gene expression in WERI-Rb-1 cells through a dual mechanism: direct, locus-specific changes in CpG islands methylation, as seen for CASP8 and BIK, and an indirect or nonspecific effect leading to the activation of downstream effectors, without altering cytosines methylation, as previously suggested in the literature [29].