A multi-omics integrative approach unravels novel genes and pathways associated with senescence escape after targeted therapy in NRAS mutant melanoma

Therapy Induced Senescence (TIS) leads to sustained growth arrest of cancer cells. The associated cytostasis has been shown to be reversible and cells escaping senescence further enhance the aggressiveness of cancers. Chemicals specifically targeting senescent cells, so-called senolytics, constitute a promising avenue for improved cancer treatment in combination with targeted therapies. Understanding how cancer cells evade senescence is needed to optimise the clinical benefits of this therapeutic approach. Here we characterised the response of three different NRAS mutant melanoma cell lines to a combination of CDK4/6 and MEK inhibitors over 33 days. Transcriptomic data show that all cell lines trigger a senescence programme coupled with strong induction of interferons. Kinome profiling revealed the activation of Receptor Tyrosine Kinases (RTKs) and enriched downstream signaling of neurotrophin, ErbB and insulin pathways. Characterisation of the miRNA interactome associates miR-211-5p with resistant phenotypes. Finally, iCell-based integration of bulk and single-cell RNA-seq data identifies biological processes perturbed during senescence and predicts 90 new genes involved in its escape. Overall, our data associate insulin signaling with persistence of a senescent phenotype and suggest a new role for interferon gamma in senescence escape through the induction of EMT and the activation of ERK5 signaling.


INTRODUCTION
Hayflick first observed in 1961 that fibroblasts stop proliferating after 50 passages and display an enlarged and flattened phenotype characteristic of senescence [1,2]. This observation was later explained by the progressive shortening of telomeres upon replication, which ultimately triggers a DNA Damage Response (DDR) leading to cell cycle arrest [3,4]. This process was called "replicative senescence". Several cellular events and environmental conditions induce senescence, and it now emerges as a generic stress response involved in central biological processes such as embryonic development, wound healing and aging [2,4,5].
Many stressors that induce senescence (e.g., oxidative stress, radiation), induce DNA damage and activate DDR. This leads to engagement of the Senescence Associated Secretory Phenotype (SASP), one of the characteristic features of senescence in which cells produce and secrete cytokines and growth factors. The SASP itself can trigger senescence in adjacent cells and contributes to the emergence of an inflammatory environment as the senescent cells, which are not eliminated by the immune system, accumulate in the organism [5]. This accumulation leads to the development of degenerative and hyperplastic pathologies and carcinogenesis, linking senescence to aging [6,7]. Research in this field is challenging due to the heterogeneous and dynamic nature of the SASP and the lack of a universal marker for senescence.
Focusing on cancer, the activation of known oncogenes results in a constant proliferative signaling which causes replication forks to trigger DNA damage and drive primary cells to undergo Oncogene-Induced Senescence (OIS) [8]. Therefore, senescence may be considered as an early evolutionary mechanism of protection against cancer leading to cell cycle arrest of overproliferating cells [9]. However, cells can escape this arrest and reestablish growth, ultimately leading to the development of cancer. Cancer cells can also undergo senescence when exposed to therapies, so called Therapy-Induced Senescence (TIS).
Thus, senescence appears as a natural barrier against cancer, but recent works suggest that senescence may be a reversible process [10]. Exit from senescence has been recently linked with the stem-like properties of cancer cells, reinforcing the idea that aside from the SASP, senescence itself may be a transient state in carcinogenesis.
Melanoma arises from the deregulated proliferation of melanocytes representing the deadliest type of skin cancer accounting for about 75% of related deaths. Standard of care involves targeted therapies and immunotherapies [11]. For the BRAF mutated subtype, representing half of the cases, specific first line inhibitors targeting both MEK and BRAF have been developed. The NRAS mutated subtype, which accounts for a quarter of all melanoma cases, lacks such a targeted approach but ongoing clinical trials are assessing the effects of combined MEK and CDK4/ 6 inhibitors (NCT01781572; NCT02065063). However, tolerance or resistance to such targeted treatments inevitably develops.
MicroRNAs have been shown to be involved in resistance to treatment. miRNAs are 20-22 nucleotides short RNAs, which in complex with the Argonaute (AGO) protein, act as posttranscriptional regulators of gene expression [12]. Canonically, miRNAs destabilise mRNAs by binding to a complementary sequence within the 3' UTR, this resulting in downregulation of the mRNA and the encoded protein [13]. Experimental methods based on immune-precipitation and sequencing have been developed to profile miRNA-mRNA interactions. The qCLASH method has the advantage of capturing direct physical interactions only, through an additional intermolecular ligation step linking mRNAs to their bound miRNAs [14].
To characterise the response of 3 NRAS mutant melanoma cell lines to a combination of MEK and CDK4/6 inhibitors, we treated the cells over 33 days and collected samples at multiple time points. Samples were analyzed by total RNA-seq, small RNA-seq, qCLASH and kinome profiling. Our results recapitulate several previous observations while displaying an unforeseen interplay between interferon signalling and senescence. Additionally, we uncover deregulated processes associated with the onset of senescence and its exit by adapting a non-negative matrix trifactorization (NMTF)-based methodology, "iCell" [15], to integrate bulk RNA-seq and single-cell data. Moreover, the dimensionality reduction and the clustering properties of the NMTF allowed us to associate novel genes to senescence escape.

