Multi-Omics Data Integration Reveals Link Between Epigenetic Modifications and Gene Expression in Sugar Beet (Beta Vulgaris Subsp. Vulgaris) in Response to Cold

DOI: https://doi.org/10.21203/rs.3.rs-949141/v1

Abstract

Background

DNA methylation is thought to influence the expression of genes, especially in response to changing environmental conditions and developmental changes. Sugar beet (Beta vulgaris ssp. vulgaris), and other biennial or perennial plants are inevitably exposed to fluctuating temperatures throughout their lifecycle and might even require such stimulus to acquire floral competence. Therefore, plants such as beets, need to fine-tune their epigenetic makeup to ensure phenotypic plasticity towards changing environmental conditions while at the same time steering essential developmental processes. Different crop species may show opposing reactions towards the same abiotic stress, or, vice versa, identical species may respond differently depending on the specific kind of stress.

Results

In this study, we investigated common effects of cold treatment on genome-wide DNA methylation and gene expression of two Beta vulgaris accessions via multi-omics data analysis. Cold exposure resulted in a pronounced reduction of DNA methylation levels, which particularly affected methylation in CHH context (and to a lesser extent CHG) and was accompanied by transcriptional downregulation of the chromomethyltransferase CMT2 and strong upregulation of several genes mediating active DNA demethylation.

Conclusion

Integration of methylomic and transcriptomic data revealed that, rather than methylation having directly influenced expression, epigenetic modifications correlated with changes in expression of known players involved in DNA (de)methylation. In particular, cold triggered upregulation of genes putatively contributing to DNA demethylation via the ROS1 pathway. Our observations suggest that these transcriptional responses precede the cold-induced global DNA-hypomethylation in non-CpG, preparing beets for additional transcriptional alterations necessary for adapting to upcoming environmental changes.

Background

Sugar beet is the main source of sugar production in Europe. Its fleshy taproot provides around 15% of the world’s annual production of sugar, solely competing with sugar cane in industrial sucrose extraction worldwide [1]. Conventionally, sugar beets are sown in spring and harvested before winter, still in the vegetative stage. In contrast, sowing sugar beet in autumn and harvesting in the following year could increase - due to the development in the winter and the acceleration of growth in the spring - the sugar content by up to 26% [2]. Therefore, one of the major goals in sugar beet breeding is the production of winter beets with the challenge to ensure the viability over the winter season, since particularly seasonal environmental stresses such as frost limit possible periods for their cultivation and thereby also restrict potential enhancement of yield. Biennial sugar beets grow vegetatively in the first season and after extended exposure to cold during the winter season, they switch to generative reproduction and acquire floral competence, as a result of vernalization and the subsequent shift towards long-day conditions in spring [3, 4]. The shoot outgrowth or stem elongation, so-called bolting, and the following flower development drastically reduce the sugar yield and the size of the beet.

Plants, as sessile organisms, are known to use several strategies to cope with changing environmental conditions through significant alterations in gene expression, or epigenetic processes to ensure phenotypic plasticity. One epigenetic key mechanism is DNA methylation, which forms 5-methylcytosine nucleotides (5-mC) by covalently linking methyl-group(s) to specific cytosines at carbon position 5. Cytosine methylation in plants is facilitated by different enzymes acting on cytosines dependent on their 3’ nucleotide environment, with distinct methyltransferases being able to modify cytosines in CpG, CHG, or CHH context, where H represents A, T, or C (Figure 1b). Sites in CpG and CHG context are considered “symmetric”, as the opposite strand inevitably carries an identical motif including a potential methylation target, which is represented by the guanine-pairing cytosine of the opposite strand [5]. Cytosines in CHH context, or methylation thereof, accordingly is “asymmetric” by nature. Apart from acting specifically towards different cytosine environments, the particular set of enzymes modulating DNA methylation is also dependent on whether methylation is being established de novo, or whether an existing methylation mark is being maintained [6]. De novo methylation is established through the RNA-directed DNA methylation (RdDM) pathway and affects cytosines in all contexts (CG, CHG, CHH). However, newly established modifications at symmetric sites would not be carried over to the daughter strand during DNA replication without an additional maintenance system. A set of distinct methyltransferases, counteract the passive loss of methylation in a context-specific manner, whereby METHYLTRANSFERASE 1 (MET1) maintains CG methylation [7, 8], whereas CHROMOMETHYLTRANSFERASE 3 (CMT3) or CMT2 maintain methylation of cytosines in CHG context [912]. CMT3 acts in presence of specific histone modifications (H3K9me2), placed by SUPPRESSOR OF VARIEGATION 3-9 HOMOLOGUE PROTEIN 4 (SUVH4), SUVH5, and SUVH6, which in turn are recruited to (CHG-) methylated DNA, thereby forming a reinforcing loop between CHG methylation and H3K9 methylation [13]. CHH methylation is mediated chromatin-specifically - either by DRM2-dependent RdDM (euchromatin), or - in H1-containing heterochromatin where RdDM is blocked - by CMT2. Cytosine methylation can also be actively removed via excision of the methylated nucleotide, followed by repair of the cleavage site (base excision repair pathway, BER). Excision of a methylated nucleotide can be facilitated by bifunctional 5-mC DNA glycosylases, such as REPRESSOR OF SILENCING 1 (ROS1), TRANSCRIPTIONAL ACTIVATOR DEMETER (DME), or DEMETER-LIKE PROTEINs (DML2 and DML3). Following further modification of the cleavage site, the gap is (filled and) repaired via DNA polymerase and ligase enzymes, e.g. LIG1 in Arabidopsis [1419].

RNA sequencing has already been used in numerous studies to investigate altered gene expression under environmental stress in plants [20]. In addition, epigenetics relating to adaptive changes in plants gain more and more importance. Especially, the application of whole-genome bisulfite sequencing (WGBS) brings expansion of knowledge about the epigenetic mechanisms controlling the development and adaptation of many plants. The treatment with bisulfite converts unmethylated cytosines into uracil while methylated cytosines stay unaffected, which allows the detection of each cytosine's methylation status. Genome-wide cytosine methylation analysis in A. thaliana already revealed methylation patterns in genes and repetitive areas [21, 22]. In crops, this technique has also been used to identify epigenetic mechanisms that have a significant impact on stress responses or developmental processes [2327]. Many other publications also have demonstrated the importance of DNA methylation in different plants concerning responses to abiotic stress, including salt, drought, cold, heat, and heavy metal exposure [2837]. Hébrard et al. [38] already reported genotype-dependent DNA methylation and their association to the expression of genes participating in bolting tolerance in sugar beet by comparing bolting-sensitive and bolting-tolerant genotypes. Further, Zakrzewski et al. [39, 40] published insights into the epigenome of sugar beet leaves and callus, focusing on methylation alterations at satellite DNAs and transposable elements. To improve our current understanding of cold-adaptation mechanisms in sugar beets, we investigated DNA methylomes of two B. vulgaris genotypes under normal conditions and after exposure to cold, focusing on conserved epigenetic alterations and their association to genome-wide gene expression. Our work provides insights into genome-wide, genotype-independent transcriptomic and epigenomic responses of Beta vulgaris to cold exposure and links the expression of putative methylation- or demethylation-related genes to alterations in DNA methylation.

Results

Experimental setup and workflow

As depicted in the overview of our workflow (Figure 1), we used plants from two different sugar beet accessions that were either grown under control conditions for 15 weeks (from now on, corresponding samples are grouped under the label CONTROL), or for 12 weeks followed by a three-week incremental cold treatment ending at a minimum temperature of 0°C (referred to as COLD). Leaf material was collected from six samples (2 genotypes x 3 biological replicates each) per condition (see Materials & Methods). From this material, DNA and RNA was extracted to subsequently analyze both the methylomic and transcriptomic cold responses (Figure 1a).

Genome-scale methylomic and transcriptomic profiling

The twelve DNA samples were subjected to directional WGBS (PE150+150) using the Illumina NovaSeq 6000 platform (Novogene, Beijing, China). After adapter and quality trimming (Q20), an average of 90 million paired-end reads (150+150) (~27 Gb) per sample were obtained, representing 36-fold coverage of the sugar beet reference genome [41]. Mapping and deduplication revealed average proportions of uniquely mapped and multi-mapped reads of 60% and 12%, respectively.

For RNA-seq, RNA extracted from all six CONTROL and all six COLD samples was used to construct poly-A selected libraries and sequenced using the Illumina Inc. HiSeq 2000 system (GATC GmbH, Konstanz, Germany). After adapter and quality trimming (Q20), an average of 88 million paired-end (151+151) reads (~13 Gb) per sample were obtained. The average proportions of uniquely mapped and multi-mapped reads were 91.06% and 5.81%, respectively (see Figure 1a and Methods).

DNA methylation patterns in sugar beet under control and cold conditions

After methylation extraction, the proportions of methylated cytosines in each context (CpG, CHG, and CHH) were analyzed per chromosome and in the genome as a whole (Figure 2). For beets grown under control conditions, unequal distribution of C methylations (mC) were detected, depending on the local sequence context, being proportionally more common in CpG (66.8 %) than in CHG (41.7 %) or CHH (9.4 %). Cold treatment reduced these percentages for mCpG (64.0 %), mCHG (40.5 %), and mCHH (8.7 %), although there was no statistical significance in any case (two-tailed Welch’s test; Figure 2). A trend towards cold-dependent demethylation was also detected on gene-feature-level. To describe these methylation patterns in the different parts of a gene, we averaged the methylation proportions for each possible context (among all 6 samples in the corresponding group) in the following regions of all annotated genes of the reference genome [41, 42]: promoter, defined as the range from 2.5 kb upstream of the transcription start site (TSS) to the TSS; 5'-untranslated region (UTR); coding sequence (CDS), intron and 3'UTR. As shown in Figure 2, methylation predominantly affected promoter and intron regions, independent of the sequence context. In contrast to CHG or CHH, CpG showed high methylation levels in CDS relative to other gene regions. The lowest methylation rates occurred along UTRs for CpG, and in CDS for CHG and CHH, respectively. Even though only 5’- and 3’- UTRs showed statistically significant methylation differences between CONTROL and COLD (two-tailed Welch’s test, p < .05) in CHH, COLD showed slightly lower average methylation than CONTROL in almost all genomic regions and cytosine contexts analyzed (Figure 2).

Global methylation levels of individual chromosomes, the entire genome (panels on the left), or within gene components (panels on the right) with light and dark blue boxes representing methylation levels of CONTROL and COLD samples, respectively. Asterisks indicate significant differences (p-value ≤ .05; one-sided t-test) between methylation levels of CONTROL and COLD.

Cold-dependent differential methylation of individual cytosines

