HMGB1 as a rheostat of chromatin topology and RNA homeostasis on the path to senescence

Spatial organization and gene expression of mammalian chromosomes are maintained and regulated in conjunction with cell cycle progression. This link is perturbed once cells enter senescence. The highly abundant HMGB1 protein, known to associate with bent and looped DNA, is actively depleted from senescent cell nuclei to act as an extracellular proinflammatory stimulus. Despite its physiological importance, we still lack understanding of the positioning and functional roles of HMGB1 on chromatin in vivo. To address this, we mapped HMGB1 binding genome-wide in different primary cells using a tailored protocol. We then integrated ChIP-seq and Hi-C data with a knot theory approach to uncover HMGB1 demarcation at the boundaries of particular topologically-associating domains (TADs). These TADs harbor genes involved in the key proinflammatory leg of the senescent transcriptional program. Moreover, we used sCLIP and siRNA-mediated knockdown to show that HMGB1 is a bona fide RNA-binding protein also affecting splicing choices. Together, our findings highlight a broader than hitherto assumed role for HMGB1 in chromatin homeostasis connected to cell cycle potency, and allow us to postulate a “rheostat” model for HMGB function with implications in cancer biology.


Introduction
The high-mobility group box 1 (HMGB1) protein, a member of the highly conserved non-histone DNA binding HMG protein family, was named after its characteristically rapid electrophoretic mobility (Štros, 2010). HMGB1, together with histone H1, are the next most abundant proteins after core histones, with one HMGB1 molecule for every ~10 nucleosomes . Despite its documented abundance in cell nuclei, HMGB1 (also known as an "alarmin") has been predominantly studied in vivo as an extracellular signaling molecule (Lohani and Rajeswari, 2016;Bianchi et al, 2017). HMGB1 is actively secreted by activated monocytes and macrophages to signal tissue damage, and passively released by necrotic and damaged cells. Once received by cells, HMGB1 molecules can be recognized with high affinity by RAGE receptors to potently signal inflammation (Scaffidi et al, 2002;Bonaldi et al, 2003;Orlova et al, 2007). In cells entering senescence, HMGB1 relocalizes from the nucleus to the cytoplasm and is then secreted to stimulate NF-κB activity via Toll-like receptor signaling. This relocalization and secretion controls and contributes to the senescence-associated secretory phenotype (SASP) of mammalian cells, representing a major paracrine contributor both in vitro and in vivo (Salminen et al, 2012;Acosta et al, 2013;Davalos et al, 2013).
Within proliferating cell nuclei, HMGB1 is thought to bind DNA in a nonspecific manner via its two HMGB-box domains, and to bend and contort the double helix -a mode of action that facilitates recruitment of various DNA-binding factors, like p53 (Štros et al, 2004;Štros, 2010;Rowell et al, 2012). HMGB1 has been studied in detail for its contribution to DNA repair (Ito et al, 2015;Mukherjee et al, 2016), V(D)J recombination (Little et al, 2013;Zagelbaum et al, 2016) or chromatin assembly (Bonaldi et al, 2002), but not so much for its transcriptional role (Calogero et al, 2009). Notably, cells lacking HMGB1 also contain reduced numbers of nucleosomes, rendering chromatin more susceptible to DNA damage, spurious transcription, and inflammatory activation (Giavara et al, 2005;El Gazzar et al, 2009;Celona et al, 2011;De Toma et al, 2014).
HMGB1 associates with its cognate DNA sites with characteristically high "on"/"off" rates, and its acidic tail is important for stabilizing binding (Pallier et al, 2003;Ueda et al, 2004;Štros, 2010;Blair et al, 2016). However, it is now understood that HMG-box DNA-binding domains are particularly sensitive to commonly-used fixatives, like formaldehyde; thus, capturing them on chromatin can be challenging (Pallier et al, 2003;Teves et al, 2016). As a result, there exist no genome-wide datasets describing HMGB binding repertoires in primary mammalian cells (see http://chip-atlas.org), and our appreciation of their on-chromatin roles remains vague. To address this, we employed here a tailored dual-crosslinking approach previously used to map HMGB2 binding genome-wide (Zirkel et al, 2018). We can now show that HMGB1 binding in proliferating primary endothelial and lung fibroblast cells is far from nonspecific, while also disparate to that by HMGB2. Once integrated with whole-genome chromosome conformation capture (Hi-C) data from proliferating and senescent cells, HMGB1-bound positions demarcate the boundaries of a considerable and specific fraction of topologically-associating domains (TADs;Dixon et al, 2012;Nora et al, 2017). This topological contribution is eliminated upon senescence entry, and knockdown experiments show that HMGB1 directly controls the expression of genes central to the senescent transcriptional program. Critically, as HMGB1 was proposed to also constitute a bona fide RNA-binding protein (Castello et al, 2016), we used sCLIP (Kargapolova et al, 2017) to show it does also influence RNA splicing. In summary, we use senescence as a model to comprehensively characterize the multifaceted in nucleo roles of HMGB1 and show how these converge to underpin transition into replicative arrest by linking chromatin homeostasis with cell cycle progression -a link with far-reaching implications for the attenuation of proliferation of cancer cells overexpressing HMGB1.

HMGB1 nuclear loss marks senescence entry coinciding with chromatin changes
To investigate the in nucleo roles of HMGB1 across cellular contexts we used primary human umbilical vein endothelial cells (HUVECs) and fetal lung fibroblasts (IMR90) that are of distinct developmental origins and have disparate gene expression programs. We defined an early-passage proliferative state and a late-passage senescent state by combining β-galactosidase staining, cell cycle staging by FACS, and MTT proliferation assays ( Figure 1A). Next, we used RNAseq data from proliferating and senescent HUVEC and IMR90 (from two different biological donors/isolates) to look into changing mRNA levels of chromatin-associated factors. Changes that were convergent between the two cell types involved genes encoding lamin B1, various histones, centrosomal proteins, cohesin and condensin complexes, as well as HMG family proteins, most of which were consistently suppressed upon senescence entry, also at the protein level (Figure 1B,C;Davalos et al, 2013;Shah et al, 2013;Rai et al, 2014;Zirkel et al, 2018). Here, we chose to focus on HMGB1 due to its high nuclear abundance ; Figure 1D) and key role in SASP induction (Davalos et al, 2013), but mostly due to its still elusive onchromatin functions, especially as regards spatial chromatin folding.
HMGB1 immunodetection in early-and late-passage cells documented a >50% decrease in its total nuclear levels in senescence-entry cell populations of HUVECs or IMR90 ( Figure 1E). However, HMGB1 nuclear depletion was most dramatic in the enlarged nuclei of senescent cells of either cell type, while smaller nuclei remained largely unaffected. FACS-sorting IMR90 based on light scattering allowed enrichment for cell populations with enlarged nuclei (i.e., ~70% of cells had larger than average nuclei, with >35% being >1.5-fold larger than the average nucleus of a proliferating cell; Figure  S1A). This showed that, within such heterogeneous populations, enlarged nuclei almost invariably represent senescent cells lacking HMGB1 and concomitantly exhibiting a drop in H3K27me3 levels that mark facultative heterochromatin -effects which would otherwise be masked ( Figure 1C and S1B,C). Last, we showed that it is those larger cells that secrete HMGB1, but not HMGB2, into their growth medium to contribute to the SASP ( Figure S1D).

