Drought-responsive dynamics of H3K9ac-marked 3D chromatin interactions are integrated by OsbZIP23-associated super-enhancer-like promoter regions in rice

Abstract

Background

In response to drought stress (DS), plants undergo complex processes that entail significant transcriptome reprogramming. However, the intricate relationship between the dynamic alterations in the three-dimensional (3D) genome and the modulation of gene co-expression in drought responses remains a relatively unexplored area.

Results

In this study, we reconstruct high-resolution 3D genome maps based on genomic regions marked by H3K9ac, an active histone modification that dynamically responds to soil water variations in rice. We discover a genome-wide disconnection of 3D genome contact upon DS with over 10,000 chromatin loops lost, which are partially recovered in the subsequent re-watering. Loops integrating promoter–promoter interactions (PPI) contribute to gene expression in addition to basal H3K9ac modifications. Moreover, H3K9ac-marked promoter regions with high affinities in mediating PPIs, termed as super-promoter regions (SPRs), integrate spatially clustered PPIs in a super-enhancer-like manner. Interestingly, the knockout mutation of OsbZIP23, a well-defined DS-responsive transcription factor, leads to the disassociation of over 80% DS-specific PPIs and decreased expression of the corresponding genes under DS. As a case study, we show how OsbZIP23 integrates the PPI cluster formation and the co-expression of four dehydrin genes, RAB16A–D, through targeting the RAB16C SPR in a stress signaling-dependent manner.

Conclusions

Our high-resolution 3D genome maps unveil the principles and details of dynamic genome folding in response to water supply variations and illustrate OsbZIP23 as an indispensable integrator of the yet unique 3D genome organization that is essential for gene co-expression under DS in rice.

Background

Rice, serving as both the world’s primary staple food and a highly susceptible crop to drought, has emerged as a pivotal model for the investigation and breeding of drought resilience in plants over the past few decades. Studies have illuminated the intricate network of coordinated expression (co-expression) that encompasses numerous drought-responsive genes, albeit with functional significance that may be comparatively minor [1, 2]. Harnessing the potential of these genes, researchers implemented strategies such as multi-gene assembly or pyramiding, leading to the development of transgenic or inbred lines that demonstrated enhanced, though not yet fully satisfactory, drought resilience [3]. Despite the rapid advances in breeding methodology, the scant understanding of the regulator mechanism of drought-responsive co-expression might be one of the main drawbacks hindering the process of rice drought-resilience improvement.

Histone post-translational modification, especially methylation and acetylation, is one of the most widely understood genomic marks in modulating gene expression [4]. Genome-wide profiling of histone H3 lysine 4 trimethylation (H3K4me3) discovered drought-triggered remodeling of H3K4me3 as well as altered expressions of the marked genes in rice [5]. A follow-up study further revealed that OsbZIP23, a master transcriptional activator in drought responses, collaborated with H3K4me3 and induced the expression of four dehydrin genes that are fundamental to drought tolerance [6]. However, only 1.6% of the H3K4me3 marked regions were specific to drought stress treatment, and the correlation between changes in H3K4me3 and differential expression levels was relatively weak, especially for genes with moderate expression [5]. Similar results were also obtained in Arabidopsis thaliana [7]. Thus, additional histone modifications might be involved to achieve a more sensitive and broad-spectrum indication for drought-responsive gene expression.

Previous studies have revealed that histone acetylation, represented by histone H3 lysine 9 acetylation (H3K9ac), generally serve as euchromatic marks that can be remodeled dynamically during developmental processes and environmental stimuli in both mammal and plant [4, 8]. In Arabidopsis thaliana, a higher level of H3K9ac can be observed at the promoters of stress-responsive genes such as AtMYB29, AtGST1, AtGSTF10, and AtANN1 when subjected to drought stress [9]. At the genome-wide level, a larger proportion of H3K9ac was found to be transiently decreased upon drought stress treatment in Maize, which marked ~ 25% of the drought stress-suppressive genes [10]. In rice, though data supporting the relationship between drought stress and H3K9ac modification was absent, H3K9ac was found to be induced by high-salt treatments and marked salt-inducible genes more significantly in seedlings than in roots [11]. These findings hint that H3K9ac possibly serves as a dynamic mark for drought-responsive genes in rice.

In addition, the modulation of transcriptome also takes place at the 3D level. Over the past few years, diverse 3D genome conformations in various organisms have been unveiled. In mammals, chromosomes are spatially organized into A/B compartments that contain active and inactive chromosome regions, respectively [12]. Each compartment is hierarchically composed of topologically associating domains (TADs), which enable interactions among regulatory elements within the boundaries while insulating interactions across boundaries [12]. In plants, although canonical TAD-like structures were also defined, they were not conserved across species [13]. Remarkably, only ~ 25% of the rice genome can be characterized by the canonical TAD concept [14]. Recent advances in high-resolution chromatin conformation capture techniques enabled closer views of rice 3D genome conformation and unveiled chromatin loops that are tethered by RNA polymerase II (RNAPII) and active/inactive histone modifications [15]. Instead of the canonical TAD, chromatin interacting domains (CIDs), i.e., genomic regions spanned by continuous chromatin loops, characterized ~ 82% of the rice genome into relatively independent spatial interacting units [15]. The absence of CTCF (CCCTC binding factor) may explain these fundamental differences in contacting domain definitions between mammals and plants [13].

Although the overall higher-order 3D genome structure is generally robust across organs and growth conditions [14, 16], spatial contact between genes and/or regulator elements can vary upon environmental stimulus. Recently, the circadian transcriptome modulation by dynamic spatial chromatin contacts was systematically profiled in rice [17, 18]. Unlike the promoter–distal element contacts identified in maize (Zea mays) and Arabidopsis thaliana [19,20,21], diurnal 3D genome remodeling triggered the formation and disassociation of time-specific promoter–promoter interacting clusters that are integrated by the promoter elements of vital circadian genes such as OsLHY and OsTOC1 [18]. However, the RNAPII-mediated ChIA-PET (chromatin interaction analysis by paired-end tag sequencing) data alone provided limited information on the regulator mechanisms of dynamic 3D genome remodeling due to the housekeeping role of RNAPII. Moreover, a study in tomato (Solanum lycopersicum) revealed that heat shock factor HSFA1a was a potential modulator for chromatin spatial re-organization under heat stress [22], implying the importance of stress-responsive transcription factors in stress-specific 3D genome remodeling. These findings emphasize the requirement for high-resolution data to elucidate the potentially unique 3D genome features in response to drought stress, which is absent for rice or other crops.

In this study, dynamic maps of the H3K9ac-marked 3D genome upon drought stress and re-watering were drawn at high-resolution using ChIA-PET. We found that genome-wide re-organizations in promoter–promoter spatial interaction (PPI) contributed to the reprogramming of gene co-expression in response to soil water variations. Such re-organizations involved disassociation and formation of PPI clusters that were integrated by promoter regions with high contact-integrating affinity, which we termed as super-promoter regions (SPRs). Moreover, the binding of OsbZIP23, a master transcription factor in rice drought stress responses, to SPRs was essential for the formation of drought-specific PPIs between functionally vital genes. These findings illustrated a novel and unique model for drought stress-driven 3D genome re-organizations in rice, in which OsbZIP23-targeting permits SPR-integrated PPI formation to promote coordinated transcription upon drought stress.

Results

H3K9ac dynamically responds to variations in water availability and is correlated with changes in the transcriptome in rice

To characterize genome-wide H3K9ac modification dynamics in response to changes in water supply, a previously established enhanced chromatin immunoprecipitation (eChIP) method was performed using MH63 (Oryza sativa L. ssp. indica cv. Minghui 63) seedlings exposed to the normal condition (NC), drought stress (DS), and subsequent re-watering (RE) (see Fig. 1a and Additional file 1: Fig. S1a, b for plant treatment and sampling procedures). Processing of the eChIP–seq reads revealed reproducible datasets with Pearson correlation coefficients (PPCs) > 0.9 between two biological replicates (Additional file 1: Fig. S1c, d). After filtering by irreproducible discovery rate (IDR) at 0.05, 18,029, 20,000, and 19,407 high-fidelity peak regions (regarded as H3K9ac-modified sites) were identified under NC, DS, and RE, respectively (Additional file 1: Fig. S1e, f; Additional file 2: Table S1). Under all three conditions, H3K9ac signals generally covered the majority of euchromatin (Additional file 1: Fig. S2a). The sequencing reads was significantly enriched around the transcription start site (TSS) (Additional file 1: Fig. S2b), and over 70% of the peaks were located in the proximal promoter regions (− 1000 – + 100 bp to TSS) (Fig. 1b), which is in consistent with the results in previous studies [8, 17, 23].

Fig. 1

