Skip to main content
Download PDF

Chromatin loops gather targets of upstream regulators together for efficient gene transcription regulation during vernalization in wheat

Genome Biology volume 25, Article number: 306 (2024) Cite this article

Abstract

Background

Plants respond to environmental stimuli by altering gene transcription that is highly related with chromatin status, including histone modification, chromatin accessibility, and three-dimensional chromatin interaction. Vernalization is essential for the transition to reproductive growth for winter wheat. How wheat reshapes its chromatin features, especially chromatin interaction during vernalization, remains unknown.

Results

Combinatory analysis of gene transcription and histone modifications in winter wheat under different vernalization conditions identifies 17,669 differential expressed genes and thousands of differentially enriched peaks of H3K4me3, H3K27me3, and H3K9ac. We find dynamic gene expression across the vernalization process is highly associated with H3K4me3. More importantly, the dynamic H3K4me3- and H3K9ac-associated chromatin-chromatin interactions demonstrate that vernalization leads to increased chromatin interactions and gene activation. Remarkably, spatially distant targets of master regulators like VRN1 and VRT2 are gathered together by chromatin loops to achieve efficient transcription regulation, which is designated as a “shepherd” model. Furthermore, by integrating gene regulatory network for vernalization and natural variation of flowering time, TaZNF10 is identified as a negative regulator for vernalization-related flowering time in wheat.

Conclusions

We reveal dynamic gene transcription network during vernalization and find that the spatially distant genes can be recruited together via chromatin loops associated with active histone mark thus to be more efficiently found and bound by upstream regulator. It provides new insights into understanding vernalization and response to environmental stimuli in wheat and other plants.

Background

Bread wheat (Triticum aestivum) is the most widely cultivated crop worldwide [1], partially benefiting from its capability to overcome cold temperatures in winter. Winter wheat requires a particular duration of cold days before flowering in warm spring, which is called vernalization, the response of temperate plants to prolonged low temperature [2,3,4,5,6]. Vernalization is critical for winter wheat, which has a longer lifespan and generally a higher yield and better grain quality than spring wheat [7,8,9].

Four genes, VERNALIZATION1 (VRN1), VRN2, VRN3, and VRN4, have been shown to control vernalization-mediated flowering in wheat [3, 10]. VRN1 encodes a MADS-box transcription factor and is the homolog of Arabidopsis floral meristem identity gene APETALA1 (AP1). The transcription of VRN1 is induced by cold temperatures and accumulates after sufficient cold exposure to promote flowering in winter wheat [11]. VRN2, whose transcription is repressed by vernalization, encodes a zinc-finger protein and acts as a flowering repressor to delay wheat flowering [12]. VRN3 is orthologous to Arabidopsis FLOWERING LOCUS T (FT) and functions as florigen in wheat [13]. VRN4 is an additional variant of VRN1 that exists in an ancient wheat in South Asia, and it is expressed to accelerate wheat flowering without the requirement of vernalization [14]. VRN1 binds to the promoters of VRN2 and VRN3 to repress the transcription of VRN2 while activating the expression of VRN3 in both wheat and barley [15, 16]. Meanwhile, VRN1 downregulates the expression of C-REPEAT BINDING FACTOR (CBF) genes, which play critical roles in freezing tolerance [16, 17]. Therefore, VRN1 plays a central role in cold resistance and vernalization response in wheat [3].

Transcription induction of vernalization-related regulators is under epigenetic regulation in both dicots and monocots [3, 18,19,20,21]. In Arabidopsis, the flowering repressor, FLOWERING LOCUS C (FLC), remains active before winter but becomes repressed by vernalization during winter [22, 23]. The silencing of FLC transcription is associated with the transition from active to repressed histone marks at the FLC locus, represented by H3K4me3 and H3K27me3, respectively [24,25,26]. In temperate grasses, including wheat, barley, and Brachypodium distachyon, transcription of VRN1, the key regulator during vernalization, is associated mainly with dynamic change of histone modification [27,28,29,30]. H3K27me3 is enriched on the VRN1 gene to inhibit its transcription before vernalization, whereas H3K4me3 is increased gradually during the prolonged cold exposure, resulting the elevated expression of VRN1 [27, 28, 30,31,32]. Besides the histone modification, the function of VRN1 is orchestrated by RNA-binding protein TaGRP2 targeting its pre-mRNA and MADS-box protein TaVRT2 binding to its promoter [31, 33, 34]. Interestingly, an alternative splicing form of VRN1, VAS, promotes its transcription by facilitating the docking of TaRF2b-TaRF2a complex at the VRN1 promoter [35]. More than the key regulators, plants respond to environmental stimulates by systematically reshaping the gene transcription network under epigenetic control [36,37,38]. The genome-wide gene transcription dynamics under different vernalization treatments in Brachypodium distachyon had been investigated, and the role of H3K4me3/H3K27me3 modifications have been demonstrated [30, 39]. Moreover, increasing evidence suggests that the three-dimensional (3D) chromatin-chromatin interaction plays a critical role in gene transcription regulation in both wheat and other species [40,41,42,43,44,45]. Nevertheless, how vernalization reshapes the genome-wide atlas of gene transcription with changed chromatin status, especially the 3D chromatin interaction is still largely unknown in wheat.

In this study, we conducted multi-omics analysis by integrating transcriptome, epigenome, and H3K4me3- and H3K9ac-associated 3D-chromatin interactome in winter wheat under different vernalized conditions to construct vernalization-mediated gene regulatory network and identify key genes that determine vernalization response in wheat. Our results not only deepen the understanding of gene transcription induced by vernalization but also provide novel targets for future adaptation breeding in wheat as well as in other temperate crops.

Results

A widely planted winter wheat cultivar, Aikang 58 (AK58) [46], was used to investigate the impact of vernalization on gene transcription in wheat. To estimate its requirement of vernalization, we measured the flowering time of AK58 plants grown under different conditions. The result showed that non-vernalized (V0) AK58 does not flower until more than 130 days after sowing. Twenty-eight-day vernalization (V28N) is sufficient to induce flowering while 14-day vernalization (V14N) leads to later flowering than 28-day vernalization (Fig. 1a, b). We then profiled transcriptome at the time points of V0 (not vernalized), V28 (vernalization at 4 °C for 28 days), and V28N6 (vernalization for 28 days followed by normal growth for 6 days), which mimics the three major phases during vernalization process in the nature (Fig. 1c). Accordingly, epigenome was generated by chromatin immune-precipitation followed by high-throughput sequencing (ChIP-seq) using antibodies against H3K4me3, H3K27me3, and H3K9ac, which are marks for active gene transcription, repressed gene transcription, and gene regulatory elements (proximal and distal) to activate gene transcription in wheat [47, 48], respectively, and by assay for transposase-accessible chromatin using sequencing (ATAC-seq) [49] and H3K4me3- and H3K9ac-associated chromatin interaction by paired-end tag sequencing (ChIA-PET) [50] (Additional file 1: Fig. S1). In total, 39 datasets were generated (Additional file 2: Table S1).

Fig. 1
figure 1

Global transcriptome dynamics during vernalization in wheat. a 28-day vernalization is more promising to initiate flowering for winter wheat cultivar, Aikang 58 (AK58), compared with 14-day vernalization. b Bar graph showing flowering time of AK58 without vernalization (V0) and with vernalization for 14 days (V14) and 28 days (V28), respectively. c Diagram showing the treatment and sampling strategy. AK58 seedlings with V0 (14 days old seedling, without vernalization), V28 (at normal temperature for 9 days and at 4 °C for 28 days), and V28N6 (normal temperature for 3 days and at 4 °C for 28 days followed by normal condition for 6 days) treatment were harvested at three-leaf stage. d PCA analysis of the transcriptome datasets. e Venn diagram showing differentially expressed genes among V0, V28, and V28N6 (p-value < 0.05, fold change > 2). f Heatmap showing seven clusters of differentially expressed genes across three treatments. g GO analysis showing the top 5 significantly enriched terms in each gene cluster. The length of light blue rectangle represents the significance of each pathway. h Sankey plot showing altered asymmetric expression pattern of homoeologs across three treatments

Full size image

Genome-wide gene transcription dynamics across the vernalization process in wheat

The principal component analysis (PCA) of RNA-seq data (Additional file 2: Table S2) uncovers significant transcriptome dynamics among treatments (Fig. 1d). We identified 13,187 (5975 upregulated and 7212 downregulated) differentially expressed genes (DEGs, fold change > 2, q-value < 0.05) when comparing the gene transcription profile in V28 samples with V0 samples (V0 vs V28), and 12,656 DEGS (7864 upregulated and 4792 downregulated) were isolated in the comparison between samples from V28N6 and V28 (V28 vs V28N6) treatments (Fig. 1e; Additional file 1: Fig. S2a, b and Additional file 2: Table S3). Furthermore, there are 5453 (3537 upregulated and 1916 downregulated) genes differentially expressed in V28N6 plants, compared with non-vernalized plants (V0 vs V28N6) (Fig. 1e; Additional file 1: Fig. S2c). Overall, 17,669 genes were dynamically expressed across the three treatments, and they were grouped into seven clusters. Among them, genes in cluster 1 and cluster 2 showed the highest expression level before vernalization (V0) and significantly reduced transcription level after prolonged cold exposure; genes in cluster 3 and cluster 4 showed elevated expression after vernalization (V28); genes in cluster 5 showed increased expression after vernalization and remain high level when the plants were transferred to normal condition, and genes in cluster 6–7 showed highest expression level in V28N6 plants (Fig. 1f; Additional file 1: Fig. S2d). Detailed Gene Ontology (GO) analysis revealed that gene related with photosynthesis, protein, and carbohydrate metabolism are active in V0 plants, while genes responsible for RNA metabolism is dominant in V28 plants and genes responding to abiotic stimulus are activated after the plants were returned to normal temperature (Fig. 1g). These well-known vernalization and flowering-related genes, including VRN1, showed a dynamic expression pattern across different vernalization conditions (Fig. 1f).

Asymmetric gene transcription contributes to plant growth and trait formation in polyploidy species [48, 51]. We investigated the asymmetric gene transcription among subgenomes during vernalization response in wheat. We identified 6218 homoeolog triads (composed of copies from A, B, and D subgenomes) from the 17,669 DEGs (Additional file 2: Table S4). We noticed that the asymmetric gene expression pattern is similar among the three conditions (Additional file 1: Fig. S2e). More than 60% of triads are balanced triads, expressing the homoeologs from A, B, and D subgenomes equally. In the unbalanced triads, the proportions of A suppressed and B suppressed are significantly higher than others (D suppressed, A/B/D dominant). Interestingly, 28.67% (1783/6218) triads changed their asymmetric gene expression pattern, and the A, B, and D homoeologs of 30.01% (1866/6218) showed dynamic gene transcription (not necessarily changed their asymmetric gene expression pattern) across the three vernalization conditions (Fig. 1h; Additional file 2: Table S4), indicating potential function differentiation of subgenome homoeologs during vernalization process in wheat. For example, A subgenome member of ODDSOC2, Ppd1, VRT2, and VRN1 and B subgenome member of TaOGT1, TaTOE1, and WPCL1 [11, 34, 52,53,54,55,56] showed distinct expression patterns compared with the other two corresponding homoeologs during vernalization (Additional file 1: Fig. S2f).

Active marks of H3K4me3 and H3K9ac are more comprehensively and highly correlated with changing gene transcription than H3K27me3 during vernalization in wheat

