MacroH2A restricts inflammatory gene expression in melanoma cancer-associated fibroblasts by coordinating chromatin looping

MacroH2A has established tumour suppressive functions in melanoma and other cancers, but an unappreciated role in the tumour microenvironment. Using an autochthonous, immunocompetent mouse model of melanoma, we demonstrate that mice devoid of macroH2A variants exhibit increased tumour burden compared with wild-type counterparts. MacroH2A-deficient tumours accumulate immunosuppressive monocytes and are depleted of functional cytotoxic T cells, characteristics consistent with a compromised anti-tumour response. Single cell and spatial transcriptomics identify increased dedifferentiation along the neural crest lineage of the tumour compartment and increased frequency and activation of cancer-associated fibroblasts following macroH2A loss. Mechanistically, macroH2A-deficient cancer-associated fibroblasts display increased myeloid chemoattractant activity as a consequence of hyperinducible expression of inflammatory genes, which is enforced by increased chromatin looping of their promoters to enhancers that gain H3K27ac. In summary, we reveal a tumour suppressive role for macroH2A variants through the regulation of chromatin architecture in the tumour stroma with potential implications for human melanoma.

MacroH2A has established tumour suppressive functions in melanoma and other cancers, but an unappreciated role in the t um our m ic ro en vi ro nment. Using an autochthonous, i mm un oc om petent mouse model of melanoma, we demonstrate that mice devoid of macroH2A variants exhibit increased tumour burden compared with wildtype counterparts. MacroH2Adeficient tumours accumulate immunosuppressive monocytes and are depleted of functional cytotoxic T cells, characteristics consistent with a compromised antitumour response. Single cell and spatial transcriptomics identify increased dedifferentiation along the neural crest lineage of the tumour compartment and increased frequency and activation of cancerassociated fibroblasts following macroH2A loss. Mechanistically, macroH2Adeficient cancerassociated fibroblasts display increased m ye lo id c he moattractant activity as a consequence of hyperinducible expression of inflammatory genes, which is enforced by increased chromatin looping of their promoters to enhancers that gain H3K27ac. In summary, we reveal a tumour suppressive role for macroH2A variants through the regulation of chromatin architecture in the tumour stroma with potential implications for human melanoma.
Among the downregulated genes, those associated with muscle and keratinization ( Fig. 2b and Extended Data Fig. 2b) indicated dis placement of normal epidermal and subcutaneous muscle structures in the dKO samples. Also downregulated were gasdermins, mediators of pyroptosis, an immunogenic form of cell death 32 , together with the CD200 axis, which inhibits myeloid cell function 33,34 . Importantly, dKO tumours downregulated markers of effector CD8 + cells, including Cd8a, Ifng and Gzmc. Overall, the transcriptomics data suggest that there is an impaired antitumour immune response in dKO animals.

Tumour-infiltrating immune cell dysfunction in dKO mice
We immunophenotyped tumours at 50 DPI (Extended Data Fig. 2c) and observed expansion of classical (CCR2 + ) and nonclassical (CCR2 − CX3CR1 + ) monocytes in dKO samples ( Fig. 2c and Extended Data Fig. 2d,e). These cells, often termed mononuclear phagocytelike myeloidderived suppressor cells, inhibit T cell function and have nega tive prognostic value in cancer 35 . Accordingly, there was a significant reduction in the relative abundance of CD8 + T cells (Fig. 2c) and their proliferation (Fig. 2d). Other CD8 cell markers and their counterparts in CD4 + cells were not affected (Extended Data Fig. 2f). Increased PD1 ligands in dKO immune cells ( Fig. 2e and Extended Data Fig. 2g), prob ably stemming from the myeloid compartment to which PDL2 is mainly restricted 36 , further indicated immunosuppression 37 .
By sorting CD8 + T cells and profiling them using RNA sequencing (RNAseq) (Supplementary Table 2), we found downregulation of the G2/M checkpoint and E2F targets, as well as genes associated with response to interferons and viral infection (Fig. 2f,g). Upregulated genes highlighted a pathological CD8 + T cell phenotype, including the RORγt receptor (Rorc), IL17F and IL23R, which are markers of T C 17 cell polarization. This response can be driven by monocytederived IL1β (the receptors of which are upregulated in the dKO CD8 + population; Fig. 2g), IL6 and IL23 (which signal though the upregulated JAK-STAT pathway) 38 (Fig. 2f), and is associated with impaired antitumour activ ity [39][40][41] . Finally, TNF signalling (upregulated in dKO mice; Fig. 2g) in activated melanomainfiltrating CD8 + T cells triggers their death 42 .
We did not observe significant changes in the relative proportions of immune cells in peripheral blood, regardless of tumourbearing status (Extended Data Fig. 2h,i), which suggested that immunophe notypic changes were TMEspecific. The increase in monocytes and decrease in CD8 + T cell abundance and functionality was recapitulated at 35 DPI (Extended Data Fig. 2j,k), which signified that the dKO immu nophenotype is not a consequence of tumour size. Together, our data indicate that without macroH2A, increased protumour inflammatory signals in the TME inhibit immunemediated tumour cell killing, which facilitates tumour growth.

Melanocyte dedifferentiation in dKO tumours
This melanoma model lacks macroH2A in a constitutive manner; thus, the dKO phenotype could stem from several cell types. We performed single cell RNAseq (scRNAseq) of three WT and three dKO melano mas, generating a dataset of ~24,000 highquality cells. We identified Mice deficient for both H2afy and H2afy2, referred to as double KO (dKO) herein, lack manifest developmental abnormalities or cancer development during ageing 10 . However, macroH2A regulates cellu lar states. For example, macroH2A variants impede reprogramming of somatic cells towards pluripotency 11-13 and mediate gene expres sion in response to stimuli, such as proinflammatory signals 7,14-16 or oncogeneinduced senescence 8 .
Melanoma incidence is rising and remains clinically challenging to treat 17 . We have previously shown that macroH2A expression is lost in advanced human melanoma and, functionally, macroH2A deple tion from melanoma cells promoted enhanced tumour growth and metastatic colonization 18 . Accordingly, overexpression of macroH2A2 induced tumour cell dormancy and suppressed the growth of dissemi nated cancer cells into overt metastasis 19 . Furthermore, high macroH2A levels correlate with favourable prognosis in lung cancer and in colon cancer 20,21 . Collectively, these data provide support for a tumour sup pressive role for macroH2A.
Epigenetic regulation of the melanoma tumour microenviron ment (TME) remains poorly characterized. Here we investigate the consequences of macroH2A deficiency on autochthonous BRAF V600E / PTENdeficient melanomas. We report an unappreciated level of intratumoral heterogeneity and an impaired antitumour immune response in dKO animals. This phenotype is driven by cancer associated fibroblasts (CAFs), which accumulate in the TME and express high levels of inflammatory genes through increased pro moter-enhancer interactions. Our study highlights a unique role for macroH2A histones in the TME by limiting the proinflammatory properties of the tumour stroma.