Regulated HMGB1 nascent transcription in single cells
Since senescence entry is idiosyncratic to each individual cell and occurs asynchronously in a given cell population, it is important to obtain a single-cell understanding of the transcriptional changes linked to the nuclear loss of HMGB1. To this end, we developed a new protocol for single-cell sequencing of nascent RNA. Typically, scRNA-seq relies on capturing and barcoding cellular mRNAs via reversetranscription of their long poly(A) tails (See et al, 2018). To capture nascent RNA, we decided to add poly(A) tails to any RNA associated with an active site of transcription. We employed the previouslyintroduced "factory" RNA-seq protocol (Caudron-Herger et al, 2015) as a basis for isolating intact nuclei rich in nascent transcripts from both proliferating and senescent HUVEC. We next removed most chromatin not attached to active transcription sites, and polyadenylated nascent transcripts in situ, before a standard library preparation via a 10X Genomics platform (overview in Figure S2A). Using this new approach, and despite low rates of non-duplet capturing (see Methods for more details), we obtained >550 single nuclei with an average of 1,650 transcripts in each nucleus, >55% of which represented introns or intron-exon junctions. Unsupervised t-SNE clustering of single-nucleus transcriptomes grouped all senescent cells in a single cluster, while subdividing the larger proliferating HUVEC population into three clusters ( Figure S2B). Reassuringly, genes differentially-expressed amongst proliferating and senescent clusters associated with GO terms relevant to senescence entry and SASP induction ( Figure S2C).
Mapping nascent RNA levels of individual gene loci onto those clusters, showed that HMGB1 and HMGB2 are actively repressed not only in senescent cells, but also in numerous proliferating cells; the senescence-induced HMGA1/A2 loci strongly produce nascent RNA in senescent cells, but are markedly upregulated also in proliferating cells showing lower HMGB1/B2 transcription (e.g., in cells of the blue cluster; Figure S2B). As expected, SASP-related genes like IL4R and MMP14 show strong transcription in cells repressing HMGB1/B2, while the p21-encoding locus, CDKN1A, has most nascent RNA detected in the senescent cluster ( Figure S2D). Notably, correlating HMGB1 and HMGB2 levels in senescent cells confirms that switching off the latter locus is an early event on the path to senescence (as recently proposed; Zirkel et al, 2018), and this trend is already obvious in proliferating cells ( Figure  S2E). Similarly, correlating HMGB1 levels with those of CDKN1A, HMGA1 or IL4R, all of which are activated upon senescence entry, we observed that their hyperactivation mostly occurs in those cells where HMGB1 nascent RNA production is diminished ( Figure S2E).

