# HSFA1a modulates plant heat stress responses and alters the 3D chromatin organization of enhancer-promoter interactions

Ying Huang, Jing An, Sanchari Sircar, Clara Bergis, Chloé Dias Lopes, Xiaoning He, Barbara da Costa, Feng-Quan Tan, Jeremie Bazin, Javier Antunez-Sanchez, et al.

<figure><img src="image_1.png" /></figure>

## To cite this version:

Ying Huang, Jing An, Sanchari Sircar, Clara Bergis, Chloé Dias Lopes, et al., HSFA1a modulates plant heat stress responses and alters the 3D chromatin organization of enhancer-promoter interactions. *Nature Communications*, 2023, 14, 10.1038/s41467-023-36227-3. http://04062805

HAL Id: hal-04062805

https://hal.inrae.fr/hal-04062805

Submitted on 7 Apr 2023

HAL is a multi-disciplinary open access archive for the deposit and dissemination of scientific research documents, whether they are published or not. The documents may come from teaching and research institutions in France or abroad, or from public or private research centers.

L'archive ouverte pluridisciplinaire HAL, est destinée au dépôt et à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

<figure><img src="image_2.png" /></figure>

Distributed under a Creative Commons Attribution 4.0 International License
<figure><img src="image_3.png" /></figure>

# HSFA1a modulates plant heat stress responses and alters the 3D chromatin organization of enhancer-promoter interactions

Received: 29 June 2022

Accepted: 19 January 2023

Published online: 28 January 2023

<figure><img src="image_4.png" /></figure>

Check for updates

Ying Huang©¹, Jing An©¹, Sanchari Sircar¹, Clara Bergis¹, Chloé Dias Lopes¹, Xiaoning He¹, Barbara Da Costa¹, Feng-Quan Tan¹, Jeremie Bazin¹, Javier Antunez-Sanchez©², Maria Florencia Mammarella³,

Ravi-sureshbhai Devani©¹, Rim Brik-Chaouche¹, Abdelhafid Bendahmane¹, Florian Frugier¹, Chongjing Xia©⁴, Christophe Rothan⁵, Aline V. Probst©⁶,

Zouine Mohamed⁷, Catherine Bergounioux¹, Marianne Delarue¹, Yijing Zhang©⁸, Shaojian Zheng©⁹, Martin Crespi¹, Sotirios Fragkostefanakis¹⁰,

Magdy M. Mahfouz©¹¹, Federico Ariel³, Jose Gutierrez-Marcos©²,

Cécile Raynaud¹, David Latrasse¹ & Moussa Benhamed©¹,¹²,¹³

<figure><img src="image_5.png" /></figure>

The complex and dynamic three-dimensional organization of chromatin within the nucleus makes understanding the control of gene expression challenging, but also opens up possible ways to epigenetically modulate gene expression. Because plants are sessile, they evolved sophisticated ways to rapidly modulate gene expression in response to environmental stress, that are thought to be coordinated by changes in chromatin conformation to mediate specific cellular and physiological responses. However, to what extent and how stress induces dynamic changes in chromatin reorganization remains poorly understood. Here, we comprehensively investigated genome-wide chromatin changes associated with transcriptional reprogramming response to heat stress in tomato. Our data show that heat stress induces rapid changes in chromatin architecture, leading to the transient formation of promoter-enhancer contacts, likely driving the expression of heat-stress responsive genes. Furthermore, we demonstrate that chromatin spatial reorganization requires HSFA1a, a transcription factor (TF) essential for heat stress tolerance in tomato. In light of our findings, we propose that TFs play a key role in controlling dynamic transcriptional responses through 3D reconfiguration of promoter-enhancer contacts.

The linear conception of the genome has been replaced by an emerging understanding of the highly complex and dynamic 3D organization of chromatin within the nucleus¹⁻³. Chromatin is a compact, multiscale-organized, and dynamic structure: in eukaryotes, DNA is wrapped around histones octamers (H3, H4, H2A, H2B)⁴, and this structure folds again to form chromatin loops,

domains and territories, organizing the genome at multiple scales⁵. Chromosomes occupy defined territories inside the nucleus, as revealed by three-dimensional -fluorescence in situ hybridization (3D-FISH) methods and by high-throughput chromosome conformation capture (Hi-C) that show preferential interactions of each chromosome with itself⁵,⁶.

<figure><img src="image_6.png" /></figure>
In mammals, chromosomes are further organized into topologically associating domains (TADs), genomic regions that preferentially form contacts with each other rather than with surrounding sequences. TADs have been proposed to allow coordinated regulation of genes situated inside the TAD by allowing contact with distal regulatory sequences with in the TAD and limiting their access to sequences outside of the TAD; although this model is still under debate⁷. TAD associations form “chromatin compartments”. Active (A) compartments represent active chromatin defined by active histone marks, a high density of genes, and an elevated transcription rate supported by high RNA polymerase II (RNAPII) activity. In contrast, B compartments represent inactive chromatin defined by repressive histone marks, DNA methylation, and an increased number of transposable elements⁸,⁹. At a smaller scale, the formation of chromatin loops brings together distal regulatory sequences with promoter regions, or isolates them from promoters¹⁰. These short and long-range interactions modulate the expression of neighbouring and distant genes¹¹. Chromatin loop may involve self-interacting transcription factors (TFs), that bring together the genome regions they bind¹²,¹³.

All these organizational levels have been defined in animals and observed in plants: A and B compartments alternate along the genome, with A compartments predominantly in telomeric regions and B compartments in pericentromeric regions in various crop species⁹,¹⁴. However, plant chromatin organization shows several differences from that in animals. First, plant TAD-like structures represent only about 25% of the rice (*Oryza sativa*) genome¹⁵ compared to 75% for humans¹⁶. However, this low number of TADs in rice be a consequence of how we currently define TADs and may not reflect the characteristics of plant genomes¹⁷. Likewise in wheat (*Triticum aestivum*), TAD-like structures correspond to large heterochromatin domains highly enriched in transposable elements, rather than to groups of co-regulated genes³. The term “intergenic condensed spacers” (ICONS) has been proposed to describe these structures, as they seem to condense silent chromatin to allow loop formation between genes located outside of the ICONS, forming transcription factories. Second, mechanisms controlling chromatin folding differ between plants and animals. In animals, cohesin complexes extrude chromatin via their ring-like shape until the chromatin binds at two convergent CTCF proteins at specific CTCF binding sites that define TAD boundaries¹⁸. However, plant genomes apparently lack CTCF homologues, therefore how TAD-like domains are structured in plants is unclear.

Although the overall structure of chromatin architecture is very robust across organs and growth conditions, detailed analyses of chromatin folding revealed it can be modulated during development¹⁹, depending on cell type²⁰ or in response to changes in environmental conditions²¹. In plants, some changes in chromatin architecture in response to various stresses likely coordinate global transcriptome modifications for appropriate cellular and physiological responses. The importance of chromatin conformation for gene regulation in response to environmental stress in *Arabidopsis thaliana* has been well documented²¹,²². However, whether and how stress induces chromatin reorganization in crops is less understood, and only a few studies have examined the plant response to temperature changes. Rice chromosomes decondense under cold stress, decreasing long-range interactions above 1Mb and increasing A-B compartment interactions¹⁵. Conversely, exposure to heat stress reduces short-range interactions between A-B compartments, as A-A and B-B compartment interactions decrease¹⁵.

Gene expression is also controlled by establishment or dissociation of interactions between distal regulatory elements (REs) and promoters, such as enhancer-promoter interactions. In plants, such long-range interactions between promoters and distal regulatory sequences are most prominent in species with large genomes, such as maize (*Zea mays*)²³,²⁴. Modulating these interactions in response to

stress is likely functionally relevant, as has been shown in Drosophila²⁵ and human cells²⁶, but remains largely unexplored in plants.

Here, we studied the chromatin-based regulation of RE-promoter contacts and how nuclear DNA reorganization affects gene expression. To this end, we captured the temporal dynamics of these contacts using a time-course analysis in tomato (*Solanum lycopersicum*) during heat stress (HS). We find that HS induces profound changes in chromatin accessibility, as well as dynamic interactions between promoters and distal REs. We further show that the Heat Shock Factor (HSF) TF HSFA1a, a master transcriptional regulator of heat stress responses in plants²⁷, plays a key role in the dynamic formation of promoter-enhancer contacts in response to heat at several loci. Together, these results allow us to propose a model in which changes in chromatin accessibility and binding of HSFA1 promote the formation of promoter-enhancer contacts to induce gene expression in response to HS.

## Results

### Tomato chromatin is highly compartmentalized and organized around transcription factories

To explore tomato chromatin architecture, we performed an immuno-detection experiment on interphase somatic nuclei using antibodies directed against two different histone modifications, histone H3 lysine 9 acetylation (H3K9ac) and histone H3 lysine 27 monomethylation (H3K27me1), associated with transcriptionally active chromatin and constitutive heterochromatin, respectively. Immuno-staining revealed these marks were not homogeneously distributed in the nucleus, indicating distinct euchromatin and heterochromatin compartments in tomato: euchromatin occurred in the centre of the nucleus while constitutive heterochromatin occurred at the nuclear periphery (Fig. 1a, Supplementary Fig. 1).