MacroH2A suppresses autochthonous melanoma growth
We crossed mice constitutively deficient for macroH2A histones (H2afy and H2afy2 dKO strain) 10 to the Braf CA Pten fl Tyr-CreERT2 trial lelic melanoma strain 22 (Extended Data Fig. 1a) and initiated tumours through the topical application of 4hydroxytamoxifen in wildtype (WT) mice and in dKO mice. Although the tumour area was similar at 25 days postinduction (DPI), by 50 DPI, dKO tumours acquired a significant >40% increase in area and a twofold increase in weight versus WT tumours. The increase was independent of sex ( Fig. 1a-c and Extended Data Fig. 1b) and involved increased vertical growth (Fig. 1d). The dKO tumours displayed accelerated development (Fig. 1e) and progressed beyond 50 mm 2 significantly earlier than WT tumours (Extended Data Fig. 1c).
Histology analyses revealed melanocytic lesions with a pig mented epicentre, loss of pigmentation on the margins and depth, and invasion through subcutaneous muscle as determined by S100 staining (Fig. 1f). MelanA staining was present in pigmented foci within tumours and normal melanocytes in the hair follicle, which suggested that most transformed cells lost melanocytic antigens (Fig. 1f). Differences were not observed in proliferation markers, melanocytic spread into the epidermis (pagetoid scatter) or pig mentation between genotypes (Extended Data Fig. 1d-f). Besides melanophages or isolated disseminated tumour cells in lymph nodes (Extended Data Fig. 1g), or rare nonpigmented S100negative lung lesions (Extended Data Fig. 1h), we did not detect overt metastases at 50 DPI. Furthermore, Pten WT dKO mice developed nevi following BRAF V600E induction 22 , but did not proliferate or transform during the lifetime of the animals (Extended Data Fig. 1i). MacroH2A1, present throughout normal skin, was retained in melanomas, whereas mac roH2A2, detected primarily at low levels in the hair follicle in normal skin, was focally present in the tumour ( Fig. 1g and Extended Data Fig. 1j). Altogether, the results indicate that macroH2A loss promotes tumour growth, which may occur through mechanisms distinct from melanocytic hyperproliferation. Article https://doi.org/10.1038/s41556-023-01208-7 33 cell clusters, including melanocytes, immune cells and CAFs, as well as rarer cell populations (Fig. 3a). Cell types and states represented by each cluster were annotated on the basis of expression of known lineage markers (Extended Data Fig. 3a), similarity (Fig. 3b) or the most signifi cant clusterspecific genes 43 (Supplementary Table 3). We performed spatial transcriptomics (ST) to visualize the distribution of tumour populations (spot clusters T1-T6) compared with normal regions (for example, epidermis) in their native tissue context ( Fig. 3c and Sup plementary Table 4). ST confirmed the identities we ascribed using scRNAseq (Extended Data Fig. 3b) and the increased tumour area and invasiveness of the dKO tumours (Fig. 3c) as observed above (Fig. 1).
In the tumour compartment, melanocytes represented only ~5% of cells (Fig. 3d); however, ~25% of cells expressed the neural crest (NC) cell marker Sox10 while lacking melanocyte markers or a MITF gene signature derived from a human melanoma scRNAseq analy sis 44 (Extended Data Fig. 3c-e). These 'NC' clusters ( Fig. 3a) expressed genes associated with developmental precursors of mouse melano cytes and Schwann cells 45 , including Ngfr (NC), Foxd3 (migratory NC and Schwann cell precursors (SCPs)) and Dhh (SCPs) (Extended Data Fig. 3c), and occupied distinct tissue niches in the ST dataset (Extended Data Fig. 3b). The NC arrested cluster displayed a transcriptional pro file (Extended Data Fig. 3c) and cell cycle distribution (Extended Data Fig. 3f) that was consistent with growth arrest. The NC clusters Hapln1 and Aqp1 expressed melanoma dissemination and NC cell migration genes. Finally, the NC Zeb2 and NC Ki67 clusters were characterized by the transcription factor (TF) Zeb2, which is involved in the proliferation of melanoma cells before and following dissemination 46 , along with Gfra3, which is associated with a NC stem cell signature and residual disease in human melanoma 47 . The NC arrested and Aqp1 clusters expressed high levels of an AXLdriven programme associated with melanoma invasion 44 , whereas the Zeb2 and Ki67 clusters expressed the highest levels of a melanocyte stem cell signature 48 (Extended Data Fig. 3d,e).
To dissect the relationship between NC cells and melanocytes, we performed cell trajectory analysis (Extended Data Fig. 3g,h). One branch incurred a growth suppressive programme in the NC arrested cluster, whereas a second branch dedifferentiated towards a state resembling the migratory and stemlike stages of NC devel opment (NC Aqp1 and Zeb2/Ki67 clusters; Extended Data Fig. 3h), which are associated with poor prognosis in human melanoma 49 . This trajectory was supported by the ST data. As expected, melanocytes were present in the dermis and in T1 spots replacing normal dermis ( Fig. 3c and Extended Data Fig. 3b). NC arrested cells were most abun dant in T2 spots, whereas migratory NC Aqp1 and Hapln1 cells defined T3 spots. Meanwhile, T4-T6 spots harboured the bulk of NC Zeb2 cells that invaded subcutaneous structures ( Fig. 3c and Extended Data Fig. 3b). Relative to the total, dKO tumours were enriched in NC Zeb2 cells to the detriment of NC arrested and Hapln1 cells (Fig. 3d).  . n WT = 17, n dKO = 18. P < 0.0001. d, Average tumour depth calculated from the tumour area and volume at the end point (50 DPI). n WT = 17, n dKO = 18. For b-d, significance was determined using Mann-Whitney twotailed test. Box plot whiskers represent the minimum to maximum range, the box plot limits the 25th to 75th percentiles, and the centre line the median. P = 0.0001. e, Tumour growth kinetics between 25 and 50 DPI. n WT = 22, n dKO = 22. Mean and 95% confidence interval error bars are shown. P values adjusted for multiple comparisons: *P < 0.05, **P < 0.01, Mann-Whitney twotailed test. Exact P values are provided as numerical source data. f, Immunohistochemical characterization of normal dorsal skin and representative tumours in a. Antigens indicated are stained pink (Vector Red substrate). g, Immunohistochemical analysis as in f, but demonstrating macroH2A1 and macroH2A2 staining in normal skin and in WT and dKO melanoma. For f and g, insets are shown at additional ×4 magnification. Staining was repeated on n WT skin = 2, n WT melanoma = 7, n dKO melanoma = 6 mice with similar results. Scale bars, 100 μm (f,g) or 1 cm (a). NS, not significant. Article https://doi.org/10.1038/s41556-023-01208-7 This result was confirmed using an unbiased approach that leveraged local cell abundance changes across conditions (Fig. 3e,f). Together, these data suggest that macroH2A deficiency promotes dedifferentia tion in the NC compartment.

Pro-tumour features of myeloid cells in the dKO TME
We identified multiple clusters of mononuclear phagocytes including macrophages (Mac), several of which expressed genes associated with protumour subtypes, including CD206 (Mrc1), Arg1, Retnla, Fn1 and C1qa 50,51 (Fig. 3a,b and Extended Data Fig. 3i). By annotating against two mouse tumour model scRNAseq datasets, we found that Mac Mrc1, which accumulated in the dKO (Fig. 3d), was fully encompassed by Mac_s2, tumourenriched macrophages that become depleted following immune checkpoint blockade 52 , and by Mac1, the human counterpart of which is associated with poor prognosis in patients with lung adenocarcinoma 53 (Extended Data Fig. 3j,k) Response  to  virus   Other  immune  e ectors   WT  dKO  WT  dKO  Il13  Rorc  Il5  Il23r  Ifngr2  Il17f  Il1r1  Il1r2  Il1rl1  Lif  Socs3  Ifnar2  Pik3r1  Ccnd3  Socs4  Jun  Junb  Tradd  Creb3l2  Pgam5   Ifit1  Ifit2  Ifit3  Eif2ak2  Adar  Star  Ifih1  Samd9l  Epsti1  Trim12c  Trim21  Cxcr4  Il12b  Hmga1b  Slfn9  Oasl2  Ifi214  F2rl1  Atp7a  Serinc3 Z-score  Immune   reinforce the finding that immunosuppressive myeloid cells accumu late in dKO tumours. Lymphoid clusters recapitulated the decrease in cytotoxic T (T c ) cells in dKO tumours (Fig. 3d), as observed by flow cytometry (Fig. 2c). Through reclustering, we refined lymphoid cells such that the T c cell cluster split into a Cd4positive CD4 circulating cluster (Extended Data Fig. 3l,m) and a bona fide CD8 cluster negative for Cd4, which was locally depleted in dKO tumours (Extended Data Fig. 3n). Overall, scRNAseq corroborated our bulk RNAseq and immunophenotyping data, providing support for the presence of increased immunosuppres sive myeloid and decreased T c cell infiltration (Extended Data Fig. 3o) in dKO melanomas.

