Introduction

Eukaryotic genomes are compartmentalized into euchromatic and heterochromatin regions (Ruiz-Velasco and Zaugg, 2017). In euchromatin, nucleosomes are more accessible, while in heterochromatin nucleosomes are more compacted and restrictive to transcription. Protein coding genes are typically euchromatic, while non-coding and repetitive elements more often reside in heterochromatin. However, transposable elements, which are targeted for DNA methylation and silencing, are found in both euchromatic and heterochromatic regions of the genome (Bourque et al., 2018). The pathways responsible for deposition of methylation and silencing of TEs in these different regions are functionally and genetically distinct (Du et al., 2015).

In plants, euchromatic repetitive regions are targeted by the RNA-directed DNA methylation pathway (RdDM) (Erdmann and Picard, 2020), in which the concerted action of non-coding RNA polymerases and small RNAs direct the DNA methyltransferase, DRM2, to specific sites of the genome. Newly invading genetic elements are targeted by the ‘non-canonical’ RdDM pathway, and involve the action of Pol II derived small RNAs (Cuerda-Gil and Slotkin, 2016). Once established, the canonical RdDM pathway takes over, whereby small RNAs are generated by the plant-specific polymerase IV (Pol IV). In both cases, a second plant-specific polymerase, Pol V, is an essential downstream component. The recruitment of Pol V to chromatin is stabilized by interaction with the DNA methylation readers SUVH2 and SUVH9 (Johnson et al., 2014). However, methylation alone is not sufficient to explain Pol V occupancy, as Pol V is restricted to the edges of long TEs and euchromatic methylated regions, despite methylation levels being high in the bodies of long TEs and in heterochromatin (Zhong et al., 2012). Therefore, additional mechanisms must exist to regulate Pol V occupancy.

Heterochromatin is targeted for DNA methylation and silencing by CMT methyltransferases (Stroud et al., 2014; Zemach et al., 2013). CMT3 and CMT2 are responsible for CHG and CHH methylation maintenance, respectively, at these regions. The CMTs form an epigenetic feedback loop with the KYP family of H3K9me2 methyltransferases, thereby maintaining high levels of repression-associated non-CG methylation and H3K9me2 at heterochromatic regions of the genome (Du et al., 2014, 2012; Li et al., 2018). As both DNA methylation and H3K9me2 are known to recruit components of the RdDM pathway (Johnson et al., 2014; Law et al., 2013), it is unclear how these regions remain functionally separated from RdDM.

The linker histone H1 associates with nucleosomes to modulate chromatin accessibility and higher order structures (Saha and Dalal, 2021). H1 binds to linker DNA as it exits from the nucleosome dyad and reduces linker DNA flexibility (Bednar et al., 2017). In plants, H1 is preferentially associated with heterochromatin, where it contributes to chromatin compaction (Bourguet et al., 2021; Choi et al., 2020; Rutowicz et al., 2019). Genetic evidence suggest that the chromatin remodeler DDM1 displaces H1 at heterochromatic TEs to facilitate access of CMT methyltransferases (Zemach et al., 2013), however, recent biochemical data suggests that this effect is likely indirect, via DDM1s interaction with the heterochromatin associated histone variant, H2A.W (Osakabe et al., 2021).

Since H1 physically locks genomic DNA to the nucleosome and is preferentially associated with heterochromatin, it is a promising candidate for preventing RdDM encroachment into these regions. Recent evidence supports this hypothesis, with RdDM associated small RNAs becoming more heterochromatically enriched in h1 knockouts (Choi et al., 2021; Papareddy et al., 2020). Here we directly test whether RdDM occupancy is affected by loss of H1, by taking advantage of an endogenous Pol V antibody (which recognizes Pol V’s largest subunit, NRPE1). We found that h1 antagonizes NRPE1 occupancy throughout the genome, particularly at heterochromatic regions, and that this effect is not limited to RdDM as it similarly impacts localization of the more recently discovered methyl-reader and activator complex (Harris et al., 2018).

Results