To explore significant alterations of DNA methylation in response to cold treatment at individual positions, we analyzed methylation for each cytosine (depth 8 in each replicate) of all chromosome-assigned scaffolds constituting the B. vulgaris genome using methylKit [43]. A cytosine site with a q-value (SLIM adjusted p-value) of < .05 was defined as differentially methylated (DMC = differentially methylated cytosine), with positions being labeled “hypermethylated” (e.g. dmCpGhyper) when methylation rates were significantly higher, or “hypomethylated” (e.g. dmCpGhypo) when methylation rates were significantly lower in cold-treated samples compared to the controls (Figure 1c). We identified 5319 DMCs distributed over the nine chromosomes (Figure 3a), with the highest DMC-density on chromosome 6 (19.7 DMCs per Mbp), and the lowest on chromosome 1 (7.5 DMCs per Mbp; Figure 3c). 3840 (72.2%), 1041 (19.6%) and 438 (8.2%) DMCs could be assigned to CpG, CHG, and CHH, respectively, with slightly higher numbers of dmCpGhypo and dmCHHhypo compared to corresponding hypermethylated Cs (Figure 3a). Of all detected DMCs, 2873 (54.0%) were located within genes (incl. promoters), with dmCpG accounting for more than 80% of all gene-associated DMCs, approximately reflecting the proportion of DMCs in CpG context to non-CpG-DMCs.

Overall, most DMCs overlapped introns, CDSs, or promoters - together accounting for more than 90% of all gene-associated DMCs (Figure 3b), except for dmCHHhypo, which showed a preference for promoters and 3’UTRs. Additionally, dmCHGhyper and dmCHHhyper showed particularly strong associations with introns.

In plants, methyltransferases and demethylases favor specific sequences for local (de)methylation [44]. To evaluate possible sequence preferences, we scanned the sequence composition within 25bp-flanking-regions (25bp up- and downstream of each DMC, i.e. 51bp in total) of all cold-dependent DMCs, separately for hyper- and hypomethylation in CpG, CHG, or CHH context, respectively (see Additional file 1). In flanking regions of dmCHG and dmCpG, we observed mainly A/T rich sequences, with a slightly higher frequency of Gs 5’ proximal to dmCHG (position -4 relative to DMC position). In contrast to dmCpG and dmCHG, Gs and Cs occurred more frequently in the neighborhood of dmCHH. Within trinucleotide-environments of hypermethylated DMCs, H positions directly next to differentially methylated cytosines were almost exclusively occupied by T or A. Some of the cytosines which had instead lost methylation in response to cold treatment preceded another (not necessarily differentially methylated) cytosine. These results are comparable to identified sequence preferences for methylation in A. thaliana [21, 22].

Cold-dependent differential methylation of genomic regions (DMRs)

Differentially methylated regions (DMRs) are stretches of DNA including multiple cytosines with collectively altered methylation levels between samples. We used HOME (Histogram of methylation), a machine learning-based tool, to identify DMRs while taking into account the sequencing depth and spatial correlation of cytosines [45]. We decided to filter all detected differentially methylated regions for a minimum methylation difference of 10% (spanning the whole DMR) between CONTROL and COLD. This left a total of 3557 DMRs (Figure 3d), with 2760 DMRs (77.6%) in CpG, 700 (19.7%) in CHG, and 97 (2.7%) in CHH context. 699 DMRs (19.7%) - all of which were assigned to either CpG or CHG - contained at least one DMC (in corresponding C context). We investigated the presence of DMRs within gene features, and – again similar to DMCs - the number of DMRs in CpG context was highest in introns, CDSs, and promoters and lowest in the 5’UTR (Figure 3e). DMRs in mCHG and mCHH showed very similar patterns. Comparable to DMCs, most of the DMRs were located on chromosome 6, while chromosome 1 showed the fewest DMRs (Figure 3f).

CHG and CHH methylation is associated with low gene expression

Integration of gene expression and DNA methylation requires consideration of sequence contexts of cytosines while discriminating between different types of genomic features affected by methylation. We used MethGet to assess associations between expression and methylation for each replicate, individually [46]. Gene expression data was normalized for transcript length to allow for the comparison of gene expression within samples. Each gene was assigned to one of six bins, based on its relative expression level. Subsequently, the methylation level along all genes within the same bin, as well as of their flanking regions (spanning about half the length of a given gene 5’ of the TSS, and half the gene-length 3’ of the TTS [transcription termination site]) was evaluated. As shown in Figure 4, the highest methylation levels were reached in flanking regions, particularly in those upstream the TSS, and this tendency was independent of the extent of expression, sequence context, or treatment. Generally, methylation of cytosines in CpG context showed a marked drop in a narrow region around the TSS, as well the TTS with methylation levels within genes reaching those observed further upstream of the gene. In CHG and CHH context, the methylation level also decreased when approaching TSS or TTS, but remained comparably low within genes. Interestingly, expression levels were not necessarily associated with the total gene methylation level of cytosines in CpG, but rather with the amplitude of the drop in methylation near the TSS, with highly expressed genes showing the steepest decline, i.e. the lowest methylation levels around the TSS and lowly expressed genes showing the lowest deviation from upstream- and within-gene regions and thus the highest methylation level in proximity to their TSS. In addition, expression of the corresponding genes seemed to be negatively correlated with the methylation levels of CHG and CHH along the entire length of the gene. There, high expression was accompanied by low methylation levels throughout, whereas methylation of lowly expressed genes was higher, approximating upstream and downstream methylation levels. In general, methylation profiles were almost identical between CONTROL and COLD. However, cold treatment uniformly decreased DNA methylation levels in CHH context for all gene-expression bins.

Effects of cold on gene expression

Analysis of differential expression [47] contrasting CONTROL vs. COLD yielded in total 2549 DEGs (Additional file 2), most of them localized on chromosome 6, fewest on chromosome 8 (Figure 5). Of them, 1244 were up-regulated, while 1305 were down-regulated in COLD (Wald-test; Bonferroni ≤ 0.05; |log2FC| ≥ 1).

Differential methylation does not globally correlate with gene expression

We identified 209 differentially expressed genes (from a total of 2549 DEGs; ≙ 8.2 %) that (including their promoter regions) were co-localized with at least one DMC. Comparing the proportion of differentially methylated genes among the DEGs (8.2 %) with the fraction of all annotated genes overlapping at least one DMC (6.6 %, corresponding to a total of 1712 DMC-associated genes), association with differential methylation was significantly enriched among our set of DEG [χ²(1, n= 26004) = 4.9091, p < .05]. However, combinations of hyper- or hypomethylation and up- or downregulation of expression did not reflect an apparent tendency (Figure 6a): from the total of 1244 transcriptionally upregulated genes, 66 were overlapped by at least one hypermethylated DMC, while similarly, 62 were overlapped by at least one hypomethylated DMC. Analogously, of the total of 1305 genes that were significantly downregulated by cold, 46 DEGs colocalized with hyper-, and 43 DEGs colocalized with hypomethylated DMCs, respectively. DMCs preferentially occupied CDSs, introns, and 3’UTRs of DEGs, independent of the sequence context of the DMC. For genes showing unaltered expression following cold exposure, this pattern shifted from 3’UTRs towards promoters (Figure 6a).

Similarly, 260 DEGs coincided with at least one of the total of 3557 regions (DMRs) where methylation levels differed by at least 10 % between CONTROL and COLD. The relative distribution of these DMRs within the affected DEGs was similar to those observed for DMCs, i.e. DMRs preferentially occurred within CDSs, introns, and 3’UTRs of DEGs (Figure 6b). Of all genes, 2287 (≙ 8.8 %) showed association with a DMR, and comparing the absolute numbers, DEGs were again enriched with DMRs [χ²(1, n = 26004) = 28.9613, p < .01]. This could lead to the conclusion that the methylation status of our set of DEGs dynamically changed not only at single positions (DMCs) but also in extended areas of methylation alterations over multiple cytosines (DMRs) as a response to cold treatment.

Differential gene expression correlates with changes in DNA methylation

Vice versa, we systematically screened our set of DEGs for known players involved in DNA (de)methylation and mapped the matching DEGs to general pathways known to affect DNA methylation. We found numerous genes whose homologs have been described to contribute to either a) de novo methylation via RdDM, b) chromomethyltransferase-mediated maintenance (or - for CHH - de novo establishment) of DNA methylation or c) active DNA demethylation and whose expression was significantly altered in response to cold (Figure 5 and Figure 7): cold treatment significantly reduced expression levels of subunits of DNA-directed RNA polymerase IV (NRPD1; Bv5_101150_qmjs) and DNA-directed RNA polymerase V (NRPE5A; Bv3_062710_wdck), possibly resulting in reduced production of siRNAs and scaffold RNAs for RdDM. In line with this, chromatin remodeler SNF2 DOMAIN-CONTAINING PROTEIN CLASSY 1 (CLSY1; Bv2_044610_uypz), an interactor of Pol IV, as well as several RNA-dependent RNA polymerases (RDR1, RDR3, RDR6, corresponding to Bv2_040860_udas, Bv8_197280_tuhe, and Bv2_030230_aisi, respectively), which are thought to produce dsRNA from mainly Pol IV-derived templates, were significantly downregulated, as was the expression of DRB4 (Bv5_109930_mrpd), which is thought to assist the processing of dsRNA fragments to siRNA via DCL4 [48, 49]. Furthermore, cold significantly decreased the expression of IDN2 (INVOLVED IN DE NOVO 2; Bv2_024880_hakm) - an interactor of DMR2 required for RdDM - as well as of DNA (cytosine-5)-methyltransferase CMT2 (Bv3_050080_yren), which mediates de novo methylation preferentially at cytosines in CHH context [50, 51]. And finally, cold-dependent transcriptional downregulation was also observed for an S-adenosylmethionine synthase (SAM2; Bv4_079640_yozf), whose product acts as a major methyl-donor in DNA methylation. Together, these transcriptional changes overall indicate a reduction in RdDM derived de novo methylation and decreased CMT2-based CHH methylation maintenance. In addition, we detected several transcriptional changes that might indicate enhancement of active DNA methylation upon cold treatment: the cold-triggered increase of NPX1 (Bv5_110670_yfki) expression, together with enhanced expression of the SWR1-components PIE1 (Bv3_065010_faku) and SWC4 (Bv5_122550_anjp) could promote the recruitment of ROS1 via deposition of H2A.Z. Furthermore, the expression of ROS1 (Bv7_160320_kstp) itself was upregulated upon cold treatment. (Note that a DNA LIGASE 1-LIKE homolog, Bv_001710_ogue, was additionally upregulated more than four-fold. However, for conciseness, we excluded this gene from global analyses due to its location on a sequence component not assigned to any chromosome.) Expression of a homolog of the putative DNA glycosylase At3g47830 | DEMETER-like protein 2 (DML2; Bv5_116310_qfjg) was significantly lowered by cold (Figure 7).

