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

ARTICLE OPEN A multi-omics integrative approach unravels novel genes and pathways associated with senescence escape after targeted therapy in NRAS mutant melanoma Vincent Gureghian 1 , Hailee Herbst 1 , Ines Kozar 2 , Katarina Mihajlovic 3 , Noël Malod-Dognin 3 , Gaia Ceddia 3 , Cristian Angeli 1 , Christiane Margue 1 , Tijana Randic 1 , Demetra Philippidou 1 , Milène Tetsi Nomigni 1 , Ahmed Hemedan 4 , Leon-Charles Tranchevent 4 , Joseph Longworth 5 , Mark Bauer 1 , Apurva Badkas 1 , Anthoula Gaigneaux 1 , Arnaud Muller 6 , Marek Ostaszewski 4 , Fabrice Tolle 1 , NatašaPržulj 3,7,8,9 and Stephanie Kreis 1,9 ✉ © The Author(s) 2023 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. Cancer Gene Therapy; https://doi.org/10.1038/s41417-023-00640-z 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 Received: 24 March 2023 Revised: 19 May 2023 Accepted: 21 June 2023 1 Department of Life Sciences and Medicine, University of Luxembourg, 6, Avenue du Swing, L-4367 Belvaux, Luxembourg. 2 Laboratoire National de Santé, Dudelange, Luxembourg. 3 Barcelona Supercomputing Center, 08034 Barcelona, Spain. 4 Luxembourg Centre for Systems Biomedicine, University of Luxembourg, Esch-sur-Alzette, Luxembourg. 5 Experimental and Molecular Immunology, Department of Infection and Immunity, Luxembourg Institute of Health, Esch-sur-Alzette, Luxembourg. 6 LuxGen, TMOH and Bioinformatics platform, Data Integration and Analysis unit, Luxembourg Institute of Health, Esch-sur-Alzette, Luxembourg. 7 Department of Computer Science, University College London, London WC1E 6BT, UK. 8 ICREA, Pg. Lluís Companys 23, 08010 Barcelona, Spain. 9 These author contributed equally: NatašaPržulj, Stephanie Kreis. ✉email: [email protected] www.nature.com/cgt Cancer Gene Therapy 1234567890();,: [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, specificfirst 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. RESULTS 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 Epithelial- to-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–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 V. Gureghian et al. 2 Cancer Gene Therapy 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. V. Gureghian et al. 3 Cancer Gene Therapy 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–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 Fig. 2 NRAS mutant melanoma trigger IFN and EMT responses upon MEKi and CDK4/6i. A Principal component analysis clusters cell lines separately reflecting distinct observed phenotypes. BGene Set Enrichment Analysis using the hallmarks gene sets from MSigDB. Color legends represent Normalised Enrichment Score (NES). All cell lines display an enrichment in interferons responses and at day 14 in EMT. CSenescence 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]. V. Gureghian et al. 4 Cancer Gene Therapy 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. 1.09 1.02 0.98 0.95 0.97 0.94 0.88 0.82 0.81 0.82 -2.22 -1.61 -1.4 2.93 2.69 2.56 2.57 2.57 2.49 2.5 2.47 2.18 2.31 -1.58 -1.05 -0.68 1.32 1.35 1.65 1.46 1.41 1.56 1.12 1.16 1.28 1.25 -1.99 -1.91 -1.48 0.86 1.16 0.84 0.89 0.87 0.76 0.73 0.72 0.71 0.66 -0.84 -0.92 -0.95 4.77 3.43 3.27 3.17 3.02 3.05 2.94 2.93 2.76 2.61 -0.59 -0.47 -0.47 -0.11 -0.14 0.1 -0.08 -0.11 -0.08 -0.02 0.04 0.03 0.01 -1.02 -1.08 -0.57 MELJUSO SKMEL30 Day1 Day4 Day33 Day1 Day4 Day33 FAK2 FAK1 ERK5 Src IRR InSR IGF1R TRKC TRKB TRKA HER3 ERK2 ERK1 Time Kinase Score a a Non -significant Significant -4 -2 0 2 4 Median Kinase Statistic SKMel30 Day 1 AAK1 Abl Arg Ack ACtR2 ACtR2B Akt1/PKBα Akt2/PKBβ Akt3/PKBγ ALK ALK1 ALK2 ALK4 ALK7 AMPKα1 AMPKα2 ANKRD3 ANPα ANPβ ARAF AurA/Aur2 AurB/Aur1 AurC/Aur3 Axl BARK1/GRK2 BARK2/GRK3 BIKE BLK BMPR1A BMPR1B BMPR2 Etk/BMX BRAF Brk BRSK1 BRSK2 BTK Bub1 BubR1 CaMK1α CaMK1β CaMK1δ CaMK1γ CaMK2α CaMK2β CaMK2δ CaMK2γ CaMK4 CAMKK1 CAMKK2 caMLCK CASK CCK4/PTK7 CCRK CDC2/CDK1 CDC7 CDK10 CDK11 CDK2 CDK3 CDK4 CDK5 CDK6 CDK7 CDK8 CDK9 CDKL1 CDKL2 CDKL3 CDKL4 CDKL5 CHED CHK1 CHK2 CK1α CK1α2 CK1δ CK1ε CK1γ1 CK1γ2 CK1γ3 CK2α1 CK2α2 CLIK1 CLIK1L CLK1 CLK2 CLK3 CLK4 COT CRIK CRK7 CSK CTK CYGD CYGF DAPK1 DAPK2 DAPK3 DCAMKL1 DCAMKL2 DCAMKL3 DDR1 DDR2 DLK DMPK DMPK2 DRAK1 DRAK2 DYRK1A DYRK1B DYRK2 DYRK3 DYRK4 EGFR EphA1 EphA10 EphA2 EphA3 EphA4 EphA5 EphA6 EphA7 EphA8 EphB1 EphB2 EphB3 EphB4 EphB6 HER2 HER3 HER4 ERK1 ERK2 ERK3 ERK4 ERK5 ERK7 FAK Fer Fes FGFR1 FGFR2 FGFR3 FGFR4 Fgr FLT1 FLT3 FLT4 FmS/CSFR FRK Fused Fyn GAK GCK GCN2 GCN2~b GRK4 GRK5 GRK6 GRK7 GSK3α GSK3β Haspin HCK HGK/ZC1 HH498 HIPK1 HIPK2 HIPK3 HIPK4 HPK1 HRI HSER HUNK ICK IGF1R IKKα IKKβ IKKε ILK InSR IRAK1 IRAK2 IRAK3 IRAK4 IRE1 IRE2 IRR ITK JAK1 JAK1~b JAK2 JAK2~b JAK3 JAK3~b JNK1JNK2 JNK3 KDR KHS1 KHS2 KIS Kit KSR KSR2 LATS1 LATS2 Lck LIMK1 LIMK2 Lkb1 Lmr1 Lmr2 Lmr3 LOK LRRK1 LRRK2 LTK Lyn LZK MAK MEK1/MAP2K1 MEK2/MAP2K2 MKK3/MAP2K3 SEK1/MAP2K4 MAP2K5 MKK6/MAP2K6 MAP2K7 MEKK1/MAP3K1 MAP3K8 MEKK2/MAP3K2 MEKK3/MAP3K3 MAP3K4 ASK/MAP3K5 MEKK6/MAP3K6 MAP3K7 MAPKAPK2 MAPKAPK3 MAPKAPK5 MARK1 MARK2 MARK3 MARK4 MAST1 MAST2 MAST3 MAST4 MASTL MELK Mer Met MINK/ZC3 MISR2 MLK1 MLK2 MLK3 MLK4 MLKL MNK1 MNK2 MOK MOS MPSK1 MRCKα MRCKβ MSK1 MSK1~b MSK2 MSK2~b MSSK1 MST1 MST2 MST3 MST4 MUSK MYO3A MYO3B MYT1 NDR1 NDR2 Nek1 Nek10 Nek11 Nek2 Nek3 Nek4 Nek5 Nek6 Nek7 Nek8 Nek9 NIK NIM1 NLK NRBP1 NRBP2 NRK/ZC4 NuaK1 NuaK2 Obscn Obscn~b OSR1 P38α P38β P38δ P38γ p70S6K p70S6Kβ PAK1 PAK2 PAK3 PAK4 PAK5 PAK6 PAS K PBK PCTAIRE1 PCTAIRE2 PCTAIRE3 PDGFRα PDGFRβ PDK1 PERK/PEK PFTAIRE1 PFTAIRE2 PhKγ1 PhKγ2 PIK3R4 Pim1 Pim2 Pim3 PINK1 PITSLRE PKAα PKAβ PKAγ PKCα PKCβ PKCδ PKCε PKCγ PKCη PKCι PKCθ PKCζ PKD1 PKD2 PKD3 PKG1 PKG2 PKN1/PRK1 PKN2/PRK2 PKN3 PKR PLK1 PLK2 PLK3 PLK4 PRKX PRKY PRP4 PRPK PSKH1 PSKH2 PYK2 QIK QSK RAF1 Ret RHODK/GRK1 RIPK1 RIPK2 RIPK3 RNAseL ROCK1 ROCK2 Ron ROR1 ROR2 Ros RSK1/p90RSK RSK1~b RSK2 RSK2~b RSK3 RSK3~b RSK4 RSK4~b RSKL1 RSKL2 RYK SBK SCYL1 SCYL2 SCYL3 SgK069 SgK071 SgK085 SGK1 SgK110 SgK196 SGK2 SgK223 SgK269 SgK288 SGK3 SgK307 SgK396 SgK424 SgK493 SgK494 SgK495 SgK496 SIK skMLCK SLK Slob smMLCK SNRK SPEG SPEG~b Src Srm SRPK1 SRPK2 SSTK STK33 STLK3 STLK5 STLK6 SuRTK106 Syk TAK1 TAO 1 TAO 2 TAO 3 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 TXK Tyk2 Tyk2~b Tyro/Sky ULK1 ULK2 ULK3 ULK4 VACAMKL VRK1 VRK2 VRK3 Wee1 Wee1B WNK1 WNK2 WNK3 WNK4 YANK1 YANK2 YANK3 Yes YSK1 ZAK ZAP70 STE CK1 AGC CAMK CMGC TK TKL Branch & Node Colour Median Kinase Statistic -4 0 4 Node Size Median Final Score 1.2 3.6 RTKs involved in enriched pathways: : ErbB signaling : Neurotrophin signaling : Insulin signaling BRAF PDGFRB FGF4 PDGFRANF1 MAPK7 HRAS FGFR3 AC D MELJUSO SKMEL30 E PDGFRA FGF4 FGFR3 BRAF MAPK7 HRAS PDGFRBNF1 B Fig. 3 MEKi and CDK4/6i lead to RTKs reactivation. A Summary of significantly deregulated kinases. BKinase activities after MEKi and CDK4/ 6i. ERK1/2 inhibition leads to RTKs activation. Rreceptors names are colored according to the enriched pathways. CPhylogenetic tree for SKMEL30 on day 1. The size of the leafs represents the “Median Final Score”, a score greater than 1.2 indicates a significant change between conditions. The color of the branches and leaves shows the “Median Kinase Statistic”, which is the difference in kinase activity. Circles highlight RTKs contributing to enriched pathways. D,EActivity of ERK5 in the senescent and adaptative cell lines. DZoom on the MAPK pathway showing the deregulation of ERK5 (MAPK7) in MELJUSO cell line. Esame for the SKMEL30 cell line. V. Gureghian et al. 5 Cancer Gene Therapy 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 Table 1. Enrichment results from PamGene data. Pathway Name Shigellosis NOD-like receptor signaling pathway Epithelial cell signaling in Helicobacter pylori infecon Progesterone-mediated oocyte maturaon RIG-I-like receptor signaling pathway ErbB signaling pathway GnRH signaling pathway Type II diabetes mellitus Fc epsilon RI signaling pathway Neurotrophin signaling pathway Toll-like receptor signaling pathway Pancreac cancer Colorectal cancer T cell receptor signaling pathway VEGF signaling pathway Acute myeloid leukemia Prostate cancer Leishmaniasis Dorso-ventral axis formaon Day1 XD-score 1.83 1.38 1.12 1.02 1.00 0.97 0.93 0.91 0.90 0.89 0.89 0.76 0.75 0.73 0.69 0.57 0.56 0.48 0.48 1.41E–14 1.99E–11 2.54E–13 1.47E–14 6.12E–20 6.69E–14 3.49E–13 p-adj 6.22E–06 1.13E–02 8.02E–07 6.70E–13 4.52E–17 5.76E–14 1.12E–08 1.52E–07 2.23E–12 6.13E–09 9.12E–08 9.34E–08 Pathway Name XD-score mTOR signaling pathway 1.13 Acute myeloid leukemia 1.05 Type II diabetes mellitus 0.89 Glioma 0.86 ErbB signaling pathway 0.86 Gap juncon 0.77 Prion diseases 0.75 Adherens juncon 0.72 Fc gamma R-mediated phagocytosis 0.71 VEGF signaling pathway 0.68 Fc epsilon RI signaling pathway 0.65 Prostate cancer 0.65 Long-term potenaon 0.65 Aldosterone-regulated sodium reabsorpon 0.64 Tight juncon 0.62 Vascular smooth muscle contracon 0.60 Insulin signaling pathway 0.59 Neurotrophin signaling pathway 0.59 Progesterone-mediated oocyte maturaon 0.58 GnRH signaling pathway 0.55 Adipocytokine signaling pathway 0.53 Hedgehog signaling pathway 0.52 Long-term depression 0.51 T cell receptor signaling pathway 0.47 Chemokine signaling pathway 0.44 Day4 Pathway Name Type II diabetes mellitus mTOR signaling pathway ErbB signaling pathway Acute myeloid leukemia Fc epsilon RI signaling pathway Progesterone-mediated oocyte maturaon Glioma Prion diseases Adipocytokine signaling pathway VEGF signaling pathway Shigellosis Neurotrophin signaling pathway Aldosterone-regulated sodium reabsorpon Gap juncon Prostate cancer Adherens juncon GnRH signaling pathway T cell receptor signaling pathway Long-term potenaon NOD-like receptor signaling pathway Pancreac cancer Fc gamma R-mediated phagocytosis Toll-like receptor signaling pathway Epithelial cell signaling in Helicobacter pylori infecon Long-term depression Vascular smooth muscle contracon Insulin signaling pathway Colorectal cancer Day33 p-adj 5.13E–08 7.58E–09 5.45E–08 1.48E–06 1.33E–10 1.50E–09 3.91E–06 8.12E–08 2.02E–10 3.34E–09 2.78E–04 7.95E–05 6.47E–08 2.89E–05 1.69E–07 2.99E–10 7.43E–07 2.18E–07 4.38E–14 3.96E–04 6.96E–11 4.24E–09 1.38E–10 2.88E–10 8.17E–07 XD-score 1.42 1.08 1.01 1.00 0.95 0.95 0.95 0.95 0.87 0.86 0.86 0.83 0.82 0.82 0.80 0.79 0.77 0.75 0.72 0.71 0.69 0.65 0.64 0.63 0.59 0.55 0.54 0.54 p-adj 1.49E–09 2.57E–11 3.31E–12 3.13E–10 4.51E–12 4.23E–12 1.69E–08 7.28E–11 1.21E–09 3.65E–08 9.13E–11 3.13E–11 2.77E–07 2.91E–05 2.60E–08 5.62E–10 1.53E–07 5.38E–13 4.81E–05 2.33E–08 9.87E–10 3.31E–05 1.45E–06 2.77E–07 1.68E–09 4.58E–09 2.77E–07 1.26E–07 Pathway enrichment results for MELJUSO cell line Pathway Name Shigellosis Fc epsilon RI signaling pathway NOD-like receptor signaling pathway Progesterone-mediated oocyte maturaon Type II diabetes mellitus RIG-I-like receptor signaling pathway Epithelial cell signaling in Helicobacter pylori infecon GnRH signaling pathway Neurotrophin signaling pathway ErbB signaling pathway Toll-like receptor signaling pathway VEGF signaling pathway Aldosterone-regulated sodium reabsorpon T cell receptor signaling pathway Leishmaniasis Adherens juncon Colorectal cancer Prion diseases Dorso-ventral axis formaon Amyotrophic lateral sclerosis (ALS) Day1 XD-score 1.27 1.17 1.12 0.95 0.94 0.91 0.88 0.87 0.86 0.80 0.75 0.73 0.69 0.69 0.65 0.65 0.64 0.54 0.52 0.48 1.09E–13 2.92E–16 p-adj 5.18E–11 2.85E–09 9.10E–05 2.45E–11 1.91E–07 1.99E–13 4.51E–07 3.10E–10 2.11E–10 2.04E–13 3.07E–16 2.12E–04 2.25E–07 1.78E–06 1.08E–03 8.88E–03 9.43E–11 1.70E–11 Pathway Name XD-score Acute myeloid leukemia 1.27 Glioma 0.91 mTOR signaling pathway 0.83 ErbB signaling pathway 0.80 Type II diabetes mellitus 0.74 Adipocytokine signaling pathway 0.71 Long-term potenaon 0.69 Fc gamma R-mediated phagocytosis 0.66 Gap juncon 0.60 Prostate cancer 0.59 Fc epsilon RI signaling pathway 0.58 Vascular smooth muscle contracon 0.57 Prion diseases 0.54 Adherens juncon 0.52 Neurotrophin signaling pathway 0.49 Apoptosis 0.49 Tight juncon 0.46 Non-small cell lung cancer 0.46 Vibrio cholerae infecon 0.46 Aldosterone-regulated sodium reabsorpon 0.46 Long-term depression 0.43 Salivary secreon 0.42 Progesterone-mediated oocyte maturaon 0.41 Gastric acid secreon 0.40 VEGF signaling pathway 0.37 Day4 Pathway Name Shigellosis NOD-like receptor signaling pathway Dorso-ventral axis formaon Pancreac cancer Type II diabetes mellitus Epithelial cell signaling in Helicobacter pylori infecon Progesterone-mediated oocyte maturaon Fc epsilon RI signaling pathway GnRH signaling pathway ErbB signaling pathway Toll-like receptor signaling pathway RIG-I-like receptor signaling pathway T cell receptor signaling pathway Leishmaniasis Colorectal cancer Endometrial cancer Neurotrophin signaling pathway VEGF signaling pathway 1.14E–04 Day33 8.19E–10 4.31E–08 3.50E–03 6.55E–06 6.05E–05 6.05E–05 2.79E–04 p-adj 2.29E–09 3.74E–07 3.74E–07 6.98E–08 7.94E–05 7.74E–04 2.42E–08 2.70E–03 1.93E–05 2.15E–08 4.20E–06 5.47E–08 1.43E–09 1.52E–09 4.02E–05 1.06E–06 7.97E–07 XD-score 1.17 1.02 0.94 0.85 0.80 0.79 0.78 0.75 0.74 0.65 0.62 0.58 0.57 0.57 0.54 0.51 0.47 0.42 p-adj 9.97E–05 3.13E–09 8.77E–06 5.44E–09 1.04E–09 6.16E–07 1.31E–09 6.62E–07 9.30E–06 1.64E–11 3.88E–10 2.69E–04 1.06E–09 2.89E–06 1.70E–08 4.11E–10 1.53E–09 1.54E–11 Pathway enrichment results for SKMEL30 cell line : Neurotrophin signaling RTKs involved in enriched pathways: : ErbB signaling : Insulin signaling Network-based enrichment was performed on the significantly deregulated kinases at different time points. MELJUSO and SKMEL30 cells show at all time points an enrichment in ErbB and neurotrophin pathway. Only MELJUSO cells show an enrichment in insulin pathway at day 4 and day 33. V. Gureghian et al. 6 Cancer Gene Therapy 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 5739 4363 5150 1550 1045 738 2268 1265 5357 9935 10522 4746 5607 443 1104 1156 13986 13999 12497 5038 1395 752 992 2803 759 420 723 723 SKMEL30 resistant sensitve IPC298 resistantsensitve ∩ ( - ) miRNA-mRNA network (nb of nodes, nb of edges) + RNA-seq ( - ) A C E B LGALS3BP TXNIP TYRP1 0 1 2 5 10 15 20 mRNAs expression across cell lines MelJuso Sensitive MelJuso Resistant SKMel30 Treatment Sensitive SKMel30 Treatment Resistant IPC298 Treatment Sensitive IPC298 Treatment Resistant 624Mel Treatment Sensitive 624Mel Treatment Resistant A375 Treatment Sensitive A375 Treatment Resistant miR-211 0 1 2 3 4 miR-211-5p expression across cell lines MelJuso Sensitive MelJuso Resistant SKMel30 Treatment Sensitive SKMel30 Treatment Resistant IPC298 Treatment Sensitive IPC298 Treatment Resistant 624Mel Treatment Sensitive 624Mel Treatment Resistant A375 Treatment Sensitive A375 Treatment Resistant 0.4-2.3 1.7 3.0 Node Fill Color: log2FoldChange -1.0 ASAP1 CDK2 CASP4 ACTG1 ANXA5 BRI3 miR-320a let-7f ASAH1 TXNIP MYC let-7a PLXNA1 let-7c GOLM1 ZNF92 ADNP WDFY1 miR-27a UGCG TOP1 UBAP2 RPS3A TPP1 RBM25 TNRC6B ILF2 TECPR2 ERC1 STAT2 EIF5B SLC20A1 CRTAP SH3PXD2A ATP6V1B2 RNFT1 miR-211 RAB5B EIF4E2 NSD1 miR-181a MYO5A TOP2A MTF2 RPS28 MLXIP NOX1 LIMD2 MT-CO3 LGALS3BP GPR143 LAMC1 miR-146a JKAMP TYRP1 IFI6 SQSTM1 HSPA8 SPARC HSP90AB1 TMF1 RNF31 GBF1 SGK1 KPNA4 FTL PPT1 KPNA2 ELOVL5 NFIA FKBP1A CCND1 MELJUSO SKMEL30 IPC298 untreatreated untreatreated untreatreated 64 128 1024 1 2 4 8 condition untreated treated miR-211-5p miRNA D GF condition count V. Gureghian et al. 7 Cancer Gene Therapy 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-specificiCell 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 andSKMEL30cells,thetop10%ofthemost“perturbed”genes were enriched in “INTERFERON ALPHA/GAMMA RESPONSE”.Atdays4and 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”,“COAGULATION”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. BTop 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. CTo 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. DRNA-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. ERelative expression of LGALS3BP, TXNIP and TYRP1 mRNAs assessed by qPCR across three NRAS and two BRAF mutant melanoma cell lines. FRelative expression of miR-211-5p across the same cell lines as determined by qPCR. GExpression of miR-211-5p in the small RNA-seq. V. Gureghian et al. 8 Cancer Gene Therapy Group stable rewired Genes eniched 10 20 30 D4 vs D0 D33 vs D4 D33 vs D0 rewired rewired rewired ADIPOGENESIS COAGULATION E2F TARGETS EMT MTORC1 SIGNALING MYC TARGETS V1 OXIDATIVE PHOSPHORYLATION D4 vs D0 D33 vs D4 D33 vs D0 rewired stable rewired stable rewired stable ADIPOGENESIS E2F TARGETS MYOGENESIS OXIDATIVE PHOSPHORYLATION D4 vs D0 D33 vs D4 D33 vs D0 rewired stable rewired stable rewired stable ADIPOGENESIS OXIDATIVE PHOSPHORYLATION PI3K AKT MTOR SIGNALING ROS PATHWAY iCELL networks: Network topology: Gene topology: MJS day0 MJS MJS day4 MJS day33 SK day0 SK day4 SK day33 IPC day0 IPC day4 IPC day33 0.00 0.35 0.31 1.06 0.59 2.04 2.38 1.50 1.56 0.35 0.00 0.26 1.02 0.85 2.24 2.59 1.66 1.61 0.31 0.26 0.00 0.96 0.77 2.13 2.49 1.57 1.51 1.06 1.02 0.96 0.00 1.04 1.64 1.98 1.02 0.83 0.59 0.85 0.77 1.04 0.00 1.56 1.91 1.07 1.36 2.04 2.24 2.13 1.64 1.56 0.00 0.71 0.69 1.20 2.38 2.59 2.49 1.98 1.91 0.71 0.00 1.09 1.63 1.50 1.66 1.57 1.02 1.07 0.69 1.09 0.00 0.87 1.56 1.61 1.51 0.83 1.36 1.20 1.63 0.87 0.00 MJS_D33 IPC_D4 IPC_D33 IPC_D0 SK_D33 SK_D0 SK_D4 MJS_D0 MJS_D4 MJS_D33 IPC_D4 IPC_D33 IPC_D0 SK_D33 SK_D0 SK_D4 MJS_D0 MJS_D4 0 0.5 1 1.5 2 2.5 IPC_D0 99 IPC_D4 95 IPC_D33 61 97 100 131 10607 MJS_D0 135 MJS_D4 158 MJS_D33 143 174 30 471 10377 SK_D0 135 SK_D4 48 SK_D33 230 48 644 103 9756 IPC298MELJUSO SKMEL30 SK IPC IPC_D33 IPC_D4 IPC_D0 SK_D33 SK_D4 SK_D0 MJS_D33 MJS_D4 MJS_D0 ADIPOGENESIS ALLOGRAFT REJECTION ANDROGEN RESPONSE APICAL JUNCTION APICAL SURFACE APOPTOSIS BILE ACID METABOLISM CHOLESTEROL HOMEOSTASIS COAGULATION DNA REPAIR E2F TARGETS EMT FATTY ACID METABOLISM G2M CHECKPOINT HEME METABOLISM IL6 JAK STAT3 SIGNALING INFLAMMATORY RESPONSE INTERFERON ALPHA RESPONSE INTERFERON GAMMA RESPONSE KRAS SIGNALING DN KRAS SIGNALING UP MITOTIC SPINDLE MTORC1 SIGNALING MYC TARGETS V1 MYC TARGETS V2 MYOGENESIS OXIDATIVE PHOSPHORYLATION P53 PATHWAY PEROXISOME PI3K AKT MTOR SIGNALING PROTEIN SECRETION ROS PATHWAY SPERMATOGENESIS TGF BETA SIGNALING TNFA SIGNALING VIA NFKB UNFOLDED PROTEIN RESPONSE UV RESPONSE DN WNT BETA CATENIN SIGNALING XENOBIOTIC METABOLISM Count 10 20 30 40 50 Enriched clusters using k-medoids Enrichment gene clustering based on topology: Clusters Number of clusters enriched with k-medois Number of clusters enriched with hierarchical clustering A B C D E F 0 500 1000 1500 1234 number of genes hierarchical clustering 0 500 1000 12345 number of clusters enriched number of clusters enriched number of genes k-medoids V. Gureghian et al. 9 Cancer Gene Therapy AUTHOR CONTRIBUTIONS SK conceived the study. SK, FT, DP, MT, TR and VG contributed to the experimental design. VG, FT, and DP acquired total RNA-seq and small RNA-seq data and AM, AG and VG processed, analysed and interpreted the data. TR, DP and CM acquired PamGene data, JL processed the data and performed Upstream Kinase Analysis. MO and AH performed network-based pathway enrichment and created xml models to overlay the inferred kinase activities; VG and CA interpreted the results. IK and HH performed qCLASH method, IK processed the data and AG, VG and HH analysed the hybrids and constructed the miRNA network. HH performed qPCR validations. MB processed the scRNA-seq data. NP, NMD, GC and KM conceived the data integration. KM, GC, and VG performed it and analysed the results under the supervision of NMD and NP. LT performed the gene clustering; AB, MB, and VG analyzed the gene clusters. VG, HH, AH, CA, AB and FT contributed to the figures. VG and SK drafted the manuscript with input from all authors. NP, NMD and MO reviewed the manuscript. All authors read and approved the final manuscript. FUNDING VG is supported by the Luxembourg National Research Fond (FNR) PRIDE DTU CanBIO [grant reference: 21/16763386]. TR is supported by the FNR PRIDE DTU CriTiCS [grant reference: 10907093]. Project-related work performed by VG, HH, CM, DP, MTN, MB, AG, FT and SK were also supported by the University of Luxembourg and the Fondation Cancer, Luxembourg (grant “SecMelPro”). KM and NP are supported by funding from the European Union’s EU Framework Programme for Research and Innovation Horizon 2020, Innovative Training Networks (MSCA-ITN-2019), funded under EXCELLENT SCIENCE - Marie Skłodowska-Curie Actions, Grant Agreement No 860895. KM, NMD, GC and NP are supported by funding from the European Research Council (ERC) Consolidator Grant 770827. NP is also supported by funding from the Spanish State Research Agency AEI 10.13039/501100011033 grant number PID2019-105500GB-I00. COMPETING INTERESTS The authors declare no competing interests. ADDITIONAL INFORMATION Supplementary information The online version contains supplementary material available at https://doi.org/10.1038/s41417-023-00640-z. Correspondence and requests for materials should be addressed to Stephanie Kreis. Reprints and permission information is available at http://www.nature.com/ reprints Publisher’s note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations. Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http:// creativecommons.org/licenses/by/4.0/. © The Author(s) 2023 V. Gureghian et al. 16 Cancer Gene Therapy