It is reported that changes in histone modification level on key regulators, like VRN1, are one of the major mechanisms in vernalization response in Triticeae grass [27, 28, 30]. To reveal the genome-wide histone modification dynamics during vernalization in wheat, we mapped loci modified by H3K4me3, H3K9ac, or H3K27me3 and the genomic regions with accessible chromatin (Fig. 2a; Additional file 1: Fig. S3 and Additional file 2: Table S5). As expected, the H3K4me3, H3K9ac, and chromatin accessibility showed positive correlation with gene transcription level revealed by RNA-seq, whereas the repressive mark H3K27me3 is negatively correlated with all the active mark and also with the gene transcription level (Additional file 1: Fig. S4a). Genes with H3K4me3, H3K9ac, or open chromatin showed significant higher expression level than genes without these modifications, while the expression of genes with H3K27me3 modification is repressed (Additional file 1: Fig. S4b).

Fig. 2
figure 2

Dynamic change of histone modifications during vernalization in wheat. a Circular diagram showing the genome-wide distribution of histone modification and chromatin accessibility in V0 plants. b Principal component analysis of H3K4me3, H3K9ac, and H3K27me3 among different treatments. c Box plots showing correlation between alteration of gene expression and differentially enriched peaks (DEPs) in its proximal region (promoter and gene body, left) and distal region (intergenic, right). d Heatmap showing DEPs of H3K4me3, H3K9ac, and H3K27me3 on the seven clusters of genes that are dynamically expressed. The red or purple color represents the peak intensity of the corresponding histone modification. e Dynamic change of histone modifications on VRN1, GRP2, and VRT2 during the vernalization process. f Heatmaps showing expression level and histone modification dynamics of published vernalization/flowering-related genes

Full size image

PCA analysis showed that the enrichment of both the active (H3K4me3 and H3K9ac) and the repressive (H3K27me3) marks is dynamically changing across different conditions (Fig. 2b). The fold change of DEGs between two adjacent stages (from V0 to V28 then to V28N6) is positively associated with H3K4me3 (R = 0.74 and 0.62) and H3K9ac (R = 0.56 and 0.46) and negatively correlated with H3K27me3 (R = − 0.25 and − 0.15) in the comparisons of V0 vs V28 and V28 vs V28N6 (Additional file 1: Fig. S5a-c), suggesting the critical role of histone modification during vernalization response, which is similar with the situation in Brachypodium distachyon [30]. We further identified peaks that are dynamically enriched across the three conditions, and we found 20,155, 38,431, and 6317 dynamically enriched peaks (DEPs) for H3K4me3, H3K9ac, and H3K27me3, respectively (Additional file 2: Table S6-8). These DEPs were clustered into seven modules for each epigenetic mark (Additional file 1: Fig. S6). To further reveal the influence of DEPs located at proximal (3.5 kb upstream of transcription start site [TSS] to transcription end site [TES], where gene promoter and gene body are included) and distal (other genomic regions out of proximal genic region) genomic regions on gene transcription activity during wheat vernalization, we further investigated the expression level of their corresponding genes. We noticed that the dynamic change of gene expression level showed higher correlation with their proximal DEPs than the distal DEPs for all of H3K4me3, H3K9ac, and H3K27me3 (Fig. 2c). We then analyzed the dynamic changes of H3K4me3, H3K9ac, and H3K27me3 solely on DEGs identified across the three vernalization treatments (Fig. 1f). Interestingly, we found that proximal DEPs of H3K4me3 and H3K9ac are highly correlated with the expression level of DEGs (Fig. 2d). While in the case of H3K27me3, only the expression of a small proportion of genes showed negative correlation with H3K27me3 (Fig. 2d–f). We focused on the 4113 H3K27me3 DEPs locating at the proximal regions of 4098 genes among conditions of V0, V28, and V28N6, and only 564 genes are differentially expressed (Fig. 2d). When we looked at these 564 differentially expressed genes, the correlation between changed gene expression and H3K27me3 modification is in generally not significant. Only the expression of genes in cluster 4 containing VRN1 is correlated with H3K27me3 change (Fig. 2d). Interestingly, wheat vernalization/flowering-related genes, such as VRN1, TaGRP2, and VRT2, are undergoing dynamic histone modifications during vernalization (Fig. 2e, f). These results demonstrated that vernalization-mediated transcriptional dynamics were associated with alteration of histone modifications, especially those in proximal regions. Moreover, it is active mark (such as H3K4me3, H3K9ac) than suppressing mark (H3K27me3) that is highly associated with the gene transcription alteration during wheat vernalization process, which may be different from the situation in Arabidopsis [19,20,21].

Chromatin loops associated with H3K4me3 and H3K9ac are highly associated with gene activation during vernalization

The above analysis suggested that the change of H3K4me3 shows the highest correlation with alteration of gene expression during wheat vernalization than H3K9ac and H3K27me3 (Fig. 2d, f; Additional file 1: Fig. S4, 5). Long-distance chromatin-chromatin interaction participates in gene transcription regulation [40,41,42]. To investigate the impact of chromatin-chromatin interaction on gene transcription during wheat vernalization process, we first performed high-resolution chromatin interaction analysis using paired-end tag sequencing (ChIA-PET) to analyze the genome-wide chromatin-chromatin interactions associated with the active mark H3K4me3 before and after vernalization (V0 and V28N6), respectively (Fig. 3a, b; Additional file 2: Table S9). We found the highest number of interactions within the DD subgenome and the second highest number of interactions between the AA and DD subgenomes in both of the cases (Fig. 3a–d). The short arm of chromosomal 1B showed a much lower interaction frequency with other genomic regions (Fig. 3a, b). This is in line with the fact that AK58 contains the 1B/1R translocation, in which chromosome 1R is introduced from rye [43, 44, 46].

Fig. 3
figure 3

Altered chromatin interactions associated with H3K4me3 during wheat vernalization. a, b ChIA-PET heatmaps showing chromatin interactions associated with histone modification H3K4me3 in V0 (a) and V28N6 (b) plants. The dashed boxes represent interactions within A, B, or D subgenome. c, d Boxplots showing interaction frequency within different subgenomes associated with H3K4me3 modification in V0 (c) and V28N6 (d) plants. e Models explaining genes and genomic loci with H3K4me3 modification and chromatin interaction. BP gene, genes with histone modification but without loops. Anchor gene, genes connected by chromatin loops. Rest gene, gene without H3K4me3 modification. PPI only, anchor genes only connected by PPI. PDI only, anchor gene with PDI. PPI & PDI, anchor genes with both PPI and PDI loops. f Statistics of basal promoter (BP), PPI interaction and PDI interaction identified in V0 (left) and V28N6 (right) plants. g Pie diagrams showing statistics of BP genes, anchor genes and rest genes in V0 and V28N6 plants. h Expression level of genes with different types of chromatin interactions. ***P < 0.001. i Boxplots showing the expression level of BP genes (basal), anchor genes with different degrees, and rest genes without H3K4me3 modification. j Column diagram showing the statistics of PPI and PDI loops in the five types of dynamic chromatin interaction (stable, up, down, appear, disappear) before and after vernalization. k, l Heatmaps showing the elevated expression level of anchor genes with increased interaction degree of PPI (k) and PDI (l) loops. m Browser screenshots showing the H3K4me3-associated chromatin interactions and gene expression level surrounding the genomic regions of the three homoeologs of VRN1 and VRT2 in V0 and V28N6 plants. inter, interchromosomal chromatin interactions. intra, intrachromosomal chromatin interactions. n Validation of the changed chromatin interactions at VRN1 locus during vernalization using DNA fluorescence in situ hybridization (FISH). Anchor regions are stained green (VRN1-A), red (VRN1-B), and carmine (VRN1-D or Chr5A:588380404), respectively. The separated red, green, and carmine channels were provided in Additional file 1: Fig. S8. Bar = 5 μm

Full size image

Depending on the location of peaks that are connected by a chromatin loop, chromatin loops were classified into proximal-proximal interaction (PPI) and proximal–distal interaction (PDI). A gene promoter enriched with H3K4me3 but has no interaction via chromatin loop with other genomic regions is defined as a basal promoter (BP). Genes with loop connected are anchor genes, including genes with only PPI, with only PDI, or with both PPI and PDI, and genes without loop linked but with H3K4me3 are BP genes [41, 57] (Fig. 3e). We identified 82,370 and 90,332 H3K4me3-associated loops in wheat before and after vernalization (Additional file 2: Table S10-11). There are 54,800 PPIs and 27,570 PDIs identified for 45,983 anchor genes and 14,286 potential distal regulatory regions in V0 plants, and 60,169 PPIs and 30,163 PDIs connecting 46,928 anchor genes and 14,803 distal regulatory genomic regions were mapped in V28N6 plants (Fig. 3f, g). The anchor genes with both PPI and PDI loops showed the highest transcription level, followed by the genes connected by PPI or PDI loops only, the anchor genes, and the BP genes in order. And both the anchor genes and BP genes are much more actively transcribed than rest genes without H3K4me3 peak (Fig. 3h). Moreover, a positive correlation between gene transcription level with chromatin interaction degree (interaction frequency) was observed in both V0 and V28N6 samples (Fig. 3i).

We wondered how the alteration of H3K4me3-associated interaction contributes to the change of gene transcription induced by vernalization. Based on the existence and number of loops, we divided anchor genes into five types, designated as “stable” (genes with the same interaction degree after vernalization), “up” (genes showed increased degree), “down” (genes showed decreased degree after vernalization), “appear” (genes had no loops before vernalization but loops appeared after vernalization), and “disappear” (genes had loops before vernalization but the loops disappeared after vernalization). The results showed that each type of the anchor genes was distributed roughly equally among the A, B, and D subgenomes with slightly less “stable” and “up” types in the B subgenome (Additional file 1: Fig. S7). The number of genes showing the “up” or “down” type is higher than the “disappear” or “appear” type after vernalization. Interestingly, most of the altered loops in the “up” and “down” types were PPI loops, while the proportion of PDI loops is predominant in the groups of “appear” or “disappear” after vernalization (Fig. 3j). These results indicated that both the quantitative and qualitative complexities of chromatin interactions contribute to the dynamic change of gene transcription triggered by vernalization, and the quantitative change is more frequently seen. We then asked which genes are affected by the changed PPI and PDI. We focused on those genes that increased the interaction degree after vernalization (“up” and “appear” types) since we are more interested in the genes that are activated by the vernalization process. Not surprisingly, these genes showed increased expression levels in V28N6 plants, including the published vernalization and flowering-related genes (Fig. 3k, l), especially the three homoeologs of VRN1 and VRT2 (Fig. 3m). To further confirm the chromatin interaction at the VRN1-A locus, we verified the interaction between VRN1-A locus and VRN1-B/D as well as an intra-chromosome interaction between VRN1-A and TraesCS5A02G392000 (probe position: Chr 5A:588380404). It was clearly shown that VRN1-A marked by green florescence got together with VRN1-B and VRN1-D marked by red and carmine florescence, respectively, after vernalization. Furthermore, VRN1-A and the genomic region of TraesCS5A02G392000 were linked together after vernalization (Fig. 3n; Additional file 1: Fig. S8).