To understand how H1 affects RdDM occupancy, we performed chromatin immunoprecipitation followed by sequencing (ChIP-seq) with a native antibody against NRPE1 (the largest catalytic subunit of Pol V, previously validated (Liu et al., 2018)). Two biological replicates of ChIP-seq in wild type control Col-0 and in h1.1-1/h1.2-1 double mutant (hereafter referred to as WT and h1, respectively) were performed. As previously reported (Liu et al., 2018; Zhong et al., 2012), NRPE1 was enriched over short TEs and at the edges of long TEs in WT (Fig. 1A-B). In the h1 mutant, however, NRPE1 enrichment was markedly increased, with NRPE1 invading the heterochromatic bodies of long TEs (Fig. 1A-B). The negative correlation between NRPE1 enrichment and TE length observed in WT was reverted in the h1 mutant (Fig. 1C). This suggests that mutation of h1 facilitates the invasion of NRPE1 to more heterochromatic regions, in which long TEs tend to reside. Consistent with this, NRPE1 increased more over heterochromatin associated CMT2 dependent hypo-CHH DMRs than at the RdDM associated DRM2 dependent CHH sites (Fig. 1D). The preferential enrichment of NRPE1 in h1 was more pronounced at TEs that overlapped with heterochromatin associated mark, H3K9me2 (Fig. 1E). Comparing NRPE1 occupancy over TEs that had previously been classified by as either ‘euchromatic’ or ‘heterochromatic’ (based on a broad range of features and small RNA expression dynamics during embryonic development (Papareddy et al., 2020)), again, we found a striking increase in NRPE1 at heterochromatic over euchromatic TEs in h1 (Fig. 1F). From a chromosomal viewpoint, NRPE1 was preferentially enriched over pericentromeric regions in h1, and showed corresponding depletion from the more euchromatic chromosomal arms (Fig. 1G, Supplementary Fig. 1). Importantly, we found no evidence for increased NRPE1 expression at the mRNA or protein level in the h1 mutant (Supplementary Fig. 2). Together, the results indicate that H1 broadly restricts NRPE1 access to chromatin, and that this effect is more pronounced in heterochromatic regions of the genome. As H1 itself is preferentially localized to heterochromatin (see Fig. 1B,D,F) (Bourguet et al., 2021; Choi et al., 2020), this suggests that H1 directly antagonises RdDM occupancy.

Figure 1.

One mechanism by which H1 could restrict RdDM access is by promoting chromatin compaction. Consistent with this idea and previous reports (Bourguet et al., 2021; Choi et al., 2020; He et al., 2019; Rutowicz et al., 2019), we observed increased accessibility in heterochromatic regions in h1, which mirrored the preferential recruitment of NRPE1 to heterochromatin (Supplementary Fig. 3AB). Loss of H1 mediated compaction in heterochromatin may be the indirect result of loss of heterochromatin associated marks such as H3K9me2. To examine this possibility, we performed H3K9me2 ChIP-seq, and although we observed some minor reductions in H3K9me2 in h1, the losses were not preferentially in heterochromatin where we see the majority of NRPE1 gain (Supplementary Fig. 3CD). This suggests that H1 primarily affects the accessibility, rather than the nature of chromatin at these heterochromatic regions.

To determine whether the antagonistic relationship between H1 and NRPE1 is reciprocal, we performed ChIP-seq using an endogenous antibody for H1, with two biological replicates in the WT and nrpe1 mutant backgrounds. H1 levels were broadly unchanged in nrpe1 at TE regions (Fig. 2A). Similarly, over NRPE1 defined peaks (where NRPE1 occupancy is strongest in WT) we observed no change in H1 occupancy in nrpe1 (Fig. 2B). The results indicate that H1 does not invade RdDM regions in the nrpe1 mutant background.

Figure 2.