HMGB1 binds active chromatin loci in a cell type-specific manner
Capturing HMGB proteins bound to chromatin has proven challenging, because their tandem HMGbox DNA-binding domains are not compatible with standard formaldehyde fixation (Pallier et al, 2003;Belmont et al, 2015;Teves et al, 2016). Here, we employed a tailored fixation strategy in ChIP assays to efficiently capture HMGB1 bound to its cognate sites genome-wide in both proliferating HUVECs and IMR90 (Figures 2A and S3A; see Methods for details). HMGB1 ChIP-seq peaks were predominantly found at the promoters and gene bodies of active genes (>75% and >65% of 810 and 593 peaks in HUVECs and IMR90, respectively; Figure 2B), overlapping regions marked by H3K27ac and oftentimes H2A.Z ( Figure 2C). HMGB1 peaks were often found clustering along chromosomes of proliferating cells ( Figure S3B), but as HMGB1 mostly binds active loci, and HUVECs have a gene expression program disparate to that of IMR90, the two repertoires only share 40 peaks ( Figure S3C). It is noteworthy that none of these peaks mark SASP-related genes, despite HMGB1's role in proinflammatory stimulation.
De novo motif discovery in accessible "footprints" within DNase I-hypersensitive sites (derived from ENCODE data) overlapping HMGB1 peaks revealed rather specific CG-rich motifs, which do not however converge into a single consensus sequence ( Figure S3D). We also surveyed these same accessible footprints for known transcription factor (TF) motifs to infer possible co-binding complexes or complexes that might replace HMGB1 on senescent chromatin. We found that HMGB1 binding sites are significantly enriched for E2F-family motifs, which are notably down-regulated upon senescence entry. Also, motifs for senescence-activated corepressors (e.g., REST and HEY2) and for the architectural ZBTB7B protein that is important for inflammatory gene activation (Nikopoulou et al, 2018) are also enriched therein ( Figure S3E). Last, we combined RNA-seq with HMGB1 ChIP-seq data to catalogue differentially-expressed genes in either cell type that are directly regulated upon HMGB1 loss in senescence. In both HUVECs and IMR90 we find >2-fold more up-rather than downregulated genes bound by HMGB1 ( Figure S3F). These senescence-induced genes are involved in ECM organization, regulation of apoptosis, as well as with NOTCH-/TGFβ-signalling that is now understood to represent a second "wave" of paracrine signalling in senescent cells (Hoare et al, 2017). On the other hand, downregulated genes in both cell types converge to cell growth and cell cycle regulation ( Figure S3F). Together, this data demonstrates how HMGB1 nuclear eviction regulates a critical leg of the program defining senescence entry.

HMGB1 demarcates a particular subset of senescence-regulated TAD boundaries
The boundaries of TADs represent genomic sites of strong local insulation of spatial interactions, and are often marked by the presence of CTCF and/or active gene promoters (Dixon et al, 2012;Nora et al, 2017). We recently showed that a considerable number of TAD boundaries in proliferating human cells are marked by HMGB2, and this demarcation is lost upon senescence entry leading to the formation of senescence-induced CTCF foci (Zirkel et al, 2018). Based on this, we reasoned that HMGB1 might also function similarly Hi-C data from proliferating and senescent HUVEC to investigate this. First, we found that >20% of the 810 HMGB1 peaks reside at TAD boundaries (called at 25-kbp resolution; Figure 2D,E). Remarkably, and much unlike HMGB2, we identified TAD boundaries marked by HMGB1, but lacking CTCF (and the converse; Figure 2E). As expected, TADs of senescent nuclei lose this demarcation (Figure 2F), and HMGB1-marked boundaries exhibit an obvious loss of insulation and rearrangement of spatial interactions ( Figure 2G). Next, we grouped TADs according to whether or not their boundaries change upon senescence entry. We found that the ~1,000 TADs remaining invariant in senescence, boundaries are more enriched for CTCF rather than HMGB1 binding. CTCF and HMGB1 demarcation is comparable in TADs shifting one boundary or collapsing into a larger TAD (~5,500 TADs in total; Figure 2H). Thus, we postulate that the loss of HMGB1 in senescence contributes to the remodeling of TADs. Last, we asked whether HMGB1 engages in long-range loop formation; the 810 HMGB1-bound genomic sites give rise to >600 interchromosomal loops, which appear to collapse upon senescence entry ( Figure 2I).

Spatial TAD co-association reveals functional chromosome compartmentalization
It is now understood that human chromosomes undergo large-scale changes upon replicative arrest (Zirkel et al, 2018), which are accentuated in "deep" senescence (Criscione et al, 2016). Looking at TADs as building blocks of chromosomes (Sexton and Cavalli, 2015), one can ask how these might spatialuly associate to give rise to "meta-TADs" (Fraser et al, 2015) and ultimately to diverse functional topologies. To address this question, we employed a bias-free and unsupervised approach based on "topologically-intrinsic lexicographic ordering" (TiLO; Johnson, 2014), whereby TADs are treated as nodes in a clustered spatial network tested for robustness via iterative network slicing (see Figure 3A and Methods for details). By applying this to TADs derived from proliferating and senescent HUVEC Hi-C data, we found that there is a marked increase in clusters that include multiple consecutive TADs in senescence, indicative of TAD merging (Figure 3B,C). This is consistent with the observation of general spatial chromatin compaction in senescence (Criscione et al, 2016;Zirkel et al, 2018) and with the fact that ~75% of TADs merge in senescence ( Figure 2G). However, we also identified individual ("singular") TADs typically positioned between multi-TAD clusters (Figure 3C,arrows). Strikingly, the boundaries of singular TADs were strongly marked by HMGB1, compared to the extremities of clusters comprised of >3 consecutive TADs ( Figure 3D). To assess the functional impact of this, we investigated expression level changes of genes embedded in different cluster types. Both singular TADs and TADs with >3 consecutive TADs harbor genes significantly more upregulated that those in all other TADs, but singular TADs also show significantly more modest gene downregulation ( Figure 3E). Notably, these two subsets harbor most of the genes differentially-regulated upon senescence entry. What further discriminates these clusters functionally, is the fact that multi-TAD clusters uniquely harbor genes involved in the regulation of chromatin organization and conformation, while those unique to singular TADs associate with SASP production and its downstream effects ( Figure 3F). Thus, TiLO identifies functional entities that involve spatial co-association of multiple TADs in cis, explaining the functional specialization of different genomic segments upon HMGB1 loss.