Discussion

Characteristics of cytosine-methylation in sugar beet

In this work, we have analyzed and compared the methylome (WGBS) and transcriptome (RNA-seq) of sugar beet under control conditions or after exposure to cold. Under control conditions, cytosine methylation proportions were 66.8% (mCpG), 41.7% (mCHG), and 9.4% (mCHH) for total sequenced CpG, CHG, and CHH, respectively (Figure 2). This deviates from the results obtained by Niederhuth et al. [52], who reported methylation rates for B. vulgaris of as high as 92.5% (mCpG), 81.2% (mCHG), and 18.7% (mCHH). It might be noteworthy that the authors highlighted that, particularly the percentage of mCHH in their study, was “unusually high” and define B. vulgaris as a “notable outlier”. The lower methylation rates for the three analyzed contexts (mCpG, mCHH and mCHG) in our analysis, in contrast, are in more concordance with the values of other angiosperm species included in Niederhuth et al. [52] and with the methylation proportion for B. vulgaris detected by Zakrzewski et al. [40]. Although data was obtained from leaf material in both studies, factors such as plant age or growth conditions [neither are documented for the comparative analysis of 52], the use of different accessions or the fact that sequencing data was generated and processed using different methods and computational tools, might to some extent account for the differences in observations of methylation levels.

Beta vulgaris alters global methylation level in response to cold

After cold treatment, methylation proportions dropped to 64.0% (mCpG), 40.5% (mCHG) and 8.7% (mCHH) which, in all cases, represented a decrease in total methylation, although none of them was statistically significant (two-tailed Welch’s test). Previous studies have also detected low to medium but consistent variations in methylation proportions due to abiotic stresses. However, the intensity and direction of the change of hypo- or hypermethylation depends on the species, type of treatment and genotype pre-adaptation (Sudan et al., 2018). In upland cotton (Gossypium hirsutum), for example, Lu et al. (2017) found a decrease of 2% in mCpG and mCHG, and an increase of 4% in mCHH, after drought stress. Gayacharan and Joel [53] observed that, after drought treatment, drought sensitive rice genotypes were generally hypermethylated, whereas the drought tolerant ones were hypomethylated. Salinity stress caused DNA demethylation in Setaria italica [25], but increased methylation in Medicago truncatula [54]. Cold stress was found to trigger DNA demethylation in roots of maize seedlings [55]. Also tartary buckwheat (Fagopyrum tataricum) tended to decrease DNA methylation following cold shock (at 0°C for several hours) or repeated short term (< 1 day) exposure to cold. However, mCHH was actually slightly denser in repeatedly cold-exposed plants compared to seedlings cultivated under control conditions [56]. In contrast, cucumber (Cucumis sativus) radicles increased DNA methylation in response to cold treatment, which was paralleled by reduced growth rates [57]. In rubber trees (Hevea brasiliensis), and fruits of sweet oranges (Citrus sinensis), or of tomato (Solanum lycopersicum), cold was shown to alter expression of specific genes involved in cold-response, volatile-, or pigment-synthesis through altered DNA methylation of corresponding promoters [37, 58, 59]. Finally, in some woody perennials like poplar, cold triggers DNA demethylation prior to bud break [60].

Large amplitudes in methylation profiles distinguish highly from lowly expressed genes

mCG, mCHG and mCHH patterning has been linked to gene expression, but methylation levels can widely vary between species. In rice, mCpG is lowest at the flanks of the TSS (approximately 100bp upstream and 500bp downstream of TSS) and the TTS (approximately 500bp upstream and 100bp downstream of TTS), mCHG and mCHH are basically absent from genes (Zemach et al., 2010). Cassava, soybean and A. thaliana have similar methylation patterns, with almost exclusive methylation of cytosines in CpG context, which become very low around the TSSs and TTSs [21, 35, 61]. Further, mCpG at the TSS and TTS tends to inversely correlate with expression levels of the corresponding gene [62, 63]. Independent of the cytosine context, we observed the lowest methylation levels in close proximity to the TSS and TTS regions of genes and much higher levels of methylation in the 2kb up- and downstream region of TSS and TTS, and - at least for mCpG - an intragenic increase in methylation. These profiles match those detected by Zakrzewski et al. [40], who reported very similar methylation patterns based on WGBS data from sugar beet leaf material.

We detected a genome-wide hypomethylation in response to cold treatment. This trend was quite evenly distributed over the whole genome of B. vulgaris with significant alterations especially in CHH context (and to some extent CHG) after cold treatment in almost all gene regions analyzed (Figure 2). Although we observed methylation in CHG and CHH context in all genes to be associated with a lower level of expression, our set of differentially expressed genes showed no obvious pattern of up- or downregulated expression and in- or decreased methylation in any sequence contexts (Figure 4 and Figure 6). This means that while absolute methylation levels in CHG and CHH context do, in fact, inversely correlate with gene expression (Figure 4), changes in methylation (i.e., DMCs or DMRs) are not necessarily associated with a significant change in expression (Figure 6).

On the other hand, we identified many homologs of known players in our set of DEGs, whose functions were already attributed to de novo- or maintenance of DNA methylation or active DNA demethylation in A. thaliana and which might explain the observed changes in methylation on a transcriptional basis.

Differential expression of chromatin modifiers and DNA (de)methylaters correlates with changes in global methylation

Cold limits precursor synthesis for DNA methylation

S-adenosylmethionine acts as co-substrate and methyl-donor for several biochemical reactions, including histone- and DNA methylation (Figure 8a). SAM (S-adenosylmethionine synthetase) catalyzes the final step in the synthesis of the molecule. Altered expression of SAMs not only affects expression of genes related to abiotic stress [6467], but also modifies the methylation status of DNA. In detail, missense mutation of SAM3/MAT4 in Arabidopsis was shown to result in decreased CHG and CHH DNA methylation [68]. While cold decreased transcription of a SAM homolog in sugar beet leaves, its product is involved in such a multitude of processes, that it is difficult to predict hypomethylation based on decreased expression of SAM2 alone.

Cold alters RNA-directed DNA methylation

However, apart from precursor synthesis, cold treatment had significant impact on expression of some enzymes considered as core components of DNA methylation-, as well as demethylation-pathways.

In plants, de novo DNA methylation of cytosine residues in all contexts is established through RdDM (Figure 8b). Canonical RdDM is initiated through the activity of DNA-directed RNA-polymerase IV [6972], which generates short, non-coding, single-stranded transcripts at sites enriched for particular histone modifications (H3K9me2). The histone marks are recognized by SHH1, which - in the presence of CLSY1 [73, 74] - triggers binding of Pol IV [7577]. Prior to acting as actual guide for DNA methylation, Pol IV transcripts are processed by a series of further enzymes: RNA-DEPENDENT RNA POLYMERASE 2 (RDR2) complements Pol IV transcripts into dsRNA [7880], which can be cleaved into siRNAs by DICER-LIKE PROTEINs [81, 82], before being loaded to ARGONAUTEs [mainly AGO4 and 6;, 83, 84]. At the target site, nascent transcripts produced by another core RdDM-polymerase, Pol V, can pair with appropriate AGO-coupled siRNAs [85]. This finally recruits DOMAINS REARRANGED METHYLASE 2 (DRM2), an interactor of AGO4, to methylate cytosines at the target site [8688].

Although we observed significantly altered expression of core RdDM components, we do not believe that the substantial global reduction of methylation levels after cold treatment can be attributed to reduction of de novo methylation by RdDM: As shown in Figure 7 and Figure 8b), cold significantly reduced expression levels of the largest subunit of Pol IV (NRPD1), as well as another polymerase-subunit specific to Pol V (NRPE5A), which together might decrease production of siRNAs and scaffold RNAs for RdDM [8991]. In addition, expression of CLSY1 was significantly decreased in cold-treated sugar beets, as was the expression of three RDRs (RDR1, 3, 6; Figure 8b) of which at least two participate in non-canonical RdDM at stages of alternative siRNA production in Arabidopsis [92]. Furthermore, cold decreased expression of DRB4, which is thought to assist processing of dsRNA fragments to siRNA via DCL4 [48, 49], as well as of IDN2 (INVOLVED IN DE NOVO 2), which binds dsRNA, associates with DRM2 and is required for recruitment of SWI/SNF chromatin remodelling complex components to RdDM target sites [93, 94]. It is thought to normally promote DRM2-mediated DNA methylation at some target loci [95, 96].

However, expression of the DRM2-homolog, DRM3, was significantly upregulated by cold. But while DRM3 is thought to regulate DRM2-mediated DNA methylation and to be required for normal maintenance of non-CG methylation, DRM3 itself lacks a conserved site required for methyltransferase activity and compared to drm2 mutants, Arabidopsis plants lacking functional DRM3, only show moderate losses of DNA methylation [9799]. Additionally, some of the most crucial (canonical) RdDM components, e.g. RDR2, DCLs, AGO4 or AGO6 and, most importantly, the major RdDM DNA methyltransferase, DRM2, remained stably transcribed in cold-treated plants (Figure 8b). And while AGO2 and Pol II – both of which were upregulated in cold-treated samples (Figure 8b) - can be involved in providing siRNAs that eventually guide DNA methylation via DRM2 through non-canonical, partially Pol IV-independent mechanisms of RdDM [92], the protein products of both genes mainly mediate other regulatory mechanisms not necessarily affecting DNA methylation.

Depending on the target and scaffold RNAs involved, RdDM can facilitate methylation of several consecutive cytosines and thus, presence of hypermethylated DMRs could hint towards increased RdDM at these positions. In our data, there was a higher number of individual CpG sites showing significant hypomethylation, than there were hypermethylated CpG positions (Figure 3a). In contrast, among longer stretches of DNA that showed considerable (> 10%) cold-dependent changes in CpG methylation (i.e. DMRs in CpG context), the majority had actually gained CpG-methylation upon cold exposure (i.e. more hyper- than hypo-methylated CpG-DMRs; see Figure 3d). DRM2 is able to methylate cytosines independent of their sequence context [86, 100]. Therefore, DMRs in a particular sequence context which have gained methylation through RdDM (i.e. through DRM2), are expected to show overlap - at least to some extent - with hypermethylated DMRs in other sequence contexts. However, although there were about 100 CpG-DMRs that coincided with CHG-DMRs and showed matching trends in their methylation change (both overlapping DMRs either gained methylation, or both exhibited loss of methylation, but not a combination of both), about half of them were hypomethylated DMRs. The remainder, i.e. overlapping hypermethylated CpG/CHG DMRs might indeed represent regions with locally enhanced RdDM activity.