Next, we asked whether ectopic recruitment of RdDM may be sufficient to evict H1 from the genome. For this, we took advantage of the ZF-DMS3 transgenic lines that we previously characterized, in which zinc finger (ZF) fused DMS3 can recruit endogenous NRPE1 to >10,000 ectopic sites throughout the genome, corresponding to ZF ‘off-target’ binding sites (Gallego-bartolome et al., 2019). With these lines, we performed NRPE1 and H1 ChIP-seq in parallel. Compared to the non-transgenic control, we confirmed that NRPE1 was highly enriched at these ‘off-target’ regions (Fig 2C). However, with the same material, we observed no depletion of H1 at the ectopically bound NRPE1 sites (Fig. 2C). To exclude the possibility that NRPE1 alone may be insufficient to evict H1, but that the fully functioning RdDM pathway may be required (which constitutes at least 32 described proteins acting in concert (Matzke and Mosher, 2014)), we focused on the subset of ZF-DMS3 off-target sites that experience hyper-CHH methylation (Gallego-bartolome et al., 2019). The gain of ectopic methylation at these regions indicates that the full RdDM pathway is recruited and functions to deposit de novo methylation. However, even at these regions where H1 is mildly enriched in WT, we saw no evidence for H1 depletion in the ZF-DMS3 transgenic lines (Fig. 2D). Together, these results indicate that H1 hierarchically and non-reciprocally restricts RdDM from overaccumulation in heterochromatin.

Given our observation that NRPE1 is redistributed in h1, we wondered whether this could help to explain some of the unusual methylation defects observed in the h1 mutant. For instance, h1 mutants experience contrasting loss and gain of DNA methylation at euchromatic and heterochromatic TEs, respectively (Bourguet et al., 2021; Zemach et al., 2013). Previous explanations to account for this suggested that the h1 afforded increased accessibility of euchromatic TEs may promote increased access by activation-associated chromatin modifying enzymes, which in turn antagonise DNA methylation (Zemach et al., 2013). However, based on our NRPE1 ChIP-seq observations, we hypothesized a simpler scenario, whereby the euchromatic loss and heterochromatic gain of methylation is caused by the relative shift of RdDM from euchromatin to heterochromatin. To directly test this possibility, we generated two independent NRPE1 CRISPR knockouts in the h1 mutant background (Supplementary Fig. 4AB), and performed whole genome bisulphite sequencing (WGBS) analysis, alongside nrpe1, h1 and WT controls. We confirmed nrpe1 loss of function in the independent CRISPR KO lines by plotting CHH methylation over NRPE1 associated high confidence DMRs (Zhang et al., 2018), finding that CHH methylation is entirely abolished over these regions in both nrpe1/h1 lines (Supplementary Fig. 4C). By analysis of chromosome level DNA methylation patterns in the different mutant backgrounds, we found that the striking pericentromeric increase in CG and CHG methylation observed in h1, is almost entirely independent of RdDM, as it was not rescued in nrpe1/h1 (Fig. 3A, Supplementary Fig. 5). In contrast, CHH methylation was mostly unaltered between the genotypes at the chromosome scale, with evidence for a slight reduction of CHH in the h1/nrpe1 double mutant as compared to h1 alone (Fig. 3A, Supplementary Fig. 5).

Figure 3.

Next, we compared methylation levels in euchromatic versus heterochromatic TEs. Heterochromatic TEs experienced major increases in CG and CHG methylation in h1 as previously reported (Zemach et al., 2013), and these were independent of NRPE1 function, as the same increases were observed in nrpe1/h1 (Fig. 3B). CHH levels were largely unchanged at heterochromatic TEs in all mutant backgrounds (Fig. 3B). At euchromatic TEs, CG and CHG methylation levels were largely unchanged in h1, while CHH levels were markedly reduced as previously reported (Fig. 3C). CHH methylation levels were further reduced in nrpe1, consistent with the well-established role of RdDM in maintaining CHH at euchromatic TEs (Stroud et al., 2013), with the additional mutation of h1 in nrpe1/h1 having a minimal effect (Fig. 3C). However, we noticed that the loss of CHG methylation in nrpe1 was largely rescued in the nrpe1/h1 double mutant (Fig. 3C). This striking effect was specific to euchromatic TEs. NRPE1 is therefore required for CHG methylation maintenance at euchromatic TEs and in the absence of H1, loss of RdDM can be functionally compensated for by other methylation pathways. This compensation is likely due to the action of CMT3 (Stroud et al., 2014; Zemach et al., 2013). The effect is reminiscent of h1’s amelioration of methylation loss in the chromatin remodeler ddm1 mutant (Zemach et al., 2013), and further supports a role for H1 in demarcating the boundary between heterochromatic and euchromatic methylation pathways.