H3K9ac modification is a dynamic histone mark in response to drought stress and subsequent re-watering. a A summary of the sample treatment and harvest workflow (see “Methods” for detailed descriptions). RSWC, relative soil water content. b H3K9ac peak categorizations in terms of the associated genomic features. c Statistics of the intersection relationships between H3K9ac peak sets under different conditions. Jaccard value equals the proportion of overlapped base pairs in the unified regions of the two peak sets. d Peak span comparisons among different conditions. The range of the hinges in the boxplot indicates the 25% and 75% quantiles; whiskers denote the 1.5 × inter-quartile range (IQR); the line inside the box indicates the median value. Statistical significances were determined by one-way analysis of variance (ANOVA) with Tukey’s Honest significant difference (HSD) test. e Genome track snapshots illustrating the dynamic changes in H3K9ac modifications and gene expressions at three representative loci under NC, DS, and RE. RPGC, reads per genomic content (per 1-bp bin); RPKM, reads per kilobase per million mapped reads (per 1-bp bin). f Volcano plots showing the quantitatively differential H3K9ac modifications. Numbers of the differential peaks are labeled. g Correlations between the fold changes of H3K9ac peak intensities and the fold changes of associated gene expression. R indicates Pearson correlation coefficient value, and P-values were determined by Kendall’s test

Upon DS and the following RE, the genomic coverages of H3K9ac-marked regions varied dynamically with Jaccard values ranging from 0.58 to 0.77 (Fig. 1c). Though the H3K9ac peak number under DS was larger than that under NC or RE, its average span was significantly narrower (Fig. 1d). This can be explained by the splitting of 980 larger NC-peaks into 2211 narrower ones under DS (Fig. 1c). In support of the statistics, we observed the DS-specific H3K9ac modification proximal to the TSS of DS-inducible genes such as OsDSSR1 [24] (Fig. 1d), as well as the DS-induced split of NC-peak proximal to OsbZIP23 [25, 26] (Fig. 1d). Besides, the narrowed peak spans under DS can also be reflected by the decreased H3K9ac modifications around DS-repressive genes like DST [27] (Fig. 1d). In addition, the DS treatment led to more intensity-increased H3K9ac modification regions than decreased ones (4979 vs 2477, Fig. 1f and Additional file 2: Table S2). Furthermore, re-watering after DS recovered the H3K9ac intensities in most regions resulted in only 1026 differentially marked regions when compared with the NC-peaks (Fig. 1f, Additional file 2: Table S2).

To evaluate the correlation between changes in H3K9ac modification and gene expression alternations, we calculated the PCC values between the H3K9ac intensity fold changes in the promoter regions and the corresponding gene expression fold changes assayed by RNA-seq (Additional file 1: Fig. S3 and Additional file 2: Table S3) upon DS treatment and re-watering. The correlations between H3K9ac intensity fold changes and gene expression fold changes were significant, with PCC > 0.6 upon DS and subsequent RE (Fig. 1g). Moreover, the correlation between the NC–RE fold changes was weak (PCC = 0.16, Fig. 1g), implying that H3K9ac was an efficient mark for dynamic expression changes upon episodic water supply transitions rather than an indicator for prolonged transcriptome influences by DS. These results indicate that H3K9ac modification may serve as an active chromatin mark that can respond to variations in water supply dynamically.

Drought stress and re-watering dynamically remodeled the landscape of H3K9ac-marked 3D genome architecture

The dynamic changes of H3K9ac in response to DS and RE prompted us to analyze the profiles of the potentially H3K9ac-associated 3D genomic interactions. Using the well-established long-read ChIA-PET method [28], we obtained H3K9ac-marked 3D genomic architectures under NC, DS, and RE with high repeatability (PCC > 0.8 based on the contact intensities between biological replicates at 100-kb resolution, Additional file 1: Fig. S4a–c). The resulting contact matrices revealed genome-wide dynamic changes in spatial interaction intensities upon DS and RE (Fig. 2a, b). If viewed at a higher resolution (i.e., 5-kb), the re-organizations of 3D contacts between genes in rice can occur in relatively shorter ranges within the traditional TAD structures (Fig. 2c). This result implicates the complexity in rice 3D genome re-programming in response to DS and RE. Thus, we performed the subsequent analysis of the H3K9ac-marked 3D genome based on chromatin loops.

Fig. 2

Overview of the H3K9ac-marked 3D genome dynamics in response to drought stress and subsequent re-watering. a Contact matrices showing the relative contact frequency determined by the normalized PET count at 100-kb resolution. b Zoomed views on the marked region in a at 100-kb resolution. c A genome track snapshot showing the spatial contacts integrated by chromatin loops and the gene expression abundances within the region marked in b. The ChIP-seq and RNA-seq signals were quantified by RPGC and RPKM, respectively, in 1-bp bin. d Genome-wide distributions of H3K9ac-marked chromatin interacting domains (CIDs) reflecting the changes in local chromatin interaction landscapes in response to drought stress and re-watering. Regions marked by cross refer to centromeres. e Comparing of the genomic spans covered by the CIDs shown in d. Statistical significances were determined by the ANOVA–HSD test. f Statistics of the overlapping and distinct CID regions. g Circos plots showing the inter-chromosome loops under NC, DS, and RE. The red quadrangles denote centromeres. h Chr01–Chr05 inter-chromosome contact matrices showing the contact intensity fold changes between different conditions (at 100-kb resolution). i The zoomed view on the contacts between the 1–4-Mb region of Chr01 and the whole Chr05. j Inter-chromosome loops connecting Chr01 and Chr05 under each condition. Black curves refer to the interaction marked in i. k Evaluation on the Chr01–Chr05 interaction labeled in i and j by FISH (fluorescence in situ hybridization). Genomic positions of the target anchors are labeled (red refers to the cy3 tag and green refers to the Alexa 488 tag). The crossbar plot denotes mean ± SD (standard deviation) of the interaction frequency based on three biological replicates. n, the number of nuclei observed. The P-value was determined by Fisher’s exact test

With a stringent cut-off of 6 on the inter-ligation pair-end tag (i-PET) count, 24,235, 9577, and 10,626 H3K9ac-marked chromatin loops were identified under NC, DS, and RE, respectively (Additional file 1: Fig. S5a and Additional file 2: Table S4). The chromatin loops showed considerable robustness with 19,199 (79.22%), 9204 (96.10%), and 8626 (81.18%) of the loops reproduced by BL-Hi-C (Additional file 1: Fig. S6a). To illustrate local chromatin interaction profiles, we constructed chromatin-interacting domains (CIDs) based on the continuously connected local loops under the three conditions via a similar method to a previous study [15] (Fig. 2d and Additional file 2: Table S5). Genomic coverage of the CIDs changed dramatically upon DS treatment with decreases in total coverage (88 Mb under DS vs 148 Mb under NC, Fig. 2d) and average CID span (137 kb under DS vs193 kb under NC, Fig. 2e), which was a result of decreased local loop number under DS (9061 under DS vs 17,601 under NC). The number and coverage of CIDs produced upon DS were limiting (16 regions, accounting for 0.88% of the total DS-CID coverage), while most DS-CIDs (624) overlapped with the NC-CID regions but with narrower spans (Fig. 2d, f; Additional file 1: Fig. S5b). This result is in accordance with the knowledge that higher-order 3D genome structures are generally robust under different conditions [14]. Of note, re-watering after DS treatment conferred limited or no recovery on the local loop number (8280) and the CID coverage (104 kb in average and 78 Mb in total) (Fig. 2d, f). Moreover, expression levels of the genes within CIDs were significantly higher than those outside CIDs under all three conditions (Additional file 1: Fig. S5c), indicating that H3K9ac-marked CIDs can generally characterize the active chromosome regions.

The DS treatment also decondensed the inter-chromosome interactions, resulting in merely 315 chromatin loops left (Fig. 2g). In contrast to the situation for intra-chromosome interaction, the number of inter-chromosome interactions recovered dramatically after re-watering (with 1939 loops, Fig. 2g). We found a hot-spot on chromosome 1 (Chr01: 2,254,083–2280089) that dynamically integrated 56, 0, and 43 inter-chromosome loops under NC, DS, and RE, respectively. Fluorescence in situ hybridization (FISH) assay probing this hot-spot and the ERF74 promoter-associating anchor on chromosome 5 (Chr05: 25,483,182–25,501,247) further supported DS-triggered spatial contact disassociation at these sites (Fig. 2h–k). In general, the DS treatment led to decondensed intra- and inter-chromosomal contacts. At the same time, subsequent re-watering considerably recovered the inter-chromosomal interactions but had a limited recovery effect on intra-chromosome interaction in terms of chromatin loop number.

Co-expression of DS-responsive genes are controlled by DS-specific chromatin loops marked by H3K9ac