CAFs produce pro-inflammatory signals in the dKO TME
We identified four CAF clusters (Extended Data Fig. 4a), which lacked distinction between inflammatory (iCAF), myofibroblastic (myCAF) and antigenpresenting CAFs (apCAF) reported in pancreatic cancer 54 (Extended Data Fig. 4b,c). This result suggests that there is distinct CAF origin or functional specialization across tumours. The CAF Meg3 cluster exhibited an almost threefold increase in dKO tumours (Fig. 3d), and with NC Zeb2, this cluster was the most relatively enriched cell type (Fig. 3e,f). By computationally assessing the weight of each cluster, CAF Meg3 was highlighted as the top driver of the dKO transcriptional profile (Fig. 4a). Moreover, we found a significant upregulation of the 'dKO tumour up' signature, which consisted of all upregulated genes in    Table 3. c, Significant hallmark pathways in a GSEA of dKO versus WT samples performed in the CAF Meg3 cluster. d, Heatmap of DEGs in CAFs sorted by flow cytometry from WT and dKO melanomas at 50 DPI, grouped under selected top gene enrichment terms defined using Homer analysis. Each column represents CAFs from an independent tumour. Expression values are normalized row wise as Zscores. e, Significant hallmark pathways in GSEA of sorted CAFs as in d. f, Expression normalized to housekeeping controls of indicated cytokine genes determined by reverse transcriptionqPCR in cultured CAFs isolated from WT and dKO tumours at the indicated times following serum stimulation. Line represents the mean of three independently performed experiments shown. Ratio paired twotailed ttest P values shown: *P < 0.05, **P < 0.01. Exact P values are provided as numerical source data. Nonsignificant differences are not labelled.  Table 2), across all dKO CAF clusters, as well as a subset of myeloid clusters (Extended Data Fig. 4d). Importantly, the upregulated cytokines in the bulk RNAseq data (for example, Ccl2, Cxcl1, Ccl11 and Il6) were significantly overexpressed in dKO CAF clusters (Fig. 4b). We additionally found upregulation of other immediateearly genes (for example, Jun and Fos), which was indicative of signal response pathway activation. This finding was confirmed by comparing our data to a signature comprising 139 immediateearly genes 55 (Extended Data Fig. 4e). GSEA revealed upregulation of inflam matory pathways, led by 'TNFα signalling via NFκB' across all CAF clusters ( Fig. 4c and Extended Data Fig. 4f), whereas NC Zeb2, the next highest cluster (Fig. 4a) did not (Extended Data Fig. 4g). Together with the increased CAF Meg3 prevalence (Fig. 3d), these data suggest that CAFs are the primary source of the abovementioned cytokines in the TME and promote the dKO immunophenotype. Of note, whereas H2afy expression was readily detected across clusters, H2afy2 was limited to CAFs (Extended Data Fig. 4h), which suggested that its loss contributes to a CAFspecific phenotype.
Next, we sorted CAFs from WT and dKO tumours by flow cytom etry on the basis of CD140a (Pdgfra) expression (Extended Data Fig. 4i). RNAseq analyses confirmed the upregulation of cytokines, among others, and activation of the TNF-NFκB pathway (Fig. 4d,e and Supplementary Table 5). Differentially expressed genes (DEGs) in the sorted CAFs significantly overlapped with those identified by scRNAseq in the combined mesenchymal populations (Extended Data Fig. 4j). We established pure cultures of sorted WT and dKO CAFs (Extended Data Fig. 4k-m). Given the inducible nature of cytokines, we stimulated primary CAF cultures with serum 56-58 , followed by quantitative PCR (qPCR) and Luminexbased quantitation ( Fig. 4f and Extended Data Fig. 4n). dKO CAFs expressed higher levels of Ccl2, Cxcl1 and Il6 at baseline, which were further induced after serum stimulation. In immortalized dermal fibroblasts (iDFs) derived from nontumour bearing WT skin and dKO skin 11 , similar results were observed (Extended Data Fig. 4o), which highlighted a conserved cellintrinsic mechanism.

Conserved macroH2A-dependent cytokine regulation in CAFs
CAFs can recruit immunosuppressive myeloid cells by secreting CCL2, IL6 and CXCL1 (reviewed in ref. 59). Thus, we performed ligand-recep tor analysis of the scRNAseq dataset. CAFs had the most prolific outgo ing interactions with other cell types both in the WT and dKO samples (Extended Data Fig. 5a), as previously described 60 . Differential analysis showed a generalized decrease in the number of interactions in the dKO samples, but an increase in the strength of communication from CAFs to NC cells (Extended Data Fig. 5a). Importantly, dKO CAFs increased signalling interactions to the mononuclear phagocyte lineage through the CCL2 and IL6 pathways, and to neutrophils and basophils through CXCL1 (Fig. 5a). Our ST data revealed that CAF Meg3 cells formed a distinct layer (Fig. 5b), partially overlapping with subcutaneous spot types (Extended Data Fig. 3b). Mac Mrc1enriched spots were enriched at the tumour periphery (Fig. 5b), and correlation analysis revealed a significant positive association between Mac Mrc1 and Meg3, Fbln1 and Lrrc15 CAF subtypes (Fig. 5c), which suggested proximity. Importantly, these cell types shared a significant negative correlation with T c cells, a characteristic consistent with local T cell exclusion (Fig. 5b,c), which was probably driven by Mrc1 + myeloid cells. We confirmed the chemoat tractant properties of CAFs in vitro by measuring the migration of WT bonemarrowderived monocytes towards WT or dKO CAFs through Transwell assays. MacroH2A dKO CAFs displayed significantly higher monocyte recruitment at later time points (Fig. 5d and Extended Data Fig. 5c), which bolstered our finding of increased monocytederived cells in the dKO TME.
Next, we predicted immune cell abundance in macroH2A high and low tumours from the melanoma cohort of The Cancer Genome Atlas (TCGA SKCM). MacroH2A1 low and macroH2A2 low primary and metastatic samples were associated with significantly reduced CD8 T cell scores ( Fig. 5e and Extended Data Fig. 5d). In primary tumours, M2 (protumour) macrophages were significantly associated with macroH2A2 low and trended for macroH2A1 low tumours (Fig. 5e). Some myeloid subtypes were negatively correlated with macroH2A2 in meta static tumours (Extended Data Fig. 5d), although macroH2A low tumours appeared overall depleted of immune cells (Extended Data Fig. 5e), which affected our ability to detect relative increases in immune sub types. We also identified key T c cell marker genes, CD8A and IFNG, and components of the tumour cytolytic activity score, GZMA and PRF1 (ref. 61), as significantly downregulated in macroH2A low tumours (Extended Data Fig. 5f). These correlations highlight antitumour immunity dysfunction in human macroH2A low melanomas.
Examining human melanomaderived primary CAF cultures (Extended Data Fig. 5g) revealed homogenous levels of macroH2A1 protein, whereas macroH2A2 spanned almost an order of magnitude ( Fig. 5f and Extended Data Fig. 5h). MacroH2A2 low CAFs secreted more CCL2, CXCL1 and IL6 when stimulated (Fig. 5g). Although not signifi cant owing to the low sample size (Extended Data Fig. 5h), we analysed a large pancancer scRNAseq dataset 60 comprising over 56,000 CAFs across 98 samples. The majority of CAFs had no detectable MACROH2A2 counts (Extended Data Fig. 5i), so pseudobulk expression values per tumour were calculated. Whereas CCL2, CXCL1 and IL6 were positively correlated to each other (Extended Data Fig. 5j), the negative correla tion between IL6 and MACROH2A2 was significant (Extended Data Fig. 5j). This result suggests that the relationship between macroH2A2 and inflammatory signalling is conserved in human CAFs. a, Comparison of signalling probability along CCL, CXCL and IL6 pathways leveraging scRNAseq data from CAF Meg3 cells to myeloid cell clusters. Dots represent significant ligand-receptor interaction pairs with increased signalling in the dKO. The dot colour represents communication probability, the dot size represents the P value of onesided permutation test, and the absence of a dot signifies a null probability of signalling. Exact P values are provided as numerical source data. b, Prediction of the spatial localization of indicated scRNAseq cell clusters in WT and dKO melanoma by label transfer onto ST data. Insets are shown at ×2 magnification. c, Correlation analysis of celltype scores for all scRNAseq clusters detected in SC data, based on the combined set of WT and dKO spots. Dots shown correspond to significant correlations (twotailed ttest adjusted for multiple comparisons, adjusted P < 0.05), heatmap colour corresponds to Pearson's correlation coefficient. Exact P values are provided as numerical source data. Black squares represent hierarchical clusters of cell types based on correlation coefficients. d, Transwell assay measuring the migration of CMFDAlabelled WT bonemarrowderived monocytes towards unlabelled WT or dKO CAFs. Monocyte counts are normalized to the CCL2 condition at 24 h. Summary of three independent experiments using different monocyte donors shown. Error bars represent s.e.m. Twotailed ttest P values shown: *P < 0.05, **P < 0.01. Exact P values are provided as numerical source data. Nonsignificant differences are not labelled. e, Comparison of deconvoluted immune cell type scores between TCGA primary melanoma tumours with MACROH2A1 and MACROH2A2 high and low expression levels. Wilcoxon ranksum test P values adjusted for multiple comparisons shown: *P < 0.05, **P < 0.01. Exact P values are provided as numerical source data. N = 35 biologically independent samples per category. The box plot centre line represents the median, the box plot limits indicate the 25th to 75th percentiles, the whiskers extend from the box limit to the most extreme value no further than 1.5× the interquartile range from the box limit, any data beyond whiskers are plotted as individual points. f, MacroH2A1 and macroH2A2 levels in a panel of 11 human melanoma CAF lines. Histone H3 was used as a control for total histone content. g, Indicated cytokine levels in CAF lines in f, stratified according to macroH2A2 levels along the median. N high = 6, n low = 5 biologically independent samples. The western blot was repeated three times.