Calling differentially methylated regions (DMRs) in the mutant backgrounds, we found that h1 hyper CHG DMRs were by far the most numerous as compared to any other context (Fig. 4A). Consistent with the average methylation chromosomal and TE plots, these hyper CHG DMRs were highly enriched over pericentromeric heterochromatin (Fig. 4B). To gain insight into the functional context of methylation in these regions, we performed overlap analysis between h1 hyper CHG DMRs and DMRs from 96 other Arabidopsis gene silencing mutants that were published previously (Stroud et al., 2013). Building similarity matrices based on pairwise overlapping scores (Co-Occurrence Statistics), we observed separation of clustering of mutants belonging to different functional pathways (Fig. 4C, Supplementary Fig. 6), consistent with previous studies (Stroud et al., 2013; Zhang et al., 2018). The h1 hyper CHG DMRs clustered primarily with components of the CG maintenance pathway, including met1 and vim1/2/3 (Fig. 4D). Interestingly, the h1 hyper CHG DMRs also clustered and correlated strongly with ddm1 hyper CHH DMRs. As heterochromatic transposons are known to transition from a state of quiescence to being actively targeted by RdDM in the ddm1 mutant (Panda et al., 2016), this again suggests that H1 controls regions of the genome that are susceptible to spurious RdDM targeting.

Figure 4.

While the majority of methylation changes observed in h1 (hyper CG and CHG) were entirely independent of NRPE1, in the chromosomal plots we noted a subtle depletion of CHH methylation in the h1/nrpe1 double mutant as compared to h1 alone (Fig. 3A), suggesting that redistribution of NRPE1 may have functional consequences on methylation patterning. To explicitly investigate these ectopically bound loci, we compared NRPE1 enrichment in WT versus h1 and identified 15,075 peaks that have significantly higher NRPE1 signal in the h1 mutant. At these regions, CG and CHG methylation changes broadly mirrored that of heterochromatic TEs, again supporting the notion that NRPE1 primarily redistributes to heterochromatic regions of the genome in h1 (Fig. 5A, Supplementary Fig. 7A).

Figure 5.

However, we also observed a significant depletion of CHH methylation at these regions in the nrpe1/h1 double mutant, as compared to either the h1 and nrpe1 alone, which was not observed when looking at heterochromatic TEs as a whole (compare Fig. 5A to Fig. 3B). This indicates that NRPE1 facilitates active deposition of CHH methylation at these newly bound locations. Reciprocally, we asked how methylation changes at regions of the genome that lose NRPE1 in h1 (1,859 peaks), and we found that methylation levels are significantly depleted in all three contexts (CG, CHG and CHH) in the h1, nrpe1 and nrpe1/h1 double mutant backgrounds (Fig. 5B, Supplementary Fig. 7B and 8). The mirrored loss and gain of methylation with changes in NRPE1 occupancy indicates that NRPE1 redistribution in h1 directly impacts methylation patterning in cis.

The results thus far support a model whereby H1 prevents RdDM encroachment into heterochromatin. As H1 likely restricts access by reducing chromatin accessibility at these regions, we reasoned that this effect would not be RdDM specific and could affect other machinery in a similar manner. Recently we identified a transcriptional activator complex, consisting of SUVH1, SUVH3, DNAJ1 and DNAJ2 (Harris et al., 2018). The SUVH1/3 methyl binding proteins directly bind CHH methylation in vitro, yet in vivo the complex shows strong preferential recruitment to RdDM sites, but not to CMT2 deposited CHH heterochromatic sites (Harris et al., 2018). We therefore wondered whether complex’s occupancy might also be affected by H1. To test this, we performed ChIP-seq on SUVH1-3xFLAG in an h1 mutant background (Fig. 6A). As with NRPE1, SUVH1 increased in the h1 background significantly more in heterochromatin, with preferential enrichment over long TEs, cmt2 dependent hypo CHH DMRs, and heterochromatic TEs (Fig. 6B). Therefore, the effect of H1 is not RdDM specific, and acts to demarcate the boundary between euchromatic and heterochromatic methylation pathways, likely through promoting nucleosome compaction and restricting heterochromatin access.