HMGB1 directly binds and regulates senescence-relevant transcripts
On top of its known ability to bind chromatin, HMGB1 was recently also classified as a bona fide RNAbinding protein (Castello et al, 2016;Trendel et al, 2018). Thus, its loss in senescence could also affect some aspect of RNA metabolism. Our suspicion was reinforced by cataloguing the proteins that coimmunoprecipitate with HMGB1 from proliferating IMR90 using mass spectrometry. This returned a diverse range of RNA-binding proteins and splicing regulators, in addition to the expected proteins involved in the regulation of chromatin conformation ( Figure 4A).
To study HMGB1 as a direct RNA binder and regulator, we performed sCLIP (Kargapolova et al, 2017) in proliferating IMR90 in two well correlated biological replicates (Figure S4A-C). These two datasets provided a compendium of 1,773 binding peaks on 866 different transcripts (e.g., ASH1L and CCNL2; Figure 4B and Table S1), which reassuringly display <25% overlap to HMGB1-bound genes in ChIP-seq. On RNA, HMGB1 was mostly found bound to exons and 5'/3' UTRs, but also to a substantial number of non-coding RNAs ( Figure 4C). These HMGB1 RNA-binding sites present a hexameric 5'-NMWGRA-3' (M=A/C, W=A/T, R=A/G) motif, irrespective of the predicted folding of the underlying RNA or the exclusion of repeat sequences from the genome build used for motif analysis (Figures 4D  and S4D). Much like the trend observed in ChIP-seq data, HMGB1 was found bound to ~3-fold more transcripts that are up-rather than downregulated upon senescence ( Figure S4E, left). Upregulated transcripts associated with senescence-relevant GO terms like "ECM organization", "wound healing" and "negative regulation of cell proliferation", while downregulated ones mostly converged on processes like RNA splicing, RNA-/miRNA-mediated gene silencing, or histone remodeling and deacetylation (Figure S4E,right). This suggests a feedback loop whereby transcripts encoding RNA regulators are affected by the senescence-induced loss of HMGB1, but also pointed to a more direct impact on splicing regulation.
We next examined how splicing patterns are altered upon senescence entry by IMR90 using Whippet (Sterne-Weiler et al, 2018). We documented significant changes involving ~4,000 transcripts, the majority of which concerned the use of alternative transcription start and polyadenylation sites (>80% of cases; Figure 4E). Interestingly, looking at the splicing changes involving 342 HMGB1-bound mRNAs, this trend remained invariant ( Figure 4E). Finally, we asked whether transcripts undergoing different types of splicing changes are involved in different functional processes. GO term analysis revealed that transcripts with alternative TSS usage encoded regulators of ECM organization and growth control, whereas those with alternative poly(A)-site usage encoded mostly splicing and RNA processing factors; HMGB1-bound transcripts, however, associated with functions from the whole of this spectrum ( Figure 4F). Thus, the nuclear loss of HMGB1 also directly affects processing of the cell's transcriptome.

HMGB1 depletion underlies induction of the senescence transcriptional program
It was previously demonstrated that transduction of human fibroblasts with shRNAs against HMGB1 sufficed to induce senescence (Davalos et al, 2013). To avoid using lentiviral vectors, we first treated HUVECs with a pool of self-delivering siRNAs targeting HMGB1. This led to a ~2-fold reduction of HMGB1 protein and RNA levels within 72 h (Figure S5A), was accompanied by a doubling of β-gal and p21-positive cells in the knockdown population, but only small changes in nuclear size and BrdU incorporation were recorded (Figure S5B-E). We subsequently turned to IMR90, where standard siRNA transfections allowed for a >10-fold decrease in HMGB1 levels, while also suppressing or inducing expression of known senescence-regulated genes (but only marginally affecting HMGB2; Figure 5A). "Deep" sequencing and analysis of mRNAs from HMGB1-knockdown and control IMR90 revealed ~900 up-and ~950 downregulated genes ( Figure 5B). GO term and gene set enrichment analyses showed that genes that were upregulated associated with the SASP and proinflammatory signalling, while those that were downregulated with changes in chromatin organization, transcriptional silencing, and the p53 pathway (Figure 5C,D).
For a more precise understanding for the role of HMGB1 in these processes, we focused on genes differentially-regulated upon HMGB1-knockdown that also carried at least one HMGB1 ChIPseq peak. This returned 44 up-and 56 downregulated genes constituting direct HMGB1 targets linked to NF-κB activation, and to chromatin organization or non-inflammatory signalling, respectively ( Figure S5F). Nonetheless, this left >850 genes in either category that could not be directly linked to an HMGB1 chromatin-binding event. To bridge this gap, we repeated the above analysis using sCLIP binding events. We found 56 and 97 HMGB1-bound mRNAs up-and downregulated upon HMGB1knockdown, respectively. Consistent to our previous observations, those upregulated could be linked to processes like ECM organization and proinflammatory activation, while downregulated ones to non-inflammatory signalling and the organization of chromatin ( Figure S5G). Notably, splicing changes deemed significant upon HMGB1-KD displayed an trend in favour of alternative TSS and poly(A)-site usage as, essentially identical to that observed in senescence, with >60% of senescence splicing events also being recorded in knockdown cells ( Figure 5E). These mRNAs encode factors linked to senescence-regulated processes like cell cycle and cell growth regulation or the p53 pathway ( Figure  5F). Taken together, our data are in support of a model whereby HMGB1 acts both at chromatin loci and on RNA transcripts to regulate cellular functions, and this interplay is disrupted in senescent cells by its nuclear depletion.