Next, we analyzed how H3K9ac-marked 3D genome remodeling affected the coordinated gene expression upon DS and the following RE. Considering H3K9ac’s main role as a transcription activation tag in the promoter regions, we checked the chromatin loops with anchors located in the promoter regions. The promoter–promoter interacting (PPI) loops, as expected, made up the largest proportion (74.1%, 73.7%, and 76.0% under NC, DS, and RE, respectively, Fig. 3a) and was focused in the following analysis. Under all three conditions, more than half of the H3K9ac-modified regions were connected by PPI loops and the corresponding promoters were classified as PPI mode (Fig. 3b). Instead, H3K9ac-marked promoters not connected by PPI loops were classified as basal mode (Fig. 3c). We then extracted the transcript abundances of the genes in PPI mode and basal mode and the statistics showed that PPI mode conferred higher expressions compared to the basal mode or the unmarked promoters (Fig. 3d). This result illustrates that chromatin interaction between the H3K9ac-marked promoter has additional induction effect on gene expression under the three conditions.

Fig. 3

Condition-specific promoter–promoter interactions (PPIs) contributed to the dynamic changes in gene expression upon drought stress and re-watering. a Proportions of chromatin loops categorized according to the associated genomic features. b Proportions of the PPI anchors and non-PPI H3K9ac peaks in the total H3K9ac peaks under each condition. c Schematic model of H3K9ac-marked PPI and basal H3K9ac modification and their impacts on gene transcription. For each model, a genomic track snapshot is illustrated as an example. Abundances of the RNA-seq and ChIP-seq signals were quantified by RPKM and RPGC (per 1-bp bin), respectively. d Comparisons of the transcription levels (log₂(TPM + 1)) of genes involved in PPI, or basal H3K9ac modification, or non-H3K9ac modification (unmarked). Statistical significances were determined by the ANOVA-HSD test. e The differential loop analysis among the three conditions. Intersections and distinctions of total loops as well as PPI loops are illustrated in the upset plot. Total loop count and PPI loop count in each intersection set are labeled in black and orange, respectively. If a loop can only be detected under a unique condition, it was considered as a specific loop to that condition. And if a loop occurred under all three conditions, it was characterized as a consensus loop. The six boxplot panels on the right shows the statistics on the transcription levels and fold changes (log₂FC) of the PPI genes (with TPM > 1) in each intersection set. Statistical significances in the boxplots were determined by the ANOVA–HSD test. f Expression dynamics of the co-expression modules categorized by Mfuzz clustering on genes with TPM > 1. g Comparisons of the co-expression memberships of the DS-specific PPI genes to those of the other genes in each of the six modules in f. Statistical significances were determined by the ANOVA–HSD test. h Schematic networks composed of PPI loops in Chr04 under each condition. Each node in the network refers to a chromatin loop anchor and each edge refers to a PPI loop. The node color denotes the co-expression module of the anchored gene. The size of a node denotes the co-expression membership of the anchored genes to module 1 (for the NC and RE networks) and module 2 (for the DS network), respectively. Condition-specific PPI loops are highlighted by the color referring to the corresponding condition. The positions of nodes are fixed among the three networks

Next, we evaluated whether the differential chromatin loops contributed to the transcriptome changes under each condition. By comparing whether any two loops, each from different conditions, connected the same pair of anchors, we classified the loops under the three conditions into seven intersection sets (Fig. 3e, Additional file 2: Table S6). Under DS, 16,989 NC-loops (including 12,484 PPI loops, Set 1 plus Set 2 in Fig. 3e) were lost, and only 2331 loops (including 1548 PPI loops, Set 6 plus Set 7 in Fig. 3e) were formed. Among them, 2107 loops (including 1432 PPI loops) were loops unique to DS. These loops were termed DS-specific loops (Set 6 in Fig. 3e). We noticed that these DS-specific PPI loops connected a group of genes vigorously induced by DS and suppressed by RE (Fig. 3e), indicating that these loops may play a vital regulatory role in the dynamic gene expression changes in response to DS. Under RE, genes connected by the recovered PPI loops (1672, Set 4 in Fig. 3e) and RE-specific PPI loops (2417, Set 3 in Fig. 3e) showed the highest expression fold changes (in average) upon RE (Fig. 3e). Furthermore, validation by BL-Hi-C reproduced 63.96–92.56% of the ChIA-PET differential chromatin loops in response to DS and RE (Additional file 1: Fig. S6b). Especially, for the disassociated loops, H3K9ac modification abundance of the unreproduced fraction was not significantly lower than that of the reproduced fraction (Additional file 1: Fig. S6b), indicating that loop disassociations identified by ChIA-PET were not mainly due to impaired H3K9ac modifications. Together, these results suggest that the formation and disassociation of H3K9ac-marked PPI loops are associated with the dynamic expression changes of the DS-responsive genes.

To further confirm if the condition-specific loop formation is correlated with global gene co-expression, we clustered all expressed genes (with TPM > 1 under at least one condition) into six non-redundant (by the similarity of 0.2, Additional file 1: Fig. S7a) co-expression modules (Fig. 3f and Additional file 2: Table S7). The membership values of the DS-specific PPI-connected genes to the six modules were compared, and the result revealed that these genes were significantly correlated with the expression patterns of modules 2 and 3 with clear DS-induction and RE-suppression trends (Fig. 3f, g; Additional file 1: Fig. S7b). Gene ontology (GO) enrichment analysis showed significant enrichment of genes related to abiotic stress-response in module 2 (Additional file 1: Fig. S7c). Similarly, PPI loops specifically existed under NC and RE connected the genes in module 1, which was represented by genes participating in plant growth or having potentially negative effects on drought response (Additional file 1: Fig. S7c). The correlation of PPI loops to gene co-expression was further illustrated by the PPI network of chromosome 4 under the three conditions (Fig. 3h). These results demonstrate that PPI loops, especially the condition-specific loops, play a crucial role in the regulation of dynamic co-expression of genes in response to water supply variations.

Dynamic changes in PPI networks under DS and RE are integrated by super-enhancer-like regions in the promoters

As illustrated by the PPI network (Fig. 3h; Additional file 2: Tables S8–S10), nodes (anchors) were not equally connected, and some nodes may obtain significantly higher degrees (number of times the node is connected in the network). This phenomenon is similar to the interactions integrated by super-enhancers in mammals and Arabidopsis [21, 29]. Therefore, we suspect that regions or elements similar to the super-enhancer may exist in rice as a composition of the proximal promoter. To characterize these regions, we classified the genome-wide PPIs into closely connected spatial clusters using a random walk algorithm and located the hub-anchor with the largest degree value within each cluster (Fig. 4a; Additional file 2: Table S8–S10). Then, a contact-integrating affinity (CIA) value (see Fig. 4a and Methods) was assigned to each hub-anchor to evaluate its PPI-integrating ability free from the influence of anchor span (Additional file 1: Fig. S8a, b). Hub-anchors with CIA value exceeding the 80% quantile were defined as super-anchors under the condition (Fig. 4b). Since the analysis was based on PPIs, the super-anchors regions were further termed super-promoter regions (SPRs).

Fig. 4

Identification of the super-promoter regions. a The working model for identifying PPI clusters, hub-anchors, and super-promoter regions (SPRs). First, closely connected PPIs were categorized into PPI clusters using the random walk algorithm. Then, the most frequently connected anchor in the cluster was defined as the hub-anchor. Next, a contact-integrating affinity (CIA) value was calculated for each hub-anchor. Finally, anchors with CIA values exceeding the 80% quantile were identified as SPRs under the condition. b Ranking of the hub-anchors according to their CIA values under each condition. The CIA thresholds for SPR are marked by dashed lines. To illustrate the variations of SPRs under different conditions, the SPRs under DS are sequentially marked by gradient colors so that variations in their CIA values can be traced in the NC and RE plots. c Heatmap showing the CIA values of all the SPRs under the three conditions. Condition-specific SPRs, whose CIA values under one condition were 1.5-fold higher than those under the other two conditions, are labeled below the heatmap. Names of the representative drought-responsive (in brown) and suppressive (in green) genes involved in the PPI clusters that were integrated by the condition-specific SPRs are as labeled. d Pearson correlation coefficients based on the CIA values of the SPRs under each condition. e Comparisons on the expression fold changes (log₂FC) of the condition-specific PPI genes and percentage condition-specific PPI (in the total PPIs) within the PPI clusters integrated by condition-specific SPRs to those within the PPI clusters integrated by other hub-anchors. P-values were determined by two-sided Wilcoxon test. f and h Comparative views on the PPI networks in Chr11 under NC and DS. Each node in the networks refers to a loop anchor and an edge refers to a PPI loop. The color of the node denotes the co-expression module of the anchored gene. The size of a node denotes the module-membership of the anchored genes to module 1 (for NC and RE networks) and module 2 (for the DS network), respectively. Condition-specific PPI loops are highlighted by the color refer to the condition. Green and brown quadrangles represent NC and DS SPRs, respectively. The node positions in the two networks are fixed and illustrated in mirror symmetry. The SPR-integrated NC-PPI clusters that disassociated under DS are marked by a dashed line. g Genome track view of the disassociated NC-PPI cluster in f and h. ChIP-seq and RNA-seq quantifications are shown in RPGC and RPKM (per 1-bp bin), respectively