Figure 6.

Discussion

Here we show that loss of H1 results in redistribution of RdDM from euchromatic to heterochromatic regions, consistent with recent findings from small RNA data (Choi et al., 2021; Papareddy et al., 2020). This antagonism, however, appears non-reciprocal, as H1 levels were not affected by RdDM loss, or artificial gain at ectopic regions of the genome. H1-mediated restriction could be in part due to its structural capacity to reduce linker DNA flexibility (Bednar et al., 2017), thereby reducing polymerase access to unwind and interact with genomic DNA. Pol V requires the DDR complex for efficient association with chromatin (Wongpalee et al., 2019), so the presence of H1 could present an additional barrier to this complex. Given that H1 has been shown to compact chromatin in part through the phase separating property of its highly charged intrinsically disordered domain (Gibson et al., 2019), another possibility is that high H1 occupancy creates a biophysical environment that excludes Pol V and related complexes. Our data together with the observed heterochromatic encroachment of Pol IV dependent small RNAs in h1 (Choi et al., 2021; Papareddy et al., 2020), along with H1’s recently described role in preventing spurious antisense transcription at protein coding genes (Choi et al., 2020), support a polymerase indiscriminate model for H1’s antagonism.

Small RNAs are relatively abundant in heterochromatin, despite Pol lV and RdDM being dispensable for methylation in these regions (Stroud et al., 2014). This suggests that Pol IV occupies heterochromatin, but that RdDM is either impotent for methylation deposition or acts redundantly with the CMT/KYP pathway. Here we detect NRPE1-dependent methylation at heterochromatic loci in h1, and combined with the finding that Pol V does not appreciably enter these regions in WT, suggests that H1 prevents RdDM from fully mobilizing in these regions by blocking Pol V. Recent evidence suggests that Pol V can be directly recruited by small RNAs to establish methylation (Sigman et al., 2021), and therefore mechanisms that preclude Pol V recruitment to heterochromatin would be particularly crucial in preventing redundant RdDM activity at regions where methylation has already been established.

In Arabidopsis thaliana, RdDM targets evolutionarily young transposable elements (Zhong et al., 2012). H1 therefore helps to improve the efficiency of RdDM by focusing it upon TEs with the greatest potential to mobilize. As transposons age, they eventually become targeted by CMTs. How this transition occurs is unclear, but our results indicate a potential role for progressive accumulation of H1 in this process. Interestingly in tomato, CMTs rather than RdDM target the evolutionarily younger elements (Wang and Baulcombe, 2020), therefore it will be important to determine whether H1 plays a similar role in maintaining the barrier between methylation pathways in this context. Recent work has shown that H1 co-operates with H2A.W to promote compaction at heterochromatic regions (Bourguet et al., 2021). However, H2A.W was found to antagonize rather than promote H1 deposition at heterochromatin. Therefore, the features that drive H1 to preferentially incorporate into heterochromatin are yet to be discovered.

Methods

Plant materials and growth conditions

Arabidopsis thaliana plants in this study were Col-0 ecotype and were grown under 16h light: 8h dark condition (on soil), or under constant light (on plates). The following plant materials were used: wild-type (WT, non-transgenic), h1 (double homozygous mutant consisting of, h1.1-1 [SALK_128430] and h1.2–1 [GABI_406H11], originally described by (Zemach et al., 2013)), nrpe1 (nrpe1-11 [SALK_029919]), DMS3-ZF (this is a homozygous transgenic line containing the pEG302-DMS3-3xFLAG-ZF construct in Col-0 WT, originally described by (Gallego-bartolome et al., 2019)), SUVH1-3xFLAG in suvh1 (this is a complementing homozygous transgenic line containing the pEG302-gSUVH1-3xFLAG construct in suvh1-1, originally described by (Harris et al., 2018)), SUVH1-3xFLAG in h1 (this line was generated for this work, by transforming the same pEG302-gSUVH1-3xFLAG construct into the h1.1-1 [SALK_128430] and h1.2–1 [GABI_406H11] double mutant background), and the h1/nrpe1 CRISPR lines (described below).