An HMGB1-dependency for lung cancer cell proliferation
The nuclear depletion of HMGB1 and its ensuing transcriptional changes are senescence hallmarks and represent an inherent antitumorigenic mechanism . However, HMGB1 is very often overexpressed in patient tumor specimens (Figure 6A), and it has been suggested that its increased titers might result in either the enhancement of cell survival or the regulation of cell death via diverse pathways that include inflammation, cell proliferation, and autophagy (Nagatani et al, 2001;Kang et al, 2013). We focused on non-small cell lung adenocarcinoma patient-derived lines, because high HMGB1 expression therein is one predictor of poor patient survival ( Figure 6B). We used three different lines, H1975, A549, and H1650, and found that cell proliferation and survival was significantly impaired via HMGB1-knockdown in 72 h in all three ( Figure 6C). Interestingly, the extent of proliferative impediment in each line was almost proportional to its proliferation rate, with the faster H1975 cells arresting completely and the slower H1650 still exhibiting some population increase upon knockdown (Figure 6C,D). Nonetheless, HMGB1-knockdown in these lines showed convergent changes in the expression of senescence markers like HMGB2, LMNB1, EZH2 or SMC2 (Figure 6E), which are actually overexpressed in correlation to HMGB1 levels in lung adenocarcinoma tumours (data from the TCGA cohort; Figure 6F). Since knocking down HMGB1 in this context also decreased HMGB2 levels, and we recently showed that reduced HMGB2 levels result in the formation of senescence-induced CTCF clusters (SICCs), we examined whether SICCs formed differentially in these lines. Indeed, we saw that the most affected H1975 cells exhibited larger and more prominent SICCs upon HMGB1-knockdown, compared to the least affected H1650 cells that displayed essentially no SICCs; the intermediately-affected A549 presented SICCs in both control and knockdown cells ( Figure  S6A-C). Thus, differences in proliferative capacity and SICC emergence correlate with the antiproliferative extent of HMGB1 depletion in these cancer lines, highlighting how targeting HMGB1 might need to be considered in therapeutic interventions for cancer treatment.