To provide a high resolution view of chromatin architecture in tomato, we performed Hi-C, a genome-wide method to detect DNA-DNA physical interactions⁸. We observed a strong signal on the main diagonal, indicating frequent interactions between adjacent loci (Fig. 1b). The Hi-C map also revealed a strong compartmentalization segregating euchromatin and constitutive heterochromatin. To assess the ability of a genomic region to form a specific type of compartmental domain with specific epigenetic features, we generated a Pearson correlation map sorted it by multiple chromatin landscapes. To that end, we reordered the rows and columns of the correlation matrix: instead of arranging according to their position along the linear sequence, we arranged bins by increasing signal of the chosen feature. Notably, when sorted according to PC1 values (discriminating between A and B compartments), the matrix displayed clear segregation of chromosome 10 into two types of compartmental domains (Fig. 1c, d). These results highlight that histone marks strongly correlate with and are predictive of compartmental organization in tomato. In addition, integrating the Hi-C interaction matrix with RNAPII chromatin immunoprecipitation combined with sequencing (ChIP-seq) at different resolutions revealed interaction hotspots between genomic bins containing actively transcribed genes, indicating that tomato chromatin is organized around transcription factories (Supplementary Fig. 2). We performed immunostaining experiments using an anti-RNAPII antibody and observed RNAPII foci (Supplementary Fig. 1b). To further validate the presence of transcription factories, we performed a HiChIP experiment, a sensitive method to analyze protein-centric chromosome conformation, using an anti-RNAPII antibody³. These analyses confirmed that active chromatin in tomato is largely organized around transcription factories (Fig. 1e), as previously described in wheat³.

### Heat stress strongly affects global chromatin organization

To understand the impact that environmental stress has on tomato chromatin 3D organization and gene expression, we focused our analysis on HS, which affects both components in different plant species²²,²⁸. To this end, we exposed tomato plants to high heat for 1
<figure><img src="image_7.png" /><figcaption>**Fig. 1** | Tomato chromatin architecture displays a strong compartmentalization and is organized around transcription factories. a Immunofluorescence detection of H3K9ac (blue) and H3K27me1 (orange) histone modifications and DAPI staining (grey) in an isolated tomato nucleus. 3 times each experiment was repeated independently with similar results. b Chromosome 10 Hi-C interaction matrix at increasing levels of resolution. H3K9ac ChIP-seq signal (blue) and H3K27me1 signal (orange) were aligned with the maps to highlight the correlation with compartmentalization in tomato. c Pearson correlations of distance-</figcaption></figure>

normalized Hi-C interaction frequency maps of chromosome 10. The ChIP-seq signal for H3K9ac (blue) and H3K27me1 (orange) as well as the -PC1 component from principal component analysis (black) were aligned to the map. d Previous interaction matrix with reorganized bins according to their PC1 value.

e Chromosome 1 interaction matrix obtained by calculating the difference between the RNAPII HiChIP and the Hi-C signal. Red dots are interactions enriched in RNAPII HiChIP compared to Hi-C, and blue dots are interactions enriched in Hi-C compared to RNAPII HiChIP.
and 6 hours (Supplementary Fig. 3), these time points were selected based on the expression of the HS early response gene HSFA2 that peaked after 1 h of exposure and returned to basal level after 4 to 8 h (Supplementary Fig. 4). We next performed RNA-seq to obtain a global view of the changes in gene expression induced by HS. We found that compared to control conditions, one hour of HS induced the expression of genes known to be involved in HS responses (Supplementary Fig. 5a). To visualize dynamic gene expression patterns, we applied an unsupervised, self-organizing map (SOM) machine-learning approach²⁹ to reveal correlations between differentially expressed genes (DEGs) in multiple treatments (Supplementary Fig. 5b). Notably, this analysis revealed that a large part of DEGs after one hour of HS reverted their expression to their original expression levels after six hours (Supplementary Fig. 5c, d). To determine the extent to which chromatin architecture is affected by HS in tomato, we measured nuclei size after 0, 1, 6 h of HS and found that nuclei increased size after heat treatment (Supplementary Fig. 6).

To test if chromatin architecture reorganization and gene expression changes under heat stress are correlated, we performed Hi-C experiments using the same conditions employed for the transcriptomic analysis. To define the effects of HS on chromatin architecture, we measured the relative interaction differences at a 100-kb resolution (Fig. 2a). Visual inspection revealed that in plants exposed to HS, the interactions within constitutive heterochromatin were reduced whereas almost no changes were observed within euchromatin. To quantify these effects, we generated scaling plots of interaction frequencies against genomic distance at a 100-kb resolution on control and heat stressed plants. We found that HS caused a decrease of chromatin interactions in constitutive heterochromatin regions whereas chromatin interactions in euchromatin were enhanced (Fig. 2b, c, Supplementary Fig. 7). These observations suggest that heat negatively affects the strength of the compartmentation in constitutive heterochromatin, and positively affects it in euchromatin. Nevertheless, we cannot rule out the possibility that the slight increase in the number of interactions detected in euchromatin could be a mere consequence of the reduction in heterochromatin interactions, which would mechanically increase the sequencing depth for euchromatic interactions. To test this, we generated saddle plots to measure the degree to which A (euchromatin) and B (constitutive heterochromatin) compartments segregate in the nucleus. We found that HS positively affected the A-A interactions and negatively affected the B-B interactions (Fig. 2d, Supplementary Fig. 8). However, we integrated RNA-seq and Hi-C data at a 5 kb resolution across the HS time course and found that HS DEGs resided in domains in which the resolution of our Hi-C does not allow to observe major changes (Fig. 2e), suggesting that regulatory dynamics may occur via intradomain contacts. To test this hypothesis, we used ATAC-seq which identifies accessible chromatin sites to identify putative REs that may physically interact with promoters of HS DEGs (Fig. 2f-k, Supplementary Fig. 9). To visualize the dynamics of chromatin accessibility patterns across the entire HS time course, we applied an unsupervised, SOM machine-learning approach²⁹ (Fig. 2f). The chromatin accessibility pattern was dynamically modified by heat stress and exhibited temporal changes in RE activation. Since these regulatory sequences can be hundreds to thousands of base pairs upstream or downstream of the gene they regulate, a simple association between chromatin accessibility and gene expression may not be always accurate. Therefore, we thus classified all the accessible chromatin peaks into two categories: those within 1.0 kb upstream of transcription start sites (TSS-proximal) and those found further away from genes (TSS-distal). We focused our analysis on TSS-proximal peaks. Using a clustering approach based on the level of chromatin accessibility, we identified five clusters (I to V) (Fig. 2g). Chromatin accessibility is necessary, but not sufficient, to induce gene expression, which is also determined by histone modifications and the availability of recruited TFs³⁰. We found that the accessibility of proximal REs correlated

positively with gene expression (Fig. 2h, i). The ATAC-seq time-course analysis can uncover sequential TF binding based on characteristic chromatin footprints. We therefore used the accessible chromatin data from the HS time course to determine the hierarchical gene regulation networks modulated by heat (Fig. 2i). We found that Heat Stress Transcription Factors (HSF) and other TF families (ERF, Myb-related and AP2) were upregulated only shortly after HS (Fig. 2j) while that TFs from the C2H2, Homeobox and TCP families were downregulated after 1 hour of HS. Moreover, we generated hierarchical TF networks by taking into account changes in gene expression as well as the presence of specific TF binding sites in their promoter region. We found that the chromatin accessibility at specific TF-binding motifs, which reflects the TFs activity, coincided with altered expression of genes harbouring an accessible binding site for this TF in their promoter (Fig. 2k). Altogether, our analyses suggest that sequential activation of TFs controls gene expression during the heat stress response via increased chromatin accessibility: a first wave of transcriptional changes could be controlled by HSF and ERF factors, and would involve early response event, and a second wave would subsequently occur under the control of WRKYs and TCPs.

## Distal and proximal REs display different chromatin signatures

Active REs are associated with specific histone modifications in mammals (i.e. H3K4me1 and H3K27ac)³¹,³². To determine if a typical chromatin signature differentiates distal and proximal REs in tomato, we performed ChIP-seq under control and HS conditions and analyzed the levels of four histone covalent modifications that have been previously associated with active transcription: H3K4me3, H3K9ac, H3K18ac and H3K27ac. We found that the accumulation of histones with these modifications did not vary significantly under HS (Supplementary Fig. 10). However, we found that nucleosomes flanking accessible proximal REs displayed a high level of H3K9ac, H3K18ac, H3K27ac and H3K4me3 (Fig. 3a, b, Supplementary Figs. 11–13). In contrast, nucleosomes flanking accessible distal REs displayed a high level of H3K9ac and H3K18ac, a moderate level of H3K27ac and low level H3K4me3 (Fig. 3a, b). These data, in combination with chromatin accessibility data, suggested that nucleosomes decorated with H3K9ac and H3K18ac and devoid of H3K4me3 may indicate active enhancer positions in the genome, and could be used to annotate enhancer regions.