Overall, while the altered expression of several RdDM- related enzymes might indicate a change in RdDM, our observations do not support the conclusion that loss of RdDM is the major determinant of the global reduction in methylation after cold treatment.

Downregulation of CMT2 suggests hypomethylation in CHH

Generally, when or if mCpG, mCHG and/or mCHH are established, these methylation marks must be actively maintained in order to retain correct DNA methylation patterning (Figure 8c). Otherwise, methylation information gets lost on the daughter-strand during DNA replication (because there is no cytosine in the complementary sequence of CHH, mCHH is, per se established de novo). In contrast to DRM2-mediated methylation, the methyltransferases involved in these pathways act sequence specific: DNA METHYLTRANSFERASE 1 (MET1) complements missing methylation in the complementary strand of hemimethylated mCpG sites [7, 8], whereas non-CpG methylation(-maintenance) relies on the CHROMOMETHYLTRANSFERASES CMT3 (mCHG and to lesser extent mCHH) and CMT2 [mCHH and mCHG;, 101]. Chromomethyltransferases preferentially act on chromatin carrying H3K9me2. Histone methylation, in turn, is established by certain H3 lysine-9-specific methyltransferases of the Suvar3-9 subfamily (SUVH4 and its homologs SUVH5 and SUVH6) which are able to bind methylated DNA and preferentially act on histones in proximity of existing mCHG and mCHH [9, 50, 102105].

Our RNA-seq data indicated contrasting transcriptional regulation of SUVH6 - which was upregulated - and CMT2 - of which one homolog (Bv3_050080_yren) was drastically downregulated in leaves of cold treated sugar beets (Figure 8c). Under the assumption that sequence- or modification-specificities are conserved in sugar beet, an increase of SUVH6 should promote H3K9me2 at sites with pre-methylated CHG and CHH [104, 106]. H3K9me2, in turn, recruits CMT3 resulting primarily in mCHG maintenance. In the absence of CMT2, (possibly increased) H3K9me2 (placed by SUVH6) thus should favor mCHG through CMT3, whereas mCHH is expected to decrease. This is in line with the high ratio of dmCHHhypo:dmCHHhyper and accordingly fits with a slightly higher number of hyper- compared to hypomethylated DMCs in CHG context. However, on a global scale, hypomethylation was observed in all cytosine contexts - including CpG, which indicates that reduction of methylation is not, or at least not exclusively due to reduced CMT-activity. Moreover, as shown in Figure 7f), at least one of the two genotypes (GT1) we included in our analysis expressed relatively high levels of another close CMT2 homolog, possibly substituting for the downregulation of Bv3_050080_yren in this genotype.

A considerable fraction of mCHH has been proposed to be dependent on CMT-, instead of RdDM-based DNA methylation. This mainly affects heterochromatin, where RdDM is blocked, whereas methylation via CMT2 can persist [11]. Accordingly, in case of severely reduced CMT2-activity, there should be an accumulation of dmCHHhypo particularly in heterochromatic regions, which tend to be rich in TEs, but are usually characterized by a low density of protein coding genes. However, cold-dependent dmCHH was rather evenly distributed over chromosomes.

Considering the relatively high expression of another, possibly redundant CMT2-homolog in one of the genotypes in our study, the even distribution of differential CHH methylation across chromosomes, and rather ubiquitous hypomethylation (not only in CHH), we see no clear evidence for a major contribution of CMT2 to the observed effects of cold on methylation.

A ROS1-homolog might drive active DNA methylation in response to cold

In addition to passive loss of DNA methylation, plants utilize a complex enzyme machinery to actively erase DNA methylation at specific positions (Figure 8d). Removal of cytosine-methylation in plants comprises excision of the entire methylated cytosine - as opposed to mere cleavage of the methyl-group - followed by repair of the cleavage site via insertion and ligation of an unmethylated cytosine [1417]. Removal of the methylated cytosine in Arabidopsis is mediated by RELEASE OF SILENCING 1 (ROS1), DEMETER (DME) and DEMETER-LIKE 2 and 3 (DML2 and DML3), which are able to demethylate DNA irrespective of the sequence context of the modified cytosine [16, 107]. Following further modification of the cleavage site, the gap is (filled and) repaired via DNA polymerase and ligase enzymes, e.g. LIG1 in Arabidopsis [19]. Again, linking histone- to DNA modification, DME requires histone linker H1 and the histone remodeling complex FACT to demethylate DNA in heterochromatin during reproduction [108110]. Similarly, the demethylase ROS1 is recruited to chromatin containing the histone variant H2A.Z. H2A.Z, in turn, is incorporated into histone octamers by the SWR1 complex, which is recruited to histone acetylation marks (via SWR1-associated MBD9 and NPX1) formerly added to the chromatin by the INCREASED DNA METHYLATION (IDM) complex [18, 111113].

We found numerous sugar beet homologs putatively mediating DNA demethylation in the described pathways. Their transcriptional changes upon cold treatment collectively support that the overall loss of DNA methylation upon cold exposure was predominantly caused by an increase of active DNA demethylation (Figure 8d): Cold triggered significant upregulation of NPX1 expression, and also significantly enhanced expression of two SWR1 components, i.e. PIE1 and SWC4, which overall might promote demethylation by ROS1 via deposition of H2A.Z. Finally, expression of ROS1 itself was significantly upregulated upon cold treatment. (Note, that expression of a DNA LIGASE 1-LIKE homolog was additionally upregulated more than four-fold. However, for conciseness, we excluded this gene from global analyses due to its location on a sequence component not assigned to any chromosome.)

Expression of a homolog of the putative DNA glycosylase At3g47830 | DEMETER-like protein 2 (DML2) was significantly lowered by cold. However, while disruption of the corresponding gene in Arabidopsis resulted in reduced methylation at sites that were heavily methylated in the wild-type, cytosine residues that were unmethylated or weakly methylated in WT, in fact, showed an increase in DNA methylation [114].

In several plant species, ROS1 expression decreases when RdDM is inhibited, for example through knock-out of a core RdDM-component [115, 116]. In contrast, despite simultaneous downregulation of several RdDM-genes including NRDP1, cold exposure resulted in upregulated ROS1 expression in sugar beets. This suggests that, either - instead of being linked to a general decrease of de novo methylation, the observed transcriptional repression of these RdDM components might rather reflect altered RdDM (for example via shift of targeted loci); or the regulatory mechanism that represses ROS1 in parallel to RdDM in other plants, is not conserved in sugar beets.

In fact, a short helitron TE in the ROS1 promoter of Arabidopsis, the so-called methylstat or MEMS (methylation monitoring sequence), which acts as a sensor for DNA methylation and (in this rather rare example) activates AtROS1 expression upon becoming methylated [115, 116], appears to be absent from the promoter of the sugar beet homolog. Besides, methylation of the cold induced ROS1 homolog of sugar beet (or of its promoter) was not significantly altered by cold.

Whereas gene promoters were overlapped by about equally many hyper- and hypomethylated DMRs in CpG context, about three quarters of all DMRs in CHH context coinciding with gene promoters were hypomethylated (Figure 3e). Hypomethylation was also clearly favored regarding promoter-associated DMCs, - particularly regarding non-CpG positions, with about twice as many promoter-associated DMCs in CHG, and more than six times as many cytosines in CHH context showing a significant reduction instead of an increase of methylation. The fact that we observed this particular association of hypomethylated DMCs with promoters, despite having a genome-wide larger number of hypermethylated DMCs in CHG, seems to imply that non-CpG hypomethylation is in some way specifically targeted towards promoters. As recently shown, ROS1 demethylates preferentially promoters of otherwise repressed genes [117]. As we detected a significantly upregulated expression of ROS1 due to cold in our set of DEGs, we consider demethylation by ROS1 as a possible explanation for the observed association of hypomethylated (non-CpG) DMCs with promoter regions. However, of the 33 genes carrying a hypomethylated DMC in non-CpG context in their promoter, only one is also differentially expressed. This gene is a homolog of the flowering pathway gene BBX32 from Arabidopsis, which in turn was shown to interact with CONSTANS-LIKE 3 (COL3) to target the promoter of FT [118], overall regulating the onset of flowering.

In contrast to RdDM-based methylation, CMTs show some specificity towards particular cytosine environments, i.e. they prefer to act on a cytosine directly preceding A and T rather than C [101]. Positions that were differentially methylated in response to cold treatment in our experiment revealed that hypermethylated CHG or CHH DMCs were almost completely devoid of another cytosine, particularly at the H directly following the cytosine that becomes (hyper)methylated (Additional file 1). In contrast, hypomethylated CHGs and CHHs were more tolerant towards additional cytosine residues. Based on Arabidopsis data, ROS1 has been proposed to counteract mainly RdDM to prevent spread of (TE-)methylation to genes [17, 119]. Under the assumption that hypomethylation is not primarily based on loss of RdDM, the increased frequency of CCG or CCH among hypo- compared to hypermethylated DMCs, fits with increased ROS1 activity counteracting methylation established through RdDM.

In summary, the transcriptional changes in sugar beet leaves in response to cold suggest an overall decrease of DNA methylation, mainly linked to enhanced active removal of methylated residues (mainly via ROS1).

Conclusion

A plants ability to adapt gene expression to changes in the environment can be crucial to their survival and/or development and it is thought to be fine-tuned by DNA methylation. Our study provided insights into conserved methylomic and transcriptomic responses of sugar beets to cold exposure. We propose that cold-dependent reduction of DNA methylation was mainly due to active removal of methylation marks through collectively upregulated expression of sugar beet homologs within the ROS1 pathway. Mechanistically, these effects seem to be common among different plants including sugar beet, tea [120], tomato [59], and poplar [60]. Strikingly, while an overall reduced DNA methylation can theoretically be facilitated either via increase of demethylation, or through passive loss, cold appears to predominantly trigger upregulation of an active DNA demethylation pathway in all of these species.

Methods

Plant material and growth conditions

The sequencing data presented in this study was obtained from two hybrid sugar beet genotypes [GT1 and GT2; corresponding to genotypes designated GT1 and GT3 in Martins Rodrigues et al. [4]]. Seed material of Beta vulgaris ssp. vulgaris GT1 and GT2 used in this study was provided by, and cultivation and collection of material was performed at KWS (Einbeck, Germany). Plants were germinated and grown on standard soil substrate ED73 (Einheitserdwerke Patzer, Germany) mixed with sand (10% v/v) under short day conditions (10 h light/14 h dark), at 60% relative humidity and 110 µmol m−2 s−1 light intensity (fluorescent tube light). GT1 and GT2 beets were either grown under control conditions (20°C during the day, 16°C at night) for 15 weeks, or - after initial cultivation under control conditions for 12 weeks – were transferred to 12°C for 7 days (acclimation phase), followed by 13 days at 4°C and finally 24h at 0°C (cold treatment). Leaf material was collected from three biological replicates per genotype and condition, transferred to liquid nitrogen and stored at -80°C until further processing for methylomic and transcriptomic analysis.

