MYCN-induced nucleolar stress drives an early senescence-like transcriptional program in hTERT-immortalized RPE cells

MYCN is an oncogenic driver in neural crest-derived neuroblastoma and medulloblastoma. To better understand the early effects of MYCN activation in a neural-crest lineage context, we profiled the transcriptome of immortalized human retina pigment epithelial cells with inducible MYCN activation. Gene signatures associated with elevated MYC/MYCN activity were induced after 24 h of MYCN activation, which attenuated but sustained at later time points. Unexpectedly, MYCN activation was accompanied by reduced cell growth. Gene set enrichment analysis revealed a senescence-like signature with strong induction of p53 and p21 but in the absence of canonical hallmarks of senescence such as β-galactosidase positivity, suggesting incomplete cell fate commitment. When scrutinizing the putative drivers of this growth attenuation, differential gene expression analysis identified several regulators of nucleolar stress. This process was also reflected by phenotypic correlates such as cytoplasmic granule accrual and nucleolar coalescence. Hence, we propose that the induction of MYCN congests the translational machinery, causing nucleolar stress and driving cells into a transient pre-senescent state. Our findings shed new light on the early events induced by MYCN activation and may help unravelling which factors are required for cells to tolerate unscheduled MYCN overexpression during early malignant transformation.

gene regulation in telomerase immortalized RPE1 cells (further referred to as RPE1-MYCN-ER) 6 . In this model, chaperonins persistently sequester the chimeric MYCN-ER protein in the cytosol. The oestrogen receptor ligand 4-OHT will displace the chaperonins, causing MYCN-ER to translocate to the nucleus where it will exert its transcriptional activity. Western blotting confirmed the expected decrease in the cytoplasmic fraction of MYCN-ER protein (Fig. 1a), and concomitant increase in nuclear levels 24 h after the addition of 4-OHT to RPE1-MYCN-ER cells (MYCN-ON) as compared to mock-treated RPE1-MYCN-ER (MYCN-OFF) and wild type RPE1 (WT) cells (Fig. 1b).
First, we evaluated the functionality of MYCN after nuclear translocation in RPE1 cells by measuring the transcriptional impact. To this end, bulk RNA sequencing was performed at three sequential time points post induction (p.i.) ( Table 1). Principal component analysis (PCA) revealed a clear separation of the MYCN-ON cells from the other conditions along the first principal component (PC) (Fig. 1c). Within the MYCN-ON samples, there was also a pronounced separation of individual time points, predominantly dictated by the second PC. This demonstrates that incubation with 4-OHT induces a specific and time-dependent transcriptional shift in MYCN-ON cells (Fig. 1c).
To assess the effect of MYCN activation on the RPE1 transcriptome, we performed gene set enrichment analysis (GSEA) using 30 publicly available MYC/MYCN target gene sets. This revealed the expected enrichment (FDR < 0.01) of MYC/MYCN targets in MYCN-ON cells at 24 h and onwards while this was not visible for the MYCN-OFF cells (Fig. 1d). In support thereof, RNA-sequencing and qPCR indeed revealed an upregulation of several selected MYCN target genes (e.g., ODC1, PTMA, WDR12, MDM2, NPM1, DKC1 7,8 (Supplementary Fig. 1a,b). When evaluating MYCN activity (translocation and signature activation) over time, the highest activity was noted at the 24 h time point, followed by sustained but attenuated activity at the later 48 h and 72 h time points (Fig. 1b, Supplementary Fig. 1c). This may be the result of a negative feedback loop as previously described 9 , a notion supported by the lower MYCN transcript levels in MYCN-ON cells as compared to MYCN-OFF cells at all time points ( Supplementary Fig. 1d).
MYCN activation induces a p53-p21 driven growth attenuation. We noticed that MYCN-ON cells did not reach the same level of confluency 72 h p.i. as compared to controls (Fig. 2a). This pointed to an attenuation in cell proliferation, which was confirmed with a focus formation assay (P-value < 0.01) (Fig. 2b). Quantitative immunofluorescence microscopy revealed that MYCN-ON cell populations harboured consistently fewer replicating (EdU-positive) cells than controls and progressively fewer cycling (Ki-67-positive) cells with time up to 72 h after MYCN induction (Fig. 2c-e). Caspase-3 and -7 activity was not significantly increased at any of the measured time points in MYCN-ON cells, suggesting there was no induction of apoptosis (Supplementary Fig. 2). DAPI-based cell cycle staging showed an increased fraction of cells in G1-phase and a significant decrease in G2-phase in MYCN-ON cells starting at 48 h p.i. (Fig. 2f).
To gain insight into the underlying cause of growth attenuation, we analysed accompanying changes in the transcriptome of MYCN-ON cells. Using gene set enrichment analysis (GSEA), we found enrichment for predicted p53 promotor-binding sites in the upregulated genes while the downregulated genes were enriched for FOXM1-and E2F1-binding targets, 48 h and 72 h after MYCN induction (Fig. 3a,b). Positive regulators of cell cycle transcription factors E2F1 and E2F2 were strongly downregulated after 24 h (Supplementary Fig. 3a). At all tested time points, transcriptional upregulation was found for Cyclin Dependent Kinase Inhibitor 1A (CDKN1A), encoding the p21 protein and a direct target gene of p53. Induction of p21 by p53 is implicated in cell cycle arrest and senescence 10,11 in accordance with the observed growth attenuation ( Fig. 3c; Supplementary  Fig. 3b). The p21 and p16 cyclin-dependent kinase inhibitors are known to play key roles in the onset of cellular senescence as their simultaneous induction cooperatively blocks the activation of cyclin D kinases 12 . While cyclin dependent kinase inhibitor 2A (CDKN2A encoding the p16 protein) expression was initially downregulated at 24 h p.i., we observed significant upregulation in MYCN-ON cells at the later 48 and 72 h p.i. time points. Of further interest, expression of cyclin dependent kinase 2 (CDK2) which has been reported to act as suppressor of oncogene induced senescence 13 is significantly and progressively downregulated on a transcriptional level ( Supplementary Fig. 3c). Quantitative immunofluorescence of MYCN-ON versus MYCN-OFF cells, revealed www.nature.com/scientificreports/ increased p53 and p21 protein (Fig. 3d,e), in the absence of notable TP53 mRNA changes ( Supplementary  Fig. 3d). Deconvolution of nuclear marker intensities along the cell cycle revealed that p21 was much more abundant in G1-phase as compared to S or G2-phase ( Supplementary Fig. 3e). This was already evident prior to MYCN induction (i.e., in MYCN-OFF), but became even more pronounced in MYCN-ON cells, suggesting this population to selectively enrich. In contrast, p53 protein levels increased in MYCN-ON cells independently of cell cycle status ( Supplementary Fig. 3e,f).
Taken together, MYCN activation induces p53-p21 transcript and protein levels along with upregulation of p53 target genes and attenuated expression of FOXM1 and E2F target genes.
Induction of a pre-senescent transcriptional program in MYCN-ON cells. The combination of a marked growth attenuation and upregulation of p53 and p21 pathways, could be indicative for the induction of MYCN-triggered (oncogene-induced) senescence. To investigate this, we queried our RNA-seq dataset for a variety of senescence-related signatures (Fig. 4a) and observed a strong enrichment as of the 48 h time point after MYCN induction for the cell cycle arrest and senescence associated gene sets (Fig. 4a). In this context, the enrichment for a YAP/TAZ gene signature 14,15 (Fig. 4a,b) is highly relevant given the reported downregulation of YAP proteins following replication-induced cellular senescence in IMR90 cells 16 . YAP/TAZ signaling has also been shown to control RRM2 transcription as a crucial effector in senescence control which is confirmed in our data showing a steep decline of RRM2 expression levels at 48 h p.i. (p < 0.0001, ANOVA) (Fig. 4c,d). In further support of a senescence inducing gene signature, we observed downregulation of the senescence-associated biomarker LMNB1 both at the transcriptional (p = 6.51e−08, ANOVA) and (at 72 h) at protein level (Fig. 4c,  www.nature.com/scientificreports/ cells, we also discovered signatures with bimodal distributions. In particular, two core senescence signatures 17,18 (p < 0.05) and a SASP-associated signature 18 (p < 0.01) showed a marked drop at the 24 h time point in MYCN-ON cells, only to return to baseline (MYCN-OFF) levels at later time points. A similar observation was made for gene sets that become upregulated upon CDK4/CCND1 knockdown (p < 0.001) (Fig. 4a).
To dissect the different patterns of temporal gene expression changes following MYCN activation in more detail, we applied the elbow method for k-means clustering which allows the identification of an unbiased number of clusters in a given data set 19 . Using this approach, a total of 10 clusters was identified that captured the overall variation in the bulk RNA-seq data. Four clusters (3, 6, 7 and 8) were connected to MYCN signatures (Fig. 5a, Supplementary Fig. 4a, Supplementary Table 2), while cluster 10 was strongly enriched for gene targets of cell cycle master regulators FOXM1 and E2F. Three clusters (3, 7 and 8) consisted of gene sets associated with ribosome biogenesis, a well-known MYC(N) driven process that is assumed to allow cells to cope with an increased demand for proteins required for cell growth and proliferation 20 (Figs. 5a, S5a). These gene sets include EBNA1BP2 (EBP2) (cluster 3), which stabilizes MYC and promotes MYC-mediated rRNA synthesis and cell proliferation, as well as genes encoding proteins involved in ribosome RNA transcription (POLR3G), processing and maturation (WDR12, RRP1B) and ribosome assembly and shaping (MDN1) (cluster 8). Genes that composed cluster 7 exhibited increased expression 24 h after MYCN activation and robust downregulation at the 48 h time point. A prominent representative of this cluster was NPM1. NPM1 controls the impaired ribosome biogenesis checkpoint which becomes activated upon unscheduled MYC(N) driven ribosome production 21 . Checkpoint activation disables MDM2, leading to p53 stabilisation and p53 pathway activation, including p21 induction. Thus, the observed NPM1 overexpression during the very early phase of MYCN activation (24 h) suggests nucleolar stress is a driving factor in the p21-driven growth delay.
Finally, in order to gain insight into the critical transcriptional changes and overall heterogeneity of transcriptional induction of CDKN1A and possible other early markers heralding the senescence transcriptome signature, we performed single cell RNA sequencing. As CDKN1A was not yet significantly upregulated 16 h p.i., but was at 24 h p.i. (Supplementary Fig. 4b), we selected two time points around the CDKN1A inflection timepoint (22 h and 24 h). Figure 5b shows the inter-cellular heterogeneity of CDKN1A expression levels compared to untreated cells further supporting our immunohistochemistry staining (Fig. 3d, Supplementary Fig. 3f). When correlating CDKN1A expression levels within cell cycle stage, the higher expression levels were proportionally increased in all stages (Fig. 5c). This indicates that CDKN1A expression is increasing overall, already 24 h p.i.

Growth attenuation in MYCN-ON cells is accompanied by cytoplasmic granules and nucleolar coalescence.
The above data prompted us to further investigate known senescence hallmarks. Unexpectedly, no evidence for ß-galactosidase staining was found at any of the time points following MYCN activation studied for transcriptome profiling. Likewise, nuclear DAPI staining did not present evidence for senescenceassociated heterochromatin foci formation (SAHFs), another canonical phenotypic senescence marker 22,23 (Supplementary Fig. 5a), and γH2AX and RPA32 staining did not reveal increased levels of DNA damage in MYCN-ON cells ( Supplementary Fig. 5b,c). However, two obvious phenotypic changes were notable in MYCN-ON cells. First, we noticed a distinct nucleolar coalescence in the DAPI images of MYCN-ON cells. Quantitative immunofluorescence for fibrillarin (a protein component of C/D box small nucleolar ribonucleoproteins), confirmed this observation together with an overall increased nucleolar size (Fig. 6a,b), independently of the cell cycle status (Fig. 6c). Furthermore, a second marked change observed by transmission microscopy was a conspicuous cytoplasmic granularity in MYCN-ON cells (Fig. 6d). We reasoned that this might reflect a change in protein turnover as a consequence of the upregulated ribosomal biogenesis and nucleolar stress as also suggested by the strong NPM1 induction. However, neither staining for lysosomes, stress granules or P-bodies 24,25 revealed colocalization with these granules, suggesting they do not represent typical relics of aberrant protein accrual and leaving their nature as yet undetermined (Supplementary Fig. 5d-g).
Taken together, we propose that unscheduled MYCN activation causes ribosome biogenesis and nucleolar stress in G2-phase, pre-priming cells for faltering in subsequent G1 and slowing down cell proliferation.

Discussion
Normal cells have tightly regulated signalling pathways that protect against insults such as unscheduled oncogene induction. Understanding how major oncogenic drivers such as MYCN activate and bypass such pathways during tumor initiation offers insights that may be of use for development of novel therapeutic strategies. Hence, we explored the use of immortalised neural crest-derived retinal pigment (RPE) cells to investigate the early effects of MYCN activation.
Unexpectedly, following MYCN activation, we observed attenuated cell growth without notable increased apoptosis but with a pronounced induction of p53 and p21. Attenuated proliferation was further accompanied by an increased number of cells in G1-phase and the induction of several previously reported senescence-induced gene signatures. While we observed a robust downregulation of LMNB1, other hallmarks of senescence (such as β-galactosidase activity and SAHFs) were not present, suggesting that cells were not undergoing full senescence. We also noticed progressive CDK2 transcriptional downregulation upon MYCN activation. CDK2 was previously shown to be essential to avoid MYC induced senescence in MEFs and required for sustained p21 and p16 induction 13 . The lack of strong p16 induction in the MYCN-ON cells might expmain the lack of full-blown induction of senescence in our model. The exact role and interaction of p21, p16 and CDK2 during MYCN activation remains to be explored, e.g. through CDK2 knock down and overexpression to further clarify the role of CDK2 in relation to p16 and p21 in the MYCN induced dynamic phenotypic changes in our model system.
Remarkably, MYCN-ON cells did show distinct phenotypic changes including nucleolar coalescence and cytoplasmic granularity, which aligned with transcriptional evidence for upregulation of ribosome biogenesis. www.nature.com/scientificreports/ Enlarged and prominent nucleoli upon increased MYCN activity have been documented before 26 , but we now connect it to a transcriptional and phenotypic switch. We propose that ribosome biogenesis and the consequent nucleolar stress drives cells into what we propose to represent a pre-senescent state. Indeed, perturbed ribosomal biogenesis and nucleolar stress is known to drive nucleolar coalescence in senescent cells [27][28][29][30] and nucleolar expansion has been observed in progeria cells as a result of increased translational throughput (ribosome biogenesis and protein synthesis) 31 . This could be mimicked by knockdown of LMNA and overexpression of progerin, suggesting a general influence of the lamin network on nucleolar plasticity. It is therefore tempting to speculate that the loss of lamin B1 maybe thus also act upstream of this process.
Oncogene activation (such as MYCN) is well known to induce alterations in ribosome biogenesis and ultimately activate the IRBC (Impaired Ribosome Biogenesis Checkpoint) as ribosomal proteins interact with the central acidic domain of MDM2 and indirectly promotes the stabilization of p53 27 . Interestingly, MYCN is implicated in upregulation of RPL (Ribosomal Protein Large) and RPS (Ribosomal Protein Small) proteins in neuroblastoma cells 20 . Based on our observations, we envision that the senescence-like response triggered by MYCN activation in RPE1-MYCN-ER cells reflects a dynamic response to the oncogenic stimulus. Intriguingly, despite robust induction of CDKN1A in the MYCN-ON cells, neither apoptosis nor complete growth arrest was observed, and cells continued to grow, albeit more slowly. In this context, the recent discovery of a so-called "Goldilocks" zone for p21 levels was observed to control the proliferation-senescence cell-fate decision after drug treatment 32 . Either a delayed or acute drug-induced p21 response led to senescence, while an intermediate p21 pulse enabled sustained proliferation. The cell-cycle dependent p21 overshoot that we witnessed in MYCN-ON cells may thus reflect an attempt to initiate cell cycle arrest, which was not completely successful, either by having insufficient intensity or improper timing.
In conclusion, we have identified a transient, pre-senescent state as early response to unscheduled MYCN activation. Further studies, including live cell imaging in this model and in vivo single cell analysis of early emerging hyperplastic lesions in MYC(N) transgenic animal tumor models can shed further light on the dynamic changes in gene activity and cellular processes that allow these cells to cope with the stress during early malignant transformation. These insights can also provide guidance towards critical dependencies and pathways as novel drug targets for MYC(N) driven malignancies, many of which are currently still associated with poor survival.

Materials and methods
Cell culture. The hTERT-immortalized MYCN retinal pigmented epithelial cell line (RPE1-MYCN-ER, kindly provided by Michael Hogarty, Children's Hospital of Philadelphia, Philadelphia, PA) carries a MYCN:ER expression construct that results in its constitutive expression as described by Maris et al. 6 . When 4-OHT is present in the culture media, MycN:ER translocates to the nucleus where it is transcriptionally active. The wilde-type (WT) hTERT-immortalized retinal pigment epithelium (hTERT-RPE1, Clontech, Palo Alto, CA) and RPE1-MYCN-ER were cultured in DMEM/F-12 HEPES (31330-038; Invitrogen), supplemented with 10% fetal bovine serum (FBS), 1% penicillin/streptomycin (15140-122; Invitrogen), 1% l-glutamin (25030-025; Invitrogen) and 0.01 mg/ml Hygromycin B (10687010; Invitrogen) according to standard procedures. Moreover, 1 µg/ ml puromycin (P7255; Sigma-Aldrich) was added to the RPE1-MYCN-ER medium composition as a selection marker for the 4-OHT inducible construct. Mycoplasma testing was done on a monthly basis. Using DNA sequencing of the KRAS coding sequence, we validated the presence of the p.G12_G13insAG variant reported by Di Nicolantonio et al. 33 in RPE-hTERT cells. We classified this as a variant of unknown significance and noted very low expression levels (not shown). In keeping with this presumed lack of functional significance Benedict et al. 34 could provide evidence for KRAS activation using a pull-down assay to specifically precipitate RAS-GTP. Therefore, we are confident that the RPE-hTERT immortalized cells used in our study are valid to test cellular and molecular effects of MYCN induction.  . MYCN activation triggers a p53-p21 driven growth attenuation. (a) Enrichr (https:// maaya nlab. cloud/ Enric hr/) analysis shows enrichment for predicted p53 promotor binding sites in the upregulated genes and enrichment for FOXM1 and E2F1 binding in the downregulated genes 48 h and 72 h p.i. Combined Z-score is depicted based on the multiplication of log p value computed with Fisher exact test and the z-score which is the deviation from the expected rank by the Fisher exact test. The top 5 ranked gene sets based on combined z-score are depicted, which are all significant (p < 0.05). Gene sets marked with an asterisk indicate significance upon multiple testing correction with ***adj p.val < 0.001 and *adj p.val < 0.05.    Nuclear and cytoplasmic fractionation. The protocol was conducted as described earlier 38       www.nature.com/scientificreports/ G3BP1 (H00010146-M01, Abnova, 1:200 dilution). As secondary antibodies goat anti-rabbit IgG (H + L) Alexa Fluor 488 (A11034, Thermo Fisher, 1:400 dilution) and goat anti-mouse IgG (H + L) CY3 (A10521, Thermo Fisher, 1:400 dilution) were used. For cytoplasmic staining HCS CellMask™ Red stain (H32712, Thermo Fisher) was added according to the manufacturer's instructions. After additional washing, the plates were maintained in 0.1% NaN 3 /PBS and stored at 4 °C for microscopy.
Lysosome staining. Live cells were stained with 5 μg/ml of the nuclear dye Hoechst 33342 (# H3570, Invitrogen) in medium at 37 °C for 5 min. After the plate was washed twice with PBS, 100 μl/well 50 nM Lysotrackerred DND-99 dye (#L7528, Invitrogen) was added in medium at 37 °C for 60 min. Thereafter, medium was replenished and prepared for imaging.
EdU incorporation assay. EdU (5-ethynyl-2′-deoxyuridine), supplied within the Click-iT EdU Alexa Fluor 647 Imaging Kit (#C10340, Thermo-Fisher Scientific, Waltham, MA, USA), was diluted in DMSO to a final concentration of 10 mM and kept at − 20 °C. EdU was added to parallel cultures growing exponentially at 37 °C in 96 well plates to final concentration of 20 μM for 30 min to label the S-phase cells before fixation with 4% paraformaldehyde. The following steps of the Click-iT reaction were performed at room temperature. Additional antibodies were added as previously described in the immunofluorescence protocol. During the secondary antibody incubation step of the immunofluorescence protocol, the Click-iT kit azide and buffer additive were removed from − 20 °C storage and allowed to thaw in a light-protected box. The buffer additive is stored as a 10× solution; once it is thawed, it must be diluted with dH2O to make a working 1× solution. During the third PBS rinse after incubation with the secondary antibody, the reaction cocktail is made, adding the azide last. The reaction cocktail must be used within 15 min of creation. Image processing and data analysis. Image processing was performed in FIJI (http:// fiji. sc), a packaged version of ImageJ freeware (W.S. Rasband, USA. National Institutes of Health, Bethesda, Maryland, USA, http:// rsb. info. nih. gov/ ij/, 1997-2014). Quantification of nuclear and spot signal intensities of (immuno-)stained cell cultures was done using a script specifically written for automated cell-based analysis (CellBlocks.ijm) 39 . In brief, the image analysis pipeline relies on the segmentation of nuclei, cells and intracellular spots and subsequent feature (number, intensity, shape…) extraction for all regions of interest. Downstream data analysis was performed in R Studio and graphs were generated by GraphPad Prism (version 8). To identify numbers (and percentages) of marker-positive and -negative cell nuclei, a manual threshold was set on the intensity values of the marker of interest. For cell cycle staging, we plotted the integrated nuclear signal intensity of the DAPI channel 40 . The results shown are averages (with standard deviation interval, tested by unpaired t-test) of a minimum of three independent experiments. Significance levels were indicated as follows: p < 0.05 (*), p < 0.01 (**), and p < 0.001 (***).
Bulk RNA-sequencing. To  www.nature.com/scientificreports/ cal value (t) and corresponding heatmaps were made using a custom R script. Signature scores were conducted using a rank-scoring algorithm 44 . For the time-series analysis, all genes were clustered according to their average counts across replicates. Counts were quantile normalized across genes and across samples. Clusters were made with k-means method from the Python sklearn package (v0. 22.2). To determine the ideal number of clusters, the elbow method was used 19 .
Single cell RNA-sequencing. Cells were harvested after treatment and washed two times with PBS with 0.04% BSA (A7030; Sigma-Aldrich). The viability and the number of cells were measured with the Countess II (Invitrogen, Carlsbad, CA, USA). Samples were prepared for loading on the 10× Genomics Chromium microfluidics controller, as described in the Chromium Single Cell 3′ protocol (10× Genomics, Pleasanton, CA, USA). Per sample, 5000 cells were loaded to end up with 2500 captured cells. During cDNA amplification, 12 PCR cycles were used. To add an index to the samples, another 12 PCR cycles were used. To check the final quality of the generated libraries, the samples were loaded on a Bioanalyzer High sensitivity chip (Agilent Technologies, Santa Clara, CA, US), and the quantity was checked with qPCR according to the Illumina protocol 'Sequencing Library qPCR Quantification protocol guide' , version February 2011. Samples were pooled equimolarly and sequenced on a Hiseq4000 lane (Illumina, San Diego, CA, USA), generating paired-end reads of 150 base pairs. Afterward, the data were trimmed to paired-end reads of 28 and 91 base pairs. Single cell analysis was performed with the Python scanpy package (v1.5.1). Cells that reported less than 1000 genes were filtered, so were genes that were reported in only 1 or 2 cells. Cells were allowed to have up to 12% mitochondrial reads. For determining the differential expression of single cell genes, the Python package diffxpy (v0.7.4) was used. To calculate the bulk RPE1-MYCN-ER signatures in the single cell uniform manifold approximation and projection (UMAP) spaces, only genes were considered with an adjusted p-value < 0.01. Single cell library sizes were adjusted to 10,000 and transformed by natural logarithm after addition of 1. UMAPs with 3 components were produced by taking 10 neighbours in a 25 component PCA space. Cell phase was determined with the S and G2M gene sets used in 45 . In figures with all treatments together, cell phase was determined across all data, in figures were only one treatment is displayed, cell phase was determined for each cell only within that treatment data. Gene senescence signature scores were calculated as the number of genes that had expression in the cells. For the bulk RPE1-MYCN-ER signatures in the single cell, the average rank of the genes in the set was calculated. Full scripts can be found at https:// github. com/ dicaso/ apety pe/ tree/ master/ examp les/ RPEMY CN.
Quantification and statistical analysis. Statistical analyses were performed in R Studio and statistical significance of differences between conditions for the immunofluorescence and caspase-glo experiments were determined by unpaired Student's t-test using GraphPad Prism (version 8). The ANOVA (analysis of variance) test followed by a post-hoc Tukey's test for multiple comparisons was used to determine differences in gene expression and signatures scores between 4 different groups with 3 (RT-qPCR) or 4 (RNA-Seq) biological replicates per condition. For RT-qPCR experiments, reference genes were excluded if GeNorm M value was greater than five and/or Coefficient of Variation greater than two, according to the qBase + software. All error bars represent the SD or 95% CI of the three or four biological replicates. For statistical testing of the focus formation (NIH-3T3) assay, cumulative growth values were modelled using a linear model with time and treatment as fixed factors and an interaction term for both (Python package statsmodels). This model explained 99% of the observed variance and performed better than a restricted model omitting the interaction term for time and treatment (LR test P-value < 0.01, interaction term P-value < 0.01 in the more complex model). The details of quantification and statistical methods used are mentioned in the figure legends. www.nature.com/scientificreports/