Experimental design and description of the data sets
To study the effects of MEKi and CDK4/6i, 3 NRAS mutant melanoma cell-lines: MELJUSO, SKMEL30 and IPC298 were selected as they showed distinct cellular responses with spindleshape (MELJUSO), augmented pigmentation (SKMEL30) and typical melanoma morphology (IPC298) phenotypes (Fig. 1). To further normalize the effect of treatment across cell lines, we determined the GI50 (drug concentration at which cell growth is reduced by half) and opted for a concentration of 110nM, 35nM, and 16nM of Binimetinib (MEKi) for MELJUSO, SKMEL30 and IPC298, respectively and 1 μM for Palbociclib (CDK4/6i) as this inhibitor led to incomplete dose-response curves (Supplementary Table 1). Upon treatment over 33 days, the 3 cell lines showed different behaviour with IPC298 suffering very little from the treatment showing constant proliferation and no morphological change. SKMEL30 on the other hand stopped proliferating and accumulated pigmentation before restoring proliferation after around 14 days. Finally, MELJUSO stopped dividing and started to stretch to acquire a spindle-shape phenotype after 4 days.
As distinct morphological changes were observed at these time points, we profiled the transcriptome of the cell lines at day 0, 4, 14 and 33 using total RNA-seq. Cellular signaling supporting the adaptation to the treatment was investigated by profiling the kinome of the senescent MELJUSO and adaptative SKMEL30 cell lines using PamGene technology. Next, we characterised the resistant miRNA interactome of the SKMEL30 and IPC298 at 33 days by combining small RNA-seq and the qCLASH method. Finally, by integrating bulk and single cell RNA-seq data, deregulated pathways and predicted genes associated with senescence escape were identified.
Cell lines display coherent transcriptomic profiles composed of interferons, EMT and senescent responses under prolonged CDK4/6i and MEKi treatment Distinct cellular responses to treatment were reflected at the transcriptomic level where the cell lines clustered separately on the PCA (Fig. 2A). As MELJUSO displayed a different evolution over time, we decided to analyze each cell line individually. Differential expression analysis, comparing later time points to day 0, revealed a greater change in the MELJUSO transcriptome. MELJUSO corresponds to 9548, 11362 and 10126 significant (p.adj <= 0.05) differentially expressed genes at day 4, day 14 and day 33, respectively when compared to SKMEL30 and IPC298 (9208, 9749, 6509 and 5436, 7575 and 3217, respectively) ( Supplementary  Fig. 1). Interestingly, across all cell lines and for all time points, we noticed the deregulation of interferon-related genes such as STAT1, IRF7, MX1, OAS2 and IFI44.
Gene Set Enrichment Analysis (GSEA) was performed on differentially expressed genes using the "hallmarks" from MSigDB, which revealed the downregulation of pathways related to proliferation such as "MYC TARGETS V1", "E2F TARGETS", and "G2M CHECKPOINT" (Fig. 2B), all affected by our inhibitors. The clustering of samples on top of the heatmap clearly separates proliferative samples (left part) comprising the IPC298 cell line and SKMEL30 day 33 from cytostatic samples (right part) with MELJUSO and early time points for SKMEL30. Interestingly, we noticed across all cell lines an enrichment in interferon responses and an enrichment in Epithelialto-Mesenchymal Transition (EMT) for the cytostatic samples (Fig. 2B). CDK4/6 inhibitors have been shown to induce cellular senescence in different cell types [16,17]. We therefore used a recently published classifier, "SENESCopedia", to predict the levels of senescence in our samples (Fig. 2C) [18]. On day 14, all cell lines display high senescence scores suggesting that they all engage in a related transcriptional programme. Of note, all samples with a high senescence score are also enriched in EMT which is coherent with previous observations of EMT in different epithelial cell types after induction of senescence (fibroblasts, colorectal and lung cancer) [19][20][21].
Altogether, these results suggest that all 3 cell lines trigger the same type of responses composed of an interferon response, a senescence and "EMT" programme but with a different intensity and temporality, which may explain the different phenotypes observed.
Kinome profiling of adapting cell lines shows reactivation of RTK pathways and suggests a key role of ERK5 To better understand how the cell lines adapt to the treatment, we acquired PamGene data ( Supplementary Fig. 2) to profile the kinome of MELJUSO and SKMEL30 cells at day 0, 1, 4 and 33 as these two cell lines showed persistent cytostasis and adaptation to treatment. MELJUSO stops proliferating and changes morphology after about 4 days while SKMEL30 acquires a dark pigmentation before restoring growth. Kinase activities were inferred from the levels of peptide phosphorylation at different time points and compared to day 0 (summarised in Fig. 3A). A phylogenetic representation shows that these two cell lines react similarly to the treatment only at early time points (Supplementary Fig. 3).
To perform pathway enrichment on a reduced set of proteins such as kinases, we used a network-based approach and obtained coherent results for the two cell lines (Table 1). Indeed, pathways like RIG-I or Toll-like receptor signaling (related to innate immunity and interferon responses) were enriched at early time points for both SKMEL30 and MELJUSO. Among the enrichment results, RTK pathways involving ErbB, neurotrophin and insulin appear (Fig. 3B, C). Some pathways were enriched before the increase in receptor Fig. 1 Experimental setup, cellular phenotypes and omics characterization. IPC298, MELJUSO and SKMEL30 cell lines display different phenotypes upon CDK4/6i and MEKi. We characterised those using RNA-seq, kinome profiling and qCLASH method. RNA-seq and qCLASH data were further combined to construct a resistant miRNA network. Bulk and single-cell RNA-seq were finally integrated into "iCell" networks. activity reached the significance threshold and those were therefore not represented on the phylogenetic tree (Supplementary Fig. 3). This signaling seems to persist after the receptor resumes normal activity (Table 1). It has previously been shown that the inhibition of ERK1/2 suppresses negative feedback on RTK expression, which can then be further activated to compensate for this inhibition [22] (Fig. 3B, Supplementary Fig. 3). Coherent with the enriched pathways, we observed in both cell lines a significant increase in the activity of RTKs such as HER3 and NTRK2 as well as insulin-related receptors (INSR, IGF1R and IRR). Noteworthy, downstream insulin signaling was enriched in the cytostatic MELJUSO only (detailed representation of pathways Supplementary Figs. [4][5][6][7][8]. Other RTKs previously related to resistance such ALK, AXL or c-MET displayed an increased activity (complete list Supplementary Fig. 9) but those were not associated with enriched downstream signaling. Within the MAPK pathway, the decrease in ERK5 (MAPK7) activity was partially relieved in SKMEL30 but not in MELJUSO (Fig. 3B). Along these lines, Tubita et al. recently showed that knock-down by shRNA or inhibition of the kinase activity of ERK5 triggers senescence in melanoma [23].
Overall, the kinome data show the reactivation of ErbB, neurotrophin and insulin pathways and enrichment of their downstream signaling. Interestingly, ERK5 activity was strongly reduced in the cytostatic cell line MELJUSO (Fig. 3B).

Contribution of miRNAs to the resistant phenotype
To see if miRNAs could participate in the adaptation to treatment and the proliferative phenotypes, we profiled the miRNome of all three cell lines by small RNA-seq at day 0 and 33. In line with their phenotypes, IPC298 showed little deregulation of miRNAs compared to SKMEL30 and MELJUSO cell lines ( Supplementary   Fig. 10). Next, we performed qCLASH analysis for the two proliferative cell lines IPC298 (early adaptation) and SKMEL30 (late adaptation) (summarised in Fig. 4A).
Although qCLASH-based miRNA interactomes do not provide a complete picture of all interactions, PCA was able to discriminate between the different cell lines and conditions ( Supplementary  Fig. 11). Previous studies reported a correlation between the number of detected hybrids and the expression levels of its miRNA and mRNA components [14,24]. These correlations hold true for our matching total and the small RNA-seq data when taking the total count across the replicates (Fig. 4B, left panels). Moreover, further summing up the hybrids for each miRNA or mRNA strengthens the correlation between RNA-seq and qCLASH data, especially for miRNAs (Fig. 4B, right panels). qCLASH data also contain transient interactions, which are unlikely to be functionally relevant for the cell and technical replicates are often used to further select most likely interactions. We explored how the detection of interactions in one, two or three replicates depends on their expression level and observed that interactions detected in triplicates were associated with a greater expression supporting the rationale of using such thresholds ( Supplementary Fig. 12). To further select the most relevant miRNAs, we constructed a network with the interactions specific to the resistant phenotype. We first selected the interactions detected in triplicates, then ensured that they were present in both resistant IPC298 and SKMEL30 but absent from their untreated counterparts (Fig. 4C). Finally, we overlaid the log fold changes from the small and total RNA-seq for the corresponding miRNAs and mRNAs onto the resulting network (Fig. 4D). Based on those representations, we selected some interactions to confirm their deregulation by qPCR. We investigated if those observations could be extended to two BRAF mutant C Senescence scores predicted by Cancer SENESCopedia. SENESCopedia webtool (https://ccb.nki.nl/publications/cancer-senescence) uses a gene expression classifier to predict senescence in cancer cell samples [18].
cell lines, A375 and 624Mel which underwent the same treatment and re-established growth ( Supplementary Fig. 13). Most of the selected interactions did not exert a robust inverse expression between miRNA and mRNA ( Supplementary Figs. 14, 15). Among mRNAs, we observed a robust upregulation of CCND1 which is a known resistance mechanism to CDK4/6i as well as a strong upregulation of IFN related gene IFI6 ( Supplementary Fig. 13). Except for the amelanotic cell line A375, we observed the upregulation of LGALS3BP, TXNIP and TYRP1 (Fig. 4E). Among miRNAs, miR-211-5p was upregulated in all proliferative cell lines except the amelanotic A375, which is known not to express it (Fig. 4F) [25], and MELJUSO (Fig. 4G), which has very low levels.   Syk   TAK1   TAO1  TAO2  TAO3   TBCK   TBK1   TEC   TESK1  TESK2   TGFβR1   TGFβR2   TIE1   TIE2   TLK1   TLK2   TNIK/ZC2   Tnk1   Trad   Trb1   Trb2   Trb3   Trio   TRKA   TRKB   TRKC   TSSK1  TSSK2   TSSK3   TSSK4   TTBK1   TTBK2   TTK   TTN  Overall, the number of interactions detected by qCLASH method is influenced by the expression level of the corresponding miRNA and mRNA. Most of the detected interactions did not exert the expected inverse expression pattern between miRNA and mRNA. miRNAs and mRNAs displayed a more robust change in expression across cell lines. Except for the amelanotic A375, upregulation of miR-211-5p is associated with resistant and proliferative cell lines.
iCell integration of bulk and single-cell RNA-seq reveals perturbed pathways and predicts genes facilitating senescence escape To gain further insights into NRAS mutant melanoma cells' adaptation to treatment, we integrated the bulk RNA-sequencing data with previous condition-matching single-cell data from our lab by adapting and applying a Non-negative Matrix Tri-Factorization (NMTF) approach called iCell [15]. We used the iCell methodology because it is a versatile data fusion framework, which has previously been used to perform integrative comparative analyses of disease and control tissue data, predicting new cancer and COVID-19 related genes [15,26] and suggesting potential drug re-purposing candidates [26].
For each condition (treated cell line and time point), we constructed Protein-Protein Interaction (PPI) networks by overlaying the experimentally validated PPI network for Homo Sapiens from BioGRID [27] with the condition-specific bulk RNA-sequencing data. Also, we derived a condition-specific co-expression (COEXP) network from the corresponding bulk RNA-seq data (see Methods section). Next, we applied Non-negative Matrix Tri-Factorization (NMTF) on the corresponding adjacency matrices of the resulting PPI, COEXP and the single-cell expression data collectively, to ensure data fusion (Methods and [15]). This resulted in the creation of nine conditionspecific iCell networks, by multiplying the G1 matrix resulting from the data fusion and its transpose (which is equivalent to computing the dot product between gene pairs) and considering as interactions the 1% highest values in each row and column (Methods and [15]). As illustrated in Supplementary Fig. 16, resulting conditionspecific iCell networks were investigated by first comparing the global topology of all resulting iCell networks. Then, these topological changes were further characterised by identifying the genes that change or conserve their local topology between each pair of conditions. We related the identified most rewired genes and the least rewired ones to MSigDB hallmarks by performing an Over-Representation Analysis (ORA) (see Methods for details). Finally, for each condition we clustered the genes based on their local topology in the iCell network and then found the MSigDB hallmarks enriched in the resulting clusters. We focused on two senescence-escape related hallmark pathways (IFN gamma and EMT) to predict novel genes potentially related based on their appearance in the clusters enriched across all conditions.
We compared the overall topology of iCell networks obtained for each condition by using thus far the most sensitive measure of global network topology, the "Graphlet Correlation Distance" (GCD-11) (see Methods and [28]). The lower the GCD-11 value, the closer the topologies of the networks [28]. We used Ward's hierarchical clustering to cluster the resulting heatmap of the GCD-11s between all iCell networks (all-to-all matrix). In line with our phenotypic observations, we showed on the heatmap (Fig. 5A) that for the resistant cell line IPC298, the conditions D4 and D33 have the lowest GCD-11 value (0.26) and hence clustered together, suggesting early adaptation to treatment of this cell line. For the two other cell lines, conditions D0 and D4 grouped together (GCD-11 values of 0.71 for SKMEL30 and 0.87 for MELJUSO), suggesting that these two cell lines adapt later.
Next, the local topology of genes in each condition-specific iCell networks was analysed to find out if some are rewired between conditions. To this end, we focused on genes that were expressed in both conditions (Venn diagrams in Fig. 5B and Supplementary  Fig. 17). Within each iCell network, the local topology around a gene in the network was quantified by using its "Graphlet Degree Vector" (GDV). We then measure the rewiring of a gene between two conditions by comparing its GDVs in the two networks using "Graphlet Degree Vector Similarity" (GDVS) (detailed in Methods and in [29]). Between two conditions, genes with high GDVSs have conserved topologies and are called "stable", while the genes with low GDVSs are rewired and are called "perturbed". We related the top 10% of the most "stable" and the top 10% of the most "perturbed" genes to MSigDB hallmarks [30] by performing ORA. For each cell line, the condition-specific iCell networks were compared between time-points (i.e., D4 vs D0, D33 vs D4, and D33 vs D0). For all cell lines, we observed that the top 10% of the most "perturbed" genes were enriched at different contrasts in "ADIPOGENESIS" (also observed in senescent fibroblasts [31]) and "OXIDATIVE PHOSPHORYLATION" (OXPHOS) (a pathway commonly dysregulated in senescence which has been proposed as a targeting strategy [32]) (Fig. 5C). Then, condition-specific iCell networks were compared between cell lines at matching time points (e.g., MELJUSO_D0 vs SKMEL30_D0). Interestingly, the top 10% of the most "stable" genes were only found to be enriched when comparing IPC298 vs MELJUSO at day 0 and when comparing MELJUSO vs SKMEL30 at day 4, which suggests that these cell lines react coherently early on and diverge later ( Supplementary Fig. 17). At day 0, and across all pairwise comparisons involving IPC298 cell line, the top 10% of the most "perturbed" genes were enriched for OXPHOS, suggesting a different initial metabolic state for IPC298. On day 4, when comparing IPC298 and SKMEL30 cells, the top 10% of the most "perturbed" genes were enriched in "INTERFERON ALPHA/GAMMA RESPONSE". At days 4 and 33, but not across all pairwise comparisons, we found that the top 10% of the most "perturbed" genes are enriched in "EPITHELIAL MESENCHYMAL TRANSITION", "ADIPOGENESIS" and "COAGULATION".
Further information on the relevance of these enriched pathways can be found in the Discussion.
Finally, to predict novel genes which could contribute to senescence escape facilitated by interferons, genes were clustered based on their GDVs signatures in each condition-specific iCell network using "Ward's hierarchical clustering" and the "kmedoids" algorithms. For each condition, we performed ORA to associate the clusters to hallmarks from MSigDB (results of ORA on the clusters obtained by using k-medoids clustering in Fig. 5D). The results for the two algorithms showed limited agreement (Fig.  5E, left), as measured by Adjusted Rand Index for each condition: 0.31, 0.38 and 0.30 for MELJUSO at day 0, day 4 and day 33; 0.32, 0.30 and 0.26 for SKMEL30 at the same time points and 0.27, 0.31 and 0.45 for IPC298. Therefore, we considered the results from the two different clustering methods (hierarchical and k-medoids) separately. Enriched clusters were filtered to exclude genes that are already known to be associated with hallmarks in MSigDB. For the remaining genes, we counted the number of times a gene was found in the clusters enriched in one of the two senescence escape-related pathways, EMT and interferon gamma (Fig. 5E  right). Then, based on the distributions of the number of times genes appeared in these enriched clusters, we predicted genes that appeared in two or more clusters as related to senescence escape. To obtain more robust predictions, we only considered genes predicted by both algorithms. This resulted in a set of 90 predicted genes, among which CCND1, a known regulator of G1 phase. Furthermore, many of our predictions are known to interact with TP53, further supporting their connection to senescence (Fig.  5F, Supplementary Table 2).
To conclude, the integrated analysis of bulk and single-cell data by using the iCell methodology further confirms our phenotypic observations and uncovers the most dysregulated processes ("OXIDATIVE PHOSPHORYLATION", "ADIPOGENESIS", "COAGULA-TION" and "EPITHELIAL MESENCHYMAL TRANSITION") during adaptation to the treatment. Furthermore, it allowed us to predict 90 new genes involved in senescence escape.

DISCUSSION
Most cancer cells inevitably develop resistance to targeted therapies. Anticipating the development of such resistance to a newly tested combination of CDK4/6i and MEKi, we characterised the transcriptomic responses in 3 NRAS mutant melanoma cell lines over an extended treatment duration of 33 days. Our RNA-seq data exhibited a strong enrichment in interferon responses, coherent with previous observations by Goel et al. [33] following CDK4/6 inhibition. This study showed that CDK4/6 inhibitors downregulate E2F2 and DNMT1, subsequently leading to a hypomethylation of the genome (Supplementary Fig. 18). The works of Chiappinelli et al. [34] and Roulois et al. [35] first showed that hypomethylation brought about by DNMT inhibitors induce the expression of otherwise silenced endogenous retroviruses (ERVs). These ERVs then form dsRNAs in the cytoplasm and trigger an anti-viral immune response associated with interferon production. Alternatively, ERVs can be reverse transcribed and produce ssDNA and dsDNA triggering the cGAS-STING pathway leading to similar responses [36]. The Fig. 4 miRNA-mRNA interactions potentially contributing to resistance. A Summary of the interactions detected through the qCLASH method. B Top part, correlations between hybrid counts (qCLASH) and miRNA/mRNA counts (RNA-seq). Bottom part, summing hybrid counts per mRNA or miRNA (qCLASH) further strengthen the correlations with the RNA-seq. C To select candidates, we considered the interactions detected in three technical replicates (colored circles). Then for the two cell lines, interactions present in the resistant state but absent in the sensitive state were considered. Finally, interactions common to the two cell lines IPC298 and SKMEL30 were used. D RNA-seq data for SKMEL30 cell line as log fold changes overlayed onto the resistant miRNA network. Selected interactions for further validation by qPCR are represented with thick lines and the corresponding mRNAs and miRNAs with bolded borders. E Relative expression of LGALS3BP, TXNIP and TYRP1 mRNAs assessed by qPCR across three NRAS and two BRAF mutant melanoma cell lines. F Relative expression of miR-211-5p across the same cell lines as determined by qPCR. G Expression of miR-211-5p in the small RNA-seq.  activation of these two pathways can be concomitant and their relative importance may be cell-type specific. This phenomenon, termed "viral mimicry", has been extended to other types of retroelements and it has been proposed to be integral part of carcinogenesis in a rheostat model [37]. It has been reported in different cancer types and in response to different epigenetic perturbations (DNMTi, HDACi, HMTi). Brägelmann et al. [38] also reported the activation of viral mimicry in the case of BRAF mutant melanoma and MAPK inhibitors (BRAFi + MEKi) further supporting the E2F2-DNMT1 axis as a mechanism of action. Melanoma cells are considered very immunogenic [39] and this partially explains why melanoma respond better to immunotherapy than other cancer types. IFN signaling is known to upregulate PD-L1, explaining why many treatments activating viral mimicry have been shown to synergise with immune checkpoint inhibitors. Here, we confirm the appearance of an IFN signature following CDK4/6 and MEK inhibition suggesting a potential induction of viral mimicry.
A "senescence score" predicted from our data with SENESCopedia [18] suggested induction of a senescent transcriptional programme in all cell lines. Cytokines composing the SASP can induce features of EMT [40] and EMT has been recently proposed to inhibit senescence in HeLa cells treated with Doxorubicin [41]. TIS and EMT share common transcription factors with opposite effects on both processes ( Supplementary Fig. 19) [42] and cancer cells escaping from TIS have been associated with a mesenchymal phenotype [43]. IFN gamma, by promoting EMT [44][45][46], could further support escape from senescence ( Supplementary Fig. 20A). To understand the relationship between IFN gamma and EMT gene sets (from the hallmarks) and a recent gene set characterizing senescence in cancer cell lines [47], we looked at changes in expression (log2 fold change) between those gene sets at day 14 where all cell lines display high senescence scores. IFN gamma shows an opposite trend in expression to the EMT and senescence gene sets ( Supplementary  Fig. 20B), suggesting that IFN may accelerate the EMT process and thereby repress senescence. Albeit interferon signaling has earlier on been linked with senescence [48], only later Katlinskaya and Yu et al. showed a role of IFNAR1 in the induction of senescence [49]. Recently, the accumulation of double stranded RNA was proposed to also contribute to the onset of senescence [50,51]. These findings suggest a link between senescence and viral mimicry, which can engage autocrine and paracrine secretion of interferons. Embryonic stem cells (ESCs) are another cellular model where transposable elements play an important role, where they are profoundly repressed during embryogenesis and leakage in their expression may impair developmental processes. Recent work by Asimi et al. in mesenchymal ESCs suggests that ERV transcription participates in the remodeling of epigenetic landscapes by competing with superenhancers for the formation of transcriptional condensates, which may explain the phenomenon of onco-exaptation (oncogene activation by a TE-derived promoter) [52]. Also recently, Liu et al. demonstrated a role of ERV expression in the induction of senescence and attributed the observed paracrine effects to the production of retrovirus like particles [53]. We propose here that the IFN gamma signalling programme activated after induction of viral mimicry may contribute to senescence escape by inducing EMT. This might represent a general cellular mechanism to cope with different extrinsic and intrinsic stressors.
To identify targetable proteins, which could explain the senescent phenotype, we profiled the kinase activities of SKMEL30 and MELJUSO cells upon CDK4/6 and MEK inhibition. We observed the reactivation of different RTKs due to ERK1/2 inhibition as previously reported by others [22,54]. Some of the corresponding pathways such as ErbB, neurotrophin and insulin were found enriched through network-based enrichment hinting at an active downstream signalling. Interestingly, all those pathways have been shown to activate ERK5. In breast cancer, ERK5 contributes to the response to EGF stimulation [55] and is involved in the resistance to ErbB-2 (HER2) inhibitors by facilitating G1-S transition [56]. In BRAF mutant melanoma resistant to BRAFi and MEKi, IGF1R activation induces ERK5 phosphorylation [57]. In senescent fibroblasts, BDNF-TrkB has been shown to activate ERK5 and was proposed to form an autocrine loop further enhancing senescent cell viability [58]. ERK5 activation has also been linked to PDGFR signalling in NRAS mutant melanoma treated with MEK and ERK1/2 inhibitors [59,60] (Fig. 3B, D, E). In other cell types, Src is also known to activate ERK5 and Focal Adhesion Kinases (FAK1 and FAK2) [61,62]. Other Src Family Kinases (SFKs) display an increased activity (Supplementary Fig. 21) and have been previously linked with EMT [63]. Finally, IFN gamma was reported to activate ERK5 and its activation is necessary for the full expression of interferon stimulated genes [64]. Linking this observation with our RNA-seq data ( Supplementary Fig. 22) further supports the idea that IFN gamma could facilitate senescence escape by activating ERK5.
Investigating the contribution of miRNAs to the establishment of resistance, we noticed the upregulation of mir-211-5p in all proliferative cell lines but the amelanotic A375. miR-211-5p is highly expressed in the melanocytic lineage as well as in the nervous tissue, as reported by the TissueAtlas2 [65]. Both cell types derive from the neurectoderm. Interestingly, miR-211-5p is known to derive from the TRPM1 transcript and both are very lowly expressed in the senescent cell line MELJUSO (Supplementary Fig.  23). miR-211-5p inhibits DUSP6, which regulates ERK5 kinase activity, and this axis has been shown to be relevant in vivo in genetic variants of A375 overexpressing mir-211 or DUSP6 and treated with BRAF/MEKi [66]. In qPCR validation experiments, mRNAs displayed a more robust change in expression with LGALS3BP, TXNIP and TYRP1 being upregulated in four out of five cell lines.
LGALS3BP and TXNIP are interferon-inducible and interestingly, LGALS3BP has been reported to promote EMT.
Integrating single-cell and bulk RNA-seq data, we observed the perturbation of OXPHOS in all cell lines. This is in line with previous works demonstrating metabolic adaptation after MAPK inhibition or CDK4/6i in BRAF mutant and uveal melanoma [67,68]. Moreover, OXPHOS metabolism has been previously reported as a generic drug resistance mechanism and as a key decision point between cell death and survival [69,70]. Upon comparison of different treatment time points, differences in processes such as "ADIPOGENESIS" and "COAGULATION" were observed. Interestingly, those have recently been associated with replicative senescence [31], cell cycle arrest caused by CDK4/6i [71] and escape from OIS [72]. EMT has been reported to inhibit senescence [19,41]. We therefore clustered the genes in the networks and enriched those clusters to report novel associations between genes and the processes of EMT and IFN gamma response, which could support senescence escape. Our Fig. 5 "iCell" integration of bulk and single-cell RNA-seq. A Comparison of iCell networks topology using Graphlet Correlation Distances (GCD). iCELL networks recapitulate phenotypic observations and early adaptation of IPC298. B Venn diagrams representing genes overlap between conditions. C Comparison of gene local topology using Graphlet Degree Signature Similarity (GDSS). Within all cell lines, the different time points were compared. ORA was performed on the top 10% most "stable" or "perturbated" genes. D Genes were clustered according to their Graphlets Degree Signatures. ORA was performed to associate gene clusters with biological function. E Left panel: number of enriched clusters detected for each gene across all conditions. Right panel, Number of detected clusters for the EMT and iFN gamma pathways (related to senescence escape) for the two algorithms. F Protein-protein interactions (BioGRID) between P53 and the genes associated with EMT and IFN gamma.
predicted genes are coherent with the results of a genetic screen aiming to identify genes contributing to senescence bypass after p16 overexpression in melanoma [73]. Interestingly, the IFNinducible gene LGALS3BP, previously reported to contribute to EMT and identified in our miRNA resistant network, was predicted by our approach to be associated with senescence escape. On the other hand, TXNIP, another IFN-inducible gene from the miRNA network, recently shown to participate in the induction of senescence upon IGF1 stimulation, was found associated only once and therefore not selected [74]. This further suggests that our approach can relate genes to function beyond their expression pattern, providing a list of 90 potential targets to overcome senescence escape in cancer cells (Fig. 5.F). For example, our topprediction NEAT1 was shown to be important for the viability of lung cancer cells [75] while its knock-down triggered interferon signaling and impaired tumor growth in vivo. Finally, NEAT1 was shown to physically interact with DNMT1 to epigenetically inhibit TP53, STING and cGAS. Among the predictions with the lowest scores, MTDH has been shown to epigenetically regulate TWIST1 and promote stemness in breast cancer [76] while AKAP12 has been suggested to protect ERK5 from PKCζ phosphorylation on residue S486A which inhibits its transcriptional activity [77].
Taken together, our data describe how NRAS mutant melanoma adapt to CDK4/6 and MEK inhibition by triggering an EMT programme. This is reminiscent of the Neural Crest Stem Cell (NCSC)-like state identified in single cell RNA-seq of patient samples after MAPKi by Rambow and colleagues [78]. This state was later characterised by FAK-dependent AKT activation [79], associated with slow or no proliferation and the expression of EMT-related transcription factors.
In sum, our work sheds new light on the intimate relationships between three interacting pathways, senescence, interferon and insulin signalling, with the EMT process (Fig. 6). Our data suggest a new role of interferon gamma, which, by triggering EMT, can facilitate senescence escape while insulin signaling is associated with a robust cytostasis and senescent transcriptome.

Total RNAseq
Transcriptome sequencing was performed using total RNA. Libraries were prepared using "TruSeq Stranded Total RNA Library Prep Gold" Kit (Illumina, USA) and sequenced as 2x75bp on a NextSeq 550 (Illumina, USA). The sequencing was carried out at Luxgen platform, LNS. Data were aligned to the reference genome (hg38) using STAR (v 2.7.9a) and differential expression analysis was performed using R (v 4.2.2) and DESeq2 (v 1.34.0). This analysis was automated in a snakemake workflow (v 0.3.2) [80] (https:// gitlab.lcsb.uni.lu/aurelien.ginolhac/snakemake-rna-seq) using docker file (v0.6) (https://hub.docker.com/r/ginolhac/snake-rna-seq/tags). Mixing all cell lines into one analysis masked the differences of each time course. Therefore, a subset per cell line was created and the differential gene expression was performed independently. To predict senescence scores with SENESCopedia [18], data were reprocessed and aligned to GENCODE (v 34) using kallisto (v 0.46.2) according to the recommendations of the authors.

Gene-list based pathway enrichment
Over-Representation Analysis (ORA) and Gene Set Enrichment Analysis (GSEA) were performed using R (v 4.2.2) and the package "clusterProfiler" (v 4.0.5). The Hallmark Gene Set Collection from MSigDB were used as reference. Differentially expressed genes used for GSEA were ranked according to the π value, which takes into consideration both the statistical significance and the biological relevance through the product of the p-value and the fold change [81]. The results were represented using the package "pheatmap" (v 1.0.12) using the default "Ward.D2" clustering algorithm. ORA results were displayed using the "ggplot2" (v 3.3.6) package.

Kinase activity (PamGene)
To assess the kinase activity within the NRAS mutant melanoma cell lines MELJUSO and SKMEL30 to a combination of CDK4/6 and MEK inhibitors, cell samples were collected at day 0, 1, 4 and 33. Cell pellets were lysed in Mper (78501, Thermo Fisher) supplemented with Halt Phosphatase Inhibitor Cocktail (78428, Thermo Fisher) and Halt Protease Inhibitor Cocktail, EDTA free (78437, Thermo Fisher), for 15 min on ice. The lysate was clarified (10000g, 15 min, 4°C) before aliquoting and freezing at −80°C. Adopting a balanced design, three replicates for each condition and cell line were Fig. 6 Graphical summary. The phenotypes observed across three different NRAS mutant melanoma cell lines can be summarised as "senescent" (left) and "resistant" (right). Insulin signalling is associated with the "senescent" phenotype and can be explained by the mechanism of action of MEKi, which downregulates ERK1/2, suppressing negative feedback on RTKs. Interferon signaling can be explained by CDK4/6i, which downregulates DNMT1 leading and lead to hypomethylation of the genome and expression of Endogenous Retroviruses (ERVs). Induced expression of intracellular ERVs trigger an innate immune response and production of interferons. IFN γ can activate ERK5, inducing EMT. This phenotype is moreover associated with upregulation of miR-211-5p, which indirectly regulates ERK5. assessed for both tyrosine kinase activity on PTK chips (86402, PamGene International B.V., 's-Hertogenbosch, the Netherlands) and serine threonine kinase activity on STK chips (87102, PamGene International B.V., 's-Hertogenbosch, the Netherlands) using the manufacturer's protocol. Each chip contains a panel of 196 (PTK) and 144 (STK) phorylatable peptides. The phosphorylation of these peptides by native kinases in the applied lysate was quantified by sequence agnostic fluorescently coupled anti phosphorylation antibodies (Suppl. Fig. 4). Image processing was achieved using the instrument manufacturer's Bionavigator software before applying COMBAT normalization and performing Upstream Kinases Analysis (UKA) using a multiple hypothesis testing approach. This identified responsible kinases for the observed phosphorylation status with a median final score for the likelihood of the kinase being responsible for observed phosphorylation to which a cutoff of 1.2 was applied and a median kinase statistic indicating the degree of and direction of kinase activity change between compared conditions. The results were compiled to construct representative kinome trees with an adaptation of the CORAL tool [82] in R (v 4.2.0).

Network-based pathway enrichment
To identify biological pathways enriched in a protein list, we applied a network-based protein set enrichment analysis approach using Enrich-Net [83]. The inputs are a list of the PamGene protein identifiers and the database of interest from which reference protein sets will be extracted. The databases include (KEGG, BioCarta, WikiPathways, Reactome, PID, and GO. Using a genome-scale molecular interaction network, we map the target and reference datasets and then perform a network analysis, which involves two steps: i) the use of a random walk with restart (RWR) algorithm to score the distance between the mapped target protein set and reference datasets in the network. ii) the comparison of the calculated scores against a background model. Extensive information on how the RWR method for node relevance scoring was developed is found in [83]. A network similarity score (XD-score) measures how close the protein set of the PamGene data is to the pathways in the databases. The XD score is compared to a classical enrichment analysis score (pvalues) which measure the significance of the overlap between pathways and the protein set. A similarity ranking table of the reference datasets is generated, including their network-based association scores (measured using the XD-score). The table ranks the pathways from a chosen database to the protein set of the PamGene data. Using an interactive graph-based visualization, we generate interoperable interaction networks for each pathway. This network (XML-based) highlights the proteins that overlap in pathways, illustrating their interactions.

Small RNAseq
Small RNAseq was performed on the same samples used for total RNAseq. Libraries were prepared using "QIASeq miRNA Library kit" (Qiagen, Germany) and sequenced as 1x75bp on a NextSeq 550 (Illumina, USA). The sequencing was carried out at Luxgen platform, LNS. Data were aligned to the miRBase (v 2.2.2.1) and piRBase (regulatoryrna.org, v 2.0) using bowtie2 (v 2.3.4.1) and differential expression analysis was performed using R (v 4.2.2) and DESeq2 (v 1.34.0). This analysis was performed using a snakemake workflow adapted from a previous workflow designed by https://github.com/ jounikuj/ for this library preparation. Briefly, the snakemake workflow applies 4 passes of trimming (including UMIs) and map the obtained reads on miRBase using bowtie2 with the parameter -very-sensitivelocal. Next, the corresponding un-mapped reads are aligned to piRBase and the successfully mapped read coming from these 2 passes of alignment are merged in a unique BAM file per sample. Finally reads are counted, after UMI collapsing (UMI_tools v. 1.0.0) using a combination of samtools (v 1.9) and bash commands. qCLASH qCLASH sample preparations were performed exactly as previously described by Kozar et al. [16] Briefly, 50 million cells were UV cross-linked at a wavelength of 250nm and 600J/cm2. AGO protein immunoprecipitation was performed using Protein G Dynabeads (Invitrogen), and 10μl 2A8 anti-Ago antibody (Millipore) per sample. qCLASH libraries were sequenced on a NextSeq500 instrument with a read length of 85bp. The raw sequences were pre-processed and analyzed as previously described [16] using the bioinformatics pipeline hyb [84], and custom scripts available at GitHub (https://github.com/RenneLab/qCLASH-Analysis). miRNA network construction and representation qCLASH data were analysed using R (v 4.2.2) and the "tidyverse" package. Interaction lists were converted to "igraph" objects with the corresponding package (v 1.3.4) and further exported to cytoscape using "RCy3" (v 2.12.4). Final representations were made using Cytoscape (v 3.8.2) and the plugin "Legend Creator" (v 1.1.6).

RT-qPCR
RNA was extracted from cell lysates using the Quick-RNA MiniPrep kit with on-column DNase I treatment (Zymo Research, USA) according to the manufacturer's protocol. RNA purity and quality were assessed using the NanoDrop2000 Spectrophotometer (Thermo Scientific, USA). 500ng of RNA were used to perform reverse transcription using the miScript RT II Kit and HiFlex Buffer Reverse Transcription (Qiagen, Netherlands). Quantitative polymerase chain reaction (qPCR) was performed in technical triplicates using the CFX384 Detection System (BioRad, USA) using 5ng and 50ng of cDNA as a template for subsequent miRNA and mRNA expression detection, respectively. Primers were purchased from Qiagen (miScript Primer Assay range) and from Eurogentec (primer list Supplementary Table 4) to respectively amplify miRNA and mRNA. Amplification curves were analyzed using CFX software (BioRad, USA) and expression was normalised using three housekeeping genes (RNU1A, RNU5A, SCARNA17 for miRNA and HPRT, PPIA, TBP for mRNA). Statistical significance was calculated for each condition considering biological replicates (n = 3) using multiple welch's t-test and an FDR of 0.01.

Single cell RNA-seq analysis
Single cell RNA-seq data from [Randic et al., in revision] were reprocessed using Seurat (v 4.2.0) and SeuratObject (v 4.1.2). Reads were aligned to human genome hg38 using STAR aligner. Transcripts were filtered to be detected in at least 3 cells and to totalise at least 100 counts across cells. Cells were filtered based on mitochondrial transcripts (<=25%), ribosomal transcripts (<= 1%), and a cell complexity (log10GenesPerUMI >=80%). Thresholds for Unique Molecular Identifiers (UMI), genes detected were different between cells lines (for MELJUSO and SKMEL30: 150 and 100, for IPC298: 350 and 150). Finally, mitochondrial, ribosomal and hemoglobin genes were removed, and genes were filtered to have at least 10 counts across the filtered cells.
iCell data-fusion framework iCell methodology [15] was applied separately to each of the three cell lines (MELJUSO, SKMEL30 and IPC298) at day 0, day 4 and day 33, to obtain integrated networks that encompass the condition-specific (i.e., defined by a cell line and day) PPI, COEX, and single-cell RNA-seq data from our lab [Randic et al., unpublished]. We used the iCell methodology because it is an intermediate data integration method that uses Non-negative Matrix Tri-factorization (NMTF) to simultaneously decompose multiple relational matrices through the inference of a single joint model [15,85]. Therefore, it does not suffer from information loss of early integration approaches that combine all datasets into a single dataset before integration, or late integration methods that first build models for each dataset and then combine them into an integrated model. Additionally, NMTF-based methods are co-clustering methods that encode the high-dimensional input data into three low-dimensional matrix factors [86]. The clustering information is contained in the so-called cluster indicator matrices that are easily interpretable (unlike the resulting matrices of other artificial intelligence algorithms), making them ideal for mining biological data.
For each condition (treated cell lines and time points), we constructed condition-specific Protein-Protein Interaction (PPI) and gene coexpression (COEXP) networks as follows. We collected the experimentally validated protein-protein interactions of Homo sapiens from BioGRID (v 4.3.195) [32], from which we excluded genes (nodes in the network) that are not expressed in our condition-specific bulk RNA-seq data (zFPKM < −3, [87]). To generate the condition-specific COEXP networks, we applied Spearman's correlation between the gene expression profiles from the condition-specific bulk RNA-seq data, using a 0.9 correlation threshold with a p-value ≤ 0.01. These two conditionspecific networks were further filtered to exclude genes (nodes in the networks) that are not expressed in any single cell according to the corresponding condition-specific single cell RNA-seq data.
For each condition, we followed the iCell methodology [15] to integrate the condition-specific PPI, COEXP, and single-cell expression data together. To this aim, we used Non-negative Matrix Tri-Factorization to simultaneously decompose the adjacency matrices of the PPI and COEXP networks, A 1 and A 2 , respectively, and the single cell RNA-seq expression data, E. More precisely, the two adjacency matrices, A 1 and A 2 , are simultaneously decomposed as the products of three matrix factors each, A 1 as the product of matrix factors G 1 , S 1 and G T 1 , and A 2 as the product of matrix factors G 1 , S 2 and G T 1 , i.e.: A 1 % G 1 S 1 G T 1 and A 2 % G 1 S 2 G T 1 , where G 1 is interpreted as the cluster indicator matrix of the genes and matrices S 1 and S 2 are interpreted as the compressed representations of the PPI and COEXP networks, respectively. The conditionspecific single-cell RNA-seq expression matrix, E, is decomposed into the product of three matrix factors G 1 , S 3 and G T 2 as: E % G 1 S 3 G T 2 , where S 3 is interpreted as the compressed representation of the single-cell expression data, G 2 is the cluster indicator matrix of single cells, and G 1 is the cluster indicator matrix of the genes that is shared across the decompositions of all input matrices, which facilitates the information flow and allows for learning from all data.
Furthermore, to ease interpretability, we constrain all matrix factors to be non-negative. We also impose orthogonality constraint to the G 2 matrix factor to obtain less ambiguous functional organization captured by the clusters [86]. This decomposition is done by minimizing the following objective function, F: min G1;G2;S1;S2;S3 ! 0 where ‖ ‖ F denotes the Frobenius norm and I is the identity matrix. Because minimizing this objective function is an NP-hard continuous optimization problem [88], we heuristically solve it by using a fixed point method that, starting from an initial solution, iteratively applies multiplicative update rules [86] to converge towards a locally optimal solution (described in the next section).
After minimizing F, we used the obtained matrix factors to create an integrated network that encompasses the condition-specific PPI, COEXP, and single cell RNA-seq data, which we call iCell network. Following the iCell methodology [15], we created this network by computing G 1 Á G T 1 and thresholding the resulting matrix to preserve only the top 1% of the strongest relationships in each row and column.
We apply this methodology for each condition separately, resulting in nine condition-specific iCell networks.

Fixed point solver
Following the iCell methodology [15], first we derive the Karush-Kuhn-Tucker (KKT) conditions, which are necessary for a solution to be optimal [88]: where ⊙ is the Hadamard (elementwise) product, I is the identity matrix, and matrices, η 1 and η 2 , are the dual variables for the primal constraints G 1 ≥ 0 and G 2 ≥ 0, respectively. Then, we derive the following multiplicative update rules for all matrix factors, G 1 , G 2 , S 1 , S 2 and S 3 , to solve the KKT conditions presented above: In our fixed point approach, we start from initial solutions of the matrix factors, which we compute using singular value decomposition (SVD) of the original input matrices [88], and iteratively use the derived update rules to compute new matrix factors, G 1 , G 2 , S 1 , S 2 and S 3 , until convergence. Initializing the matrix factors with SVD reduces the number of iterations that are needed to achieve convergence and makes our solver deterministic [89]. To satisfy the non-negativity constraint of NMTF, we generate the initial solutions by taking the absolute values of the entries from the matrices computed using SVD.
The iterative process ends when the objective function does not decrease anymore, which we measure every ten iterations with: F À1 À F 0 F 0 10 À3 , where F 0 is the value of the objective function at the current iteration, i, and F −1 is the value computed at iteration i-10.

Analyzing the topology of iCells
To mine the biological knowledge hidden in the iCell networks, we applied graphlet-based methods because they are the most sensitive measures of network topology to date [28,29]. Graphlets are small, connected, nonisomorphic, induced subgraphs of a large network that appear at any frequencies in the network [88]. Within a graphlet, the symmetries between groups of nodes are called automorphism orbits and are used to characterise the different topological positions that a node participates in [90]. These orbits are used to generalise the notion of the node degree (the number of edges that the node touches in the network): for each orbit, the "graphlet degree of the node" in the network is the number of times the node touches (is found at the position of) the particular orbit [28]. Following the methodology of Yaveroğlu et al. (2014), we characterised the topology of each node in a network by using the 11 nonredundant orbits of 2-to 4-node graphlets. Therefore, the topology of every node is captured by an 11dimensional vector called "Graphlet Degree Vector" (GDV).
We compared the global topological dissimilarity between two condition-specific iCell networks by using "Graphlet Correlation Distance" (GCD-11), because it is the most sensitive network distance measure [28]. To do this, first we characterised the global topology of a network with its "Graphlet Correlation Matrix" (GCM), which is an 11 × 11 symmetric matrix that contains the Spearman's correlations between GDVs over all nodes of the network. Then, the GCD-11 between two networks is the Euclidean distance between the upper triangle values of their GCMs. The lower the GCD-11 value, the closer the topologies of the networks [28].
Additionally, we measure the change in the wiring of a gene between two condition-specific iCell networks by using the "Graphlet Degree Vector Similarity" (GDVS) between the GDVs of the node in the two networks, which is computed as follows. Given the GDV vector of a gene (node) in one network, a, and the GDV vector of the same gene (node) in another network, b, the distance between the i th coordinates of vectors a and b is defined as: where w i is the weight of orbit i that accounts for dependencies between the orbits [29]. Then, GDVS between vectors a and b is defined as: GDVS is the similarity measure in (0,1] range, such that the similarity equal to 1 means that the two GDV vectors are identical, and thus, that the gene has the same local topology in both networks. We use GDVS to identify the genes whose topology changed between two conditionspecific iCell networks, which we termed "perturbed" (i.e., whose GDVSs were close to 0), as well as the genes whose local topology did not change, which we termed "stable" (i.e., whose GDVSs were close to 1). These comparisons were made between different time points for the same cell line, as well as between different cell lines for the same time point.

Clustering
For each cell line and each time point, we performed clustering of the Graphlet Degree Vectors (GDVs). Prior to clustering, the features were weighted by using previously defined weights [29] and then logtransformed. In addition, both the NMTF and GCM matrices were scaled to a 0 mean and unit variance. For each input matrix, we first computed the Hopkins statistics to estimate cluster tendency. Then, to identify the numbers of clusters, we estimate the cluster validity of potential solutions reporting between 50 and 150 clusters, using 21 different metrics. We rank all clustering solutions for each metric. Then, for each input matrix, the number of clusters is selected based on the solution associated with the minimal median rank across all 21 metrics. We then performed hierarchical and k-medoids clustering using a Pearson correlation-based distance and validated the obtained repartitions using silhouette width and the between/within average distance ratio. All scripts are available on the associated Git repository.

Network representation of gene predictions
Protein-Protein interactions between TP53, a known regulator of senescence, and genes associated with "EMT" and "INTERFERON GAMMA RESPONSE" which have been linked with senescence escape were extracted from BioGRID (v 4.4.217). A network representation was created by using Cytoscape (v 3.9.1). The representation was further refined by using stringApp (v 2.0.0) and clusterMaker2 app (v 2.3.2) and modified by using Inkscape (v 1.2.1).