Sample preparation, sequencing and data availability

Sample extraction (total genomic DNA from leaves), library preparation and whole genome bisulfite sequencing (Illumina NovaSeq 6000 platform) was provided as a custom service (Novogene, Beijing, China). For WGBS, the PBAT protocol for paired-end sequencing (PE150) was used. Total RNA for RNA-seq was extracted from leaf material as described in Martins Rodrigues et al. [4]. Poly-A selection, library preparation and stranded, paired-end sequencing (151+151; Illumina HiSeq 2000) was provided as a custom service (GATC GmbH, Konstanz, Germany). Raw sequencing data (RNA-seq and Whole Genome Bisulfite Sequencing) have been deposited to NCBI’s Sequence Read Archive under BioProject ID PRJNA74855, accessible via https://www.ncbi.nlm.nih.gov/sra/PRJNA748559 upon acceptance of this article (also see paragraph “Availability of Data and Materials” in section “Declarations”).

Pre-processing and mapping

Raw sequence reads from both sequencing methods were inspected using fastQC [121] and multiQC [122]. Using BBDuk version 38.69 [123] the reads were filtered for N-reads, common contaminants, and known adaptor sequences from both ends of the reads using a k-mer length between 11 – 23. Reads were additionally filtered based on read quality (Q20) and length (minimum 35 bp). WGBS reads had an additional step after adaptor trimming in which the first 18 nucleotides of each read were removed (ftl = 18) using BBDuk version 38.69 [123] as recommended by the company, to avoid bias due to the sequencing technique. Trimming success was confirmed based on quality reports generated for trimmed data using fastQC and multiQC again. Adapter- and quality-trimmed RNA-seq reads were mapped to RefBeet1.2.2 (43, [RefBeet-1.2.fna.gz] downloaded from: 44), using STAR version 2.5.0a [125]. Cleaned WGBS reads were mapped to the same reference genome using BISMARK version 22.3 [126]. Duplicated reads within the WGBS data were removed by BISMARK’s deduplication function.

Methylation extraction and detection of DMCs and DMRs

The methylation status together with its positional coverage was evaluated based on BISMARK’s mapping alignments and methylation extractor function output (genome-wide cytosine report) for cytosines with a minimum read coverage of eight per position in each sample. Corresponding bam-files generated by BISMARK were used as input for detection of DMCs using the methylKit package [43] within Bioconductor [127] in R [128], which was used to detect differentially methylated cytosines (comparisons of COLD vs CONTROL, for each cytosine context separately). A cytosine site with a (by default SLIM adjusted) q-value < .05 was defined as differentially methylated. Coverage-filtered BISMARK methylation extraction outputs were used for the detection of DMRs via HOME, a histogram-based machine learning approach version 1.0.0 [45].

Analysis of differential gene expression

Mapped transcriptomic reads were quantified on gene-level for all predicted protein coding genes [ published in 42, BeetSet-2, available from 124 as BeetSet-2.unfiltered_genes.1408.gff3] using featureCounts version 1.5.0. provided with the Subread package [129]. The output was used to analyze differential expression between CONTROL and COLD using DESeq2 version 1.32.0 [47] with a reduced design formula, correcting for genotype specific differences in the transcriptomic cold response. Genes with an absolute Log2FoldChange ≥ 1 and an adjusted p-value (Bonferroni correction) of ≤ .05 were defined as differentially expressed (see Additional file 2).

Pseudochromosome construction and visualization

For plots showing distribution of detected cold-responsive elements - i.e. of differentially expressed genes (DEGs), differentially methylated cytosines (DMCs), or differentially methylated regions (DMRs) - along chromosomes, all localized scaffolds (designated BvchrX.scaYYY), followed by all unplaced scaffolds (BvchrX_un.scaYYY, where X denotes the corresponding chromosome and YYY indicates the [incremental] number of the scaffold) assigned to the same chromosome were strung together with 10kb pseudogaps inserted between individual sequence components (i.e. scaffolds). Based on the (cummulative) lengths of (preceding) scaffolds and pseudogaps, a position-adjustment table was generated, carrying - for each scaffold - a constant, from which relative positions of a given feature with respect to the corresponding (pseudo-)chromosome was calculated (based on the feature’s coordinates as given in the original annotation file, i.e. relative to its scaffold) directly prior to visualizing the data. If not otherwise specified, plots were generated using ggplot2 [130] in R [128] and finalized using Inkscape version 1.0.2-2 [131].

Functional classification and phylogenetic analysis

A fasta file containing protein sequences for all annotated BeetSet-2 genes [42] was obtained from the Beta vulgaris Resource [124]. The file was checked for consistency (valid format, valid character usage etc.) using the Mercator4 Fasta Validator and sequences were assigned to MapMan bins using Mercator4 v3.0, both available from the plaBi dataBase [132]. Data for all 285 gene products assigned to MapMan category 12 (“Chromatin organization”) was extracted. Genes not assigned to any chromosome were discarded, leaving 267 genes. Of those, 38 were differentially expressed in response to cold. Following closer inspection of the DEGs assigned to MapMan category “Chromatin organization”, literature search on those genes (with a focus on those with predicted effects on cytosine modification) revealed eleven additional DEGs with homologs related to DNA methylation or demethylation. For the complete list of genes (MapMan-based as well as manually assigned) see Additional file 3.

Phylogenetic analysis was done based on amino acid sequences of DME/ROS- and chromomethyltransferase homologs from Arabidopsis thaliana and Beta vulgaris, as compiled in PLAZA HOM04D001046, or HOM04D001291, respectively [133]. Additionally, the peptide sequence of another putative nuclease that we found among our DEGs, and which in NCBI annotation data was described as “DEMETER-LIKE 2”, was included with the DME/ROS-set. Analyses were done based on these sequences for each set, independently, using default parameters of PhyML+SMS/OneClick-workflow available via NGphylogeny [134, 135], and subsequent visualization in iToL [136].

Declarations

Ethics approval and consent to participate

Sugar beet hybrid lines GT1 and GT2 used in this study were selected by KFW and are proprietary to KWS SAAT SE & Co. KGaA (Einbeck, Germany). Seed material of Beta vulgaris ssp. vulgaris GT1 and GT2 was provided by, and cultivation as well as collection of sample material was performed at KWS SAAT SE & Co. KGaA (Einbeck, Germany) in compliance with institutional, national, and international guidelines and legislation. This article did not contain any studies with human participants or animals and did not involve any endangered or protected species.

Consent for publication

Not applicable.

Availability of data and materials

Beta vulgaris subsp. vulgaris lines used in this study can be obtained from KWS SAAT SE & Co. KGaA for non-commercial research purposes. Raw sequencing data (RNA-seq and Whole Genome Bisulfite Sequencing) have been deposited to NCBI’s Sequence Read Archive under BioProject ID PRJNA74855 and will be made available via https://www.ncbi.nlm.nih.gov/sra/PRJNA748559 upon acceptance of this manuscript. Reviewers can access all metadata via https://dataview.ncbi.nlm.nih.gov/object/PRJNA748559?reviewer=8u2fcim9s740odhk4tql5tc8tauntil data is released to the public upon acceptance. Refbeet1.2.2 [41] and corresponding annotations (BeetSet-2; 42) were downloaded from ‘The Beta vulgaris Resource’ [124] and served as reference genome during mapping and any downstream analyses.

Competing interests

OC, KF-W, WK, and FL are employed by KWS SAAT SE & Co. KGaA. This study further received (partial) funding from KWS SAAT SE & Co. KGaA. The funder was involved in selection and propagation of germplasm material and plant growth. All authors declare no other competing interests.

Funding

This work was supported by a research grant to HEN and US by the Federal Ministry of Education and Research, Germany (BMBF project “Betahiemis”; FKZ: 031B0185). IK is a PhD candidate funded by KWS SAAT SE & Co. KGaA. 

Authors' contributions

US, HEN, FL and WK acquired funding; JMC, US, FL, WK, US, HEN, BP, WZ, KH and CM designed the research. KFW selected sugar beet accessions used in this study. IK and CMR harvested the material. SG, WZ, IK, CMR and CM performed experiments; SG and CM analyzed the data; and SG, JMC and CM wrote the manuscript. All authors proof‐read and approved the manuscript.

Acknowledgements

We thank Michaela Brock, David Pscheidt and Jörg Hofmann (FAU) for excellent technical assistance and Wiebke Rettberg, Wenke Rettberg, and Joern Koch (KWS) for help and support with sugar beet harvesting. We further acknowledge Yun Jaeyoung (Novogene) for support and advice on WGBS sequencing. 

References

1. OECD/FAO. Sugar. In: OECD-FAO Agricultural Outlook 2021-2030. Paris: OECD Publishing; 2021. doi:10.1787/19428846-en.

2. Hoffmann CM, Kluge-Severin S. Light absorption and radiation use efficiency of autumn and spring sown sugar beets. F Crop Res. 2010;119:238–44.

3. Hoffmann CM, Kluge-Severin S. Growth analysis of autumn and spring sown sugar beet. Eur J Agron. 2011;34:1–9. doi:10.1016/j.eja.2010.09.001.

4. Martins Rodrigues C, Müdsam C, Keller I, Zierer W, Czarnecki O, Corral JM, et al. Vernalization alters sink and source identities and reverses phloem translocation from taproots to shoots in sugar beet. Plant Cell. 2020;32:3206–23.

5. Zhang X, Yazaki J, Sundaresan A, Cokus S, Chan SW-L, Chen H, et al. Genome-wide high-resolution mapping and functional analysis of DNA methylation in Arabidopsis. Cell. 2006;126:1189–201. doi:10.1016/j.cell.2006.08.003.

6. Matzke MA, Mosher RA. RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nat Rev Genet. 2014;15:394–408. doi:10.1038/nrg3683.

7. Kankel MW, Ramsey DE, Stokes TL, Flowers SK, Haag JR, Jeddeloh JA, et al. Arabidopsis MET1 cytosine methyltransferase mutants. Genetics. 2003;163:1109–22.

8. Kim J, Kim JH, Richards EJ, Chung KM, Woo HR. Arabidopsis VIM proteins regulate epigenetic silencing by modulating DNA methylation and histone modification in cooperation with MET1. Mol Plant. 2014;7:1470–85. doi:10.1093/mp/ssu079.