Generation of CRISPR lines

Two independent NRPE1 knockout lines were generated in the background of h1 (h1.1-1 h1.2-1 double, see above), designated as h1/nrpe1 cKO_49 and h1/nrpe1 cKO_63. The lines were generated using the previously published pYAO::hSpCas9 system (Yan et al., 2015) as described by (Ichino et al., 2021). Briefly, the two guides (ATTCTTGACGGAGAGATTGT, and TCTGGCACTGACAAACAGTT) targeting hSpCas9 to NRPE1 (AT2G40030) exons were sequentially cloned in the SpeI linearized pYAO::hSpCas9 construct by In-Fusion (Takara, cat #639650). The final vector was electroporated into AGL0 agrobacteria and transformed in h1 mutant plants by agrobacterium-mediated floral dipping. T1 plants were selected on ½ MS agar plates with hygromycin B and were PCR genotyped to identify plants with deletions spanning the guide region. Selected lines were taken to T2 and genotyped to identify those that had segregated out the pYAO::hSpCas9 transgene. Cas9 negative plants that were homozygous for an NRPE1 deletion (h1/nrpe1 cKO_49) or an indel-induced frameshift to premature stop codon (h1/nrpe1 cKO_63) were confirmed in the T3 generation (see Supplementary Fig. 4).

Chromatin immunoprecipitation (ChIP) and western blot

All the material for ChIP used in this paper were from pooled 10-12 days old seedlings grown on Murashinge and Skoog agar plates (1/2X MS, 1.5% agar, pH5.7). The ChIP protocol used has been previously described (Ichino et al., 2021; Villar and Kohler, 2018), with minor modifications. Briefly, 2-4g of seedling tissue was used for each sample. Samples were ground in liquid nitrogen and crosslinked for 10 or 12 minutes at room temperature in Nuclei Isolation buffer (50mM Hepes pH8, 1M sucrose, 5mM KCl, 5mM MgCl2, 0.6% triton X-100, 0.4mM PMSF, 5mM benzamidine hydrochloride, cOmplete EDTA-free Protease Inhibitor Cocktail (Roche)) containing 1% formaldehyde. The crosslinking reactions were stopped with 125mM glycine by incubation at room temperature for 10 minutes. Crosslinked nuclei were filtered through one layer of miracloth, washed with Extraction Buffer 2, centrifuged through a layer of Extraction Buffer 3 and lysed with Nuclei lysis buffer (buffer compositions are described in the published protocol (Villar and Kohler, 2018)). Chromatin was sheared using a Bioruptor Plus (Diagenode) (20 cycles of 30 seconds on and 30 seconds off) and immunoprecipitated overnight at 4°C with antibody. anti-NRPE1 (endogenous antibody, described in (Liu et al., 2018)), anti-H1 (Agrisera, AS11 1801), anti-H3K9me2 (abcam, 1220), anti-H3 (abcam, 1791), anti-FLAG (for both the SUVH1-3xFLAG and the DMS3-ZF ChIPs, we used Anti-FLAG M2 (Sigma, F1804)). 25μl each of Protein A and Protein G magnetic Dynabeads (Invitrogen) were added to each sample and incubated for 2 more hours at 4°C. The immunoprecipitated chromatin was washed with Low Salt (2X), High Salt, LiCl and TE buffers, for 5 minutes each at 4°C (buffer compositions are described in the published protocol (Villar and Kohler, 2018)) and eluted twice with 250μl of elution buffer (1% SDS and 0.1M NaHCO3) in a thermomixer at 65°C and 1000 rpm, 20 minutes for each elution. Reverse crosslinking was performed overnight at 65°C in 0.2M NaCl, and proteins were degraded by proteinase K treatment at 45°C for 5h. The DNA fragments were purified using phenol:chloroform and ethanol precipitated overnight at -20°C. Libraries were prepared using the Ovation Ultralow System V2 kit (NuGEN, 0344NB-A01) following the manufacturer’s instructions, with 15 cycles of PCR. Libraries were sequenced on a HiSeq 4000 or NovaSeq 6000 instrument (Illumina) using either single- or paired-end 50bp reads.

For the western blot shown in Supplementary Fig. 2B, 100ul of sheared chromatin was used as the input material for SDS-PAGE (as previously described (Harris et al., 2016)), using the anti-NRPE1 and anti-H1 antibodies described above.

ChIP-seq data analysis

ChIP-seq data were aligned to the TAIR10 reference genome with Bowtie2 (v2.1.0) (Langmead and Salzberg, 2012) allowing only uniquely mapping reads with zero mismatch. Duplicated reads were removed by Samtools. ChIP-seq peaks were called by MACS2 (v2.1.1) and annotated with ChIPseeker (Yu et al., 2015). Differential peaks were called by bdgdiff function in MACS2 (Zhang et al., 2008). ChIP-seq data metaplots were plotted by deeptools (v2.5.1) (Ramírez et al., 2016).

Whole-genome bisulfite sequencing (BS-seq) library preparation

Whole genome bisulfite sequencing libraries were generated as previously described (Harris et al., 2016), with minor modifications. The DNeasy Plant Mini Kit (Qiagen #69104) was used to isolate genomic DNA. 150μg of genomic DNA was used as input, a Covais S2 instrument was used for shearing (2 minutes), the Epitect Kit (QIAGEN #59104) was used for bisulfite conversion, and the Ultralow Methyl Kit (NuGEN) was used for library preparation and paired-end libraries were sequenced on a NovaSeq 6000 instrument.

Whole-genome bisulfite sequencing (BS-seq) data analysis

Previously published whole-genome bisulfite sequencing data for mutants and wild type were reanalyzed from previous paper (Stroud et al., 2013). Briefly, Trim_galore (http://www.bioinformatics.babraham.ac.uk/projects/trim_galore/) was used to trim adapters. BS-seq reads were aligned to TAIR10 reference genome by BSMAP (v2.90) and allowed 2 mismatches and 1 best hit (-v 2 -w 1) (Xi and Li, 2009). Reads with three or more consecutive CHH sites were considered as unconverted reads and filtered. DNA methylation levels were defined as #C/ (#C + #T). DMR overlapping analysis were conducted by mergePeaks (-d 100) of Homer (Heinz et al., 2010) with WGBS data published previously (Stroud et al., 2013). The estimated conversion rates for all WGBS libraries are provided in Supplementary Table 1.

Data and code availability

Data supporting the findings of this work are available within the paper and its Supplementary Information files. All high-throughput sequencing data generated in this study are accessible at NCBI’s Gene Expression Omnibus (GEO) via GEO Series accession number GSE225480 (https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE225480). The customized codes used in this study are available upon reasonable request. The following secure token has been created to allow review of record GSE225480 while it remains in private status: XXX

Acknowledgements

We thank Suhua Feng and Mahnaz Akhavan for support with high-throughput sequencing at the UCLA Broad Stem Cell Research Center BioSequencing Core Facility. C.J.H is supported by a Royal Society University Research Fellowship (URF\R1\201016). This work was supported by NIH R35 GM130272 to S.E.J.. S.E.J is an Investigator of the Howard Hughes Medical Institute.

Author contributions

C.J.H. and S.E.J. designed the research. C.J.H., Z.Z. and S.E.J. wrote the manuscript. L.I. helped to isolate the h1/nrpe1 CRISPR knockout lines. S.F. generated the WGBS libraries from genomic DNA. C.J.H. performed all other experiments. Z.Z. and C.J.H. performed the bioinformatic analyses.