Article MacroH2As regulate enhancer-promoter contacts affecting enhancer activity and sensitivity to inflammatory cytokines Graphical abstract Authors David Corujo, Roberto Malinverni, Juan Carrillo-Reixach, ..., Josep Manyé, Carolina Armengol, Marcus Buschbeck Correspondence carmengol@igtp.cat (C.A.), mbuschbeck@carrerasresearch.org (M.B.) In brief Corujo et al. show that the histone variants macroH2A1.2 and macroH2A2 regulate the response of cancer cells to inflammatory cytokines on the level of chromatin structure. MacroH2A-depen- dent changes in enhancer-promoter contacts are associated with increased sensitivity to tumor necrosis factor alpha and reduced sensitivity to interferon gamma. Highlights d MacroH2As facilitate the response of cancer cells to TNF- a-NF-kB but repress IFNg-STAT1 d Changes in promoter-enhancer contacts co-occur or precede changes in enhancer activity d MacroH2A2 is a more potent repressor of inflammatory response genes than macroH2A1.2 Corujo et al., 2022, Cell Reports 39, 110988 June 21, 2022 ª 2022 The Author(s). https://doi.org/10.1016/j.celrep.2022.110988 ll ll OPEN ACCESS Article MacroH2As regulate enhancer-promoter contacts affecting enhancer activity and sensitivity to inflammatory cytokines David Corujo,1,2 Roberto Malinverni,1 Juan Carrillo-Reixach,2,3 Oliver Meers,1,4 Arce Garcia-Jaraquemada,5 Marguerite-Marie Le Pannérer,1,6 Vanesa Valero,1 Ainhoa Pérez,1 Álvaro Del Rı́o-Álvarez,2,3 Laura Royo,2,3 Beatriz Pérez-González,7 Helena Raurell,7 Rafael D. Acemel,8 José M. Santos-Pereira,8 Marta Garrido-Pontnou,9 José Luis Gómez-Skarmeta,8 Lorenzo Pasquali,7 Josep Manyé,5,10 Carolina Armengol,2,3,10,* and Marcus Buschbeck1,2,11,* 1Cancer and Leukaemia Epigenetics and Biology Program, Josep Carreras Leukaemia Research Institute (IJC), Campus ICO-GTP-UAB, Badalona, Barcelona 08916, Spain 2Program for Predictive and Personalized Medicine of Cancer, Germans Trias i Pujol Research Institute (PMPPC-IGTP), Badalona, Barcelona 08916, Spain 3Childhood Liver Oncology Group, Germans Trias i Pujol Research Institute (IGTP), Badalona, Barcelona 08916, Spain 4Doctoral Programme in Biomedicine, Universitat de Barcelona, Facultat de Farmàcia i Ciències de l’Alimentació, Barcelona 08028, Spain 5IBD Research Group, Germans Trias i Pujol Research Institute (IGTP), Badalona, Barcelona 08916, Spain 6PhD Programme in Biomedicine, Universitat Pompeu Fabra, Barcelona 08003, Spain 7Department of Medicine and Life Sciences, Universitat Pompeu Fabra, Barcelona 08003, Spain 8Centro Andaluz de Biologı́a del Desarrollo (CABD), CSIC-Universidad Pablo de Olavide-Junta de Andalucı́a, Sevilla 41013, Spain 9Pathology Department, Hospital Vall d’Hebron, Barcelona 08035, Spain 10Liver and Digestive Diseases Networking Biomedical Research Centre (CIBEREHD), Madrid 28029, Spain 11Lead contact *Correspondence: carmengol@igtp.cat (C.A.), mbuschbeck@carrerasresearch.org (M.B.) https://doi.org/10.1016/j.celrep.2022.110988 SUMMARY MacroH2A histone variants have a function in gene regulation that is poorly understood at the molecular level. We report that macroH2A1.2 and macroH2A2 modulate the transcriptional ground state of cancer cells and how they respond to inflammatory cytokines. Removal of macroH2A1.2 and macroH2A2 in hepatoblastoma cells affects the contact frequency of promoters and distal enhancers coinciding with changes in enhancer activity or preceding them in response to the cytokine tumor necrosis factor alpha. Although macroH2As regulate genes in both directions, they globally facilitate the nuclear factor kB (NF-kB)-mediated response. In contrast, macroH2As suppress the response to the pro-inflammatory cytokine interferon gamma. MacroH2A2 has a stronger contribution to gene repression than macroH2A1.2. Taken together, our results suggest that macroH2As have a role in regulating the response of cancer cells to inflammatory signals on the level of chromatin structure. This is likely relevant for the interaction of cancer cells with immune cells of their microenvironment. INTRODUCTION MACROH2A2 (previously H2AFY2) encodes macroH2A2 (Posa- vec et al., 2013). A large number of loss-of-function studies sug- Histone variants can replace replication-coupled histones in the gest that macroH2A proteins promote and stabilize differenti- nucleosome and endow chromatin regions with specific ated cell states (Buschbeck and Hake, 2017). MacroH2A biochemical and biophysical characteristics (Martire and Banas- histone variants promote the differentiation of stem cells (Barrero zynski, 2020). Compared with replication-coupled histones, his- et al., 2013a; Creppe et al., 2012) and inhibit the induction of plu- tone variants differ in sequence, expression timing, and mRNA ripotency by nuclear transfer or by transduction of Yamanaka processing (Buschbeck and Hake, 2017). Some of the most sub- factors (Barrero et al., 2013b; Gaspar-Maia et al., 2013; Pasque stantial sequence changes are found in three histone H2A vari- et al., 2011). In cancer, the role of macroH2As is context- and ants known as macroH2As that share a peculiar tripartite struc- isoform dependent (Ghiraldini et al., 2021). Reported observa- ture consisting of the histone fold, a linker, and a macrodomain tions include a tumor-suppressive function of macroH2A2 in (Chakravarthy et al., 2005; Timinszky et al., 2009). Through alter- melanoma (Kapoor et al., 2010) and an association of mac- native splicing, the MACROH2A1 gene (previously H2AFY) gives roH2A1.1 with poor survival in triple-negative breast cancer (Lav- rise to the isoforms macroH2A1.1 and macroH2A1.2, while igne et al., 2014). Cell Reports 39, 110988, June 21, 2022 ª 2022 The Author(s). 1 This is an open access article under the CC BY-NC-ND license (http://creativecommons.org/licenses/by-nc-nd/4.0/). ll OPEN ACCESS Article The molecular basis for how macroH2A histone variants genes, we describe how macroH2A sensitivity is associated mediate these physiological and pathophysiological functions with differences in chromatin architecture and enhancer activity. is not known. The field generally assumes that it is mediated through an influence on gene regulation. However, how mac- RESULTS roH2A proteins impact on gene regulation is not understood. So far, we know that the presence of macroH2A on transcribed Loss of macroH2A1.2 and macroH2A2 globally regions is associated with a repressed state (Gamble et al., modulates the transcriptional ground state 2010). This is explained by the fact that reincorporation of mac- The loss of macroH2A histone variants has drastic conse- roH2A relative to replication-coupled H2A is disfavored after quences for nuclear organization and global three-dimensional passing of the polymerase complex in transcribed regions (Sun (3D) heterochromatin architecture (Douet et al., 2017). To gain et al., 2018). On inactive genes, the presence of macroH2A pro- a better understanding of whether these changes are associated vides a layer of repression in addition to histone deacetylation with changes in gene regulation, we have decided to study the and DNA methylation (Hernández-Muñoz et al., 2005). However, loss of function of macroH2A on the transcriptional ground state deletion of macroH2As affects gene expression in both direc- and the response to innate immune signals. Specifically, we have tions (Chen et al., 2014; Hurtado-Bagès et al., 2020), suggesting used the hepatoblastoma-derived cell line HepG2 (López-Ter- that the role of macroH2As is more complex than those of purely rada et al., 2009) that expresses macroH2A1.2 and macroH2A2 repressive chromatin components. In addition, macroH2A was but not macroH2A1.1 (Figure 1A). Short hairpin RNA (shRNA)- found to be required for the proper induction of gene expression mediated knockdown allowed us to efficiently deplete both mac- in response to different signals (Creppe et al., 2012; Dell’Orso roH2A proteins without notable compensation by macroH2A1.1 et al., 2016; Ouararhni et al., 2006). One of the best-studied ex- (Figure 1A). We reduced the serum concentration from the cell amples is the macroH2A1.1-dependent response of fibroblasts culture media to minimize the effect of stimuli coming from the to a cocktail of cytokines secreted by senescent cells, which is complex composition of fetal bovine serum and analyzed the essential for the propagation of the paracrine signal (Chen transcriptome in control and macroH2A-deficient HepG2 cells et al., 2015). Moreover, macroH2As have an important role in nu- treated or not with TNF-a. Principal-component analysis clear organization and are required for maintaining the charac- confirmed the similarity of biological replicates and that both teristic three-dimensional architecture of heterochromatin depletion of macroH2A and TNFa treatment affected the tran- (Douet et al., 2017), a function that can be mediated by the linker scriptome (Figure 1B). Focusing first on the untreated condition, (Kozlowski et al., 2018). At present, it is not clear how the function we found that a substantial number of genes were differentially in chromatin architecture ties in with transcriptional regulation. expressed when comparing macroH2A-depleted cells with their To gain a better understanding of how macroH2A affects gene controls (Figure 1C). The numbers of those up- and downregu- regulation in cancer, we have studied the transcriptional lated as well as their effect sizes were in a similar range. Top-en- response to signals from innate immune cells. We have taken riched Gene Ontologies of upregulated genes were related to cell two different approaches. To study the response to a specific motility, metabolism, and the inflammatory response, while signal, we have chosen tumor necrosis factor alpha (TNF-a), those of downregulated genes were related to morphogenesis, which is one of the major pro-inflammatory cytokines involved development, and adhesion (Figures 1D and 1E). Pathway anal- in the innate immune response (Zelová and Ho sek, 2013). Acti- ysis indicated downregulation of the PI3K pathway and upregu- vated macrophages are the primary source of TNF-a, but other lation of several pathways including WNT and JAK-STAT (Fig- cells of the immune system can also secrete it. The role of ure 1F). We observed similar macroH2A-dependent differences TNF-a in cancer is complex. On one hand, TNF-a is an important in gene expression when comparing macroH2A-depleted and effector of cytotoxic T cells and natural killer cells, but on the control cells treated with TNF-a (Figures S1A–S1C). 60% of up- other, it can favor tumor growth and immune evasion (Balkwill, regulated and 53% of downregulated genes overlapped in un- 2009). Signal transduction after TNF-a is mainly mediated by treated and TNFa-treated cells, while the others were only differ- the canonical pathway involving the IKK complex and nuclear entially expressed in one of the two conditions (Figure S1D). In factor kB (NF-kB) transcription factors (Taniguchi and Karin, conclusion, the stable depletion of macroH2A1.2 and mac- 2018). The NF-kB family is composed of 5 DNA-binding proteins roH2A2 in hepatoblastoma cells changes their transcriptional that include RelA and RelB and that can form different homo- and ground state. heterodimeric transcription factors (Perkins, 2012). To study the impact of macroH2A on the inflammatory response in a more Lack of macroH2A affects the transcriptional response complex and natural environment, we have studied cells grown to TNF-a as xenografts in mice possessing intact innate immunity and As expected, TNF-a prominently induced genes representing analyzed previously generated datasets on primary patient the inflammatory response, while repressed genes were related samples. to metabolic pathways (Figures 2A and S2A–S2C). 663 (43.7%) Here, we describe how the levels of macroH2As modulate the of the total of 1,448 TNF-a-induced genes were further differen- ground state of transcription in a way that determines the quality tially expressed in macroH2A-depleted cells (Figure 2B). This of the response to inflammatory signals. Specifically, we present proportion was significantly higher than for TNF-a-insensitive results showing that macroH2A1.2 and macroH2A2 affect the genes, defined as those that were readily detected but neither transcriptional responses to TNF-a and interferons in a complex repressed nor induced by TNF-a treatment. 462 of these 633 and, to some extent, opposing manner. For TNF-a-response genes were significantly less expressed in macroH2A-deficient 2 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS ) KD 1/2 ) 2 A TL (D m H l (C s h ontro FLAG-mH2A C A B C sh 1.1 1.2 2 DKD vs CTL (untreated) mH2A1.2 40kDa 1633 1580 PC2: 32% variance DOWN UP −log(adjusted pv) 10 200 CTL mH2A1.1 40kDa unt. DKD MACROH2A2 mH2A2 40kDa 0 CTL 100 +TNFα MACROH2A1 FLAG 40kDa DKD −10 H3 15kDa 0 −20 −10 0 10 −8 −4 0 4 8 PC1: 55% variance log2FC D E F UP DOWN cell morphogenesis involved in PROGENy monocarboxylic acid metabolic process differentiation Pathway wound healing heart development PI3K Score small molecule catabolic process cell part morphogenesis WNT 1 response to wounding ribonucleoprotein complex biogenesis 0.5 p53 0 positive regulation of cell motility cellular component morphogenesis n. genes Androgen −0.5 organic hydroxy compound metabolic regulation of cell projection process organization 80 JAK−STAT −1 positive regulation of locomotion ncRNA metabolic process 90 positive regulation of cellular brain development CTL DKD component movement 100 inflammatory response regulation of plasma membrane bounded cell projection organization positive regulation of hydrolase regulation of cell adhesion activity 0 5 10 15 0 5 10 15 −log10(adjusted pv) −log10(adjusted pv) Figure 1. Loss of macroH2A modulates the transcriptional ground state of HepG2 cells (A) Total cell extracts from double-knockdown (DKD) HepG2 cells stably transduced with shRNAs for macroH2A1 and macroH2A2 and control (CTL) cells were analyzed by immunoblotting using the indicated antibodies. Polyclonal DKD cells stably expressing equal amounts of FLAG-tagged macroH2A proteins were used as reference. All samples were detected on the same blotted membrane (per antibody), and white spaces indicate a separation where lanes of the blot have been re-arranged. The antibody used for macroH2A1.2 also detects macroH2A1.1, but as it is not expressed, the signal can be interpreted as macroH2A1.2. (B) Principal-component (PC) plot of RNA sequencing (RNA-seq) samples from CTL and DKD cells treated or not with TNF-a at 20 ng/mL for 24 h. Three different replicates were seeded and treated on different days. (C) Volcano plot of differentially expressed genes between DKD and CTL in untreated conditions plotted as log2 fold change (log2FC; Benjamini-Hochberg adjusted p value < 0.05). (D) Top enriched Gene Ontologies in upregulated genes. (E) Top enriched Gene Ontologies in downregulated genes. (F) Inferred pathway activity using PROGENy. cells (Figure 2C). A smaller proportion of TNF-a-induced genes tested upregulated genes, while the individual knockdown of (171 of 633) were further upregulated by loss of macroH2A. Us- macroH2A1.2 only affected one (Figures S3C and S3D). In ing independent samples and semiquantitative PCR, we have contrast to the double knockdown, the individual knockdowns confirmed the bidirectional influence of macroH2A depletion did not downregulate any of the 5 tested genes and in some on a panel of TNF-a-induced genes (Figures 2D and 2E). Many cases even caused the opposite effect (Figure S3D). of these genes were also differentially expressed in untreated Taken together, here we identify an enrichment between TNF- conditions. Next, we analyzed the same panel of genes in a a- and macroH2A-regulated genes. The influence of macroH2As rescue setting in which we individually expressed exogenous is bidirectional and requires the loss of both macroH2A1.2 and macroH2A1.2 and macroH2A2 in macroH2A-depleted cells. macroH2A2. From the immunoblot, we estimated that the achieved levels were close to the endogenous levels of macroH2A1.2 and mac- MacroH2A affects the activity and TNF-a sensitivity of roH2A2 in parental HepG2 cells (Figure 1A). Interestingly, these some enhancers levels of the individual macroH2A proteins were sufficient to To get a better understanding of how macroH2A affects gene partially rescue a subset of upregulated, but not downregulated, regulation on the molecular level, we took a closer look at the genes (Figures 2F, S3A, and S3B). In this, macroH2A2 was more extended genomic loci and in particular enhancers of three potent than macroH2A1.2. In line with this observation, the indi- genes that were sensitive to both macroH2A depletion and vidual knockdown of macroH2A2 upregulated five out of five TNF-a treatment. These included DKK1, whose expression Cell Reports 39, 110988, June 21, 2022 3 ll OPEN ACCESS Article A +TNFα vs untreated (CTL) B C TNFα-induced genes Repressed Induced 200 462 171 1170 OPTN 1448 TNFα-Insensitive TNFα-Induced TNFα-Repressed DOWN UP −log(adjusted pv) −log(adjusted pv) EFNB2 200 (13225) (1488) (1170) UBD 17.8% 43.7% 44.4% TTN KLHL13 PIK3AP1 NFKBIA DKD vs CTL TTC39B IFIH1 100 100 RELB UP ITIH5 SYTL2 PCDH9 CAVIN2 CLDN19 CCL20 DOWN PLAU FBLIM1 HTR1D ITGAM STK17B CXCL8 TNF Not affected RBP1 SLC6A14 DPEP1 EREG DKK1 0 log2(OR) = 1.85 log2(OR) = 1.89 0 −8 −4 0 4 8 −6 −4 −2 0 2 4 6 log2FC * log2FC DKD vs CTL (+TNFα) * D E F TNFα TNFα DOWN UP 8 * 20 * 150 - + - + Rel. expression Rel. expression CTL CTL * 15 100 % Rescue 6 * DKD DKD * DKD vs CTL 4 10 * 50 effect size * * * * * * * * 0 2 * 5 * * * * * * * * * * * * * * * * * * Non−significant * * −50 Significant 0 0 −100 ER U EG R 9 ST 1 FB 7B 1 5 EF N SY 2 2 1 C P1 C 19 TT IN2 KL 9B D 3 K3 1 H P1 C D 14 H BP M IH B TL IH L1 PI K A TT SL 1 N 6A K1 N PE 3 K A D TR PL LI IT IF AV H 2 2 2 2 C LD PC 1. 2A 1. 2A D 2A 2A H H −m −m H H −m −m AG AG AG AG FL FL FL FL Figure 2. An enrichment between TNF-a- and macroH2A-regulated genes (A) Volcano plot of the log2FC of gene expression comparing TNF-a-treated CTL cells with untreated CTL cells. All genes with an adjusted p value <0.05 were considered to be differentially expressed. (B) Pie charts showing the enrichment of macroH2A-regulated genes among TNF-a-induced and -repressed genes relative to TNF-a-insensitive genes (detected at the same levels in treated and untreated conditions). A Fisher’s exact test on a 2 3 2 contingency table has been used to make the indicated comparisons (*adjusted p value < 0.05; log2(OR), log2 odds ratio). (C) Volcano plot of the log2FC of 633 macroH2A-sensitive, TNF-a-induced genes comparing CTL and DKD in TNF-a-treated conditions. The number of induced and repressed genes is indicated. (D) Relative expression of a panel of TNF-a-induced genes downregulated in DKD cells compared with CTL cells assayed by qRT-PCR. TNF-a treatment was 20 ng/mL for 24 h (n = 3 biological replicates). A two-tailed Student’s t test was used to make the indicated pairwise comparison. *p < 0.05. (E) Same as in (E) for TNF-a-induced genes upregulated in DKD cells compared with CTL cells. (F) FLAG-tagged macroH2A1.2 and macroH2A2 were individually and exogenously expressed in DKD cells. The resulting protein levels are shown in Figure 1A. The effect on the mRNA levels of 20 macroH2A and TNF-regulated genes (gene panels from D and E) was plotted for each gene relative to the difference between DKD and CTL (set to 100%). Significant rescues are indicated by dark bullets (p < 0.05, two-tailed Student’s t test, n = 3 biological replicates). The full data used to generate this summary are shown in Figures S3A and S3B. See also Figures S1–S3. was increased already by macroH2A depletion and further by and the regulatory annotation of the latest version of the TNF-a treatment in a dose- and IKK-dependent manner, leading GeneHancer database (Fishilevich et al., 2017). Taken together, to increased secretion of WNT-antagonistic DKK1 protein these datasets suggest that the DKK1 gene is regulated by (Figures S4A–S4C). Using our previously generated chromatin several 30 enhancers, of which the strongest coincides with the immunoprecipitation sequencing (ChIP-seq) data (Douet et al., MBL2 gene located at a distance of 0.45 Mb. Knockdown of 2017) and a peak-calling algorithm optimized for broader re- macroH2A increased the expression of MBL2 in a similar manner gions, we found that both macroH2A1.2 and macroH2A2 are en- as DKK1 but did not affect the much closer 50 gene PRKG1 or its riched around the DKK1 gene and several 30 prime regions (Fig- antisense transcript PRKG1-AS (Figure S4D). According to an ure 3A). This broad signal is typical for macroH2A histone available ENCODE dataset in HepG2 cells, the NF-kB transcrip- variants and makes the functional interpretation of their genomic tion factor RELA has a single annotated binding site at the DKK1 profile challenging (Douet et al., 2017; Gamble et al., 2010; Rat- gene promoter. Next, we tested whether depletion of macroH2A nakumar et al., 2012). Here, we have used epic2 (Stovner and would affect the level of H3K27ac on 30 enhancers as a proxy for Sætrom, 2019), a peak caller particularly aimed at identifying their activity. We found that several of the tested enhancers, broad enrichment regions in diffuse ChIP-seq signals, to define including those around the MBL2 gene, increased their basal ac- domains of significant macroH2A1.2 and macroH2A2 enrich- tivity and their sensitivity to TNF-a (Figure 3B). On a second up- ment. To identify and visualize regulatory enhancers of the regulated gene, PIK3AP1, macroH2A depletion affected the DKK1 gene, we have used the ROADMAP chromatin state anno- basal activity of two out of three tested enhancers located at 3’ tation for HepG2 (Roadmap Epigenomics Consortium et al., (Figures 3C and 3D). EREG is another TNF-a response gene 2015), the histone H3 K27 acetylation profile it is based on, but downregulated in macroH2A-depleted cells. Two 30 4 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS A DKK1 C PIK3AP1 E EREG 40 40 40 mH2A1.2 0 0 0 40 40 40 mH2A2 0 0 0 15 15 15 RELA 0 0 0 H3K27Ac 100 100 100 Enhancers 0 0 0 Amplicons Enhancer 1 2 34 56 7 8 910 Associations DKK1 MBL2 TLL2 PIK3AP1 LCOR SLIT1 MTHFD2L EREG AREG BTC PRKG1 PRKG1−AS1 TM9SF3 EPGN chr10 chr10 chr4 54Mb 54.2Mb 54.4Mb 54.6Mb 98.2Mb 98.4Mb 98.6Mb 98.8Mb 75Mb 75.2Mb 75.4Mb 75.6Mb 75.8Mb B D F DKK1 PIK3AP1 EREG 20 * 20 CTL * CTL CTL DKD DKD DKD 15 10 15 * * * ChIP / H3 ChIP / H3 ChIP / H3 * * 10 10 * 5 * * 5 * 5 * * * * * * * * 0 0 0 TNFα - + - + - + - + - + - + - + - + TNFα - + - + - + - + - + - + - + - + TNFα - + - + - + - + - + - + IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac IgG K27Ac 1 2 3 4 neg. 5 6 7 8 9 10 Figure 3. MacroH2A depletion affects the activity of enhancers of both up- and downregulated genes (A) Genome browser view of the extended DKK1 locus showing enrichment profiles for macroH2A1.2 macroH2A2, RELA, and H3K27ac in HepG2 cells. Peaks and regions of significant enrichment are indicated for macroH2A1.2, macroH2A2, and RELA below the enrichment profile. Enhancer regions are indicated using the annotation provided by the ROADMAP chromatin state of HepG2. Orange arches indicate functional gene-enhancer associations from a large compendium of cell types annotated in the GeneHancer database. The numbers indicate the position of amplicons used for ChIP-qPCR analysis in (B). (B) ChIP-qPCR data for H3K27ac in the regions indicated in (A). Values have been normalized to H3 ChIP signal. Immunoglobulin G (IgG) is used as a negative CTL. Cells were cultured in the absence or presence of TNF-a at 20 ng/mL for 24 h (n = 4 biological replicates, *p < 0.05, two-tailed Student’s t test on the indicated comparisons, K27Ac stands for H3K27Ac). (C) Same as in (A) but for the extended locus for the upregulated gene PIK3AP1. (D) Same as in (B) but for loci indicated in (C) and a negative CTL region (n = 3 biological replicates, *p < 0.05, two-tailed Student’s t test on the indicated comparisons, K27Ac stands for H3K27Ac). (E) Same as in (A) but for the downregulated gene EREG. (F) Same as in (B) but for the loci indicated in (E) (n = 3 biological replicates, *p < 0.05, two-tailed Student’s t test on the indicated comparisons, K27Ac stands for H3K27Ac). See also Figure S4. enhancers lost their sensitivity to TNF-a in macroH2A-depleted we have analyzed DKK1-MBL2 and EREG as examples of genes cells (Figures 3E and 3F). This was not the case for a third and that are up- or downregulated by macroH2A depletion after TNF- neighboring enhancer that we tested and that happened to be a treatment. We found that the DKK1 gene was positioned at the bound by RELA in parental HepG2 cells. Importantly, the level 50 end of a topologically associated domain (TAD) that spans a of macroH2A1.2 and macroH2A2 enrichment at tested en- little more than 0.5 Mb and the MBL2 gene at the 30 end (Fig- hancers was low to intermediate and not affected by TNF-a ure 4A). This was indicated by contact probability maps from treatment (Figures S4E and S4F). HiC data generated by ENCODE (Yardımcı et al., 2017) and Taken together, these results suggest that the absence and our own HiChIP data generated with an antibody directed presence of macroH2A1.2 and macroH2A2 can affect both the against H3 modified by trimethylation on lysine 4 (H3K4me3), a activity of enhancers in their ground state as well as the sensi- mark of active promoters. The H3K4me3 HiChIP data further tivity of enhancers to incoming inflammatory signals. showed that the DKK1 promoter primarily established contacts with the downstream region inside the same TAD, including the Changes in enhancer activity and sensitivity are regulatory distal enhancers located at the MBL2 gene. Next, associated with macroH2A-dependent changes in we used semiquantitative UMI-4C to probe for macroH2A- chromatin architecture dependent changes in contact frequency of this regulatory From our previous study, we know that loss of macroH2A has a region. In macroH2A-depleted cells, we observed increased drastic impact on large-scale chromatin organization in the nu- contact with regions containing the DKK1 promoter and other cleus. To get a better understanding of whether this also occurs macroH2A-sensitive enhancers 65 kb 30 of DKK1 (Figure 4A, bot- in macroH2A-sensitive inflammatory response genes, we turned tom panel). In control cells, these regions were enriched in mac- to chromosome-conformation-capture methods. Specifically, roH2A1.2 and macroH2A2. H3K4me3 HiChIP further showed Cell Reports 39, 110988, June 21, 2022 5 ll OPEN ACCESS Article DKK1-MBL2 EREG A HiC (HepG2 untreated) B (HepG2 untreated) HiC H3K4me3 H3K4me3 HiChIP HiChIP chr10 54Mb 54.2Mb 54.4Mb 54.6Mb 54.8Mb chr4 75Mb 75.2Mb 75.4Mb 75.6Mb 75.8Mb PRKG1 AS1 MTHFD2L PRKG1 AREG DKK1 MBL2 EPGN EREG BTC UMIs normalized trend UMIs normalized trend 15 30 10 CTL * CTL * * DKD 20 * DKD 5 10 0 mH2A1 0 mH2A1 mH2A2 mH2A2 Enhancers Enhancers 1 2 34 8 9 10 chr10 54Mb 54.2Mb 54.4Mb 54.6Mb 54.8Mb chr4 75Mb 75.2Mb 75.4Mb 75.6Mb 75.8Mb Figure 4. MacroH2A-regulated loci display changes in chromatin architecture (A) Three-dimensional architecture of the DKK1-MBL2 locus. The top panel shows the contact frequency maps of HepG2 HiC data generated by ENCODE. The middle panel shows the contact frequency map of the same region using H3K4me3 HiChIP data generated by us in untreated HepG2 cells. The bottom panel shows UMI-4C data that was used to quantify changes in contact frequency in untreated CTL and DKD HepG2 cell lines. The UMI-4C profile contains aggregated sequencing data from four biological replicates. The black arrow indicates the position of the viewpoint at the predicted regulatory region at the MBL2 gene locus. As in Figure 3A, enhancers and macroH2A-enriched regions are annotated. Red arrowheads indicate the regions analyzed by ChIP-qPCR in Figure 3B. The dotted line indicates the position of the DKK1 promoter. (B) Same as in (A) but for the extended EREG locus. The UMI-4C profile contains aggregated sequencing data from three biological replicates. that the promoter of the downregulated gene EREG established study of cancer and autoimmune diseases (Bancroft and Kelly, contacts with regions both at 30 and 50 (Figure 4B). In the absence 1994) https://www.jax.org/news-and-insights/2004/october/nod- of macroH2A, two regions, one 50 and one 30 , reduced contact cb17-prkdc-scid-j-strain-an-ideal-model-for-cancer-andautoimmune- frequency with the EREG promoter. Interestingly, both regions disea#. While we did not observe any difference in growth or overall contained annotated regulatory enhancers and were enriched tissue composition (Figures S5A and S5B), we confirmed that the in macroH2A1.2 and macroH2A2. repression of MACROH2A1 and MACROH2A2 transcripts was sta- Taken together, these results show that depletion of mac- ble and still prominent in xenografts after resection (Figure 5A). roH2A1.2 and macroH2A results in changes in chromatin archi- Next, we performed mRNA sequencing of five tumors of each tecture and promoter-enhancer contacts. On the upregulated experimental group. Principal-component analysis indicated that gene DKK1, increased expression and enhancer activity is asso- the overall expression pattern of xenografts was similar within the ciated with increased contacts. In the case of EREG, changes in same group but greatly differed between groups (Figure S5C). Us- chromatin architecture, specifically reduced contacts, preceded ing an adjusted p value of 0.05 as cutoff, we identified 2,113 genes changes in gene expression and enhancer activity that were only to be downregulated and 2,293 to be upregulated when comparing observed after TNF-a treatment (Figure 3F). macroH2A-deficient xenografts with their control counterparts (Figure 5B). Both MACROH2A1 and MACROH2A2 were among MacroH2A has a complex function in inflammatory the most downregulated genes, thus validating the performance responses globally facilitating NFkB responses while of our mRNA sequencing. Enriched Gene Ontologies included ri- inhibiting responses to interferon gamma (IFNg)-STAT1 bosomal pathways in downregulated genes and catabolic pro- Our results so far have established that the presence and absence cesses in upregulated genes (Figure S5D). Gene set enrichment of macroH2A1.2 and macroH2A2 in HepG2 cells affects genes analysis indicated an enrichment of genes related with the inflam- involved in the inflammatory response to TNF-a when treated in matory response (Figure 5C) among the upregulated genes. In culture. To test whether macroH2A affects cells exposed to cyto- addition, genes related to hypoxia were enriched among the upre- kines in a physiologic setting, we have decided to perform xeno- gulated genes (Figure S2E), which resonated previous findings grafts. From one of our previous studies, we know that HepG2 relating macroH2A1 isoforms to redox metabolism in breast cancer cells lacking macroH2A1.2 but retaining macroH2A2 are able to cells (Dardenne et al., 2012). To gain a better understanding of the grow as xenografted tumors in immunosuppressed mice (Lo Re pathways mediating the observed transcriptional changes, we et al., 2017). Specifically, we injected macroH2A-deficient have used a method to infer transcription factor activity (Alvarez HepG2 cells and their control counterparts into partially immuno- et al., 2016; Garcia-Alonso et al., 2019). Among the downregulated suppressed NOD.CB17-PRKDCSCID mice (Figure 5A). These mice regulons in macroH2A-deficient xenografts was the NF-kB tran- lack functional B and T cells but retain their innate immune system, scription factor RELB (Figure 5D), which is part of the canonical making them a useful and frequently used model system for the signaling pathway of TNF-a (Zelová and Ho sek, 2013). This is in 6 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS A HepG2 Resected xenografts B C CTL X.CTL DKD X.DKD Inflammatory response Enrichment score 300 2113 2293 Enrichment MACROH2A1 MACROH2A2 MACROH2A2 0.4 profile DOWN UP −log(adjusted pv) 2.0 2.0 Hits * 200 MACROH2A1 0.2 Rel. expression 1.5 1.5 0 * 100 DKK1 1.0 1.0 0.5 0.5 0 UP log2FC X.DKD vs X.CTL DOWN −8 −4 0 4 8 0.0 0.0 NES = 1.67 log2FC X.DKD vs X.CTL pv = 0.007 FDR = 0.198 KD KD TL TL C C D D X. X. X. X. TNFα IFNγ IFNγ D E NFκB (RELB) IFNγ response genes PRDM14 SREBF2 SREBF1 POU2F2 TCF7L2 STAT5A FOXO4 CEBPA FOSL1 EPAS1 GATA2 STAT6 STAT4 STAT1 USF2 ETV4 ELK1 TP63 MITF ATF6 IRF1 WT1 * YY1 AR 75 Rel. expression X.CTL * CTL 50 DKD X.DKD 25 * * * 0 - + - + - + - + - + IFNγ Regulator activity 1 F1 0 0 7 AT X6 L1 I2 IR IF ST XC -2 -1 0 1 2 D D C Figure 5. MacroH2A proteins affect the quality of the inflammatory response in vivo (A) Relative expression of macroH2A1- and macroH2A2-encoding transcripts in resected xenografts of macroH2A-deficient HepG2 cells (X.DKD) and CTL cells (X.CTL) was assayed by qRT-PCR. Rhombi indicate the mean values and thick horizontal lines the median values per group (n = 10 tumors per condition, a two- tailed Student’s t test was used to make the indicated pairwise comparison, *p < 0.05). (B) Volcano plot of differentially expressed genes between X.DKD and X.CTL in untreated conditions plotted as log2FC (Benjamini-Hochberg adjusted p value < 0.05). The number of up- and downregulated genes are given. (C) Gene set enrichment profile using all significantly differentially expressed genes ranked by their log2FC. The M5932 - HALLMARK_INFLAMMATORY_RESPONSE gene set from the Molecular Signatures Database is shown. NES, normalized enrichment score; FDR, false discover rate; pv, p value. (D) Gene expression count matrices were used to infer changes in regulon activity using VIPER and the DoRothEA regulon sets. Mediators of IFNg (STAT1 and IRF1) and TNF-a signaling (NF-kB) are indicated. Values in the heatmap are scaled by column. (E) Expression of STAT1, IRF1, and three IFNg response genes assayed by qRT-PCR in HepG2 CTL and DKD cultured in the absence or presence of IFNg at 100 ng/mL for 48 h (n = 3 biological replicates, a two-tailed Student’s t test was used to make the indicated pairwise comparison, *p < 0.05). See also Figure S5 and Table S1. line with our observation that loss of macroH2A led to a substan- response in a more natural environment of cytokines. The tran- tially larger number of TNF-a-induced genes being down- than up- scriptomic analysis points toward a complex role of macroH2A regulated (Figure 2B). Upregulated regulons included IFN regulato- in inflammatory signaling. The net output of macroH2A’s pres- ry factor 1 (IRF1) and STAT1, two well-known mediators of another ence is, on one hand, the promotion of TNF-a-NF-kB signaling type of pro-inflammatory signaling (Figure 5D). Specifically, IRF1 and, on the other, a suppression of IFNg-STAT1 signaling. and STAT1 are the main transcriptional mediators of IFNs, which are central cytokines of the innate immune response (Feng et al., MacroH2A1 expression levels are related to 2021). To test an influence of macroH2A on IFN signaling, we inflammatory transcription factor activity in turned again to cell culture and analyzed known IFNg response hepatoblastoma tumor samples genes in macroH2A-deficient HepG2 and control cells. The The association of macroH2A variants with cancer is context- expression of the transcription factors STAT1 and IRF1 as well and isoform dependent. In cancers of the liver tissue, both can- as the IFNg-response genes CXCL10, DDX60, and IFI27 was cer-promoting and -repressing functions have been attributed to strongly increased in cells lacking macroH2A upon treatment macroH2A histone variants (Hsu et al., 2021). Hepatoblastoma is with IFNg (Figure 5E). the most common type of pediatric liver cancer and usually de- Taken together, the growth of macroH2A-deficient and control velops in infants and young children during the first 3 years of life HepG2 cells as xenografts allowed us to investigate their (Schnater et al., 2003). We have recently described a detailed Cell Reports 39, 110988, June 21, 2022 7 ll OPEN ACCESS Article HB low vs high A B MACROH2A1 C HB HepG2 14 * 2078 1150 low vs high DKD vs CTL 20 DOWN UP Expression (log2) −log(adjusted pv) 12 * MACROH2A1 15 NT 10 HB 10 2391 838* 2233 5 8 0 −4 0 4 1 2 2A 2A log2FC H H RO RO AC AC M M D E HB F Pan-Cancer (TCGA) low vs high 5.0 64 283 STAT1 RELA HNF4A 100 100 log2FC HB low vs high 2.5 STAT1 % of samples 75 75 CEBPA R2 = 0.193* Regulon TF activity MYC 0.0 50 50 UP RELA DOWN Chi-Square E2F1 25 25 −2.5 pv < 2.2e-16 E2F4 HIF1A 0 0 −5.0 378 113 −4 −2 0 2 4 6 1 2 3 4 1 2 3 4 −5.0 −2.5 0.0 2.5 5.0 NES log2FC HepG2 DKD vs CTL (unt.) (FDR < 0.05) low high low high MACROH2A1 expression Figure 6. MacroH2A1 expression levels are related to inflammatory transcription factor activity in hepatoblastoma primary tumor samples (A) Boxplot showing the log2 expression values of MACROH2A1 and MACROH2A2 for 32 hepatoblastoma (HB) tumor samples and matched non-tumor (NT) tissue samples. An unpaired t test with Welch correction was used to compare the two groups, *p < 0.05. (B) Volcano plot of differentially expressed genes between the tercile of HB tumors with lowest and highest MACROH2A1 expression (log2FC; Benjamini-Hoch- berg adjusted p value < 0.05). For tercile separation, please see Figure S6A. (C) Venn diagram showing the significant overlap of differentially expressed genes in low versus high MACROH2A1 and DKD versus CTL (both datasets have been reduced to the genes that have been detected in both experiments). A Fisher’s exact test was used to evaluate the significance of the overlap, *p < 0.05. (D) Scatterplot of log2FC values of the differentially expressed genes between low versus high MACROH2A1 and DKD versus CTL. The blue line represents the regression line obtained by fitting a linear model on the data points, and R2 is the coefficient of determination of the resulting linear regression model (p < 2.16 3 10-16). A chi square test of independence has been performed on the categories defined by the four quadrants, and the p value is reported. (E) Regulon activity as annotated in DoRothEA has been inferred with VIPER, and normalized enrichment scores (NESs) with a false discovery rate (FDR) < 0.05 are plotted. (F) The samples from the TCGA Pan-Cancer cohort with up- or downregulated inferred activities of the transcription factors STAT1 (n = 1,000) and NF-kB (RELA; n = 256) by the RABIT algorithm were ranked in four quartiles according to their MACROH2A1 expression values. The proportion in percentage of samples with up- or downregulated activities for each transcription factor is plotted for each quartile. See also Figure S6. molecular characterization of one of the largest cohorts of pa- ure 6C). The correlation between both datasets was positive, as tients with hepatoblastoma (HB) (Carrillo-Reixach et al., 2020). the large majority were differentially expressed in a concordant This analysis included the transcriptome-wide analysis of 32 manner (Figure 6D). However, this was not the case for paired tumor and adjacent normal liver tissues. While the data MACROH2A2 (Figures S6B–S6E). Pathway analysis of the HB did not allow us to resolve differences between the two mac- tumor samples indicated that low MACROH2A1 levels were roH2A1 splice isoforms, it showed that total levels of macroH2A1 associated with increased STAT1 signaling and reduced activity and macroH2A2 transcripts were upregulated in tumor tissues of the NF-kB transcription factor RELA (Figure 6E). Low compared with adjacent normal tissue and that MACROH2A1 MACROH2A2 levels were also associated with increased was more prominently expressed than MACROH2A2 (Figure 6A). STAT1 signaling but, in contrast to MACROH2A1, also with We have separated the tumor samples in low, medium, and high higher levels of RELA activity (Figure S6F). To test whether this MACROH2A1-expressing tumors (Figure S6A). Comparison of differential role of macroH2A could be more general and not low with high expressing tumors identified 3,229 differentially ex- restricted to HB, we turned to data of the Pan-Cancer cohort pressed genes (Figure 6B). A significant proportion of these of The Cancer Genome Atlas (TCGA). We have divided all sam- genes (838) were overlapping with genes differentially expressed ples with an up- or downregulated inferred transcription factor in macroH2A-depleted HepG2 compared with control cells (Fig- activity for STAT1 and the NF-kB transcription factor RELA in 4 8 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS quartiles of macroH2A1 expression (Figures S6B and S6C). As we observed a reduction of enhancer-promoter contacts that shown in Figure 6F, high STAT1 activity was preferentially found preceded changes in gene expression and enhancer activity af- in low MACROH2A1-expressing samples, while high NF-kB ac- ter TNF-a treatment. Taken together, results are compatible with tivity was associated with high MACROH2A1 levels. For the hypothesis that the influence of macroH2A on enhancers is MACROH2A2, we observed the same inverse correlation with downstream of its influence on 3D chromatin architecture. Inter- STAT1 activity and, to a lesser extent, a positive correlation estingly, identified sites of gained contacts in DKK1-MBL2 locus with NF-kB activity (Figure S6G). as well as sites that lost contact in the EREG locus were enriched In conclusion, the transcriptomic analyses of cell cultures and in macroH2A. In our future study, we will test whether this func- xenografts after macroH2A depletion and the comparison of tional association exists genome wide. high and low macroH2A-expressing HB samples suggest that Others have previously put forward the hypothesis that highly macroH2A variants have a role in modulating the inflammatory localized macroH2A-containing nucleosomes might affect the response of liver cancer cells. Taken together, our results sug- binding of transcription factors to their regulatory elements. In gest that the presence of macroH2A facilitates NF-kB responses support of this hypothesis, macroH2A1.2-containing nucleo- while inhibiting IFNg-STAT1 responses. The similar association somes interfered with NF-kB binding in vitro (Angelov et al., between transcription factor activity and macroH2A levels in 2003) and with ATF2 binding in B cells (Agelopoulos and Thanos, the large Pan-Cancer cohort of TCGA suggests that this bidirec- 2006), while in myogenic cells, macroH2A1.2 promoted the tional function of macroH2A is more general. recruitment of Pbx1 to enhancers (Dell’Orso et al., 2016). In a more general manner, singular macroH2A-containing nucleo- DISCUSSION somes positioned on promoters were suggested to facilitate the binding of some transcription factors while repelling others, Here, we provide evidence that macroH2A modulates the tran- thereby contributing to transcriptional robustness (Lavigne scriptional response of cancer cells to inflammatory signaling et al., 2015). However, in most cell types, including the here- on the level of chromatin. MacroH2A histone variants do so in used HepG2 cells, macroH2A is widely distributed in the an unprecedented manner, facilitating some responses while genome, and this distribution profile is not compatible with highly dampening others. This involves the regulation of enhancers positioned nucleosomes on promoters or enhancers (Changol- and 3D chromatin architecture of sensitive loci. We will focus kar et al., 2010; Douet et al., 2017; Gamble et al., 2010; Ratnaku- our discussion on two major points. First, we will discuss mar et al., 2012). Attempts to quantify the distribution of mac- possible molecular mechanisms by which macroH2A could roH2A have estimated that up 30% of the genome contains impact on enhancer elements. Second, we will discuss an some enrichment of macroH2A (Gamble et al., 2010), which emerging role of macroH2A in regulating autocrine and paracrine needs to be put into context with the fact that, in most cell types, signaling and its relevance for cancer. macroH2A1 levels make up 1% of the histone H2A pool (Busch- How macroH2A affects gene regulation is poorly understood. beck and Hake, 2017; Leroy et al., 2013). This means that most Here, we found that depletion of both macroH2A1.2 and mac- genomic regions scored as macroH2A bound by ChIP-seq or roH2A2 affected the transcriptional ground state of cancer cells similar epigenomic methods are likely not saturated for mac- involving both the up- and downregulation of a large number of roH2A binding. A plausible interpretation of these observations genes. When focusing on inflammatory response genes defined is that the exact positioning of macroH2A proteins is less rele- by the upregulation after treatment with TNF-a, we found that vant for their function and that their presence in few nucleo- macroH2A-regulated genes were strongly overrepresented. somes possibly interspersed with nucleosomes containing other Interestingly, the regulation of these genes required the stable H2As might suffice for macroH2A to exert its molecular function. depletion of both macroH2A1.2 and macroH2A2, since neither An important lesson comes from heterochromatin where inter- depletion nor re-expression of macroH2A1.2 or macroH2A2 mediate levels of macroH2A enrichment are sufficient to mediate alone was sufficient to elicit or rescue, respectively, the full effect 3D compaction in the nuclear space and its association with the observed in the double knockdown cells. Interestingly, we found nuclear lamina (Douet et al., 2017). Future studies will be needed that upregulated genes were generally more sensitive to pertur- to resolve the causality of macroH2A-dependent alterations in bations of individual macroH2As than downregulated genes, 3D chromatin architecture and gene regulation. possibly suggesting that repression is more directly linked to In our study, we have focused on the response of HB cells to the presence of macroH2A. The effect of perturbing macroH2A2 the acute inflammatory signal induced by the cytokine TNF-a. alone was greater than macroH2A1.2 alone. We have further Although ablation of macroH2A led to the increased expression studied the genomic loci of three macroH2A-regulated TNF-a of some TNF-a response genes such as DKK1 and PIK3AP1, response genes in greater detail. MacroH2A depletion led to globally, the effect was a dampening of the TNF-a response. changes in gene expression associated with changes in Most interestingly, we found that macroH2A had the opposite H3K27ac levels on distal enhancers. In several cases, these role for the response to IFNg mediated by JAK-STAT signaling, changes already occurred in the absence of TNF-a and were which is also observed in cancer-associated fibroblasts (Fili- associated with changes in 3D chromatin architecture including pescu et al., 2022). IFNg is a pleiotropic cytokine that can have enhancer-promoter contacts. In cases where changes in gene cytotoxic effects on tumor cells and promote the recruitment expression, enhancer activity measured by H3K27ac, and and action of immune cells to target tumor cells but, at the enhancer-promoter activity coincide, we are unable to distin- same time, can induce proliferation of cancer cells and other guish cause and consequence. However, in the case of EREG, pro-tumorigenic effects (Gocher et al., 2021). To our knowledge, Cell Reports 39, 110988, June 21, 2022 9 ll OPEN ACCESS Article such a chromatin-regulated shift in the relative response to two performed on representative gene loci. The limited number of major cytokines of the innate immunity has not been reported these examples limits our capacity to generalize our observation so far, and its physiological consequences remain elusive. that macroH2A-dependent changes in 3D chromatin architec- Furthermore, it will be interesting to assess to what extent the ture co-occur or precede changes in enhancer activity. The com- role of macroH2A extends to other paracrine stimuli from the bination of genome-wide approaches with faster perturbation microenvironment and in particular immune cells. methods will be needed to overcome these limitations in the Several lines of evidence indicate that macroH2A in cancer future. cells might also play a role in the opposite direction shaping the microenvironment through its impact on secretory pathways. STAR+METHODS In breast and prostate cancers, macroH2A inhibits genes encod- ing secreted factors promoting osteoclastogenesis (Kim et al., Detailed methods are provided in the online version of this paper 2018a, 2018b). Inversely, others have reported a positive influ- and include the following: ence. In hepatocellular carcinoma cells, for instance, macroH2A depletion promotes the secretion of cytokines and chemokines d KEY RESOURCES TABLE affecting the activation state of exposed lymphocytes (Lo Re d RESOURCE AVAILABILITY et al., 2020). This included an increased secretion of TNF-a, B Lead contact which, combined with our results, indicates that macroH2A B Materials availability can regulate cancer cells through autocrine loops. In human fi- B Data and code availability broblasts, macroH2A is required for the response to senes- d EXPERIMENTAL MODEL AND SUBJECT DETAILS cence-inducing factors and promoted the expression of genes B Cell lines encoding them, thereby propagating the signal in a paracrine B Animals manner (Chen et al., 2015). The here-identified macroH2A- d METHOD DETAILS repressed gene DKK1 encodes a paracrine factor that is a B Cloning of plasmids well-described negative regulator of the WNT signaling pathway B Cell lines, culture condition, gene transduction and (Huang et al., 2018). The literature describes a number of cancer- treatments promoting functions for DKK1 that include a contribution to im- B Xenografts mune evasion of disseminated cancer cells (Malladi et al., B Immunoblotting 2016). Interestingly, a previous study has assigned opposing B RNA isolation and gene expression analysis by RT- roles for macroH2A1 and DKK1 in affecting colorectal cancer qPCR progression through senescence (De Barrios et al., 2017). B Transcriptomic analysis Taken together, a picture emerges placing macroH2A in a cen- B Enrichment analysis of gene sets and gene ontologies tral role in intercellular communication by affecting, on one hand, B Inference of pathway and transcription factor activity the response to inflammatory signals and, on the other, the base- from transcriptomic data line and inducible expression of genes encoding autocrine and B Chromatin immunoprecipitation paracrine factors. Depending on the signal and its intracellular B UMI-4C pathway, macroH2A can both enhance or suppress the intensity B HiChIP of the inflammatory response. This has implications for under- B Use of published datasets standing the interaction of cancer cells with cells of their micro- d QUANTIFICATION AND STATISTICAL ANALYSIS environment and, in particular, immune cells. An inflammatory microenvironment has been recognized as a hallmark of cancer SUPPLEMENTAL INFORMATION (Hanahan and Weinberg, 2011). Specifically, after failed elimina- Supplemental information can be found online at https://doi.org/10.1016/j. tion of a nascent tumor by immune cells, a switch in the quality of celrep.2022.110988. the immune reaction including changes in the polarization of tu- mor-associated macrophages can promote tumor progression ACKNOWLEDGMENTS and immune evasion (Galdiero et al., 2018). The evasion of the immune system is recognized as another hallmark of cancer (Ha- We thank François Le Dily (CRG, Spain) for his help with sequencing the nahan and Weinberg, 2011). In the future, it will be exciting to UMI4C samples, members of the M.B. lab for valuable discussion of the re- sults, and the staff of IJC and IGTP core facilities for excellent support. further dissect how macroH2A-regulated paracrine signaling Some of our results have been obtained using data generated by the TCGA modulates the interaction of cancer cells with their microenviron- Research Network (https://www.cancer.gov/tcga). This research project was ment and to evaluate the relevance of this regulatory axis for im- supported by the national grants RTI2018-094005-B-I00 and BFU2015- mune therapy. 66559-P from FEDER/Ministerio de Ciencia e Innovación - Agencia Estatal de Investigación (to M.B.); PI09/00751, PI10/02082, and PI13/02340 from Limitations of the study the Instituto de Salud Carlos III (to C.A.); the MECD fellowship FPU14/06542 (to D.C.); AGAUR 2019 FI_B01024 and 2022 FI_B00528 fellowships (to In this study, we have studied transcriptional changes after J.C.-R. and A.D.R.-A., respectively); and predoctoral fellowships BES-2016- stable depletion of macroH2As in HepG2 cells upon TNF-a 077251 (to M.-M.L.P.) and PRE2019-088529 (to O.M.). Research in the M.B. treatment for 24 h. These time frames limit our capability of dis- lab is further supported by the following grants: the Marie Sk1odowska Curie tinguishing direct and indirect effects. The detailed characteriza- Training network ‘‘INTERCEPT-MDS’’ H2020-MSCA-ITN-2020-953407 (to tion of enhancer activity and chromatin conformation has been M.B.); MINECO-ISCIII PIE16/00011 (to M.B.); the Deutsche José Carreras 10 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS Leukämie Stiftung DJCLS 14R/2018 (to M.B.); AGAUR 2017-SGR-305 (to Buschbeck, M., and Hake, S.B. (2017). Variants of core histones and their roles M.B.); and Fundació La Marató de TV3 257/C/2019 (to M.B.). C.A.’s research in cell fate decisions, development and cancer. Nat. Rev. Mol. Cell Biol. 18, has received funding from the European Union’s Horizon 2020 research and 299–314. https://doi.org/10.1038/nrm.2016.166. innovation program under grant agreement nos. 668596 (ChiLTERN) and Buschbeck, M., Uribesalgo, I., Wibowo, I., Rué, P., Martin, D., Gutierrez, A., 826121 (iPC) as well as from CIBERehd (CB06/04/0033) and AGAUR (2017- Morey, L., Guigó, R., López-Schier, H., and Di Croce, L. (2009). The histone SGR-490). C.A. was supported by Ramón y Cajal (RYC-2010-07249) of the variant macroH2A is an epigenetic regulator of key developmental genes. Ministry of Science and Innovation of Spain. Research at the IJC is supported Nat. Struct. Mol. Biol. 16, 1074–1079. https://doi.org/10.1038/nsmb.1665. by the ‘‘La Caixa’’ Foundation, the Fundació Internacional Josep Carreras, Cantariño, N., Fernández-Figueras, M.T., Valero, V., Musulén, E., Malinverni, Celgene Spain, and the CERCA Programme/Generalitat de Catalunya. R., Granada, I., Goldie, S.J., Martı́n-Caballero, J., Douet, J., Forcales, S.V., et al. (2016). A cellular model reflecting the phenotypic heterogeneity of mutant AUTHOR CONTRIBUTIONS HRAS driven squamous cell carcinoma. Int. J. Cancer 139, 1106–1116. https:// doi.org/10.1002/ijc.30139. Conceptualization, D.C. and M.B.; data curation, R.M. and J.C.-R.; formal Carrillo-Reixach, J., Torrens, L., Simon-Coma, M., Royo, L., Domingo-Sàbat, analysis, R.M., D.C., A.D.R.-A., J.C.-R., and M.G.-P.; funding acquisition, M., Abril-Fornaguera, J., Akers, N., Sala, M., Ragull, S., Arnal, M., et al. (2020). C.A. and M.B.; investigation, D.C., V.V., A.P., A.G.-J., M.-M.L.P., L.R., H.R., Epigenetic footprint enables molecular risk stratification of hepatoblastoma and O.M.; methodology, B.P.-G., R.D.A., J.M.S.-P., J.L.G.-S., and L.P.; re- with clinical implications. J. Hepatol. 73, 328–341. https://doi.org/10.1016/j. sources, C.A.; supervision, J.M., J.L.G.-S., C.A., L.P., and M.B.; writing – orig- jhep.2020.03.025. inal draft, D.C. and M.B.; writing – review and editing, C.A. Chakravarthy, S., Gundimella, S.K.Y., Caron, C., Perche, P.-Y., Pehrson, J.R., Khochbin, S., and Luger, K. (2005). Structural characterization of the histone DECLARATION OF INTERESTS variant macroH2A. Mol. Cell Biol. 25, 7616–7624. https://doi.org/10.1128/ mcb.25.17.7616-7624.2005. The authors declare no competing interests. Changolkar, L.N., Singh, G., Cui, K., Berletch, J.B., Zhao, K., Disteche, C.M., and Pehrson, J.R. (2010). Genome-wide distribution of macroH2A1 histone Received: October 19, 2021 variants in mouse liver chromatin. Mol. Cell Biol. 30, 5473–5483. https://doi. Revised: April 4, 2022 org/10.1128/mcb.00518-10. Accepted: May 31, 2022 Chen, H., Ruiz, P.D., Novikov, L., Casill, A.D., Park, J.W., and Gamble, M.J. Published: June 21, 2022 (2014). MacroH2A1.1 and PARP-1 cooperate to regulate transcription by pro- moting CBP-mediated H2B acetylation. Nat. Struct. Mol. Biol. 21, 981–989. REFERENCES https://doi.org/10.1038/nsmb.2903. Chen, H., Ruiz, P.D., McKimpson, W.M., Novikov, L., Kitsis, R.N., and Gamble, Abugessaisa, I., Noguchi, S., Hasegawa, A., Kondo, A., Kawaji, H., Carninci, M.J. (2015). MacroH2A1 and ATM play opposing roles in paracrine senes- P., and Kasukawa, T. (2019). refTSS: a reference data set for human and cence and the senescence-associated secretory phenotype. Mol. Cell 59, mouse transcription start sites. J. Mol. Biol. 431, 2407–2422. https://doi.org/ 719–731. https://doi.org/10.1016/j.molcel.2015.07.011. 10.1016/j.jmb.2019.04.045. Creppe, C., Janich, P., Cantariño, N., Noguera, M., Valero, V., Musulén, E., Agelopoulos, M., and Thanos, D. (2006). Epigenetic determination of a cell- Douet, J., Posavec, M., Martı́n-Caballero, J., Sumoy, L., et al. (2012). Mac- specific gene expression program by ATF-2 and the histone variant mac- roH2A1 regulates the balance between self-renewal and differentiation roH2A. EMBO J. 25, 4843–4853. https://doi.org/10.1038/sj.emboj.7601364. commitment in embryonic and adult stem cells. Mol. Cell Biol. 32, 1442– 1452. https://doi.org/10.1128/mcb.06323-11. Alvarez, M.J., Shen, Y., Giorgi, F.M., Lachmann, A., Ding, B.B., Hilda Ye, B., and Califano, A. (2016). Functional characterization of somatic mutations in Dardenne, E., Pierredon, S., Driouch, K., Gratadou, L., Lacroix-triki, M., Espi- cancer using network-based inference of protein activity. Nat. Genet. 48, noza, M.P., Zonta, E., Germann, S., Mortada, H., Villemin, J.P., et al. (2012). 838–847. https://doi.org/10.1038/ng.3593. Splicing switch of an epigenetic regulator by RNA helicases promotes tu- mor-cell invasiveness. Nat. Struct. Mol. Biol. 19, 1139–1146. https://doi.org/ ^ té, J., Hans, F., Co, J., Angelov, D., Molla, A., Perche, P.Y., Khochbin, S., Co 10.1038/nsmb.2390. Bouvet, P., Dimitrov, S., Merci, D.D., et al. (2003). The histone variant MacroH2A interferes with transcription factor binding and SWI/SNF nucleo- Dell’Orso, S., Wang, A.H., Shih, H.-Y., Saso, K., Berghella, L., Gutierrez-Cruz, some remodeling. Mol. Cell 11, 1033–1041. https://doi.org/10.1016/s1097- G., Ladurner, A.G., O’Shea, J., O’Shea, J.J., Sartorelli, V., and Zare, H. (2016). 2765(03)00100-x. The histone variant MacroH2A1.2 is necessary for the activation of muscle en- hancers and recruitment of the transcription factor Pbx1. Cell Rep. 14, 1156– Balkwill, F. (2009). Tumour necrosis factor and cancer. Nat. Rev. Cancer 9, 1168. https://doi.org/10.1016/j.celrep.2015.12.103. 361–371. https://doi.org/10.1038/nrc2628. Douet, J., Corujo, D., Malinverni, R., Renauld, J., Sansoni, V., Posavec- Barrero, M.J., Sese, B., Martı́, M., and Izpisua Belmonte, J.C. (2013a). Macro , M., Cantariño, N., Valero, V., Mongelard, F., Bouvet, P., et al. Marjanovic histone variants are critical for the differentiation of human pluripotent cells. (2017). MacroH2A histone variants maintain nuclear organization and hetero- J. Biol. Chem. 288, 16110–16116. https://doi.org/10.1074/jbc.m113.466144. chromatin architecture. J. Cell Sci. 130, 1570–1582. Bancroft, Gregory, J, and Kelly, John, P. (1994). Macrophage Activation and Durand, N.C., Shamim, M.S., Machol, I., Rao, S.S.P., Huntley, M.H., Lander, Innate Resistance to Infection in SCID Mice. Immunobiology 191 (4), E.S., and Aiden, E.L. (2016). Juicer provides a one-click system for analyzing 424–431. https://doi.org/10.1016/S0171-2985(11)80448-1. loop-resolution Hi-C experiments. Cell Syst. 3, 95–98. https://doi.org/10. Barrero, M.J., Sese, B., Kuebler, B., Bilic, J., Boue, S., Martı́, M., Martı́, M., and 1016/j.cels.2016.07.002. Izpisua Belmonte, J.C. (2013b). Macrohistone variants preserve cell identity by ENCODE Project Consortium; Dunham, I., Kundaje, A., Aldred, S.F., Collins, preventing the gain of H3K4me2 during reprogramming to pluripotency. Cell P.J., Davis, C.A., Doyle, F., Epstein, C.B., Frietze, S., Harrow, J., et al. Rep. 3, 1005–1011. https://doi.org/10.1016/j.celrep.2013.02.029. (2013). An integrated encyclopedia of DNA elements in the human genome. De Barrios, O., Gyo} rffy, B., Fernández-Aceñero, M.J., Sánchez-Tilló, E., Sán- Nature 489, 57–74. chez-Moral, L., Siles, L., Esteve-Arenys, A., Roué, G., Casal, J.I., Darling, D.S., Fellmann, C., Hoffmann, T., Sridhar, V., Hopfgartner, B., Muhar, M., Roth, M., et al. (2017). ZEB1-induced tumourigenesis requires senescence inhibition via Lai, D.Y., Barbosa, I.A.M., Kwon, J.S., Guan, Y., et al. (2013). An optimized mi- activation of DKK1/mutant p53/Mdm2/CtBP and repression of macroH2A1. croRNA backbone for effective single-copy RNAi. Cell Rep. 5, 1704–1713. Gut 66, 666–682. https://doi.org/10.1136/gutjnl-2015-310838. https://doi.org/10.1016/j.celrep.2013.11.020. Cell Reports 39, 110988, June 21, 2022 11 ll OPEN ACCESS Article Feng, H., Zhang, Y.B., Gui, J.F., Lemon, S.M., and Yamane, D. (2021). Inter- Kapoor, A., Goldberg, M.S., Cumberland, L.K., Ratnakumar, K., Segura, M.F., feron regulatory factor 1 (IRF1) and anti-pathogen innate immune responses. Emanuel, P.O., Menendez, S., Vardabasso, C., Leroy, G., Vidal, C.I., et al. PLoS Pathog. 17, e1009220. https://doi.org/10.1371/journal.ppat.1009220. (2010). The histone variant macroH2A suppresses melanoma progression Filipescu, D., Hasson, D., Tung, N., Carcamo, S., Humblin, E., Goldberg, M., through regulation of CDK8. Nature 468, 1105–1109. https://doi.org/10. Vyas, N., Beaumont, K., Ghiraldini, F., Salmon, H., et al. (2022). MacroH2A re- 1038/nature09590. stricts melanoma progression via inhibition of in ammatory gene expression in Kim, J., Shin, Y., Lee, S., Kim, M., Punj, V., Lu, J.F., Shin, H., Kim, K., Ulmer, cancer-associated broblasts. Res. Square. https://doi.org/10.21203/rs.3.rs- T.S., Koh, J., et al. (2018a). Regulation of breast cancer-induced osteoclasto- 1578023/v1. genesis by MacroH2A1.2 involving EZH2-mediated H3K27me3. Cell Rep. 24, Fishilevich, S., Nudel, R., Rappaport, N., Hadar, R., Plaschkes, I., Iny Stein, T., 224–237. https://doi.org/10.1016/j.celrep.2018.06.020. Rosen, N., Kohn, A., Twik, M., Safran, M., et al. (2017). GeneHancer: genome- Kim, J.M., Shin, Y., Lee, S., Kim, M.Y., Punj, V., Shin, H.I., Kim, K., Koh, J.M., wide integration of enhancers and target genes in GeneCards. Database 2017, Jeong, D., and An, W. (2018b). MacroH2A1.2 inhibits prostate cancer-induced 1–17. https://doi.org/10.1093/database/bax028. osteoclastogenesis through cooperation with HP1a and H1.2. Oncogene 37, Galdiero, M.R., Marone, G., and Mantovani, A. (2018). Cancer inflammation 5749–5765. https://doi.org/10.1038/s41388-018-0356-3. and cytokines. Cold Spring Harb. Perspect. Biol. 10, a028662. https://doi. Kozlowski, M., Corujo, D., Hothorn, M., Guberovic, I., Mandemaker, I.K., org/10.1101/cshperspect.a028662. Blessing, C., Sporn, J., Gutierrez-Triana, A., Smith, R., Portmann, T., et al. Gamble, M.J., Frizzell, K.M., Yang, C., Krishnakumar, R., and Kraus, W.L. (2018). MacroH2A histone variants limit chromatin plasticity through two (2010). The histone variant macroH2A1 marks repressed autosomal chro- distinct mechanisms. EMBO Rep. 19, e44445. https://doi.org/10.15252/ matin, but protects a subset of its target genes from silencing. Genes Dev. embr.201744445. 24, 21–32. https://doi.org/10.1101/gad.1876110. Lavigne, A.-C., Castells, M., Mermet, J., Kocanova, S., Dalvai, M., and By- Garcia-Alonso, L., Holland, C.H., Ibrahim, M.M., Turei, D., and Saez-Rodri- stricky, K. (2014). Increased macroH2A1.1 expression correlates with poor guez, J. (2019). Benchmark and integration of resources for the estimation of survival of triple-negative breast cancer patients. PLoS One 9, e98930. human transcription factor activities. Genome Res. 29, 1363–1375. https:// https://doi.org/10.1371/journal.pone.0098930. doi.org/10.1101/gr.240663.118. Lavigne, M.D., Vatsellas, G., Polyzos, A., Mantouvalou, E., Sianidis, G., Mara- Gaspar-Maia, A., Qadeer, Z.A., Hasson, D., Ratnakumar, K., Adrian Leu, N., ziotis, I., Agelopoulos, M., and Thanos, D. (2015). Composite macroH2A/ Leroy, G., Liu, S., Costanzi, C., Valle-Garcia, D., Schaniel, C., et al. (2013). NRF-1 nucleosomes suppress noise and generate robustness in gene expres- MacroH2A histone variants act as a barrier upon reprogramming towards plu- sion. Cell Rep. 11, 1090–1101. https://doi.org/10.1016/j.celrep.2015.04.022. ripotency. Nat. Commun. 4, 1565. https://doi.org/10.1038/ncomms2582. Leroy, G., Dimaggio, P.A., Chan, E.Y., Zee, B.M., Blanco, M.A., Bryant, B., Fla- Gel, B., and Serra, E. (2017). KaryoploteR: An R/Bioconductor package to plot niken, I.Z., Liu, S., Kang, Y., Trojer, P., et al. (2013). A quantitative atlas of his- customizable genomes displaying arbitrary data. Bioinformatics 33 (19), 3088– tone modification signatures from human cancer cells. Epigenet. Chromatin 6, 3090. https://doi.org/10.1093/bioinformatics/btx346. 20. https://doi.org/10.1186/1756-8935-6-20. Ghiraldini, F.G., Filipescu, D., and Bernstein, E. (2021). Solid tumours hijack the López-Terrada, D., Cheung, S.W., Finegold, M.J., and Knowles, B.B. (2009). histone variant network. Nat. Rev. Cancer 21, 257–275. https://doi.org/10. Hep G2 is a hepatoblastoma-derived cell line. Hum. Pathol. 40, 1512–1515. 1038/s41568-020-00330-0. https://doi.org/10.1016/j.humpath.2009.07.003. Gocher, A.M., Workman, C.J., and Vignali, D.A.A. (2021). Interferon-g: team- Love, M.I., Huber, W., and Anders, S. (2014). Moderated estimation of fold mate or opponent in the tumour microenvironment? Nat. Rev. Immunol. 22, change and dispersion for RNA-seq data with DESeq2. Genome Biol. 15, 550. 158–172. https://doi.org/10.1038/s41577-021-00566-3. Liberzon, A., Birger, C., Thorvaldsdóttir, H., Ghandi, M., Mesirov, J.P., and Goldman, M.J., Craft, B., Hastie, M., Repec ka, K., McDade, F., Kamath, A., Tamayo, P. (2015). The Molecular Signatures Database Hallmark Gene Set Banerjee, A., Luo, Y., Rogers, D., Brooks, A.N., et al. (2020). Visualizing and in- Collection. Cell Syst. 1 (6), 417–425. https://doi.org/10.1016/j.cels.2015.12. terpreting cancer genomics data via the Xena platform. Nat. Biotechnol. 38, 004. 675–678. https://doi.org/10.1038/s41587-020-0546-8. Lo Re, O., Fusilli, C., Rappa, F., Van Haele, M., Douet, J., Pindjakova, J., Ro- Hanahan, D., and Weinberg, R.A. (2011). Hallmarks of cancer: the next gener- cha, S.W., Pata, I., Valcı́ková, B., Uldrijan, S., et al. (2017). Induction of cancer ation. Cell 144, 646–674. https://doi.org/10.1016/j.cell.2011.02.013. cell stemness by depletion of macrohistone H2A1 in hepatocellular carcinoma. Hepatology 67, 636–650. https://doi.org/10.1002/hep.29519. Hernández-Muñoz, I., Lund, A.H., van der Stoop, P., Boutsma, E., Muijrers, I., Verhoeven, E., Nusinow, D.A., Panning, B., Marahrens, Y., and van Lohuizen, Lo Re, O., Mazza, T., Giallongo, S., Sanna, P., Rappa, F., Vinh Luong, T., Li M. (2005). Stable X chromosome inactivation involves the PRC1 Polycomb Volti, G., Drovakova, A., Roskams, T., van Haele, M., et al. (2020). Loss of his- complex and requires histone MACROH2A1 and the CULLIN3/SPOP ubiquitin tone macroH2A1 in hepatocellular carcinoma cells promotes paracrine-medi- E3 ligase. Proc. Natl. Acad. Sci. U S A 102, 7635–7640. https://doi.org/10. ated chemoresistance and CD4+CD25+FoxP3+ regulatory T cells activation. 1073/pnas.0408918102. Theranostics 10, 910–924. https://doi.org/10.7150/thno.35045. Hsu, C.J., Meers, O., Buschbeck, M., and Heidel, F.H. (2021). The role of mac- Malladi, S., MacAlinao, D.G., Jin, X., He, L., Basnet, H., Zou, Y., De Stanchina, roh2a histone variants in cancer. Cancers 13, 3003. https://doi.org/10.3390/ E., and Massagué, J. (2016). Metastatic latency and immune evasion through cancers13123003. autocrine inhibition of WNT. Cell 165, 45–60. https://doi.org/10.1016/j.cell. 2016.02.025. Huang, Y., Liu, L., and Liu, A. (2018). Dickkopf-1: current knowledge and related diseases. Life Sci. 209, 249–254. https://doi.org/10.1016/j.lfs.2018. Mallona, I., Dı́ez-Villanueva, A., Martı́n, B., and Peinado, M.A. (2017). Chainy: 08.019. an universal tool for standardized relative quantification in real-time PCR. Bio- Hurtado-Bagès, S., Posavec Marjanovic, M., Valero, V., Malinverni, R., Corujo, informatics 33, 1411–1413. https://doi.org/10.1093/bioinformatics/btw839. D., Bouvet, P., Lavigne, A.C., Bystricky, K., and Buschbeck, M. (2020). The his- Martire, S., and Banaszynski, L.A. (2020). The roles of histone variants in fine- tone variant MacroH2A1 regulates key genes for myogenic cell fusion in a tuning chromatin organization and function. Nat. Rev. Mol. Cell Biol. 21, splice-isoform dependent manner. Cells 9, 1109. https://doi.org/10.3390/ 522–541. https://doi.org/10.1038/s41580-020-0262-8. cells9051109. Morgenstern, J.P., and Land, H. (1990). Advanced mammalian gene transfer: Jiang, P., Freedman, M.L., Liu, J.S., and Liu, X.S. (2015). Inference of transcrip- high titre retroviral vectors with multiple drug selection markers and a comple- tional regulation in cancers. Proc. Natl. Acad. Sci. U S A 112, 7731–7736. mentary helper-free packaging cell line. Nucleic Acids Res. 18, 3587–3596. https://doi.org/10.1073/pnas.1424272112. https://doi.org/10.1093/nar/18.12.3587. 12 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS Mumbach, M.R., Rubin, A.J., Flynn, R.A., Dai, C., Khavari, P.A., Greenleaf, genes reveal signaling footprints in cancer gene expression. Nat. Commun. 9, W.J., and Chang, H.Y. (2016). HiChIP: efficient and sensitive analysis of pro- 20. https://doi.org/10.1038/s41467-017-02391-6. tein-directed genome architecture. Nat. Methods 13, 919–922. https://doi. Schwartzman, O., Mukamel, Z., Oded-Elkayam, N., Olivares-Chauvet, P., Lu- org/10.1038/nmeth.3999. bling, Y., Landan, G., Izraeli, S., and Tanay, A. (2016). UMI-4C for quantitative Ouararhni, K., Hadj-Slimane, R., Ait-Si-Ali, S., Robin, P., Mietton, F., Harel-Bel- and targeted chromosomal contact profiling. Nat. Methods 13, 685–691. lan, A., Dimitrov, S., and Hamiche, A. (2006). The histone variant mH2A1.1 in- https://doi.org/10.1038/nmeth.3922. terferes with transcription by down-regulating PARP-1 enzymatic activity. Soneson, C., Love, M.I., and Robinson, M.D. (2016). Differential analyses for Genes Dev. 20, 3324–3336. https://doi.org/10.1101/gad.396106. RNA-seq: transcript-level estimates improve gene-level inferences. Pasque, V., Gillich, A., Garrett, N., and Gurdon, J.B. (2011). Histone variant F1000Research 4, 1521. https://doi.org/10.12688/f1000research.7563.2. macroH2A confers resistance to nuclear reprogramming. EMBO J. 30, Stovner, E.B., and Sætrom, P. (2019). Epic2 efficiently finds diffuse domains in 2373–2387. https://doi.org/10.1038/emboj.2011.144. ChIP-seq data. Bioinformatics 35, 4392–4393. https://doi.org/10.1093/bioin- Patro, R., Duggal, G., Love, M.I., Irizarry, R.A., and Kingsford, C. (2017). formatics/btz232. Salmon provides fast and bias-aware quantification of transcript expression. Subramanian, A., Tamayo, P., Mootha, V.K., Mukherjee, S., and Ebert, B.L. Nat. Methods 14, 417–419. https://doi.org/10.1038/nmeth.4197. (2005). Gene set enrichment analysis : a knowledge-based approach for inter- Perkins, N.D. (2012). The diverse and complex roles of NF-kB subunits in can- preting genome-wide. Proc. Natl. Acad. Sci. U S A 102, 15545–15550. cer. Nat. Rev. Cancer 12, 121–132. https://doi.org/10.1038/nrc3204. Sun, Z., Filipescu, D., Andrade, J., Gaspar-Maia, A., Ueberheide, B., and Bern- Perkins, J.R., Dawes, J.M., McMahon, S.B., Bennett, D.L.H., Orengo, C., and stein, E. (2018). Transcription-associated histone pruning demarcates mac- Kohl, M. (2012). ReadqPCR and NormqPCR: R packages for the reading, qual- roH2A chromatin domains. Nat. Struct. Mol. Biol. 25, 958–970. https://doi. ity checking and normalisation of RT-qPCR quantification cycle (Cq) data. org/10.1038/s41594-018-0134-5. BMC Genom. 13, 296. https://doi.org/10.1186/1471-2164-13-296. Taniguchi, K., and Karin, M. (2018). NF-kB, inflammation, immunity and can- Posavec, M., Timinszky, G., and Buschbeck, M. (2013). Macro domains as cer: coming of age. Nat. Rev. Immunol. 18, 309–324. https://doi.org/10. metabolite sensors on chromatin. Cell. Mol. Life Sci. 70, 1509–1524. https:// 1038/nri.2017.142. doi.org/10.1007/s00018-013-1294-4. Timinszky, G., Till, S., Hassa, P.O., Hothorn, M., Kustatscher, G., Nijmeijer, B., Ramos-Rodrı́guez, M., Subirana-Granés, M., and Pasquali, L. (2021). UMI4- Colombelli, J., Altmeyer, M., Stelzer, E.H.K., Scheffzek, K., et al. (2009). A mac- Cats: an R package to analyze chromatin contact profiles obtained by UMI- rodomain-containing histone rearranges chromatin upon sensing PARP1 acti- 4C. Bioinformatics, 1–3. https://doi.org/10.1093/bioinformatics/btab392. vation. Nat. Struct. Mol. Biol. 16, 923–929. https://doi.org/10.1038/nsmb. Ratnakumar, K., Duarte, L.F., LeRoy, G., Hasson, D., Smeets, D., Vardabasso, 1664. C., Bönisch, C., Zeng, T., Xiang, B., Zhang, D.Y., et al. (2012). ATRX-mediated Weinstein, J.N., Collisson, E.A., Mills, G.B., Shaw, K.R.M., Ozenberger, B.A., chromatin association of histone variant macroH2A1 regulates a-globin Ellrott, K., Shmulevich, I., Stuart, J.M., Chang, K., Creighton, C.J., et al. expression. Genes Dev. 26, 433–438. https://doi.org/10.1101/gad.179416. (2013). The cancer genome atlas pan-cancer analysis project. Nat. Genet. 111. 45, 1113–1120. https://doi.org/10.1038/ng.2764. Ritchie, M.E., Phipson, B., Wu, D., Hu, Y., Law, C.W., Shi, W., and Smyth, G.K. Wickham, H. (2016). ggplot2: Elegant Graphics for Data Analysis (Springer- (2015). Limma powers differential expression analyses for RNA-sequencing Verlag New York). and microarray studies. Nucleic Acids Res. 43, e47. https://doi.org/10.1093/ Wingett, S., Ewels, P., Furlan-Magaril, M., Nagano, T., Schoenfelder, S., nar/gkv007. Fraser, P., and Andrews, S. (2015). HiCUP: pipeline for mapping and process- Ritz, C., and Spiess, A.N. (2008). qpcR: an R package for sigmoidal model se- ing Hi-C data. F1000Research, 1310. https://doi.org/10.12688/f1000research. lection in quantitative real-time polymerase chain reaction analysis. Bioinfor- 7334.1. matics 24, 1549–1551. https://doi.org/10.1093/bioinformatics/btn227. Wu, T., Hu, E., Xu, S., Chen, M., Guo, P., Dai, Z., Feng, T., Zhou, L., Tang, W., Roadmap Epigenomics Consortium; Kundaje, A., Meuleman, W., Ernst, J., Bi- Zhan, L., et al. (2021). clusterProfiler 4.0: a universal enrichment tool for inter- lenky, M., Yen, A., Heravi-Moussavi, A., Kheradpour, P., Zhang, Z., Ziller, M.J., preting omics data. Innovation 2, 100141. https://doi.org/10.1016/j.xinn.2021. Wang, J., et al. (2015). Integrative analysis of 111 reference human epige- 100141. nomes. Nature 518, 317–330. https://doi.org/10.1038/nature14248. Yardımcı, G.G., Ozadam, H., Sauria, M., Ursu, O., Yan, K.-K., Yang, T., Chak- Schnater, J.M., Köhler, S.E., Lamers, W.H., Von Schweinitz, D., and Aronson, raborty, A., Kaul, A., Lajoie, B., Song, F., et al. (2017). Measuring the reproduc- D.C. (2003). Where do we stand with hepatoblastoma? A review. Cancer 98, ibility and quality of Hi-C data. Genome Biol. 20, 57. 668–678. https://doi.org/10.1002/cncr.11585. Zelová, H., and Hosek, J. (2013). TNF-a signalling and inflammation: interac- €nemann, M., Sieber, A., Uhlitz, F., Sauer, S., Gar- Schubert, M., Klinger, B., Klu tions between old acquaintances. Inflamm. Res. 62, 641–651. https://doi. €thgen, N., and Saez-Rodriguez, J. (2018). Perturbation-response nett, M.J., Blu org/10.1007/s00011-013-0633-0. Cell Reports 39, 110988, June 21, 2022 13 ll OPEN ACCESS Article STAR+METHODS KEY RESOURCES TABLE REAGENT or RESOURCE SOURCE IDENTIFIER Antibodies anti-macroH2A1 (Buschbeck et al., 2009; Douet et al., 2017) N/A anti-macroH2A2 (Douet et al., 2017) N/A anti-DKK1 Santa Cruz Biotechnology Cat#sc-374574; RRID: AB_10989416 anti-FLAG Sigma-Aldrich Cat#F1804; RRID:AB_262044 anti-H3 Abcam Cat#ab1791; RRID:AB_302613 anti-H3K4me3 Abcam Cat#ab8580; RRID:AB_306649 anti-H3K27Ac Abcam Cat#ab4729; RRID:AB_2118291 anti-IgG Abcam Cat#ab46540; RRID:AB_2614925 IRDye 680RD Goat anti-Rabbit IgG LI-COR Biosciences Cat# 926-68071, RRID:AB_10956166 Secondary Antibody IRDye 680RD Goat anti-Mouse IgG LI-COR Biosciences Cat# 926-68070, RRID:AB_10956588 IRDye 800CW Goat anti-Rabbit IgG LI-COR Biosciences Cat# 926-32211, RRID:AB_621843 IRDye 800CW Goat anti-Mouse IgG LI-COR Biosciences Cat# 926-32210, RRID:AB_621842 Goat Anti-Rabbit Immunoglobulins/HRP Dako P0448 Goat Anti-Mouse Immunoglobulins/HRP Dako P0447 Chemicals, peptides, and recombinant proteins Dulbecco’s modified Eagle medium Thermo Fisher Scientific Cat# 11960085 (DMEM) Fetal Bovine Serum Thermo Fisher Scientific Cat# 10270106 L-Glutamine (200mM) Thermo Fisher Scientific Cat# 25030024 Penicillin-Streptomycin (10,000 U/mL) Thermo Fisher Scientific Cat# 15140122 Recombinant Human TNF-a Peprotech Cat# 300-01A Recombinant Human IFN-g Peprotech Cat# 300-02 BMS-345541 Abcam Cat# ab144822 Matrigel Matrix Corning Cat# 354234 DpnII New England Biolabs Cat# 174R0543M T4 ligase Promega Cat# M1804 Critical commercial assays Maxwell RSC simplyRNA Cells Kit Promega Cat# AS1390 Pierce ECL Western Blotting Substrate Thermo Fisher Scientific Cat# 32106 TRIzol Reagent Thermo Fisher Scientific Cat# 15596026 First Strand cDNA Synthesis Kit Thermo Fisher Scientific Cat# K1612 LightCycler 480 SYBR Green I Master Roche Cat# 04707516001 Magna ChIPTM Protein A+G Magnetic Merck Cat# 16-663 Beads ChIP DNA Clean & Concentrator Zymo Research Cat# D5205 QubitTM dsDNA HS Thermo Fisher Scientific Cat# Q32851 Deposited data RNA-Seq, HiChIP and UMI-4C This paper GSE161859 RNA-Seq from 32 patients with (Carrillo-Reixach et al., 2020) GSE133039 hepatoblastoma ChIP-Seq of macroH2A1 and macroH2A2 (Douet et al., 2017) GSE58175 in HepG2 ChIP-Seq of H3K27Ac in HepG2 ENCODE (ENCODE Project Consortium GEO:GSM733743 et al., 2013) (Continued on next page) e1 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS Continued REAGENT or RESOURCE SOURCE IDENTIFIER ChIP-Seq of RELA in HepG2 ENCODE (ENCODE Project Consortium GSE170794 et al., 2013) HiC in HepG2 ENCODE (ENCODE Project Consortium ENCSR194SRI et al., 2013) GeneHancer Database (Fishilevich et al., 2017) N/A refTSS (Abugessaisa et al., 2019) N/A HMM Chromatin Segmentation (Core (Roadmap Epigenomics Consortium et al., N/A 15-state model) 2015) TCGA Pan-Cancer Batch effects (Weinstein et al., 2013), accessed through EB++AdjustPANCAN_IlluminaHiSeq_ normalized mRNA data the UCSC Xena Browser (Goldman et al., RNASeqV2.geneExp.xena 2020) TCGA Pan-Cancer RABIT transcription (Jiang et al., 2015), accessed through the RABIT/pancan/RABIT_pancan.HiSeq.V2 factor regulatory impact, HiSeqV2 UCSC Xena Browser (Goldman et al., 2020) Molecular Signatures Database (Liberzon et al., 2015) MSigDB v7.2 Experimental models: Cell lines HepG2 ATCC Cat#HB-8065; RRID:CVCL_0027 Experimental models: Organisms/strains NOD.CB17-PRKDCSCID Mouse Charles River Laboratories Strain code:394 Oligonucleotides For a list of all primers used see Table S2 This paper N/A Recombinant DNA pBabe.puro (Morgenstern and Land, 1990) Addgene# 1764 pBabe.puro – FLAG-mH2A1.2 This paper N/A pBabe.puro – FLAG-mH2A2 This paper N/A SGEP (Fellmann et al., 2013) Addgene# 111170 SGEP - sh Renilla This paper N/A SGEP - sh mH2A1 This paper N/A SGEP - sh mH2A2 This paper N/A Software and algorithms Salmon (Patro et al., 2017) RRID:SCR_017036 DESeq2 (Love et al., 2014) RRID:SCR_015687 tximport (Soneson et al., 2016) RRID:SCR_016752 limma (Ritchie et al., 2015) RRID:SCR_010943 GSEA v4.03 for Mac (Subramanian et al., 2005) RRID:SCR_003199 clusterProfiler (Wu et al., 2021) RRID:SCR_016884 PROGENy (Schubert et al., 2018) N/A VIPER (Alvarez et al., 2016) N/A DoRothEA (Garcia-Alonso et al., 2019) RRID:SCR_005569 Juicer (Durand et al., 2016) RRID:SCR_017226 RESOURCE AVAILABILITY Lead contact Further information and requests for resources and reagents should be directed to and will be fulfilled by the lead contact, Marcus Buschbeck (mbuschbeck@carrerasresearch.org). Materials availability All unique/stable reagents generated in this study are available from the lead contact with a completed Materials Transfer Agreement (MTA). Cell Reports 39, 110988, June 21, 2022 e2 ll OPEN ACCESS Article Data and code availability d The RNA-Seq, HiChIP and UMI-4C data generated and analyzed for this study has been deposited in the Gene Expression Omnibus and are publicly available as of the date of publication. Accession numbers are listed in the key resources table. d This paper analyzes existing, publicly available data. The accession numbers for the datasets are listed in the key resources table. d This paper does not report original code. d Any additional information required to reanalyze the data reported in this paper is available from the lead contact upon request. EXPERIMENTAL MODEL AND SUBJECT DETAILS Cell lines HepG2 cells (ATCC, HB-8065) were routinely cultured in Dulbecco’s modified Eagle medium (DMEM) containing 4.5 g/L glucose (Gibco) supplemented with 10% v/v fetal bovine serum (FBS) (Gibco), 2 mM L-glutamine (Gibco), 50 U/mL penicillin (Gibco) and 50 mg/mL streptomycin (Gibco) at 37 C in 5% CO2. Animals Four weeks old NOD.CB17-PRKDCSCID immunosuppressed mice were purchased from Charles River Laboratories (strain code 394). The animals were kept in specific pathogen free (SPF) conditions for the duration of the experiment. The animals were kept in the stall for two weeks before the experimental procedure. Three male and two female animals were used for each experimental group. The Animal Experimentation Commission of the Catalan Government (Comissó d’Experimentació Animal de la Generalitat de Catalunya) has approved the animal experimentation procedures performed in this study with project reference 9500. METHOD DETAILS Cloning of plasmids All cloning procedures were performed following standard cloning techniques involving PCR amplification, restriction enzyme diges- tion and ligation reactions. Full-length FLAG-tagged mouse macroH2A1.2 and macroH2A2 were generated by cloning the PCR amplified cDNA of mouse macroH2A2 into a pCMV-Tag2A and then subcloning into a pBabe.puro backbone. Short-hairpin RNA se- quences were cloned into pRRL-SFFV-GFP- miRE-PGK-Puro (SGEP) vectors as described (Fellmann et al., 2013). Cell lines, culture condition, gene transduction and treatments HepG2 cell lines with a stable shRNA-mediated double knock-down of macroH2A1 and macroH2A2 (DKD) and their random sequence control counterparts (CTL) have been described previously (Douet et al., 2017). HepG2 cells (ATCC, HB-8065) were routinely cultured in Dulbecco’s modified Eagle medium (DMEM) containing 4.5 g/L glucose (Gibco) supplemented with 10% v/v fetal bovine serum (FBS) (Gibco), 2 mM L-glutamine (Gibco), 50 U/mL penicillin (Gibco) and 50 mg/mL streptomycin (Gibco) at 37 C in 5% CO2. Stable cell lines were generated by transduction of viral vectors selected as previously described (Cantariño et al., 2016). pRRL- SFFV-GFP- miRE-PGK-Puro (SGEP, Fellmann et al., 2013)) vectors encoding either a shRNA cassette against transcripts of Renilla (non-targeting control), mH2A1 or mH2A2 were used for transduction to generate stable single knock-downs in HepG2. The following procedure was followed to treat the cells with TNFa or IFNg: 500.000 HepG2 cells were seeded in six-well plates and incubated overnight in standard conditions. Cells were starved of FBS over night by reducing its concentration to 0.5%. After that, TNFa or IFNg (PeproTech, resuspended in PBS 0.1% BSA) was added to the media to a final concentration of 20 ng/mL or 100ng/mL, respectively (or at different concentrations if indicated in the figure). An equal volume of PBS 0.1% BSA was added as the untreated control condition. The cells were treated for 24h with TNFa and 48h with IFNg, then collected and subject to protein and RNA expres- sion analysis. The same procedure was followed when evaluating the effect of the IKK inhibitor BMS-345541 (Abcam, resuspended in DMSO) with the following modifications: 150.000 cells were seeded in 12-well plates, the inhibitor was added at a final concentration of 2 mM 1h prior to adding TNFa and an equal volume of DMSO was used as the untreated control condition. To obtain samples for ChIP experiments, the same treatment protocol was followed but 10 million cells were seeded in P15 plates and grown for two days before applying the o/n serum starvation, subsequent TNFa treatment and fixation after 24h. Xenografts The Animal Experimentation Commission of the Catalan Government (Comissó d’Experimentació Animal de la Generalitat de Cata- lunya) has approved animal experimentation procedure used for the xenograft experiment with project reference 9500. Four weeks old NOD.CB17-PRKDCSCID immunosuppressed mice were purchased from Charles Rivers Laboratories (strain code 394). The an- imals were kept in specific pathogen free (SPF) conditions for the duration of the experiment. The animals were kept in the stall for two weeks before the experimental procedure. CTL and DKD HepG2 cells were resuspended in a Matrigel (Corning) and PBS mix 1:1. The animals were anesthesiated with isofluorane and two subcutaneous injections were performed in each lower flank of the animal consisting of 100 mL Matrigel:PBS 1:1 containing 5 million cells. A total of 5 mice (three males, two females) per condition was e3 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS used, effectively resulting in 10 injections for each cell line. The weight and overall condition of the mice were regularly monitored. The tumor size was monitored every 2-3 days and measured with a digital caliper. Two measurements were performed for each tumor approximating the longer (a) and shorter (b) axis of an ellipse fitted to the tumor and the tumor volume was estimated by calculating ða3bÞ2 =2. 21 days after the injection, tumors were resected with the mice under anesthesia and all animals were sacrificed. Tumors were cleaned with PBS and cut in half. One half was fixed in formalin and embedded in paraffin for histological analysis, while the other half was immediately frozen in liquid nitrogen for later RNA extraction. Immunoblotting Protein extraction was performed by resuspending cell pellets Laemmli’s buffer. After heating at 95 C for 10 min, the samples were analyzed by Western blot. Samples were loaded on 12% polyacrylamide gels and run at 36 mA for 45-60 min to achieve optimal separation and then transferred at 220 mA for 90 min onto a nitrocellulose membrane (GE healthcare). After transfer, protein load was checked by Ponceau staining and the membranes were blocked using 5% non-fat milk (Nestlé) in TBST for 15 min. Mem- branes were incubated with primary antibodies o/n at 4 C. The day after, membranes were washed thrice with TBST for 10 min and incubated for 1 h at RT with either Fluorophore-conjugated secondary antibodies (Figure 5B, LI-COR Biosciences IRDye 680RD and IRDye 800CW) or HRP-conjugated antibodies (Figure 1A, DakoCytomation). Membranes were then washed again thrice with TBST for 10min. For the detection of fluorescent signals, the dried membranes were scanned with an Odyssey CLx Imager. For the detection of chemiluminescent signals, the membranes were incubated 1 min with PierceTM ECL Western Blotting Substrate (Thermo Scientific) chemiluminescent reagent mix in a 1:1 ratio. Dried membranes were overlaid with photo- graphic film (GE Healthcare) for a few seconds up to 5 min to achieve proper exposures and the films were developed with a FujiFilm FPM-100A developer. For the detection of secreted proteins, 1 mL of culture media was collected and used for acetone precipitation of extracellular pro- teins. First, the sample was briefly centrifuged to remove any detached cells that might be in suspension. Then, four volumes of cold (20 C) acetone were added to the tube. After vortexing, the sample was incubated o/n at 20 C. The sample was centrifuged 10min at 13000g and the supernatant was discarded. The protein pellet was air-dried at room temperature for 30 min and. RNA isolation and gene expression analysis by RT-qPCR Cultured cells were collected by trypsinization and washed twice with PBS. Total RNA was isolated using the Maxwell RSC sim- plyRNA Cells Kit (Promega) together with the Maxwell RSC Instrument (Promega) according to the manufacturer’s instructions. To- tal RNA from xenograft samples was isolated using the TRIzolTM Reagent (ThermoFisher Scientific). Concentration and quality of the extracted RNA was checked using a Nanodrop instrument (Thermo Scientific). For quantitative PCR after reverse transcription (RT-qPCR), 1 mg of total RNA was used for cDNA synthesis using the First Strand cDNA synthesis kit (Thermo Scientific) according to the manufacturer’s instructions. Resulting cDNAs were diluted 10 fold in auto- claved ddH2O previous to real-time PCR analysis using a LightCycler 480 II instrument (Roche). A 10 mL reaction mix containing 5 mL LightCycler 480 SYBR Green I Master (Roche), 0.5 mL of forward and reverse primers (for a final concentration of 1 mM), 2 mL of diluted cDNA and 2 mL of autoclaved ddH2O was loaded in 384-well plates. The median efficiency of each primer pair was determined in every RT-qPCR plate run using the software chainy (Mallona et al., 2017) that is based on the R packages qpcR and NormqPCR (Perkins et al., 2012; Ritz and Spiess, 2008). The expression of a target RNA in a sample relative to a control sample was calculated using a variation of the DCt method that accounts for the efficiency of the ðCtcontrol  Ctsample Þ primers as follows: efficiencytarget . All values were normalized by the average expression of two reference genes: GAPDH and RPLP0. The normalized relative expression value is used in all the plots and statistics. All quantifications of RNA expression are relative, so that the chosen control sample will have an expression of 1 for all target genes and the rest of the samples will have an expression that represents a multiple or fraction of the expression in the control. Transcriptomic analysis One microgram of total RNA was submitted to the Genomics Unit of the Centre for Genomic Regulation (Barcelona, Spain) who prepared the libraries with a poly(A) selection and sequenced them in an Illumina HiSeq2500 using 50bp single reads obtaining >30M reads per sample. Five different samples of each X.CTL and X.DKD and three replicates of CTL and DKD untreated or treated with TNFa were sequenced. To quantify the transcript expression we used Salmon (Patro et al., 2017) version 0.14.0, using libtype A and the validateMappigs option active. The raw reads were ‘‘quasi-aligned’’ to the human Ensembl genome annotation GRCh38.12. The quantification results obtained from Salmon were imported to the DESeq2 R package (Love et al., 2014) using the tximport R package (Soneson et al., 2016). For the xenograft samples, the design formula for DE- Seq2 included only the type of tumor as a variable (X.CTL or X.DKD) and the contrast X.DKD vs X.CTL was applied using the results function to obtain log2FC and p-values. For the cultured HepG2 cells samples, the experimental groups were defined as a combination of cell line and treatment, resulting in four groups: CTL untreated, DKD untreated, CTL + TNFa and DKD + TNFa. The design formula for DESeq2 included the experimental group plus the experimental replicate. The variance-stabilized log2 gene expression values were extracted using the rlog or vst functions (for HepG2 and hepatoblastoma datasets, respectively), which minimizes differences between samples for rows with small counts. For the analysis of normal liver tissue and hepato- blastoma tumor expression, the variance-stabilized log2 expression values of MACROH2A1 were used to classify the tumor Cell Reports 39, 110988, June 21, 2022 e4 ll OPEN ACCESS Article samples in 3 quantiles. The quantile annotation was used as design formula for DESeq2 and the contrast between the lowest and highest mH2A1 expression quantiles was performed. Specific contrasts where extracted as indicated in each figure to obtain log2FC and p-values using the lfcShrink function with the ashr type of log2FC shrinkage. The removeBatchEffect function from the limma R package (Ritchie et al., 2015) was used on the log2 expression matrix to reduce the variability explained by different experimental replicates for the purposes of Principal Component Analysis, sample distances calculation and log2 gene expression plotting. Ensembl Gene IDs where mapped to HGNC symbols whenever possible, otherwise the Ensembl Gene ID was kept as annotation. In all cases, genes with a Benjamini-Hochberg adjusted p-value lower than 0.05 were selected as being significantly differentially expressed. Enrichment analysis of gene sets and gene ontologies Differentially expressed gene lists with an HGNC annotation were ranked by log2FC value for a Gene Set Enrichment Analysis (GSEA) using the pre-ranked functionality of the GSEA software (Subramanian et al., 2005) (v4.03 for Mac). As gene database we selected the hallmark gene sets from the Molecular Signatures Database (MSigDB v7.2 updated September 2020) (Liberzon et al., 2015) which contains gene expression signatures derived by aggregating many MSigDB gene sets to represent well-defined biological states or processes. For gene ontology enrichment analysis, the function enrichGO from the R package clusterProfiler (Wu et al., 2021) was used, with minGSS = 20 and an adjusted p-value cut-off of 0.05. Enriched categories were ranked according to the number of differentially ex- pressed genes and a subset of the top categories listed was plotted in the figures. Inference of pathway and transcription factor activity from transcriptomic data We used the R package PROGENy (Schubert et al., 2018) on the rlog variance-stabilized gene expression matrices to calculate pathway activity scores per sample according to the top 100 genes of each pathway. After scaling the pathway scores to the control samples, we applied a linear model to identify the pathways significantly different between experimental groups (pv of the estimate coefficient < 0.05). We used the R packages VIPER (Alvarez et al., 2016) and DoRothEA (Garcia-Alonso et al., 2019) to evaluate reg- ulon activity from gene expression data. For Figures 6 and S2, the ranked log2FC of differentially expressed genes was used as input for msVIPER with regulons of the highest confidence (A). For Figure 5, The normalized gene expression counts were extracted using DESeq2 function counts and used as input for VIPER as implemented in the wrapper function run_viper in DoRothEA with the param- eters method = ‘‘scale’’, minsize = 4 and eset.filter = FALSE and the two regulons of most confidence (A, B). A two-tailed Student’s t-test was performed and TFs with an adjusted p-value (Bonferroni-Hochberg) less than 0.05 were considered to have significantly different inferred activity between the experimental groups. Chromatin immunoprecipitation Between 20 and 30 million HepG2 cells were fixed with a PBS solution containing 1% FA for 10 min at RT. The fixation was quenched by the addition of glycine to a concentration of 0.125M and incubation for 5 min at RT. Cells were then washed with cold PBS and scraped off the plate. Nuclei extraction was achieved by resuspending the fixed cells in 500 mL lysis buffer I (5 mM PIPES pH 8, 85 mM KCl, 0.5% NP-40, 1 mM PMSF, 50 mg/mL leupeptin) and incubating for 30min at 4 C in rotation. Lysed cells were centrifuged at 5000 rpm for 5 min at 4 C, the supernatant was discarded and the pellet resuspended in 500 mL lysis buffer II (1% SDS, 10 mM EDTA pH 8, 50 mM TRIS pH 8, 1 mM PMSF, 50 mg/mL leupeptin) and incubated for 30min at 4 C in rotation. Chromatin shearing was achieved by sonication in a Bioruptor Pico (Diagenode, 10 cycles 30sec ON/30 sec OFF at high intensity). Samples were then centrifuged for 10 min at 10000rpm to remove insoluble debris and a 10 mL aliquot of lysate was purified, quantified and run on an agarose gel to test the size range after sonication. 35 mg of sheared chromatin were used for each immunoprecipitation (IP) reaction diluted with 9 volumes of IP dilution buffer (1% Triton X-100, 150 mM NaCl, 2 mM EDTA pH 8, 20 mM TRIS pH 8, 1mM PMSF). Samples were precleared with 20 mL of Magna ChIPTM protein A/G magnetic beads (Merck Millipore) for 2 h at 4 C in rotation. 10% of the precleared lysate was taken as input. Antibodies were added to each IP and samples were incubated o/n at 4 C in rotation. The following day 20 mL of Magna ChIPTM protein A/G magnetic beads (Merck Millipore) were added to capture the antibody-bound chromatin fraction for 2 h at 4 C in rotation. Beads were washed twice at RT with each of the following buffers in this order: mixed micelle was buffer (140 mM NaCl, 20 mM Tris-HCl pH 8, 5 mM EDTA pH 8, 5% sucrose w/v, 1% Triton X-100 v/v, 0.2% SDS, 0.02% NaN3), LiCl/detergent wash (0.5% deoxycholic acid so- dium salt w/v, 1 mM EDTA, 250 mM LiCl, 0.5% NP-40 v/v, 10 mM Tris-HCl pH 8), buffer 500 (0.1% deoxycholic acid sodium salt w/v, 1mM EDTA, 50mM HEPES pH 7.5, 500 mM NaCl, 1% Triton X-100 v/v, 0.02% NaN3) and TE buffer (10mM Tris-HCl pH8, 1 mM EDTA pH 8). The enriched chromatin was eluted by resuspending the beads in 200 mL of elution buffer (1%SDS, 100 mM NaHCO3), 30 sec of vortexing and incubation for 30 min at RT in rotation. Two sequential elutions were performed and the el- uates were combined. The input samples were included again in the experimental procedure at this point. Crosslinking reversal was achieved by adding NaCl to a final concentration of 200 mM and incubating the samples o/n at 65 C shaking. The following day the protein in the samples was digested by adding 2 mL of Proteinase K 10 mg/mL, 16 mL 1M TRIS pH 6.5 and 8 mL 500 mM EDTA and incubating for 2 h at 45 C. ChIP DNA was purified using ChIP DNA Clean & Concentrator (Zymo Research) and analyzed by qPCR in a similar way to the ‘‘Gene expression analysis by RT-qPCR’’ description. The input of each sample was used as the ‘‘control’’ sample to obtain a relative enrichment over input, using the same formula than for RNA expression. e5 Cell Reports 39, 110988, June 21, 2022 ll Article OPEN ACCESS The relative enrichment was then converted to a percentage. The percentage of input value divided by the signal of the H3 ChIP is used in the ChIP plots. UMI-4C UMI-4C was performed essentially following a described protocol (Schwartzman et al., 2016). More specifically, 2.5 million HepG2 cells were seeded in 6-well plates. The day after, the cells were fixed by addition of PBS 1% FA and incubation for 10 min at RT. The reaction was stopped by adding glycine to a final concentration of 0.125M and incubation for 5 min on ice. The cells were then collected by scraping, washed with PBS and resuspended in 500 mL of lysis buffer (10 mM TRIS pH 8, 10 mM NaCl, 0.2% NP-40) and incubated for 30 min on ice. Nuclei were isolated by centrifugation for 5 min at 760 g and resuspended in 500 mL of DpnII buffer (NEB). 15 mL of SDS 10% were added to the tube and incubated for 60min at 37 C shaking at 900rpm. Then, 150 mL of 10% Triton X-100 were added and the sample was incubated for 60 min at 37 C. Afterward, digestion of the nuclei was achieved by addition of 400U DpnII (NEB) and o/n incubation at 37 C. The next morning, 200U DpnII were added again and the sample was incubated for 4 more hours at 37 C, afterward DpnII was inactivated for 20min at 65 C. Digested nuclei were pelleted by centrifugation for 5 min at 600 g and resuspended in 5.7mL milli-Q H2O. Then, 700 mL T4 ligation buffer (Promega) and 20 mL T4 ligase (Promega, 3 U/mL) were added. The samples were incubated overnight at 16 C. The samples were then digested with 5 mL RNAse A (10 mg/mL) and 10 mL proteinase K (10 mg/mL) for 30min at 37 C and decrosslinked o/n at 65 C. DNA was purified by phenol-chloroform extraction and ethanol precipitation and quantified with the QubitTM dsDNA HS Assay Kit (Thermo Fisher). A total of 5 mg of DNA were sonicated with a Covaris S2 instrument with the following parameters: duty cycle 10%, intensity 4, cycles per burst 200 and time 60 s. Sequencing libraries were prepared as described (Schwartzman et al., 2016). The procedure is based in a nested PCR using oligonucleotides aligned to a specific region functioning as ‘‘bait’’ or ‘‘point of view’’ for capturing contacts. The sequences of the primers used can be found in Table S1. Libraries were sequenced in an Illumina NovaSeq instrument to achieve a depth of >1 million reads per library PE 150bp. UMI-4C analysis was performed using the UMI4Cats R package (Ramos-Rodrı́guez et al., 2021). The reads were aligned on hg19 Human genome assembly. For differential testing the function fisherUMI4C (included in the package) was applied, performing a Fisher’s Exact Test in the experiment region using a bin size of 10kb. HiChIP HiChIP was performed mainly as described (Mumbach et al., 2016) starting with 20 million HepG2 CTL cells. After In situ contact library generation, samples were sonicated in a Covaris M220 with the following settings: duty cycle 10%, PIP 75W, 100 cycles/burst, time 50 . A total of 20ug of anti-H3K4me3 antibody was used to immunoprecipitate chromatin. The experiment was replicated twice and libraries sequenced at the BGI in an Illumina HiSeq 4000 obtaining 200M paired-end 50bp reads per replicate. Sequencing data from HiChIP was analyzed using the HiCUP (Wingett et al., 2015) pipeline version 0.7.3. The result was transformed into a .hic file using juicer-pre software from the juicer (Durand et al., 2016) suite. The code with the exact parameters and specifications of the analysis pipeline are available at the following repository URLs: https://github.com/RMalinverni/Omics_miscellanea/blob/master/HiCUP_script_and_configuration.sh https://github.com/RMalinverni/Omics_miscellanea/blob/master/Hic.pipe.clean_by_chrom.sh https://github.com/RMalinverni/Omics_miscellanea/blob/master/HiC_post_analysis.sh. Use of published datasets The ChIP-Seq profile of macroH2A1 and macroH2A2 (GEO accession GSE58175) has been previously described (Douet et al., 2017). epic2 was used to call enriched regions (Stovner and Sætrom, 2019). In addition, we have used the following datasets available from ENCODE (ENCODE Project Consortium et al., 2013): HepG2 ChIP-Seq profile of H3K27Ac (GEO accession GSM733743), RELA (GEO GSE170794), HepG2 HiC (ENCSR194SRI). The ENCODE HiC datast was processed with Juicer (Durand et al., 2016) and loops identified with the included tool HiCCUPS. From the TCGA Pan-Cancer datasets (Weinstein et al., 2013) we used ‘‘Batch effects normalized mRNA data’’ and samples with a positive (upregulation) or negative (downregulation) score of ‘‘RABIT transcription factor regulatory impact, HiSeqV2’’ (Jiang et al., 2015). TCGA data sets were accessed and downloaded from the UCSC Xena platform (Goldman et al., 2020). The GeneHancer database (Fishilevich et al., 2017) and the TSS positions compiled in refTSS (Abugessaisa et al., 2019) were used to identify and plot gene- enhancer functional interactions. We used the 15 chromatin states as defined by the Roadmap Epigenomics Consortium (Roadmap Epigenomics Consortium et al., 2015). The generation and analysis of RNA-Seq data for the 32 paired samples of hepatoblastoma tumor and adjacent normal liver tissue included in the present study has been described (Carrillo-Reixach et al., 2020, GSE133039). All patients have given their informed consent and the Ethical Committee for Clinical Research of the Hospital Germans Trias I Pujol (Spain) has approved the procedures. QUANTIFICATION AND STATISTICAL ANALYSIS In all bar plots, the height of the bar corresponds to the mean value and the bars indicate the Standard Error. In all box plots, the box signifies the upper (75th) and lower quartiles (25th), the median is represented by a horizontal line within the box and the mean is Cell Reports 39, 110988, June 21, 2022 e6 ll OPEN ACCESS Article represented by a rhombus within the box. The upper whisker extends from the upper hinge to the largest value no further than 1.5 * IQR from the hinge (where IQR is the inter-quartile range, or distance between the first and third quartiles). The lower whisker extends from the hinge to the smallest value at most 1.5 * IQR of the hinge. Data beyond the end of the whiskers are plotted individually as ‘‘outliers’’. The statistical test and comparison used to calculate p-values, number of replicates and their nature, p-values set as sig- nificance level and the number of replicates are reported in each figure and/or figure legend. Most plots were generated using ggplot2 (Wickham, 2016). All heatmaps were generated using pheatmap (https://CRAN.R-project.org/package=pheatmap). In all hetmaps, values are scaled by row unless stated otherwise. ChIP-Seq profiles and genomic features were plotted using karyoploter (Gel and Serra, 2017). e7 Cell Reports 39, 110988, June 21, 2022