9. Du J, Johnson LM, Jacobsen SE, Patel DJ. DNA methylation pathways and their crosstalk with histone methylation. Nat Rev Mol Cell Biol. 2015;16:519–32. doi:10.1038/nrm4043.

10. Yaari R, Noy-Malka C, Wiedemann G, Auerbach Gershovitz N, Reski R, Katz A, et al. DNA METHYLTRANSFERASE 1 is involved in mCG and mCCG DNA methylation and is essential for sporophyte development in Physcomitrella patens. Plant Mol Biol. 2015;88:387–400. doi:10.1007/s11103-015-0328-8.

11. Zemach A, Kim MY, Hsieh P-H, Coleman-Derr D, Eshed-Williams L, Thao K, et al. The Arabidopsis nucleosome remodeler DDM1 allows DNA methyltransferases to access H1-containing heterochromatin. Cell. 2013;153:193–205. doi:10.1016/j.cell.2013.02.033.

12. Stroud H, Greenberg MVC, Feng S, Bernatavichute Y V, Jacobsen SE. Comprehensive analysis of silencing mutants reveals complex regulation of the Arabidopsis methylome. Cell. 2013;152:352–64. doi:10.1016/j.cell.2012.10.054.

13. Wendte JM, Schmitz RJ. Specifications of targeting heterochromatin modifications in plants. Mol Plant. 2018;11:381–7. doi:10.1016/j.molp.2017.10.002.

14. Agius F, Kapoor A, Zhu J-K. Role of the Arabidopsis DNA glycosylase/lyase ROS1 in active DNA demethylation. Proc Natl Acad Sci. 2006;103:11796–801. doi:10.1073/pnas.0603563103.

15. Gong Z, Morales-Ruiz T, Ariza RR, Roldán-Arjona T, David L, Zhu J-K. ROS1, a repressor of transcriptional gene silencing in Arabidopsis, encodes a DNA glycosylase/lyase. Cell. 2002;111:803–14. doi:10.1016/s0092-8674(02)01133-9.

16. Morales-Ruiz T, Ortega-Galisteo AP, Ponferrada-Marín MI, Martínez-Macías MI, Ariza RR, Roldán-Arjona T. DEMETER and REPRESSOR OF SILENCING 1 encode 5-methylcytosine DNA glycosylases. Proc Natl Acad Sci U S A. 2006;103:6853–8.

17. Zhu J-K. Active DNA demethylation mediated by DNA glycosylases. Annu Rev Genet. 2009;43:143–66. doi:10.1146/annurev-genet-102108-134205.

18. Liu R, Lang Z. The mechanism and function of active DNA demethylation in plants. J Integr Plant Biol. 2020;62:148–59.

19. Li Y, Duan C-G, Zhu X, Qian W, Zhu J-K. A DNA ligase required for active DNA demethylation and genomic imprinting in Arabidopsis. Cell Res. 2015;25:757–60. doi:10.1038/cr.2015.45.

20. Martin L, Fei Z, Giovannoni J, Rose JKC. Catalyzing plant science research with RNA-seq. Front Plant Sci. 2013;4. doi:10.3389/fpls.2013.00066.

21. Cokus SJ, Feng S, Zhang X, Chen Z, Merriman B, Haudenschild CD, et al. Shotgun bisulphite sequencing of the Arabidopsis genome reveals DNA methylation patterning. Nature. 2008;452:215–9. doi:10.1038/nature06745.

22. Lister R, O’Malley RC, Tonti-Filippini J, Gregory BD, Berry CC, Millar AH, et al. Highly integrated single-base resolution maps of the epigenome in Arabidopsis. Cell. 2008;133:523–36. doi:10.1016/j.cell.2008.03.029.

23. Boyko A, Kovalchuk I. Epigenetic control of plant stress response. Environ Mol Mutagen. 2008;49:61–72. doi:https://doi.org/10.1002/em.20347.

24. Chinnusamy V, Zhu JK. Epigenetic regulation of stress responses in plants. Curr Opin Plant Biol. 2009;12:133–9.

25. Pandey G, Yadav CB, Sahu PP, Muthamilarasan M, Prasad M. Salinity induced differential methylation patterns in contrasting cultivars of foxtail millet (Setaria italica L.). Plant Cell Rep. 2017;36:759–72. doi:10.1007/s00299-016-2093-9.

26. Annacondia ML, Magerøy MH, Martinez G. Stress response regulation by epigenetic mechanisms: changing of the guards. Physiol Plant. 2018;162:239–50. doi:https://doi.org/10.1111/ppl.12662.

27. Thiebaut F, Hemerly AS, Ferreira PCG. A role for epigenetic regulation in the adaptation and stress responses of non-model plants. Front Plant Sci. 2019;10. doi:10.3389/fpls.2019.00246.

28. Al-Harrasi I, Al-Yahyai R, Yaish MW. Differential DNA methylation and transcription profiles in date palm roots exposed to salinity. PLoS One. 2018;13:e0191492. doi:10.1371/journal.pone.0191492.

29. Cheng J, Niu Q, Zhang B, Chen K, Yang R, Zhu J-K, et al. Downregulation of RdDM during strawberry fruit ripening. Genome Biol. 2018;19:212. doi:10.1186/s13059-018-1587-x.

30. Eichten SR, Springer NM. Minimal evidence for consistent changes in maize DNA methylation patterns following environmental stress. Front Plant Sci. 2015;6. doi:10.3389/fpls.2015.00308.

31. Fan SK, Ye JY, Zhang LL, Chen HS, Zhang HH, Zhu YX, et al. Inhibition of DNA demethylation enhances plant tolerance to cadmium toxicity by improving iron nutrition. Plant Cell Environ. 2019;43:275–91. doi:10.1111/pce.13670.

32. Li R, Hu F, Li B, Zhang Y, Chen M, Fan T, et al. Whole genome bisulfite sequencing methylome analysis of mulberry (Morus alba) reveals epigenome modifications in response to drought stress. Sci Rep. 2020;10:8013. doi:10.1038/s41598-020-64975-5.

33. Liu T, Li Y, Duan W, Huang F, Hou X. Cold acclimation alters DNA methylation patterns and confers tolerance to heat and increases growth rate in Brassica rapa. J Exp Bot. 2017;68:1213–24. doi:10.1093/jxb/erw496.

34. Lu X, Wang X, Chen X, Shu N, Wang J, Wang D, et al. Single-base resolution methylomes of upland cotton (Gossypium hirsutum L.) reveal epigenome modifications in response to drought stress. BMC Genomics. 2017;18:297. doi:10.1186/s12864-017-3681-y.

35. Song Q-X, Lu X, Li Q-T, Chen H, Hu X-Y, Ma B, et al. Genome-wide analysis of DNA methylation in soybean. Mol Plant. 2013;6:1961–74. doi:10.1093/mp/sst123.

36. Sun S, Zhu J, Guo R, Whelan J, Shou H. DNA methylation is involved in acclimation to iron-deficiency in rice (Oryza sativa). Plant J. 2021. https://onlinelibrary.wiley.com/doi/abs/10.1111/tpj.15318.

37. Tang X, Wang Q, Yuan H, Huang X. Chilling-induced DNA demethylation is associated with the cold tolerance of Hevea brasiliensis. BMC Plant Biol. 2018;18:70. doi:10.1186/s12870-018-1276-7.

38. Hébrard C, Peterson DG, Willems G, Delaunay A, Jesson B, Lefèbvre M, et al. Epigenomics and bolting tolerance in sugar beet genotypes. J Exp Bot. 2016;67:207–25. doi:10.1093/jxb/erv449.

39. Zakrzewski F, Schubert V, Viehoever P, Minoche AE, Dohm JC, Himmelbauer H, et al. The CHH motif in sugar beet satellite DNA: a modulator for cytosine methylation. Plant J. 2014;78:937–50. doi:https://doi.org/10.1111/tpj.12519.

40. Zakrzewski F, Schmidt M, Van Lijsebettens M, Schmidt T. DNA methylation of retrotransposons, DNA transposons and genes in sugar beet (Beta vulgaris L.). Plant J. 2017;90:1156–75. doi:10.1111/tpj.13526.

41. Dohm JC, Minoche AE, Holtgräwe D, Capella-Gutiérrez S, Zakrzewski F, Tafer H, et al. The genome of the recently domesticated crop plant sugar beet ( {Beta} vulgaris ). Nature. 2014;505:546–9. doi:10.1038/nature12817.

42. Minoche AE, Dohm JC, Schneider J, Holtgräwe D, Viehöver P, Montfort M, et al. Exploiting single-molecule transcript sequencing for eukaryotic gene prediction. Genome Biol. 2015;16:184. doi:10.1186/s13059-015-0729-7.

43. Akalin A, Kormaksson M, Li S, Garrett-Bakelman FE, Figueroa ME, Melnick A, et al. methylKit: a comprehensive R package for the analysis of genome-wide DNA methylation profiles. Genome Biol. 2012;13:R87. doi:10.1186/gb-2012-13-10-r87.

44. Vanyushin BF, Ashapkin V V. DNA methylation in higher plants: Past, present and future. Biochim Biophys Acta - Gene Regul Mech. 2011;1809:360–8. doi:10.1016/j.bbagrm.2011.04.006.

45. Srivastava A, Karpievitch Y V, Eichten SR, Borevitz JO, Lister R. HOME: A histogram based machine learning approach for effective identification of differentially methylated regions. BMC Bioinformatics. 2019;20:253. doi:10.1186/s12859-019-2845-y.

46. Teng C-S, Wu B-H, Yen M-R, Chen P-Y. MethGET: web-based bioinformatics software for correlating genome-wide DNA methylation and gene expression. BMC Genomics. 2020;21:375. doi:10.1186/s12864-020-6722-x.

47. Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15. doi:10.1186/s13059-014-0550-8.

48. Hiraguri A, Itoh R, Kondo N, Nomura Y, Aizawa D, Murai Y, et al. Specific interactions between Dicer-like proteins and HYL1/DRB- family dsRNA-binding proteins in Arabidopsis thaliana. Plant Mol Biol. 2005;57:173–88. doi:10.1007/s11103-004-6853-5.

49. Nakazawa Y, Hiraguri A, Moriyama H, Fukuhara T. The dsRNA-binding protein DRB4 interacts with the Dicer-like protein DCL4 in vivo and functions in the trans-acting siRNA pathway. Plant Mol Biol. 2007;63:777–85. doi:10.1007/s11103-006-9125-8.

50. Stroud H, Do T, Du J, Zhong X, Feng S, Johnson L, et al. Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis. Nat Struct Mol Biol. 2014;21:64–72. doi:10.1038/nsmb.2735.