MacroH2A-regulated genes are enriched in super-MCDs
CUT&RUN analyses of macroH2A1 was performed in WT cultured CAFs (Extended Data Fig. 6a). As previously reported 62 , macroH2A was excluded from the bodies of highly expressed genes and retained at lowly expressed ones (Extended Data Fig. 6b). We identified mac roH2A1 chromatin domains (MCDs) 62 , which we partitioned into 'super' and 'standard' classes on the basis of enrichment and size (Extended Data Fig. 6c-e). Genomewide, macroH2A1 was enriched at signifi cant DEGs compared with a control set of static genes with matched expression (Fig. 6a). These DEGs also significantly overlapped MCDs (Fig. 6b). Notably, the 39 significantly upregulated inflammatory genes (Supplementary Table 5) showed higher average macroH2A1 occupancy compared with static genes (Fig. 6c).
Intergenic enrichment of macroH2A suggested that it may regulate cis-regulatory elements. Traditional enhancers (TEs) and superenhanc ers (SEs) are emerging regulators of NFκBdriven inflammatory gene transcription 63-67 . Moreover, macroH2A regulates specific promoterenhancer contacts 68 and suppresses a subset of enhancers 69 . There fore, we performed chromatin immunoprecipitation with sequencing (ChIP-seq) for H3K27ac, which marks active promoters (Extended Data
We next utilized iDFs to query whether macroH2A regulates induc ible genes beyond CAFs. We performed RNAseq before and after serum stimulation and found that TNF signalling through the NFκB pathway was significantly upregulated in the dKO samples following stimulation (Extended Data Fig. 8a and Supplementary Table 6). The upregulated serumresponsive genes (n = 139; Extended Data Fig. 8b,c) were pref erentially enriched in macroH2A1 and macroH2A2, as measured in WT DFs 62 , relative to static genes (Fig. 6f). Generally, DEGs and DAEs were enriched in macroH2A variants in WT fibroblasts compared with static regions, regardless of the change in direction (Extended Data

MacroH2A loss leads to rewired chromatin looping
Emerging evidence suggests that macroH2A regulates 3D genome organization 68,73-76 . To assess this possibility, and to annotate functional promoter-enhancer pairs, we performed MicroC 77 coupled with pro moter capture (pcMicroC) 78 in CAFs (Extended Data Fig. 9a). At 10 kb resolution, we identified a similar number of promoteroriginating loops in WT and dKO CAFs, with similar size and score distributions ( Fig. 7a and Extended Data Fig. 9b). Distal loop ends were enriched for open chromatin and active enhancers (Extended Data Fig. 9c), which validated the functional nature of the contacts.
A comparison of loop coordinates revealed moderate overlap between WT and dKO CAFs (Fig. 7a), which suggested that there was genomewide reorganization of promoter contacts. However, these changes were also driven by small shifts in the distal ends of loops to an adjacent 10 kb bin (Extended Data Fig. 9d). By examining the activity of enhancers located at distal loop ends, we found that gain of H3K27ac in the dKO samples was more frequently associated with unique dKO loops (and loss of H3K27ac with unique WT loops) compared with static enhancers ( Fig. 7b and Extended Data Fig. 9e). This association in the dKO samples was highly significant for enhancers within superMCDs (Fig. 7b), which suggested that macroH2A instructs functional looping properties of the chromatin fibres it is highly enriched in. Similarly, by   Article https://doi.org/10.1038/s41556-023-01208-7 associating changes in the total number of loops originating in each gene with DEGs, the overlap reached the highest significance and level of directional correlation in macroH2A superMCDs ( Fig. 7c and Extended Data Fig. 9f). Because loop changes often did not match those in gene or enhancer activity (Fig. 7b,c), we intersected genes upregulated in dKO CAFs with bona fide target genes of enhancers gaining H3K27ac, as determined by looping data. Although only 20 genes fit these stringent criteria, 17 had a net increase in the number of loops in dKO CAFs, and 7 were part of inflammatory signalling (Fig. 7d), with significant enrich ment for NFκB targets (Supplementary Table 7). Loci such as Ccl2, Il6, Cxcl1, Ptgs2 (the predominant prostaglandinendoperoxide synthase in CAFs) and Kitl (which encodes Kit ligand/stem cell factor) were located within superMCDs ( Fig. 7e- Altogether, these results demonstrate that macroH2A represses inflam matory genes in CAFs by restricting enhancer contacts and/or activity.

Discussion
Our study of a macroH2Adeficient melanoma mouse model revealed an unappreciated role for this histone variant in the TME. By profiling macroH2Adependent chromatin looping, we identified widespread changes in the promoter-enhancer interaction landscape. Together with changes in enhancer activity, this finding suggests that macroH2A may enforce the position of enhancers relative to nuclear compart ments or the ability of enhancers to interact with their cognate pro moters. Accordingly, previous studies have suggested that macroH2A affects 3D chromatin organization at several scales, including stabi lizing nucleosome-DNA contacts to limit the mobility of the chro matin fibre 73,76 , associating with the boundaries of laminaassociated domains and promoting heterochromatin 74,75 , blocking BRD4 binding at macroH2Abound enhancers 69 , and changing contact frequency of promoter-enhancer pairs, the activity of which is affected by mac roH2A1 or macroH2A2 depletion 68 . Interestingly, however, changes in chromatin accessibility were minimal, which suggests that macroH2A loss does not affect chromatin remodelling, as previously reported 79 .
In our model, macroH2Adependent regulatory mechanisms converged on a small set of inflammatory genes, which underscores a role for macroH2A in limiting inflammatory signalling in vivo. Nota bly, CAFs hijack these inflammatory genes as a mechanism of tumour immune escape 80-83 , and our studies suggest that dKO CAFs promote an immunosuppressive environment that leads to increased melanoma burden (Fig. 7h).
MacroH2A deficiency in cancer promotes tumour growth through multiple mechanisms 1 . Our previous report that macroH2A blocks the proliferative and metastatic capacity of melanoma cells 18 aligns with the increased size of primary tumours observed in dKO mice. Here, we revealed that macroH2A suppresses dedifferentiation along the NC lineage towards a state associated with advanced disease and poor prognosis 44,47,49 , immune evasion and immunotherapy resistance 84-87 . We uncovered a conserved role for macroH2A in mouse and human melanoma CAFs, which produce increased cytokines when macroH2A levels are low. This phenotype appears CAF intrinsic, although we can not exclude the possibility that a hyperinducible response to stimuli occurs in other macroH2Adeficient cells of the TME. Furthermore, given that macroH2A loss decreases the frequency of the CAF Fbln1 cluster (Fig. 3d), which expresses myofibroblast markers (Extended Data Fig. 4a), and downregulates the myofibroblastassociated genes Lrrc15 and Fbln1 (refs. 88,89) in sorted CAFs (Fig. 4d), we propose that the increased inflammation observed is a consequence of skewed der mal fibroblast polarization towards iCAFs at the expense of myCAFs 54 .
Inflammatory signalling was among the first identified microenvi ronmental cues that induce melanoma dedifferentiation 84,90,91 , which raised the possibility that CAFs not only recruit immunosuppressive myeloid cells but may also promote melanocyte dedifferentiation. Such crosstalk occurs in colorectal cancer, in which Ptgs2 expression in CAFs drives the expansion of tumourinitiating stem cells in a paracrine man ner 92 . We speculate that the convergence of these mechanisms would predict poor response of macroH2A low tumours to immunotherapy, with potential to stratify patients.

Online content
Any methods, additional references, Nature Portfolio reporting sum maries, source data, extended data, supplementary information, acknowledgements, peer review information; details of author contri butions and competing interests; and statements of data and code avail ability are available at https://doi.org/10.1038/s41556023012087. Bars under macroH2A CUT&RUN tracks indicate 'super' (purple) and 'standard' (magenta) macroH2A chromatin domains. Below H3K27ac tracks, bright and dark bars indicate TEs and SEs, respectively; red, blue and green denote gain, loss and no change, respectively of H3K27ac level in dKO versus WT CAFs. Chromatin loops, originating at the promoter of the highlighted gene, are shown in red if specific for the dKO, blue for the WT, and black if shared. f, As in e, but for the Il6 locus. g, As in e, but for the Ptgs2 locus. h, Model of the impact of macroH2A loss on the melanoma TME. In the absence of macroH2A, inflammatory genes in CAFs become intrinsically hyperinducible owing to increased enhancer activity and promoter-enhancer looping. This leads to an increased production of proinflammatory cytokines by CAFs, which in turn attract Mrc1 + myeloid cells with a protumour phenotype. Accumulating myeloid cells inhibit CD8 + Tcellmediated tumour cell killing, which results in increased tumour size in dKO animals. CAFdriven inflammatory signalling could also promote melanoma dedifferentiation through mechanisms that are yet to be determined (dashed lines). Illustration in h by Jill K. Gregory, reproduced with permission from © Mount Sinai Health System.    Fig. 1b). For tumour growth and histology experiments, Cremediated recombination of Braf and Pten alleles was induced in 7-11weekold mice through the application of 1 μl of 5 mM 4hydroxytamoxifen (70% Zisomer, Sigma H6278, dis solved in ethanol) on the medial dorsal skin 24 h after hair removal with depilatory cream. Tumour length and width were measured from 25 DPI onwards using calipers, and the area was calculated assuming an ellipti cal tumour shape. Tumour thickness was approximated from the end point tumour area and weight, assuming an elliptic cylinder shape and tissue density = 1. Tumours were collected no later than 50 DPI before reaching a humane end point (tumour length over 1 cm, presence of ulceration). Five WT and 10 dKO male and 7 WT and 5 dKO female mice aged 9-10 weeks at induction were used for tumour immunopheno typing at 50 DPI. Three WT and 5 dKO male and 5 WT and 4 dKO female mice aged 9-10 weeks at induction were used for peripheral blood immune cell counts at 50 DPI. Two WT and 1 dKO male and 4 WT and 5 dKO female mice aged 10-11 weeks at induction were used for tumour immunophenotyping at 35 DPI. Two WT and 3 dKO male and 4 WT and 3 dKO female mice aged 7 weeks were used for peripheral blood immune cell counts in the tumournaive setting. Two females and 2 males each for WT and dKO mice aged 9 weeks at induction were used for CD8 T cell sorting at 50 DPI for RNAseq. Two females and 2 males each for WT and dKO mice aged 10 weeks at induction were used for CAF sort ing at 50 DPI for RNAseq. Two WT and 2 dKO females aged 10 weeks at induction were used for CAF sorting at 50 DPI for ATAC-seq. One WT and 1 dKO male aged 13 weeks at induction were used for CAF sorting at 50 DPI for culturing. Females were used for scRNAseq to avoid the potential impact of sexspecific transcripts on integration. Pairs of WT and dKO agematched females were induced at 10, 12 and 9 weeks in 3 independent experiments. One WT and 1 dKO female were induced at 9 weeks and processed at 35 DPI for ST analyses.

Histology
Resected tissue was fixed in neutral buffered formalin for 24 h, paraffin embedded and sectioned at the ISMMS Biorepository and Pathology Core. Antigens were retrieved in citrate buffer (pH 6) in a domestic pressure cooker for 5 min. Immunodetection was performed using ImmPRESS polymer, ImmPACT Vector Red and NovaRed kits (Vector). Primary antibodies are listed in Supplementary Table 8. Ten random ×40 objective fields within the tumour were scored by boardcertified dermatopathologists at ISMMS (M.S.G. and N.S.V.) for mitotic cells (H3S10ph) and Ki67. Pagetoid scatter was measured as the number of melanoma cells within the epidermis across ten ×40 fields. The degree of pigmentation was estimated on sections stained with haematoxylin and eosin as the fraction of tissue area containing pigment using a ×2 objective.

Tumour dissociation
Resected BRAF V600E /PTENdeficient melanomas were cut into 1-2 mm fragments. For immunophenotyping and CD8 + T cell sorting, tumour fragments were digested in RPMI containing 400 U ml -1 collagenase IV (Gibco), 100 U ml -1 hyaluronidase (Sigma) and 100 μg ml -1 DNAse (Roche) at 37 °C for 1 h 96 , aspirated 5 times through a 14 G needle, and filtered through a 70 μm cell strainer. Immune cells were enriched by centrifugation through a discontinuous 40/90 Percoll (GE Healthcare Life Sciences) gradient 97 . For scRNAseq and CAF sorting, tumour frag ments were digested using a Tumour Dissociation kit, mouse (Miltenyi) in DMEM using the soft/medium protocol according to the manufac turer's instructions. For scRNAseq specifically, red blood cell lysis was performed in ACK (ammonium-chloride-potassium) buffer for 5 min on ice, followed by a wash in 1× DPBS containing 0.04% BSA.

RNA extraction
Snapfrozen tumours were disrupted in QIAzol reagent (Qiagen) by milling with zirconia beads (Fisher). Sorted cells were centrifuged and resuspended in QIAzol. Cultured cells were lysed directly in their culture vessel. After the addition of chloroform (Sigma) and centrifu gation, RNA was isolated from the aqueous phase using an RNEasy

RNA-seq analysis
Reads were quasimapped to the Gencode M25 (GRCm38.p6) gene set using salmon (v. by filtering for an independent hypothesis weighting 100 adjusted P value of <0.05 and a log 2 fold change of >0.75 or < -0.75. GSEA was performed with the fgsea (v.1.22.0) package 101 on the entire expressed gene set preranked based on the Wald test statistic computed using DESeq2. Significant pathways were reported using an adjusted P value cutoff of 0.05 and ordered on normalized enrichment scores. Gene ontology enrichment was performed using HOMER (v.4.10) 102 . Heatmaps were generated for visualization purposes using counts normalized with the variance stabilization transformation of DESeq2 and corrected for library preparation batch or sexspecific effects when applicable using limma (v.3.54.0) 103 .

Reverse transcription-qPCR
A First Strand cDNA Synthesis kit (OriGene), FastStart Universal SYBR Green Master (Rox) mix (Roche) and primers listed in Supplementary  Table 7 were used for RTqPCR. Amplification was performed on a CFX384 instrument (BioRad), and target genes were quantified rela tive to Hprt and Sdha housekeeping genes using CFX Manager (v.3.1) software (BioRad).

scRNA-seq library preparation
Dropletbased scRNAseq (Chromium, 10x Genomics) was applied to singlecell suspensions of BRAF V600E /PTENdeficient melanomas from 3 WT and 3 dKO mice; 6 × 10 3 cells were loaded per well for each of the 6 samples, and partitioning and library preparation were performed according to the manufacturer's protocol for the 3′ v.3 chemistry.

scRNA-seq analysis
Filtered gene-barcode matrices generated using cellranger (v.7.0.1) with the mm102020A transcriptome (10x Genomics) were further analysed using Seurat (v.4.0) 104 . Lowquality cells (fewer than 500 unique molecular identifiers or 250 detected genes, or at least 20% mitochondrial transcripts) were removed, and data were normalized using SCTransform (v.2), regressing out cell cycle scores (determined using the CellCycleScoring function) and mitochondrial ratio. Data sets were integrated using reciprocal principal component analysis on the first 50 principal components with an alignment strength of ten. Within this combined dataset, clusters were determined using a variety of resolutions and annotated using top clusterspecific genes conserved between WT and dKO cells, as well as known lineage genes for expected cell types. At final resolution (0.6), selected to avoid overclustering while distinguishing rare cell types, five clusters with a low number of significant genes, also displaying a low nuclear tran scriptome complexity, were considered to represent lowquality cells and were discarded. The least abundant cluster, characterized by genes of multiple lineages, probably contained doublets and was also discarded. dKO versus WT DEGs were identified within each cluster on the RNA slot using the FindMarkers function in Seurat. GSEA was performed using SCPA (v.1.5.1) 105 . Reclustering of related cell types was performed by subsetting the Seurat object to the relevant clusters, then performing normalization, integration and clustering as described above. Pseudotime trajectories for NC lineage clusters were calculated using Monocle 3 (v.1.3.1) 106 , ordering cells in pseudotime under the assumption that they shared a common precursor and transformation initiated in the highest Tyrexpressing cluster. Unbiased cell identity mapping to existing scRNAseq datasets was performed by subsetting Seurat objects to myeloid cell clusters, then determining the highest similarity in terms of gene expression to the reference cell types using singleR (v.1.10.0) 107 . Local changes in cell abundance were profiled using miloR (v.1.5.0) 108 with the following parameters for the entire dataset: 50 principal components, 40 nearest neighbours, sampling proportion 0.1 and sampling refinement algorithm. Thirty nearest neighbours and a sampling proportion of 0.2 were used for the lower number of cells in the lymphoid reclustered dataset. Celltype prioritization to evaluate the contribution to changes in gene expression between WT and dKO samples was performed using Augur (v.1.0.3) 109 using default parameters except for a minimum cell number of 50. Cell signalling through ligand-receptor interaction analysis was performed using CellChat (v.1.6.1) 110 using the comparison workflow on a merged object.

Detailed cluster annotation
In addition to those described in the text, the following genes were used as markers during cluster annotation. Melanocytes were identified through the expression of genes associated with pigment production, such as Mitf, Mlana and Dct 45 . The NC arrested cluster expressed high levels of the cell cycle inhibitors p21(CIP1) (encoded by Cdkn1a) and p19(ARF) (encoded by Cdkn2a), the EGFlike ligands amphiregulin (Areg) and epiregulin (Ereg), and the histone variant H2A.J (H2afj). Amphiregulin expression is associated with BRAF V600E induced senes cence in melanocytes 111 and H2A.J accumulates in senescent fibro blasts 112 . Another NC cluster expressed Hapln1, an ECM crosslinker the downregulation of which in aged fibroblasts promotes melanoma dis semination 113,114 . The NC Aqp cluster was characterized by aquaporin 1 (Aqp1) and Tfap2b, two factors that orchestrate NC cell migration 115 [131][132][133] . Two clusters with high levels of the myofibroblast marker Acta2 expressed Lrrc15 and Fbln1, respectively, both genes associated with immunoregulatory myofibroblast popula tions in pancreatic cancer and in breast cancer 88,89 .

ST analysis
Tissue was prepared according to 10x Genomics Visium V1 Slide-3′ Spatial guidelines. Tumours were frozen in a bath of isopentane and liquid nitrogen, stored at -80 °C in a sealed container, then embedded in OCT. Tissue was sectioned at temperatures of -20 °C for the cryostat chamber and -10 °C for the specimen head at 10 μm thickness. An optimal permeabilization time of 18 min was determined using a Visium Spatial Tissue Optimization kit (10x Genomics). Sequencing data were mapped using spaceranger (v.2.0.0) using the mm102020A transcrip tome (10x Genomics) to generate spatial gene expression matrices, which were processed according to the Seurat spatial vignette (https:// satijalab.org/seurat/articles/spatial_vignette.html). We integrated data normalized with SCTransform (v.2) across the two ST samples. For label transfer, we then integrated this dataset on the first 30 principal components with scRNAseq data normalized using SCTransform (v.2).
The probabilistic distribution of cell types identified by scRNAseq in each ~50 μm spot of the spatial array was calculated using the Transfer Data function. To assess celltype colocalization within the tissue 134 136,137 . All humanderived cell cultures were stripped of any patient identifiers before they were provided to our laboratory, and their use is not subjected to Institutional Review Board approval. Ethics information relevant to their collection is provided in the noted refer ences. Cells were maintained on plates coated with Cultrex Basement Membrane Extract, PathClear (BioTechne) in Advanced DMEM/F12 (Gibco) supplemented with 5% heatinactivated FBS, 2 mM lglutamine (Gibco), 10 ng ml -1 recombinant human EGF (Gibco), 400 ng ml -1 hydro cortisone (Sigma), 24 μg ml -1 adenine (Sigma), 100 μg ml -1 Primocin and 25 μg ml -1 Plasmocin (Invivogen); 10 μM Y27632 dihydrochlo ride (Tocris) was added to the medium during the initial expansion of cultures. For serum starvation before cytokine induction, hydrocor tisone, adenine and EGF were omitted, and the serum concentration was reduced to 0.25%. After 24 h, starvation medium was replaced with complete human CAF medium without Y27632 dihydrochloride. Supernatant was collected, centrifuged to remove cells and debris, and flashfrozen. Cells were counted using a Guava Muse Cell Analyzer Count & Viability kit (Luminex).

Protein quantitation
Western blotting was performed as previously described 62 using the antibodies listed in Supplementary Table 8. For cytokine quantitation in mouse samples, total protein was extracted using RIPA buffer and normalized to the lowest sample concentration. Analytes of inter est were quantified using a Mouse Cytokine Array/Chemokine Array 31Plex by Eve Technologies. For cytokine quantitation in human CAFs, supernatant was diluted 1:5, 1:10 and 1:20 with assay diluent and sub jected to ELISA in duplicate against human CCL2, IL6 (BioLegend) and CXCL1 (R&D Systems) according to the manufacturers' instructions. Data points within range of the standard curve were normalized by cell number and averaged for each sample.

Monocyte Transwell migration assay
Monocyte migration was performed as previously described 138 but with modifications. In brief, 7.5 × 10 4 WT or dKO CAFs were seeded per well of a 24well cell culture insert companion plate (Corning), in triplicate, allowed to grow for 24 h reaching confluency, then starved for 24 h as described above. Monocytes were isolated from femur and tibia bone marrow of tumournaive WT 129S mice using a MojoSort Mouse Mono cyte Isolation kit (BioLegend) and labelled with 3 μM CMFDA green dye in serumfree RPMI 1640 medium for 30 min at 37 °C. After washing with R10 medium (RPMI 1640 supplemented with 10% heatinactivated FBS, 20 mM HEPES, 0.5 mM sodium pyruvate, 1% penicillin-streptomycin, 1× MEM amino acids without lglutamine, 1× MEM nonessential amino acids, Gibco), unbound dye was allowed to diffuse out by another incubation for 30 min at 37 °C in R10. Next, 1 × 10 5 monocytes in 300 μl R10 were added to a FluoroBlok insert with 3 μm pore size (Corning), which functioned as the upper chamber. Concomitantly, starvation medium on CAFs in the companion plate, which functioned as a lower chamber, was replaced with 800 μl R10. One well each of R10 and R10 supplemented with 10 ng ml -1 recombinant mouse CCL2 (BioLegend) were used as negative and positive controls, respectively. Plates were imaged for green widefield fluorescence at the Z position correspond ing to the CAF layer, every 6 h for a 72 h time course, in a Cytation 7 automated microscope (Agilent). Image stitching, background subtrac tion, segmentation and cell count were performed using the onboard Gen5 (v.3.12) software (Agilent). To account for differences in monocyte reactivity among donor mice, cell numbers were normalized to the number of monocytes migrated at 24 h in the CCL2positive control.

TCGA data analysis
TCGA melanoma cohort (SKCM) 139 RNAseq raw counts and sample annotations were downloaded using the TCGAbiolinks (v.2.25.3) pack age 140 . Gene counts were normalized using the TMM method from the edgeR (v.3.40.1) package 141 to account for biases arising from library size and gene length. Primary and metastatic lesions were separately analysed, and comparisons for TME composition and gene expression were performed between the first and third tercile of MACROH2A1 and MACROH2A2 expression. Deconvolution of immune populations 142 and tumour/stroma/immune microenvironmental composition 143 was performed using the IOBR (v.0.99.9) package 144 . Differential gene expression analysis was performed using DESeq2 (ref. 99) with raw counts as input.

Human CAF scRNA-seq reanalysis
scRNA count matrix of 855,271 highquality cells and associated meta data annotations from a previous study 60 (GSE210347) were reanalysed. Counts were scaled using the Seurat LogNormalize function with a scale factor of 10,000. The top 2,000 most highly variable genes were identified using the vst selection method of FindVariableFeatures. The RunFastMNN function from the SeuratWrappers (v.0.3.0) pack age was utilized to perform batch correction using the 'SampleID' metadata column. The top 30 principal components were used in the downstream analysis. To investigate the correlation between MAC-ROH2A1-MACROH2A2 and IL6-CXCL1-CCL2 expression, we generated pseudobulk counts from fibroblasts. In brief, the filtered Seurat object containing annotated fibroblasts was converted into a SingleCellEx periment (v.1.12) object using the raw counts per cell. The function aggregate.Matrix was used to sum and collapse the raw counts of each cell by their 'SampleID' annotation. Samples with <100 cells and nontumour tissue were filtered out. Collapsed raw counts were nor malized using the vst function from DESeq2 and the Pearson's correla tion between genes was calculated using the chart.Correlation function from the PerformanceAnalytics (v.2.0.4) package.

ATAC-seq library preparation
5 × 10 4 sorted CAFs per sample were processed for ATAC-seq as pre viously described 145 . The optimal number of library amplification cycles was determined 146 . Libraries were subjected to doublesided size selection using SPRIselect beads (Beckman Coulter) at ratios of 0.55 and 1.2× before sequencing on a NextSeq 500 (Illumina) in 75 bp pairedend mode.

ChIP-seq analysis
Adapters were trimmed using Trimmomatic (v.0.36) 155 , and alignment, read filtering and genome coverage calculation were performed as for ATAC-seq. The bam files of WT and dKO samples were concat enated into a master bam file, which was used to call significant peaks with matching input controls using MACS2 (v.2.1.0) 152 . Cutoff values for qvalue significance were determined posthoc, testing several qvalues on the basis of the signaltobackground ratio. Quantifica tion of reads in significant peaks for all samples was performed using bedtools (v.2.29.2) 153 multicov. TEs and SEs were called on the basis of H3K27ac enrichment using the ROSE algorithm (rank ordering of superenhancers) (v.0.11) 156,157 with a stitching distance of 12.5 kb and a TSS exclusion zone size of 2.5 kb. The ROSE algorithm was also used to extract H3K27ac levels at TEs and SEs for WT and dKO samples individually. The average H3K27ac signal value across all elements was calculated for each sample and further used for normalization between samples. The log 2 fold change ratio of the normalized signal (-0.75 < log 2 fold change < 0.75) was used to call differential TEs and SEs. All other enhancers were considered static.

CUT&RUN
2 × 10 5 shortterm cultured CAFs before or after 30 min of serum stim ulation were processed for CUT&RUN 158 for macroH2A1 (antibody ab37264, Abcam, lot number GR2780201, 1 μg per reaction). Cells were permeabilized with 0.0085% digitonin and incubated with antibody. DNA was cleaved with the CUTANA pAGMNase (Epicypher). Released DNA was processed using a NEBNext Ultra II Library Preparation kit for Illumina (NEB), including a 14cycle amplification step. Libraries were sequenced and reads were processed as for ATAC-seq.

CUT&RUN analysis
MacroH2A1 enrichment determined by CUT&RUN was benchmarked by correlation analysis with published macroH2A1 ChIP-seq data in dermal fibroblasts 62 . Enrichment at the level of read pileups had a correlation coefficient of 0.79 (Extended Data Fig. 6a). MacroH2A1 chromatin domains were called separately for stimulated and unstimu lated samples using SEACR (v.1.3) 159 without a control in stringent mode with a threshold of 0.05, and significant peaks were merged using bedtools (v.2.29.2) 153 if within 25 kb distance. Merged peak and alignment files for unstimulated and stimulated cells were then concat enated to generate master files for the ROSE algorithm 156,157 . A stitching distance of 12.5 kb, no TSS exclusion and normalization to element size were used in ROSE to rank macroH2A1 domains. The equivalent of SEs called by ROSE based on the macroH2A1 signal were termed superMCDs. The equivalent of TEs was then divided into groups of standard and low macroH2A1 signals using a cutoff above 3,780 (signal density × length units) determined posthoc on the basis of the macroH2A1 enrichmenttobackground ratio. Most domains with low macroH2A1 signal were also below 1.5 kb in length, represented less than 0.5% of the genome, showed higher enrichment of IgG over the macroH2A signal, and were therefore excluded from further analysis.
For comparing macroH2A1 enrichment levels, static genes were selected as follows: for each differentially expressed gene, all genes within 25% of its averaged transcripts per million (TPM) level in WT sam ples were identified, and 10 (inflammatory up genes comparison) or 3 (all DEG comparison) genes that were neither differentially expressed nor less than 1 kb in length were randomly selected, added to a running list, and duplicates were removed.

pcMicro-C
MicroC was performed using a Dovetail MicroC kit (Cantata Bio) according to the manufacturer's protocol. Chromatin from 1 × 10 6 WT and dKO cultured CAFs, crosslinked after serum stimulation, was digested with 3 μl MNase to obtain a mononucleosomal fraction within specifications. Each chromatin preparation was subjected to proximity ligation and subsequent library preparation in duplicate. Next, 400 ng of each library was pooled and subjected to promoter capture using a Dovetail Target Enrichment kit and a Mouse Pan Promoter Panel (Cantata Bio). The manufacturer's protocol was adjusted by halving the streptavidin bead elution volume to enable the use of the entire postcapture material for the amplification reaction and performing only seven PCR cycles. The resulting 4plex library was sequenced on a NextSeq 500 high output lane in a 75 bp pairedend configuration.

Chromatin looping analysis
Interaction calling and significance thresholding was performed based on the workflow developed by Dovetail Genomics (https://microc. readthedocs.io/en/latest/) and the CHiCAGO tool recommendations 160 . Sequenced reads were aligned using BWA (v.0.7.17r1188) 161 and fil tered using PairTools (v.1.0.2) 162 and SAMTools (v.1.9) 149 to identify nonduplicated read pairs. The technical replicates were merged, and significant interactions (chromatin loops) were called using CHiCAGO (v.1.2.0) 163 with default parameters at 10 kb resolution 160 . Significant Interactions were defined as those with a CHiCAGO score of ≥5 per condition. For visualization on the UCSC genome browser, output loop files were converted from the WashU EpiGenome Browser interactBED format to the bigInteract format using BEDTools (v.2.30.0) 153,164 and UCSCUtils (v.2.9) 165,166 . Bait map plots depicting all called interac tions per bait were generated using the plotBaits function native to CHiCAGO. Histograms showing the enrichment of functional elements (ATAC-seq peaks and H3K27ac peaks) at the distal end of loops were generated as a default output of the CHiCAGO pipeline. Interactions were defined as shared if both ends overlapped between conditions, and unique otherwise. Promoters and enhancers at proximal and distal ends of loops, respectively, were divided into classes on the basis of the type of macroH2A1 domain they overlapped with-super, standard and absent-and the nature of the loops underlying them-shared, WTspecific or dKOspecific. For the enhancerspecific association, enhancers were further divided on the basis of the direction of their change in H3K27ac levels by macro level to generate a count matrix of changes in enhancer activity versus changes in interactions, stratified by macroH2A1 level. We performed a Chisquare test of independence to determine associations between enhancer and looping deregulation in the absence of macroH2A. The mouse promoter panel uses multiple baits to capture alternative promoters of the same gene; therefore, we summed loop counts for all baits associated with the same gene name for overlap with gene expression changes. We compared the loop counts per gene in dKO versus WT samples and further stratified genes on the basis of their differential expression status and macroH2A1 occupancy. We performed a Chisquare test of independence to deter mine associations between gene expression and looping deregulation in the absence of macroH2A.

Statistics and reproducibility
Information on the number of biologically independent samples ana lysed and the number of times experiments were performed is included Nature Cell Biology Article https://doi.org/10.1038/s41556-023-01208-7 in the figure legends. No assumption of data normality was made, and nonparametric statistical tests were performed except when n = 3 replicates (comparison of celltype frequency in scRNAseq data, CAF RTqPCR) for which data distribution was assumed to be normal but this was not formally tested. All statistical tests performed were twosided except when noted. No statistical method was used to predetermine sample sizes, but were similar to those reported in previous publica tions 22,97,167 . No data were excluded from the analyses, and biological samples were excluded from the study only if sample preparation or data acquisition failed. The experiments were not randomized, and the investigators were not blinded to allocation during experiments and outcome assessment.

Reporting summary
Further information on research design is available in the Nature Port folio Reporting Summary linked to this article.

Data availability
The transcriptomics and epigenomic datasets, including raw and pro cessed sequencing data generated and analysed during the current study, are available in the Gene Expression Omnibus (GEO) reposi tory under accession number GSE200751 (https://www.ncbi.nlm.nih. gov/geo/query/acc.cgi?acc=GSE200751) and in this article's table files. TCGA melanoma data are publicly available through the NCI Genomic Data Commons (GDC) data portal under project ID TCGA SKCM (https://portal.gdc.cancer.gov/projects/TCGASKCM). The human pancancer scRNAseq dataset mined in this study is available in the GEO repository under accession number GSE210347 (https://www. ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE210347). The mouse M25 (GRCm38.p6) genome assembly and gene set used for transcriptom ics and epigenomic analyses are available at Gencode (https://www. gencodegenes.org/mouse/release_M25.html). All other data sup porting the findings of this study is available from the corresponding authors upon request. Source data are provided with this paper.

Code availability
All packages used for data analysis are publicly available. No custom code was generated for this study. All scripts used for bulk RNAseq, scRNAseq, ATAC-seq, ChIP-seq and CUT&RUN data analyses in this study are available from the corresponding authors upon request.  Fig. 1 | Characterization of macroH2A WT and dKO murine melanomas. a) Breeding strategy to obtain WT and dKO mice used for melanoma induction. b) Tumor area in males and females, n WTM = 17, n WTF = 21, n dKOM = 25, n dKOF = 14. No significant intersex differences were observed within genotypes (KruskalWallis with Dunn's multiple comparisons test). Box plot whiskers represent minmax range, box plot limits -25 th to 75 th percentiles, center line -median. c) KaplanMeier analysis of diseasefree survival. Events represent melanoma growth beyond an area of 50 mm 2 , n WT = 22, n dKO = 22. Pvalue = 0.0009, logrank (MantelCox) test. d) Immunohistochemical scoring of H3S10ph and Ki67 at 50 DPI in ten random fields per tumor. e) Pagetoid spread of melanocytic cells into epidermal structures scored in ten epidermiscontaining fields per tumor on H&Estained sections. f) Relative pigmented area of the tumor, estimated on a 2X magnification ensemble view of H&Estained sections. df, n WT = 11, n dKO = 13, center line represents median, significance determined using a MannWhitney twotailed test. Exact Pvalues are provided as numerical source data. g) Tumordraining axillary lymph nodes stained for S100 (reddishbrown NovaRed substrate color). h) Lung sections of tumorbearing mice at 50 DPI. gh, 4X (scale bar = 1 mm) and 20X (scale bar = 100 μm) magnification shown; inserts shown at 2.5X additional magnification. Staining and analysis were performed in 2 animals per genotype with similar results. i) KaplanMeier analysis of diseasefree survival following nevus induction, n at least 1 WT allele = 7, n dKO = 6. Progression to melanoma, defined as radial and/or vertical lesion growth, was not observed during the lifespan of the mice. j) Validation of antibody specificity for macroH2A variants in immunohistochemistry, using liver of either macroH2A1 or macroH2A2 single KO mice. Scale bar = 100 μm, inserts shown at 5X additional magnification. Experiment was performed twice with similar results.  Peripheral blood, 50 DPI tumor-bearing mice  Trdc Tcf7 S1pr1 Gramd3 Cd4 Ifng Cd8b1 Cd8a Cd3e Itgax Itgam          Arrows depicting signaling direction are colored according to emitting cell type. Arrow weight is proportional to interaction strength. b) Differential interaction strength and number between WT and dKO cell type pairs. Arrows are colored according to direction of change (WT -blue, dKO -red) and weight represents amplitude of change. c) Transwell assay measuring the migration of CMFDAlabeled WT bone marrowderived monocytes towards unlabelled WT or dKO CAFs. Average cell counts of 3 technical replicates and error bars representing SEM are shown for WT and dKO CAF conditions. No replicates are performed for negative and positive controls. Individual experiments using different monocyte donors shown. d) Comparison of deconvoluted immune cell type scores 142 between TCGA metastatic melanoma tumors with MACROH2A1/2 high and low expression levels. n = 123 biologically independent samples per category. e) As in (d) for estimated immune, stromal and tumor purity scores 143 . n primary = 35, n metastatic = 123 biologically independent samples per category. de, Wilcoxon ranksum test Pvalues adjusted for multiple comparisons shown: * < 0.05, ** < 0.01, *** < 0.001, **** < 0.0001. Exact Pvalues are provided as numerical source data. Box plot center line represents the median, box plot limits -25 th to 75 th percentiles, whiskers extend from the box limit to the most extreme value no further than 1.5 * interquartile range from the box limit, any data beyond whiskers is plotted as individual points. f) Expression of genes associated with antitumor cytotoxic activity in human primary and metastatic melanoma samples from the TCGA cohort, segregated by MACROH2A1 or MACROH2A2 gene expression levels. High and low terciles are compared. Independent hypothesis weighted Wald test Pvalues adjusted for multiple comparisons, computed by DESeq2, shown. g) Detection of CAF markers by western blot in a panel of 11 human melanoma primary CAF cultures. Membrane stain for total protein shown for loading control. h) Correlation between macroH2A2 protein levels and indicated cytokine secretion in CAF cultures from (g). n = 11 biologically independent samples. Error bars correspond to SEM of 3 technical replicates of blotbased quantification for macroH2A2 and up to 3 dilutions each in 2 technical replicates for cytokine ELISA. Spearman correlation statistics, calculated on the average values without considering individual technical replicates, are shown. i) UMAP plots of a human pancancer scRNAseq dataset 60 , together with macroH2A gene expression in fibroblasts. j) Correlation analysis between pseudobulk expression levels of selected cytokines and macroH2A genes in CAFs from (i). Pearson correlation coefficients and significance shown above the diagonal, pairwise scatter plots shown below. Histograms of individual gene expression comprise the diagonal. Pvalues shown: ** < 0.01, *** < 0.001. Exact Pvalues are provided as numerical source data.   Fig. 6 | Genome-wide macroH2A2 occupancy in cultured WT CAFs. a) Correlation of genomewide enrichment of macroH2A between ChIPseq and CUT&RUN methodologies. ChIPseq and associated inputs were previously generated in dermal fibroblasts 62 , CUT&RUN was performed in serumstarved unstimulated and serumstimulated WT CAFs, with associated IgG control. Pearson's correlation coefficients shown, conditions are clustered based on Euclidean distance. b) Metagene profile of macroH2A CUT&RUN signal in cultured WT CAFs before and after serum stimulation across genes binned into quartiles according to their expression levels. Genes not detected are either not expressed or not mappable. n = 3000 randomly selected genes within each group. c) Metagene profile of macroH2A CUT&RUN signal in cultured WT CAFs before and after serum stimulation at different classes of MCDs and genomic regions where MCDs are absent. d) Size distributions of regions in (c). Box plot center line represents the median, box plot limits -25 th to 75 th percentiles, whiskers extend from the box limit to the most extreme value no further than 1.5 * interquartile range from the box limit. cd, n super = 1560, n standard = 5440, n low = 1556, n absent = 23615 regions. e) Proportion of the genome comprised within different classes of MCDs.