Discussion
Unlike its well documented proinflammatory role, the intracellular positioning and gene expression control exerted by HMGB1 on mammalian chromosomes is poorly understood. Here, we were able to assign a multifaceted role to HMGB1 -as an on-chromatin regulator of active gene loci, and as a bona fide RNA-binding protein recognizing a distinct subset of transcripts. As a result, we can deduce the following. First, we that the loss of HMGB1 from the nuclei of senescent cells mostly triggers upregulation of its previously-bound target loci and mRNAs, suggesting that HMGB1 tends to act a "buffering" factor thereon. Second, that ~1/5 of HMGB1-bound positions mark the boundary of a TAD, and many of these TADs specifically harbor SASP-related genes induced upon senescence entry. Last, that the loss of HMGB1 correlates essentially only with alternative usage of TSS and polyadenylation sites in its bound transcripts. These observations come to substantiate a previous hypothesis that low nuclear titers of HMGB1 are necessary for the fully-fledged development of the SASP (Davalos et al, 2013). This is because there is a need for alleviating the regulatory effects that HMGB1 exerts on active promoters, on mRNAs being processed, as well as on TAD boundaries. This is a rather unique example of a regulatory circuit, where the programmed deregulation in one cellular compartment (the nucleus) is in direct and quantitative crosstalk with the signaling deployed in another (in extracellular space to initiate paracrine activation). Thus, the senescent regulatory program has a strong, if not hierarchical, dependency on the nuclear events preceding SASP deployment.
Recently, we characterized the function of the sister protein to HMGB1, HMGB2, for the entry into replicative senescence (Zirkel et al, 2018). The loss of HMGB2 is an event preceding the loss of HMGB1, and leads to the formation of large senescence-induced CTCF clusters (SICCs). This has an apparent effect on the spatial architecture of chromosomes, and concomitantly on gene expression. Intriguingly, direct HMGB2 target loci are also typically upregulated when relieved of HMGB2 binding; however, this is the only similarity between the functions of HMGB1 and HMGB2. The loss of HMGB1 does not trigger SICC formation, the same way that the loss of HMGB2 does not trigger immediate senescence entry. Also, HMGB1 and HMGB2 bind non-overlapping targets and also demarcate TADs in two distinct modes -HMGB2 marks the extremities of TADs that shift one boundary in senescence, while HMGB1 is mostly found at the boundaries of TADs that collapse together, in line with the overall compaction observed in senescent chromosomes (Criscione et al, 2016;Zirkel et al, 2018). Finally, it is important to note that the HMGB1 knockdown does not reduce HMGB2 levels in primary human cells, nor does the converse hold true, meaning that the two pathways these related factors control do not really cross one another, but are rather deployed in parallel.
Interestingly, the knockdown of HMGB1 in primary lung fibroblasts exhibits differential gene expression patterns that are partially inversed in the same cells upon senescence entry (e.g., the negative regulation of RNAPII transcription is suppressed in the knockdown, but not in senescence). This hints towards a coordinated counter-regulation of HMGB1 effects by the rest of the program of cells entering senescence. This can be explained by the fact that the nuclear presence of HMGB1 is linked to favorable proautophagic effects that enhance cell survival and limit programmed cell death (Tang et al, 2010). This might also be a simple way to explain the strong overexpression of HMGB1 (and often also HMGB2) in various cancer types (Li et al, 2014;Zirkel et al, 2018). This overexpression seems to favor increased cell proliferation (Li et al, 2014), and it is only reasonable to assume that its targeting might be an effective anti-cancer strategy. Here, we used three lung adenocarcinoma lines to show that indeed a simple siRNA-mediated inhibition of HMGB1 suffices for replicative arrest and cell death. However, we observed that the response of each of these three lines correlated inversely to their proliferation rates, suggesting that higher rates come with stronger addiction to HMGB1 presence. Moreover, the formation (or strengthening) of SICCs in these cancer cells also aligned well with their response to HMGB1 knockdown -i.e., obvious SICC emergence signified replicative arrest.
In summary, the above allows to propose a simple model by which HMGB1 titers can be seen as a "rheostat" of cell cycle potency of a given cell. Primary proliferating cells maintain normal nuclear HMGB1 levels, cells entering senescence arrest upon nuclear depletion of HMGB1, while aberrantly proliferating cancer cells actively overexpress HMGB1 and are addicted to it for propagation ( Figure  6G). Thus, in a next step, elucidating the exact molecular dependencies of cancer cells to HMGB1 (and probably also HMGB2) nuclear overrepresentation for proliferation may lead to new ideas for combinatorial cancer treatments.

Primary cell culture and senescence markers
HUVECs from single, apparently healthy, donors (passage 2; Lonza) were continuously passaged to replicative exhaustion in complete Endopan-2 supplemented with 2% FBS under 5% CO2. Cells were constantly seeded at ~10,000 cells/cm 2 , except for late passages when they were seeded at ~20,000 cells/cm 2 . Single IMR90 isolates (I90-10 and -79, passage 5; Coriell Biorepository) were continuously passaged to replicative exhaustion in MEM (M4655, Sigma-Aldrich) supplemented with non-essential amino acids and 10% FBS under 5% CO2. Senescence-associated β-galactosidase assay (Cell Signaling) was performed according to the manufacturer's instructions to evaluate the fraction of positivelystained cells at different passages. Cell proliferation was monitored by MTT assays at different passages. In brief, ~5,000 cells are seeded in 96-well format plates in quadruplicates. On the next day, the medium is replaced with 100 ml fresh medium plus 10 ml of a 12 mM MTT stock solution (Invitrogen), and cells are incubated at 37 o C for 4 h. Subsequently, all but 25 mL of the medium is removed from the wells, and formazan dissolved in 50 mL DMSO, mixed thoroughly and incubated at 37 o C for 10 min. Samples are then mixed again and absorbance read at 530 nm. Measurements are taken at 24, 48 and 72 h post-seeding, background subtracted, and normalized to the 24 h time point. DNA methylation at six selected CpG islands (Franzen et al, 2017) was measured by isolating genomic DNA at different passages and performing targeted pyrosequencing (Cygenia GmbH). Finally, nascent DNA synthesis was monitored by EdU incorporation and subsequent labelling with Click-iT chemistry (Click-iT EdU Imaging Kit; Invitrogen). In brief, cells were incubated in 10 mM EdU for 7 h, fixed using 3.7% PFA/PBS for 15 min at room temperature, permeabilized, and labelled as per manufacturer's instructions, before imaging on a widefield Leica microscope.

Protein extraction and western blotting
For assessing protein abundance at different passages, ~4 x 10 6 cells were gently scraped off 15-cm dishes, and pelleted for 5 min at 600 x g. The supernatant was discarded, and the pellet resuspended in 100 mL of ice-cold RIPA lysis buffer (20 mM Tris-HCl pH 7.5, 150 mM NaCl, 1 mM EDTA pH 8.0, 1 mM EGTA pH 8.0, 1% NP-40, 1% sodium deoxycholate) containing 1x protease inhibitor cocktail (Roche), incubated for 20 min on ice, and centrifuged for 15 min at >20,000 x g to pellet cell debris and collect the supernatant. The concentration of the nuclear extracts was determined using the Pierce BCA Protein Assay Kit (Thermo Fisher Scientific), before extracts were aliquoted and stored at -70 o C to be used for western blotting. Resolved proteins were detected using the antisera mentioned above, plus a mouse monoclonal anti-H3K9me3 (1:200; Active motif 39286).