We noticed that the PPI-integrating ability of SPRs changed dramatically under different conditions, illustrated by the DS-SPRs as landmarks (Fig. 4b). Fifty five, 68, and 35 SPRs specifically existed under NC, DS, and RE, respectively, with a CIA fold change cutoff of 1.5 (Fig. 4c). Of note, 67 of the DS-specific SPRs also obtained higher CIA values under DS than under NC and/or RE when validated by the BL-Hi-C PPI networks (Additional file 1: Fig. S6c), illustrating the robustness of the SPR identifications. Besides, the recovery of 3D genome after RE was more reflected by PPCs of the CIAs of SPRs than by CIDs or unclassified loops (Fig. 4d), suggesting that the contacts integrated by the SPRs may play a priority role in the remodeling of 3D genome. Indeed, the PPI clusters integrated by the DS-specific SPRs involved many important drought-inducible genes such as SNAC1 [30], SAPK6 [31], RAB16A–D [6], and OsTZF5 [32], while those by NC-specific SPRs involved typical drought-suppressive genes such as OMTN3 [33] (Fig. 4c). Furthermore, PPI clusters integrated by the specific SPRs under each condition showed a higher proportion of condition-specific PPI loops, as well as higher expression fold changes in genes connected by the condition-specific PPIs compared to those by non-specific hub-anchors (Fig. 4e). The result indicated that PPI clusters integrated by the condition-specific SPRs may have greater contribution to the dynamic changes on 3D genome and transcriptome in response to DS and RE. The dynamics of the PPI cluster integrated by the SPR upstream OsMH_11G0045700 (Chr11: 2,500,675–2502560) was illustrated as an example (Fig. 4f–h). In addition, motif enrichment analysis was performed to investigate transcription factors associated with the DS-specific SPRs. The result revealed the existence of binding sites for transcription factors with potential functions in DS response, such as bZIP and WRKY transcription factors (Additional file 1: Fig. S8c). These results highlight the existence of super-enhancer-like elements in the proximal promoter regions in rice and outline the reprogramming of the rice 3D genome upon episodic DS and RE through the alternations in the PPIs integrated by the SPRs.

OsbZIP23 is essential for the formation of DS-specific PPIs

Previous works implicated that OsbZIP23 not only functions as a direct activator of gene transcription but also acts as a potential mediator for chromatin modifications [6, 25]. The enriched bZIP motif in the DS-specific SPRs prompted us to investigate the potential role of OsbZIP23 in 3D genome remodeling under DS. The expression of OsbZIP23 was significantly induced under DS and reduced upon RE (Fig. S9a), in accordance with the previous result [26]. Using the long-read ChIA-PET, we identified the H3K9ac-marked chromatin architecture under DS in a previously reported osbzip23 T-DNA mutant [25] with high reproducibility (PCC = 0.93, based on the 100-kb resolution contact matrices of the two biological replicates, Additional file 1: Fig. S9b). Compared to the MH63 wild-type (WT) data, mutation of OsbZIP23 significantly reshaped the landscape of the 3D genome under DS (Fig. 5a). Based on the 50-kb contact matrices, the PCC between osbzip23 and the other MH63-WT samples were relatively low (0.71–0.76, Additional file 1: Fig. S9c). Identification of chromatin loops using the same criteria as the MH63-WT samples revealed fewer loops (8366, Additional file 2: Table S11) in osbzip23, even when the total i-PET count in osbzip23 was almost twofold of that in the MH63-WT DS sample (Additional file 2: Table S12). In addition, CID construction using the local loops (8007) revealed less (439) but larger (198 kb in average) CIDs in osbzip23 compared to the DS-treated MH63-WT (Fig. 5b, c; Additional file 2: Table S13). The median intra-chromosome loop span was also larger than that in MH63-WT under DS (Additional file 1: Fig. S9d). However, the average gene expression within the CIDs in osbzip23 was significantly lower than that in MH63-WT under DS (Fig. 5d), indicating that the genes within the active chromosome were less induced by DS in the osbzip23 mutant.

Fig. 5

Knockout of OsbZIP23 altered the 3D genome landscape under DS. a The contact matrix showing the relative contact intensities in terms of osbzip23 vs MH63-WT under DS (at 50-kb resolution). b The genome-wide distribution of H3K9ac-marked CIDs in MH63-WT and osbzip23 under DS. c Comparisons of the CID spans in osbzip23 under DS with those in MH63-WT under NC and DS. Statistical significances were determined by the ANOVA–HSD test. d Statistics on the CID gene expression levels (log₂(TPM + 1)) under DS in MH63-WT and osbzip23. The P-value was calculated by two-sided Wilcoxon test. e Intersection analysis on the total loops as well as the PPI loops between MH63-WT and osbzip23. The proportion of DS-specific loops in each intersection set is marked in brown. f Boxplots illustrating the statistical analysis on the expression levels (log₂(TPM + 1)) and fold changes (log₂FC) of the PPI genes in each intersection group. Statistical significances were determined by the ANOVA–HSD test. g Ranking of the hub-anchors in terms of their CIA values. The DS SPRs in MH63-WT are marked by gradient colors for tracing in osbzip23. h Relative intensities of anti-OsbZIP23 ChIP-seq reads along the genomic span of DS PPI loops. i The proportion of OsbZIP23-binding PPIs in the DS-specific PPIs that disassociated in osbzip23 under DS. j The heatmap showing the CIA values of SPRs in MH63-WT and osbzip23 under DS. The DS-specific SPRs with 1.5-fold or more CIA decreases in osbzip23 were termed suppressed SPRs in osbzip23. And the suppressed SPRs with OsbZIP23 binding site(s) were regarded as SPRs targeted by OsbZIP23. Names of the representative DS-responsive genes involved in the PPI clusters integrated by the OsbZIP23-targeted SPRs are highlighted. The boxplot shows the relative OsbZIP23-binding intensities (IP vs input, under DS) within the total DS-specific SPRs and the suppressed ones in osbzip23. The P-value was determined by two-sided Wilcoxon test. k The expression comparison of genes involved in the OsbZIP23-targeted SPR-integrated PPI clusters to those involved in clusters integrated by other DS-specific SPRs. The P-value was determined by two-sided Wilcoxon test. l, m The Chr11 PPI networks in MH63-WT and osbzip23 under DS. The OsbZIP23-targeted SPR-integrated PPI cluster involving four dehydrin genes are highlighted by dashed circles. n Genome track showing the spatial contacts between the promoters of the four dehydrin genes (RAB16A–D), which were integrated by the OsbZIP23-targeted SPR in the RAB16C promoter. The ChIP-seq and RNA-seq signals are shown in RPGC and RPKM (per 1-bp bin), respectively

Next, we compared the loops in osbzip23 to those in MH63-WT under DS. The result showed 4923 disassociated MH63-WT loops (including 3409 PPI loops, Set 1 in Fig. 5e) and 3712 osbzip23-specific loops (including 2689 PPI loops, Set 3 in Fig. 5e). Strikingly, about 80% (1179) of the DS-specific PPI loops in the MH63-WT were lost in osbzip23 (Set 1 in Fig. 5e; Additional file 2: Table S14), leading to decreases in average expression of drought-responsive genes in osbzip23 (illustrated by the expression fold changes for Set 1 in Fig. 5f). We analyzed the SPRs in osbzip23 and found that the CIA values of DS-SPRs in MH63-WT also changed in osbzip23 (Fig. 5g). Such a dramatic influence on 3D genome conformation by OsbZIP23 mutation led us to suspect if OsbZIP23 directly participates in the formation of DS-specific loops. Therefore, we examined the OsbZIP23 binding sites by ChIP-seq using a previously approved antibody [26, 34]. The data revealed 4026 (NC) and 8961 (DS) high-fidelity OsbZIP23 binding sites with IDR < 0.05 by two biological replicates (Additional file 1: Fig. S10 and Additional file 2: Table S15). We observed a sharp enrichment of OsbZIP23 binding signals under DS, but not under NC, in the anchor regions of DS-PPI loops (Fig. 5h). We also found that ~ 75% (881) DS-specific loops that disassociated in osbzip23 involved at least one anchor directly bound by OsbZIP23 (Fig. 5i; Additional file 2: Table S14). In addition, the bZIP motif identified in DS-specific SPRs was significantly enriched in the center of the OsbZIP23 binding sites specific to DS (n = 7254, with Fold change > 1 and FDR < 0.05 in terms of DS vs NC) and aligned with the OsbZIP23 motif at the ACGT core element (Fig. S8d). Together, these results suggest that OsbZIP23 binding at the anchor regions may play a predominant role in the formation of DS-specific loops.