In mammals, it has been reported that distal enhancers recruit RNAPII³³. We postulated that this feature could also be present at distal active enhancers in plants. To test this hypothesis, we performed RNAPII ChIP-seq on tomato plants subjected to HS and analyzed its accumulation on distal REs specifically active at 1 hour after HS. We identified a large number of RNAPII binding sites near TSSs, but also present at the distal REs identified by our chromatin analyses (Fig. 3c, Supplementary Figs. 14, 15). We then determined motifs associated with these HS-activated distal REs and found that both proximal and distal elements contained motifs for the same TFs, suggesting that the same TFs could bind both proximal and distal REs (Fig. 3d).

## 3D chromatin reorganization transiently induces promoter-RE interactions in response to HS

To determine whether distal REs interact with promoters, we further combined Hi-C with a capture enrichment step (C-Hi-C), a powerful tool to delineate spatial and functional relationships at specific chromatin regions. To this end, we generated a biotinylated RNA bait library specifically targeting 212,735 promoter regions (Supplementary Data 1) and performed a C-Hi-C across the HS time course to reveal the different chromatin associations during HS. As expected, C-Hi-C enables the identification of a large number of chromatin contacts, increasing the resolution of our analysis (Supplementary Fig. 16). Interestingly, we found that in addition to the classical 1 vs. 1 RE-promoter contacts, we observed two types of interactions implicating multiple loci: (i) one promoter interacted with several REs, called
<figure><img src="image_8.png" /></figure>

"promoter-centric hubs", and (ii) one RE interacted with several genes, named "RE-centric hubs" (Fig. 3e). To further dissect the properties of RE-promoter contacts, we carried out a differential analysis integrating interaction data from control and HS treatments. This analysis revealed notable dynamic changes between the different types of RE-promoter contacts (Fig. 3e, Supplementary Fig. 17, Supplementary Data 2). Distal REs are thought to be brought into physical proximity to

their target promoters via the three-dimensional looping of chromatin mediated by structural proteins such as TFs³⁴,³⁵. Therefore, we hypothesized that the strength of a RE-promoter chromatin interaction could be linked to the accessibility of both the RE and the promoter to TFs. To test this, we classified the promoter-associated interactions into three categories: (i) interactions for which neither anchor is not accessible, (ii) interactions for which only one anchor is not accessible,
**Fig. 2** | The high-order tomato chromatin organization is strongly affected by heat stress. a Chromosome 10 Hi-C interaction frequency heatmap at a 100-kb resolution in control and heat stress conditions (1 h). Each pixel denotes all interactions between any to 100 kb genomic loci from the linear genome. Intensity represents log₂ normalized contact frequencies. b Averaged scaling plot of interaction frequencies against increasing genomic distance for all tomato pericentromeric regions (gene poor region). The genomic bin size is 100 kb. c Averaged scaling plot of interaction frequencies against increasing genomic distance for all tomato chromosome arms (gene-rich regions). The genomic bin size is 100 kb. d Saddle plots of Hi-C data representing the mean observed interaction frequencies divided by expected interaction frequencies between 20 kb bins. Interactions between A compartments are in the top right, and interactions between B compartments are in the bottom left.

e Hi-C interaction heat maps at a 5 kb resolution surrounding a genomic region with differential gene expression along the time course analysis. f Self-organizing map (SOM) of chromatin accessibility profiles at 0, 1 and 6 h of heat stress. Red areas mark over-accessible chromatin sites, blue areas mark under-accessible chromatin sites. g K-means clustering of differentially accessible peaks. h Example of proximal regulatory elements (REs) showing correlation between RNA-seq signals and ATAC-seq signals along the time course. i Correlation between ATAC-seq and RNA-seq data shown by heatmaps of proximal accessible chromatin sites and corresponding gene expression levels. j HOMER DNA-motif enrichment analyses of proximal accessible chromatin peaks. k TF network built by integration of RNA-seq and ATAC-seq data. Source data of Fig. 2g are provided as a Source Data file.

and (iii) interactions for which both anchors are accessible. The strength of interactions in which both anchors are accessible was stronger than the interactions in which only one anchor or no anchors are accessible (Supplementary Fig. 18). Thus, HS induces transient promoter/RE interactions that correlate with an increase in chromatin accessibility at both sites. Together with the observation that both proximal and distal REs are enriched in the same TF binding motifs points these results point to a putative role of DNA binding proteins such as TFs to promote or stabilize these interactions.

Our integrative analyses revealed that the HSF binding motif was over-represented among the REs involved in the promoter-enhancer interactions only after 1 h of HS treatment. The first exhaustive overview of HSF family in plants was described in *Arabidopsis*, where 21 representatives were identified, belonging to three classes (A, B, and C) according to the structural features of their oligomerization domains³⁶. Class A HSFs are essential for transcriptional activation, whereas class B and C HSFs have no activator function because they lack the acidic amino acid residue containing motifs³⁷. Previous studies have shown that HSFA1a is the master regulator of the transcriptional HS primary response²⁷,³⁸,³⁹. Interestingly, previous work in *S. lycopersicum* showed that HSFA1a is both a master regulator of the HS response and basal thermo-tolerance³⁸. We thus hypothesized that HSFA1a may be implicated in RE-Promoter interactions during early HS responses. To test this, we first determined HSFA1a binding sites using a TF-DNA binding assay called DNA affinity purification and sequencing (DAP-seq)⁴⁰. This approach revealed a large number (5034) of putative HSFA1a binding sites (Fig. 4a, Supplementary Fig. 19). As expected, these binding sites were preferentially found in the immediate vicinity of the TSS of annotated genes (Supplementary Fig. 20) which were associated with 1,035 HS DEGs (Fig. 4b, Supplementary Fig. 21). Interestingly, a significant part of these potential binding sites was found in chromatin regions that became transiently accessible 1 h after HS (Supplementary Fig. 22). Implementing a de novo motif discovery analysis we found that the HSFA consensus sequence (TCTAGAANNTTCT) was over-represented within the HSFA1a DAP-seq peaks (Fig. 4c). Integrating DAP-seq and transcriptomic data revealed that 65% of the putative HSFA1a target genes that were differentially expressed after 1 h of HS were up-regulated whereas 35% were down-regulated supporting the view that HSFA1a acts primarily as a transcriptional activator (Fig. 4b). Next, we explored the relationship between chromatin accessibility and HSFA1a-enriched loci by integrating ATAC-seq and DAP-seq data (Fig. 4d). A significant proportion of chromatin regions accessible after 1 h of HS displayed HSFA1 binding, whereas this was not the case in chromatin regions that were accessible in control conditions. We then tested whether the strength of spatial contacts between HSFA1a binding sites and distal REs differed in heat stress vs. control conditions. To test this idea, we generated aggregate plots (aggregated C-Hi-C matrices) to quantify the mean of aggregated/stacked C-Hi-C submatrices between HSFA1a-distal RE loops after an observed/expected transformation of the C-Hi-C matrix, for control and HS conditions. These analyses revealed an increase of the strength of the HSFA1a-distal RE contacts after 1 h of HS compared to the control (Fig. 4e, left). Since we

previously observed that many distal REs were bound by RNAPII, we performed a HiChIP experiment with an anti-RNAPII antibody to detect and measure the strength of the HSFA1a-distal RE interactions. This analysis revealed an increase of the number and strength of the distal RE-HSFA1a contacts after 1 h of HS compared to the control (Fig. 4e, f). Taken together, our data suggest that HSFA1a plays a major role in RE-promoter interactions.

### HSFA1a plays a major role in the dynamic formation of promoter-enhancer contacts in response to HS

To determine if the dynamic formation of chromatin loops induced by HS is dependent on HSFA1a activity, we undertook a genetic approach. We used two independent *hsfa1a* knock-down transgenic lines (*hsfa1a-line1* and *hsfa1a-line2*) carrying a hairpin RNA construct that down-regulates *HSFA1a* expression²⁷ (Supplementary Fig. 23). In accordance with previous studies, we found that compared to wild type (WT) plants, *hsfa1a* knock-down transgenic plants displayed strong defects in growth when exposed to HS (Fig. 5a). To determine the role of HSFA1a on HS-mediated gene expression, we performed a transcriptomic analysis on control and *hsfa1a* knock-down transgenic lines subjected to 1 h of HS (Fig. 5b). Gene ontology analysis on down-regulated genes in *hsfa1a* knock-down transgenic lines compared to WT plants revealed significant enrichment in genes involved in the HS response (Fig. 5c). We found that 247 of the HSFA1a targets identified by DAP-seq were misregulated in *hsfa1a* knock-down transgenic lines, 186 (75.3%) being downregulated and 61 (24.7%) upregulated (Fig. 5d and Supplementary Fig. 24). This data reinforces the view that HSFA1a acts as a positively regulator of gene expression in response to HS.