Chromatin immunoprecipitation (ChIP) sequencing and analysis
For each batch of ChIP experiments ~25 million proliferating cells, cultured to > 80% confluence in 15cm dishes, were crosslinked in 1.5 mM EGS/PBS (ethylene-glycol-bis-succinimidyl-succinate; Thermo) for 20 min at room temperature, followed by fixation for 40 min at 4 o C in 1%PFA. From this point onward, cells were processed via the ChIP-IT High Sensitivity kit (Active motif) as per manufacturer's instructions. In brief, chromatin was sheared to 200-500 bp fragments on a Bioruptor Plus (Diagenode; 2x 9 cycles of 30 sec on and 30 sec off at the highest power setting), and immunoprecipitation was carried out by adding 4 mg of a monoclonal HMGB1 antiserum (Developmental Studies Hybridoma Bank; PCRP-HMGB1-4F10) to ~30 mg of chromatin and rotating overnight at 4 o C in the presence of protease inhibitors. Following addition of protein A/G agarose beads and washing, DNA was purified using the ChIP DNA Clean & Concentrator kit (Zymo Research) and used in next-generation sequencing on a HiSeq4000 platform (Illumina) to obtain at least 25 million reads were obtained of both sample and its respective ''input''. Raw reads (typically 100 bp-long) were mapped to the reference human genome (hg19) using BWA (Li and Durbin, 2010), and the resulting .BAM files were processed using Picard tools (http://broadinstitute.github.io/picard/) before MACS2 software (Zhang et al, 2008) was used to identify signal enrichment over input. Thresholded HMGB1 ChIP-seq peaks per each cell type were annotated using Chipseeker (Yu et al, 2015) and are listed in Table S2A; .BAM files were used in ngs.plot (Shen et al, 2014) for plotting signal coverage over particular genomic positions for different conditions/cell types. Finally, transcription factor recognition motif enrichments within DHS footprints under HMGB1 ChIP-seq peaks were calculated using the Regulatory Genomics Toolbox (Gusmao et al, 2014). Note that all other ChIP-seq datasets used here come from previous work (Zirkel et al, 2018).

Whole-genome chromosome conformation capture (Hi-C) and TiLO analysis
Hi-C data from proliferating and senescent HUVEC were generated previously (Zirkel et al, 2018), and the HiTC Bioconductor package was used to annotate, correct data for biases in genomic features (Servant et al, 2012), and visualize 2D heatmaps with a maximum resolution of 20-kbp. For plotting insulation heatmaps and "loop-o-grams", normalized interactions values in the twenty 20-kbp bins around each HMGB1 peak were added up, normalized to the median value in each matrix and plotted provided the local maxima are higher than the third quantile of Hi-C data in the matrix. All R scripts are available on request; HMGB1-anchored loops are listed in Table S2B.
For Topologically-intrinsic Lexicographic Ordering (TiLO), we directly applied an algorithm from mathematical knot theory that makes zero assumptions about the structure, shape or number of clusters in the data (Johnson, 2014). In brief, topologically-intrinsic ordering was used to permutate the linear order of TADs (as a starting organization level in the Hi-C matrices) until a certain "robustly irreducible" topological condition is satisfied. Then, the "pinch ratio" algorithm is used (Heisterkamp and Johnson, 2013) is applied to heuristically slice the network at connections between TADs were local interaction minima are, while also considering noise in the matrices. Finally, this analysis returns a list of TADs grouped into multiple clusters in cis, also via its built-in measure for network robustness defining the end-point.

Single-cell nascent RNA sequencing and analysis
Proliferating (p. 4) and senescent HUVEC (p. 16) were washed once in an isotonic near-physiological buffer (PB) that maintains the cells' transcriptional activity and subjected immediately to the first steps of the "factory RNA-seq" protocol (Melnik et al, 2016). In more detail, cell nuclei are gently isolated using PB+0.4% NP-40, DNase I-treated at 33oC for 25 min to detach most chromatin, pelleted and washed once in ice-cold PB, before polyadenylation of nascent RNA as described (Kargapolova et al, 2017). Next, ~2,500 cells from each state were loaded onto the Chromium 10X Genomics platform for encapsulation in oil droplets and generation of barcoded cDNA libraries from individual nuclei as per manufacturer's instructions. Despite the documented 0.8% chance of capturing a cell duplet on this platform, HUVEC nuclei are particularly prone to aggregation. As a result, 494 proliferating and 129 senescent single nuclei were efficiently captured and processed. Following sequencing on a HiSeq4000 platform (Illumina), and mapping to the reference genome (hg38) using STAR (Dobin et al, 2013) and filtering via UMI-tools (Smith et al, 2017), ~45,000 and ~60,000 reads were generated per each proliferating or senescent cell, respectively. Poor quality cells were excluded (i.e. cells with <300 or >5,000 expressed genes), as were genes expressed <10 cells. This returned 1,650 robustly captured transcripts per cell on average, with >55% of reads mapping to introns or exon-intron junctions, 575 genes being expressed in at least 25% of all cells, and with the 50 most-expressed genes taking over ~22% of all sequencing reads. Mapped and filtered data were then processed and visualized using a compilation of ZINB-WaVE (ver. 1.3.4; Risso et al, 2018) and Seurat (ver. 2.3.4;Butler et al, 2018) for clustering and t-SNE visualization. ZINB-WaVE was used to create a low-dimensional representation of the gene expression matrix, and factors inferred in the ZINB-WaVE model were added as one of the low-dimensional data representations in the Seurat object in order to identify cell subpopulations via a shared nearest neighbour (SNN) modularity optimization-based clustering algorithm as applied in the FindClusters function of Seurat. Visualization was performed using the t-SNE map by applying the Rtsne function on the ZINB-WaVE output. This map was integrated into the Seurat object and used to plot gene expression. For differential gene expression analysis between the clusters, we applied a combination of ZINB-WaVE and DESeq2 (Love et al, 2014), where the posterior probabilities that counts were generated from the negative binomial count component of the ZINB-WaVE model were used as observation-level weights in the estimation of regression parameters in DESeq2 ( Van den Berge et al, 2018). Differentially-expressed genes identified via this method were filtered using a threshold of log2FC>±1 plus P-value<0.05 and are listed in Table S3A.