To further evaluate the role of OsbZIP23 in DS-triggered 3D genome remodeling, we investigated OsbZIP23’s effect on DS-specific SPR formation under DS. Among the 68 DS-specific SPRs, over 50% (37) showed more than 1.5-fold decreases in CIA values in the osbzip23 compared to MH63-WT, and these SPRs were regarded as suppressed SPRs in osbzip23 (Fig. 5j; see CIA values in Additional file 2: Table S9 and Table S16). Compared with the total DS-specific SPRs, the suppressed SPRs in osbzip23 were more intensively bound by OsbZIP23 in MH63-WT (P = 0.033 by two-sided Wilcoxon test) and 18 of these 37 SPRs, termed as OsbZIP23-targeted SPRs, were directly bound by OsbZIP23 (Fig. 5j). The PPI clusters integrated by the OsbZIP23-targeted SPRs involved genes with higher expression levels under DS (Fig. 5j, k; in Additional file 2: Table S9). As a representative of these PPI clusters, we illustrate the DS-specific PPIs, integrated by an SPR at RAB16C, among four YSK2 dehydrin genes (RAB16A–D, Fig. 5l–n) that play vital roles in rice drought resistance [6, 35]. The PPIs as well as expressions of these four genes were hardly detected under NC or RE but dramatically induced under DS (the DS-specific PPIs were confirmed by both ChIA-PET and BL-Hi-C, Fig. 5n and Additional file 1: Fig. S11a), which is in accordance with previous reports [6, 26]. In the DS-treated osbzip23 mutant, however, H3K9ac modification and PPIs of the four genes were almost wiped out and the expressions of at least three of the genes were significantly decreased (Fig. 5n and Additional file 2: Table S17). DS-specific binding of OsbZIP23 to the promoter of RAB16A–D was also confirmed by ChIP-qPCR (Fig. S11b). These results suggest the existence of a potential 3D co-expression hub mediated by OsbZIP23, which is summarized in a schematic model presented in Fig. 6a.

Fig. 6

OsbZIP23 was crucial for the formation of the DS-specific PPI cluster and co-expression among the RAB16 genes. a A schematic model illustrating the OsbZIP23-mediated PPI cluster integrating the co-expression of four dehydrin genes (RAB16A–D). In the absence of stress stimulus or OsbZIP23, the four RAB16s were spatially distant to each other, which led to impaired co-expression. Under DS, OsbZIP23 can mediate the formation of PPIs among RAB16s, thus providing structural basis for the co-expression of RAB16s by potential activation complex(es). The undetected but potential binding of OsbZIP23 at the promoter of RAB16D was marked in grey. b RT-qPCR evaluation of the stress-specific co-expression that modulated by OsbZIP23-mediated PPIs among the RAB16s. Five µM of exogenous ABA (abscisic acid) was used as an equivalent for DS signal. The relative expression levels were determined by 2^−ddCt, and the ABA-treated MH63-WT protoplast without the RAB16C promoter-targeting gRNA was used as the reference sample (quantifications for this sample were set as 1). Schematic models were exhibited to help understand the result under each transfection condition. The crossbar plot denotes mean ± SD (standard deviation) of the relative expression level. Statistical significances were determined by Fisher’s least-significant difference (LSD) test. c A schematic diagram summarizing the principles in dynamic 3D genome re-organizations under DS and RE in rice. In brief, DS caused the disassociation of chromatin loops and the suppression of NC-favorable gene co-expression. On the other hand, DS also triggered the formation of PPIs among DS-responsive genes via the mediation of OsbZIP23 at the SPRs, which led to activated co-expression of these genes. In the following RE, the 3D genome landscape was partially recovered, and the co-expression of DS-responsive genes was closed

We further validated the essentiality of OsbZIP23 for the formation of the RAB16 PPI cluster and co-expression of the four genes using a dead-Cas9 (dCas9)-guided activation system in rice protoplasts (see Methods). Considering that OsbZIP23 is a key component in the abscisic acid (ABA) signal pathway and its expression as well as transcriptional activity are induced under exogenous ABA treatment [26], and the expression of the RAB16s was also highly responsive to ABA treatment [6], we used exogenous ABA as an equivalent for DS treatment in this assay. The result showed that the RAB16C SPR-targeting dCas9–VP16 can simultaneously induce the expression of the four RAB16s only in the presence of OsbZIP23 and exogenous ABA (Fig. 6b). In the osbzip23-derived protoplast, the RAB16C-targeting dCas9–VP16 only functioned as a basal activator as indicated by the induced expression of RAB16C but not the other three RAB16s, even under ABA treatment (Fig. 6b). Moreover, basal but not coordinated activation was also observed in the MH63-WT transfection without ABA treatment (Fig. 6b), which supported the sequencing result that formation of the spatial cluster was a stress signal-dependent process (Fig. S11a). In all, this case study elucidates the essentiality of OsbZIP23-mediated DS-specific PPI clusters to the co-expression of DS-responsive genes.

Discussion

The ability to react to and survive water-insufficient conditions is a landmark in the evolutionary route of plant landing. Previous studies have unveiled diverse features in the transcription programming upon drought stress, while the mechanism(s) modulating the coordinated changes in such a great number of genes remains elusive [2]. Suggestions have been made that, like the situations in developmental stage transitions, re-organization in chromatin spatial contacts may play vital roles in the modulation of drought-responsive gene expression. However, genome-wide evidence for this hypothesis is still absent. Furthermore, previous rice 3D genome studies detected limited spatial reprogramming at the TAD or A/B compartment levels in response to cold or heat stress and failed to draw direct links between the 3D genome re-organization and the co-expression alternations of stress-responsive genes, especially those with vital functions [14, 36].

With the high-resolution chromatin interaction maps by the ChIA-PET technology in this study, we managed to discover the comprehensive connections between genes beyond the detection of TAD or compartment-based methodology (Figs. 2c and 5n). In contrast to the limited influences of cold or heat stress in terms of the layout of TADs or compartments, drought stress drastically reshaped the rice 3D genome with over 10,000 disassociated PPI loops as well as significant decreases in average loop range and CID span (Figs. 2d and 3e; Additional file 1: Fig. S9d). Most strikingly, over 90% of inter-chromosome loops were lost under DS. Such a huge decondensation in chromatin spatial interaction is counter-intuitive since dehydration is expected to lead to more condensed cell components. However, the greater number of downregulated genes (than the upregulated ones) under DS (Additional file 1: Fig. S3b) and the low DS vs NC expression fold changes of NC-specific PPI-connected genes (Fig. 3e) indicate such phenomenon is in accordance with the transcriptome alternation. In addition, ~ 1500 PPI loops emerged upon DS and connected the highly inducible genes under DS (Fig. 3e). These results imply that the reshaping of rice 3D genome under the drought condition was an organized transcriptome regulatory process rather than merely damages caused by the unfavorable physiological condition.

It is also noticeable that, in contrast to the nearly complete recovery of H3K9ac modification and gene expression upon RE, the 3D genome conformation was only partially recovered (Figs. 2e and 3e). One possible explanation is that the restorations on H3K9ac modification is highly dynamic to water supply and happens prior to the rebuilt of chromatin interactions, and the sampling timing in our design (4 days after RE) may be insufficient for complete recovering of the chromatin loops. Nevertheless, re-formed PPI loops after RE did connect the genes with high induction levels by RE (Set 4 in Fig. 3e), implying that PPIs between genes with potential functions in RE may be priorly recovered. In addition, the DS PPI-loops preserved after RE connected highly upregulated genes under DS (Set 7 in Fig. 3e). This may account for certain “memory” mechanisms for rapid re-activation of these DS-responsive genes in case of repeated water shortage.

The existence of super-enhancer or super-enhancer-like elements in plants has not been systematically identified until a recent study in Arabidopsis that discovered intergenic super-enhancers with similar attributes to those in mammals [21]. Our results, as well as recent rice 3D genome studies, show that chromatin loops more intensively connect certain anchors in the proximal promoter regions, and they may function as super-enhancer-like regions in integrating chromatin spatial contacts (Fig. 4b) [29]. Like the super-enhancers in other species, when analyzed using a network approach, the degrees of such anchors were strongly correlated with their genomic spans (Additional file 1: Fig. S8a) [29]. Nevertheless, the concept of super-enhancer has been challenged by the fact that the spatial contact-integrating ability of a super-enhancer may be a collective result of all the elements within the large region and thus may not be regarded as an attribute of a single unit [37]. In this regard, we invented the super-promoter regions (SPRs) by considering the region size and network scale (see Methods). Our analysis shows that SPRs defined in this manner are informative since the PPI clusters they integrated generally contains higher percentage of condition-specific loops and genes with greater expression changes in response to the corresponding treatments (Fig. 4e). Some of the SPRs may integrate few interactions within relatively tiny PPI clusters, but these clusters can also involve vital genes in DS response, as represented by the interactions among the RAB16A–D genes (Fig. 5n), and are likely to be ignored if judged by the traditional super-enhancer concept [29]. Thus, the definition of SPRs may be more practically useful than super-enhancer or super-enhancer-like elements, at least for studying rice 3D genome dynamics in DS response. Our result demonstrates the enhancer-like role of certain rice promoters, which is, to date, a unique 3D genome regulatory mechanism in rice. However, the existence of distal intergenic enhancers in rice cannot be excluded since the H3K9ac modification did not efficiently mark these regions (Additional file 1: Fig. S2b).