51. Shen X, De Jonge J, Forsberg SKG, Pettersson ME, Sheng Z, Hennig L, et al. Natural CMT2 variation is associated with genome-wide methylation changes and temperature seasonality. PLoS Genet. 2014;10:e1004842. doi:10.1371/journal.pgen.1004842.

52. Niederhuth CE, Bewick AJ, Ji L, Alabady MS, Kim K Do, Li Q, et al. Widespread natural variation of DNA methylation within angiosperms. Genome Biol. 2016;17:194. doi:10.1186/s13059-016-1059-0.

53. Gayacharan C, Joel AJ. Epigenetic responses to drought stress in rice (Oryza sativa L.). Physiol Mol Biol Plants. 2013;19.

54. Yaish MW, Al-Lawati A, Al-Harrasi I, Patankar HV. Genome-wide DNA methylation analysis in response to salinity in the model plant caliph medic (Medicago truncatula). BMC Genomics. 2018;19:78. doi:10.1186/s12864-018-4484-5.

55. Steward N, Ito M, Yamaguchi Y, Koizumi N, Sano H. Periodic DNA methylation in maize nucleosomes and demethylation by environmental stress. J Biol Chem. 2002;277:37741–6.

56. Song Y, Jia Z, Hou Y, Ma X, Li L, Jin X, et al. Roles of DNA methylation in cold priming in tartary buckwheat. Front Plant Sci. 2020;11.

57. Chen B, Saltveit ME, Beckles DM. Chilling-stress modifies DNA methylation level in cucumber (Cucumis sativus L.) seedling radicle to regulate elongation rate. Sci Hortic (Amsterdam). 2019;252:14–9. doi:10.1016/j.scienta.2019.03.023.

58. Sicilia A, Scialò E, Puglisi I, Lo Piero AR. Anthocyanin biosynthesis and DNA methylation dynamics in sweet orange fruit [Citrus sinensis L. (Osbeck)] under cold stress. J Agric Food Chem. 2020;68:7024–31. doi:10.1021/acs.jafc.0c02360.

59. Zhang B, Tieman DM, Jiao C, Xu Y, Chen K, Fei Z, et al. Chilling-induced tomato flavor loss is associated with altered volatile synthesis and transient changes in DNA methylation. Proc Natl Acad Sci. 2016;113:12580–5. doi:10.1073/pnas.1613910113.

60. Conde D, Le Gac A-L, Perales M, Dervinis C, Kirst M, Maury S, et al. Chilling-responsive DEMETER-LIKE DNA demethylase mediates in poplar bud break. Plant Cell Environ. 2017;40:2236–49.

61. Wang H, Beyene G, Zhai J, Feng S, Fahlgren N, Taylor NJ, et al. CG gene body DNA methylation changes and evolution of duplicated genes in cassava. Proc Natl Acad Sci. 2015;112:13729–34. doi:10.1073/pnas.1519067112.

62. Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science (80- ). 2010;328:916–9.

63. Li X, Zhu J, Hu F, Ge S, Ye M, Xiang H, et al. Single-base resolution maps of cultivated and wild rice methylomes and regulatory roles of DNA methylation in plant gene expression. BMC Genomics. 2012;13:300. doi:10.1186/1471-2164-13-300.

64. Ezaki B, Higashi A, Nanba N, Nishiuchi T. An S-adenosyl Methionine Synthetase (SAMS) gene from Andropogon virginicus L. confers aluminum stress tolerance and facilitates epigenetic gene regulation in Arabidopsis thaliana. Front Plant Sci. 2016;7. doi:10.3389/fpls.2016.01627.

65. Heidari P, Mazloomi F, Nussbaumer T, Barcaccia G. Insights into the SAM synthetase gene family and its roles in tomato seedlings under abiotic stresses and hormone treatments. Plants. 2020;9:586. doi:10.3390/plants9050586.

66. Yu J-G, Lee G-H, Park Y-D. Physiological role of endogenous S-adenosyl-L-methionine synthetase in Chinese cabbage. Hortic Environ Biotechnol. 2012;53:247–55. doi:10.1007/s13580-012-0021-7.

67. Ma C, Wang Y, Gu D, Nan J, Chen S, Li H. Overexpression of S-adenosyl-L-methionine synthetase 2 from sugar beet M14 increased arabidopsis tolerance to salt and oxidative stress. Int J Mol Sci. 2017;18:847. doi:10.3390/ijms18040847.

68. Meng J, Wang L, Wang J, Zhao X, Cheng J, Yu W, et al. METHIONINE ADENOSYLTRANSFERASE4 mediates DNA and histone methylation. Plant Physiol. 2018;177:652–70. doi:10.1104/pp.18.00183.

69. Eamens A, Vaistij FE, Jones L. NRPD1a and NRPD1b are required to maintain post-transcriptional RNA silencing and RNA-directed DNA methylation in Arabidopsis. Plant J. 2008;55:596–606. doi:10.1111/j.1365-313X.2008.03525.x.

70. He X-J, Hsu Y-F, Pontes O, Zhu J, Lu J, Bressan RA, et al. NRPD4, a protein related to the RPB4 subunit of RNA polymerase II, is a component of RNA polymerases IV and V and is required for RNA-directed DNA methylation. Genes Dev. 2009;23:318–30. doi:10.1101/gad.1765209.

71. Herr AJ, Jensen MB, Dalmay T, Baulcombe DC. RNA polymerase IV directs silencing of endogenous DNA. Science (80- ). 2005;308:118–20. doi:10.1126/science.1106910.

72. Onodera Y, Haag JR, Ream T, Nunes PC, Pontes O, Pikaard CS. Plant nuclear RNA polymerase IV mediates siRNA and DNA methylation- dependent heterochromatin formation. Cell. 2005;120:613–22. doi:10.1016/j.cell.2005.02.007.

73. Yang D-L, Zhang G, Wang L, Li J, Xu D, Di C, et al. Four putative SWI2/SNF2 chromatin remodelers have dual roles in regulating DNA methylation in Arabidopsis. Cell Discov. 2018;4:1–14. doi:10.1038/s41421-018-0056-8.

74. Zhou M, Palanca AMS, Law JA. Locus-specific control of the de novo DNA methylation pathway in Arabidopsis by the CLASSY family. Nat Genet. 2018;50:865–73. doi:10.1038/s41588-018-0115-y.

75. Law JA, Du J, Hale CJ, Feng S, Krajewski K, Palanca AMS, et al. Polymerase IV occupancy at RNA-directed DNA methylation sites requires SHH1. Nature. 2013;498:385–9.

76. Law JA, Vashisht AA, Wohlschlegel JA, Jacobsen SE. SHH1, a homeodomain protein required for DNA methylation, as well as RDR2, RDM4, and chromatin remodeling factors, associate with RNA polymerase IV. PLOS Genet. 2011;7:e1002195. doi:10.1371/journal.pgen.1002195.

77. Zhang H, Ma Z-Y, Zeng L, Tanaka K, Zhang C-J, Ma J, et al. DTF1 is a core component of RNA-directed DNA methylation and may assist in the recruitment of Pol IV. Proc Natl Acad Sci. 2013;110:8290–5. doi:10.1073/pnas.1300585110.

78. Greenberg MVC, Ausin I, Chan SWL, Cokus SJ, Cuperus JT, Feng S, et al. Identification of genes required for de novo DNA methylation in Arabidopsis. Epigenetics. 2011;6:344–54. doi:10.4161/epi.6.3.14242.

79. Singh J, Mishra V, Wang F, Huang HY, Pikaard CS. Reaction mechanisms of Pol IV, RDR2, and DCL3 drive RNA channeling in the siRNA-directed DNA methylation pathway. Mol Cell. 2019;75:576-589.e5. doi:10.1016/j.molcel.2019.07.008.

80. Haag JR, Ream TS, Marasco M, Nicora CD, Norbeck AD, Pasa-Tolic L, et al. In vitro transcription activities of Pol IV, Pol V, and RDR2 reveal coupling of Pol IV and RDR2 for dsRNA synthesis in plant RNA silencing. Mol Cell. 2012;48:811–8. doi:10.1016/j.molcel.2012.09.027.

81. Liu Q, Feng Y, Zhu Z. Dicer-like (DCL) proteins in plants. Funct Integr Genomics. 2009;9:277–86.

82. Xie Z, Johansen LK, Gustafson AM, Kasschau KD, Lellis AD, Zilberman D, et al. Genetic and functional diversification of small RNA pathways in plants. PLoS Biol. 2004;2:e104.

83. Zilberman D, Cao X, Jacobsen SE. ARGONAUTE4 control of locus-specific siRNA accumulation and DNA and histone methylation. Science (80- ). 2003;299:716–9. doi:10.1126/science.1079695.

84. Havecker ER, Wallbridge LM, Hardcastle TJ, Bush MS, Kelly KA, Dunn RM, et al. The Arabidopsis RNA-directed DNA methylation Argonautes functionally diverge based on their expression and interaction with target loci. Plant Cell. 2010;22:321–34. doi:10.1105/tpc.109.072199.

85. Carbonell A. Plant ARGONAUTEs: Features, Functions, and Unknowns. Methods Mol Biol. 2017;1640:1–21.

86. Cao X, Jacobsen SE. Role of the Arabidopsis DRM methyltransferases in de novo DNA methylation and gene silencing. Curr Biol. 2002;12:1138–44. doi:10.1016/S0960-9822(02)00925-9.

87. Zhong X, Du J, Hale CJ, Gallego-Bartolome J, Feng S, Vashisht AA, et al. Molecular mechanism of action of plant DRM de novo DNA methyltransferases. Cell. 2014;157:1050–60. doi:10.1016/j.cell.2014.03.056.

88. Zhang H, Lang Z, Zhu JK. Dynamics and function of DNA methylation in plants. Nature Reviews Molecular Cell Biology. 2018;19:489–506.

89. Haag JR, Pikaard CS. Multisubunit RNA polymerases IV and V: purveyors of non-coding RNA for plant gene silencing. Nat Rev Mol Cell Biol. 2011;12:483–92. doi:10.1038/nrm3152.

90. Ream TS, Haag JR, Wierzbicki AT, Nicora CD, Norbeck AD, Zhu JK, et al. Subunit compositions of the RNA-silencing enzymes Pol IV and Pol V reveal their origins as specialized forms of RNA polymerase II. Mol Cell. 2009;33:192–203. doi:10.1016/j.molcel.2008.12.015.

91. Zhou M, Law JA. RNA Pol IV and V in gene silencing: Rebel polymerases evolving away from Pol II’s rules. Curr Opin Plant Biol. 2015;27:154–64. doi:10.1016/j.pbi.2015.07.005.

92. Cuerda-Gil D, Slotkin RK. Non-canonical RNA-directed DNA methylation. Nat Plants. 2016;2:1–8. doi:10.1038/nplants.2016.163.