HMGB1 sCLIP and analysis
sCLIP was performed on ~25 million UV-crosslinked nuclei from proliferating IMR90 as previously described (Kargapolova et al, 2017) using the same the monoclonal HMGB1 antiserum (DSHB; PCRP-HMGB1-4F10) as for ChIP. Following sequencing of strand-specific libraries on a HiSeq4000 platform (Illumina), raw reads were mapped to the human reference genome (hg19). Consistent peaks were identified by overlapping intervals of peaks with a P-value <0.05 from 2 biological replicates to obtain 1773 peaks. This peak annotation was used to count reads uniquely aligned to each peak region using HTSeq, HMGB1-bound transcript coordinates were retrieved via Ensembl (GRCh37, ver. 78) and annotated using HOMER (http://homer.ucsd.edu), while Gene Ontology analysis was performed using Metascape (www.metascape.org). Finally, the final merged peak list was use for de novo motif analysis using ssHMM (Heller et al, 2017) and significantly enriched motifs were next compared to existing RBP motifs to predict proteins potentially recognizing similar sequences using Tomtom (http://memesuite.org/tools/tomtom). All HMGB1-bound mRNAs are listed in Table S1.

siRNA-mediated HMGB1 knockdown
HUVECs were seeded at ~20,000 cells/cm 2 the day before transfection. Self-delivering Accell-siRNA pools (Dharmacon) targeting HMGB1, plus a non-targeting control (NTC; fluorescently-tagged to allow transfection efficiency to be monitored), were added to the cells at a final concentration of 1 mM. Knockdown efficiency was assessed 72 h after transfection using RT-qPCR and immunofluorescence. For the knockdown in lung adenocarcinoma lines and IMR90 cells, standard siRNA transfections were carried out using RNAiMAX (Invitrogen) as previously described (Zirkel et al, 2018).

Total RNA isolation, sequencing, and analysis
Control and HMGB1-knockdown were harvested in Trizol LS (Life Technologies) and total RNA was isolated and DNase I-treated using the DirectZol RNA miniprep kit (Zymo Research). Following selection on poly(dT) beads, barcoded cDNA libraries were generated using the TruSeq RNA library kit (Illumina) and were paired-end sequenced to at least 50 million read pairs on a HiSeq4000 platform (Illumina). Raw reads were mapped to the human reference genome (hg19) using default settings of the STAR aligner (Dobin et al, 2013), followed by quantification of unique counts using featureCounts (Liao et al, 2014). Counts were further normalized via the RUVs function of RUVseq (Risso et al, 2014) to estimate factors of unwanted variation using those genes in the replicates for which the covariates of interest remain constant and correct for unwanted variation, before differential gene expression was estimated using DESeq2 (Love et al, 2014). Genes with an FDR <0.01 and an absolute (log2) foldchange of >0.6 were selected as differentially-expressed and are listed in Table S3B. For splicing analysis, a reference index on the basis of hg19 annotation was first constructed, combined with all splice sites contained in the mapped RNA-seq reads. Raw reads were then aligned using Whippet (Sterne-Weiler et al, 2018) to the constructed index in order to quantify and annotate alternative splicing events. Subsequent plots were plotted using BoxPlotR (http://shiny.chemgrid.org/boxplotr/), and GO term/pathway enrichment bar plots using Metascape (http://metascape.org/gp/index.html; Tripathi et al, 2015).

Data availability
The sequencing data generated in this study can be found in the NCBI Gene Expression Omnibus (GEO) repository under the accession number GSE98448, except HMGB1 sCLIP data that will be found under a GSE# accession number upon manuscript acceptance. Figures S1-S6 and Tables S1-S3.

Author contributions
SK, AZ, NS, NJ, YK, TG, and AM performed experiments; MN, NJ, YK, and EGG performed computational analyses; AP and AB analyzed single-cell RNA-seq data; MK and RTU provided cancer cell lines and performed knockdowns; JA and PN provided library preps and sequencing; AP conceived the study and wrote the manuscript with input from all coauthors.

Conflicts of interest
The authors declare that they have no conflict of interest.