To assess whether HSFA1a influences RE-promoter loop formation, we focused our analysis on *Solyc09g074475* and *HSFA2*. These two genes were selected because they are both heat-responsive in an HSFA1-dependent manner, they display an HSFA1a binding site in their promoter region, and dynamic RE-promoter contacts in response to HS. Furthermore, HSFA2 is well-known as a key component of thermotolerance in *Arabidopsis*⁴¹. C-Hi-C and HiChIP-RNAPII revealed the presence of a discrete RE that interacts with *Solyc09g074475* and an unannotated gene, and that this interaction was enhanced in response to HS (Fig. 5e, Supplementary Fig. 25a). Notably, the formation of these two chromatin loops correlated with a strong activation of both genes in response to HS, thus suggesting that this element may act as a transcriptional enhancer. To test whether this putative RE has the potential to function as an enhancer, we generated a construct containing this RE sequence fused to a minimal mosaic virus (CaMV) 35 S promoter to drive the expression of a luciferase (LUC) reporter after transformation of tobacco *N. benthamiana* leaves. We found that compared to the minimal promoter, the addition of the *Solyc09g074475* RE sequence significantly increased the expression of the LUC reporter (*P* = 8.40e-05, Student's *t*-test) (Fig. 5e). These results supported the view that this genomic region can act as a transcriptional enhancer. To assess whether HSFA1a is necessary for the formation of chromatin loops with this regulatory sequence, we used chromosome conformation capture (3C) followed by qPCR (3C-qPCR)
<figure><img src="image_9.png" /><figcaption>**Fig. 3** | Distal and proximal REs displayed different chromatin signatures.</figcaption></figure>

a Heatmaps of ATAC-seq and ChIP-seq (H3K4me3, H3K9ac, H3K18ac and H3K27ac) profiles over proximal (grey) and distal (purple) regions at 0, 1 and 6 h of heat stress. b Example of a distal regulatory element (RE) and of the histone modifications surrounding it after 0, 1 and 6 h of heat stress. c Heatmaps of ATAC-seq signal and RNAPII occupancy over distal REs at 0, 1 and 6 h after heat stress. d K-means

clustering of distal differentially accessible chromatin regions and the corresponding top enriched motif identified by HOMER (v4.11) for each cluster. Binomial distribution test was used in detected ATACseq peaks. e Example of two types of transient chromatin loop formation: a promoter-centric hub and a RE-centric hub. Significant interactions for C-HiC and HiChIP-RNAPII data were detected with HOMER (v4.11) using cumulative binomial distribution with p value 0.05.
a

<figure><img src="image_10.png" /></figure>

b

<figure><img src="image_11.png" /></figure>

c

<figure><img src="image_12.png" /></figure>

P-value = 1e-12.017

d

<figure><img src="image_13.png" /></figure>

Proximal

<figure><img src="image_14.png" /></figure>

Distal

<figure><img src="image_15.png" /></figure>

e

<figure><img src="image_16.png" /></figure>

f

<figure><img src="image_17.png" /></figure>

on plants grown under control conditions or subjected to 1 h of HS. We found that in WT plants, the strength of the two loops, named A and B, increased after 1 h of HS, whereas in *hsfa1a* knock-down transgenic lines the strength of the two loops did not increase but instead was significantly reduced (Loop A: *PI* = 3.20e-05, *P2* = 3.20e-05 and Loop B: *PI* = 8.15e-05, *P2* = 4.23e-05, Student's *t*-test) (Fig. 5f). With respect to the *HSFA2* locus, C-Hi-C and HiChIP-RNAPII revealed that it is potentially regulated by a chromatin loop with two distal REs (named distal

RE-A and distal RE-B) (Fig. 5g, Supplementary Fig. 25b). The formation of this chromatin loop correlated with the activation of *HSFA2*, thus suggesting that these elements could act as transcriptional enhancers. To test this, we followed the previously described approach and found that compared to the control minimal promoter construct, the signal strength of either distal REs was significantly increased (RE-A, *P* = 3.21e-04; RE-B, *P* = 1.87e-04) (Fig. 5g). We used 3C-qPCR to determine whether HSFA1a influences this distal RE-promoter loop formation, and
**Fig. 4** | HSFA1a binding sites correlate with anchors of transient promoter-enhancer chromatin loops formed in response to heat stress. a Heatmaps of control Halo-tag and HSFA1a density in the region of ±1 kb around the HSFA1a peaks. b Proportion of HSFA1a targets that were mis-regulated after 1 h of heat stress compared to the control condition; 65% of HSFA1a targets were upregulated after 1 h of heat stress. c Top enriched motif in HSFA1a peaks identified by HOMER (v4.11). Binomial distribution test was used in detected DAPseq peaks. d Heatmaps of DAP-seq and ATAC-seq profiles surrounding proximal (top) or distal (bottom) regions accessible only in the control condition (left) or accessible only after 1 h of heat stress (right). e Aggregate plots quantifying the mean of aggregated/stacked

C-Hi-C submatrices (left) or RNAPII HiChIP (right) between HSFA1a and distal regulatory element (RE) loops after an observed/expected transformation of the matrix, for control and heat stress conditions. f A connection between gained contacts and induction of heat-related genes is exemplified. DAP-seq Halo-tag signal is represented in black and HSFA1a in red. RNA-seq signals are represented by orange peaks, ATAC-seq signals by blue and ChIP-seq signal of RNAPII by purple peaks. Chromatin interactions identified by C-Hi-C are represented by red lines and HiChIP RNAPII by purple lines. Significant interactions for C-Hi-C and HiChIP-RNPII data were detected with HOMER (v4.11) using cumulative binomial distribution with *p*-value 0.05.

found that depleting HSFA1a decreased the strength of HS-specific contacts between the distal putative enhancer and the HSFA2 locus (P1 = 3.26e-08, P2 = 3.26e-08, Student's *t*-test) (Fig. 5h). Similar experiments were conducted with two additional candidate enhancers interacting with genes involved in heat stress response, leading to the same results (Supplementary Fig. 26). Collectively, our results suggest that HSFA1a plays a major role in the dynamic formation and/or stability of promoter-enhancer interactions that underpin the initial transcriptional responses to HS (Fig. 6).

## Discussion

This study presents a deep analysis of the tomato nuclear architecture by integrating capture-Hi-C, ATAC-seq, ChIP-seq, HiChIP, and RNA-seq data. Results obtained highlight the specific functionalities of the different levels of chromatin spatial organization underscoring the regulatory importance of nuclear organization as a whole. Chromatin conformation in tomato is strongly compartmentalized and organized around focalized sites of transcription in the active chromatin that resemble transcription factories$^{3,42}$. Although the precise mechanisms underlying transcription factories organization remain poorly understood, they are thought to play a critical role in clustering co-regulated genes$^{43}$. The presence of transcription factories in the tomato nucleus reveals their importance in transcriptional regulation and in particular on the optimization of transcription in response to environmental cues by concentrating TFs in discrete nuclear compartments.

Temporal and spatial gene expression is influenced by dynamics chromatin changes and is implicated in the control of developmental programs as well as in response to environmental stresses$^{14,21}$. In metazoan interphase nuclei, chromatin is organized into TADs, which are thought to bring together co-regulated genes. In Drosophila, HS induces dramatic rearrangement of the 3D nuclear organization$^{44}$, suggesting that TAD organization of metazoan genomes is plastic and is reconfigured in response to environmental changes$^{25}$. However, a recent study in human and Drosophila cells has shown that HS only affects a subset of enhancer-promoter interactions whereas TADs and compartment structures remain mostly stable upon HS$^{45}$. In plants, HS induces high-order reorganization of the constitutive heterochromatin compartment$^{22,28}$. In line with these results, we observed that HS in tomato negatively affects the chromatin interactions between constitutive heterochromatin regions, thus affecting the strength of constitutive heterochromatin compartmentalization.

REs are non-coding DNA sequences containing binding motifs for one or several TFs involved in transcriptional regulation of target genes, and are also key components of the 3D chromatin organization. These sequences act as proximal or distal RE and their activity is reflected by their accessibility, which is likely controlled by the action of pioneer TFs and chromatin modifiers$^{46,47}$. In animals, distal REs are associated with the specific histone marks H3K27ac and H3K4me1, defining a specific chromatin state$^{48}$. By contrast, in plants, studies in maize and wheat have proposed that H3K27ac, H3K4me3 and H3K9ac are histone marks preferentially associated with distal REs$^{49–51}$. Here, we focused on differentially accessible REs as a proxy for active REs and found that active distal and proximal REs displayed different chromatin signatures. We have found that in tomato distal REs are

characterized by the presence of nucleosomes rich in H3K9ac and H3K18ac and devoid of H3K4me3, flaking the region with maximal accessibility. Importantly, we demonstrated that some of these RE are capable to act as transcriptional enhancers.

3D chromatin architecture dynamically controls the access of REs to their target genes by promoting or inhibiting RE-promoter interactions. RE sequences can be hundreds to thousands of base pairs upstream or downstream from the gene they regulate, making it difficult to define which REs control which genes. The rapid development sequencing technologies has enabled multiple high-throughput experimental methods to measure chromatin interactions, including Capture-Hi-C and HiChIP, which are particularly useful to detect enhancer-promoter interactions. In animals, promoter C-Hi-C has been used for the analysis of long-range RE-promoter interactions and their control of gene expression using different cell types. These analyses have revealed a massive change in RE-promoter interactions taking place during cellular differentiation$^{52–56}$. In plants, RE-promoter interactions were analyzed in maize using HiChIP for H3K4me3, H3K27me3 and H3K27ac, which identified long range RE-promoter interactions that contribute to the control of gene expression$^{57,58}$. This analysis uncovered the existence of long-range RE-promoter interactions. In line with this observation, using both capture-Hi-C and HiChIP approaches in tomato, several thousands of RE-promoter interactions were discovered. Notably, we have found examples for a promoter interacting with several REs as well as examples where one RE interacts with multiple genes. These two types of chromatin could be defined as promoter-centric hubs or RE-centric hubs respectively$^{59}$.

One of the mechanisms proposed to shape the 3D genome structure in animals is the loop extrusion model, that involves the structural maintenance of chromosomes (SMC) cohesin complex and the zinc-finger DNA-binding protein CTCF$^{60,61}$. Plants genomes appear to lack CTCF homologues, but SMC complexes are highly conserved. Therefore, this loop extrusion may be conserved in plants and involve yet unknown factors that act as insulators. In addition to the action of SMC and CTCF, TFs serve as protein anchors and determine 3D genome organization$^{59}$. TFs may mediate RE-promoter interactions through different mechanisms, including direct protein homodimerization$^{62}$ and co-factor protein recruitment to generate oligomers$^{35}$. In *Arabidopsis*, the CONSTANS TF plays a crucial role in forming of an enhancer-promoter loop to control the expression of *FLOWERING LOCUS T*$^{63}$, thus suggesting that TFs also impact 3D chromatin organization in plants. However, most of studies in plants have focused on deciphering the role of TFs in 3D chromatin architecture and gene regulation under steady-state conditions, thus providing clues on the mechanisms involved in the maintenance of RE-promoter interaction rather than in their formation. To decipher the molecular mechanisms involved in establishing distal RE-promoter interactions, we carried out a multidimensional study combining stress at different time points. Our data show that 3D chromatin interactions are dynamically modified by HS to allow a temporally evolving connection between a diversity of activated REs and promoters. Furthermore, our integrative analyses revealed that the binding motif of HSF was over-represented in the sequences involved in RE-promoter contacts at the onset of HS treatment. These data suggest that
<figure><img src="image_18.png" /><figcaption>b</figcaption></figure>

<figure><img src="image_19.png" /><figcaption>c</figcaption></figure>

<figure><img src="image_20.png" /><figcaption>d</figcaption></figure>

<figure><img src="image_21.png" /><figcaption>e</figcaption></figure>

<figure><img src="image_22.png" /><figcaption>f</figcaption></figure>

<figure><img src="image_23.png" /><figcaption>g</figcaption></figure>

<figure><img src="image_24.png" /><figcaption>h</figcaption></figure>

<figure><img src="image_25.png" /></figure>

transcription factors of the HFS family play a major role in RE-promoter contact formation. Studies in yeast have shown that HSF1 mediates the clustering of its targets and dynamically influences the reconfiguration of the genome³⁵,⁶⁴. Our data for tomato indicates that HSFA1a plays a critical role in the dynamic formation of promoter-enhancer contacts and in controlling the transcriptional response at the onset of HS. Future work shall help determine whether loop

formation is a pre-requisite for transcription activation, or is rather a consequence of it. Few studies have assessed the dynamic role of TFs in the formation of chromatin loops in response to a stimulus, even in animals. One study revealed that KLF4 and ZNF750 TFs, rather than cohesins, are key to establishing of RE-promoter contacts to promote gene expression during epidermal cell differentiation³⁵. Collectively, our findings shed light on the complexity of distal and proximal RE
**Fig. 5** | HSFA1a plays a major role in the dynamic formation of promoter-enhancer contacts in response to a heat stress. a Phenotypes of WT and *hsfa1a* knock-down (via RNA interference) lines after 0 and 6 h of heat stress. b Expression level heatmap of genes that were down-regulated in two independent *hsfa1a* lines after 1 h of heat stress compared to the WT. c Gene ontology enrichment of down-regulated genes (n = 280). P-values calculated using Fisher's exact test. d Pie chart representing the 247 HSFA1a target genes that were mis-regulated in *hsfa1a* lines after 1 h of heat stress compared to the WT. Among them, 75.3% were down-regulated in *hsfa1a* compared to the WT. e, g Screenshot of the *Solyc09g074475* and *HSFA2* locus. DAP-seq against Halo-tag and HSFA1 are represented in black and red. The RNA-seq, ATAC-seq and RNAPII ChIP-seq signal is represented in orange, blue and purple respectively. Chromatin interactions from C-Hi-C data and RNAPII HiChIP are represented by red and purple lines. The right green box of e shows a zoom of a distal RE that interacts with two heat stress-responsive genes (*Solyc09g074475* and

one non-annotated gene). The purple box of **g** shows a zoom of two distal REs named A and B that interact with **HSFA2**. The histogram represents the luciferase activity from a transient luciferase assay using a **RE-mini35S::LUC** construct or the **mini35S::LUC** alone as a control. The mean LUC/REN activity levels were normalized to **mini35S::LUC** as control (n = 3 biological replicates), **f, h** (Top) Design of the 3C-qPCR assays used to analyze the distal RE loops in the WT and **hsfa1a** lines in response to heat stress. (Bottom) Quantitative 3C; the relative interaction frequencies were calculated as described in the Materials and Methods. Data are the average of three biological replicates each with three technical replicates (n = 3). Circles in **e, g** denote relative LUC/REN activity values. Circles in **e, f, g, h** denote relative interaction frequency. Bars indicate mean values +/- SD of three replicates. Exact *P*-values are shown (two-tailed, two-sample Student's *t*-test). Source data are provided as a Source Data file.

<figure><img src="image_26.png" /></figure>

**Fig. 6** | Model recapitulating the role of HSFA1a in enhancer-promoter loop formation to regulate gene expression in response to heat stress in tomato. The diagram was created using BioRender.
contact loop formation and provides evidence for a key role of TFs in
controlling transcriptional response through the 3D reconfiguration of
the chromatin structure in animals and plants.

Methods

Plant material and growth conditions

Tomato (Solanum lycopersicum) lines used in this study were in the cv. M82 and cv. Moneymaker backgrounds. The *hsfa1a* mutants line1 and line2 were obtained from Dr.Fragkostefanakis's lab²⁷. Seeds were directly sown on soil and plants grown in growth chambers at 24 °C under long-day (16 h light) conditions. For heat treatment, 4-week-old plants were treated at 45 °C for 1 h and 6 h in a climatic chamber (Aralab). The fourth leaf was used for all experiments in this research.

Immunostaining

Immunostaining experiments were performed according to Latrasse et al⁶⁵. Briefly, fourth leaves of 4-week-old tomato plants were fixed and then nuclei were isolated, placed on a poly-lysine slide, and incubated overnight at 4 °C with primary antibodies (400X diluted) of H3K9me2 (Abcam, ab1220), H3K9ac (Millipore, 07-352), H3K27me1 (Millipore, 07-448) and RNAPII (Abcam, ab26721). Slides were washed and then incubated for 1 h at room temperature in the dark with Goat anti-Rabbit Alexa Fluor Plus 488 (A11034 Invitrogen) and Goat anti-Mouse Alexa Fluor Plus 555 (A32727 Invitrogen) or with Goat anti-Rabbit Alexa Fluor Plus 555 (A32732 Invitrogen) and Goat anti-Mouse Alexa Fluor Plus 488 (A32723 Invitrogen) secondary antibodies (400X diluted). DNA was counterstained with 4,6 diamidino-2-phenylindole (DAPI) in Slow Fade Diamond Antifade mounting media. Slides were directly imaged on a confocal microscope (Zeiss Microsystems).

Assay for transposase-accessible chromatin with high-
throughput sequencing

Approximately 100 mg of 4-weeks-old tomato fourth leaves were ground, passed through a 100 µm filter and the nuclei were isolated using nuclei isolation buffer (0.25 M sucrose, 10 mM Tris-HCl pH8, 10 mM MgCl2, 1% Triton X-100, 5 mM 2-mercapto-ethanol, and 0.1 mM protease inhibitors). The nuclei were re-suspended in 2x TD buffer (20 mM Tris-HCl, pH 7.6, 10 mM MgCl2, and 20% dimethyl formamide) and 2.5 µl of Tn5 transposase (Illumina FC-121-1030). The transposition reaction was performed at 37 °C for 30 min, and DNA was purified using a Qiagen MinElute Kit (QIAGEN, Cat.No.28004). DNA libraries were amplified for a total of 10 cycles as described by Buenrostro et al⁶⁶ and Jegu et al⁶⁷. DNA libraries were checked for quality and quantified using a 2100 Bioanalyzer (Agilent) and subjected to 2 × 75 bp high-throughput sequencing by NextSeq 500 (Illumina). Two independent biological replicates were generated.

RNA-seq assay

Total RNA was extracted from the fourth leaves of 4-week-old tomato with Nucleospin RNA kit (Macherey-Nagel), according to the manufacturer's instructions. RNA-seq libraries were prepared from 2 µg of total RNA using NEBNext Ultra II Directional RNA library Preparation Kit (NEB) according to the manufacturer's instructions. RNA libraries were checked for quality and quantified using a 2100 Bioanalyzer (Agilent) and subjected to 1 × 75 bp high-throughput sequencing by NextSeq 500 (Illumina). Two independent biological replicates were generated.

In situ Hi-C assay

In situ Hi-C experiments were performed according to Concia et al.³ using DpnII enzyme (New England Biolabs). DNA libraries were prepared using NEBNext Ultra II DNA library preparation kit (NEB) according to the manufacturer's instructions (10 cycles for the PCR amplification step). DNA libraries were checked for quality and quantified using a 2100 Bioanalyzer (Agilent) and the libraries were

subjected to 2 × 75 bp high-throughput sequencing by NextSeq 500 (Illumina). Two independent biological replicates were generated.

C-Hi-C assay

For C-Hi-C, the in situ Hi-C libraries were used for the capture step. Capture was performed by using SureSelect XT Target Enrichment System for Illumina Paired-End Multiplexed Sequencing Library (Agi- lent) according to the manufacturer's recommendations. All tomato cultivar M82 promoters were selected1.5 kb in length and 212,735 probes (length of each probe = 120 bp) were designed to catch all promoters. The capture step was performed using 1 µg of in situ Hi-C libraries following the manufacturer's recommendations. The quality of the libraries was assessed using a 2100 Bioanalyzer (Agilent), and the libraries were subjected to 2 × 75 bp paired-end high-throughput sequencing by NextSeq 500 (Illumina). Two independent biological replicates were generated

HiChIP assay

Nuclei were isolated from the fourth leaves of 4-week-old tomato plants using the same procedure as for the in situ Hi-C experiments. The Hi-ChIP protocol from Mumbach et al.⁶⁸ was then applied using the DpnII restriction enzyme (NEB) and 3ug of anti-RNAPII antibody (Abcam, ab26721). The quality of the libraries was assessed using a Bioanalyzer (Agilent), and the libraries were subjected to 2 × 75 bp paired-end high-throughput sequencing by NextSeq 500 (Illumina).

ChIP-seq assay

ChIP-seq assays were performed on the fourth leaves of 4-week-old tomato plants according to Bio-protocol of Ramirez-Prado et al.⁶⁹ using 3ug of anti-H3K27me1 (Millipore, 07-448), anti-H3K9ac (Millipore, 07-352), anti-H3K14ac (Millipore, 07-353), anti-H3K18ac (Millipore, 07-354), anti-H3K27ac (Abcam, ab4729), anti-H3K4me3 (Millipore, 07-473), and anti-RNAPII (Abcam, ab26721) antibodies. ChIP-seq libraries were prepared from 10 ng of DNA using NEBNext Ultra II DNA Library Prep Kit for Illumina (NEB) according to the manufacturer's instructions. Two independent biological replicates were generated for each time point of heat stress. DNA libraries were checked for quality and quantified using a 2100 Bioanalyzer (Agilent) and subjected to 1 × 75 bp high-throughput sequencing by NextSeq 500 (Illumina).

Amplified DNA affinity purification sequencing assay (Amp-DAP-seq)

AmpDAP-seq assays were performed according to Bartlett et al.⁴⁰ with modifications. Briefly, gDNA was extracted from 4-week-old tomato fourth leaves, fragmented, and a library was generated using NEBNext Ultra II DNA library preparation kit (NEB, #E7645S/L, #E7103S/L). *SolycHsfa1a* open reading frames were transferred into the Gateway-compatible pIX-HALO expression vector and expressed using Promega TNT Coupled Wheat Germ Extract Systems (Promega, cat. no. L4130). The expressed proteins were immobilized on Magne HALO-Tag beads (Promega, G728A), washed, and incubated with the DNA library. After bead washing, DNA was eluted and amplified with indexed NEBNext primers (20 PCR cycles were used). DNA libraries were checked for quality and quantified using a 2100 Bioanalyzer (Agilent) and subjected to 1 × 75 bp high-throughput sequencing by NextSeq 500 (Illumina). Two independent biological replicates were generated.

3C-qPCR

Four-week-old tomato fourth leaves were cross-linked in 1% (v/v) formaldehyde at room temperature for 20 min. Cross-linked plant material was ground and nuclei were isolated and treated with 0.5% SDS at 62 °C for 5 min. The SDS was then supplemented with 2% Triton X-100. Digestions were performed overnight at 37 °C using 150U of DpnII enzyme (New England Biolabs). Restriction enzymes
were inactivated by adding 1.6% SDS and incubated at 62°C for 20 min. SDS was supplemented with 1% Triton X-100. DNA was ligated by incubating at 22°C for 5 h using 4000U of T4 DNA ligase (Fermentas). Reverse crosslinking was performed by overnight treatment at 65°C. DNA was recovered after Proteinase K treatment by phenol-chloroform extraction and ethanol precipitation. Relative interaction frequencies were calculated by quantitative PCR using 15 ng of DNA. A region uncut by DpnII was used to normalize the amount of DNA. Details for primers used for 3C-qPCR are listed in Supplementary Data 3. Two independent biological replicates were generated.

## Dual-luciferase assay

Luciferase activity was quantified using reagents from the Dual Luciferase Reporter Assay System (Promega, E2920) according to the manufacturer's instructions.

## Analysis of ChIP-seq data

Trimmomatic-0.38⁷⁰ was used for trimming with the following parameters: minimum length of 36 bp; mean Phred quality score >30; leading and trailing bases removal with base quality <5. The reads were mapped onto the M82 V1.0⁷¹ assembly using Bowtie2 V2.3.5⁷² with mismatch permission of 1 bp (Supplementary Data 4). To identify significantly enriched regions, we used MACS2 V2.2.7.1⁷³ with the following peak-calling parameters: number of duplicate reads at a location: 1; mfold of 5:50; q-value cutoff: 0.05; extsize 200; broad peak. To extract the average scores across the genomic regions, multiBigwigSummary command from the deepTools package⁷⁴ was used with default parameters on the RPGC normalized bigWig files.

## Differential expression analysis

Single-end reads of RNA-seq samples were trimmed using Trimmomatic-0.38 with the parameters: minimum length of 30 bp; mean Phred quality score >30; leading and trailing base removal with base quality <5. STAR aligner⁷⁵ was used to map the reads to the M82 V1.0⁷¹ genome assembly (Supplementary Data 4). Raw read counts were then extracted using the featureCounts V2.0.0 utility from the R Subread package based on the M82_v1.1.0 gene annotations⁷¹. Finally, we used DESeq2 V1.38.0⁷⁶ to identify differentially expressed genes. Genes with read counts ≥ 50 in at least 2 samples were considered for differential expression analysis.

## Analysis of HiChIP and C-Hi-C data

Raw FASTQ files were preprocessed with Trimmomatic-0.38 to remove Illumina sequencing adaptors. The 5' and 3' ends with a quality score below 5 (Phred+33) were trimmed and reads <30 bp after trimming were dropped. The trimmed files were then processed with HiC-Pro V2.11.4⁷⁷. The reads were aligned using Bowtie2 V2.3.5 onto the M82V1.0 assembly with default settings, except for the parameter "-score-min L, -0.6, -0.8" (Supplementary Data 4). Invalid ligation products (such as dangling ends, fragments ligated on themselves, and ligations of juxtaposed fragments) were discarded. Valid pairs were used to produce raw interaction matrixes at various resolutions. Finally, ".hic" files were generated with the software Juicer Tools and visualized with the tool Juicebox⁷⁸.

## Identification and analysis of genomic interactions

Valid pairs generated were further analyzed using HOMER V4.11⁷⁹ for different resolutions (500 bp for short-range interactions and 20 kb for long-range interactions). For the short range interactions, both anchors of the interactions (genomic bins) were annotated with genes using bedtools intersect⁸⁰. We removed interactions without any gene annotations, self-loops, and duplicates.

## Reporting summary

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

## Data availability

The data that support this study are available from the corresponding author upon reasonable request. Raw sequencing data generated in the course of this study have been deposited to the Gene Expression Omnibus (GEO) database under accession number GSE206365. Source data are provided with this paper.

## References

1. Rodriguez-Granados, N. Y. et al. Put your 3D glasses on: plant chromatin is on show. *J. Exp. Bot.* **67**, 1–17 (2016).

2. Huang, Y. et al. The matrix revolutions: towards the decoding of the plant chromatin three-dimensional reality. *J. Exp. Bot.* **71**, 5129–5147 (2020).

3. Concia, L. et al. Wheat chromatin architecture is organized in genome territories and transcription factories. *Genome Biol.* **21**, 104 (2020).

4. Doenecke, D. et al. Histones: genetic diversity and tissue-specific gene expression. *Histochem. Cell Biol.* **107**, 1–10 (1997). 1997 1071.

5. Misteli, T. The self-organizing genome: principles of genome architecture and function. *Cell* **183**, 28–45 (2020).

6. Kempfer, R. & Pombo, A. Methods for mapping 3D chromosome architecture. *Nat. Rev. Genet.* **21**, 207–226 (2020).

7. Misteli, T. & Finn, E. H. Chromatin architecture is a flexible foundation for gene expression. *Nat. Genet.* **53**, 426–427 (2021). 2021 534.

8. Lieberman-Aiden, E. et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. *Science* **326**, 289–293 (2009).

9. Dong, P. et al. 3D chromatin architecture of large plant genomes determined by local A/B compartments. *Mol. Plant.* **10**, 1497–1509 (2017).

10. Hansen, A. S. et al. Distinct classes of chromatin loops revealed by deletion of an RNA-binding region in CTCF. *Mol. Cell* **76**, 395 (2019).

11. Di Giannartino, D. C., Polyzos, A. & Apostolou, E. Transcription factors: building hubs in the 3D space. *Cell Cycle* https://doi.org/10.1080/15384101.2020.1805238 (2020).

12. Khanna, N., Zhang, Y., Lucas, J. S., Dudko, O. K. & Murre, C. Chromosome dynamics near the sol-gel phase transition dictate the timing of remote genomic interactions. *Nat. Commun.* **10**, 2771 (2019).

13. Sabari, B. R. et al. Coactivator condensation at super-enhancers links phase separation and gene control. *Sci* **361**, eaar3958 (2018).

14. Zhao, L. et al. Chromatin loops associated with active genes and heterochromatin shape rice genome architecture for transcriptional regulation. *Nat. Commun.* **10**, 1–13 (2019).

15. Liu, C., Cheng, Y. J., Wang, J. W. & Weigel, D. Prominent topologically associated domains differentiate global chromatin packing in rice from Arabidopsis. *Nat. Plants* **3**, 742–748 (2017).

16. An, L. et al. OnTAD: Hierarchical domain structure reveals the divergence of activity among TADs and boundaries. *Genome Biol.* **20**, 1–16 (2019).

17. Grob, S. Plants are not so different. *Nat. Plants* **3**, 690–691 (2017). 2017 39.

18. Kim, S., Yu, N. K. & Kaang, B. K. CTCF as a multifunctional protein in genome regulation and gene expression. *Exp. Mol. Med.* **47**, e166–e166 (2015).

19. Shang, J. Y. & He, X. J. Chromatin-remodeling complexes: conserved and plant-specific subunits in *Arabidopsis*. *J. Integr. Plant. Biol.* **64**, 499–515 (2022).

20. Zhou, S., Jiang, W., Zhao, Y. & Zhou, D. X. Single-cell three-dimensional genome structures of rice gametes and unicellular zygotes. *Nat. Plants* **5**, 795–800 (2019).
21. Bhadouriya, S. L., Mehrotra, S., Basantani, M. K., Loake, G. J. & Mehrotra, R. Role of chromatin architecture in plant stress responses: an update. *Front. Plant. Sci.* **11**, 2131 (2021).

22. Sun, L. et al. Heat stress-induced transposon activation correlates with 3D chromatin organization rearrangement in Arabidopsis. *Nature* https://doi.org/10.1038/s41467-020-15809-5 (2020).

23. Peng, Y. et al. Chromatin interaction maps reveal genetic regulation for quantitative traits in maize. *Nat. Commun.* **10**, 1–11 (2019). 2019 101.

24. Ricci, W. A. et al. Widespread long-range cis-regulatory elements in the maize genome. *Nat. Plants* **5**, 1237–1249 (2019).

25. Li, L. et al. Widespread rearrangement of 3D chromatin organization underlies polycomb-mediated stress-induced silencing. *Mol. Cell* **58**, 216–231 (2015).

26. Vertii, A. Stress as a chromatin landscape architect. *Front. Cell Dev. Biol.* **9**, 790138 (2021).

27. Mishra, S. K. et al. In the complex family of heat stress transcription factors, HsfA1 has a unique role as master regulator of thermotolerance in tomato. *Genes. Dev.* **16**, 1555–1567 (2002).

28. Liang, Z. et al. Reorganization of the 3D chromatin architecture of rice genomes during heat stress. *BMC Biol.* **19**, 53 (2021).

29. Löffler-Wirth, H., Kalcher, M. & Binder, H. OposSOM: R-package for high-dimensional portraying of genome-wide expression landscapes on bioconductor. *Bioinformatics* **31**, 3225–3227 (2015).

30. Klemm, S. L., Shipony, Z. & Greenleaf, W. J. Chromatin accessibility and the regulatory epigenome. *Nat. Rev. Genet.* **20**, 207–220 (2019).

31. Bannister, A. J. & Kouzarides, T. Regulation of chromatin by histone modifications. *Cell Res.* **21**, 381–395 (2011). 2011 213.

32. Kimura, H. Histone modifications for human epigenome analysis. *J. Hum. Genet.* **58**, 439–445 (2013).

33. Kim, T. K. et al. Widespread transcription at neuronal activity-regulated enhancers. *Nat* **465**, 182–187 (2010).

34. Furlong, E. E. M. & Levine, M. Developmental enhancers and chromosome topology. *Sci* **361**, 1341–1345 (2018).

35. Kim, S. & Shendure, J. Mechanisms of interplay between transcription factors and the 3D genome. *Mol. Cell* **76**, 306–319 (2019).

36. Nover, L. et al. *Arabidopsis* and the heat stress transcription factor world: How many heat stress transcription factors do we need? *Cell Stress. Chaperones* **6**, 177–189 (2001).

37. Scharf, K. D., Berberich, T., Ebersberger, I. & Nover, L. The plant heat stress transcription factor (Hsf) family: Structure, function and evolution. *Biochim. Biophys. Acta Gene Regul. Mech.* **1819**, 104–119 (2012).

38. El-shershaby, A. et al. Functional diversification of tomato HsfA1 factors is based on DNA binding domain properties. *Gene* **714**, 143985 (2019).

39. Zhao, J., Lu, Z., Wang, L. & Jin, B. Plant Responses to Heat Stress: Physiology, Transcription, Noncoding RNAs, and Epigenetics. *Int. J. Mol. Sci.* https://doi.org/10.3390/ijms22010117 (2020).

40. Bartlett, A. et al. Mapping genome-wide transcription-factor binding sites using DAP-seq. *Nat. Protoc.* **12**, 1659–1672 (2017).

41. Charng, Y. Y. et al. A heat-inducible transcription factor, HsfA2, is required for extension of acquired thermotolerance in Arabidopsis. *Plant. Physiol.* **143**, 251–262 (2007).

42. Dong, P., Tu, X., Liang, Z., Kang, B. H. & Zhong, S. Plant and animal chromatin three-dimensional organization: Similar structures but different functions. *J. Exp. Bot.* **71**, 5119–5128 (2020).

43. Rieder, D., Trajanoski, Z. & McNally, J. G. Transcription factories. *Front. Genet.* **3**, 221 (2012).

44. Chathoth, K. T. et al. The role of insulators and transcription in 3D chromatin organization of flies. *Genome Res.* **32**, 682–698 (2022).

45. Ray, J. et al. Chromatin conformation remains stable upon extensive transcriptional changes driven by heat shock. *PNAS* **116**, 19431–19439 (2019).

46. Zaret, K. S. & Carroll, J. S. Pioneer transcription factors: establishing competence for gene expression. *Genes. Dev.* **25**, 2227 (2011).

47. Lai, X. et al. The LEAFY floral regulator displays pioneer transcription factor properties. *Mol. Plant.* **14**, 829–837 (2021).

48. Heintzman, N. D. et al. Distinct and predictive chromatin signatures of transcriptional promoters and enhancers in the human genome. *Nat. Genet.* **39**, 311–318 (2007).

49. Oka, R. et al. Genome-wide mapping of transcriptional enhancer candidates using DNA and chromatin features in maize. *Genome Biol.* **18**, 1–24 (2017).

50. Li, Z. et al. The bread wheat epigenomic map reveals distinct chromatin architectural and evolutionary features of functional genetic elements. *Genome Biol.* **20**, 1–16 (2019).

51. Xie, Y. et al. Enhancer transcription detected in the nascent transcriptomic landscape of bread wheat. *Genome Biol.* **23**, 1–15 (2022).

52. Mifsud, B. et al. Mapping long-range promoter contacts in human cells with high-resolution capture Hi-C. *Nat. Genet.* **47**, 598–606 (2015).

53. Schoenfelder, S. et al. The pluripotent regulatory circuitry connecting promoters to their long-range interacting elements. *Genome Res.* **25**, 582–597 (2015).

54. Javierre, B. M. et al. Lineage-specific genome architecture links enhancers and non-coding disease variants to target gene promoters. *Cell* **167**, 1369–1384.e19 (2016).

55. Rubin, A. J. et al. Lineage-specific dynamic and pre-established enhancer-promoter contacts cooperate in terminal differentiation. *Nat. Genet.* **49**, 1522–1528 (2017).

56. Siersbæk, R. et al. Dynamic rewiring of promoter-anchored chromatin loops during adipocyte differentiation. *Mol. Cell* **66**, 420–435.e5 (2017).

57. Tao, S. et al. Single-cell transcriptome and network analyses unveil key transcription factors regulating mesophyll cell development in maize. *Genes* **13**, 374 (2022).

58. Ouyang, W., Xiong, D., Li, G. & Li, X. Unraveling the 3D genome architecture in plants: present and future. *Mol. Plant.* **13**, 1676–1693 (2020).

59. Campigli, D., Giammartino, D., Polyzos, A. & Apostolou, E. Transcription factors: building hubs in the 3D space. *Cell Cycle* https://doi.org/10.1080/15384101.2020.1805238 (2020).

60. Sanborn, A. L. et al. Chromatin extrusion explains key features of loop and domain formation in wild-type and engineered genomes. *Proc. Natl Acad. Sci. USA* **112**, E6456–E6465 (2015).

61. Fudenberg, G. et al. Formation of chromosomal domains by loop extrusion. *Cell Rep.* **15**, 2038–2049 (2016).

62. Kyrchanova, O. & Georgiev, P. Mechanisms of enhancer-promoter interactions in higher eukaryotes. *Int. J. Mol. Sci.* **22**, 1–20 (2021).

63. Cao, S. et al. A Distal CCAAT/NUCLEAR FACTOR Y complex promotes chromatin looping at the FLOWERING LOCUS T promoter and regulates the timing of flowering in *Arabidopsis*. *Plant. Cell* **26**, 1009–1017 (2014).

64. Chowdhary, S., Kainth, A. S., Pincus, D. & Gross, D. S. Heat shock factor 1 drives intergenic association of its target gene loci upon heat shock. *Cell Rep.* **26**, 18–28.e5 (2019).

65. Latrasse, D. et al. MAPK-triggered chromatin reprogramming by histone deacetylase in plant innate immunity. *Genome Biol.* **18**, 131 (2017).

66. Buenrostro, J. D., Giresi, P. G., Zaba, L. C., Chang, H. Y. & Greenleaf, W. J. Transposition of native chromatin for fast and sensitive epigenomic profiling of open chromatin, DNA-binding proteins and nucleosome position. *Nat. Methods* **10**, 1213–1218 (2013).

67. Jégu, T. et al. The *Arabidopsis* SWI/SNF protein BAF60 mediates seedling growth control by modulating DNA accessibility. *Genome Biol.* **18**, 114 (2017).
68. Mumbach, M. R. et al. HiChIP: Efficient and sensitive analysis of protein-directed genome architecture. *Nat. Methods* **13**, 919–922 (2016).

69. Ramirez-Prado, J. S., Latrasse, D. & Benhamed, M. Histone modification ChIP-seq on *Arabidopsis* thaliana plantlets. *Bio-Protoc.* **11**, e4211 (2021).

70. Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. *Bioinformatics* **30**, 2114–2120 (2014).

71. Alonge, M. et al. Automated assembly scaffolding elevates a new tomato system for high-throughput genome editing. *bioRxiv* https://doi.org/10.1101/2021.11.18.469135 (2021).

72. Langmead, B. & Salzberg, S. L. Fast gapped-read alignment with Bowtie 2. *Nat. Methods* **9**, 357–359 (2012).

73. Gaspar, J. Improved peak-calling with MACS2. *bioRxiv* https://doi.org/10.1101/496521 (2018).

74. Ramírez, F. et al. deepTools2: a next generation web server for deep-sequencing data analysis. *Nucleic Acids Res.* **44**, W160–W165 (2016).

75. Dobin, A. et al. STAR: Ultrafast universal RNA-seq aligner. *Bioinformatics* **29**, 15–21 (2013).

76. Love, M. I., Huber, W. & Anders, S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. *Genome Biol.* **15**, 550 (2014).

77. Servant, N. et al. HiC-Pro: An optimized and flexible pipeline for Hi-C data processing. *Genome Biol.* **16**, 259 (2015).

78. Durand, N. C. et al. Juicebox provides a visualization system for Hi-C contact maps with unlimited zoom. *Cell Syst.* **3**, 99–101 (2016).

79. Heinz, S. et al. Simple combinations of lineage-determining transcription factors prime cis-regulatory elements required for macrophage and B cell identities. *Mol. Cell* **38**, 576–589 (2010).

80. Quinlan, A. R. & Hall, I. M. BEDTools: a flexible suite of utilities for comparing genomic features. *Bioinformatics* **26**, 841–842 (2010).

## Acknowledgements

This work was supported by the European Research Council ERC (Project 101044399-3Dwheat), Agence National de la Recherche ANR (ANR-21-CE20-0036-4D Heat Tomato) and by the Institut Universitaire de France (IUF) to M.B. Moreover, this work was supported by China Scholar Council fellowships (201806690005) to Y.H. Figure 6 was created using BioRender (www.biorender.com) with publication permissions.

## Author contributions

M.B. designed the research. Y.H., M.F.M., and R.B.-C. collected the plant materials. Y.H. performed the ATAC-seq, RNA-seq, AmpDAP-seq, Capture Hi-C, Hi-ChIP, ChIP-seq, immunofluorescence staining

experiments, and the ATAC-seq, RNA-seq, AmpDAP-seq, C-Hi-C, Hi-ChIP, ChIP-seq libraries with D.L.'s help. D.L. performed the Illumina sequencing. Y.H. performed chromatin conformation capture. Y.H. performed dual-luciferase assays with C.R., F.-Q.T. and R.D.'s help. J.A., S.S., X.H., J.A.-S., D.L., C.X., Y.Z., and M.B. performed the bioinformatics analyses. Y.H., J.A., and M.B. made the figures and supplemental data. M.B., C.B., C.D.L., and Y.H. wrote the paper with support from all authors. B.D.C., A.B., C.B., J.B., F.F., M.D., Z.M., Y.Z., S.Z., C.X., F.A., F.S., M.C., M.M.M., J.G.-M., C.R., and D.L. analyzed the data and provided critical feedback. The authors read and approved the final manuscript.

## Competing interests

The authors declare no competing interests.

## Additional information

**Supplementary information** The online version contains supplementary material available at https://doi.org/10.1038/s41467-023-36227-3.

**Correspondence** and requests for materials should be addressed to Moussa Benhamed.

**Peer review information** *Nature Communications* thanks Junpeng Shi, and the other, anonymous, reviewer(s) for their contribution to the peer review of this work.

**Reprints and permissions information** is available at http://www.nature.com/reprints

**Publisher's note** Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

**Open Access** This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons license, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons license and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this license, visit http://creativecommons.org/licenses/by/4.0/.

© The Author(s) 2023

<footnote n="1">¹Université Paris-Saclay, CNRS, INRAE, Univ Evry, Institute of Plant Sciences Paris-Saclay (IPS2), 91405 Orsay, France. ²School of Life Science, University of Warwick, Coventry CV4 7AL, UK. ³Instituto de Agrobiotecnología del Litoral, CONICET, Universidad Nacional del Litoral, Colectora Ruta Nacional 168km 0, 3000 Santa Fe, Argentina. ⁴Wheat Research Institute, School of Life Sciences and Engineering, Southwest University of Science and Technology, Mianyang, Sichuan 621010, China. ⁵INRA and University of Bordeaux, UMR 1332 Biologie du Fruit et Pathologie, F-33140 Villenave d'Ornon, France. ⁶GReD, Université Clermont Auvergne, CNRS, INSERM, BP 38, 63001 Clermont-Ferrand, France. ⁷Laboratoire Génomique et Biotechnologie du Fruit (GBF), UMR990, INRA/INP-ENSAT, Castanet-Tolosan, France. ⁸State Key Laboratory of Genetic Engineering, Collaborative Innovation Center of Genetics and Development, Department of Biochemistry, Institute of Plant Biology, School of Life Sciences, Fudan University, Shanghai 200438, China. ⁹State Key Laboratory of Plant Physiology and Biochemistry, College of Life Sciences, Zhejiang University, Hangzhou 310058, China. ¹⁰Department of Biosciences, Molecular Cell Biology of Plants, Goethe University Frankfurt am Main, Max-von-Laue Str. 9, 60438 Frankfurt am Main, Germany. ¹¹Division of Biological and Environmental Sciences and Engineering, King Abdullah University of Science and Technology, Thuwal 23955-6900, Kingdom of Saudi Arabia. ¹²Université de Paris, Institute of Plant Sciences Paris-Saclay (IPS2), F-75006 Paris, France. ¹³Institut Universitaire de France (IUF), Orsay, France.✉e-mail: moussa.benhamed@u-psud.fr</footnote>