As plants lack a CTCF homolog, transcription factors are suspected to function as potential mediators or regulators of the plant 3D genome. Recent studies provided support to this hypothesis in plants with TCP1 in Marchantia and HSFA1a in tomato [22, 38]. However, their impact on 3D genome structure was either limited or not convincing with genome-wide evidence, and most importantly, the factor(s) modulating the DS-specific 3D genome re-organization have not been identified. In this study, we profiled the genome-wide effect of OsbZIP23 mutation on the H3K9ac-marked fraction of rice 3D genome under DS and discovered the dissociation of over 80% DS-specific PPI loops (Fig. 5e). So far, no other report has shown stronger genome-wide regulatory effect of a plant transcription factor on 3D genome structure. Furthermore, as the DS-specific loops were uniquely formed upon DS, we can also conclude from the result that OsbZIP23 may have a vital role in the loop formation rather than maintaining under DS (considering that 61.24%, i.e., 3353 of the PPIs in Set 2 and 5 in Fig. 3e were preserved in osbzip23; Additional file 2: Table S9 and Table S14).

We also confirmed the binding of OsbZIP23 to at least one anchor in ~ 75% DS-specific loops that were lost in osbzip23 (Fig. 5h, i) by ChIP-seq. These results suggested that OsbZIP23 participates directly in the process of DS-specific loop formation. However, as OsbZIP23 does not hold the structural basis like CTCF in mammals to handle the loop structure directly, other factors or mechanisms in loop shaping may be involved in an OsbZIP23-dependent manner. A previous report implied OsbZIP23’s role as a H3K4me3 modification mediator in this region [6]. It is also noticeable that the H3K9ac modification in regions like the RAB16 gene cluster was almost impaired in the osbzip23 mutant (Fig. 5n). Moreover, a recent study showed that OsbZIP23 can interact with ADA2 in rice [39]. Thus, we speculate that OsbZIP23 possibly acts as a pioneer factor to mediate the H3K9ac modification and then loop formation in the promoters of DS-responsive genes by directing the histone acetyltransferases GCN5 through the bridging of ADA2 [40]. It is promising that OsbZIP23 and these co-factors can be utilized to optimize 3D genome structure and fine-tune DS-responsive gene expression in rice breeding.

It is also noticeable that not all DS-specific SPRs with OsbZIP23-binding (n = 33) were suppressed in osbzip23 in terms of their CIA values (Fig. 5j). A possible explanation might be that other drought-responsive transcription factors (like the activation form of OsbZIP46) may share the binding sites of OsbZIP23 and perform redundant functions [41]. Thus, OsbZIP23 might not be the sole decisive factor for rice 3D genome remodeling under DS. Together with the recent findings on other 3D genome-associated transcription factors, it is promising that the next-stage studies focusing on the potentially diverse mechanisms behind these transcription factors will give fascinating insight into 3D genome organization in rice.

Conclusions

In summary, our study provides high-resolution 3D genome interaction maps highlighting the dynamic changes in promoter–promoter interactions of genes in response to drought stress and subsequent recovery. For the first time, we discovered promoter-localized enhancer-like regions in rice, i.e., the super-promoter regions that integrate PPI clusters heavily loaded with condition-specific PPIs and genes. Strikingly, we characterized OsbZIP23 as a possible regulator for rice 3D genome re-organization under drought stress. These findings emphasize the 3D genome dynamics as an indispensable layer in transcriptome reprogramming upon water supply extremes. A schematic summary of the results of this study is shown in Fig. 6c.

Methods

Plant materials and drought stress treatment

Wild-type (WT) cultivar MH63 (Oryza sativa L. ssp. indica cv. Minghui 63) or osbzip23 T-DNA mutant seeds were soaked for 3 days until germination. Seedings with similar growth states were then planted in 12 buckets (30 plants per bucket) with equal amounts of soil and water. The buckets were randomly placed in a growth chamber with an air temperature of 32 °C (day) or 29 °C (night), 12 h of light (constant light intensity 15,000 lx) per day, and a relative humidity of 80%. The weight of soil in each bucket was 2.1 kg, and the initial water was 800 mL. After 14 days of normal growth with ample and equal watering, 8 buckets of the 4-leaf stage seedlings were randomly chosen for drought stress treatment while the other 4 buckets continued normal growth. The drought stress treatment was performed by ceasing irrigation for a total of 15 days. The weights of the 8 buckets were adjusted daily with water to keep equal drought progress. During days 11–15, the relative soil water content (RSWC) was kept ~ 5% for each bucket. The RSWC was determined as follows:

$${\text{RSWC}}_{i}=\frac{{m}_{i}-{m}_{dry}}{{m}_{0}-{m}_{dry}} \times 100\%$$

Where ${m}_{i}$ represent the weight of the bucket at day i of the treatment, ${m}_{0}$ is the bucket weight on the day of water cutoff, and ${m}_{dry}$ is the bucket weight with dry soil measured prior to the experiment. After that, half of the 8 buckets were randomly selected for re-watering, while the normally grown and drought-treated samples were harvested. Then, after 4 days of equal and ample re-watering (RSWC = 100%), the samples for recovery were harvested. The timing of each harvest was fixed at 08:00 in the morning (see Additional file 1: Fig. S1a, b for settings on growth conditions and seedling morphology at sampling). For osbzip23, only drought-treated samples were harvested. Samples from the four buckets under the same treatment were combined and regarded as one biological replicate. The above treatment process was repeated twice with different batches of seedlings. The MH63-WT seedlings used in this study were the self-hybrid descendants of the plants used for gapless-genome assembly in a pervious study [42]. The osbzip23 T-DNA insertion mutant is the self-hybrid descendant of the homozygous mutant used in previous studies, and its identity has been approved by Northern blot, Southern blot, and genetic complementation [25, 26].

Whole transcriptome sequencing (RNA-seq) library construction

Total RNA was isolated using the RNeasy Plant Mini Kit according to the manufacturer’s instructions (QIAGEN, 74,904). Two milligrams total RNA was used for library construction following the manufacturer’s recommendations of TruSeq RNA kit (Illumina, 20,020,599). Sequencing was performed using an Illumina HiSeq X Ten system with the 150-bp paired-end strategy (sequencing service was provided by Annoroad Gene Technology, Beijing, China).

Chromatin immunoprecipitation sequencing (ChIP-seq) library construction

ChIP-seq libraries were performed as reported with minor modifications [34]. Rice seedling strips were cross-linked with 1% formaldehyde (SIGMA, F8775) for 10–15 min at room temperature in a vacuum chamber, and 0.2 M glycine was applied to quench the reaction. For each reaction, ~ 0.1 g sample was ground in liquid nitrogen into fine powder and then lysed in 300 μL of Buffer S (50 mM HEPES–KOH (pH7.5), 150 mM NaCl, 1 mM Ethylene Diamine Tetraacetic Acid [EDTA], 1% Triton X-100, 0.1% sodium deoxycholate, 1% SDS) for 10 min at 4 °C. The homogenate was mixed with 1.2 mL of Buffer F (50 mM HEPES–KOH (pH 7.5), 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate). After that, the chromatin in the homogenate was fragmented to 200–600 bp by sonication with a Bioruptor sonication system (Diagencode, UCD-200). Then, the mixture was centrifuged at 12,000 × g for 10 min at 4 °C to recover the chromatin in the supernatant. Twenty microliters of the supernatant was reserved as input, and the remaining was used for ChIP.

To coat antibodies to magnetic beads, 50 μL Dynabeads protein G beads (Life Technologies, 10003D) were balanced with PBST buffer (phosphate-buffered saline (PBS) with 0.1% Tween 20) twice before incubating with 1–2 μg anti-histone H3 acetyl K9 antibody (Abcam, ab10812) or previously approved anti-OsbZIP23 polyclonal antibody [34] for 6–8 h on a rotator at 4 °C. After removing excessive antibodies by washing with PBST twice, the antibody-conjugated beads were mixed with the sheared chromatin and incubated at 4 °C overnight with gentle rotation. The immunoprecipitated chromatin was washed successively with low-salt ChIP buffer (50 mM HEPES–KOH (pH7.5), 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, and 0.1% SDS) (three times), high-salt ChIP buffer (identical to low-salt ChIP buffer expect for the use of 350 mM NaCl) (twice), ChIP wash buffer (10 mM Tris–HCl (pH 8.0), 250 mM LiCl, 0.5% NP-40, 1 mM EDTA, 0.1% sodium deoxycholate) (once), and TE buffer (10 mM Tris–HCl, pH 8.0, and 1 mM EDTA) (twice). The protein–DNA complex was eluted from the magnetic beads with 100 μL Elution Buffer (50 mM Tris–HCl pH 7.5, 10 mM EDTA, 1% SDS) at 65 °C for 15 min and then reverse-cross-linked by overnight incubation with 20 mg/mL proteinase K and 5 M NaCl at 55 °C. ChIP DNA libraries were prepared using the NEBNext® UltraTM II DNA library prep kit for Illumina sequencing (New England BioLabs, E7645). Next, library DNA fragments with 250–650 bp length were preserved using AMPure XP beads (Beckman, A63881). Finally, the libraries were sequenced using the Illumina HiSeq X Ten system with the 150-bp paired-end strategy (Annoroad Gene Technology, Beijing, China).