Except for H3K4me3, H3K9ac also showed close correlation with dynamic gene transcription during vernalization. We then uncovered the chromatin interactions associated with H3K9ac (Additional file 2: Table S9, 12, 13). It revealed that the interaction frequency within DD is the highest, and the patterns of interaction frequency before and after vernalization are similar (Additional file 1: Fig. S9a, b). There are 44,170 and 45,010 anchor genes linked by H3K9ac-associated chromatin interactions in V0 and V28N6 plants respectively (Additional file 1: Fig. S9c). Not surprisingly, gene with both PPI and PDI showed the highest gene expression level, and the more interactions one gene was involved, the higher expression of the gene can be detected (Additional file 1: Fig. S8d, e). We also found that the chromatin interactions are dynamically changing before and after vernalization, and the interaction degree of vernalization-related genes including VRN1 and VRT2 is significantly increased (Additional file 1: Fig. S9f-h). Interestingly, all the above results are highly comparable with the results obtains from ChIA-PET data using H3K4me3 antibody (Fig. 3c–m), which is line with the fact that the anchor genes marked by H3K4me3 and H3K9ac are largely overlapped (Additional file 1: Fig. S9i). These results supported the pivotal role of H3K4me3- and H3K9ac-associated 3D chromatin interaction on gene transcription activation during vernalization in wheat. Since H3K4me3 and H3K9ac seems to play very comparable role in indicating active gene transcription, we mainly used the H3K4me3-related data for further analysis.

H3K4me3-associated chromatin loops gather targets bound by upstream regulator together for efficient gene activation

In animal cells, the CCCTC-binding factor (CTCF) mediates the formation of topologically associated domains for efficient gene transcription regulation [58]. We were curious how the H3K4me3-associated chromatin interactions activate gene transcription during vernalization. We first tried to figure out the possible upstream regulator (especially transcription factors, TFs) of the anchor genes that showed increased interaction degree and enhanced expression level (Fig. 3k, l) with the information of TF footprints revealed by ATAC-seq in V28N6 plants (Additional file 2: Table S5). We constructed a gene regulatory network, and strikingly, VRN1 appeared to be the core regulator with the highest number of outdegree (to regulate others) to regulate a group of anchor genes (Fig. 4a). This result indicated that the targets of the upstream regulator were gathered together by chromatin loops. To further confirm this hypothesis, we mapped the binding targets of VRN1 by DNA affinity purification sequencing (DAP-seq) [59] in combination with the information of TF footprints by ATAC-seq in V28N6. A total of 18,033 possible binding targets were identified, and the classical CArG box was identified as the binding motif of VRN1 (Fig. 4b, c). Intriguingly, we found that most of the predicted VRN1 target genes were anchor genes that are involved in chromatin looping regardless of transcriptional upregulated (81.1%, 3866/4765) or downregulated genes (86.8%, 3685/4245) (Fig. 4d). There are more VRN1 targets (2977 vs 1836) showing increased interaction degree (up or appear) than decreased degree genes (down or disappear) in both transcriptional upregulated (1459 vs 987) and downregulated anchor genes (1518 vs 849) (Fig. 4e). In addition to VRN1, we also tested whether other hub regulators involved in vernalization behavers in the same manner. We did not detect significant chromatin loops at both VRN2 and VRN3 loci, which may be due to the fact that H3K4me3 is hardly seen but H3K27me3 is highly enriched at these loci (Additional file 1: Fig. S10). We focused on VRT2 that had been reported to interact with VRN1 to play key role in vernalization [34]. We applied the publicly available DAP-seq data for VRT2 (Additional file 1: Fig. S11a, b) [60] and analyzed the dynamic change of chromatin interactions of the targets. The result showed that chromatin interaction degree was increased at VRT2 locus after vernalization (Fig. 3m). Moreover, similar with the case of VRN1-A, a much higher proportion (2880 vs 671; 2730 vs 440) of VRT2 targets are involved in chromatin loops and the degrees are dynamically changing before and after vernalization (Additional file 1: Fig. S11c, d). These results support the notion that chromatin loops gathered vernalization-related genes together for efficient co-transcription of these genes regulated by master regulators, such as VRN1 and VRT2.

Fig. 4
figure 4

