1. Introduction
Diurnal adaptation in plants involves a coordinated set of responses to optimize growth and resource allocation across the day-night cycle, as well as the vegetative-to-reproductive transition (Coen & Prusinkiewicz, Reference Coen and Prusinkiewicz2024; Nozue & Maloof, Reference Nozue and Maloof2006; Venkat & Muneer, Reference Venkat and Muneer2022; Wang et al., Reference Wang, Han and Jeffrey Chen2024). At the physiological level, plants adjust key processes such as photosynthesis, respiration and stomatal movement, water and ion transport (Matthews et al., Reference Matthews, Vialet-Chabrand and Lawson2018; Venkat & Muneer, Reference Venkat and Muneer2022). At the molecular level, diurnal variations are detected in the proteome, the transcriptome and the epigenome (Nozue & Maloof, Reference Nozue and Maloof2006; Wang et al., Reference Wang, Han and Jeffrey Chen2024 and references below). While some of these rhythms are regulated by circadian rhythms and remain stable despite temporary disturbances, others respond directly to environmental cues that vary throughout the day-night (diurnal) cycle (Nozue & Maloof, Reference Nozue and Maloof2006).
In Arabidopsis, diurnal variations of the proteome influence a variety of processes including light perception, photosynthesis, metabolism and ion transport (Uhrig et al., Reference Uhrig, Echevarria-Zomeno, Schlapfer, Grossmann, Roschitzki, Koerber and Gruissem2021). Diurnal variations extend to protein isoforms, such as acetylated and phosphorylated forms, which are crucial for the proteins’ functional activity (Uhrig et al., Reference Uhrig, Schlapfer, Roschitzki, Hirsch-Hoffmann and Gruissem2019; Uhrig et al., Reference Uhrig, Echevarria-Zomeno, Schlapfer, Grossmann, Roschitzki, Koerber and Gruissem2021). Protein translation, which, along with protein degradation, plays a critical role in controlling the abundance of proteins available to the cell, also shows diurnal variation (Duncan & Millar, Reference Duncan and Millar2022; Li et al., Reference Li, Duncan, Ganguly, Lee, Crisp, Wijerathna-Yapa and Millar2022; Seaton et al., Reference Seaton, Graf, Baerenfaller, Stitt, Millar and Gruissem2018; Uhrig et al., Reference Uhrig, Echevarria-Zomeno, Schlapfer, Grossmann, Roschitzki, Koerber and Gruissem2021). Protein abundance also correlates with transcriptional fluctuations, a process called translational coincidence (Seaton et al., Reference Seaton, Graf, Baerenfaller, Stitt, Millar and Gruissem2018). Moreover, diurnal transcriptome dynamics also include the occurrence of different transcript isoforms because of diurnal variation in splicing activity controlled by the circadian clock (Yang et al., Reference Yang, Li, Sancar and Oztas2020). Further, the pattern of diurnal transcriptional changes varies with age and development, suggesting a cross-talk between environmental signals and intrinsic cues on the physiological state of the plant (Jung et al., Reference Jung, Kim, Lee, Kang, Kim, Kim and Nam2024; Redmond et al., Reference Redmond, Ronald, Davis and Ezer2024).
In mice, transcriptional oscillations can be largely explained by rhythmic RNA Polymerase II (RNA Pol II) recruitment at promoters rather than rhythmic transition from paused to productive transcript elongation (Le Martelot et al., Reference Le Martelot, Canella, Symul, Migliavacca, Gilardi, Liechti and Cycli2012). In plants, RNA Pol II also exhibits genome-wide rhythmic occupancy, also shown in Rice to correlate with changes in the spatial organization of the genome (Deng et al., Reference Deng, Gao, Zhao, Zhang, Zhang, Guo and Li2022). Core clock genes show high spatial proximity in the nuclear space (chromatin connectivity) in the morning and weaker connectivity in the evening, indicating a shift from coordinated transcription in the morning to discretely regulated transcriptional regions in the evening (Deng et al., Reference Deng, Gao, Zhao, Zhang, Zhang, Guo and Li2022). Likewise, chromatin accessibility differs in Arabidopsis in response to photoperiod (Tian et al., Reference Tian, Li, Wang, Xu, Zhang, Zeb and Liu2021). Such diurnal changes in chromatin organization are also reflected at the epigenetic level as shown by studies conducted on Arabidopsis thaliana and Brassica napus (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016; Deng et al., Reference Deng, Gao, Zhao, Zhang, Zhang, Guo and Li2022; Fung-Uceda et al., Reference Fung-Uceda, Gomez, Rodriguez-Casillas, Gonzalez-Gil and Gutierrez2024; Maric & Mas, Reference Maric and Mas2020; Perales & Mas, Reference Perales and Mas2007; Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019). Notably, morning-phased or evening-phased genes show concomitant changes in the accumulation of transcriptionally permissive marks such as H3K4me3, H3K9ac, H3K27ac and H3S28p (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016; Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019). Among those genes are transcription factors, circadian regulators and components of the starch catabolic pathway (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016). By contrast, diurnal fluctuations in H3K27me1 or DNA methylation, associated with a transcriptionally repressive state, affect a different subset of genes, including cell cycle and DNA damage response genes (Balcerowicz, Reference Balcerowicz2024; Fung-Uceda et al., Reference Fung-Uceda, Gomez, Rodriguez-Casillas, Gonzalez-Gil and Gutierrez2024).
Spatial reorganization of the genome, changes in chromatin accessibility and epigenetic composition along the diurnal rhythm are likely involving general chromatin architects to either enable smooth transitions, stabilize new states, or both. Here, we explored the role of H1 linker histones (H1) in diurnal transcriptome fluctuations and discovered that they are crucial for coordinating the expression of genes with diurnal rhythms. Notably, genes that lose diurnal synchronization in the absence of H1 display binding sites for NAC transcription factors (TF), suggesting that H1 may influence their regulation by modulating accessibility to these TFs. Using large-scale quantitative imaging of nuclei isolated at different time points during the photoperiod, we observed subtle variations in nuclear size and chromatin distribution that were influenced by H1. We also found a decrease in global H3K4me3 and active RNA Pol II in the evening, accompanied by an increase in H3K27me3, with a role for H1 in H3K27me3 accumulation in the evening and RNA Pol II (Ser2P) in the late night and morning in the mutant.
Collectively, this study revealed a role for linker histone (H1) in regulating diurnal transcriptome and chromatin variations in Arabidopsis, focusing on how they might influence the coordinated expression of genes with rhythmic patterns. These findings pave the way to investigate the role of general chromatin architects in transcription and chromatin organization oscillations.
2. Results
2.1. H1 depletion affects transcriptional diurnal dynamics
The Arabidopsis genome encodes three H1 variants: H1.1, H1.2 and H1.3, the latter being stress-inducible and expressed only in guard cells under regular growth conditions (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015). Loss-of-function mutant lines were first obtained by expressing long interference RNAs (Wierzbicki & Jerzmanowski, Reference Wierzbicki and Jerzmanowski2005). Thereafter, T-DNA insertion lines were isolated and combined to create a 3h1 - triple homozygous mutant line (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019). This line was created by introgressing mutant alleles from different accessions (Ler and Col_0, respectively) and backcrossed four times to Col_0 (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015) possibly leaving 6.2% residual Ler genome (Collard et al., Reference Collard, Jahufer, Brouwer and Pang2005). While potential confounding effects in phenotypic analyses, arising from this small fraction of the hybrid genome, can be controlled using complementation approaches (Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019), having access to a triple mutant within a homogeneous genetic background would benefit future studies. We thus generated a new mutant line in a Col_0 background using CRISPR-Cas9 gene editing (Bieluszewski et al., Reference Bieluszewski, Szymanska-Lejman, Dziegielewski, Zhu and Ziolkowski2022). 3h1crispr harbours a small deletion in each of the three H1 genes, leading to frameshift mutations (Figure 1a, Supplementary File 1). Efficient depletion of the whole H1 complement was confirmed by the loss of detectable H1 in immunostained 3h1crispr nuclei (Figure 1b). The 3h1crispr knock-out also expectedly reproduces the loss of conspicuous chromocenters in which heterochromatin is compacted, early flowering phenotype under long days and increased lateral root number described previously (Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019) (Figure 1b, Supplementary File 1, and below). Consequently, all experiments in this study were conducted using this new mutant line, thereafter referred to as 3h1.
Figure 1. H1 depletion affects the diurnal control of gene expression. (a) Schematic representation of the CRISPR Cas9-mediated deletions in the H1.1, H1.2 and H1.3 coding regions, respectively (thick line, exon; thin line, intron). All deletions create frameshift mutations (see Supplementary File 1). (b) H1 Immunostaining (green) and DNA counterstaining (grey) of 3h1crispr leaf nuclei confirm H1 depletion and heterochromatin decondensation, respectively, as previously shown for the T-DNA-based 3h1 mutant (Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019). Scale bar: 2 μm (c). Representative images of 12-day-old wt and 3h1 mutant seedlings grown under continuous light or under long day conditions as indicated. (d) Number of lateral roots per seedling (dot) grown under continuous light (7-day-old seedlings, left) or under long day photoperiod conditions (12-day-old seedlings, right), for wild-type segregant seedlings (white) and 3h1 crispr mutant seedlings (grey). ***, P < 0.001, Kruskal–Wallis test with Bonferroni correction. The numbers of seedlings scored are indicated on the X axis. (e) Genes with co-regulated expression in wild-type (wt) seedlings were partitioned into six clusters (n, number of genes in each cluster). The graphs show their expression in 19 dag seedlings, wt (green traces) and mutant (3h1, pink traces), sampled at six time points as indicated. The overlay (right) allows visualizing genes losing synchronised, co-expression in the mutant (DSG, De-Synchronised Genes) in contrast to the bundled traces corresponding to genes with synchronised co-expression (SG, Synchronised Genes). (f) Deviation score -measuring the cumulative differences in expression of diurnally regulated genes in mutant seedlings across all time points – for all genes (left) and for the top and bottom 5% defining the DSG and SG classes, respectively. (g) Example of DSG expression in 3h1 and wt of the functional categories ‘response to heat (GO:0009408)’ and ‘regulation of innate immune response (GO:0045088)’. The values correspond to averaged, normalised expression levels. See also Supplementary Data 1 and Supplementary Figure S1.
Under continuous light, the 3h1 mutant displays an increased number of lateral roots, while overall root length remains comparable to the wild type (Supplementary Figure S1A). However, this excess of lateral roots is no longer apparent under a long-day photoperiod (Figure 1c–d, Supplementary Figure S1A, B). Notably, when both lateral root primordia and elongated lateral roots were scored, 3h1 still showed a higher initiation rate of lateral roots (Supplementary Figure S1B). These findings suggest that although 3h1 retains the capacity to initiate more lateral roots, their subsequent growth may be impaired under long-day conditions compared to under continuous light. The contrasting root behaviour – enhanced lateral root growth under continuous light but reduced growth under a long-day photoperiod – led us to hypothesize that H1 may play a broader role in regulating diurnal transcriptional dynamics. To test this, we analysed the transcriptome of wild-type (wt) and 3h1 mutant seedlings grown for 19 days under long-day (16h light–8h night) conditions, with the temperature fluctuating between 20°C at night and 22°C during the day. We profiled triplicate samples at six time points: ZT4, ZT13, ZT16, ZT20, ZT22, ZT24. These time points correlated with epigenomic changes at morning- or evening-phase loci (H3K9Ac at ZT4 and ZT16; H3K4me3 at ZT13) (data from Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019) and peak expression of chromatin modifiers such as JMJ14 at ZT20 and SDG2 at ZT22 (Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019). The replicate transcriptome profiles showed a high correlation at each time point (Supplementary Figure S1C). We then focused on genes with diurnal fluctuations and differential expressions between 3h1 and wt plants (LFC>1, FDR<0.01, n = 1040, Time and Treatment Contrast, TTC genes, Supplementary Data1). Among them, we identified 34 genes specifically upregulated during lateral root development (Gala et al., Reference Gala, Lanctot, Jean-Baptiste, Guiziou, Chu, Zemke and Nemhauser2021) (Supplementary Data1), which could explain the difference in lateral root number under continuous light. To gain a broader understanding of possible alteration in the diurnal kinetics of gene expression in 3h1, we clustered the previously selected TTC genes according to their kinetics in the wt and compared their diurnal expression profile in the mutant (Figure 1e, six clusters; Supplementary Figure S1D, eight clusters). In the wild type, gene clusters showed synchronised expression, whereas in the mutant, their expression was more variable and desynchronized for a subset of genes, peaking at different times compared to the wild type. In wild-type plants, H1-encoding genes do not exhibit diurnal changes in expression levels (Supplementary Figure S1E), suggesting that the impact of H1 depletion on diurnal gene regulation is not directly related to global H1 protein abundance. To further analyse genes most dramatically affected in their diurnal kinetics in the mutant, we used a deviation score that measures the cumulative discrepancy between wt and mutant profiles (Supplementary Figure S1F). Using this score, we selected genes with deviation scores in the top 5% range and bottom 5% range, respectively. This generated two gene groups, corresponding to de-synchronized genes (DSG, n = 197) and synchronized genes (SG, n = 773), respectively (Figure 1f).
2.2. H1 depletion mostly affects PRC2 targets enriched in NAC binding motifs
To determine whether DSG genes, which show altered diurnal kinetics upon H1 depletion, are clock-regulated, we verified their expression levels in conditions where seedlings were first entrained by a photo- and a thermoperiod before being transferred to continuous light (DIURNAL data resource, Mockler et al., Reference Mockler, Michael, Priest, Shen, Sullivan, Givan and Chory2007). The analysis showed that none of the DSGs are regulated by the clock, albeit two genes encoding subunits of the light-harvesting complex (Supplementary Figure S2A). Instead, about one-third of DSGs are light-responsive according to a recent analysis of Arabidopsis seedling development (Schivre et al., Reference Schivre, Wolff, Mirasole, Armanet, Davidson, Vidal and Barneche2025) (Supplementary Figure S2A).
Next, to answer the question whether DSG represent a random or specific subset of diurnally regulated processes, we analysed the GO terms associated with DSG and SG. Nearly 80% of the biological processes enriched in DSG were linked to responses to biotic and abiotic stresses, compared to only 10% in SG (Figure 2a). In contrast, SG were enriched in fundamental cellular metabolism and nuclear functions (Figure 2a). This suggests that H1 regulates the timing of stress-response gene expression. DSG were generally expressed at lower levels than stably expressed genes (SGs) across all time points (Figure 2b, Supplementary Figure S2B). This lower expression correlated with higher average levels of H3K27me3 at DSG, at least at day start, based on available data (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016; Figure 2c, Supplementary Figure S2C). While DSG showed similar average H1 enrichment to SGs (based on seedling data; time point not specified, Bourguet et al., Reference Bourguet, Picard, Yelagandula, Pelissier, Lorkovic, Feng and Mathieu2021, Figure 2c, Supplementary Figure S2C), H3K27me3 levels at DSG may still be influenced by H1 occupancy. Supporting this, approximately half of the DSG overlap with regions depleted of H3K27me3 in cotyledon tissues of the h1 double mutant (h1.1 h1.2) (Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023; Supplementary Figure S2C). Conversely, DSG show a significant depletion of H3K4me3, not only at day start (Figure 2c) but also across the diurnal phase (Supplementary Figure S2D, data from Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019). Since H1 deposition interplays with DNA methylation (Wierzbicki & Jerzmanowski, Reference Wierzbicki and Jerzmanowski2005; Zemach et al., Reference Zemach, Kim, Hsieh, Coleman-Derr, Eshed-Williams, Thao and Zilberman2013), we also asked whether DSG were distinct in their DNA methylation levels compared to SG. Using available whole genome bisulfite sequencing data from 2-week-old seedlings grown under a long-day photoperiod (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015), we found that DSG were, in average, 8x less methylated in CpG context than SG (Figure 2c). In the CHG and CHH contexts, DSG exhibit similar methylation levels than SG in average albeit for a few loci showing elevated mCHG levels (Supplementary Figure S2E). Yet, DNA methylation levels do not change significantly more in DSG compared to SG in a 3h1 mutant background (Supplementary Figure S2E). This suggests that DNA methylation may play little role in the deregulation of DSG diurnal expression in 3h1 mutant seedlings.
Figure 2. DSG are mostly response genes likely regulated by NAC transcription factors. (a) Genome ontology analysis of DSG and SG showing biological processes that are significantly enriched (P < 0.05). Processes are grouped in major categories indicated in the legend. The X axis shows the % of genes falling into each category. (b) Averaged, normalised expression levels of SG and DSG in the wild-type at each time point. (c) Relative enrichment levels of H3K27me3 at day start (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016), H3K4me3 at ZT3 (Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019), of H1 (Bourguet et al., Reference Bourguet, Picard, Yelagandula, Pelissier, Lorkovic, Feng and Mathieu2021) and mCpG (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015) at unknown time points, in DSG vs SG. P values are from a Kruskal–Wallis test with Bonferroni correction. (d) Top two most enriched motifs among DSG and their occurrence in % among DSG vs SG (pie charts). See also Supplementary Figure S2 and Supplementary Data 2.
Next, to investigate whether DSG regulation shares a common control mechanism, we searched for binding motifs of transcriptional regulators. We found a significant recurrence of motifs corresponding to NAC protein binding sites (Figure 2d, Supplementary Data 2). NAC proteins (NO APICAL MERISTEM, ATAF1/2 and CUP-SHAPED COTYLEDON) are part of a large, ancient, plant-specific transcription factor family involved in plant growth, senescence, flowering, fruit ripening, immunity and abiotic stress responses (Kim et al., Reference Kim, Alvarez, Rana, Preciado, Liu and Begcy2024) (Olsen et al., Reference Olsen, Ernst, Leggio and Skriver2005). A role for NAC transcription factors in diurnal regulation has not been reported so far. While some of them can be diurnally or clock-regulated (Ma et al., Reference Ma, Zhou, Ma, Tan, Huang, Segbo and Gao2025; Peng et al., Reference Peng, Phung, Zhai and Neff2020), it is not the case for our candidates whose binding motifs are enriched in DSGs (Supplementary Figure S2E). Whether NAC TFs interact with or compete against H1 proteins at binding sites to influence DSG rhythmic transcription remains to be determined. Remarkably, DSG did not contain binding sites for circadian clock regulators, further supporting the idea that their rhythmic expression may be independently regulated or only indirectly controlled by the clock.
2.3. H1 contributes to diurnal chromatin dynamics
H1 is broadly distributed across the genome (Bourguet et al., Reference Bourguet, Picard, Yelagandula, Pelissier, Lorkovic, Feng and Mathieu2021; Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015) and acts as a general chromatin organizer (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020; He et al., Reference He, Yu, Wang, Zhang, Bai, Li and Feng2024; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019; Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023). Despite this broad role, it can influence the transcriptional timing of specific genes (results above). We therefore asked whether H1 contributes to large-scale chromatin reorganization that could influence or accompany transcriptional reprogramming.
H1 visibly affects chromatin organization at the cytological level (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020; He et al., Reference He, Yu, Wang, Zhang, Bai, Li and Feng2024; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019) and 3D genome topology (Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023). Using quantitative imaging, we measured several nuclear and chromatin features. Nuclei were isolated from wild-type and 3h1 mutant seedlings at the same time points as before in replicate experiments (see Methods and Supplementary Data 3). In total, ~1530 DNA-stained nuclei were imaged (Supplementary Figure S3A). Using Nucl.Eye.D (Johann To Berens et al., Reference Johann To Berens, Schivre, Theune, Peter, Sall, Mutterer and Molinier2022), we performed automated segmentation to quantify nuclear size (area), shape (roundness, circularity, solidity, aspect ratio) and DNA intensity distribution. While analysing the mean values of these features indicated variation across time points (Supplementary Figure S3B, C), a principal component analysis (PCA) best captured the diurnal changes, showing distinct trajectories during the day and night in wt nuclei (Figure 3a). Clearly, 3h1 mutant nuclei occupied a separate region in the PC space, likely due to their larger size and reduced heterochromatin content (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019), which also affects DNA intensity distribution (Supplementary Figure S3A). Yet, diurnal cytological changes showed a different trajectory in the mutant, indicating that H1 plays a role in chromatin reorganization dynamics during the photoperiod (Figure 3a). We measured specifically the relative heterochromatin fraction (RHF) present in conspicuous chromocenters (Fransz et al., Reference Fransz, De Jong, Lysak, Castiglione and Schubert2002). We observed that the RHF remained constant on average across time points and as expected (Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019), the 3h1 mutant showed lower RHF and fewer chromocenters (Figure 3b and Supplementary Figure S3D).
Figure 3. Diurnal changes in nuclear size, chromatin density, global epigenetic and transcriptional levels in Arabidopsis seedlings. Seedling nuclei were analysed by quantitative imaging following isolation, immunostaining and DNA counterstaining, before batch-segmentation using Nucl.Eye.D (Johann To Berens et al., Reference Johann To Berens, Schivre, Theune, Peter, Sall, Mutterer and Molinier2022). (a) PCA analysis of nuclear features, including area, shape descriptors (roundness, circularity, solidity, aspect ratio) and DNA density distribution (intensity mean and standard deviation -std). The distinct PC space occupied by wt and 3h1 samples, as well as distinct trajectories along the day-night cycle, were manually coloured and indicated. PC plot: columns with similar annotations were collapsed by taking the median inside each group. Unit variance scaling was applied to rows and SVD with imputation was used to calculate principal components (PC) 1 and 2 (X and Y axes) that explain 78% and 11% of the total variance, respectively. The weights of each descriptor for PC1 and PC2 are shown on the right. The number of nuclei analysed is given in Supplementary Figure S3A and ranges between 40 and 73 per sample. (b–e) Quantification of the relative heterochromatin content (RHF, B) and relative levels of H3K27me3 (c), H3K4me3 (d) and RNA Pol II Ser2P (e) in wt (blue) and 3h1 (red) nuclei. Relative levels: for each nucleus, the immunostaining signal intensity was normalised by the DNA signal intensity. Differences between genotypes were evaluated using a Mann–Whitney U test. ns, not significant. *, P < 0.05. **, P < 0.01. ***, P < 0.001. ****, P < 0.0001. The right panels show representative images for each experiment: a DNA-stained nucleus segmented by nucl.Eye.D (B, green contours), nuclei immunostained (green) and counterstained with propidium iodide (magenta). Scale bar: 2μm. See also Supplementary Figure S3, Supplementary Data 3 and BioStudies accession S-BIAD1915.
Despite its strong influence on the epigenetic landscape, H1 depletion has a modest effect on transcription in Arabidopsis (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020; Hsieh et al., Reference Hsieh, He, Buttress, Gao, Couchman, Fischer and Feng2016; Lyons & Zilberman, Reference Lyons and Zilberman2017; Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019; Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023; Wierzbicki & Jerzmanowski, Reference Wierzbicki and Jerzmanowski2005; Zemach et al., Reference Zemach, Kim, Hsieh, Coleman-Derr, Eshed-Williams, Thao and Zilberman2013). Since the photoperiod is associated with rhythmic shifts in histone marks such as H3K4me3, H3K9Ac, H3K27Ac, H3K27me3 and H3K27me1 at diurnally regulated loci (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016; Fung-Uceda et al., Reference Fung-Uceda, Gomez, Rodriguez-Casillas, Gonzalez-Gil and Gutierrez2024; Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019), we asked whether H1 might contribute to some of these changes. To test this, we used quantitative immunostaining for measuring the nuclear abundance of H3K4me3 and H3K27me3 – marking transcriptionally permissive and repressive states, respectively – as well as productive RNA Polymerase II (Ser2-phosphorylated, Ser2P). In nuclei of wild-type seedlings, all three markers showed diurnal variations (Supplementary Figure S3E–G). Notably, H3K27me3 levels peaked at the end of the day (ZT16) and decreased before the end of the night (ZT24). Interestingly, a subset of nuclei at ZT16 and ZT20 showed a nearly 4-fold increase compared to ZT4, suggesting that the amplitude of this response might be cell-type specific. The 3h1 mutant also showed diurnal fluctuations in H3K27me3 but with less sharp transitions than in the wt, leading to significantly elevated levels towards the end of the day (ZT13) and end of the night (ZT22) (Figure 3c).
Conversely, both H3K4me3 and RNA Pol II Ser2P reached their lowest levels at ZT16-ZT20, rising slowly during the night until ZT4 in wt plants (Supplementary Figure S3F, G). The reduction in productive RNA Pol II (Ser 2P) at the end of the day correlated with a decrease in paused RNA Pol II (Ser 5P) (Supplementary Figure S3H), suggesting a global reduction in transcriptional activity as the night approaches. In comparison, 3h1 nuclei also showed diurnal fluctuations, but – as for H3K27me3 - with a strong attenuation of the diurnal variations. This led to higher levels of productive RNA Pol II (Ser2P) towards the end of the night (ZT22, ZT24) and lasting until the morning (ZT4) in conjunction with slightly elevated levels of H3K4me3 (Figure 3d, e).
To determine whether the increased RNA Pol II activity in 3h1 was a result of a reduction in the inactive fraction of RNA Pol II or reflected a general increase in available RNA Pol II, we quantified RNA Pol II Ser5P and total RNA Pol II (NP) levels, respectively, by immunostaining as before and from ZT4 to ZT16. The results supported the latter scenario (Supplementary Figure 1).
In conclusion, our quantitative immunostaining analyses revealed large-scale, diurnal changes in chromatin organization and epigenetic composition, with fluctuations in global levels of H3K4me3 and H3K27me3 in wild-type seedlings. These changes are coincidental with a global reduction in RNA Pol II toward the end of the day, followed by an increase in the morning. In this process, linker histones appear to play a role, particularly in suppressing transcriptional activity during the night and early morning.
3. Discussion
H1 linker histones are essential components of chromatin. They regulate its structure, compaction and epigenetic landscape in both plants and animals (Fyodorov et al., Reference Fyodorov, Zhou, Skoultchi and Bai2018; He et al., Reference He, Yu, Wang, Zhang, Bai, Li and Feng2024; Over & Michaels, Reference Over and Michaels2014; Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019; Zemach et al., Reference Zemach, Kim, Hsieh, Coleman-Derr, Eshed-Williams, Thao and Zilberman2013). In Arabidopsis, H1 depletion causes mild developmental changes (Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019). Our study revealed that H1 plays a role in synchronizing the expression of a subset of rhythmically regulated genes, particularly those involved in stress responses. This aligns with several studies highlighting the involvement of H1 variants in responses to biotic and abiotic stresses in Arabidopsis. For example, compared to wild-type plants, 3h1T-DNA mutant seedlings exhibit enhanced resistance to Pseudomonas syringae and Botrytis cinerea, as well as increased defence priming in response to flg22 (Sheikh et al., Reference Sheikh, Nawaz, Tabassum, Almeida-Trapp, Mariappan, Alhoraibi and Hirt2023). They also show greater resistance to heat stress, particularly when combined with DNA methylation inhibition (Liu et al., Reference Liu, de Jonge, Trejo-Arellano, Santos-Gonzalez, Kohler and Hennig2021), but display hypersensitivity to salt stress at germination (Perrella et al., Reference Perrella, Fasano, Donald, Daddiego, Fang, Martignago and Amtmann2024). Additionally, plants lacking the H1.3 variant demonstrate increased tolerance to drought and low-light conditions (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015) and seedlings with altered H1.3 expression show slightly enhanced resistance to salt stress (Wu et al., Reference Wu, Xu, Meng, Fang, Xia, Zhang and Fan2022). Stress responses are subject to natural variation in Arabidopsis (Miller et al., Reference Miller, Song, Shi, Juenger and Chen2015; Zhao et al., Reference Zhao, Fu, Wang, Zhang, Liang, Cui and Liu2024). Whether the residual Ler/Col hybrid genomic fraction in the 3h1T-DNA mutant contributes – albeit possibly modestly – to the altered stress responses reported remains an open question. This issue could now be addressed using the 3h1crispr mutant in the Col-0 background in future studies. Nonetheless, our finding that diurnally regulated biotic and abiotic stress-response genes show altered rhythmic expression in the 3h1crispr mutant strongly supports a central role for H1 in fine-tuning biotic and abiotic stress responses.
While it is already known that a significant fraction of diurnally regulated genes are associated with shifts in various histone marks (Baerenfaller et al., Reference Baerenfaller, Shu, Hirsch-Hoffmann, Futterer, Opitz, Rehrauer and Gruissem2016; Fung-Uceda et al., Reference Fung-Uceda, Gomez, Rodriguez-Casillas, Gonzalez-Gil and Gutierrez2024; Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019), our findings extend this knowledge by identifying rhythmic changes beyond these loci. Specifically, we observed changes in chromatin organization and nuclear levels of the key epigenetic marks H3K27me3 and H3K4me3, as well as in RNA Pol II levels. Additionally, we found that H1 contributes to these global changes. Given that H1 interacts closely with DNA methylation (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020; Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015; Wierzbicki & Jerzmanowski, Reference Wierzbicki and Jerzmanowski2005; Zemach et al., Reference Zemach, Kim, Hsieh, Coleman-Derr, Eshed-Williams, Thao and Zilberman2013), it remains to be determined whether global DNA methylation levels also fluctuate in a diurnal manner.
Several mechanisms, potentially cooperating, could explain the disrupted transcription timing in the absence of H1: a direct effect on RNA Pol II progression, an indirect effect through modulation of H3K4me3, H3K27me3 or other marks accumulation, such as DNA methylation, or a combination thereof. Indeed, H1 regulates RNA Polymerase II (RNA Pol II) progression during key phases of transcription. In both yeast and mammals, linker histones cause RNA Pol II to pause near the nucleosome entry site (Hirano et al., Reference Hirano, Ehara, Kujirai, Uejima, Takizawa, Sekine and Kurumizaka2022). This is due to the close spatial proximity of the entry and exit DNA strands constrained by H1 on the nucleosome that creates a steric clash between the DNA exit strand and RNA Pol II, which is engaged with the entry strand. This conflict is resolved by transcription elongation factors, which remove H1 from the entry strand (Hirano et al., Reference Hirano, Ehara, Kujirai, Uejima, Takizawa, Sekine and Kurumizaka2022). In Arabidopsis, the absence of H1 leads to irregular nucleosome spacing, with increased occurrences of both short (<170bp) and long (>300nt) nucleosome repeat lengths (NRL) (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019). Short nucleosome spacing may intensify steric clashes with RNA Pol II, which elongation factors cannot efficiently resolve. Conversely, longer NRL regions may promote extended RNA Pol II progression, where H1 would normally limit it. Thus, H1 absence could disrupt the regulation of transcriptional elongation during diurnal fluctuations. Also, intergenic transcripts, which emerge in the absence of H1 (Choi et al., Reference Choi, Lyons, Kim, Moore and Zilberman2020), might contribute to transcriptional timing. Additionally, in Arabidopsis, H1 is known to influence the deposition and maintenance of DNA methylation as well as histone marks such as H3K4me3 and H3K27me3 (Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015; Rutowicz et al., Reference Rutowicz, Lirski, Mermaz, Teano, Schubert, Mestiri and Baroux2019; Wierzbicki & Jerzmanowski, Reference Wierzbicki and Jerzmanowski2005; Zemach et al., Reference Zemach, Kim, Hsieh, Coleman-Derr, Eshed-Williams, Thao and Zilberman2013). Hence, the disrupted transcriptional timing observed in H1-deficient plants may be linked to altered epigenetic states that are poorly regulated in the absence of H1. Supporting this hypothesis, De-Synchronised Genes (DSG) in the 3h1 mutant are largely enriched in H3K27me3 and half display reduced H3K27me3 level upon H1 depletion, while their expression increases in the absence of H1. DNA methylation on the contrary may not contribute much to DSG regulation since mCpG levels are much lower at DSG loci than SG and barely change upon H1 depletion. However, this alone does not fully explain changes in transcription kinetics. Temporally resolved, diurnal profiles of H3K27me3 and DNA methylation in both wild-type and 3h1 mutant plants are needed. Likely, however, additional factors may act in concert with H1 to integrate epigenetic states and enable diurnal transcriptional dynamics.
Linker histones (H1) are not essential for plant growth or viability under controlled conditions. However, under stress – both biotic and abiotic – H1 plays a role in modulating the stress response (Liu et al., Reference Liu, de Jonge, Trejo-Arellano, Santos-Gonzalez, Kohler and Hennig2021; Perrella et al., Reference Perrella, Fasano, Donald, Daddiego, Fang, Martignago and Amtmann2024; Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015; Sheikh et al., Reference Sheikh, Nawaz, Tabassum, Almeida-Trapp, Mariappan, Alhoraibi and Hirt2023), influencing physiological and growth adaptations possibly through their influence on RNA Pol II and/or PRC2 activity (Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023). It is interesting to note that genes whose diurnal regulation is disrupted upon H1 depletion are primarily involved in biotic and abiotic stress responses, heavily marked by H3K27me3 during 24-h cycles and not regulated by the circadian clock. Clock-controlled genes, in contrast, are not subject to H3K27me3 deposition (Malapeira et al., Reference Malapeira, Khaitova and Mas2012) and often relate to core metabolic processes, light perception and growth regulation, although an interplay with abiotic responses is identified (Panter et al., Reference Panter, Muranaka, Cuitun-Coronado, Graham, Yochikawa, Kudoh and Dodd2019). The specificity of functional categories impacted by the absence of H1 is intriguing, especially considering the widespread deposition of H1 variants across the genome (Bourguet et al., Reference Bourguet, Picard, Yelagandula, Pelissier, Lorkovic, Feng and Mathieu2021; Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015). One plausible mechanism enabling selective effects of H1 at DSG is a competitive interaction between H1 and other transcriptional regulators, such as NAC transcription factors, which control stress responses and other processes (Olsen et al., Reference Olsen, Ernst, Leggio and Skriver2005). This model is supported by the observation that NAC binding motifs are enriched at loci with desynchronized expression in the 3h1 mutant. In this context, H1 may restrict NAC factor access to these motifs. However, NAC target genes are not simply overexpressed throughout the photoperiod in the absence of H1 and their corresponding candidate NAC transcription factors are themselves not diurnally regulated. This suggests that other factors modulate the competitive binding between H1 and NAC in a diurnal manner. One possibility is the role of H1 post-translational modifications (PTMs), which influence H1 binding kinetics (Catez et al., Reference Catez, Ueda and Bustin2006). Phosphorylation is a prime candidate for destabilizing H1 from chromatin, which is especially relevant considering the significant fluctuations in the Arabidopsis phospho-proteome during the photoperiod (Uhrig et al., Reference Uhrig, Echevarria-Zomeno, Schlapfer, Grossmann, Roschitzki, Koerber and Gruissem2021). Competitive binding between H1 and specific chromatin regulators has been demonstrated, which explains the impact of H1 depletion at TRB target loci (Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023). Thus, this model of competitive binding suggests that H1, a generic chromatin factor, could regulate specific subsets of loci at certain times, potentially in response to environmental or developmental signals.
In conclusion, H1 may facilitate transcriptional reprogramming by coordinating chromatin accessibility, RNA Pol II progression and competition with transcriptional regulators. However, the mechanisms by which these processes achieve precise temporal control at specific loci remain unclear and further research is required to elucidate these questions. Photoreceptors, metabolic feedback loops, or both may directly or indirectly influence the interaction between H1, epigenetic marks, NAC transcription factors and RNA Pol II, ultimately synchronizing the expression of gene subsets with similar transcriptional kinetics.
4. Material and Methods
4.1. Plant growth
The Arabidopsis mutant line 3h1(T-DNA) was obtained from Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015. The newly created 3h1(CRISPR) is described in Supplementary File 1, along with the new and previously reported single H1 mutant alleles. Seedlings were grown on 0.5x Murashige & Skoog (MS) salt base (Carolina Biological Supply, US) with 1% PhytoAgar (Duchefa Biochemie, Netherlands) under continuous light or long-day conditions (16h/8h), depending on the experiment. The temperature was maintained at 21-22°C during the day and 18°C at night. For transcriptomic experiments, all seedlings were grown simultaneously in the same growth chamber, with six replicate plates per genotype and time point, randomly distributed on the shelves for sample collection. For cytological experiments, seedlings were grown in two separate incubators with an inverted day-night cycle to allow tissue collection for all time points within a day. Samples were collected in two replicate experiments (A and B) per genotype and time point.
4.2. Quantification of root length and lateral root number
Seedlings were grown as before, but vertically on square petri dishes and scanned 7 and 12 days after germination, before measurements of root length using Fiji (Schindelin et al., Reference Schindelin, Arganda-Carreras, Frise, Kaynig, Longair, Pietzsch and Cardona2012) and manual scoring of lateral roots.
4.3. RNA-seq data processing
Short reads generated in this study were deposited at NCBI Sequence Read Archive (SRA, www.ncbi.nlm.nih.gov/sra) and are accessible through the accession number PRJNA1244618. Reads were trimmed and quality-checked with fastp (Chen et al., Reference Chen, Zhou, Chen and Gu2018) (version 0.20.1). Transcripts were quantified with Salmon (Patro et al., Reference Patro, Duggal, Love, Irizarry and Kingsford2017) (version 1.4.0) using the cDNA and gene annotation available from araport.org (Araport 11). Variation in read counts was analysed with a general linear model in R (version 3.6.1) with the package DESeq2 (Love et al., Reference Love, Huber and Anders2014) (version 1.24.0) according to a factorial design with the two explanatory factors ‘time’ and ‘treatment’ combined into a single factor. Specific conditions were then compared with linear contrasts. Within each comparison, P-values were adjusted for multiple testing (Benjamini-Hochberg). Regions with an adjusted P-value (false discovery rate, FDR) below 0.01 and a minimal log2 fold-change (i.e., the difference between the log2-transformed, normalized sequence counts) of 1 were differentially expressed. Normalized sequence counts were calculated accordingly to DESeq2 and log2(x+1) transformed.
4.4. Definition of gene sets
To identify genes with differences in diurnal expression between the 3h1 mutant and the wild-type, we used 1) linear, quadratic and cubic contrast (Rosenthal & Rosnow, Reference Rosenthal and Rosnow1985); 2) contrasts that compare one time point to all other time points (one test for each other time point); and 3) contrasts that compare two adjacent time points to all other time points (also one test for each other time point). Genes that showed significant differences in any of these contrasts using the wild-type data were considered ‘genes that change over time’.
To identify genes with differences in expression between 3h1 and wild-type, the treatments were compared at each and across all time points. Genes that showed significant differences in any of these comparisons were defined as ‘genes that change by treatment’. However, to focus on genes that show the most ‘de-synchronized’ pattern, we calculated an index of change as the sum of the absolute LFCs between 3h1 and wt at all individual time points. Large values indicate that genes show a 3h1-specific pattern that is clearly distinct from the wild-type and not just shifted in intensity. From all genes, the top 2.5% were intersected with the genes that changed by treatment, resulting in 216 de-synchronized genes. Likewise, we devised a set of genes that were synchronized between 3h1 and wt seedlings by intersecting the bottom 2.5% with the genes that do not show any differences in expression between treatments (578 synchronized genes).
4.5. Time course visualization
To visualize the time-course behaviour of genes of interest, we first grouped genes into clusters using wild-type data. For this, replicates from each time and treatment combination were averaged. Expression of each gene was standardized and the wild-type data were clustered into 4–8 clusters using Mfuzz (Kumar, Reference Kumar and Futschik2007). 3h1 data were added to the visualization but not used while clustering. Only genes belonging to the cluster core are shown. Those were extracted with the function acore() from the Mfuzz package (Kumar, Reference Kumar and Futschik2007).
4.6. Functional gene set characterization
To functionally characterize a gene set of interest (e.g., synchronized and de-synchronized genes), we tested for enrichment of gene ontology (GO) terms with topGO (version 2.32.0) (Alexa et al., Reference Alexa, Rahnenfuhrer and Lengauer2006). Analysis was based on gene counts (genes in the set of interest compared to all annotated genes) using the ‘weight’ algorithm with Fisher’s exact test (both implemented in topGO). A term was identified as significant if the P-value was below 0.05.
4.7. Motif discovery
To search for motifs in promoters (2 kb upstream of TSS) of genes of interest, we used findMotifs.pl from the HOMER suite (version 4.11) (Heinz et al., Reference Heinz, Benner, Spann, Bertolino, Lin, Laslo and Glass2010) with the synchronized genes as the control group. As we observed many NAC matches, we extracted all matches and generated a consensus motif with the R-package ‘universalmotif’ (https://github.com/bjmt/universalmotif). The consensus motif was then searched in the promoter sequences with FIMO (MEME Suite, version 5.0.2, 2011) (Bailey et al., Reference Bailey, Boden, Buske, Frith, Grant, Clementi and Noble2009).
4.8. Cross-analyses with ChIP-seq and RNA-seq data
Publicly available ChIP-seq raw data were retrieved from SRA/NCBI. Datasets from seedling tissues included H3K4me3 and H3K9me3 (Song et al., Reference Song, Huang, Yu, Ando, Mas, Ha and Chen2019), H3K27me3 (Veluchamy et al., Reference Veluchamy, Jegu, Ariel, Latrasse, Mariappan, Kim and Benhamed2016) and H1 (Bourguet et al., Reference Bourguet, Picard, Yelagandula, Pelissier, Lorkovic, Feng and Mathieu2021). Reads were trimmed and quality-checked with fastp (version 0.20.1) (Chen et al., Reference Chen, Zhou, Chen and Gu2018) and aligned to the Arabidopsis reference genome (TAIR10) with bowtie 2 (version 2.3.5.1) (Langmead & Salzberg, Reference Langmead and Salzberg2012) keeping unique alignments. Duplicate reads were removed with Picard (version 1.140; broadinstitute.github.io/picard/). Coverage was normalized and visualized with DeepTools (version 3.3.0) (Ramirez et al., Reference Ramirez, Dundar, Diehl, Gruning and Manke2014); bamCoverage with --normalizeUsing CPM, computeMatrix scale-regions and plotHeatmap). In addition, determination of differential H3K27me3 enrichment in wt and 2h1 seedlings for DSG and SG was done according to a DESeq2 analysis reported in [Teano et al., Reference Teano, Concia, Wolff, Carron, Biocanin, Adamusova and Barneche2023, Supplementary Table S1). The values reported in Supplementary Table 2 represent Log2 fold change, with NA corresponding to genes without significant H3K27me3 enrichment upon peak detection with MACS2. The cross-analysis of DSG and SG with light-induced genes during cotyledon de-etiolation was done using a Spike-In RNAseq dataset (Schivre et al., Reference Schivre, Wolff, Mirasole, Armanet, Davidson, Vidal and Barneche2025). Values represent normalised transcript levels, with NA corresponding to genes with no detectable reads.
4.9. WGBS data processing
Short reads from whole genome bisulfite sequencing for Col_0 wt seedlings and 3h1 (T-DNA allele, Rutowicz et al., Reference Rutowicz, Puzio, Halibart-Puzio, Lirski, Kotlinski, Kroten and Jerzmanowski2015) mutant seedlings were obtained from ArrayExpress (E-MTAB-2807). Reads were (quality) trimmed with fastp (version 0.20.0 with the options – trim_front1 2, (Chen et al., Reference Chen, Zhou, Chen and Gu2018) and aligned to the reference genome with Bowtie2 (version 2.3.5.1, Langmead and Salzberg (Reference Langmead and Salzberg2012)) in combination with Bismark (version 0.22.3 with the option – non_directional (Krueger & Andrews, Reference Krueger and Andrews2011)). The Col-0 reference genome was obtained from www.arabidopsis.org (TAIR10). Alignments were deduplicated with bismark and methylation tables were extracted with MethylDackel (version 0.5.0, github.com/dpryan79/MethylDackel). Files were merged and only cytosines with a coverage between 5 and 100 within all four samples were kept. 11.1 million cytosines passed this filter. Cytosines were mapped to the Araport 11 annotation to obtain average methylation levels per gene and methylation context. For each gene and context, we averaged the methylation levels for the two groups (wt and 3h1) and calculated the difference between the groups (delta). Genes that were identified as SG and DSG were extracted and tested within each context for different DNA methylation levels in wt, comparing SG and DSG groups, or differential response in 3h1 vs wt (delta) using a two-sided t-test.
To assess differences in DNA-methylation in more detail, the two groups (3h1 and wt) were compared to each other with the R-package DSS (version 2.32.0, (Park & Wu, Reference Park and Wu2016) using the functions DMLfit.multiFactor() and DMLtest.multiFactor(). P-values were corrected for multiple testing and thereby converted to false discovery rates (FDRs, Benjamini and Hochberg (Reference Benjamini and Hochberg1995)). Cytosines were identified as significantly differentially methylated (DMC) if the FDR was below 0.05. 25,218 cytosines exhibited significant differential methylation (0.23%). We extracted the genes for which there were at least 10 cytosines in one of the contexts tested (34,626) and subset them to the ones with at least 20% differentially methylated cytosines: 468 genes were left but only 10 of them were in the list of SGs and DSGs. However, all of them were SGs. The genes were AT1G31440, AT2G03350, AT2G04540, AT2G04880, AT2G14120, AT2G14680, AT2G17410, AT2G18876, AT2G19430 and AT2G19470.
4.10. Nuclei extraction and immunostaining
Nuclei were extracted essentially as described previously (Kracik-Dyer & Baroux, Reference Kracik-Dyer and Baroux2025), leaving out the embedding step and adding a clearing procedure. In brief, leaves from 19-day-old A. thaliana seedlings were fixed, washed then finely chopped with a sterile razor blade in nuclear isolation buffer. The suspension was briefly sonicated and filtered before being distributed as 5 μL drops on cleaned slides (EprediaTm SuperFrost PlusTM Adhesion Slides, ThermoFisher Scientific). Nuclei were air-dried until used or stored at 4°C. Prior to immunostaining, slides underwent a clearing treatment to enhance downstream imaging quality. Specifically, slides were placed in a Coplin jar and treated for 5 minutes in a 1:1 solution of 100% ethanol and xylene, followed by a 5-min wash in 100% methanol, a 5-min wash in a 1:1 solution of methanol and 1× PBS and a final 10-minute wash in 1× PBS.
Immunostaining was carried out following a 45-min blocking step at room temperature using 4% BSA in 1×PBS,143 mM NaCl, followed by a 10-min wash in 1× PBS. The primary and secondary antibodies were each applied sequentially at a 1:1000 dilution after having verified that this allowed a quantitative analysis (Supplementary Figure S3J) in 1× PBS containing 1% BSA, 143 mM NaCl, 0.05% Tween20 and incubated at 37°C for 30 min each, with a 10-min wash in 1× PBS between applications. Antibodies used in this study: Rabbit-anti-H3K27me3 (Active Motif cat# 39155), Rabbit-anti-H3K4me3 (Abcam cat# ab8580), Rabbit-anti-RNA Polymerase II CTD Phospho S2 (Abcam cat# ab5095), Rabbit-anti-RNA Polymerase II CTD Phospho S5 (Abcam cat# ab5131), Goat-anti-Rabbit Alexa 488 conjugate (ThermoFisher Scientific, cat# A-11008). Nuclei were stained with 2 μg/mL 4′,6-diamidino-2-phenylindole (DAPI) or 2 μg/mL Propidium Iodide (PI) in PBS for 10 min, washed 1 min in PBS and mounted in VECTASHIELD® Antifade Mounting Medium complemented with DAPI (Vector Laboratories H-1100) or PI (Vector Laboratories H-1300), respectively, before being covered with 18 mm × 18 mm, 0.17 mm thick cover glasses, sealed with transparent nail polish and stored in the dark at 4°C until imaging.
4.11. Fluorescence microscopy imaging
To minimize selection bias, nuclei were consistently selected based on DNA staining with DAPI rather than antibody staining. Nuclei that appeared damaged or morphologically abnormal were excluded from imaging. Images were acquired using a Leica DM6000B microscope (Leica microsystems GmbH, Germany) equipped with a 63×/1.4 NA oil immersion objective and an Andor Neo 5.5 sCMOS camera (Andor, Oxford Scientific Instruments, USA). Images were captured in a 16-bit format using the ‘low-noise gain’ setting, and 15-200ms exposure, depending on the antibody and fluorescence filters (Chroma Technology Corp, USA: A4, L5 and Tx2 for DAPI, Alexa488 and PI detection, respectively). For quantifying immunostaining signal, single-plane images were acquired, while short z-stacks were acquired – using system’s optimized settings – for downstream quantification of the relative heterochromatin fraction (RHF). Acquisition parameters were set for each fluorochrome to achieve a good pixel intensity distribution while avoiding saturation and kept constant throughout the slides.
4.12. Image analysis
For RHF analyses, images were segmented and quantified using Nucl.eye.D using the author’s instructions (Johann To Berens et al., Reference Johann To Berens, Schivre, Theune, Peter, Sall, Mutterer and Molinier2022) and using 2D projections of the z-stacks as input images. Max projections were created in batch using a macro in Fiji (Schindelin et al., Reference Schindelin, Arganda-Carreras, Frise, Kaynig, Longair, Pietzsch and Cardona2012). The quantification of immunostaining signals was carried out using Fiji (Schindelin et al., Reference Schindelin, Arganda-Carreras, Frise, Kaynig, Longair, Pietzsch and Cardona2012). For this, a macro was written to (i) separate the channels for each image (Image/Colour/Split Channels), (ii) stack all images in one file (Image/Stacks/Images to Stack/Copy (centre); Use Titles as Labels), create one file per channel, (iii) create a montage for each stack file (Image/Stacks/Make Montage using Increment = 2, 6 Columns, 3 Rows – adjust to needs, First slice = 1 for the first montage, first channel; First slice = 2 for the second montage, second channel). Then, Regions of Interest (ROIs) were manually drawn around each nucleus in the montage containing the DNA signal, with ROIs saved in the ROI manager and used to extract signal intensities. The saved ROIs were opened and applied to the other montage to extract signal intensities from the second channel. Signal intensities were expressed as a ratio (antibody to DNA) for creating the graphs, plotted with RStudio (RStudio [2022.12.0+353] http://www.rstudio.com/). Unless otherwise stated, the data distribution across timepoints and genotypes were statistically evaluated using an ANOVA with Tukey HSD test. PCA analyses were conducted with the online tool ClustVis (Metsalu & Vilo, Reference Metsalu and Vilo2015).
Open peer review
To view the open peer review materials for this article, please visit http://doi.org/10.1017/qpb.2025.10033.
Acknowledgements
The authors thank Valeria Gagliardini, Christoph Eichenberger, Arturo Bolanos and Peter Kopf (University of Zurich) for general lab support, Professor Ueli Grossniklaus (University of Zürich) for valuable discussions.
Competing interest
The authors declare no conflict of interest.
Data availability statement
RNAseq data are available at NCBI Sequence Read Archive (SRA) with the reference number PRJNA1244618. Image data are available at Biostudies accession number S-BIAD1915.
Author contributions
KR and CBa conceived and designed the study, PZ and AJ designed elements of the study. KR, LP, JS, DZ-P, TB, MK and RSR performed the experiments. MWM performed bioinformatic analyses. KR, LP, CBo, FB and CBa performed data analyses and statistical analyses. CB wrote the manuscript. KR, LP, PZ, TB, MK, CBo and FB contributed manuscript amelioration. AI was used for language and grammar corrections.
Funding statement
The project has received funding from the Swiss National Science Foundation (SNSF Grant #310030_185186 and IZCOZ0_182949 to CBa), University of Zürich (Stiftung für Wissenschaftliche Forschung #STWF-19-023 to CBa), the National Science Center, Poland (UMO-2018/31/D/NZ2/02974 to MK), the Velux Foundation (Project 1107 to CBa and FB), from Agence Nationale de la Recherche (ANR-20-CE13-0028 to CBo, ANR-24-CE12-1113-01, ANR-24-CE20-2108 and ANR-22-CE20-0001 to FB) and by the CNRS program EpiPlant, France.
Supplementary material
The supplementary material for this article can be found at http://doi.org/10.1017/qpb.2025.10033.
Open access
Comments
No accompanying comment.