KDM5-driven transcriptional noise fuels plasticity-led awakening and relapse in paediatric cancer

KDM5-driven transcriptional noise fuels plasticity-led awakening and relapse in paediatric cancer | bioRxiv /* */ /* */ <!-- <!-- /*! * yepnope1.5.4 * (c) WTFPL, GPLv2 */ (function(a,b,c){function d(a){return"[object Function]"==o.call(a)}function e(a){return"string"==typeof a}function f(){}function g(a){return!a||"loaded"==a||"complete"==a||"uninitialized"==a}function h(){var a=p.shift();q=1,a?a.t?m(function(){("c"==a.t?B.injectCss:B.injectJs)(a.s,0,a.a,a.x,a.e,1)},0):(a(),h()):q=0}function i(a,c,d,e,f,i,j){function k(b){if(!o&&g(l.readyState)&&(u.r=o=1,!q&&h(),l.onload=l.onreadystatechange=null,b)){"img"!=a&&m(function(){t.removeChild(l)},50);for(var d in y[c])y[c].hasOwnProperty(d)&&y[c][d].onload()}}var j=j||B.errorTimeout,l=b.createElement(a),o=0,r=0,u={t:d,s:c,e:f,a:i,x:j};1===y[c]&&(r=1,y[c]=[]),"object"==a?l.data=c:(l.src=c,l.type=a),l.width=l.height="0",l.onerror=l.onload=l.onreadystatechange=function(){k.call(this,r)},p.splice(e,0,u),"img"!=a&&(r||2===y[c]?(t.insertBefore(l,s?null:n),m(k,j)):y[c].push(l))}function j(a,b,c,d,f){return q=0,b=b||"j",e(a)?i("c"==b?v:u,a,b,this.i++,c,d,f):(p.splice(this.i++,0,a),1==p.length&&h()),this}function k(){var a=B;return a.loader={load:j,i:0},a}var l=b.documentElement,m=a.setTimeout,n=b.getElementsByTagName("script")[0],o={}.toString,p=[],q=0,r="MozAppearance"in l.style,s=r&&!!b.createRange().compareNode,t=s?l:n.parentNode,l=a.opera&&"[object Opera]"==o.call(a.opera),l=!!b.attachEvent&&!l,u=r?"object":l?"script":"img",v=l?"script":u,w=Array.isArray||function(a){return"[object Array]"==o.call(a)},x=[],y={},z={timeout:function(a,b){return b.length&&(a.timeout=b[0]),a}},A,B;B=function(a){function b(a){var a=a.split("!"),b=x.length,c=a.pop(),d=a.length,c={url:c,origUrl:c,prefixes:a},e,f,g;for(f=0;f<d;f++)g=a[f].split("="),(e=z[g.shift()])&&(c=e(c,g));for(f=0;f<b;f++)c=x[f](c);return c}function g(a,e,f,g,h){var i=b(a),j=i.autoCallback;i.url.split(".").pop().split("?").shift(),i.bypass||(e&&(e=d(e)?e:e[a]||e[g]||e[a.split("/").pop().split("?")[0]]),i.instead?i.instead(a,e,f,g,h):(y[i.url]?i.noexec=!0:y[i.url]=1,f.load(i.url,i.forceCSS||!i.forceJS&&"css"==i.url.split(".").pop().split("?").shift()?"c":c,i.noexec,i.attrs,i.timeout),(d(e)||d(j))&&f.load(function(){k(),e&&e(i.origUrl,h,g),j&&j(i.origUrl,h,g),y[i.url]=2})))}function h(a,b){function c(a,c){if(a){if(e(a))c||(j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}),g(a,j,b,0,h);else if(Object(a)===a)for(n in m=function(){var b=0,c;for(c in a)a.hasOwnProperty(c)&&b++;return b}(),a)a.hasOwnProperty(n)&&(!c&&!--m&&(d(j)?j=function(){var a=[].slice.call(arguments);k.apply(this,a),l()}:j[n]=function(a){return function(){var b=[].slice.call(arguments);a&&a.apply(this,b),l()}}(k[n])),g(a[n],j,b,n,h))}else!c&&l()}var h=!!a.test,i=a.load||a.both,j=a.callback||f,k=j,l=a.complete||f,m,n;c(h?a.yep:a.nope,!!i),i&&c(i)}var i,j,l=this.yepnope.loader;if(e(a))g(a,0,l,0);else if(w(a))for(i=0;i (function(w,d,s,l,i){w[l]=w[l]||[];w[l].push({'gtm.start':new Date().getTime(),event:'gtm.js'});var f=d.getElementsByTagName(s)[0];var j=d.createElement(s);var dl=l!='dataLayer'?'&l='+l:'';j.src='//www.googletagmanager.com/gtm.js?id='+i+dl;j.type='text/javascript';j.async=true;f.parentNode.insertBefore(j,f);})(window,document,'script','dataLayer','GTM-M677548'); Skip to main content Home About Submit ALERTS / RSS Search for this keyword Advanced Search New Results KDM5-driven transcriptional noise fuels plasticity-led awakening and relapse in paediatric cancer Alejandro Allo Anido , Cecilia Roux , Emilia Chen , Siân Hamer , Abigail Shea , Ayeh Sadat Sadr , Christie English , Charlotte Butterworth , Harvey Che , Angella Bellini , Birgit Geoerger , Gudrun Schleiermacher , Louis Chesler , Michael David Morgan , Alejandra Bruna doi: https://doi.org/10.1101/2025.10.13.682004 Alejandro Allo Anido 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Cecilia Roux 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Emilia Chen 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Siân Hamer 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Abigail Shea 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: alejandra.bruna{at}icr.ac.uk Ayeh Sadat Sadr 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Christie English 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Charlotte Butterworth 4 SIRiC RTOP Laboratory (Recherche Translationelle en Oncologie Pediatrique), U1330 INSERM and SIREDO Integrated Pediatric Oncology Center, PSL research University, Institut Curie ; Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Harvey Che 3 Division of Clinical Studies, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Angella Bellini 4 SIRiC RTOP Laboratory (Recherche Translationelle en Oncologie Pediatrique), U1330 INSERM and SIREDO Integrated Pediatric Oncology Center, PSL research University, Institut Curie ; Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Birgit Geoerger 5 Department of Pediatric and Adolescent Oncology, Gustave Roussy Cancer Campus, Université Paris-Saclay , Villejuif, France ; INSERM U1015 Find this author on Google Scholar Find this author on PubMed Search for this author on this site Gudrun Schleiermacher 4 SIRiC RTOP Laboratory (Recherche Translationelle en Oncologie Pediatrique), U1330 INSERM and SIREDO Integrated Pediatric Oncology Center, PSL research University, Institut Curie ; Paris, France Find this author on Google Scholar Find this author on PubMed Search for this author on this site Louis Chesler 3 Division of Clinical Studies, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Michael David Morgan 6 Institute of Medical Sciences, School of Medicine, Medical Sciences and Nutrition, University of Aberdeen , Scotland, United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site Alejandra Bruna 1 Centre for Evolution and Cancer, Division of Cancer Biology, The Institute of Cancer Research; Sutton , United Kingdom 2 Centre for Children and Young People’s Cancer, The Institute of Cancer Research; Sutton , United Kingdom Find this author on Google Scholar Find this author on PubMed Search for this author on this site For correspondence: alejandra.bruna{at}icr.ac.uk Abstract Full Text Info/History Metrics Supplementary material Preview PDF ABSTRACT How drug-tolerant persister (DTP) cells escape quiescence to drive tumour relapse is a central unresolved question in cancer evolution. Here, we identify transcriptional noise (TN), defined as the stochastic variability in gene expression, as a latent property of paediatric cancer cells that becomes a driver of adaptive regrowth after treatment withdrawal. Using functional assays, lineage tracing, single-cell transcriptomics, and multiscale landscape modelling, we show that therapy enriches mesenchymal-like tolerant states in neuroblastoma without clonal selection, while post-treatment awakening is a stochastic process fuelled by noise-enabled plasticity in cell-identity programmes. The histone demethylase KDM5A relocates to noisy cell-state genes during awakening, promoting H3K4me3 removal and chromatin remodelling at these loci. KDM5 inhibition abrogates this process, and suppresses transcriptional noise, halts DTP exit, and prevents tumour recovery in both neuroblastoma and hepatoblastoma models. These results establish DTP as an exploitable evolutionary bottleneck, positioning KDM5-mediated transcriptional noise as an actionable therapeutic target to limit cancer adaptation and relapse Download figure Open in new tab INTRODUCTION Cancer is the leading cause of death in children, accounting for more than one in five deaths among individuals under 14 1 . Unlike adult cancers, which often arise from an accumulation of somatic mutations, paediatric cancers typically emerge in utero due to disrupted differentiation and uncontrolled growth of transformed embryonic cells 2 . These malignancies are driven by relatively few genetic alterations 3 , 4 , yet they also frequently evade treatment due to their adaptive evolution capabilities. The embryonic cells from which paediatric cancers arise are inherently plastic during development, a property essential for their ability to differentiate into diverse cell types 5 – 9 . The developmental origin of childhood cancers, coupled with their low mutational burden, suggests that non-genetic mechanisms, rather than the extensive acquisition of novel mutations, play a key role in treatment resistance and relapse. As such, paediatric cancers offer a unique model for studying plasticity-driven adaptive evolution. Neuroblastoma is one of the most common paediatric cancers, affecting approximately 100 children every year in the United Kingdom 1 . Half of these children will present with “high-risk” disease and half of them will succumb to neuroblastoma due to treatment resistance or relapse 10 , 11 . As a paradigm of paediatric cancer, neuroblastoma exhibits limited genetic drivers, usually restricted to MYCN amplification and/or structural alterations in chromosome 17 and 1 (accounting to up to 65% high-risk patients). Other common drivers include alterations in TERT, ATRX and ALK 12 – 14 . Genomic profiling of relapsed tumours has revealed further acquisition of alterations in these loci, but lack of shared novel drivers 15 , highlighting the contribution of non-genetic mechanisms to neuroblastoma evolution in patients. Despite this, neuroblastoma displays remarkable phenotypic heterogeneity. It remains unclear whether this arises from an imprinting of developmental plasticity and/or different cells of origin within the sympatho-adrenal lineage – the tissue of origin of these tumours 16 – 19 . Two cell identities have been reported and well characterized in neuroblastoma pre-clinical models 20 . The noradrenergic identity (ADRN) – resembling differentiated neuroblasts – and the mesenchymal state (MES) – associated with treatment resistance and which some groups have linked to dedifferentiated Schwann cell precursors 16 . These ADRN and MES states are regulated by exclusive sets of transcription factors that form core regulatory circuitries (CRC) 21 – 23 . Recent studies support plasticity between ADRN and MES identities in pre-clinical models 24 , yet its impact on patients remains debatable. While some works question the presence of MES cells in patient tumours 19 , 25 , 26 , others suggest they only emerge in response to treatment 21 , 23 . Studying the phenotypic dynamics and mechanisms driving phenotypic plasticity post-treatment is essential to understand the processes fuelling resistance and relapse. Although specific signalling pathways could have a role in promoting plasticity 24 , the mechanism by which these seemingly spontaneous changes in transcriptional programs are carried out in relatively homogenous populations remains largely elusive. Transcriptional noise is the stochastic variability in gene expression among individual cells in a genetically homogenous population. It mainly arises from the discontinuous nature of gene expression, with bursts of transcription followed by valleys of no expression. The asynchrony of this process generates stochastic expression differences amongst cells at any given timepoint 27 , 28 . Transcriptional noise has been associated with phenotype diversification and cell-fate determination in multiple organisms, and it is considered a “bet-hedging” strategy to ensure survival in rapidly changing environments. It allows HIV to exit the latent phase 29 , it promotes a competent phenotype in bacteria 30 , it has a role in the development of organs in mouse and the fruit fly 31 , 32 , and more recently it has been linked to the development of resistance in neuroblastoma 33 . We hypothesized that transcriptional noise drives plasticity and adaptive evolution in neuroblastoma and other cancers by increasing the likelihood of cells to spontaneously rewire their transcriptional programs towards fitter phenotypes. Here, we apply a preclinical treatment framework in neuroblastoma cell lines and patient-derived models to study post-treatment adaptation. Integrating lineage barcoding and single-cell transcriptomic data we show that treatment enriches tolerant states with MES-like phenotype with limited clonal selection, whereas post-treatment awakening and regrowth is driven by stochastic transitions into proliferative ADRN identities. This process is enabled by a surge in transcriptional noise in cell-identity genes during awakening, coinciding with an increased binding of the histone demethylase KDM5A and a contraction of H3K4me3 marks at noisy loci. Notably, we extended these observations to hepatoblastoma models and show that pharmacological inhibition of KDM5 activity reduces noise and delays awakening and recovery in both neuroblastoma and hepatoblastoma. Collectively, our findings establish transcriptional noise-driven plasticity as a key evolutionary force underlying neuroblastoma progression after treatment and propose KDM5A as and actionable target to halt post-treatment adaptation in paediatric cancers. RESULTS Highly plastic bridge cells mediate heterogeneity in treatment-naïve populations Phenotypic heterogeneity is a well-established determinant of cancer evolution, observed in both preclinical and clinical models across many tumour types 34 – 36 , including genomically silent paediatric cancers 4 . However, its origins remain debated: whether it arises from branching clonal evolution, distinct cells of origin, or phenotypic plasticity. Here, we asked whether plasticity is intrinsic to naïve cancer cell populations and relevant to early adaptation. We focused on neuroblastoma, the most common extracranial solid paediatric cancer, known for extensive phenotypic heterogeneity 21 , 22 , 37 . To investigate intrinsic plasticity, we profiled treatment-naïve neuroblastoma populations using single cell RNA-sequencing (scRNA-seq) across four models: the SK-N-SH cell line and three patient-derived models (the xenograft GR-NB5 and the organoid lines NB039 and NB067). Clustering analyses confirmed heterogeneity across all models (Figure S1A-D). However, previous studies in preclinical and clinical samples have yielded conflicting results, due to the complexity of the neuroblastoma phenotypic landscape and the limitations of the traditional clustering method 16 , 18 , 25 , 38 , 39 . While most evidence supports the presence of at least three major cell phenotypic populations reminiscent of different stages of sympathoadrenal development 38 , 40 , reconciling these findings across patient tumours, patient-derived models, and long-term cultured cell lines remains challenging. We first annotated our models based on canonical adrenergic (ADRN) and mesenchymal (MES) states, which have been previously characterised as established states in neuroblastoma cell line models 21 , 41 (Figure S1E-H). ADRN cells are highly proliferative and governed by core regulatory circuits involving transcription factors such as PHOX2B, GATA3, and HAND2, which drive neuronal differentiation programs. MES cells are regulated by AP-1, PRRX1, and RUNX1, promoting a mesenchymal-like, migratory, and therapy-resistant phenotype. A population of cells expressing mixed adrenergic and mesenchymal features was also observed, hereon in referred to as bridge (BRI) cells. Patient-derived organoids only contained cells of the ADRN identity, highlighting the limitations of cell line-derived signatures to recapitulate clinical heterogeneity 37 , 38 , 40 . Thus, we used annotations derived from primary neuroblastomas 40 and, besides ADRN cells, we identified populations expressing proliferation markers (PROL) as well as rare cells expressing both canonical ADRN markers and patient-derived MES markers, which we classified as BRI (Figure S1I-J). Traditional clustering methods impose discrete boundaries that limit the assessment of phenotypic heterogeneity and plasticity as a continuum. To resolve this, we applied MuTrans, a multiscale dynamical manifold approach, to construct three-dimensional phenotypic landscapes and capture the dynamic properties of cell state transitions 42 ( Figure 1 ). In this approach, each cell was arranged within the 3D landscape and their position within attractor wells provide insight into both the stability of cell states and the dynamics of transitions within the landscape. Cells within deep attractor wells display low transcriptomic entropy and are committed to a cell state, while those traversing over shallower valleys show high entropy and are predicted to represent transitioning plastic states ( Figure 1A-H ). We observed that BRI cells in SK-N-SH and GR-NB5 display higher entropy than MES and ADRN states (Figure S1K-L), suggesting a higher transitioning potential for these cells. In line with this, plastic wells generally contained the highest percentage of BRI cells, while committed attractor states were mainly composed of MES, ADRN or PROL cells. For instance, in SK-N-SH, clusters 1, 2 and 4 were annotated as plastic based on their entropy levels, and contained 35.5%, 9.5% and 62.8% BRI cells, respectively, while all other clusters contained only 1.71-4.8% of cells in this state. Similarly, in GR-NB5 plastic clusters contained 10.3-19.7% BRI cells, compared to 2.3-6.6% in committed wells. Also in NB039 and NB065 BRI cells accounted for 5.8-27.7% and 9.6-59.2% of the cells in plastic wells, respectively ( Figure 1I-L ). Interestingly, these two models exhibited very divergent landscapes despite containing similar cell identities (Figure SF1). NB039, with wild-type ALK , was characterised by a prototypical landscape including plastic and committed cells, whereas cell states in NB067, with mutant ALK , were highly unstable, coinciding with an enrichment of BRI cells in this model ( Figure 1 ). This may suggest a contribution of ALK status to plasticity potential in neuroblastoma. Download figure Open in new tab Figure 1. Phenotypic heterogeneity in neuroblastoma models is mediated by highly plastic BRI cells. (A-D) Dynamic manifold landscapes in untreated neuroblastoma models. Well depth correlated with stability of the attractor state. Cells are coloured by attractor well and were annotated as “plastic” or “committed” based on the entropy levels in (E-H). (E-H) Entropy values per attractor state, estimated from transcriptomic data using MuTrans. (I-L) Percentage cells in each cell state (ADRN/BRI/MES for SK-N-SH and GR-NB5; ADRN/BRI/PROL for NB039 and NB067) across attractor wells. To functionally validate our computational inferences, we prospectively isolated and studied distinct cell states in SK-N-SH cells using the cell surface marker CD44 24 (Figure S2A). CD44 low cells expressed ADRN markers and exhibited the highest proliferation rates, while CD44 high cells displayed a MES phenotype and were more resistant to standard chemotherapies, consistent with prior observations 21 , 22 , 24 (Figure S2B-G). Isolation of pure ADRN, MES and BRI (CD44 int ) populations demonstrated bidirectional plasticity between these cell states (Figure S3A-B). BRI cells were highly unstable and displayed progenitor-features. Significantly, expression of BRI markers in neuroblastoma patients correlated with poorer survival, underscoring the malignant potential of this transient, plastic phenotype (Figure S3C-E). Collectively, these results provide functional evidence that phenotypic interconversion occurs in treatment-naïve populations, a process we define as “intrinsic plasticity”. They also identify BRI cells as transitional states actively mediating shifts between more stable cellular identities, acting as biomarkers of plasticity. Importantly, the presence of pre-existing drug-tolerant phenotypes in untreated populations supports the view that phenotypic diversity, driven by plasticity rather than genetic mutation, functions as an evolutionary bet-hedging strategy, safeguarding a reservoir of pre-adapted states capable of surviving future therapeutic pressures. Therapy selects phenotypes, but awakening is driven by plasticity and stochastic clonal dynamics To investigate plasticity at single-cell resolution in naïve, unperturbed populations and during the early stages of adaptive evolution under therapy, we applied an expressed DNA barcoding approach within a preclinical treatment framework ( Figure 2A ). Download figure Open in new tab Figure 2. Lineage tracing reveals stochastic awakening of DTP cells. (A) SK-N-SH cells were tagged with Cellecta barcodes in a lentiviral system and expanded. The population (referred to as POT) was divided into triplicates and from each of these, the cells were distributed equally into six flasks. Among the six replicates, three were treated with IC50 doses of cisplatin, and the remaining three replicates served as untreated controls. After one week of treatment, 10% of the cells from each replicate were collected and cryopreserved (T1). The remaining cells in all replicates were re-seeded and left untreated. The collection process was repeated one week (T2) and four weeks (T4) later. (B) Proportion of lineages containing one unique cell state, two or three cell states in untreated samples. Plots display mean ± SD of three biological replicates. (C). Distribution of cells across cell states and barcode lineages in POT and a representative replicate of the untreated samples, for the top 100 most frequent lineages. (D) Shannon diversity index per experimental timepoint. Bars represent average of three replicates. (E) Cumulative frequency of each barcode lineage across experimental timepoints. Barcodes were ranked by increased frequency. (F) Area under the curve for the cumulative frequencies shown in (E). based on comparisons Statistical significance was determined using pairwise t-tests and Bonferroni multiple-comparison correction methods between each experimental condition and untreated samples (*P ≤ 0.05, **P ≤ 0.01, ***P ≤ 0.001). (G) Frequency of the 100 most highly abundant lineages within POT visualized across all other samples and replicates, ordered most to least frequent. Dote size reflects percentage representation of each barcode within a sample. Individual SK-N-SH cells were marked with unique genomic barcodes detectable in both DNA and RNA sequencing platforms, enabling simultaneous tracking of clonal and phenotypic dynamics (Materials and Methods). In our experimental setup, one barcode represents one clone and its progeny, as the DNA barcode. In this system, each barcode uniquely marks a clone and its progeny, as the barcode is stably integrated into the genome and inherited across cell divisions. To ensure high lineage resolution and experimental reproducibility in single-cell analyses, the barcode library was designed to trace approximately 100,000 potential lineages and was expanded 850-fold prior to experimentation to maximise the recovery of identical barcodes across independent replicates. At this point, a sample of the population – comprised of an average of 16,124 unique barcodes – was taken, hereafter referred to as the pooled population of transduced cells (POT), and the remaining cells were subsequently used to establish biological replicates for downstream analyses. To investigate intrinsic plasticity at clonal resolution, we integrated barcode and transcriptomic data, linking lineage identity with transcriptionally defined cell states. This analysis enabled the assignment of phenotypic composition to each clone ( Figure 2B and 2C ). We observed that approximately 50% of the clones in the naïve population comprised cells occupying more than one cell state, indicative of intrinsic plasticity. For example, in the POT sample, 40.5% of cellular clones were linked to a single phenotypic state, 17.6% to two states, and 41.9% to three states, demonstrating that phenotypic transitions occur in a substantial fraction of unperturbed cells. Having established a method to quantify clonal plasticity, we next investigated how plasticity interacts with treatment-induced selection using a preclinical treatment framework that features on- and off-treatment cycles ( Figure 2A ). Barcoded SK-N-SH cells were exposed to cisplatin for seven days (T1) and subsequently cultured in drug-free medium for one week (persistence; T2) or four weeks (awakening; T3). Single-cell transcriptomic analysis revealed that treatment enriched for drug-tolerant persister (DTP) cells displaying MES-like phenotypes, while ADRN-like cells were depleted. Upon drug withdrawal, these DTP-MES cells re-entered the cell cycle and regrew by transitioning to highly proliferative ADRN-like states, restoring phenotypic heterogeneity (Figure S4A-C). The existence of intrinsic plasticity in untreated populations ( Fig. 2B–C ), suggests that treatment may predominantly select pre-existing MES-like phenotypes rather than inducing new ones. To validate this, we sorted ADRN and MES populations using CD44-based flow cytometry and exposed them to cisplatin. Cisplatin preferentially eliminated ADRN cells, while MES cells persisted, subsequently undergoing MES-to-ADRN transitions upon drug withdrawal driving regrowth and phenotypic diversification (Figure S4D). These results establish that treatment selects phenotypes, but post-treatment recovery depends on plasticity-led transitions. We next investigated how these phenotypic transitions are reflected at the clonal level. Barcode sequencing revealed minimal clonal selection during treatment (T1), as indicated by highly similar clonal identities and frequency distributions across replicates. In contrast, awakening and regrowth (T3) were accompanied by a marked reshaping of clonal diversity and cumulative frequency distributions ( Figure 2D-F ). Only a small, replicate-specific subset of clones contributed disproportionately to regrowth, indicating a stochastic recovery process ( Figure 2G ). Ranking of lineages by counts and logarithmic frequency confirmed that low-frequency clones at recovery (T3) followed a linear, Poisson-like distribution, consistent with stochastic awakening, whereas persistent clones at T1 displayed non-linear distributions, consistent with treatment selection (Figure S4E). At higher frequencies, this relationship was lost due to the exponential expansion of early-awakening lineages. Altogether, these results reveal a two-step adaptive strategy under strong therapeutic pressure: (i) phenotypic selection seeds persistence through survival of tolerant MES-like states and (ii) plasticity-driven awakening, fuelled by stochastic processes, re-establishes phenotypic heterogeneity and drives tumour regrowth. Transcriptional noise concentrates in transitional bridge states Transcriptional noise (TN), defined as stochastic variability in gene expression, is a well-described source of phenotypic diversification in homogeneous populations 29 , 30 , 43 , 44 . Given the stochastic nature of awakening of neuroblastoma DTPs, we asked whether TN contributes to intrinsic and post-treatment plasticity in neuroblastoma preclinical models. To address this, we analysed scRNA-seq data from barcoded SK-N-SH cells as described in Figure 1 . Using the BASiCS framework 45 , 46 to separate technical from biological variability, we focused on genes with stable mean expression across populations, thereby avoiding variability independent of noise. When comparing the highly plastic BRI cells to the more stable MES and ADRN cells, we consistently observed a higher level of noise (i.e. expression overdispersion) in BRI cells. Of the total number of genes that were classified as highly variable, 64% and 89% were noisier in the BRI compared with the MES and ADRN cells, respectively ( Figure 3A ; Table S2) Download figure Open in new tab Figure 3. Transcriptional noise promotes phenotypic plasticity in cell lines and patient-derived models. (A) Differential transcriptional noise (overdispersion) analysis between BRI and MES cells (left) or BRI and ADRN cells (right). Dots represent individual genes, dashed lines statistical thresholds and numbers percentage of all noisy genes in each cell state. (B) GSEA of BRI-specific noisy genes in (A). ( C) Representative flow cytometry plots of SK-N-SH cells treated for five days with DMSO, IdU or Cisplatin. Coloured boxes represent CD44 gating strategy and numbers represent the percentage of cells in that gate. (D) Representative flow cytometry histograms across conditions and timepoints (Iso, isotype; Cis, cisplatin). (E) Quantification of the percentage of ADRN (CD44low), BRI (CD44int) and MES (CD44high) cells during IdU or cisplatin treatment. (F) Wound healing assay in SK-N-SH treated with the noise inducer IdU or the plasticity inducer TGF-ß. Lines represent average wound area across two replicates and shaded area represents standard error (G) Representative images of patient-derived xenograft (PDX) before and after being treated ex vivo with IdU, TGF-ß or a combination of both. (H) Quantification of the area, roundness and number of fusions between PDX spheroids in (G). (I) Gene Set Enrichment Analysis (GSEA) of genes differentially expressed between spheroids treated with a combination of IdU+TGF-ß and any other treatment. (J) GSEA of public MES and ADRN signatures in the genes differentially expressed spheroids treated with a combination of IdU+TGF-ß and any other treatment. Stars represent statistical significance resulting from a two-way (E) or one-way (F, H) ANOVA test with Tukey’s HSD (‘****’ p <0.0001 ‘***’ p <0.001, ‘**’ p<0.01, ‘*’ p 0.05’). Gene set enrichment analysis (GSEA) further showed that noisy genes in BRI versus MES were enriched for biological processes related to cell motility and locomotion, consistent with MES-associated features ( Figure 3B ). Conversely, genes with higher noise in BRI compared with ADRN were enriched for neurogenesis, neuronal development, and proliferation pathways, reflecting ADRN-associated programmes. Together, these data demonstrate that intrinsic plasticity in neuroblastoma is linked to elevated TN, and that noise disproportionately affects genes associated with cell-state identity. This suggests that TN weakens state boundaries and provides a transcriptional substrate for phenotypic transitions. Noise induction promotes plasticity in cells and patient-derived models To test whether TN directly enhances plasticity, we employed a validated TN readout by transducing a destabilised GFP (d2GFP) reporter into SK-N-SH cells 43 . We treated d2GFP-expressing SK-N-SH cells with the TN inducer, IdU 43 , and analysed noise and phenotypic changes by flow cytometry. IdU increased transcriptional noise (d2GFP expression variability) in MES cells (Figure S5A-B) and enriched the population for this phenotype ( Figure 3C-E ), similar to the effect observed in response to cisplatin during the adaptation process. Remarkably, noise-induction also enriched the population for BRI cells. We also observed a similar enrichment of BRI cells when we generated noise through transient pausing of transcription using the transcriptional elongation inhibitor, 5,6-dichloro-1-β-D-ribofuranosylbenzimidazole (DRB) 47 , 48 , confirming the association with noise and not DNA damage (Figure S5C). Altogether, these results demonstrate a direct link between the induction of transcriptional noise and an increase in phenotypic plasticity, reflected by enrichment of plastic BRI states and transitions toward MES cells. We next evaluated the functional consequences of TN induction. In wound-healing assays, IdU treatment increased the migratory potential of neuroblastoma cells in a dose-dependent manner ( Figure 3F and Figure S5D). The magnitude of migration after 48h was comparable to that induced by TGF-β, a well-established driver of the most well-studied plasticity process known as epithelial-to-mesenchymal transition 49 , 50 . Based on the prior knowledge that MES cells display higher migratory and invasive potential than ADRN cells 24 , these results further confirmed that TN induced an increase in resistant and aggressive MES-like cells within neuroblastoma populations. We next tested whether TN could induce plasticity in patient-derived neuroblastoma cells grown ex vivo in 3D cultures. Cells were treated with IdU, TGF-β, or a combination of both. Single treatments had only modest effects, likely reflecting an increased phenotypic fidelity in established tumour cells 24 (Figure S5E). However, a significantly different phenotype was observed when we combined the effects of both plasticity (TGF-ß) and TN (IdU) enhancers ( Figure 3G-H ). Treated spheroids became larger, more irregularly shaped, and showed a tendency to fuse with neighbouring spheroids, although the latter did not reach statistical significance. This morphological change was accompanied by transcriptional evidence of plasticity, including activation of EMT-related pathways and enrichment of MES-associated gene signatures, as demonstrated by RNA sequencing of treated spheroids ( Figure 3I-J ). Altogether, this evidence reveals that elevating TN boosts plasticity, enriching transitional states and tolerant-like phenotypes in both cell lines and patient-derived models. Transcriptional noise peaks at persister exit and increases variability in cell-identity genes We have shown that treatment of phenotypically diverse neuroblastoma populations with cisplatin selects for slow-cycling DTP cells with MES-like features. Additionally, we demonstrated that upon cisplatin withdrawal plasticity drives their awakening and reconstitution of the original heterogeneous tumour population (Figure S4). Given the close association between TN and plasticity, we next asked whether plasticity-led recovery from cisplatin similarly depends on changes in TN. To address this, we quantified TN using BASiCS on scRNA-seq data from barcoded SK-N-SH cells treated with cisplatin ( Figure 2 ). Analyses were restricted to MES-classified cells, thereby eliminating confounding effects of phenotypic heterogeneity during recovery, as MES represents the only cell state detected immediately after treatment removal. In untreated populations, 2.7% of genes were classified as highly variable, consistent with intrinsic plasticity. This fraction decreased during cisplatin persistence (1.3%) but rose sharply upon treatment withdrawal during awakening and regrowth (8.0%), when plasticity is at its maximum. Once the population recovered and heterogeneity was restored, TN levels dropped back to baseline (3.0%) ( Figure 4A-B ; Table S3). Importantly, TN fluctuations were not correlated with cell-cycle state or viability, indicating that changes in noise reflect regulatory variability rather than proliferative differences (Figure S6A-D). These results establish that TN and plasticity are tightly coupled during therapy response, positioning TN as a key regulator of plasticity-led adaptive evolution in neuroblastoma. Download figure Open in new tab Figure 4. Phenotypic plasticity during awakening of DTPs is associated with transcriptional noise in cell identity genes. (A) High Variable Genes (HVG) analysis in MES cells across stages of cisplatin treatment. Dots represent individual genes, dashed lines statistical threshold and numbers percentage of the whole transcriptome that corresponds to HVG genes. Boxes label important MES (orange) and ADNR (purple) master regulators. (B) Percentage of noisy genes across treatment stages. (C) Overlap of high noise genes across stages of cisplatin treatment. (D) GSEA of high noise genes in all treatment stages (“Core”) or only during awakening (“Awakening”). Highlighted text relates to cell-state identity pathways in MES and ADRN cells. To identify genes underlying this phenomenon, we characterised the noisy gene pool across treatment stages. We uncovered a core set of 432 genes that remained highly variable before and after treatment ( Figure 4C ). These genes were enriched for plasticity-associated processes (e.g. EMT) and cell identity programs (e.g. cell adhesion, motility, neurogenesis) ( Figure 4D ). During awakening, TN expanded to include hundreds of additional genes involved in similar hallmarks and biological processes, including important regulators of the ADRN and MES identities, namely, PHOX2A , PHOX2B , KLF7 and PRRX1 , RUNX1 , RUNX2 , YAP1 , respectively 21 , 22 ( Figure 4 ). Notably, these genes exhibited minimal changes in mean expression and no consistent directionality, underscoring that TN is the primary driver of expression variability (Figure S4E-G). Together, these findings suggest that TN precedes and facilitates plasticity-led recovery by increasing expression variability in lineage regulators, thereby enhancing the likelihood that MES DTPs transition into proliferative ADRN states once therapy pressure is relieved. Transcriptional noise dynamics during therapy response are conserved across patient-derived models We then sought to evaluate whether the link between TN and plasticity-led recovery observed in preclinical cell lines also applies to patient-derived neuroblastoma models. To this end, we analysed scRNA-seq data of two neuroblastoma patient-derived organoid models (PDOs; NB039 and NB067) and a patient-derived xenograft (PDX; GR-NB5) cultured ex vivo ( Figure 1 and Figure S1). Despite differences in the gene sets defining cellular states across models, treatment-induced phenotypic dynamics were highly conserved ( Figure 5A and S7A-B). To facilitate cross-model comparisons, we annotated cell states based on their dynamics in response to treatment rather than by gene identity. Cisplatin reduced highly proliferating cells (ADRN in PDX; PROL in PDOs; hereafter “PROL”) and increased treatment-resistant reservoir cells (MES in PDX; ADRN in PDOs; collectively “RES”). Upon drug withdrawal, the awakening phase was characterised by the enrichment of previously rare bridge cells (BRI), marking phenotypic transitions towards repopulation and restoration of heterogeneity ( Figure 5B ). Consistent with this functional annotation, BRI cells exhibited the highest TN levels, further supporting their role as biomarkers of plasticity (Figure S7C). Download figure Open in new tab Figure 5. Transcriptional noise enables plasticity in patient-derived neuroblastoma models. (A) Phenotypic dynamics of PDO (NB039, NB067) and PDX (GR-NB5) models during cisplatin treatment, derived from scRNA-seq analysis. (B) Schematic of the functional annotation of cell states during cisplatin treatment across neuroblastoma models and their corresponding cell identity (as in Figure S1) (C) High Variable Genes (HVG) analysis in RES cells across stages of cisplatin treatment. Dots represent individual genes, dashed lines statistical threshold and numbers percentage of the whole transcriptome that corresponds to HVG genes. (D) Percentage of noisy genes across treatment stages and models. (E) GSEA of high noise genes before treatment and during awakening (“Core”) or only during awakening (“Awakening”) in PDOs. To study the involvement of TN in this process, we used BASiCS to estimate the levels of noise in RES cells during the different stages of treatment, as this is the dominant DTP phenotype upon cisplatin induction, akin to the MES phenotype in pre-clinical cell lines ( Figure 5B-D ; Table S4). In both organoid models, TN increased markedly during the awakening phase preceding recovery, with 4.5% (NB039) and 4.9% (NB067) of genes classified as highly variable, compared with 2.4% and 1.5%, respectively, in the pre-treatment condition. In contrast, RES cells from the GR-NB5 model did not display significant changes in TN, consistent with the lack of recovery observed upon drug withdrawal ( Figure 5A ). Shared noisy genes between cell lines and patient-derived organoids were limited and mainly involved in cell cycle regulation (Figure S7D). Consistent with our findings in cell lines, we identified a core gene set that maintain high noise levels before and after treatment, encompassing genes involved in cell-cycle progression, translation, and neuronal phenotype regulation. During awakening, this noisy gene pool expanded substantially, incorporating a broader cohort of genes associated with neuroblastoma cell identity programmes ( Figure 5E and S7E) Together, these results demonstrate that increased TN during the awakening phase is a conserved feature of plasticity-led recovery across neuroblastoma models. This elevated transcriptional variability affects key regulators of neuronal and adrenergic identity, suggesting that noise in lineage-specifying genes serves as a molecular driver of persister exit and tumour repopulation following therapy withdrawal. We hypothesized that the expression fluctuations generated by TN must be transiently stabilised across cell divisions to meaningfully contribute to plasticity-led evolution. To test this, we applied MemorySeq analysis 51 to SK-N-SH cells to distinguish between heritable and non-heritable gene expression variability. A total of 27 single-cell clones were derived and expanded from parental SK-N-SH cells and bulk RNA-seq was performed on each clone. If the high expression states generated by transcriptional noise were heritable, we would expect variability in the proportion of high-expressing cells across clones, depending on the timing of state acquisition. Conversely, if expression differences were short-lived, the proportion of high-expressors would remain consistent across clones ( Figure 6A ). Download figure Open in new tab Figure 6. Expression states generated by transcriptional noise are heritable. (A) Schematic of the MemorySeq experiment. A total of 27 single-cell clones were established from SK-N-SH, expanded and RNA-seq was performed to distinguish between heritable and non-heritable patterns of expression. Additionally, 48 samples of a mixed population were used as controls to account for technical noise. (B) Overlap between genes classified as noisy during awakening from cisplatin treatment and genes with heritable expression patterns in the MemorySeq experiment. (C) Genes identified as highly variable in the MemorySeq experiment. (D) Visualization of expression variability between clones for selected heritable genes (related to neuroblastoma cell-state identity) and non-heritable genes. MemorySeq identified 250 genes with heritable expression patterns in the clonal populations, compared to only 26 in a mixed-population control accounting for technical noise. Remarkably, 55.6% of the heritable genes overlapped with those highly variable during awakening from cisplatin treatment, including canonical MES and ADRN cell identity markers ( Figure 6B-D ; Table S5). These results indicate that the expression patters generated by transcriptional noise in cell identity genes can be captured and stabilised over several generations, supporting a role for stochastic gene expression in driving sustained phenotypic plasticity. A targeted screen identifies KDM5 as a suppressor of noise This observed link between TN and phenotypic plasticity suggested that modulating noise could represent a novel strategy to limit adaptive evolution in neuroblastoma. We therefore aimed to identify regulators of noise that could reduce the likelihood of phenotypic transitions and, in turn, limit cancer cell adaptability in response to treatment. To this end, we measured TN levels in SK-N-SH cells expressing a d2GFP reporter following treatment with a focused panel of 14 compounds, including drugs currently used in neuroblastoma management and agents targeting pathways of interest ( Figure 7A ; Table S1). After 5 days of treatment, we observed a strong correlation between TN (quantified as the standard deviation of d2GFP expression) and plasticity potential (measured as the proportion of BRI cells) ( Figure 7B ). Importantly, no correlation was detected between mean expression and variability of d2GFP, excluding confounding effects of mean expression on noise (Figure S8A). Download figure Open in new tab Figure 7. KDM5 inhibition reduces transcriptional noise and halts awakening of persister cells after cisplatin treatment. (A) Schematic of drug screening pipeline in reporter SK-N-SH. (B) Correlation between the change in standard deviation of the reporter protein d2GFP (proxy for transcriptional noise) and the percentage of bridge cells (CD44int, marker for plasticity) for each compound in (A). Both variables were assessed by flow cytometry. (C) Distribution of d2GFP expression in each treatment condition compared to vehicle-treated control cells; representative of two replicates. Both treatment groups represent the same number of cells and individual expression values are normalized to the average across all cells. Numbers represent the fold change in d2GFP between treated and untreated cells. (D) Change in ADRN and MES cells during stages of cisplatin treatment in the absence of the KDM5 inhibitor KDM5-C70 or when it is added only once or weekly after cisplatin removal. Data points were normalized to the percentage of each phenotype on the day of cisplatin removal. Stars represent statistical significance resulting from a one-way ANOVA test with Tukey’s HSD (‘****’ p <0.0001 ‘***’ p <0.001, ‘**’ p<0.01, ‘*’ p 0.05’). Most compounds had little or no effect on TN. Four drugs increased noise (fold change in standard deviation >1.5): the positive control IdU, the reactive oxygen species inducer auranofin, and the WNT pathway modulators XAV939 (inhibitor) and CHIR99021 (activator), suggesting these pathways may act as noise enablers. The only compound that reduced TN was KDM5-C70, a pan-KDM5 histone demethylase inhibitor ( Figure 7C ). KDM5A and KDM5B are epigenetic regulators that remove histone 3 lysine 4 trimethylation (H3K4me3) marks at gene promoters. The impact on d2GFP variability in response to their inhibition, therefore, points to a mechanistic link between chromatin regulation and control of TN. Notably, KDM5A/B have been shown to have critical roles during development and cancer 52 . In neuroblastoma, KDM5B is one of the MES signature genes 21 , highlighting the key role of this epigenetic regulator in MES cell state identity. Consistent with this, we detected an increase in KDM5A/B protein expression upon cisplatin withdrawal, indicating that KDM5 activity is dynamically upregulated during recovery and may facilitate adaptive reprogramming following treatment (Figure S8B). KDM5 inhibition blocks awakening and regrowth To test whether KDM5-mediated control of TN could regulate plasticity-led awakening and regrowth post-treatment, we used our preclinical treatment framework encompassing persistence, awakening, and recovery phases ( Figure 2A ). KDM5-C70 was administered either once or weekly following cisplatin withdrawal. In the control arm, cells followed the expected dynamics of DTP entry and exit, with MES-like persisters giving rise to proliferative ADRN phenotypes during awakening and recovery. By contrast, KDM5 inhibition markedly delayed this recovery, with weekly dosing producing the strongest effect ( Figure 7D ). Pre-treatment with KDM5-C70 before cisplatin withdrawal produced a similar outcome (Figure S8C), highlighting the potential benefit of reducing noise both during chemotherapy and in the adjuvant setting to prevent adaptation. Importantly, in untreated SK-N-SH cells, KDM5-C70 showed minimal toxicity and did not alter phenotypic composition, excluding selective cytotoxicity as an explanation (Figure. S8D–E). Together, these findings indicate that KDM5 activity promotes plasticity-led exit from drug-tolerant persister states, facilitating tumour regrowth. Inhibiting KDM5, and thereby reducing transcriptional noise, impedes the phenotypic transitions required for awakening, highlighting a potential therapeutic strategy to target cancer’s adaptive evolution and prevent recurrence and relapse. KDM5A primes noisy identity genes to enable awakening To dissect the mechanism by which KDM5 activity regulates plasticity-led awakening and regrowth, we analysed the chromatin-binding landscapes of KDM5A, KDM5B, and H3K4me3 in noisy loci during awakening. KDM5B binding remained low and largely unchanged throughout drug-response phases ( Figure 8C and S9A). In sharp contrast, KDM5A binding increased following cisplatin exposure and rose further during the awakening phase after drug withdrawal. This increase was significantly reduced by KDM5-C70 treatment, and the same patterns were observed in ADRN, BRI and MES identity genes ( Figure 8A , 8C and S9B-C). These results indicate that KDM5A directly contributes to plasticity-led evolution through chromatin engagement. Download figure Open in new tab Figure 8. KDM5A binding drives noise during awakening by compacting H3K4me3 peaks in SK-N-SH cells. (A) KDM5A and (B) H3K4me3 occupancy at noisy genes during awakening across stages of cisplatin treatment. Genes are ordered and coloured by signal intensity, normalised to the untreated condition. (C) Signal intensity at noisy genes during awakening across cisplatin treatment stages for KDM5A, KDM5B and H3K4me3. (D) Width of H3K4me3 peaks in noisy genes during awakening across cisplatin treatment stages. We next investigated whether these binding dynamics affected H3K4me3 deposition, the canonical substrate of KDM5A. Although global H3K4me3 levels at noisy genes increased during persistence, they remained largely unchanged during awakening and upon KDM5 inhibition ( Figure 8B-C ). However, we observed a significant decrease in H3K4me3 peak width during awakening, which was rescued by KDM5 inhibition ( Figure 8D ). This mechanism, rather than overall methylation, has been previously associated with increased transcriptomic heterogeneity 53 . Altogether, these results indicate that KDM5A primes chromatin at noisy, cell identity genes to facilitate plasticity-led phenotypic transitions. By selectively binding and reshaping H3K4me3 architecture at plasticity-linked loci, KDM5A orchestrates transcriptional noise and enables the awakening and regrowth of drug-tolerant persisters following treatment. Transcriptional noise drives phenotypic plasticity in hepatoblastoma models Finally, we wanted to investigate whether transcriptional noise could mediate plasticity and treatment recovery in other paediatric cancers. To this end, we used similar analytical approaches in scRNA-seq data of the hepatoblastoma cell lines HuH6 and HepG2 previously generated by our lab 54 . Hepatoblastomas are amongst the paediatric cancers with the lowest mutational burden 4 and recent studies have suggested that phenotypic plasticity drives adaptation to therapy in this tumour type 55 , 56 . Previous work from our group has shown that hepatoblastoma cell lines exist in a phenotypic continuum (similar to neuroblastoma) and that phenotypic plasticity between these states drives adaptation to treatment 54 . In hepatoblastoma models, however, these semi-stable phenotypes recapitulate transcriptomic features of normal liver development. In the HuH6 cell line this included a slow-proliferating precursor-like reservoir state (RES) resembling early progenitors, two differentiated and proliferative states (PROL-I, PROL-II) resembling differentiated hepatocytic cell states and a population of bridge (BRI) cells between these two more stable phenotypes with progenitor features 54 (Figure S10A). Quantification of the levels of noise using BASiCS, showed a significant increase in highly variable genes in BRI cells (3.2% of all genes), when compared to RES (0.8%), PROL-I (0.8%) and PROL-II states (0.6%) ( Figure 9A ). Notably, most noisy genes identified in the other states were also classified as highly variable in BRI cells, and were enriched for functions related to plasticity processes, development and hepatoblastoma cell identity (Figure S10B-C). These findings suggest that, as in neuroblastoma, BRI cells exhibit the highest levels of TN, consistent with their higher phenotypic plasticity. In addition, we show that TN in these transitioning cells particularly affects genes associated with dynamic cell state regulation. Download figure Open in new tab Figure 9. Transcriptional noise regulates plasticity in hepatoblastoma models. (A) High Variable Genes (HVG) analysis across untreated cell states or (B) in RES cells across stages of cisplatin treatment in HuH6. Dots represent individual genes, dashed lines statistical threshold and numbers percentage of the whole transcriptome that corresponds to HVG genes. Boxed highlight genes involved in plasticity-led treatment adaptation. (C) Percentage of noisy genes across treatment stages. (D) Number of clones after awakening in HepG2 cells treated with cisplatin for seven days and KDM5-C70 once after cisplatin removal (E) Confluency of HepG2, HuH6 and HB303 cells treated with cisplatin for two days. Lines represent average of two or three replicates, shaded areas represent standard deviation of the mean. Dashed lines highlight day of cisplatin treatment and removal. Confluency was normalized to the percentage on the day of cisplatin removal. Stars represent statistical significance resulting from a t-test (‘****’ p <0.0001 ‘***’ p <0.001, ‘**’ p<0.01, ‘*’ p 0.05’). We further studied the role of noise in driving recovery from cisplatin in hepatoblastoma. In this model, and similarly to neuroblastoma, plasticity drives the wakening and regrowth of RES-DTP cells post cisplatin treatment, leading to a recovered population that captures the original phenotypic heterogeneity 54 (Figure S10D). We followed the changes in noise in the RES cells during the different stages of treatment to discern whether noise was having a role in this recovery process. Noise quantification showed similar dynamics to those observed in the neuroblastoma drug response framework. The highest levels of TN were observed during awakening, with levels of highly variable genes resetting to pre-treatment values after recovery ( Fig. 9B-C ; Table S6). As we had uncovered in neuroblastoma, we observed a core of 152 noisy genes at baseline, which, during awakening expanded to include thousands of additional genes, primarily involved in cell cycle, differentiation and plasticity regulation. These included genes known to drive plasticity during cisplatin treatment in hepatoblastoma, as well as markers of RES and PROL cells 54 ( Figure 9B and S10E-F). In contrast, HepG2 cells did not display high levels of noise, neither did these levels differ between cell states or treatment phases (Figure S11A-D). We hypothesise that the low levels of TN stem from limited plasticity under these experimental conditions, as we observed a weaker bottleneck on treatment and an incompletely recovered population post-treatment (Figure S11C). Despite the limited noise in HepG2, we do observe an overlap between noisy genes in both hepatoblastoma models as well as the biological processes they are involved in (Figure S11E-F). To test the actionability of targeting noise to prevent plasticity in hepatoblastoma, we treated HuH6, HepG2 and the PDX line HB303 with KDM5-C70 after cisplatin withdrawal. KDM5 inhibition halted awakening in HepG2 and HB303 cells, and moderately delayed this process in HuH6 cells ( Figure 9D-E ), suggesting a common mechanism in neuroblastoma and hepatoblastoma. Taken together, our data suggests that the cell state diversification preceding and fuelling phenotypic evolution in paediatric cancers is driven by a regulated increase in transcriptional noise, specifically targeting genes associated with cell state identity. Whether this mechanism is conserved across other paediatric and adult cancers remains an open question for future investigation. DISCUSSION Childhood cancers, such as neuroblastoma, are characterised by disruptions in developmental processes 2 , 57 , allowing tumours to exploit embryonic mechanisms such as phenotypic plasticity to adapt without requiring additional genetic mutations. However, direct evidence linking plasticity to cancer evolution has been limited due to technological constraints. Here, through the integration of experimental evolution, expressed DNA barcoding, single-cell multi-omics, dynamic manifold modelling, and chromatin profiling, we provide a framework that resolves this long-standing question. We show that plasticity, rather than genetic selection, drives early adaptive evolution in paediatric tumours, allowing them to acquire high-fitness phenotypes without new mutations. A substantial fraction of cancers are thought to evade treatment through non-genetic mechanisms of adaptation, often involving entry into drug-tolerant persister (DTP) states with quiescent features 58 . This process may be particularly relevant to paediatric cancers, whose genomic silence limits the availability of mutational routes for adaptation 4 . In neuroblastoma, previous studies have linked ADRN-to-MES transitions with treatment resistance and invasiveness 21 – 23 . However, whether such transitions arise through selection of pre-existing subpopulations or true phenotypic interconversion has been unclear due to technological constraints. Our lineage tracing and functional analyses show that plasticity is an intrinsic property of neuroblastoma cell populations. This intrinsic plasticity likely serves as an evolutionary bet-hedging strategy, maintaining a small reservoir of cells that – by chance – exist in tolerant, MES-like states before treatment. Upon exposure to chemotherapy, phenotypic selection enriches for these pre-adapted MES-DTPs. Yet, crucially, tumour regrowth after treatment withdrawal is not driven by clonal selection but by plasticity-led awakening, where dormant DTPs stochastically transition back to proliferative ADRN-like states. Similar dynamics of dormancy followed by stochastic reactivation have been observed in adult cancers such as breast cancer, supporting the idea that plasticity-first evolution is a general principle of tumour relapse 59 . Despite plasticity being widely studied in the context of neuroblastoma, the molecular mechanisms enabling such widespread transcriptional rewiring are poorly understood. The seemingly stochastic awakening of MES DTPs after treatment, suggests that plasticity is driven by spontaneous transcriptomic changes rather than coordinated signalling pathways. Transcriptional noise – the stochastic variability in gene expression – has a well-described role in promoting phenotypic diversification in homogeneous populations 29 , 30 , 43 , 60 , 61 . Our data show that transcriptional noise is highest in bridge (BRI) cells, the most plastic and least stable state in neuroblastoma, and that inducing noise experimentally (via IdU or DRB treatment) is sufficient to increase phenotypic transitions and the emergence of BRI cells. In addition, post-treatment stochastic awakening is accompanied by a sharp rise in transcriptional noise, peaking just before population regrowth. This enhances expression variability in cell identity genes, increasing the likelihood that dormant MES-DTPs stochastically transition into proliferative phenotypes once the drug is removed. These results provide direct experimental evidence in clinically relevant models that transcriptional noise drives plasticity-led adaptation, representing a mechanistic bridge between stochastic gene expression and evolutionary fitness. While previous studies proposed that noise facilitates stochastic fate transitions in stem or cancer cells, they have largely relied on mathematical models or limited experimental systems. For example, Hang et al. described that noise-mediated heterogeneity in Sca-I expression influences fate decisions in murine hematopoietic stem cells 60 and Gupta et al. applied Markov modelling to show stochastic state transitions in breast cancer cells 61 . Our current study however directly links noise dynamics to tumour recovery and relapse, establishing TN as a functional determinant of evolutionary potential. Importantly, we extend these observations to hepatoblastoma models, suggesting that noise-mediated plasticity may represent a more generalized mechanism across pediatric cancers. Although the specific noisy genes differed between models, we identified conserved sets of transcriptional patterns: a shared core of noisy genes, which was consistently present across treatment stages, and an awakening-specific set that emerged post-treatment. These gene sets were both associated with the particular cell identity of the model: plasticity (i.e. EMT), cell migration, and neurogenesis in heterogeneous neuroblastoma cell lines; neurogenesis in adrenergic-enriched patient-derived organoids; and liver development and plasticity-related processes in hepatoblastoma. Together, these data support a model in which transcriptional noise affects lineage-determining gene networks to promote cell state transitions during both intrinsic and adaptive plasticity. Our mechanistic analyses identify the histone demethylase KDM5A as a key regulator of noise-induced plasticity. We show that KDM5A binding increases at noisy, cell-identity genes after treatment, coinciding with a narrowing of H3K4me3 domains, a chromatin signature associated with increased transcriptional variability. Pharmacological inhibition of KDM5 with KDM5-C70 prevented this H3K4me3 peak sharpening, reduced transcriptional noise, suppressed phenotypic plasticity, and blocked tumour regrowth in both neuroblastoma and hepatoblastoma models. These findings uncover a direct mechanistic link between chromatin remodelling and transcriptional noise regulation during cancer evolution. This aligns with prior work showing that histone methylation modulates gene-expression variability in embryonic stem cells 62 , and with observations that KDM5 activity enhances transcriptomic heterogeneity and treatment resistance in ER⁺ breast cancer 53 . Together, these studies suggest a conserved role for KDM5-driven noise regulation as a facilitator of plasticity and adaptive evolution across developmental and cancer contexts. Finally, we observed that a substantial subset of noisy genes retained their expression variability across multiple cell generations, implying a form of transcriptional memory. These findings support a model in which initially stochastic gene expression states can become stabilized, allowing noise-induced expression patterns to persist and influence future cell fate decisions. We propose that this stabilization is favored by the architecture of neuroblastoma core regulatory circuitries (CRC), where key transcription factors regulate themselves and each other through multilayered positive feedback loops 63 . This hypothesis is supported by prior studies showing that such feedback loops can stabilize noise-induced high expression of fitness-associated genes in rare subpopulations, thereby sustaining specific cell fate decisions 29 , 30 . Specifically to neuroblastoma, transcriptional noise within the JNK signaling network, amplified by feedback regulation, has been shown to drive resistance to treatment and create a “memory effect” of transient variability 33 . However, whether similar feedback mechanisms within CRCs underpin the heritability and stabilization of transcriptional noise and its consequences in our models remains to be determined. Future work will be essential to dissect how noise-induced transcriptional states are selected, maintained, or erased across therapeutic cycles and how this contributes to cancer evolution. Collectively, our findings demonstrate that neuroblastoma cells leverage plasticity to generate phenotypic diversity in treatment-naïve conditions and to enable awakening from DTP states after therapy. This study provides, to our knowledge, the first experimental evidence directly linking transcriptional noise to phenotypic plasticity and relapse in cancer models. Beyond advancing our mechanistic understanding of plasticity and adaptive evolution, these findings position noise as a novel non-genetic driver of plasticity, a key cancer hallmark. This has important implications for the identification of biomarkers associated with increased adaptive potential and malignant progression. Moreover, our results open new therapeutic avenues by targeting transcriptional noise and plasticity to prevent relapse and limit tumour evolution. Although the clinical application of KDM5 inhibitors has been constrained by toxicity concerns 64 , recent advances in targeted protein degradation technologies – such as proteolysis-targeting chimeras (PROTACs) – offer a promising strategy to overcome these limitations and translate our findings into clinical benefit 65 – 67 . Ultimately, targeting the regulatory mechanisms underpinning transcriptional noise may represent a novel approach to reducing cancer plasticity and improving long-term treatment outcomes. MATERIALS AND METHODS Cell lines and treatment SK-N-SH human neuroblastoma cells were originally provided by L. Chesler, ICR, London. HuH6 human hepatoblastoma cells were provided by R. Kappler, LMU, Germany. HepG2 human hepatoblastoma cells were purchased from ATCC. All cell lines were cultured in DMEM (Thermo Fisher Scientific, USA) supplemented with 10% fetal bovine serum (FBS) (Thermo Fisher Scientific, USA), 100 IU penicillin and 100μg/mL streptomycin (Sigma-Aldrich, USA) at 37°C and 5% CO 2 . Unless specified otherwise, cells were treated with 50µM IdU (I7125, Sigma-Aldrich, USA), 100nM TGF-ß (240-B-010, Biotechne, USA), 30µM KDM5-C70 (HY-120400, MedChem Express, USA) or an IC50 concentration of cisplatin (CAY13119, Cambridge Bioscience, UK): 3.18µM for SK-N-SH, 11.1µM for HuH6 and 14.1µM HepG2. To achieve transient pausing of transcription, cells were treated with 200µM DRB (Sigma-Aldrich, USA), which was removed after 6 hours to promote normal transcription. Wound healing assay SK-N-SH cells were plated in a tissue-culture 48 well-plate (Thermo Fisher Scientific, USA) at a density of 20,000 cells per cm 2 . After 24 hours, a scratch was done with a pipette tip and wells were washed twice with PBS to remove detached cells and provided with fresh media containing the desired drug. DMSO (D2650, Sigma-Aldrich, USA) was added at a final concentration of 0.25% v/v. IdU (I7125, Sigma-Aldrich, USA) was titrated at 10, 20 and 50µM; if not stated otherwise, the concentration used was 50µM. TGF-ß (240-B-010, Biotechne, USA) was supplied at a final concentration of 100nM. Cell migration was monitored using the IncuCyte S3 instrument (Sartorius, Germany) acquiring images every hour at 10X magnification. Measurements of wound area were done using an ImageJ plug-in 68 and normalized to wound area of each condition at 0 hours. Flow cytometry and sorting Treated cells were detached using TrypLE Express (Thermo Fisher Scientific, USA) and collected in 5ml polystyrene FACS tubes. Cells were washed in in 1XPBS + 1% Bovine Serum Albumin (Thermo Fisher Scientific, USA), here FB buffer, pelleted at 400xg for 5 minutes at 4°C and resuspended in 100μL FB. For experiments not including KDM5-C70, cells were stained with 2μL of anti-CD44 APC/Fire 750 antibody (103062, BioLegend, USA) for 30 min at 4°C. Cells were washed by centrifugation in FB and resuspended in 300μL FB containing 25ng/μL DAPI (130-111-570, Miltenyi, Germany). In experiments including KDM5-C70, washed cells were first stained with 1μL Zombie Red viability dye (423109, BioLegend, USA) for 10 min at room temperature and 2μL of anti-CD44 APC/Fire 750 antibody for 30 min at 4°C. Cells were then washed by centrifugation in FB and resuspended in 300μL FB. Cells were analyzed using a BD LSR II or Symphony A5 (BD Bioscience, USA) or sorted using a Symphony S6 (BD Bioscience, USA). Data was analyzed using the FlowJo software. Colony forming assay For the colony forming assay cells were seeded at 1000 cells/well of a 6-well plate in complete media. Colony growth was monitored, and media was replenished once a week. Plates were imaged using the Incucyte S3 system (Sartorius), and images were analysed using the confluence tool, which generates the yellow segmentation mask visible in the output images. Segmentation parameters, including hole filling and area thresholds, were optimised per experiment and then held constant throughout the analysis. For colony number quantification, the Analyze Particles tool in Fiji (ImageJ2) was used. Immunofluorescence For immunofluorescence staining, 2× 10^5 cells were seeded in a 24-well plate containing glass coverslips. The day after, cells were treated with Cisplatin for the indicated amount of time. After treatment, cells were fixed with 4% PFA for 30 minutes at 4°C. Cells were then permeabilized with 0.2% Triton X-100 for 10 minutes at RT. Cells were washed twice and incubated in PBS buffer containing 1% BSA and 2% FCS (IFF buffer) for 30 minutes at room temperature (RT). Cells were stained overnight at 4 °C with primary antibody diluted in IFF buffer (anti-CD44 IM7, 1:500, eBioscience; anti-gH2A.XpPhospho (Ser139) 2F3, 1:200, Biolegend). The following day, cells were washed three times with PBS and subjected to conjugated secondary antibody diluted in IFF buffer (anti-Mouse-AF488 A21202 Invitrogen, anti-Rat-AF647 A78947 Invitrogen) at 1:2000 for 2h at RT. Coverslips were mounted with DAPI Fluoromount-G mounting medium and imaged the with a confocal microscope (Zeiss LSM700). Barcoding of SK-N-SH cells The CloneTracker library was a gift from Andrea Sottoriva (BCXP1M3RP-XSP, Cellecta). SK-N-SH cell lines were cultured in normal growth media and barcoded by lentiviral transduction according to manufacturer’s recommendations. Following lentiviral titration, 1 million neuroblastoma cells were transduced with a lentiviral library using a multiplicity of infection (MOI) of 0.1 (corresponding to 10% of cells transduced) to attain a single barcode per cell. 100,000 barcoded cells were FACS-sorted based on RFP-positivity and expanded. This population constitutes the “original population” or POT. Experimental design of single cell RNA sequencing of barcoded SK-N-SH cell line Barcoded SK-N-SH cells were seeded equally into 2 treatment groups and divided across three replicates. The seeding density was adjusted based on each treatment and a total of 5M cells were seeded for each sample. Untreated cells were collected at confluency. The remaining three samples were exposed to cisplatin treatment (3.18 uM). After 1 week (T1), the treatment was removed and cells were detached, 10% of the total population was cryopreserved in four aliquots of 2.5% and the remaining cells were reseeded at a fixed density and grown under normal growth media for 1 week (T2) or 4 weeks (T3). Cell counts were determined via the Countess II Automatic Cell Counter (Thermo Fisher). Single cell RNA profiling of barcoded SK-N-SH cell line After thawing, single cells were washed with PBS and were resuspended with a buffer (Ca++/Mg++ free PBS + 0.04% BSA). Cells were tagged with the CMO oligos (10X Genomics) as per manufacturer’s instructions. Samples were multiplexed in groups of 6, each containing 5000 cells/samples for a total of 30000 cells/pool. Pooled single cell suspensions were loaded on a Chromium Single Cell 3′ Chip (10X Genomics) and were run in the Chromium Controller to generate single-cell gel bead-in-emulsions using the 10X genomics 3′ Chromium v2.0 platform as per manufacturer’s instructions. Single-cell RNA-seq and CMO libraries were prepared according to the manufacturer’s protocol and the library quality was confirmed with a Bioanalyzer High-Sensitivity DNA Kit (Agilent, 5067-4627). Cellecta barcode libraries were prepared using 25% of the adaptor-ligated cDNA to amplify only the barcoded cDNA product. Barcodes were pre-amplified using the P5-Read1 primer (ACACTCTTTCCCTACACGACGCTCTTCCGATCT) and the P7-Adapter-FBPI Cellecta primer (GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCTCCGACCACCGAACGCAACGC ACGCA), following which Illumina indexes were added. Samples were sequenced on an Illumina NovaSeq 6000 according to standard 10X Genomics protocol, aiming for approximately 50,000 reads per cell for gene expression libraries and 1000 reads per cell for CMO/Cellecta barcode libraries. Single cell data pre-processing 10x Genomics CellRanger (v6.1.2) was used to process raw sequencing FASTQs, using cellranger count (or multi in the case of samples multiplexed with CellPlex) pipeline to conduct read alignment to human genome reference (GRCh38-2024-A), sample demultiplexing and generate feature-barcode matrix. For PDX models, read alignment was conducted using combined human and mouse genome reference (GRCh38_and_GRCm39-2024-A), to observe any mouse gene expression contamination within the sample. EmptyDrops method from DropletUtils (v1.20.0) package 69 , 70 was implemented on cellranger output H5 matrix to remove empty droplets with false discovery rate of 0.1%. Data files were merged and then normalised using logNormCounts from scuttle (v1.10.1) package 54 for downstream QC and processing. Cell-level filtering was carried out with the following thresholds: mitochondrial gene content 5%, and UMI counts >1000 for cell lines, mitochondrial gene content 0.3%, and UMI counts > 1000 for PDO models, mitochondrial gene content 0.2%, and UMI counts > 1000 for PDX models. Filtering at the gene-level was implemented to remove genes with fewer than five reads. Downstream normalisation and scaling were performed using Seurat (v4.3.0) 71 with the SCTransform function, using the top 1000 variable genes, and regressing out cell cycle and mitochondrial genes. The SCTransform data output was used to perform a principal component analysis (PCA). The batch correction algorithm Harmony (v0.1.1) was applied across CMO multiplexed pools (where appropriate), and the output used for Uniform Manifold Approximation and Projection (UMAP)-based dimensional reduction. Clustering was carried out using the Seurat FindClusters function and gene expression markers for the clusters were identified by a Wilcoxon Rank Sum test, using the Seurat FindAllMarkers function. Gene set enrichment analysis using clusterProfiler (v4.6.2) 72 , 73 was then carried out on the cluster markers. To compare changes in cell clusters, statistical testing was conducted using student t-test between each sample and T0 (untreated), with adjusted p-value calculated using Bonferroni method for multiple-comparison correction. Cell state score generation Single-cell expression signatures for the ADRN and MES cell states were created using the AddModuleScore function in Seurat and published gene sets 21 . To calculate the cell state score we used a method based on that used by Tan et al. 74 . Statistical testing (t-tests) was used to determine whether differences between ADRN and MES gene set expression were significant. The t-statistic was taken as the cell state score and the statistical significance was further used to segregate cells into phenotypic groups: ADRN (AMT score < 0; P < 0.05), intermediate ADRN (AMT score 0; P ≥ 0.05) and MES (AMT score > 0; P < 0.05). The two intermediate classifications were grouped together as INT. To compare changes in phenotypic cell state assignment, statistical testing was conducted using student t-test with adjusted p-value calculated using Bonferroni method for multiple-comparison correction. Cellecta barcodes amplification and next generation library preparation SK-N-SH cells were harvested and pelleted at the relevant timepoints. Genomic DNA was isolated using AllPrep kit (Qiagen) according to manufacturer’s instructions and quantified using Qubit (Life Technologies). Amplicon PCR was conducted using 50ng of DNA (or 20ng for samples with lower DNA yield) along with 2x Accuzyme mix (Bioline) to amplify barcodes using the following primer sequences (Forward - ACTGACTGCAGTCTGAGTCTGACAG, Reverse - CTAGCATAGAGTGCGTAGCTCTGCT). After identifying and purifying PCR product (consisting of 80-bp product + 30-bp semi-random barcode), NGS libraries were prepared using NEBnext Ultra II DNA library preparation kit for Illumina (New England Biolabs) and quantified with using a TapeStation (Agilent Genomics). Samples were balanced and pooled together to have final molarity of 4nM, NGS was then performed by ICR Genomics facility (TPU) using NovaSeq SP PE150. Cellecta barcodes bioinformatics analysis FASTQ files from bulk experiments were mapped to a custom index reference generated from white-list (all possible sequences) of Cellecta barcodes, using Burrows-Wheeler aligner (BWA-MEM tool; version 0.7.17) 75 . Once mapped BAM files were generated using Sequence Alignmnet/Map (SamTools; v1.18) 76 . Gene counts were quantified using FeatureCounts from Rsubread. Parameters were set to filter out reads with mapping quality less than 30 (minMQS = 30). To observe lineage dynamics, representative percentage was calculated for each lineage within each sample (frequency of Cellecta barcode / frequency of all Cellecta barcodes × 100). Shannon diversity index 77 of all DNA barcodes from each experimental samples was calculated between conditions to observe changes in lineage diversity as per following equation: H=−∑[(pi)×log(pi)], where pi is proportion of individual lineage within the population. Statistical significance of results was determined using pairwise t-tests . Cumulative proportion quantification of all DNA barcodes from each sample was calculated by ranking all barcodes from most to least abundant (sum of frequency). Cumulative proportion was then calculated as cumulative sum of barcode abundance / sum of abundance. For statistical quantification an area under the curve (AUC) value was calculated for each sample using trapz() (from R package prcama) of cumulative proportion and abundance from cumulative proportion plot. Statistical significance of results was determined using pairwise t-tests between untreated and each other samples. For top 100 POT barcode analysis, barcodes in POT sample were ranked most to least frequent and the top 100 most abundant barcodes were extracted. Remaining experimental samples were filtered for these barcodes and visualised. For top 100 recovery barcode analysis, barcodes in each T3 sample were ranked most to least frequent and the top 100 most abundant barcodes were extracted. Samples from the same experimental treatment arm were filtered for these barcodes and visualised. Common colour scheme for barcodes of interest at T3 was generated across all top 100 barcodes lists from all experimental arms, to allow visualisation of common barcodes across experimental arms. For scRNA-seq experiments Cellecta barcodes from 10x Chromium libraries were isolated using STAR solo 78 with the same white list of barcodes used for BWA for bulk experiments. Where multiple barcodes were observed within a single cell, cells were assigned the barcode with the higher proportion of representation from read counts. To account for sparse data, lineages were filtered to ensure each barcodes were observed at least 3 times within the entire dataset. For plasticity score analysis, each lineage was given a plasticity score based on the number of cell states occupied by that lineage. Statistical significance of results was determined using pairwise t-tests between untreated and each other samples. To assess the relationship between cell frequency and plasticity group assignment we conducted spearman’s correlation coefficient in each condition. P-value represented on each plot corresponds to correlation coefficient significance across all replicates within each condition. Plasticity transition dynamics analysis was conducted by assigning each cell within each lineage a plasticity score based on the number of cell states occupied; either =1 for a single phenotypic state or > 1 for multiple phenotypic states. For analysis frequency of cells in each group were aggregated across replicates. Chi-squared statistical testing was conducted in a pairwise fashion to observe whether there was an association between sample and plasticity based on frequency of each plasticity group. Statistical test was conducted between untreated and each other sample and subject to multiple comparison correction using Bonferroni’s method. Directionality of statistical significance was calculated as (frequency of cells in >1 plasticity group / total cell frequency in sample) - (frequency of cells in >1 plasticity group from untreated / total cell frequency in untreated). Positive was termed change value > 0 and negative < 0. To identify and track cellular lineages that successfully persisted throughout the treatment, we implemented a matched clone analysis focusing on clones present in both untreated and recovery conditions. We identified winner clones by filtering criteria. Only clones with a minimum of 3 cells in both the untreated condition and the 4-week off-treatment timepoint were included in the analysis. For each matched clone, we tracked their AMT state transitions across all time points. The dominant cellular state at each timepoint was determined by the most frequent annotation among cells within each clone. To maintain consistency and statistical reliability, the minimum cell threshold of 3 cells per condition was applied across all timepoints. Clone trajectories were visualized using alluvial flow diagrams implemented with the ggalluvial R package, enabling tracking of individual clone over time. Patient-derived xenograft establishment PDX (IC-pPDX-112 and GR-NB5) were established from neuroblastoma tumor biopsies obtained at diagnosis, as described previously, following informed consent of parents or guardians 79 , 80 . This study was performed in accordance with the recommendations of the European Community (2010/63/UE) for the care and use of laboratory animals. Experimental procedures were approved by the ethics committee of the Institut Curie CEEA-IC #118 (Authorization APAFIS#11206-2017090816044613-v2) in compliance with the international guidelines. The establishment of PDX received approval by the Institut Curie institutional review board OBS170323 CPP ref 3272; n° de dossier 2015-A00464-45. In brief, tumor samples were engrafted subcutaneously in immunocompromised mice. Upon successful establishment (> P2) and upon reaching ethical size, PDX tumors were cryopreserved for banking, morphological and molecular analyses according to procedures described previously. Spheroid morphology quantification IC-pPDX-112 was dissociated as previously described ( 63 ) and cultured in DMEM supplemented with 10% FBS. 80,000 cells/well were seeded in a tissue-culture flat-bottom 96-well plate (Thermo Fisher Scientific, USA). After 10 days in culture, cells were treated with DMSO (0.25% v/v) or IdU (50µM) in the presence or absence of TGF-ß (100nM) and monitored using the IncuCyte ZOOM instrument (Sartorius, Germany) at a 10X magnification every hour. To analyze their morphology, spheroids were manually segmented, and area, roundness and total number of fusion events were calculated using ImageJ. RNA sequencing of PDX IC-pPDX-112 was dissociated as previously described 81 and cultured in DMEM supplemented with 10% FBS. 1×106 cells/well were seeded in tissue-culture 6-well plates. After 10 days in culture, cells were treated with DMSO (0.25% v/v) or IdU (50µM) in the presence or absence of TGF-ß (100nM). 5 days after treatment cells were collected and RNA was isolated using the RNeasy Micro Kit (Qiagen, Germany). RNA quality was assessed using the Bioanalyzer System (Agilent, USA) and 600ng of RNA were used for library preparation with the Illumina Stranded mRNA Prep, Ligation kit (Illumina, USA). Sequencing was performed in a NovaSeq instrument (Illumina, USA) using 30 million clusters per sample and PE100 reads. Pre-processing of the raw sequencing reads was performed using an in-house pipeline available via https://zenodo.org/records/13744441 . Briefly, reads were trimmed using TrimGalore v0.6.2 to remove adapters and quality-checked using FastQC v.0.11.9. Reads aligning to the mm10 murine reference genome were removed using xengsort ( 64 ) to avoid contamination from host mouse cells. Remaining reads were aligned to the GRCh37/hg19 human reference genome using STAR aligner v2.7.6 and gene-level counts were quantified with featureCounts using GENCODE annotation. Differential gene expression analysis was performed using DESeq2 (v1.40.2) in R (v4.3.0). Patient-derived models and treatment for scRNA-seq analysis The PDX line GR-NB5 (MYCN AMP , previously named MAPGR-B25-NB-117 82 ) was originally provided by G. Schleiermacher (Curie Institute, Paris, France). Cells were obtained from enzymatical dissociation of the PDX and cells were cultures in in low-Glucose GlutaMAX DMEM (Thermo Fisher) supplemented with 25% Ham’s F-12 Nutrient Mix, 2% FBS (Thermo Fisher), B-27 Supplement (50X) (Thermo Fisher), minus vitamin A N-2 Supplement (100X) (Thermo Fisher), 100 U/mL penicillin-streptomycin (Sigma), 20 ng/mL Animal-Free Recombinant Human EGF (Peprotech), 40 ng/mL Recombinant Human FGF-basic (Peprotech), 200 ng/mL Recombinant Human IGF-I (Peprotech), 10 ng/mL Recombinant Human PDGF-AA (Peprotech) and 10 ng/mL Recombinant Human PDGF-BB (Peprotech). NB039 (MYCN AMP ALK WT ) and NB067 (MYCN AMP ALK R1275Q ) patient-derived spheroids were originally provided by J. Molenaar (Princess Maxima Centrum, Netherlands). PDOs were cultured in low-Glucose GlutaMAX DMEM (Thermo Fisher) supplemented with 25% Ham’s F-12 Nutrient Mix, B-27 Supplement (50X) (Thermo Fisher), minus vitamin A N-2 Supplement (100X) (Thermo Fisher), 100 U/mL penicillin-streptomycin (Sigma), 20 ng/mL Animal-Free Recombinant Human EGF (Peprotech), 40 ng/mL Recombinant Human FGF-basic (Peprotech), 200 ng/mL Recombinant Human IGF-I (Peprotech), 10 ng/mL Recombinant Human PDGF-AA (Peprotech) and 10 ng/mL Recombinant Human PDGF-BB (Peprotech). After thawing, single cells were washed with PBS and were loaded onto a Chromium Single Cell 3′ Chip (10X Genomics) and were run in the Chromium Controller to generate single-cell gel bead-in-emulsions using the 10X genomics 3′ Chromium v2.0 platform as per manufacturer’s instructions. Single-cell RNA-seq libraries were prepared according to the manufacturer’s protocol and the library quality was confirmed with a Tapestation High-Sensitivity DNA Kit (Agilent, 5067-5585). Samples were sequenced on an Illumina NovaSeq 6000 according to standard 10X Genomics protocol, aiming for approximately 50,000 reads per cell for gene expression libraries and 1000 reads per cell for CMO/Cellecta barcode libraries. MuTrans analysis Pre-processing of the scRNA-seq data, split by drug condition and experimental timepoint, was carried out analogous to as outlined above using Scanpy (v1.9.3) 83 . Using the filtered data as inputs, the count matrices were normalised and log transformed with a pseudo-count of 1. The number of transcripts per cell and the mitochondrial content were regressed out and the data was scaled, prior to PCA. Dimensionality reduction using UMAP and Leiden clustering was carried out. Cluster markers were then found using Wilcoxon rank-sum tests. A neighbourhood graph was computed for the processed data using scanpy.pp.neighbors (metric = “cosine”). Construction of cell-fate dynamical manifolds and downstream analysis of cell-state transition was then performed using the MUltiscale method for TRANSient cells 84 . BRI Cell Signature Differential expression analysis was used to identify genes that were upregulated in INT cells compared to ADRN and MES cells in the single cell data (FindMarkers in Seurat package). Genes were considered differentially expressed if adjusted p-value was less than 0.05 and log fold change > 1. The robustness of the signature was validated across different datasets, including patient-derived sample and public datasets 85 (Great Ormand Street Hospital with 5 patients and Princess Maxima Center with 16 patients). The expression of the BRI signature was assessed by AddModulScore function from Seurat package. This function calculates a module score for each cell based on the expression levels of signature genes. Survival analysis In the analysis of a large bulk RNA-sequencing dataset comprising 498 primary neuroblastoma patient samples (GSE62564/GSE49711), the enrichment score of BRI signature genes were calculated using Gene Set Variation Analysis (GSVA) 86 with the single sample Gene Set Enrichment Analysis (ssGSEA). Based on the enrichment score, samples were categorised in to “high” and “low” BRI groups. Kaplan-Meier survival analysis was then conducted and the log-rank test was used to determine the correlation between high BRI scores and overall survival (OS). Transcriptional noise quantification Single-cell RNA-seq data on HuH6 and HepG2 was obtained from https://dataview.ncbi.nlm.nih.gov/object/PRJNA1247955?reviewer=g4kgeej9of8fb9sf5ng8hdepnl and pre-processed as previously describede. Processed expression matrices were then analyzed in R studio (v4.3.0) using the BASiCS package (v2.12.3) 45 , 46 . For intrinsic plasticity analysis, untreated cells were split into equally sized groups based on cell state. Group size was determined by the smallest cell state. For analysis of plasticity during cisplatin recovery, mesenchymal/reservoir cells were split into equally sized groups based on timepoint. Markov chain Monte Carlo simulations were performed using BASiCS_MCMC function to estimate expression parameters (N = 4000, Thin = 10, Burn = 2000). Chain convergence was visually inspected. For differential noise analysis between conditions BASiCS_TestDE function was used with default parameters. Identification of highly variable genes (HVG) within specific conditions was performed using BASiCS_DetectHVG function with a VarThreshold of 0.6. MemorySeq SK-N-SH cells were plated at a density of 0.5 cells/well in 96-well plates to achieve single-cell seeding. Cells were cultured in 50 μL DMEM supplemented with 10% FBS and 50 μL conditioned media from parental SK-N-SH. The cultures were maintained for several weeks and expanded to 24-well and 6-well plates. Upon reaching confluency in a 6-well plate, cells were harvested, and 100,000 cells were flash-frozen for RNA-seq. A total of 27 MemorySeq clones were generated following this methodology. Additionally, 48 technical noise controls were generated by sampling 100,000 cells from a parental SK-N-SH population. RNA extractions were performed using AllPrep Micro (Qiagen, Germany). Library preparations were performed by Azenta Life Sciences using NEBNext Ultra II RNA Library Prep Kit (New England Biolabs, USA). Sequencing was performed in an Illumina NovaSeq X plus instrument (Illumina, USA) using a 2×150bp Paired End (PE) configuration. Raw sequencing reads were quality controlled to remove low-quality bases, then aligned to the human GRCh38/hg38 reference genome. Gene-level counts were generated with featureCounts. Raw counts were converted to Transcripts Per Million (TPM) values to normalize for sequencing depth and gene length. Low-quality samples with fewer than 500,000 mapped reads were excluded. Genes were filtered based on median TPM and gene length thresholds to reduce technical noise, following the criteria previously described 51 . Batch effects were identified through Principal Component Analysis (PCA) and corrected using limma’s removeBatchEffect function (v3.58.1). To assess expression variability, the coefficient of variation (CV), mean expression, skewness, and kurtosis were calculated for each gene across control and sorted samples. To identify heritable gene expression patterns, a Poisson regression model relating CV to mean expression was fitted and highly expressed genes (mean TPM > 1.5) with residuals above the 98th percentile were selected, indicating significantly higher variability than expected based on their expression level. All analysis was performed using R v4.3.2. Generation of d 2 GFP cell lines To construct the destabilized GFP (d 2 GFP) reporter, the d 2 GFP gene was amplified by PCR from the Lentiviral-TOP-dGFP-reporter plasmid (Addgene plasmid #14715) and cloned downstream of an EF1a promoter in the lentiviral plasmid pLX307 Luciferase (Addgene plasmid #117734) using the NEBuilder HiFi DNA Assembly Cloning kit (E5520S, New England Biolabs, USA). Lentiviral-TOP-dGFP-reporter was a gift from Tannishtha Reya and pLX307 Luciferase was a gift from William Hahn. Lentiviral particles were produced in HEK293T cells using psPAX2 and pMD2.G as packaging vectors and diluted 1:1 with fresh DMEM to infect SK-N-SH cells. Transduced cells were selected with 1µg/mL puromycin (P33020, Melford Labs, UK) until complete eradication of control untransduced cells was achieved. The remaining d 2 GFP-positive cells were then sorted using a Symphony S6 instrument (BD Biosciences, USA) to increase purity. Drug screen SK-N-SH d 2 GFP were seeded in duplicated in a 6-well plate at a density of 10,000 cells per cm 2 . After 24h, cells were treated with the appropriate concentration of each drug (Table S1). Five days after treatment, the expression of CD44 and d 2 GFP was assessed by flow cytometry. For analysis, 2000 CD44 high cells were sampled and d 2 GFP expression in each cell was normalized to the average expression across the sample. KDM5-C70 treatment during recovery from cisplatin For SK-N-SH, control untreated cells were plated at a seeding density of 10,000 cells per cm 2 , while treated cells were seeded at 20,000 cells per cm 2 in a 6-well plate in triplicate. After 24h, cells were treated with cisplatin (3.18μM) for seven days, after which cisplatin was washed and KDM5-C70 (30μM) was added only once or weekly for a month. The phenotypic composition of the population was assessed by CD44 flow cytometry at 0, 7, 14 and 35 days post cisplatin treatment. For live imaging and growth assessment in hepatoblastoma models, 50,000 cells were seeded in duplicate or triplicates in a 24-well plate. After 24h, cells were treated with an IC50 dose of cisplatin for 48h. After that period, cisplatin was washed, KDM5-C70 (30μM) or DMSO was added once and cells were followed until normal growth was achieved. Cell growth was monitored and quantified using the IncuCyte S3 system (Sartorius, Germany). For the colony-forming assay 200,000 HepG2 cells were seeded in duplicates in a 6-well plate. After 24h, cells were treated with cisplatin (14.1μM) for seven days. After that period, cisplatin was washed, KDM5-C70 (30μM) or DMSO was added once and cells were kept in culture with weekly media changes until day 35 post cisplatin treatment. On day 35, cells were imaged with the IncuCyte S3 system to quantify colony diameter. After this, media was aspirated, cells were washed with 1xPBS, stained with GIEMSA dye (Sigma-Aldrich, USA) and imaged to quantify number of colonies. RT-qPCR Between 1×10 5 and 4×10 5 cells were harvested and their RNA isolated using the RNeasy Micro Kit (74004, QIAGEN, Germany). 1µg of RNA was used to synthesize cDNA using the RevertAid First Strand cDNA Synthesis Kit (K1622, Thermo Fisher Scientific, USA). qPCR was performed using PoweUp SYBR Master Mix (A25780, Thermo Fisher Scientific, USA) in a C100 Touch Thermal Cycler supplied with a CFX96 Real-Time System (Bio-Rad, USA). Gene expression was calculated with the ΔΔCt method using as internal control the average expression of two housekeeping genes, ACTB and RPS9. Western blotting 1×10 6 cells were harvested and lysed using RIPA buffer (Thermo Fisher Scientific, USA) supplemented with proteinase and phosphatase inhibitors. Lysates were centrifuged at 13,000xg for 20 min to isolate the proteome. Protein lysates were quantified using the DC Protein Assay kit (Bio-Rad, USA). 30µg of protein were denaturalized by incubating at 95°C for 5min in the presence of 2-mercaptoethanol. Samples were loaded in a Bolt 4-12% Bis-Tris Plus Gel (Thermo Fisher Scientific, USA), run using MOPS buffer and transferred into a PVDP membrane. Membranes were blotted using anti-KDM5A (Jarid1A/RBBP2 (EPR18651) Rabbit mAb, ab194286, Abcam), anti-KDM5B (JARID1B (E2X6N) Rabbit mAb, 15327, Cell Signalling Technology) or anti-GAPDH (GAPDH (14C10) Rabbit mAb, 2118, Cell Signalling Technology) primary antibodies and anti-Rabbit IgG-HRP (Goat Anti-Rabbit Immunoglobulins/HRP, P044801-2, Agilent Technologies, UK) secondary antibody at the recommended concentrations. Dose response curve Cell lines were seeded at 2000 cells/well in 384-well plates and cultured for 24 hours prior to dosing. Cells were dosed using a 7-dose response curved with three technical replicates per concentration. Cell viability was assessed using CellTiterGlo-3D (Promega, USA) at day 3 after dosing. Cell viability was normalized to background and solvent controls, and dose response curves were fitted using isotonic regression, as previously described 87 . CUT&Tag sample preparation and analysis SK-N-SH cells were seeded in duplicate in a 6-well plate and treated with an IC50 concentration of cisplatin or left untreated. After one week, untreated and treated cells were collected and viably frozen (“Untreated” and “Persistence” timepoints). For the remaining cells cisplatin treatment was removed. One arm was left to recover untreated and a second arm was treated with 30μM KDM5-C70. Cells were collected and frozen a week later (“Awakening” and “Awakening + KDM5-C70”, respectively). Samples were processed using the Cut & Tag-IT Assay Kit (Active Motif, USA) as per manufacturer’s instructions with 1μg of antibody: anti-KDM5A (Jarid1A/RBBP2 (EPR18651) Rabbit mAb, ab194286, Abcam), anti-KDM5B (JARID1B (E2X6N) Rabbit mAb, 15327, Cell Signalling Technology), anti-H3K4me3 (Tri-Methyl-Histone H3 (Lys4) (C42D8) Rabbit mAb, 9751S, Cell Signalling Technology). Sequencing was performed in a NovaX instrument using 50 paired-end cycles. CUT&Tag sequencing data were processed following the computational pipeline described by Zheng et al. 88 , with modifications for time series analysis. Raw paired-end sequencing reads were aligned to the human reference genome (GRCh38) using Bowtie2 (v2.4.2) with CUT&Tag-optimized parameters (--local --very-sensitive --no-unal --no-mixed --no-discordant -I 10 -X 700). Aligned reads were coordinate-sorted and indexed using SAMtools (v1.11). Quality control metrics including alignment rates, duplication rates, and library complexity were assessed for all samples using custom scripts adapted from the tutorial. Duplicate reads were removed using Picard MarkDuplicates (v2.23.8) with REMOVE_DUPLICATES=true. Peak calling was performed using MACS2 (v2.2.7.1) in paired-end mode (-f BAMPE) with a q-value threshold of 0.1, using untreated samples as matched controls for each histone mark. Reproducible peaks between biological replicates were identified using Homer mergePeaks with a 7.5 kb merging distance for consensus peak set generation. For quantitative visualization, normalized coverage tracks were generated using bamCoverage with RPGC normalization (effective genome size: 2,913,022,398 bp) and 10 bp binning. Signal profiles and heatmaps were created using deepTools computeMatrix and plotHeatmap, with analyses centred on transcription start sites (TSS) and peak summits using ±3 kb flanking regions. Gene-specific chromatin dynamics were assessed by assigning peaks to the nearest gene within 5 kb using bedtools window to evaluate lineage-specific chromatin remodeling across treatment time points. Data availability The raw single-cell RNA sequencing data used in this study can be accessed through the National Centre for Biotechnology Information - Sequence Read Archive (NCBI - SRA) using the following links: https://dataview.ncbi.nlm.nih.gov/object/PRJNA1071370?reviewer=m1girilboh9soufg9bhoio9d9k ; https://dataview.ncbi.nlm.nih.gov/object/PRJNA1247955?reviewer=g4kgeej9of8fb9sf5ng8hdepnl . The raw bulk RNA sequencing data on the PDX cells is available under restricted access, through the European Genome-Phenome Archive (EGA) with ID EGAD50000001770.The raw bulk RNA sequencing used for MemorySeq analysis can be accessed here: https://dataview.ncbi.nlm.nih.gov/object/PRJNA1256291?reviewer=ockms98q7ckoirf9hj0vmgtnn6 Code availability Code used to process the data in this study is available in the following GitHub repository: https://github.com/PCEBrunaLab . AUTHOR CONTRIBUTIONS Methodology, A.A.A., C.R., S.H., A.S., E.C., C.E., A.S.S., C.B., A.Be.; experimental design, A.A.A., C.R., S.H., A.S., M.D.M., A.B..; analytical investigation, A.A.A., S.H., E.C., A.S.S., H.C., A.B.; writing–original draft, A.A.A., C.R., A.B.; writing–review and editing, A.A.A., C.R., S.H., A.S., E.C., C.E., M.D.M., A.B.; resources, L.C., G.S., B.G.; conceptualisation, A.B.; project administration, A.B.; funding acquisition, G.S., L.C., A.B. ACKNOWLEDGEMENTS We thank Katerina Rekopoulou and Nik Matthews for assistance with sample and library preparation for single-cell RNA-sequencing; Jan Molenaar for providing the patient-derived spheroid lines; Andrea Sottoriva for sharing the Cellecta barcode library and expertise in barcoding the cells as well as Salvatore Milite and Haider Tari for assistance and technical advice in analysis and barcode sequencing processing; Yura Grabovska and Chris Jones for technical expertise in the analysis of scRNA-seq data; Leor Weinberger, Melanie Beckett, Elizabeth R. Tucker, Evon Poon and Karen Barker for helpful discussion and technical advice; Peijie Zhou for help with the MuTrans analysis; Oscar M. Rueda and Maurizio Callari for technical advice regarding WGS and WES analyses. The PDX model was developed within the MAPPYACTS study following informed consent from patients and legal representatives. The authors thank all patients and parents/caregivers who participated in the Mappyacts trial and consented to the development of preclinical models. We acknowledge the ICR facilities: Light Microscopy and Confocal, Flow Cytometry, Genomics and Scientific Computing. We acknowledge the RSE Group at The Institute of Cancer Research for providing software enhancements, specifically Rachel Alcraft. Their support was instrumental in achieving the results presented in this paper. A.B. is supported by ICR London. This study was supported by Little Princess Trust Innovation Grant LPTINS2\9 (to A.B. and L.C.) and VAGABOND Marie Sklodowska-Curie Actions ITN (Horizon 2020) Grant 956285 (A.B., L.C. and G.S.). Funder Information Declared The Little Princess Trust , LPTINS2\9 Horizon 2020 , 956285 REFERENCES 1. ↵ Stiller , C. , Irvine , L. & Welham , C. Children, Teenagers and Young Adults UK Cancer Statistics Report 2021 . https://www.qub.ac.uk/research-centres/nicr/FileStore/PDF/UKReports/Filetoupload,1024819,en.pdf ( 2021 ). 2. ↵ Marshall , G. M. et al. The prenatal origins of cancer . Nat Rev Cancer 14 , 277 – 289 ( 2014 ). OpenUrl CrossRef PubMed Web of Science 3. ↵ Ma , X. et al. Pan-cancer genome and transcriptome analyses of 1,699 paediatric leukaemias and solid tumours . Nature 555 , 371 – 376 ( 2018 ). OpenUrl CrossRef PubMed 4. ↵ Gröbner , S. N. et al. The landscape of genomic alterations across childhood cancers . Nature 555 , 321 – 327 ( 2018 ). OpenUrl CrossRef PubMed 5. ↵ Gore , L. , DeGregori , J. & Porter , C. C . Targeting developmental pathways in children with cancer: what price success? The Lancet Oncology 14 , e70 – e78 ( 2013 ). OpenUrl PubMed 6. Moore , S. W . Developmental genes and cancer in children . Pediatr Blood Cancer 52 , 755 – 760 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 7. Shah , N. & Sukumar , S . The Hox genes and their roles in oncogenesis . Nat Rev Cancer 10 , 361 – 371 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 8. Pulford , K. et al. The emerging normal and disease-related roles of anaplastic lymphoma kinase . Cell Mol Life Sci 61 , 2939 – 2953 ( 2004 ). OpenUrl CrossRef PubMed Web of Science 9. ↵ Filbin , M. & Monje , M . Developmental origins and emerging therapeutic opportunities for childhood cancer . Nat Med 25 , 367 – 376 ( 2019 ). OpenUrl CrossRef PubMed 10. ↵ Maris , J. M . Recent advances in neuroblastoma . N Engl J Med 362 , 2202 – 2211 ( 2010 ). OpenUrl CrossRef PubMed Web of Science 11. ↵ Basta , N. O. et al. Factors associated with recurrence and survival length following relapse in patients with neuroblastoma . Br J Cancer 115 , 1048 – 1057 ( 2016 ). OpenUrl CrossRef PubMed 12. ↵ Pugh , T. J. et al. The genetic landscape of high-risk neuroblastoma . Nat Genet 45 , 279 – 284 ( 2013 ). OpenUrl CrossRef PubMed 13. Bosse , K. R. & Maris , J. M . Advances in the translational genomics of neuroblastoma: From improving risk stratification and revealing novel biology to identifying actionable genomic alterations . Cancer 122 , 20 – 33 ( 2016 ). OpenUrl CrossRef PubMed 14. ↵ Schramm , A. et al. Mutational dynamics between primary and relapse neuroblastomas . Nat Genet 47 , 872 – 877 ( 2015 ). OpenUrl CrossRef PubMed 15. ↵ Gundem , G. et al. Clonal evolution during metastatic spread in high-risk neuroblastoma . Nat Genet 55 , 1022 – 1033 ( 2023 ). OpenUrl CrossRef PubMed 16. ↵ Jansky , S. et al. Single-cell transcriptomic analyses provide insights into the developmental origins of neuroblastoma . Nat Genet 53 , 683 – 693 ( 2021 ). OpenUrl CrossRef PubMed 17. Villalard , B. et al. Neuroblastoma plasticity during metastatic progression stems from the dynamics of an early sympathetic transcriptomic trajectory . Nat Commun 15 , 9570 ( 2024 ). OpenUrl CrossRef PubMed 18. ↵ Bedoya-Reina , O. C. et al. Single-nuclei transcriptomes from human adrenal gland reveal distinct cellular identities of low and high-risk neuroblastoma tumors . Nat Commun 12 , 5309 ( 2021 ). OpenUrl CrossRef PubMed 19. ↵ Kildisiute , G. et al. Tumor to normal single-cell mRNA comparisons reveal a pan-neuroblastoma cancer cell . Sci Adv 7 , eabd3311 ( 2021 ). OpenUrl FREE Full Text 20. ↵ Ross , R. A. , Spengler , B. A. & Biedler , J. L . Coordinate Morphological and Biochemical Interconversion of Human Neuroblastoma Cells . JNCI: Journal of the National Cancer Institute 71 , 741 – 747 ( 1983 ). OpenUrl PubMed Web of Science 21. ↵ van Groningen , T. et al. Neuroblastoma is composed of two super-enhancer-associated differentiation states . Nat Genet 49 , 1261 – 1266 ( 2017 ). OpenUrl CrossRef PubMed 22. ↵ Boeva , V. et al. Heterogeneity of neuroblastoma cell identity defined by transcriptional circuitries . Nat Genet 49 , 1408 – 1413 ( 2017 ). OpenUrl CrossRef PubMed 23. ↵ van Groningen , T. et al. A NOTCH feed-forward loop drives reprogramming from adrenergic to mesenchymal state in neuroblastoma . Nat Commun 10 , 1530 ( 2019 ). OpenUrl CrossRef PubMed 24. ↵ Thirant , C. et al. Reversible transitions between noradrenergic and mesenchymal tumor identities define cell plasticity in neuroblastoma . Nat Commun 14 , 2575 ( 2023 ). OpenUrl CrossRef PubMed 25. ↵ Dong , R. et al. Single-Cell Characterization of Malignant Phenotypes and Developmental Trajectories of Adrenal Neuroblastoma . Cancer Cell 38 , 716 – 733 .e6 ( 2020 ). OpenUrl CrossRef PubMed 26. ↵ Jansky , S. et al. Single-cell transcriptomic analyses provide insights into the developmental origins of neuroblastoma . Nat Genet 53 , 683 – 693 ( 2021 ). OpenUrl CrossRef PubMed 27. ↵ Chalancon , G. et al. Interplay between gene expression noise and regulatory network architecture . Trends Genet 28 , 221 – 232 ( 2012 ). OpenUrl CrossRef PubMed Web of Science 28. ↵ Urban , E. A. & Johnston , R. J . Buffering and Amplifying Transcriptional Noise During Cell Fate Specification . Front Genet 9 , 591 ( 2018 ). OpenUrl CrossRef PubMed 29. ↵ Weinberger , L. S. , Burnett , J. C. , Toettcher , J. E. , Arkin , A. P. & Schaffer , D. V . Stochastic gene expression in a lentiviral positive-feedback loop: HIV-1 Tat fluctuations drive phenotypic diversity . Cell 122 , 169 – 182 ( 2005 ). OpenUrl CrossRef PubMed Web of Science 30. ↵ Maamar , H. , Raj , A. & Dubnau , D . Noise in gene expression determines cell fate in Bacillus subtilis . Science 317 , 526 – 529 ( 2007 ). OpenUrl Abstract / FREE Full Text 31. ↵ Tsuboi , A. et al. Olfactory neurons expressing closely linked and homologous odorant receptor genes tend to project their axons to neighboring glomeruli on the olfactory bulb . J Neurosci 19 , 8409 – 8418 ( 1999 ). OpenUrl Abstract / FREE Full Text 32. ↵ Wernet , M. F. et al. Stochastic spineless expression creates the retinal mosaic for colour vision . Nature 440 , 174 – 180 ( 2006 ). OpenUrl CrossRef PubMed Web of Science 33. ↵ Hastings , J. F. et al. Memory of stochastic single-cell apoptotic signaling promotes chemoresistance in neuroblastoma . Sci. Adv . 9 , eabp8314 ( 2023 ). OpenUrl CrossRef PubMed 34. ↵ Marusyk , A. & Polyak , K. Tumor heterogeneity: causes and consequences . Biochim Biophys Acta 1805 , 105 – 117 ( 2010 ). OpenUrl CrossRef PubMed 35. Marine , J.-C. , Dawson , S.-J. & Dawson , M. A . Non-genetic mechanisms of therapeutic resistance in cancer . Nat Rev Cancer 20 , 743 – 756 ( 2020 ). OpenUrl CrossRef PubMed 36. ↵ Lenz , G. et al. The Origins of Phenotypic Heterogeneity in Cancer . Cancer Res 82 , 3 – 11 ( 2022 ). OpenUrl Abstract / FREE Full Text 37. ↵ Bonine , N. et al. NBAtlas: A harmonized single-cell transcriptomic reference atlas of human neuroblastoma tumors . Cell Reports 43 , 114804 ( 2024 ). OpenUrl PubMed 38. ↵ Chapple , R. H. et al. An integrated single-cell RNA-seq map of human neuroblastoma tumors and preclinical models uncovers divergent mesenchymal-like gene expression programs . Genome Biology 25 , 161 ( 2024 ). OpenUrl PubMed 39. ↵ Kameneva , P. et al. Single-cell transcriptomics of human embryos identifies multiple sympathoblast lineages with potential implications for neuroblastoma origin . Nat Genet 53 , 694 – 706 ( 2021 ). OpenUrl CrossRef PubMed 40. ↵ Patel , A. G. et al. A spatial cell atlas of neuroblastoma reveals developmental, epigenetic and spatial axis of tumor heterogeneity . bioRxiv 2024.01.07.574538 ( 2024 ) doi: 10.1101/2024.01.07.574538 . OpenUrl Abstract / FREE Full Text 41. ↵ Thirant , C. et al. Reversible transitions between noradrenergic and mesenchymal tumor identities define cell plasticity in neuroblastoma . Nat Commun 14 , 2575 ( 2023 ). OpenUrl CrossRef PubMed 42. ↵ Zhou , P. , Wang , S. , Li , T. & Nie , Q . Dissecting transition cells from single-cell transcriptome data through multiscale stochastic dynamics . Nat Commun 12 , 5609 ( 2021 ). OpenUrl CrossRef PubMed 43. ↵ Desai , R. V. et al. A DNA repair pathway can regulate transcriptional noise to promote cell fate transitions . Science 373 , eabc6506 ( 2021 ). OpenUrl Abstract / FREE Full Text 44. ↵ Thanawala , S. U. et al. Regional Modulation of a Stochastically Expressed Factor Determines Photoreceptor Subtypes in the Drosophila Retina . Developmental Cell 25 , 93 – 105 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 45. ↵ Vallejos , C. A. , Marioni , J. C. & Richardson , S . BASiCS: Bayesian Analysis of Single-Cell Sequencing Data . PLoS Comput Biol 11 , e1004333 ( 2015 ). OpenUrl CrossRef PubMed 46. ↵ O’Callaghan , A. , Eling , N. , Marioni , J. C. & Vallejos , C. A . BASiCS workflow: a step-by-step analysis of expression variability using single cell RNA sequencing data . F1000Res 11 , 59 ( 2022 ). OpenUrl 47. ↵ Chodosh , L. A. , Fire , A. , Samuels , M. & Sharp , P. A . 5,6-Dichloro-1-β-D-ribofuranosylbenzimidazole inhibits transcription elongation by RNA polymerase II in vitro . Journal of Biological Chemistry 264 , 2250 – 2257 ( 1989 ). OpenUrl Abstract / FREE Full Text 48. ↵ Bartman , C. R. et al. Transcriptional Burst Initiation and Polymerase Pause Release Are Key Control Points of Transcriptional Regulation . Mol Cell 73 , 519 – 532 .e4 ( 2019 ). OpenUrl CrossRef PubMed 49. ↵ Xu , J. , Lamouille , S. & Derynck , R . TGF-beta-induced epithelial to mesenchymal transition . Cell Res 19 , 156 – 172 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 50. ↵ Katsuno , Y. & Derynck , R . Epithelial plasticity, epithelial-mesenchymal transition, and the TGF-β family . Developmental Cell 56 , 726 – 746 ( 2021 ). OpenUrl CrossRef PubMed 51. ↵ Shaffer , S. M. et al. Memory Sequencing Reveals Heritable Single-Cell Gene Expression Programs Associated with Distinct Cellular Behaviors . Cell 182 , 947 – 959 .e17 ( 2020 ). OpenUrl CrossRef PubMed 52. ↵ Yoo , J. et al. Drawing a line between histone demethylase KDM5A and KDM5B: their roles in development and tumorigenesis . Exp Mol Med 54 , 2107 – 2117 ( 2022 ). OpenUrl CrossRef PubMed 53. ↵ Hinohara , K. et al. KDM5 Histone Demethylase Activity Links Cellular Transcriptomic Heterogeneity to Therapeutic Resistance . Cancer Cell 34 , 939 – 953 .e9 ( 2018 ). OpenUrl CrossRef PubMed 54. ↵ Chen , E. et al. The Role of Phenotypic Plasticity in Adaptation to Treatment and Prospective Plasticity Drivers in Hepatoblastoma . 2025.04.16.648783 Preprint at doi: 10.1101/2025.04.16.648783 ( 2025 ). OpenUrl Abstract / FREE Full Text 55. ↵ Hirsch , T. Z. et al. Integrated Genomic Analysis Identifies Driver Genes and Cisplatin-Resistant Progenitor Phenotype in Pediatric Liver Cancer . Cancer Discov 11 , 2524 – 2543 ( 2021 ). OpenUrl Abstract / FREE Full Text 56. ↵ Roehrig , A. et al. Single-cell multiomics reveals the interplay of clonal evolution and cellular plasticity in hepatoblastoma . Nat Commun 15 , 3031 ( 2024 ). OpenUrl PubMed 57. ↵ Ponzoni , M. et al. Recent advances in the developmental origin of neuroblastoma: an overview . J Exp Clin Cancer Res 41 , 92 ( 2022 ). OpenUrl CrossRef PubMed 58. ↵ Bell , C. C. & Gilan , O . Principles and mechanisms of non-genetic resistance in cancer . Br J Cancer 122 , 465 – 472 ( 2020 ). OpenUrl CrossRef PubMed 59. ↵ Rosano , D. et al. Long-term Multimodal Recording Reveals Epigenetic Adaptation Routes in Dormant Breast Cancer Cells . Cancer Discov 14 , 866 – 889 ( 2024 ). OpenUrl CrossRef PubMed 60. ↵ Chang , H. H. , Hemberg , M. , Barahona , M. , Ingber , D. E. & Huang , S . Transcriptome-wide noise controls lineage choice in mammalian progenitor cells . Nature 453 , 544 – 547 ( 2008 ). OpenUrl CrossRef PubMed Web of Science 61. ↵ Gupta , P. B. et al. Stochastic state transitions give rise to phenotypic equilibrium in populations of cancer cells . Cell 146 , 633 – 644 ( 2011 ). OpenUrl CrossRef PubMed Web of Science 62. ↵ García-Blay , Ó. et al. Multimodal screen identifies noise-regulatory proteins . Dev Cell 60 , 133 – 151 .e12 ( 2025 ). OpenUrl PubMed 63. ↵ Shendy , N. A. M. , Zimmerman , M. W. , Abraham , B. J. & Durbin , A. D . Intrinsic transcriptional heterogeneity in neuroblastoma guides mechanistic and therapeutic insights . Cell Reports Medicine 3 , 100632 ( 2022 ). OpenUrl PubMed 64. ↵ Jose , A. et al. Histone Demethylase KDM5B as a Therapeutic Target for Cancer Therapy . Cancers 12 , 2121 ( 2020 ). OpenUrl PubMed 65. ↵ Chirnomas , D. , Hornberger , K. R. & Crews , C. M . Protein degraders enter the clinic - a new approach to cancer therapy . Nat Rev Clin Oncol 20 , 265 – 278 ( 2023 ). OpenUrl CrossRef PubMed 66. Iida , T. et al. Identification of Proteolysis Targeting Chimeras (PROTACs) for Lysine Demethylase 5 and Their Neurite Outgrowth-Promoting Activity . Chem Pharm Bull (Tokyo) 72 , 638 – 647 ( 2024 ). OpenUrl PubMed 67. ↵ Guan , T. et al. Discovery of an efficacious KDM5B PROTAC degrader GT-653 up-regulating IFN response genes in prostate cancer . European Journal of Medicinal Chemistry 272 , 116494 ( 2024 ). OpenUrl PubMed 68. ↵ Suarez-Arnedo , A. et al. An image J plugin for the high throughput image analysis of in vitro scratch wound healing assays . PLoS ONE 15 , e0232565 ( 2020 ). OpenUrl CrossRef PubMed 69. ↵ Griffiths , J. A. , Richard , A. C. , Bach , K. , Lun , A. T. L. & Marioni , J. C . Detection and removal of barcode swapping in single-cell RNA-seq data . Nat Commun 9 , 2667 ( 2018 ). OpenUrl CrossRef PubMed 70. ↵ Lun , A. T. L. et al. EmptyDrops: distinguishing cells from empty droplets in droplet-based single-cell RNA sequencing data . Genome Biol 20 , 63 ( 2019 ). OpenUrl CrossRef PubMed 71. ↵ Satija , R. , Farrell , J. A. , Gennert , D. , Schier , A. F. & Regev , A . Spatial reconstruction of single-cell gene expression data . Nat Biotechnol 33 , 495 – 502 ( 2015 ). OpenUrl CrossRef PubMed 72. ↵ Wu , T. et al. clusterProfiler 4.0: A universal enrichment tool for interpreting omics data . The Innovation 2 , 100141 ( 2021 ). OpenUrl PubMed 73. ↵ Yu , G. , Wang , L.-G. , Han , Y. & He , Q.-Y . clusterProfiler: an R Package for Comparing Biological Themes Among Gene Clusters . OMICS: A Journal of Integrative Biology 16 , 284 – 287 ( 2012 ). OpenUrl CrossRef PubMed 74. ↵ Tan , T. Z. et al. Epithelial-mesenchymal transition spectrum quantification and its efficacy in deciphering survival and drug responses of cancer patients . EMBO Mol Med 6 , 1279 – 1293 ( 2014 ). OpenUrl Abstract / FREE Full Text 75. ↵ Li , H. & Durbin , R . Fast and accurate short read alignment with Burrows–Wheeler transform . Bioinformatics 25 , 1754 – 1760 ( 2009 ). OpenUrl CrossRef PubMed Web of Science 76. ↵ Danecek , P. et al. Twelve years of SAMtools and BCFtools . GigaScience 10 , giab008 ( 2021 ). OpenUrl CrossRef PubMed 77. ↵ Shannon , C. E . A Mathematical Theory of Communication . Bell System Technical Journal 27 , 379 – 423 ( 1948 ). OpenUrl CrossRef 78. ↵ Dobin , A. et al. STAR: ultrafast universal RNA-seq aligner . Bioinformatics 29 , 15 – 21 ( 2013 ). OpenUrl CrossRef PubMed Web of Science 79. ↵ Berlanga , P. et al. The European MAPPYACTS Trial: Precision Medicine Program in Pediatric and Adolescent Patients with Recurrent Malignancies . Cancer Discov 12 , 1266 – 1281 ( 2022 ). OpenUrl CrossRef PubMed 80. ↵ Marques Da Costa , M. E. , et al. A biobank of pediatric patient-derived-xenograft models in cancer precision medicine trial MAPPYACTS for relapsed and refractory tumors . Commun Biol 6 , 949 ( 2023 ). OpenUrl PubMed 81. ↵ Stewart , E. et al. Orthotopic patient-derived xenografts of paediatric solid tumours . Nature 549 , 96 – 100 ( 2017 ). OpenUrl CrossRef PubMed 82. ↵ Boeva , V. et al. Heterogeneity of neuroblastoma cell identity defined by transcriptional circuitries . Nat Genet 49 , 1408 – 1413 ( 2017 ). OpenUrl CrossRef PubMed 83. ↵ Wolf , F. A. , Angerer , P. & Theis , F. J . SCANPY: large-scale single-cell gene expression data analysis . Genome Biol 19 , 15 ( 2018 ). OpenUrl CrossRef PubMed 84. ↵ Zhou , P. , Wang , S. , Li , T. & Nie , Q . Dissecting transition cells from single-cell transcriptome data through multiscale stochastic dynamics . Nat Commun 12 , 5609 ( 2021 ). OpenUrl CrossRef PubMed 85. ↵ Kildisiute , G. et al. Tumor to normal single-cell mRNA comparisons reveal a pan-neuroblastoma cancer cell . Sci Adv 7 , eabd3311 ( 2021 ). OpenUrl FREE Full Text 86. ↵ Hänzelmann , S. , Castelo , R. & Guinney , J . GSVA: gene set variation analysis for microarray and RNA-Seq data . BMC Bioinformatics 14 , 7 ( 2013 ). OpenUrl CrossRef PubMed 87. ↵ Bruna , A. et al. A Biobank of Breast Cancer Explants with Preserved Intra-tumor Heterogeneity to Screen Anticancer Compounds . Cell 167 , 260 – 274 .e22 ( 2016 ). OpenUrl CrossRef PubMed 88. ↵ Zheng , Y. , Ahmad , K. & Henikoff , S . CUT&Tag Data Processing and Analysis Tutorial v1 . Preprint at doi: 10.17504/protocols.io.bjk2kkye ( 2020 ). OpenUrl CrossRef View the discussion thread. Back to top Previous Next Posted October 14, 2025. Download PDF Supplementary Material Email Thank you for your interest in spreading the word about bioRxiv. NOTE: Your email address is requested solely to identify you as the sender of this article. Your Email * Your Name * Send To * Enter multiple addresses on separate lines or separate them with commas. You are going to email the following KDM5-driven transcriptional noise fuels plasticity-led awakening and relapse in paediatric cancer Message Subject (Your Name) has forwarded a page to you from bioRxiv Message Body (Your Name) thought you would like to see this page from the bioRxiv website. Your Personal Message CAPTCHA This question is for testing whether or not you are a human visitor and to prevent automated spam submissions. Share KDM5-driven transcriptional noise fuels plasticity-led awakening and relapse in paediatric cancer Alejandro Allo Anido , Cecilia Roux , Emilia Chen , Siân Hamer , Abigail Shea , Ayeh Sadat Sadr , Christie English , Charlotte Butterworth , Harvey Che , Angella Bellini , Birgit Geoerger , Gudrun Schleiermacher , Louis Chesler , Michael David Morgan , Alejandra Bruna bioRxiv 2025.10.13.682004; doi: https://doi.org/10.1101/2025.10.13.682004 Share This Article: Copy Citation Tools KDM5-driven transcriptional noise fuels plasticity-led awakening and relapse in paediatric cancer Alejandro Allo Anido , Cecilia Roux , Emilia Chen , Siân Hamer , Abigail Shea , Ayeh Sadat Sadr , Christie English , Charlotte Butterworth , Harvey Che , Angella Bellini , Birgit Geoerger , Gudrun Schleiermacher , Louis Chesler , Michael David Morgan , Alejandra Bruna bioRxiv 2025.10.13.682004; doi: https://doi.org/10.1101/2025.10.13.682004 Citation Manager Formats BibTeX Bookends EasyBib EndNote (tagged) EndNote 8 (xml) Medlars Mendeley Papers RefWorks Tagged Ref Manager RIS Zotero Tweet Widget Facebook Like Google Plus One Subject Area Cancer Biology Subject Areas All Articles Animal Behavior and Cognition (7629) Biochemistry (17660) Bioengineering (13881) Bioinformatics (41912) Biophysics (21436) Cancer Biology (18578) Cell Biology (25482) Clinical Trials (138) Developmental Biology (13372) Ecology (19889) Epidemiology (2067) Evolutionary Biology (24302) Genetics (15599) Genomics (22483) Immunology (17728) Microbiology (40365) Molecular Biology (17163) Neuroscience (88540) Paleontology (666) Pathology (2830) Pharmacology and Toxicology (4821) Physiology (7637) Plant Biology (15130) Scientific Communication and Education (2045) Synthetic Biology (4290) Systems Biology (9818) Zoology (2269)