93. Ausin I, Mockler TC, Chory J, Jacobsen SE. IDN1 and IDN2 are required for de novo DNA methylation in Arabidopsis thaliana. Nat Struct Mol Biol. 2009;16:1325–7. doi:10.1038/nsmb.1690.

94. Finke A, Kuhlmann M, Mette MF. IDN2 has a role downstream of siRNA formation in RNA-directed DNA methylation. Epigenetics. 2012;7:950–60.

95. Greenberg MVC, Deleris A, Hale CJ, Liu A, Feng S, Jacobsen SE. Interplay between active chromatin marks and RNA-directed DNA methylation in Arabidopsis thaliana. PLoS Genet. 2013;9:e1003946. doi:10.1371/journal.pgen.1003946.

96. Jiang D, Yang W, He Y, Amasino RM. Relatives of the human Lysine-Specific Demethylase1 repress the expression of and and thus promote the floral transition. Plant Cell. 2007;19:2975–87. doi:10.1105/tpc.107.052373.

97. Henderson IR, Deleris A, Wong W, Zhong X, Chin HG, Horwitz GA, et al. The de novo cytosine methyltransferase DRM2 requires intact UBA domains and a catalytically mutated paralog DRM3 during RNA-directed DNA methylation in Arabidopsis thaliana. PLoS Genet. 2010;6:e1001182.

98. Zhong X, Hale CJ, Nguyen M, Ausin I, Groth M, Hetzel J, et al. Domains rearranged methyltransferase3 controls DNA methylation and regulates RNA polymerase V transcript abundance in Arabidopsis. Proc Natl Acad Sci U S A. 2015;112:911–6.

99. Costa-Nunes P, Kim JY, Hong E, Pontes O. The cytological and molecular role of DOMAINS REARRANGED METHYLTRANSFERASE3 in RNA-dependent DNA methylation of Arabidopsis thaliana. BMC Res Notes. 2014;7:721. doi:10.1186/1756-0500-7-721.

100. Zhang Y, Harris CJ, Liu Q, Liu W, Ausin I, Long Y, et al. Large-scale comparative epigenomics reveals hierarchical regulation of non-CG methylation in Arabidopsis. Proc Natl Acad Sci. 2018;115:E1069--E1074. doi:10.1073/pnas.1716300115.

101. Gouil Q, Baulcombe DC. DNA methylation signatures of the plant chromomethyltransferases. PLoS Genet. 2016;12:e1006526. doi:10.1371/journal.pgen.1006526.

102. Ebbs ML, Bartee L, Bender J. H3 lysine 9 methylation is maintained on a transcribed inverted repeat by combined action of SUVH6 and SUVH4 methyltransferases. Mol Cell Biol. 2005;25:10507–15.

103. Jackson JP, Johnson L, Jasencakova Z, Zhang X, PerezBurgos L, Singh PB, et al. Dimethylation of histone H3 lysine 9 is a critical mark for DNA methylation and gene silencing in Arabidopsis thaliana. Chromosoma. 2004;112:308–15. doi:10.1007/s00412-004-0275-7.

104. Johnson LM, Bostick M, Zhang X, Kraft E, Henderson I, Callis J, et al. The SRA methyl-cytosine-binding domain links DNA and histone methylation. Curr Biol. 2007;17:379–84. doi:10.1016/j.cub.2007.01.009.

105. Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11:204–20. doi:10.1038/nrg2719.

106. Li X, Harris CJ, Zhong Z, Chen W, Liu R, Jia B, et al. Mechanistic insights into plant SUVH family H3K9 methyltransferases and their binding to context-biased non-CG DNA methylation. Proc Natl Acad Sci U S A. 2018;115:E8793--E8802.

107. Penterman J, Zilberman D, Huh JH, Ballinger T, Henikoff S, Fischer RL. DNA demethylation in the Arabidopsis genome. Proc Natl Acad Sci U S A. 2007;104:6752–7.

108. Frost JM, Kim MY, Park GT, Hsieh PH, Nakamura M, Lin SJH, et al. FACT complex is required for DNA demethylation at heterochromatin during reproduction in Arabidopsis. PNAS. 2018;115:E4720–9. doi:10.1073/pnas.1713333115.

109. He S, Vickers M, Zhang J, Feng X. Natural depletion of histone H1 in sex cells causes DNA demethylation, heterochromatin decondensation and transposon activation. Elife. 2019;8.

110. Rea M, Zheng W, Chen M, Braud C, Bhangu D, Rognan TN, et al. Histone H1 affects gene imprinting and DNA methylation in Arabidopsis. Plant J. 2012;71:776–86. doi:10.1111/j.1365-313X.2012.05028.x.

111. Lang Z, Lei M, Wang X, Tang K, Miki D, Zhang H, et al. The Methyl-CpG-Binding Protein {MBD7} facilitates active DNA demethylation to limit DNA hyper-methylation and transcriptional gene silencing. Mol Cell. 2015;57:971–83. doi:10.1016/j.molcel.2015.01.009.

112. Li Y, Kumar S, Qian W. Active DNA demethylation: mechanism and role in plant development. Plant Cell Rep. 2018;37:77–85. doi:10.1007/s00299-017-2215-z.

113. Nie W-F, Lei M, Zhang M, Tang K, Huang H, Zhang C, et al. Histone acetylation recruits the SWR1 complex to regulate active DNA demethylation in Arabidopsis. Proc Natl Acad Sci. 2019;116:16641–50. doi:10.1073/pnas.1906023116.

114. Ortega-Galisteo AP, Morales-Ruiz T, Ariza RR, Roldán-Arjona T. Arabidopsis DEMETER-LIKE proteins DML2 and DML3 are required for appropriate distribution of DNA methylation marks. Plant Mol Biol. 2008;67:671–81. doi:10.1007/s11103-008-9346-0.

115. Lei M, Zhang H, Julian R, Tang K, Xie S, Zhu J-K. Regulatory link between DNA methylation and active demethylation in Arabidopsis. Proc Natl Acad Sci. 2015;112:3553–7. doi:10.1073/pnas.1502279112.

116. Williams BP, Pignatta D, Henikoff S, Gehring M. Methylation-sensitive expression of a DNA demethylase gene serves as an epigenetic rheostat. PLOS Genet. 2015;11:e1005142. doi:10.1371/journal.pgen.1005142.

117. Halter T, Wang J, Amesefe D, Lastrucci E, Charvin M, Singla Rastogi M, et al. The Arabidopsis active demethylase ROS1 cis-regulates defence genes by erasing DNA methylation at promoter-regulatory regions. Elife. 2021;10:e62994. doi:10.7554/eLife.62994.

118. Tripathi P, Carvallo M, Hamilton EE, Preuss S, Kay SA. Arabidopsis B-BOX32 interacts with CONSTANS-LIKE3 to regulate flowering. Proc Natl Acad Sci. 2017;114:172–7. doi:10.1073/pnas.1616459114.

119. Tang K, Lang Z, Zhang H, Zhu J-K. The DNA demethylase ROS1 targets genomic regions with distinct chromatin modifications. Nat Plants. 2016;2:1–10. doi:10.1038/nplants.2016.169.

120. Tong W, Li R, Huang J, Zhao H, Ge R, Wu Q, et al. Divergent DNA methylation contributes to duplicated gene evolution and chilling response in tea plants. Plant J. 2021;106:1312–27. doi:10.1111/tpj.15237.

121. Andrews S, Krueger F, Seconds-Pichon A, Biggins F, Wingett S. FastQC. A quality control tool for high throughput sequence data. Babraham Institute. 2015;1:1. https://www.bioinformatics.babraham.ac.uk/projects/ fastqc/.

122. Ewels P, Magnusson M, Lundin S, Käller M. MultiQC: summarize analysis results for multiple tools and samples in a single report. Bioinformatics. 2016;32:3047–8. doi:10.1093/bioinformatics/btw354.

123. Bushnell B. BBMap. https://sourceforge.net/projects/bbmap/.

124. Himmelbauer H, Dohm JC, Minoche AE. The Beta vulgaris Resource - sugar beet genetics and genomics. https://bvseq.boku.ac.at/. Accessed 12 Feb 2021.

125. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, Jha S, et al. STAR: Ultrafast universal RNA-seq aligner. Bioinformatics. 2013;29:15–21. doi:10.1093/bioinformatics/bts635.

126. Krueger F, Andrews SR. Bismark: A flexible aligner and methylation caller for Bisulfite-Seq applications. Bioinformatics. 2011;27:1571–2. doi:10.1093/bioinformatics/btr167.

127. Huber W, Carey VJ, Gentleman R, Anders S, Carlson M, Carvalho BS, et al. Orchestrating high-throughput genomic analysis with Bioconductor. Nat Methods. 2015;12:115–21. doi:10.1038/nmeth.3252.

128. R Core Team. R: A Language and Environment for Statistical Computing. 2021;2. http://www.r-project.org.

129. Liao Y, Smyth GK, Shi W. featureCounts: an efficient general-purpose read summarization program. Bioinformatics. 2014;30:923–30.

130. Wickham H. ggplot2: elegant graphics for data analysis. Springer-Verlag New York; 2016. https://ggplot2.tidyverse.org.

131. Inkscape Project. Inkscape. 2021. https://inkscape.org.

132. Schwacke R, Y. Ponce-Soto G, Krause K, M. Bolger A, Arsova B, Hallab A, et al. MapMan4: a refined protein classification and annotation framework applicable to multi-omics data analysis. Mol Plant. 2019;12:879–92. doi:https://doi.org/10.1016/j.molp.2019.01.003.

133. Van Bel M, Diels T, Vancaester E, Kreft L, Botzki A, de Peer Y, et al. PLAZA 4.0: an integrative resource for functional, evolutionary and comparative plant genomics. Nucleic Acids Res. 2018;46:D1190--D1196. doi:10.1093/nar/gkx1002.

134. Dereeper A, Guignon V, Blanc G, Audic S, Buffet S, Chevenet F, et al. Phylogeny.fr: robust phylogenetic analysis for the non-specialist. Nucleic Acids Res. 2008;36 Web Server issue:W465–9. doi:10.1093/nar/gkn180.

135. Lemoine F, Correia D, Lefort V, Doppelt-Azeroual O, Mareuil F, Cohen-Boulakia S, et al. NGPhylogeny.fr: new generation phylogenetic services for non-specialists. Nucleic Acids Res. 2019;47:W260–5. doi:10.1093/nar/gkz303.

136. Letunic I, Bork P. Interactive Tree Of Life (iTOL) v4: recent updates and new developments. Nucleic Acids Res. 2019;47:W256--W259. doi:10.1093/nar/gkz239.