Long-read ChIA-PET library construction

Rice long-read ChIA-PET libraries were constructed with a modified ChIA-PET method [28]. Tissues were dual cross-linked with 1% formaldehyde for 10 min and 1.5 mM ethylene glycol bis (succinimidylsuccinate) (Thermo Fisher Scientific, 21,565) for 30 min at room temperature. After that, 0.2 M glycine was applied to halt the cross-linking reaction. Then, samples were ground in liquid nitrogen to a fine powder, and approximately 5 g powder was lysed in 2 × volume of buffer S. Next, each sample was applied with 8 × volume of buffer F and sonicated in a Bioruptor sonication system (Diagencode, UCD-200) to yield chromatin fragments with 1–3 kb in length. The chromatin fragments were collected by centrifuging at 12,000 × g for 10 min at 4 °C. Eight hundred microliters protein G with 25 μg anti-Histone H3 acetyl K9 antibody (Abcam, ab10812) was applied to the chromatin fragments and incubated at 4 °C overnight with gentle rotation. The supernatant was discarded, and the beads were sequentially washed with 5 mL low-salt ChIP buffer three times, 5 mL high-salt ChIP buffer twice, 5 mL ChIP wash buffer once, and 5 mL 1 × TE buffer twice. ChIP DNA was used for end-repair and A-tailing using T4 DNA polymerase (Promega, M421F) and Klenow enzyme (NEB, M0212L). Proximity ligation of ChIP DNA was performed using a biotinylated bridge-linker (forward strand: 5′-[5Phos] CGCGATATC/iBIOdT/TATCTGACT-3′, reverse strand: 5′-[5Phos] GTCAGATAAGATATCGCGT-3′). Protein-DNA complex was reverse cross-linked, and the ChIA-PET libraries were prepared using Tn5 transposase (VAHTS, TD501). The libraries were sequenced (2 × 150 bp) using Illumina Hiseq X-Ten (Annoroad Gene Technology, Beijing, China).

RNA-seq data processing

The raw reads were filtered and trimmed by fastp (v0.23.0) [43] with parameters “-q 15” and then mapped to the MH63RS3 assembly [42] using STAR software (v2.7.10a) with the default parameters [44]. Transcripts were assembled and quantified using StringTie (v2.1.7) [45] with the reference-guide mode (“-e”). After removing silent genes with average TPM (transcript per million) < 1 under each condition, the read counts of the rest of the genes were then submitted to R package DEseq2 [46] to determine the differentially expressed genes (DEGs) among the treatments based on its negative binormal model with a threshold of |log₂(fold change)|> 1 and P-adjust (Wald test P-value adjusted by Benjamini–Hochberg method) < 0.05. See Additional file 2: Table S18 for RNA-seq sample summary.

ChIP-seq data processing