Efficient gene transcription achieved by H3K4m3-associated chromatin loops during vernalization. a VRN1 is the core upstream regulator of genes from Fig. 3k, l, with the highest outdegree in the network. b Read density of VRN1 DAP-seq data surrounding the transcription start site (TSS). c CArG box was identified from VRN1 DAP-seq data. d Column diagram showing the number of putative VRN1 target genes connected by loops during vernalization process both in transcriptional upregulated (left) and downregulated (right) gene sets. e Number of loop-connected VRN1 target genes with different altered types (stable, up, down, appear, disappear) of interaction degree both in transcriptional upregulated (left) and downregulated (right) gene sets. f, g Vernalization-related regulatory networks based on TF footprint and chromatin interactions in V0 (f) and V28N6 (g) plants. The light green line represents specific link in V0 plants; the light purple line indicates specific link in V28N6 plants. h Validation of the chromatin interactions at ZNF10 loci after vernalization by DNA FISH. Anchor regions are stained carmine (ZNF10-5A or ZNF10-5B) and green (Chr5A:593761121 or Chr5B:581709899 respectively. Bar = 5 μm

Full size image

Hence, the systematic impact of dynamic chromatin interaction on vernalization response was investigated. The regulatory networks constructed by combining TF footprint with chromatin interactions contain 9324 links in the V0 network connecting 6779 nodes (4396 DEGs and 2383 distal sites) and 9871 links in the V28N6 network connecting 6248 nodes (4797 DEGs and 1451 distal sites) (Fig. 4f, g). Significant differences in V0 and V28N6 regulatory networks with limited overlaps can be observed (Additional file 1: Fig. S12a, b). We further identified 364 TFs with dynamic gene expression as upstream regulators to drive expression change of downstream genes before and after vernalization (Additional file 2: Table S14). Among them, we found key regulators for vernalization, such as VRN1, VRT2, TaFDL2, WAG1, and WAG2, and genes including TaAGL17-B1, TaAGL17-B3, and TaSTK-D1 that belong to MADS-box family, which have been reported to play an important role in plants flowering or flower development [61]. Notably, a group of C2H2-type zinc-finger (ZNF) genes is targeted by VRN1 and other proteins. More ZNF genes are gathered together by chromatin loops after vernalization (Fig. 4g) than in V0 plants (Fig. 4f), and this result was further confirmed by cellular evidences from both ZNF10-5A and ZNF10-5B with their corresponding interacting genomic locus (Fig. 4h). This result again supported the hypothesis that target genes of upstream key regulators are recruited by long-distance chromatin loops for efficient regulation.

The ZNF genes are involved in vernalization-related flowering time control in wheat

Subsequently, by examining the role of ZNF genes, we investigated the biological significance of the long-distance loop mediated gene activation. We identified 39 ZNF genes showing differential expression during wheat vernalization (Fig. 5a). Notably, the change in H3K4me3 level at 25 loci is correlated with the dynamic change of the target gene expression, while for the rest 14 ZNFs, the correlation between dynamic change of gene expression and H3K4me3 level is pretty weak (Additional file 1: Fig. S13). Most of these ZNF genes were suppressed by vernalization treatment, indicating the negative role of ZNF genes in wheat vernalization. To further support the function of ZNF genes in wheat flowering time control, we conducted haplotype analysis of all the differentially expressed ZNFs by considering the natural variation of flowering time from the published data [62]. We found that there are 12 ZNFs showing significant differences between haplotypes (Fig. 5b), indicating the putative function of these ZNF genes.

Fig. 5
figure 5

Validating the function of ZNF genes in flowering time control. a Heatmap showing differently expressed ZNF genes during vernalization process. b Genetic effect by haplotype analysis for ZNF genes using wheat flowering time data. c Positive TaZNF10-5B-OE transgenic plants showing delayed flowering compared with wildtype (WT) and negative transgenic plants. Bar = 10 cm. d Boxplot showing the statistics of flowering time of positive and negative TaZNF10-5B-OE transgenic plants and the wildtype. e Homozygous TaZNF10-5B-KO plants showing earlier flowering compared with wildtype plants. Bar = 10 cm. f Boxplot showing the statistics of flowering time data collected from TaZNF10-5B-KO and wildtype plants. g Boxplots showing the statistics of flowering time of TaZNF10-5B-OE and TaZNF10-5B-KO and negative plants without vernalization (V0) and with vernalization for 14 days (V14) and 28 days (V28), respectively

Full size image

To determine the function of ZNFs, we picked TaZNF10 (Fig. 5a, b) as target gene to generate transgenic lines for further analysis. Among the 12 ZNFs with significant difference among haplotypes, TaZNF10 showed less degree of significance between haplotypes than most of ZNF genes. If this gene shows clear function in flowering time control, other ZNFs should have even higher chance to regulate flowering time in wheat. Notably, the positive TaZNF10-5B-OE lines showed 10–13 days’ delay in flowering time than wild type and negative control (Fig. 5c, d). We also tried to produce a knockout mutant. Unfortunately, we failed to create a mutant with all of the three homoeologs in A, B, and D subgenomes knocked out but only obtained a mutant for TaZNF10-5B (TraesCS5B02G406000). Nevertheless, slight but significant difference could be seen between the earlier flowering homozygous mutant (TaZNF10-5B-KO, with 15-bp deletion in the exon to result 5 amino acid loss in the conserved C2H2 domain [Additional file 1: Fig. S14]) and the wild type (Fig. 5e, f). Furthermore, to demonstrate that the TaZNF10 is involved in vernalization, we treated all of the transgene negative, knock-out mutant, and the overexpression lines without vernalization (V0) or with vernalization at 4 °C for 2 weeks (V14) or 4 weeks (V28), respectively. The results showed that without vernalization, the knock-out mutant flowered about 1 week earlier than the transgene negative plants, and the overexpressed plants did not flower 160 days after sowing. With 2 weeks’ vernalization treatment, there are only less than 5 days’ difference between the knockout mutant and the negative control, and the overexpression lines started to flower at around 140 days after sowing. When treated by 4 weeks’ cold, a significant difference in flowering time among knock-out mutant, overexpressed line, and the negative control can be still observed, but the difference gets much smaller than the situation under 2 weeks’ vernalization treatment (Fig. 5g). All these results approved that TaZNF10 is involved in the vernalization-related flowering time control in wheat, which indicated that the recruitment of targets via chromatin loops is of great importance during wheat vernalization response and flowering time control.

Discussion

Winter wheat is always sown in autumn, and the seedlings stop to grow (or grow very slowly) in cold winter. When the spring comes with the warm temperature, the seedlings start to develop, and soon, the transition to reproductive growth happens for future spike and seed development. During this life span, vernalization is crucial because it promotes the transition from the vegetative to the reproductive phase. Requirement of long-time cold exposure (vernalization) before flowering ensures proper flowering time in warm spring [3, 9, 63].

Several genes have been identified to participate in vernalization response and flowering time control in wheat [9, 64]. To study gene transcription regulation during vernalization response in wheat, we mimic the overwinter growth of wheat and picked the three representative phases for investigation (Fig. 1c; Additional file 1: Fig. S1). Plants respond to environmental stimulus by changing gene transcription promptly. It would make sense to trace the dynamic trajectory of gene transcription in a time-series manner starting from the first few minutes after cold treatment begins. However, cold exposure triggers two different types of responses in plants that need to pass winter: cold adaptation, a rapid response to short-term cold and vernalization, a long-term response to long-term coldness [21]. Plants establish freezing tolerance through cold acclimation after a short period of exposure to low temperatures [65], but it takes weeks or even months to reach the required vernalized state. Only prolonged exposure to cold promotes the plants to enter the reproductive growth [66]. It is challenging to distinguish vernalization from cold response, a response to stress condition, to define the starting point of vernalization. Thus, our study here focuses on the impact of sufficient vernalization on gene transcription and the underlying regulatory mechanism in wheat. Although the initial gene transcription events at the beginning of vernalization might be missed, the information required for maintaining the sufficient vernalized status can be captured with the sampling strategy we applied. Of course, more detailed studies on the difference and connection between prompt cold response and vernalization would significantly enhance our understanding of plants responding to cold temperature.

We found that in total, 17,669 genes significantly altered their expression during vernalization, including a set of published genes involved in vernalization, such as VRN1, VRN2, VRT2, TaGRP2, and ODDSOC2 (Fig. 1e, f; Additional file 1: Fig. S2a-d). Interestingly, the A, B, and D subgenomes contribute unequally during the process. For instance, the expression of VRN1-A is higher than VRN1-B and VRN1-D (Additional file 1: Fig. S2f). The asymmetric gene expression contributes to function divergence among homoeologs, which ultimately leads to enhanced environmental adaptability [48]. Histone modification has been proven to have a fundamental impact on gene transcription [67]. During vernalization response, the two functionally antagonistic histone marks, H3K4me3 and H3K27me3, are highly correlated with transcription regulation on key vernalization regulators like FLC and VRN1 [24,25,26]. However, although we saw the negative correlation between H3K27me3 and gene transcription (Additional file 1: Fig. S4,5c), the genome-wide alteration of gene expression during vernalization in wheat showed strong correlation with active markers of especially H3K4me3 and H3K9ac but very weak correlation with repressive mark H3K27me3 (Fig. 2d–f), which has been shown to play critical role in vernalization-mediated gene suppression in Arabidopsis. These findings suggest that gene activation is the major mechanism during wheat vernalization response.

Gene transcription is known to be influenced by chromatin structure [40,41,42]. We uncovered H3K4me3- and H3K9ac-associated chromatin interactions before and after vernalization using a long-read ChIA-PET strategy in wheat for the first time (Fig. 3a, b; Additional file 1: Fig. S9). Different from high-throughput chromosome conformation capture (Hi-C), which probes physical interactions among all genomic loci [68, 69], ChIA-PET is known to be mediated by specific protein tethering linear genomic elements to high-order chromatin architecture [50, 70]. The three-dimensional interactome showed that chromatin interactions between gene proximal regions is a more potent contributor to gene expression than interactions between gene proximal and distal regions (Fig. 3h; Additional file 1: Fig. S9d). Interestingly, histone modification in the proximal regions (intragenic and promoter) showed stronger correlation with the altered gene expression than that in distal regions (intergenic) (Fig. 2c). These findings suggest that regulation on proximal elements may be the main driving force in response to vernalization. Gene transcription network is usually leaded by master regulators targeting a group of downstream targets [71]. To achieve fast and efficient binding and gene transcription regulation, one ideal situation is that the targets of an upstream regulator are spatially closed thus the protein can easily find the targets, like a shepherd managing a flock of sheep with sheepfold. We found that binding targets of upstream regulators, represented by VRN1 and VRT2, were gathered together by H3K4me3-associated chromatin loops (Fig. 4d; Additional file 1: Fig. S11). Here, we proposed the shepherd model to describe the efficient gene transcription activation upon vernalization vis chromatin looping. Compared with the transcription factory model that proposes efficient transcription of transcriptionally active genes in discrete regions in the nucleus by RNA polymerase II [72], the shepherd model proposed here emphasized that genes activated by common upstream regulatory TFs are spatially clustered in 3D space, enabling synchronized activation (or repression) of their activities. While transcription factors are generally considered to be part of transcription factories, it does not explain stage- and tissue-specific regulation of gene expression; the shepherd model addresses the question how transcription factories attain specificity in space and time. This “shepherd” working model provided the chance for efficient gene transcription activation or repression.

By combing chromatin interaction and TF binding-based gene regulatory network and genetic variation, we identified TaZNFs that encode C2H2-type zinc finger proteins, as flowering repressors in wheat (Fig. 5a, b). C2H2-type zinc finger proteins have been reported to be associated with flowering time control in higher plants, such as SUF4 [73] and LATE [74] as negative regulators in Arabidopsis and SIP1 [75] and RID1 [76] as positive regulators in rice. Although higher transcription level of TaZNF10 was detected in V28N6 than V0 plants, much-delayed flowering time of TaZNF10-5B-OE plant (Fig. 5c, d) suggests the negative role of TaZNF10-5B in wheat flowering time control. The increased expression does not necessarily mean positive effect. TaSOC1, a negative regulator in wheat flowering [64], also showed an increased expression level after vernalization. We speculated that the transcription of TaZNF10 in V28N6 plants may be co-regulated by multiple upstream regulators, and some of them are negative regulators of flowering time in wheat.

Conclusions

We performed a multi-omics analysis to dissect the gene transcription regulation mechanism of vernalization response in wheat. We discovered a “shepherd” model for gene transcription activation during response to vernalization, in which the targets of upstream regulator are gathered together via chromatin loops for efficient recognition and binding by upstream regulator. We provide insights into gene transcription regulation in wheat vernalization response. Further studies under more biotic/abiotic conditions in more species would help to answer whether this is general mechanism in responding to environmental stimulus.

Methods

Plant materials and treatment

The winter wheat cultivar AK58 was used in this study. After being sterilized with 75% ethanol for 45 s, the seeds were soaked in water containing 1.6% H2O2 for 24 h at room temperature and pre-germinated on wet filter paper. The germinated seeds were transferred to soil and grown in a growth chamber with 16 h light and 8 h dark cycle at 22 °C for normal growth and at 4 °C for vernalization treatment. For V0 treatment, plants were grown at 22 °C for 14 days. For V28 treatment, plants were grown at 22 °C for 9 days and were transferred to 4 °C for 28 days. For V28N6 treatment, plants were grown at 22 °C for 3 days, followed by cold exposure at 4 °C for 28 days and then normal growth at 22 °C for 6 days. Then, we harvested all the entire aerial parts of AK58 seedlings (including leaves, shoot, and shoot apical meristem) from different treatments but with synchronized developmental status at three-leaf stages.

RNA-seq library preparation

The fresh seedlings of AK58 were harvested and ground into powder with liquid nitrogen. Total RNA was extracted using TRIzol reagent (Invitrogen, 15596026) according to the manufacturer’s instructions. In brief, 1 mL TRIzol was added to 0.2 g powder, vortex for 1 min and incubated at room temperature for 10 min. After centrifuging at 13,000 rpm for 10 min at 4 °C, all the supernatant was collected and mixed with 200 μL chloroform by vigorous shaking. The mixture was incubated at room temperature for 20 min and was then centrifuged at 13,000 rpm for 15 min at 4 °C. The supernatant was collected carefully and was intensively mixed with 500 μL isopropanol to precipitate the total RNA. After washing twice with 75% alcohol prepared in DEPC water, RNA was collected, air-dried, and dissolved with 40 μL DEPC water. The MGIEasy RNA Library Prep Kit was used to construct RNA libraries, and a MGI DNBseq-T7 (paired-end 150 bp reads) platform was used for sequencing.

Sample preparation for ChIP-seq and ChIA-PET

Samples for ChIP-seq and ChIA-PET experiments were cross-linked immediately after harvested according to published protocol with slight modifications [41, 42]. Briefly, samples for ChIA-PET were cut into pieces and dual cross-linked with 1.5 mM ethylene glycol bis (succinimidyl succinate) (EGS, Thermo Scientific, 21565) for 20 min and 1% formaldehyde (Thermo Scientific, 28906) for another 20 min in MC buffer (10 mM sodium phosphate, PH7.0, 50 mM NaCl and 0.1 M sucrose) [77] under vacuum at 37 °C. The fixation was stopped by adding glycine solution to a final concentration of 0.125 M. The fixed tissues were washed with MC buffer three times and with sterilized water for another three times. After residual water was absorbed using filter paper, the fixed tissues were quick-frozen in liquid nitrogen and stored at –80 °C until use. For ChIP-seq, the fixation with EGS is omitted.

Immunoprecipitation and ChIP-seq library preparation

ChIP-seq was conducted according to the published protocols with minor modification [77, 78]. About 0.2 g crosslinked tissue was ground into power in liquid nitrogen, and 1.5 × volume (300 μL) Buffer S (50 mM HEPES–KOH (pH 7.5), 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, 1% SDS, 1 tablet protease inhibitor cocktail/50 mL) was added to lyse the tissue with rotation/incubation at 4 °C for 30 min. After adding 2 × volume (600 μL) Buffer F (50 mM HEPES–KOH (pH 7.5), 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, 1 tablet protease inhibitor cocktail/50 mL) and incubation at 4 °C for another 15 min, the chromatin was fragmented by sonication using a sonicator (Covaris S220, peak power: 140; duty factor: 6.0; cycles/burst: 200; duration: 500 s) for 10 min and was centrifuged at maximum speed for 10 min at 4 °C twice and then transfer the supernatant containing fragmented chromatin to a new tube and incubated with 30 μL Dynabeads™ protein A (Invitrogen, 10002D) for pre-clearing with rotation at 4 °C for 1 h. After separating from the beads with magnetic rack, 50 μL supernatant was transferred as the input sample. The remaining part of chromatin supernatant was incubated with 40 μL Dynabeads™ protein A and 5 μL antibodies against H3K4me3 (ABclonal, A2357), H3K9ac (Cell Signaling Technology, 9649S), and H3K27me3 (Cell Signaling Technology, 9733S) at 4 °C with rotation for immunoprecipitation overnight. Then, the beads were washed with 1 mL low-salt ChIP buffer (50 mM HEPES–KOH (PH7.5), 150 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, 0.1% SDS) three times, 1 mL high-salt ChIP buffer (50 mM HEPES–KOH (PH7.5), 350 mM NaCl, 1 mM EDTA, 1% Triton X-100, 0.1% sodium deoxycholate, 0.1% SDS) twice, 1 mL ChIP Wash/LiCl Buffer (10 mM Tris–HCl (PH8.0), 250 mM LiCl, 0.5% NP-40, 1 mM EDTA, 0.1% sodium deoxycholate) once, and 1 mL TE buffer (10 mM Tris–HCl (PH 8.0), 1 mM EDTA) twice in order. The protein-DNA complexes were eluted from beads by addition of 100 μL fresh prepared ChIP Elution buffer (8.5 mM Tris–HCl (PH 8.0), 0.85 mM EDTA, 1% SDS, 5 mM DTT) and was incubated at 75 °C for 10 min with 900 rpm. The supernatant was mixed with 15 μL 5 M NaCl and 12–13 μL proteinase K (20 mg/mL) and was incubated at 65 °C overnight for decrosslinking. The ChIP DNA was precipitated with solution containing 750 μL absolute ethanol, 30 μL 3 M sodium acetate, and 1 μL glycogen (50 mg/mL). After centrifugation at 13,000 rpm for 30 min at 4 °C, the air-dried ChIP DNA was eluted in 100 μL TE buffer. The QIAquick ® PCR Purification Kit (Qiagen, 28106) was used to purify the ChIP DNA. The ThruPLEX DNA-Seq Kit (Takara, R400675) and DNA Single Index Kit (Takara, R400695) were used to prepare ChIP-seq libraries following the manufacturer. Then, fragments ranging from 250 to 700 bp in the library were selected using SPRIselect beads (Beckman Coulter, B23318) and Illumina NovaSeq 6000 platform was applied for sequencing (paired-end 150 bp reads).

ATAC-seq library preparation

For ATAC-seq library preparation, fresh wheat tissues were chopped into homogenate immediately after harvesting in pre-cooled D-Sorbitol buffer (0.6 M D-sorbitol, 20 mM MES, 20 mM KCl, 10 mM MgCl2, 0.2% Triton X-100, 1 mM β-mercaptoethanol, 1 × PMSF, 0.1 mg/mL BSA, 1 tablet protease inhibitor cocktail/50 mL). The homogenate was filtered twice with 30-μm cell strainer. The generated crude nuclei were stained with 4′,6-diamidino-2-phenylindole (DAPI), and 50,000 nuclei were sorted in 400 μL pre-cooled D-Sorbitol buffer using a flow cytometer (BD FacsAria II SORP). Centrifugation at 1000 g for 10 min at 4 °C was performed, and the supernatant was removed carefully. The sorted nuclei were collected and was mixed with TD buffer (10 mM Tris–HCl, 5 mM MgCl2, 10% DMF) and 1.5 μL TTE Mix (Vazyme, TD501) immediately. After incubation at 37 °C for 30 min, the reaction was terminated, and DNA was purified using MinElute PCR Purification Kit (Qiagen, 28006). NEBNext® High-Fidelity 2X PCR Master Mix (NEB, M0541L) and TruePrep® Index Kit V2 for Illumina (Vazyme, TD202) were used to prepare the ATAC-seq library. The library fragments were selected with SPRIselect beads (Beckman Coulter, B23318) and sequenced with Illumina NovaSeq 6000 system (paired-end 150 bp).

Long-read ChIA-PET library preparation

ChIA-PET experiment and library preparation were performed according to the published method with modifications [41, 42, 50, 78], and the ChIP part was similar to that of ChIP-seq described above. In brief, about 4.5 g cross-linked tissues were ground into fine powder in liquid nitrogen, mixing thoroughly with 6.75 mL Buffer S. After incubation for 30 min at 4 °C with rotation, 27 mL Buffer F were added, incubating for another 15 min at 4 °C. The chromatin was sonicated in a BioRuptor (Diagenode) with high-level intensity (30 s ON and 50 s OFF for 20 cycles) into 1–3-kb fragments. After centrifugation at 12,000 rpm for 10 min at 4 °C, the supernatant was mixed with 100 μL Dynabeads™ protein A and rotated 3 h at 4 °C for pre-cleaning. To generate antibody-coated beads, another 800 μL Dynabeads™ protein A was mixed with 100 μL H3K4me3 antibody (ABclonal, A2357) and was incubated for 8 h at 4 °C with rotation. After being separated from pre-cleaning beads using a magnetic rack, the supernatant containing fragmented chromatin was mixed with antibody-coated beads and was incubated at 4 °C overnight with rotation for immunoprecipitation. Then, the beads were washed with 5 mL low-salt ChIP buffer three times, 5 mL high-salt ChIP buffer twice, 5 mL ChIP wash buffer, and 5 mL TE buffer in order. The end-repairing and A-tailing of ChIP DNA was conducted using T4 DNA polymerase (Promega, M4215) and Klenow fragment enzyme (NEB, M0212L), respectively, and the proximity ligation of ChIP DNA was performed using biotinylated bridge linker [50]. The protein-DNA complex was reverse cross-linked using 12 μL 5 M NaCl and 10 μL proteinase K; then, ChIP DNA was extracted using phenol–chloroform-isoamyl alcohol (pH 7.9) and precipitated by addition of 650 μL isopropanol, 60 μL 3 M sodium acetate, and 2 μL 15 mg/mL GlycoBlue (Invitrogen, AM9515). Tn5 transposase from TruePrep DNA Library Prep Kit V2 for Illumina (Vazyme, TD501) was used to construct the ChIA-PET library. Fragments from library were selected using AMPure XP beads (Beckman Coulter, A63880) and sequenced using Illumina NovaSeq 6000 system (paired-end 150 bp).

DAP-seq library preparation

DAP-seq experiment was conducted as previously described [59, 79] with minor modifications. In brief, VRN1 was fused with HA tags in pTnT vector (pTnT-C-3 × HA). The VRN1 protein was expressed using Quick Coupled Transcription/Translation System (Promega). The genomic DNA of Chinese spring was extracted and was sonicated to the average fragment size of 300 bp. Then, the Amp-DAP initial DNA library was prepared by adaptor ligation. The expressed VRN1 protein and Amp-DAP initial library were co-incubated at 4 °C and then was washed with TBS for 3 times followed by PBS-NP40 wash for 5 times. Subsequently, the protein-DNA complex was incubated at 95 °C; then, the DNA was purified using MinElute PCR Purification Kit (Qiagen, 28006). The VRN1-DAP library was sequenced using Illumina NovaSeq 6000 platform (paired-end 150 bp).

RNA-seq data analysis

Raw sequencing reads were trimmed by Trimmomatic v0.32 (https://github.com/usadellab/Trimmomatic) with parameters of “TruSeq3-PE.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36.” High-quality clean reads were used to quantified the transcript abundance via Kallisto v0.48.0 (https://github.com/pachterlab/kallisto) using IWGSC v1.0 assembly and v1.1 annotation as reference [1]. The output pseudoalignment bam files were converted to bigwig format file using deepTools v3.5.0 (https://github.com/deeptools/deepTools). Differentially expressed genes (DEGs) were identified using Sleuth v0.30.0 (http://pachterlab.github.io/sleuth); genes with |log2 (fold change)|> 1, P adjust value < 0.05, and at least expressed (TPM > 0.5) in one sample were regarded as DEGs.

ChIP-seq data analysis

Raw sequencing reads were trimmed as RNA-seq. High-quality clean reads were aligned to IWGSC v1.0 reference genome using Bowtie2 v2.4.4 (https://github.com/BenLangmead/bowtie2). Aligned reads with MAPQ lower than 5 were filtered using SAMTools (https://github.com/samtools/samtools), and duplicated alignments were removed using Picard (v2.23.9) (http://broadinstitute.github.io/picard/). The narrow peaks for H3K4me3 and H3K9ac were called by the ChIP-seq data analysis pipeline [80], which developed based on the ENCODE project recommend guidelines. The IDR (irreproducible discovery rate) value greater than 0.05, 0.05, and 0.01 for the original replicates, the pseudo-replicates split from original replicates, and pooled pseudo-replicates were used to obtain consistent peaks among replicates. Specifically, H3K27me3 peak calling was performed using model-based analysis of ChIP-seq data by MACS (v 2.2.7.1) as described previously [48, 81], applying a stringent P-value threshold of less than 1 × 10−5 and utilizing the broad-cutoff option set to 0.05. Peaks with at least one nucleotide of overlap between the two replicates were identified, retained, and merged for subsequent analysis. The peak annotation was performed by BEDTools v2.27 (https://github.com/arq5x/bedtools2), the promoter region was set from 3500 bp upstream to 1500 bp downstream of the TSS, the downstream 1500 bp of TSS to TES were defined as gene body, and the rest regions were called distal regions. To identify DEPs, we first merged the peaks from three different time points using “merge” command of BEDTools. Subsequently, we utilized “multicov” of BEDTools to calculate the count for each peak across different time points. A matrix with peaks as rows and count numbers for different time points as columns was created. This matrix was then used as the input file for Deseq2 to identify differentially modified peaks. To quantify the reads density on each peak, the fragment per million (FPM, similar to TPM in RNA-seq) value was calculated using number of read per fragment divided by sum of all reads per fragment in millions.

ATAC-seq data analysis

Raw sequencing reads were trimmed by Trimmomatic v0.32 with parameters of “NexteraPE-PE.fa:2:30:10:8:TRUE LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 -threads 32 MINLEN:36.” The high-quality clean reads were processed using WheatATAC pipeline, which consist of reads mapping, peak calling, IDR calculation, combination of consistent, peak annotation, and data quality evaluation [82].

ChIA-PET data analysis

Raw sequencing reads after merging biological replicates were trimmed by referencing the parameters and standards applied in ATAC-seq data analysis. High-quality clean reads were processed using CHIA-PET computational package for automatic linker filtering, reads mapping to IWGSC v1.0 reference genome based on BWA, mapped reads purifying, dividing the reads into different categories, interaction calling, and data evaluation (https://github.com/GuoliangLi-HZAU/ChIA-PET_Tool). Specially, the anchor files for H3K4me3 and H3K9ac were generated using the 3-kb upstream and downstream region of the middle site of the peak summits of V0 and V28N6. The interaction frequency was calculated by HiCExplorer v3.5.1 (https://github.com/deeptools/HiCExplorer).

Detection of the VRN1 target genes

Tobias (https://github.com/loosolab/TOBIAS) was applied to predict TF footprints by analyzing ATAC-seq data and the DAP-seq data filtered with OCRs. The BEDTools “closest” command was utilized to identify DAP-seq peaks that was located within the OCRs. VRN1-binding targets were predicted based on the combination of TF footprint and the genomic regions defined by OCRs located DAP-seq peaks.

Hybridization chain reaction (HCR)-based DNA FISH

HCR-based DNA FISH was conducted following the published method [83]. In brief, the specific probe for target chromatin interaction regions were designed and labeled with cy3, cy5, and fam dyes. Nuclei were isolated from crosslinked (1.5 mM EGS and 1% formaldehyde orderly) AK58 seedlings (aerial parts) and were placed on slide for air-dry. After denaturation in PBS buffer containing 0.2% Triton X-100, dehydration in ethanol gradient (70%, 90%, 100%), RNA digestion with RNase, and post-fixation with 4% formaldehyde, the hybridization buffer containing 0.2 pmol of each specific targeting probe was added to the nuclei for incubation at 45 °C overnight. Then, the slides were washed with 75%, 50%, and 25% probe wash buffer gradient in order to remove excess probe. After pre-amplification with amplification buffer, the nuclei were incubated with amplification buffer containing 6 pmol of each fluorescently labeled hairpin at room temperature (25 °C) overnight. After washing with 5 × SSCT and PBS buffer orderly to remove excess hairpins, nuclei were stained with DAPI and observed under a Nikon SIM Super-resolution microscope.

Generation of transgenic lines

The full-length cDNA of TaZNF10-5B (TraesCS5B02G406000) was isolated from AK58. For the TaZNF10-5B-OE (over expression), the full entire coding sequence without stop codons was recombined into the destination vector, modified pCAMBIA3300 in which the maize ubiquitin promoter is placed before the multiple cloning site (MCS), using the one step cloning kit (Vazyme, C112). To generate the mutant of TaZNF10, two gRNAs targeting the coding region were designed (Additional file 1: Fig. S14), and modified CRISPR-Cas9 vector pBUE414, in which TaU3 replaces the OsU3 promoter, was used to express Cas9 and the gRNAs. The vector was constructed using the Golden Gate Assembly (NEB, E1601). The plasmids were transformed into Agrobacterium tumefaciens strain EHA105 (WEIDI, AC1010). Wheat transformation was performed by infecting immature embryos of KN199 with EHA105 carrying the constructs. Primer targeting ubiquitin promoter and gene-specific primer were combined to identify positive TaZNF10-5B-OE lines, and gene-specific primers were designed to genotype mutation created by CRISPR-Cas9. All the primers used for construction and identification are listed in Additional file 2: Table S15.

Data availability

All data harvested from RNA-seq, ChIP-seq, ATAC-seq, and ChIA-PET experiments in this study are available at the Genome Sequence Archive at the Beijing Institute of Genomics (BIG) Data Center, Chinese Academy of Sciences, under accession number CRA015030 [84]. The VRN1 DAP-seq data generated in this study is available at the Genome Sequence Archive at the Beijing Institute of Genomics (BIG) Data Center, Chinese Academy of Sciences, under accession number of CRA020327 [85]. The VRT2 DAP-seq data was downloaded from NCBI Gene Expression Omnibus under the accession number of GSM6019668 [86]. The datasets used for the natural variation analysis were downloaded from China National GeneBank DataBase under the accession number of CNP0003712 [87]. The code used in this study is available in GitHub under GPL-3.0 license (https://github.com/hcph/Wheat_Vernalization_Regulation_Network) [88] and in Zenodo (https://zenodo.org/records/14060044) [89].

References

  1. International Wheat Genome Sequencing Consortium (IWGSC). Shifting the limits in wheat research and breeding using a fully annotated reference genome. Science. 2018;361(6403):eaar7191.

  2. Finnegan EJ. Vernalization. Curr Biol. 2012;22(12):R471–2.

    Article  CAS  PubMed  Google Scholar 

  3. Xu S, Chong K. Remembering winter through vernalisation. Nat Plants. 2018;4(12):997–1009.

    Article  PubMed  Google Scholar 

  4. Kim DH, Doyle MR, Sung S, Amasino RM. Vernalization: winter and the timing of flowering in plants. Annu Rev Cell Dev Biol. 2009;25:277–99.

    Article  CAS  PubMed  Google Scholar 

  5. Amasino RM. Vernalization and flowering time. Curr Opin Biotechnol. 2005;16(2):154–8.

    Article  CAS  PubMed  Google Scholar 

  6. Trevaskis B, Hemming MN, Dennis ES, Peacock WJ. The molecular basis of vernalization-induced flowering in cereals. Trends Plant Sci. 2007;12(8):352–7.

    Article  CAS  PubMed  Google Scholar 

  7. Penrose L. Yield of early dryland sowing of wheat with winter and spring habit in southern and central New South Wales. Aust J Exp Agr. 1993;33(5):601–8.

    Article  Google Scholar 

  8. Cann DJ, Schillinger WF, Hunt JR, Porker KD, Harris FAJ. Agroecological advantages of early-sown winter wheat in semi-arid environments: a comparative case study from southern Australia and Pacific Northwest United States. Front Plant Sci. 2020;11:568.

    Article  PubMed  PubMed Central  Google Scholar 

  9. Xiao J, Liu B, Yao Y, Guo Z, Jia H, Kong L, Zhang A, Ma W, Ni Z, Xu S, et al. Wheat genomic study for genetic improvement of traits in China. Sci China Life Sci. 2022;65(9):1718–75.

    Article  PubMed  Google Scholar 

  10. Trevaskis B. Wheat gene for all seasons. Proc Natl Acad Sci U S A. 2015;112(39):11991–2.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  11. Yan L, Loukoianov A, Tranquilli G, Helguera M, Fahima T, Dubcovsky J. Positional cloning of the wheat vernalization gene VRN1. Proc Natl Acad Sci U S A. 2003;100(10):6263–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  12. Yan L, Loukoianov A, Blechl A, Tranquilli G, Ramakrishna W, SanMiguel P, Bennetzen JL, Echenique V, Dubcovsky J. The wheat VRN2 gene is a flowering repressor down-regulated by vernalization. Science. 2004;303(5664):1640–4.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  13. Yan L, Fu D, Li C, Blechl A, Tranquilli G, Bonafede M, Sanchez A, Valarik M, Yasuda S, Dubcovsky J. The wheat and barley vernalization gene VRN3 is an orthologue of FT. Proc Natl Acad Sci U S A. 2006;103(51):19581–6.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  14. Kippes N, Debernardi JM, Vasquez-Gross HA, Akpinar BA, Budak H, Kato K, Chao S, Akhunov E, Dubcovsky J. Identification of the VERNALIZATION 4 gene reveals the origin of spring growth habit in ancient wheats from South Asia. Proc Natl Acad Sci U S A. 2015;112(39):E5401–10.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  15. Chen A, Dubcovsky J. Wheat TILLING mutants show that the vernalization gene VRN1 down-regulates the flowering repressor VRN2 in leaves but is not essential for flowering. PLoS Genet. 2012;8(12): e1003134.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  16. Deng W, Casao MC, Wang P, Sato K, Hayes PM, Finnegan EJ, Trevaskis B. Direct links between the vernalization response and other key traits of cereal crops. Nat Commun. 2015;6:5882.

    Article  PubMed  Google Scholar 

  17. Dhillon T, Pearce SP, Stockinger EJ, Distelfeld A, Li C, Knox AK, Vashegyi I, Vagujfalvi A, Galiba G, Dubcovsky J. Regulation of freezing tolerance and flowering in temperate cereals: the VRN-1 connection. Plant Physiol. 2010;153(4):1846–58.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  18. Simpson GG. The autonomous pathway: epigenetic and post-transcriptional gene regulation in the control of Arabidopsis flowering time. Curr Opin Plant Biol. 2004;7(5):570–4.

    Article  CAS  PubMed  Google Scholar 

  19. Sung S, Amasino RM. Vernalization and epigenetics: how plants remember winter. Curr Opin Plant Biol. 2004;7(1):4–10.

    Article  CAS  PubMed  Google Scholar 

  20. Song J, Angel A, Howard M, Dean C. Vernalization - a cold-induced epigenetic switch. J Cell Sci. 2012;125(16):3723–31.

    CAS  PubMed  Google Scholar 

  21. Zografos BR, Sung S. Vernalization-mediated chromatin changes. J Exp Bot. 2012;63(12):4343–8.

    Article  CAS  PubMed  Google Scholar 

  22. Michaels SD, Amasino RM. FLOWERING LOCUS C encodes a novel MADS domain protein that acts as a repressor of flowering. Plant Cell. 1999;11(5):949–56.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  23. Sheldon CC. The molecular basis of vernalization: the central role of FLOWERING LOCUS C (FLC). Proc Natl Acad Sci U S A. 2000;97(7):3753–8.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  24. Whittaker C, Dean C. The FLC locus: a platform for discoveries in epigenetics and adaptation. Annu Rev Cell Dev Bi. 2017;33:555–75.

    Article  CAS  Google Scholar 

  25. Hepworth J, Dean C. Flowering Locus C’s lessons: conserved chromatin switches underpinning developmental timing and adaptation. Plant Physiol. 2015;168(4):1237–45.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Bastow R, Mylne JS, Lister C, Lippman Z, Martienssen RA, Dean C. Vernalization requires epigenetic silencing of FLC by histone methylation. Nature. 2004;427(6970):164–7.

    Article  CAS  PubMed  Google Scholar 

  27. Oliver SN, Finnegan EJ, Dennis ES, Peacock WJ, Trevaskis B. Vernalization-induced flowering in cereals is associated with changes in histone methylation at the VERNALIZATION1 gene. Proc Natl Acad Sci U S A. 2009;106(20):8386–91.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  28. Diallo AO, Ali-Benali MA, Badawi M, Houde M, Sarhan F. Expression of vernalization responsive genes in wheat is associated with histone H3 trimethylation. Mol Genet Genomics. 2012;287(7):575–90.

    Article  CAS  PubMed  Google Scholar 

  29. Oliver SN, Deng W, Casao MC, Trevaskis B. Low temperatures induce rapid changes in chromatin state and transcript levels of the cereal VERNALIZATION1 gene. J Exp Bot. 2013;64(8):2413–22.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  30. Huan Q, Mao Z, Chong K, Zhang J. Global analysis of H3K4me3/H3K27me3 in Brachypodium distachyon reveals VRN3 as critical epigenetic regulation point in vernalization and provides insights into epigenetic memory. New Phytol. 2018;219(4):1373–87.

    Article  CAS  PubMed  Google Scholar 

  31. Xiao J, Xu S, Li C, Xu Y, Xing L, Niu Y, Huan Q, Tang Y, Zhao C, Wagner D, et al. O-GlcNAc-mediated interaction between VER2 and TaGRP2 elicits TaVRN1 mRNA accumulation during vernalization in winter wheat. Nat Commun. 2014;5:4572.

    Article  CAS  PubMed  Google Scholar 

  32. Niu D, Gao Z, Cui B, Zhang Y, He Y. A molecular mechanism for embryonic resetting of winter memory and restoration of winter annual growth habit in wheat. Nat Plants. 2024;10(1):37–52.

    Article  CAS  PubMed  Google Scholar 

  33. Xu S, Xiao J, Yin F, Guo X, Xing L, Xu Y, Chong K. The protein modifications of O-GlcNAcylation and phosphorylation mediate vernalization response for flowering in winter wheat. Plant Physiol. 2019;180(3):1436–49.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  34. Xie L, Zhang Y, Wang K, Luo X, Xu D, Tian X, Li L, Ye X, Xia X, Li W, et al. TaVrt2, an SVP-like gene, cooperates with TaVrn1 to regulate vernalization-induced flowering in wheat. New Phytol. 2021;231(2):834–48.

    Article  CAS  PubMed  Google Scholar 

  35. Xu S, Dong Q, Deng M, Lin D, Xiao J, Cheng P, Xing L, Niu Y, Gao C, Zhang W, et al. The vernalization-induced long non-coding RNA VAS functions with the transcription factor TaRF2b to promote TaVRN1 expression for flowering in hexaploid wheat. Mol Plant. 2021;14(9):1525–38.

    Article  CAS  PubMed  Google Scholar 

  36. Lloyd JPB, Lister R. Epigenome plasticity in plants. Nat Rev Genet. 2022;23(1):55–68.

    Article  CAS  PubMed  Google Scholar 

  37. He G, Elling AA, Deng XW. The epigenome and plant development. Annu Rev Plant Biol. 2011;62(1):411–35.

    Article  CAS  PubMed  Google Scholar 

  38. Feil R, Fraga MF. Epigenetics and the environment: emerging patterns and implications. Nat Rev Genet. 2012;13(2):97–109.

    Article  CAS  PubMed  Google Scholar 

  39. Huan Q, Mao Z, Zhang J, Xu Y, Chong K. Transcriptome-wide analysis of vernalization reveals conserved and species-specific mechanisms in Brachypodium. J Integr Plant Biol. 2013;55(8):696–709.

    Article  CAS  PubMed  Google Scholar 

  40. Tang Z, Luo OJ, Li X, Zheng M, Zhu JJ, Szalaj P, Trzaskoma P, Magalska A, Wlodarczyk J, Ruszczycki B, et al. CTCF-mediated human 3D genome architecture reveals chromatin topology for transcription. Cell. 2015;163(7):1611–27.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  41. Zhao L, Wang S, Cao Z, Ouyang W, Zhang Q, Xie L, Zheng R, Guo M, Ma M, Hu Z, et al. Chromatin loops associated with active genes and heterochromatin shape rice genome architecture for transcriptional regulation. Nat Commun. 2019;10(1):3640.

    Article  PubMed  PubMed Central  Google Scholar 

  42. Deng L, Gao B, Zhao L, Zhang Y, Zhang Q, Guo M, Yang Y, Wang S, Xie L, Lou H, et al. Diurnal RNAPII-tethered chromatin interactions are associated with rhythmic gene expression in rice. Genome Biol. 2022;23(1):7.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  43. Concia L, Veluchamy A, Ramirez-Prado JS, Martin-Ramirez A, Huang Y, Perez M, Domenichini S, Rodriguez Granados NY, Kim S, Blein T, et al. Wheat chromatin architecture is organized in genome territories and transcription factories. Genome Biol. 2020;21(1):104.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  44. Jia J, Xie Y, Cheng J, Kong C, Wang M, Gao L, Zhao F, Guo J, Wang K, Li G, et al. Homology-mediated inter-chromosomal interactions in hexaploid wheat lead to specific subgenome territories following polyploidization and introgression. Genome Biol. 2021;22(1):26.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  45. Yuan J, Sun H, Wang Y, Li L, Chen S, Jiao W, Jia G, Wang L, Mao J, Ni Z, et al. Open chromatin interaction maps reveal functional regulatory elements and chromatin architecture variations during wheat evolution. Genome Biol. 2022;23(1):34.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  46. Jia J, Zhao G, Li D, Wang K, Kong C, Deng P, Yan X, Zhang X, Lu Z, Xu S, et al. Genome resources for the elite bread wheat cultivar Aikang 58 and mining of elite homeologous haplotypes for accelerating wheat improvement. Mol Plant. 2023;16(12):1893–910.

    Article  CAS  PubMed  Google Scholar 

  47. Li Z, Wang M, Lin K, Xie Y, Guo J, Ye L, Zhuang Y, Teng W, Ran X, Tong Y, et al. The bread wheat epigenomic map reveals distinct chromatin architectural and evolutionary features of functional genetic elements. Genome Biol. 2019;20(1):139.

    Article  PubMed  PubMed Central  Google Scholar 

  48. Wang M, Li Z, Zhang Y, Zhang Y, Xie Y, Ye L, Zhuang Y, Lin K, Zhao F, Guo J, et al. An atlas of wheat epigenetic regulatory elements reveals subgenome divergence in the regulation of development and stress responses. Plant Cell. 2021;33(4):865–81.

    Article  PubMed  PubMed Central  Google Scholar 

  49. Bajic M, Maher KA, Deal RB. Identification of open chromatin regions in plant genomes using ATAC-seq. Methods Mol Biol. 2018;1675:183–201.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  50. Li X, Luo OJ, Wang P, Zheng M, Wang D, Piecuch E, Zhu JJ, Tian SZ, Tang Z, Li G, Ruan Y. Long-read ChIA-PET for base-pair-resolution mapping of haplotype-specific chromatin interactions. Nat Protoc. 2017;12(5):899–915.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  51. Zhang L, He C, Lai Y, Wang Y, Kang L, Liu A, Lan C, Su H, Gao Y, Li Z, et al. Asymmetric gene expression and cell-type-specific regulatory networks in the root of bread wheat revealed by single-cell multiomics analysis. Genome Biol. 2023;24(1):65.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  52. Greenup AG, Sasani S, Oliver SN, Talbot MJ, Dennis ES, Hemming MN, Trevaskis B. ODDSOC2 is a MADS box floral repressor that is down-regulated by vernalization in temperate cereals. Plant Physiol. 2010;153(3):1062–73.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  53. Chen A, Li C, Hu W, Lau MY, Lin H, Rockwell NC, Martin SS, Jernstedt JA, Lagarias JC, Dubcovsky J. Phytochrome C plays a major role in the acceleration of wheat flowering under long-day photoperiod. Proc Natl Acad Sci U S A. 2014;111(28):10037–44.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  54. Fan M, Miao F, Jia H, Li G, Powers C, Nagarajan R, Alderman PD, Carver BF, Ma Z, Yan L. O-linked N-acetylglucosamine transferase is involved in fine regulation of flowering time in winter wheat. Nat Commun. 2021;12(1):2303.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  55. Zikhali M, Wingen LU, Leverington-Waite M, Specel S, Griffiths S. The identification of new candidate genes Triticum aestivum FLOWERING LOCUS T3–B1 (TaFT3-B1) and TARGET OF EAT1 (TaTOE1-B1) controlling the short-day photoperiod response in bread wheat. Plant Cell Environ. 2017;40(11):2678–90.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  56. Mizuno N, Kinoshita M, Kinoshita S, Nishida H, Fujita M, Kato K, Murai K, Nasuda S. Loss-of-function mutations in three homoeologous PHYTOCLOCK 1 genes in common wheat are associated with the extra-early flowering phenotype. PLoS ONE. 2016;11(10): e0165618.

    Article  PubMed  PubMed Central  Google Scholar 

  57. Peng Y, Xiong D, Zhao L, Ouyang W, Wang S, Sun J, Zhang Q, Guan P, Xie L, Li W, et al. Chromatin interaction maps reveal genetic regulation for quantitative traits in maize. Nat Commun. 2019;10(1):2632.

    Article  PubMed  PubMed Central  Google Scholar 

  58. Krijger PHL, de Laat W. Regulation of disease-associated gene expression in the 3D genome. Nat Rev Mol Cell Biol. 2016;17(12):771–82.

    Article  CAS  PubMed  Google Scholar 

  59. Bartlett A, O’Malley RC, Huang SC, Galli M, Nery JR, Gallavotti A, Ecker JR. Mapping genome-wide transcription-factor binding sites using DAP-seq. Nat Protoc. 2017;12(8):1659–72.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  60. Zhang Y, Li Z, Liu J, Zhang Y, Ye L, Peng Y, Wang H, Diao H, Ma Y, Wang M, et al. Transposable elements orchestrate subgenome-convergent and -divergent transcription in common wheat. Nat Commun. 2022;13(1):6940.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  61. Schilling S, Kennedy A, Pan S, Jermiin LS, Melzer R. Genome-wide analysis of MIKC-type MADS-box genes in wheat: pervasive duplications, functional conservation and putative neofunctionalization. New Phytol. 2020;225(1):511–29.

    Article  CAS  PubMed  Google Scholar 

  62. Gao J, Hu X, Gao C, Chen G, Feng H, Jia Z, Zhao P, Yu H, Li H, Geng Z, et al. Deciphering genetic basis of developmental and agronomic traits by integrating high-throughput optical phenotyping and genome-wide association studies in wheat. Plant Biotechnol J. 2023;21(10):1966–77.

    Article  PubMed  PubMed Central  Google Scholar 

  63. Zhang J, Li XM, Lin HX, Chong K. Crop improvement through temperature resilience. Annu Rev Plant Biol. 2019;70:753–80.

    Article  CAS  PubMed  Google Scholar 

  64. Luo X, Liu B, Xie L, Wang K, Xu D, Tian X, Xie L, Li L, Ye X, He Z, et al. The TaSOC1-TaVRN1 module integrates photoperiod and vernalization signals to regulate wheat flowering. Plant Biotechnol J. 2023;22(3):635–49.

    Article  PubMed  PubMed Central  Google Scholar 

  65. Thomashow MF, Gilmour SJ, Stockinger EJ, Jaglo-Ottosen KR, Zarka DG. Role of the Arabidopsis CBF transcriptional activators in cold acclimation. Physiol Plantarum. 2001;112(2):171–5.

    Article  CAS  Google Scholar 

  66. Amasino R. Vernalization, competence, and the epigenetic memory of winter. Plant Cell. 2004;16(10):2553–9.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  67. Xiao J, Lee US, Wagner D. Tug of war: adding and removing histone lysine methylation in Arabidopsis. Curr Opin Plant Biol. 2016;34:41–53.

    Article  CAS  PubMed  Google Scholar 

  68. Liang Z, Li G, Wang Z, Djekidel MN, Li Y, Qian M-P, Zhang MQ, Chen Y. BL-Hi-C is an efficient and sensitive approach for capturing structural and regulatory chromatin interactions. Nat Commun. 2017;8(1):1622.

    Article  PubMed  PubMed Central  Google Scholar 

  69. Lieberman-Aiden E, van Berkum NL, Williams L, Imakaev M, Ragoczy T, Telling A, Amit I, Lajoie BR, Sabo PJ, Dorschner MO, et al. Comprehensive mapping of long-range interactions reveals folding principles of the human genome. Science. 2009;326(5950):289–93.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  70. Li G, Ruan X, Auerbach Raymond K, Sandhu Kuljeet S, Zheng M, Wang P, Poh Huay M, Goh Y, Lim J, Zhang J, et al. Extensive promoter-centered chromatin interactions provide a topological basis for transcription regulation. Cell. 2012;148(1–2):84–98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  71. Kaufmann K, Pajoro A, Angenent GC. Regulation of transcription in plants: mechanisms controlling developmental switches. Nat Rev Genet. 2010;11(12):830–42.

    Article  CAS  PubMed  Google Scholar 

  72. Shaban HA, Barth R, Bystricky K. Navigating the crowd: visualizing coordination between genome dynamics, structure, and transcription. Genome Biol. 2020;21(1):278.

    Article  PubMed  PubMed Central  Google Scholar 

  73. Kim S, Choi K, Park C, Hwang H-J, Lee I. SUPPRESSOR OF FRIGIDA4, encoding a C2H2-type zinc finger protein, represses flowering by transcriptional activation of Arabidopsis FLOWERING LOCUS C. Plant Cell. 2006;18(11):2985–98.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  74. Weingartner M, Subert C, Sauer N. LATE, a C2H2 zinc-finger protein that acts as floral repressor. Plant J. 2011;68(4):681–92.

    Article  CAS  PubMed  Google Scholar 

  75. Jiang P, Wang S, Zheng H, Li H, Zhang F, Su Y, Xu Z, Lin H, Qian Q, Ding Y. SIP1 participates in regulation of flowering time in rice by recruiting OsTrx1 to Ehd1. New Phytol. 2018;219(1):422–35.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  76. Wu C, You C, Li C, Long T, Chen G, Byrne ME, Zhang Q. RID1, encoding a Cys2/His2-type zinc finger transcription factor, acts as a master switch from vegetative to floral development in rice. Proc Natl Acad Sci U S A. 2008;105(35):12915–20.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  77. Kaufmann K, Muino JM, Osteras M, Farinelli L, Krajewski P, Angenent GC. Chromatin immunoprecipitation (ChIP) of plant transcription factors followed by sequencing (ChIP-SEQ) or hybridization to whole genome arrays (ChIP-CHIP). Nat Protoc. 2010;5(3):457–72.

    Article  CAS  PubMed  Google Scholar 

  78. Zhao L, Xie L, Zhang Q, Ouyang W, Deng L, Guan P, Ma M, Li Y, Zhang Y, Xiao Q, et al. Integrative analysis of reference epigenomes in 20 rice varieties. Nat Commun. 2020;11(1):2658.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  79. Lai X, Vega-Léon R, Hugouvieux V, Blanc-Mathieu R, van der Wal F, Lucas J, Silva CS, Jourdain A, Muino JM, Nanao MH, et al. The intervening domain is required for DNA-binding and functional identity of plant MADS transcription factors. Nat Commun. 2021;12(1):4760.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  80. Chen D, Yan W, Fu LY, Kaufmann K. Architecture of gene regulatory networks controlling flower development in Arabidopsis thaliana. Nat Commun. 2018;9(1):4534.

    Article  PubMed  PubMed Central  Google Scholar 

  81. Liu Y, Yuan J, Jia G, Ye W, Jeffrey Chen Z, Song Q. Histone H3K27 dimethylation landscapes contribute to genome stability and genetic recombination during wheat polyploidization. Plant J. 2020;105(3):678–90.

    Article  PubMed  Google Scholar 

  82. Wang X, Chen C, He C, Chen D, Yan W. Mapping open chromatin by ATAC-seq in bread wheat. Front Plant Sci. 2022;13:1074873.

    Article  PubMed  PubMed Central  Google Scholar 

  83. Ouyang W, Zhang X, Guo M, Wang J, Wang X, Gao R, Ma M, Xiang X, Luan S, Xing F, et al. Haplotype mapping of H3K27me3-associated chromatin interactions defines topological regulation of gene silencing in rice. Cell Rep. 2023;42(4): 112350.

    Article  CAS  PubMed  Google Scholar 

  84. Liu Y, Xu X, He C, Jin L, Zhou Z, Gao J, Guo M, Wang X, Chen C, Ayaad MH, Li X and Yan W. Chromatin loops gather targets of upstream regulators together for efficient gene transcription regulation during vernalization in wheat. Datasets. Genome Sequence Archive. https://ngdc.cncb.ac.cn/gsa/browse/CRA015030 (2024).

  85. Liu Y, Xu X, He C, Jin L, Zhou Z, Gao J, Guo M, Wang X, Chen C, Ayaad MH, Li X and Yan W. Chromatin loops gather targets of upstream regulators together for efficient gene transcription regulation during vernalization in wheat. Datasets. Genome Sequence Archive. https://ngdc.cncb.ac.cn/gsa/browse/CRA020327 (2024).

  86. Zhang Y, Li Z, Liu J, Zhang Y, Ye L, Peng Y, Wang H, Diao H, Ma Y, Wang M, Xie Y, Tang T, Zhuang Y, Teng W, Tong Y, Zhang W, Lang Z, Xue Y and Zhang Y. Transposable elements orchestrate subgenome-convergent and -divergent transcription in common wheat. Datasets. NCBI Gene Expression Omnibus. https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE192815 (2022).

  87. Gao J, Hu X, Gao C, Chen G, Feng H, Jia Z, Zhao P, Yu H, Li H, Geng Z, Fu J, Zhang J, Cheng Y, Yang B, Pang Z, Xiang D, Jia J, Su H, Mao H, Lan C, Chen W, Yan W, Gao L, Yang W and Li Q. Deciphering genetic basis of developmental and agronomic traits by integrating high-throughput optical phenotyping and genome-wide association studies in wheat. Datasets. China National Genebank. https://db.cngb.org/search/project/CNP0003712/(2023).

  88. Liu Y, Xu X, He C, Jin L, Zhou Z, Gao J, Guo M, Wang X, Chen C, Ayaad MH, Li X and Yan W. Chromatin loops gather targets of upstream regulators together for efficient gene transcription regulation during vernalization in wheat. Github. https://github.com/hcph/Wheat_Vernalization_Regulation_Network (2024).

  89. Liu Y, Xu X, He C, Jin L, Zhou Z, Gao J, Guo M, Wang X, Chen C, Ayaad MH, Li X and Yan W. Chromatin loops gather targets of upstream regulators together for efficient gene transcription regulation during vernalization in wheat. Zenodo. https://zenodo.org/records/14060044 (2024).

Download references

Acknowledgements

We thank Prof. Mingming Xin (China Agricultural University) for kindly providing the modified pBUE414 vector for generating CRISPR-Cas9 mutant, Prof. Kerstin Kaufmann (Humboldt-Universität zu Berlin) for valuable suggestion in improving the manuscript, and Miss Pei Chen from the wheat transformation platform for help with the generating the transgenic plants. We sincerely thank the computing platform of the National Key Laboratory of Crop Genetic Improvement in Huazhong Agricultural University for providing the computational resources.

Peer review information

Wenjing She was the primary editor of this article and managed its editorial process and peer review in collaboration with the rest of the editorial team.

Review history

The review history is available as Additional file 3.

Funding

This work was supported by grants from Biological Breeding-Major Projects (2023ZD04076), the National Key Research and Development Program of Hubei Province (2022BBA154), the Foundation of Hubei Hongshan Laboratory (2021hszd010) and the National Key Laboratory of Crop Genetic Improvement Self-Research Program (ZW19A0201), and the China Postdoctoral Science Foundation (2019M662665).

Author information

Author notes

  1. Yanyan Liu, Xintong Xu and Chao He contributed equally to this work.

Authors and Affiliations

  1. National Key Laboratory of Crop Genetic Improvement, Hubei Hongshan Laboratory, Huazhong Agricultural University, Wuhan, 430070, China

    Yanyan Liu, Xintong Xu, Chao He, Liujie Jin, Ziru Zhou, Jie Gao, Minrong Guo, Xin Wang, Chuanye Chen, Xingwang Li & Wenhao Yan

  2. Plant Research Department, Nuclear Research Center, Egyptian Atomic Energy Authority, Cairo, 13759, Egypt

    Mohammed H. Ayaad

Authors
  1. Yanyan Liu
    View author publications

    You can also search for this author inPubMed Google Scholar

  2. Xintong Xu
    View author publications

    You can also search for this author inPubMed Google Scholar

  3. Chao He
    View author publications

    You can also search for this author inPubMed Google Scholar

  4. Liujie Jin
    View author publications

    You can also search for this author inPubMed Google Scholar

  5. Ziru Zhou
    View author publications

    You can also search for this author inPubMed Google Scholar

  6. Jie Gao
    View author publications

    You can also search for this author inPubMed Google Scholar

  7. Minrong Guo
    View author publications

    You can also search for this author inPubMed Google Scholar

  8. Xin Wang
    View author publications

    You can also search for this author inPubMed Google Scholar

  9. Chuanye Chen
    View author publications

    You can also search for this author inPubMed Google Scholar

  10. Mohammed H. Ayaad
    View author publications

    You can also search for this author inPubMed Google Scholar

  11. Xingwang Li
    View author publications

    You can also search for this author inPubMed Google Scholar

  12. Wenhao Yan
    View author publications

    You can also search for this author inPubMed Google Scholar

Contributions

W.Y. conceived the study and designed the experiment; Y.L. and X.X. performed the experiments and generated the data with assistance from L.J., Z.Z., M.G., X.W., C.C., and X.L.; C.H. and Y.L. conducted the data analysis with support from J.G.; Y.L., X.X., and C.H. drafted the manuscript with input from M.A.; W.Y. organized the data and finalized the manuscript.

Corresponding author

Correspondence to Wenhao Yan.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare that they have no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Supplementary Information

13059_2024_3437_MOESM1_ESM.docx

Additional file 1: Fig. S1-14. Fig. S1 Schematic outline of multi-omics analysis of wheat vernalization. Samples were harvested at V0 (no vernalization), V28 (vernalization for 28 days) and V28N6 (vernalization for 28 days and recovered for 6 days). Fig. S2 Transcription analysis of AK58 vernalization process. a-c, Volcano plots showing differentially expressed genes in V0 vs V28 (a), V28 vs V28N6 (b) and V0 vs V28N6 (c) analysis. d, Expression trends of DEGs in the seven different clusters. e, Statistics of different homoeologs asymmetric expression pattern in V0, V28 and V28N6 plants. Triads were used to analysis when one or more of the subgenome homoeologs belongs to DEGs. f, Expression patterns of vernalization/flowering-related genes at V0, V28 and V28N6. Different colors indicating different homoeologs in A, B, D subgenomes. Fig. S3 Circular diagram showing the genome-wide change of histone modification and chromatin accessibility in V28 and V28N6 plants. Fig. S4 The relationship of histone modification and chromatin accessibility with gene expression. a, Heatmap showing correlation coefficient of RNA-seq dataset with ATAC-seq and ChIP-seq datasets of H3K4me3, H3K9ac and H3K27me3 harvested in this study. b, Box plots showing the expression level of genes with or without (none) corresponding histone modification or chromatin accessibility. ***P < 0.001. Fig. S5 The relationship of DEGs fold change with H3K4me3 (a), H3K9ac (b) and H3K27me3 (c) modification fold change in V0 vs V28 (top) and V28 vs V28N6 (bottom) comparations. Fig. S6 Heatmap showing seven clusters of differentially enriched peaks (DEPs) of H3K4me3, H3K9ac and H3K27me3. The red or purple color represents the peak intensity of the corresponding modification. Fig. S7 Statistics of anchor genes in different subgenomes that with (up, down, appear, disappear) or without (stable) altered interaction degree during wheat vernalization process. Fig. S8 Validation of changed chromatin interactions of VRN1-A (green), VRN1-B (red) and VRN1-D (carmine) locus during wheat vernalization process. Bar = 5 μm. Fig. S9 Altered chromatin interactions associated with H3K9ac during wheat vernalization. a-b, Boxplots showing interaction frequency within different subgenomes associated with H3K9ac modification in V0 (a) and V28N6 (b) plants. *P < 0.05, **P < 0.01, ***P < 0.001. c, Pie diagrams showing statistics of BP genes, anchor genes and rest genes in V0 and V28N6 plants. d, Expression level of genes with different types of chromatin interactions. e, Boxplots showing the expression level of BP genes (basal), anchor genes with different degrees, and rest genes without H3K9ac modification. f, Column diagram showing the statistics of PPI and PDI loops in the five types of dynamic chromatin interaction (stable, up, down, appear, disappear) before and after vernalization. g, h, Heatmaps showing the elevated expression level of anchor genes with increased interaction degree of PPI (g) and PDI (h) loops associated with H3K9ac. i, Venn diagram showing the overlaps of anchor genes identified using H3K4me3 and H3K9ac antibodies respectively in V0 (top) and V28N6 (bottom) plants. Fig. S10 Browser screenshots showing no H3K4me3-associated chromatin interaction surrounding the genomic regions of the three homoeologs of VRN2 (top) and VRN3 (bottom) in V0 and V28N6 plants. Fig. S11 Targets of VRT2 are gathering by H3K4me3-associated chromatin loops during vernalization. a, Read density of VRT2 DAP-seq data surrounding the transcription start site (TSS). b, CArG box was identified from publicly VRT2 DAP-seq data. c, Column diagram showing the number of putative VRT2 target genes connected by loops during vernalization process both in transcriptional upregulated (left) and downregulated (right) gene set. d, Number of loop-connected VRT2 target genes with different altered types (stable, up, down, appear, disappear) of interaction degree both in transcriptional upregulated (left) and downregulated (right) gene sets. Fig. S12 Venn diagram showing overlaps of links (a) and nodes (b) in V0 and V28N6 regulatory networks presented in Fig. 4f, g. Fig. S13 The dynamics of gene expression (left) and H3K4me3 histone modification (right) level of identified ZNF genes during vernalization. Fig. S14 Genotype of TaZNF10-5B-KO plants. a, Diagram showing the construct used for CRISPR-Cas9 gene editing of TaZNF10-5B. b, Scheme showing the structure of TaZNF10-5B gene and the target sites of designed sgRNAs. Black and white rectangles represent the exon and untranslated region of TaZNF10-5B respectively. Green or highlight letters represent the PAM site. The number -15 indicates the deleted bases in TaZNF10-5B-KO lines. g1, sgRNA1; g2, sgRNA2.

13059_2024_3437_MOESM2_ESM.xlsx

Additional file 2: Table S1-S15. Table S1 Statistics of datasets generated in this study. Table S2 Summary of RNA-seq libraries harvested in this study. Table S3 Differentially expressed genes during AK58 vernalization. Table S4 The reshaping of homoeolog asymmetric expression pattern during wheat vernalization. Table S5 Summary of ChIP-seq and ATAC-seq libraries harvsted in this study. Table S6 Differentially enriched peaks of H3K4me3 during AK58 vernalization. Table S7 Differentially enriched peaks of H3K9ac during AK58 vernalization. Table S8 Differentially enriched peaks of H3K27me3 during AK58 vernalization. Table S9 Summary of ChIA-PET libraries harvested in this study. Table S10 Histone modification H3K4me3-associated interactions in V0 plants. Table S11 Histone modification H3K4me3-associated interactions in V28N6 plants. Table S12 Histone modification H3K9ac-associated interactions in V0 plants. Table S13 Histone modification H3K9ac-associated interactions in V28N6 plants. Table S14 Gene list of transcription factor identified in TF footprint analysis. Table S15 Primers used in this study.

Additional file 3. Review history.

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, 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 licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence 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 licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Liu, Y., Xu, X., He, C. et al. Chromatin loops gather targets of upstream regulators together for efficient gene transcription regulation during vernalization in wheat. Genome Biol 25, 306 (2024). https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-024-03437-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doiorg.publicaciones.saludcastillayleon.es/10.1186/s13059-024-03437-x

Keywords

Genome Biology

ISSN: 1474-760X

Contact us