The raw reads were filtered and trimmed by fastp (v0.23.0) [43] with parameters “-q 15” and mapped to the MH63RS3 assembly [42] using Bowtie2 (v2.4.5) [47] with parameters “–sensitive, –no-unal”. Uniquely mapped nonredundant reads were preserved by picard (https://github.com/broadinstitute/picard) and were passed to macs2 (v2.2.7.1) [48] for peaks calling (parameters “–call-summits, -f BAMPE, -p 0.01”). Peaks called by the two biological replicates were filtered by an irreproducible discovery rate (IDR) framework (https://github.com/nboley/idr) with a threshold of 0.05 to yield high-fertility peaks. Differentially-covered or shared (with at least 1-bp overlap) H3K9ac peak regions among the conditions were identified using bedtools (v2.29.1) [49]. For quantitative analysis of differential occupation, uniquely mapped nonredundant reads with the union peak regions of all conditions were quantified and normalized (by library size) by R package DiffBind (v3.0) [50]. Then, differential occupation was determined using the DEseq2 method with a threshold of |log₂(fold change)|> 1 and FDR < 0.05. Annotation of the peak-flanking genes was done using the R package ChIPseeker [51] with the − 1000 to + 100 region relative to TSS as the promoter region. The sequencing signal tracks were generated using bamCoverage in the deeptools kit (v3.5.1) [52] with the RPGC (reads per genomic content) normalization and a bin size of 1 bp. The read coverage profiles over genomic regions were calculated and plotted by deeptools with a bin size of 10 bp. See Additional file 2: Table S19 for the ChIP-seq sample summary.

ChIA-PET data processing

The linker filtering, read mapping, redundancy removal, and chromatin interaction identification of the ChIA-PET data was processed using ChIA-PET tools v3 [53], an updated version of the automatic ChIA-PET data processing pipeline previously applied in rice ChIA-PET data processing [15]. We use the union H3K9ac ChIP-seq peak regions under all conditions as seed anchor regions for loop calling (including the osbzip23 samples) so that loops under different conditions can be directly compared based on their anchors. Since the reproducibility between biological replicates was generally good (with PCC > 0.8 based on the intensity of inter-ligation PETs (i-PETs) in the 100-kb window, Additional file 1: Fig. S1d), we perform the loop calling based on the pooled data of the two replicates. We defined the uniquely mapped, nonredundant PETs covering a genome span > 8000 bp as i-PETs. Considering the technical noise and the sequencing depth, candidate chromatin loops with false discovery rate (FDR) < 0.05 and unique i-PET count > 6 were preserved as confident loops and passed to subsequent analysis (see Additional file 2: Table S12 for data summary). To generate a contact matrix, i-PETs and self-ligation PETs (s-PETs) determined by ChIA-PET tools v3 were pooled and submitted to the bedpe2Matrix program in ChIA-PET2 software [54]. The raw matrices were normalized to the one with the smallest data scale among all samples by HiCExplorer [55]. Then, matrices were compared and visualized using the hicCompareMatrices and hicPlotMatrix functions in HiCExplorer, respectively. To determine chromatin-interacting domains (CIDs), we set loops with genomic spans < 200 kb as local-interacting loops. We then identified continuously linked local-interacting loops based on shared anchors by bedtools [49], and the genomic regions covered by at least three continuous local-interacting loops were defined as a CID, representing a potential topologically associated structure of local chromatin interactions.

BL-Hi-C library preparation and data analysis

The BL-Hi-C library was constructed as previously described [17]. The MH63 seedlings were cross-linked with 1% formaldehyde (SIGMA, F8775) for 15 min and halted by 0.2 M glycine, followed by 1.5 mM ethylene glycol-bis (EGS, Thermo Fisher Scientific, 21,565) treatment for 40 min at 25 °C. Isolated rice nuclei were digested with MboI restriction enzyme (NEB, R0147). The resulting chromatin fragments were subjected to end-repair and A-tailing using T4 DNA polymerase (Promega, M421F) and Klenow enzyme (NEB, M0212L), respectively. In situ proximity ligation of DNA fragments was performed using a biotinylated bridge linker (forward strand: 5′-[5Phos] CGCGATATC/iBIOdT/TATCTGACT-3′, reverse strand: 5′-[5Phos] GTCAGATAAGATATCGCGT-3′). Then, the protein-DNA complex was reverse cross-linked, and BL-Hi-C libraries were prepared using Tn5 transposase (VAHTS, TD501). The libraries were sequenced (2 × 150 bp) using an Illumina Hiseq X-Ten (Beijing AnnoRoad Company).

Raw BL-Hi-C data were analyzed using ChIA-PET tools v3 [53]. Chromatin loops were called based on the same seed anchors as ChIA-PET (the pooled H3K9ac ChIP-seq peak regions under all three conditions). An i-PET ( with > 8 kb genomic span) ≥ 6 was set as the threshold for a significant chromatin loop. See Additional file 2: Table S20 for data summary.

Analysis of differential chromatin loop and condition-specific interacting gene pairs

The “pairtopair” function in bedtools (v2.29.1) [49] was used to determine the identical loops between two conditions under the criteria that a loop is regarded as co-occurred only when the pair of anchors match each other. If a loop is unique to a condition in all comparisons, it is termed a condition-specific loop. The loop anchor-flanking genes were annotated by the R package ChIPseeker [51] with the same parameters for ChIP-seq analysis. If the anchor overlapped with multiple genes, only the one with the highest TPM under the corresponding condition was used as a representative. If a pair of genes obtained interacting anchors in their proximal promoter regions (− 1000 to + 100 bp relative to TSS), it was termed a promoter-promoter interacting (PPI) gene pair. If the chromatin loop is condition-specific, then the PPI genes were defined as condition-specific PPI genes.

Hybridization chain reaction (HCR)-based FISH for inter-chromosomal interacting sites

The HCR-FISH assay was conducted following a reported protocol with a few modifications [56]. MH63 seedlings dual cross-linked with 1% formaldehyde and 1.5 mM EGS were chopped in LB01 lysis buffer (15 mM Tris–HCl, pH 7.5, 2 mM Na₂-EDTA, 0.5 mM spermine tetrahydrochloride, 80 mM KCl, 20 mM NaCl, 0.1% Triton X-100) on ice to a fine suspension. The suspension was filtered through the nyclon mesh (pore size 35–50 μm) twice, and then loaded on the surface of equal volume of sucrose buffer (20 mM Tris–HCl pH8.0, 2 mM MgCl₂, 2 mM EDTA, 15 mM 2-ME, 1.7 M sucrose, 0.2% TritonX-100). After centrifugation at 2200 × g for 20 min at 4 °C, the supernatant was discarded and the pellet (nuclei) was resuspended in 1 ml LB01 lysis buffer. The isolated nuclei were spread on microscope slides and permeabilized with PBS buffer containing 0.2% Triton X-100 in Coplin jars. Next, ethanol gradient (70%, 90%, 100%) was applied to the slides for dehydration. After that, RNA was digested with 100 μg/mL RNase solution and the nuclei was post-fixed with 4% formaldehyde. DNA was denatured with denaturing solution (70% deionized formamide, 2 × SSC) at 80 °C for 3 min, and dehydrated in a precooled ethanol gradient immediately. Then, the slides containing nuclei were pre-incubated at 45 °C with the hybridization buffer (50% formamide, 5 × SSC, 9 mM citric acid (pH 6.0), 0.1% Tween 20, 50 μg/mL heparin, 1 × Denhardt’s solution, 10% dextran sulfate) for 30 min followed by an overnight incubation at 45 °C with the hybridization buffer containing 0.2 pmol probe (see Additional file 2: Table S21 for probe sequences). Then excessive probes were removed by incubation with a 75%, 50%, 25% probe wash buffer gradient (50% formamide, 5 × SSC, 9 mM citric acid (pH 6.0), 0.1% Tween-20, and 50 mg/ml heparin) at 45 °C. The slides were pre-amplified in amplification buffer (5 × SSC, 0.1% Tween 20, 10% dextran sulfate) at 25 °C for 30 min. Hairpin solution with fluorescent tags (cy3 and Alexa 488) was then added and incubated overnight at 25 °C. The excessive hairpins were washed with 5 × SSCT buffer (5 × SSC, 0.1% Tween 20) at 25 °C. Finally, slides were stained with 4',6-diamidino-2-phenylindole (DAPI) (Vector, H-1220) before observation under a structured illumination microscope (Nikon, N-SIM). The assay was repeated three times with different batches of seedlings.

Co-expression analysis for drought-responsive genes

Co-expression clustering of the RNA-seq data (with gene with TPM ≥ 1) was performed using the R package Mfuzz based on the Euclidean distance and c-means objective functions [57]. In this soft-clustering process, we pre-defined a series of module members and calculated the membership (ranging from 0 to 1) of each gene to each of the modules based on its TPM values under the three conditions. A gene was categorized into a module if its membership to that module was larger than the other modules. The similarities of the modules were assessed using the “overlap” function in the Mfuzz package. Once a pre-defined module number yielded modules with a coupling coefficient > 0.2, the corresponding categorizations were chosen for further analysis. The co-expression module’s function interpretation was based on the gene ontology (GO) enrichment result generated using AgriGO v2.0 [58]. Significant enrichment of the GO terms was determined by a hypergeometric test with Benjamini–Hochberg adjustment.

Identification of super-promoter regions

We determined the super-promoter regions (SPRs), a concept derived from the super-enhancers, under each condition via a network-based method using the R package tidygraph (https://CRAN.R-project.org/package=tidygraph). We constructed the spatial PPI networks using PPI genes as nodes and chromatin loops connecting the genes as edges. Then, a random walk algorithm was applied to the networks to identify spatial clusters (or subnetworks) containing closely connected PPIs. This was achieved using the group_walktrap function in the tidygraph package [59] with a walking step of 4, and the i-PET counts of PPI loops were used as weight. In each PPI cluster, the nodes (promoter-located chromatin loop anchors) obtaining the highest degree (the number of times a node is connected) were considered as the hub-anchor, and in this regard, the promoter region corresponding to the hub-anchor was termed as hub-promoter regions. Less centralized PPI clusters (degree of the hub-anchor was less than 3) with few PPIs (less than 6 loops and 4 genes) were removed from further analysis. To evaluate the loop-integrating abilities of the hub-promoter regions free from the influences of anchor region size (larger regions generally have higher chances to be connected by more interactions, Additional file 1: Fig. S8a) and the total network scale, we defined contact-integrating affinity (CIA) of a promoter region (i) as follows:

$${\text{CIA}}_{i}=\frac{{degree}_{i} \times {10}^{7}}{\sum_{j=1}^{n}{degree}_{j} \times {anchor}_{i} span in bp}$$

Where degree means the number of times the promoter region i is connected in the PPI cluster, and n is the total number of anchors in this subnetwork. After the CIA values were quantified for the hub-anchors, we ranked them and defined those with CIA values exceeding the 80% quantile as SPRs. If the CIA of an SPR under a condition is 1.5-fold higher than that in the other two conditions, it is regarded as an SPR specific to that condition. TF binding motifs within the DS-specific SPRs (with ± 500 bp spans) were discovered using the AME program in MEME-suit (v5.5.4) with the default settings [60]. In this analysis, non-DS-specific SPRs were set as background, and the 2022 version of the JASPAR core plant motif dataset [61] was used as the candidate motif dataset. The cutoff P-value (determined by Fisher’s exact test) for significant motif enrichment was set to 10⁻². DS-specific SPRs with OsbZIP23 binding site(s) under DS (defined by anti-OsbZIP23 ChIP-seq) were termed OsbZIP23-targeted SPRs in response to drought stress.

Rice protoplast transformation and verification of SPRs by RT–qPCR

To validate the SPR identified by the sequencing data, a CRISPR (clustered regularly interspaced short palindromic repeats)/Cas9-mediated activator was constructed by fusing a dead-Cas9 (dCas9) with the transcription activator VP16 [62]. The CDS of VP16 was amplified from the GAL4–VP16 gene in pFX-E24.2-15R [63] and fused to the ORF (open reading frame) of a previously reported dCas9 [64]. Then, the dCas9–VP16 CDS was cloned into the pRGEB32 vector to replace the original Cas9 (between the BstBI and XbaI sites) [65]. Afterward, the gRNA cassette carrying the RAB16C-targeting spacer sequence was cloned into the pRGEB32 following the developer’s guideline [65]. The working model of the activator is illustrated in Fig. 6b. MH63-WT and osbzip23 protoplast was isolated and transfected as previously described [66]. For each protoplast fraction, 5 µg plasmids carrying the dCas9–VP16 ORF with/without the RAB16C-targeting gRNA were transformed. After overnight incubation at 25 °C, 5 µM ABA was applied to the treatment groups for another 6-h incubation. After that, all samples were harvested, and RNA was extracted using TRIzol™ reagent (Invitrogen™, 15,596,026). The RT-qPCR was performed using a QuantStudio™ 6 Flex system (Applied Biosystems™, 4,485,697) and SYBR™ Green PCR mix (Applied Biosystems™, 4,309,155). The result was quantified using the 2^−ddCt [67] based on the mock MH63-WT sample with dCas9–VP16 targeting at RAB16C (quantified as 1). The assay was repeated three times with different batches of seedlings. Primers used for vector construction and RT–qPCR are listed in Additional file 2: Table S21.

Data visualization

If not specified, data visualization was performed in R (v4.2.2) (https://www.r-project.org/). Genome tracks were visualized using IGV [68] and pyGenomeTracks [69]. Network analysis and visualization were performed by R packages tidygraph (https://CRAN.R-project.org/package=tidygraph) and ggraph (https://CRAN.R-project.org/package=ggraph).

Statistical analysis

If not specified, statistics significances were determined by either two-sided Wilcoxon test (for paired comparisons) or one-way analysis of variance (ANOVA) with Tukey’s Honest significant difference (HSD) test (for comparisons among multiple samples). All statistical analyses were performed in R (v4.2.2).