- Open Access
Multi-omics data integration reveals link between epigenetic modifications and gene expression in sugar beet (Beta vulgaris subsp. vulgaris) in response to cold
BMC Genomics volume 23, Article number: 144 (2022)
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.
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.
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.
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 . 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% . 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 (Fig. 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 . 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 . 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 [9,10,11,12]. 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 . CHH methylation is mediated chromatin-specifically - either by DRM2-dependent RdDM (euchromatin), or - in H1-containing heterochromatin where RdDM is blocked - by CMT2 [11, 14, 15]. 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 [16,17,18,19,20,21].
RNA sequencing has already been used in numerous studies to investigate altered gene expression under environmental stress in plants . 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 [24, 25]. In crops, this technique has also been used to identify epigenetic mechanisms that have a significant impact on stress responses or developmental processes [26,27,28,29,30]. 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 [31,32,33,34,35,36,37,38,39,40]. Hébrard et al.  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. [42, 43] 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.
Experimental setup and workflow
As depicted in the overview of our workflow (Fig. 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 × 3 biological replicates each) per condition (see Methods). From this material, DNA and RNA was extracted to subsequently analyze both the methylomic and transcriptomic cold responses (Fig. 1a).
Cold treatment reduces DNA methylation levels
The 12 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 . Mapping and deduplication revealed average proportions of uniquely mapped and multi-mapped reads of 60 and 12%, respectively (see Fig. 1a and Methods).
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 (Fig. 2a). For beets grown under control conditions, the highest methylation levels were detected for cytosines in CpG context (mCpG: 66.8%), followed by mCHG (41.7%) and mCHH (9.4%). Cold treatment slightly decreased global methylation levels of cytosines in all sequence contexts, with methylation levels reduced by 2.8% (mCpG: 64%), 1.2% (mCHG: 40.5%), or 1.3% (mCHH: 8.7%) in the cold treated group, although differences were not statistically significant in any case (two-tailed Welch’s test; Fig. 2a). 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 [22, 44]: an upstream region spanning 2.5 kb in 5′ direction of the transcription start site (TSS); 5′-untranslated region (UTR); coding sequence (CDS), intron and 3’UTR. As shown in the right panel of Fig. 2, methylation levels generally reached maxima upstream of genes and within introns. In detail, mCHG and mCHH levels peaked extragenically, i.e. in upstream regions, while maximum mCpG was detected along introns. Overall, CDSs and 5’UTRs showed comparatively low methylation levels in all contexts, with mCHG and mCHH reaching a minimum in CDSs, and lowest mCpG detected for 5’UTRs. Although COLD showed slightly lower average methylation than CONTROL in almost all gene-related features (i.e. in all categories including 5′ flanks) and for all contexts analyzed, these differences were mostly rather subtle, while statistically significant differences were detected only between 5′- or 3′- UTR mCHH levels of COLD vs CONTROL (two-tailed Welch’s test, p < .05).
Exposure to cold alters methylation of individual cytosines in all sequence contexts
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 . 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 (Fig. 1c). We identified 5319 DMCs distributed over the nine chromosomes (Fig. 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; Fig. 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 (Fig. 3a). Of all detected DMCs, 2873 (54.0%) were either located directly within, or occurred within 2.5 kb 5′ of TSSs, with dmCpG accounting for more than 80% of all DMCs overlapping a gene or its upstream flank, approximately reflecting the proportion of DMCs in CpG context to non-CpG-DMCs.
Overall, most DMCs overlapped introns, CDSs, or occurred within 2.5 kb upstream of TSSs - together accounting for more than 90% of all gene-associated or 5′-proximal DMCs (Fig. 3b), except for dmCHHhypo, which showed a preference for upstream regions and 3’UTRs. Additionally, dmCHGhyper and dmCHHhyper showed particularly strong associations with introns. Note, that while this was beyond the scope of our study, we additional provide corresponding data summarizing overlaps of DMCs (and DMRs, see next chapter) with transposable elements (TE) in Additional File 1, Supplementary Fig. S2 and S3. In plants, methyltransferases and demethylases favor specific sequences for local (de)methylation . To evaluate possible sequence preferences, we scanned the sequence composition within 25 bp-flanking-regions (25 bp up- and downstream of each DMC, i.e. 51 bp in total) of all cold-dependent DMCs, separately for hyper- and hypomethylation in CpG, CHG, or CHH context, respectively (see Additional file 1, Supplementary Fig. S4). 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 [24, 25].
Differentially methylated regions in CHH context preferentially become hypomethylated in response to cold
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 . 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 (Fig. 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 CDSs and introns and lowest in 5’UTRs (Fig. 3e). DMRs in CHG context were similarly distributed over gene features, but peaked within introns, both contrasting CHH DMRs, which predominantly overlapped upstream regions but – apart from some introns - were basically absent from other gene features. Comparable to DMCs, most of the DMRs were located on chromosome 6, while chromosome 1 showed the fewest DMRs (Fig. 3f).
CHG and CHH methylation is associated with low gene expression
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 Fig. 1a and Methods).
Integration of gene expression and DNA methylation additionally 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 . 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 groups, based on its relative expression level. Subsequently, the methylation level along all genes within the same group, 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 Fig. 4 (also see Additional File 1, Supplementary Fig. S5), 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, closer inspection revealed small but noticeable differences between 5′ mCHH peak levels of CONTROL and COLD samples (indicated in Fig. 4, CHH panel), with gene groups of COLD samples reaching slightly lower maximum methylation levels than gene groups with comparable expression levels of CONTROL samples. But because gene groups are reconstituted for each sample individually and thus their composition partially differs between samples, an effect of 5′ mCHH levels on gene expression cannot be inferred from this observation.
Cold triggers genotype-independent changes in gene expression
Analysis of differential expression  contrasting CONTROL vs. COLD yielded in total 2549 DEGs (Additional File 2), most of them localized on chromosome 6, fewest on chromosome 8 (Fig. 5). Of them, 1244 were up-regulated, while 1305 were down-regulated in COLD (Wald-test; padj ≤.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 were either directly overlapped by one or more DMCs, or had at least one DMC co-localizing with their 2.5 kb upstream regions. Comparing the proportion of DMC-associated genes among the DEGs (8.2%) with the fraction of DMC-associated genes (again with gene-ranges extended to comprise an additional 2.5 kb 5′ of TSSs) among all annotated genes (6.6%, corresponding to a total of 1712 DMC-associated genes), association with differential methylation was significantly enriched among our set of DEGs [χ2(1, n = 26,004) = 4.9091, p < .05]. However, combinations of hyper- or hypomethylation and up- or downregulation of expression did not reflect an apparent tendency (Fig. 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 upstream regions (Fig. 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 (Fig. 6b). Of all genes, 2287 (≙ 8.8%) showed association with a DMR, and comparing the absolute numbers, DEGs were again enriched with DMRs [χ2(1, n = 26,004) = 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.
Altered expression of genes involved in DNA (de)methylation correlates with observed methylation changes
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 (see Additional File 2 and section Functional classification and phylogenetic analysis in Methods for details on the classification, annotation and pathway assignments of DEGs) and whose expression was significantly altered in response to cold (see Fig. 7 a-c for expression levels of these genes in CONTROL or COLD, respectively; and see Fig. 8 for a schematic model depicting pathways, in which the highlighted DEGs are predicted to interact with additional [non-deregulated] genes to eventually drive alterations of DNA methylation in response to cold). 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 [51, 52]. Furthermore, cold significantly decreased the expression of IDN2 (INVOLVED IN DE NOVO 2; Bv2_024880_hakm) - an interactor of DRM2 required for RdDM - as well as of a homolog of DNA (cytosine-5)-methyltransferase CMT2 (Bv3_050080_yren, also see Fig. 7d and f), which mediates de novo methylation preferentially at cytosines in CHH context [15, 53]. However, phylogenetic analysis suggests that another gene, whose expression was not cold-dependent nor consistent between genotypes (Fig. 7d, f), but which showed partially high expression (comparable to Bv3_050080_yren), might act as the actual functional homolog of CMT2 at least in one of the genotypes, or act redundantly to the DEG candidate, Bv3_050080_yren. 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 demethylation 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. Phylogenetic analyses and the almost complete lack of expression of other sugar beet homologs from the same gene family (Fig. 7e and g) contributed to the functional classification of this gene as the main ROS1-homolog. (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 (Fig. 7).
The empirical results reported here should be considered in the light of some limitations. A large proportion of the genome of Beta vulgaris (42.3% of the Refbeet assembly, ) consists of repetitive sequences, rendering it a challenging target particularly during read mapping. Furthermore, the same reference was used to map reads to - irrespective of the treatment. This means that potential genomic rearrangements (e.g. transpositions) that had occurred as a response to the applied stress, are inevitably missed during the analysis. This is mainly relevant for analyses of WGBS data, which (as opposed to RNA-seq) theoretically covers the entire genome. Retrotransposons are the most abundant type of repetitive elements identified in Refbeet1.2.2 , and activation as well as transposition, in fact, has been reported to be triggered by abiotic stresses, including cold [54, 55]. Although this was beyond the scope of our study, and has been examined in detail in earlier studies [42, 43], we detected significant decrease of mCHH along several retrotransposon subtypes (mainly DNA and LTR, Additional File 1 Supplementary Fig. S2 and S3), which could potentially be linked to transcriptional activation of transposons or nearby genes. The fact that we used Refbeet1.2.2 as reference, which is based on an independent sugar beet accession (KWS2320, ), further complicates proper alignment of reads, as well as interpretation of WGBS data. Galewski and McGrath  recently examined lineage-specific variation (LSV), comparing closely related Beta vulgaris crop-types and different Beta vulgaris accessions of the same subspecies between or within groups. With about 0.23%, the described LSV between accessions of the same subspecies (sugar beet) was low compared to the high variation between crop types (99.37%) - as was the variation observed within other crop type lineages (approximately 0.40%). A comparatively large fraction of variation distinguishing different sugar beet accessions, however, was concentrated on chromosome 6, which could be due to drift or divergent selection on this chromosome. It is conceivable that the ability to change methylation status at this chromosome in response to the environment could contribute to a fine-tuned regulation of breeding-relevant genes. And this would be consistent with the observed high DMC density on chromosome 6 (Fig. 3c) following cold exposure. However, genetic variation (distinguishing the accessions used in this study from the reference) can affect the outcome of WGBS analysis, for one, because it negatively impacts mapping efficiency, and secondly, because methylation is eventually determined based on comparison of the read to the reference. This means that a C (reference genome) to T mutation present in both accessions, is eventually scored as a non-methylated cytosine (which are converted to T during bisulfite treatment). In turn, for positions representing a T to C mutation, where the C could be methylated or non-methylated in the analyzed accessions, no methylation information will be extracted, resulting in loss of information for such positions. In the present study we mainly focused on the investigation of genotype-independent adaptation of Beta vulgaris subsp. vulgaris to cold exposure. By eventually comparing control grown plants with cold treated plants of the same accession(s), genotype-dependent responses and variation should be cancelled out, at least partially. While this comes with a price, as the combination of two accessions eventually dilutes effects that are pronounced only in one of the genotypes, it should also decrease artifacts arising from variation between our accessions and the reference, for example during the examination of differential methylation between COLD and CONTROL. Yet, we did detect several low-difference DMCs (Fig. 3c) and found that the majority of them occurred highly concentrated within rather narrow sequence segments. Upon further inspection, we could assign most of these regions to areas, where coverage was highly above average. The seemingly high coverage, in turn, could be attributed to comparatively large numbers of falsely mapped reads in those regions, which represented highly repetitive sequences. More precisely, low-difference-high-coverage “DMC” clusters occurred almost exclusively within sequences annotated as or predicted to code for rRNAs (see also Additional File 1, Supplementary Fig. S2 and S3).
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 (Fig. 2). This deviates from the results obtained by Niederhuth et al. , 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.  and with the methylation proportion for B. vulgaris detected by Zakrzewski et al. . 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 Niederhuth et al. ), 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 fine-tunes global methylation levels 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 . In upland cotton (Gossypium hirsutum), for example, Lu et al.  found a decrease of 2% in mCpG and mCHG, and an increase of 4% in mCHH, after drought stress. Gayacharan and Joel  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 , but increased methylation in Medicago truncatula . Cold stress was found to trigger DNA demethylation in roots of maize seedlings . 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 . In contrast, cucumber (Cucumis sativus) radicles increased DNA methylation in response to cold treatment, which was paralleled by reduced growth rates . 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 [40, 64, 65]. Finally, in some woody perennials like poplar, cold triggers DNA demethylation prior to bud break .
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 100 bp upstream and 500 bp downstream of TSS) and the TTS (approximately 500 bp upstream and 100 bp downstream of TTS), mCHG and mCHH are basically absent from genes . 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 [24, 38, 68]. Further, mCpG at the TSS and TTS tends to inversely correlate with expression levels of the corresponding gene [67, 69]. 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 2 kb 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. , 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 (Fig. 2; Additional File 1, Supplementary Fig. S1). 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 (Fig. 4 and Fig. 6). This means that while absolute methylation levels in CHG and CHH context do, in fact, inversely correlate with gene expression (Fig. 4), changes in methylation (i.e., DMCs or DMRs) are not necessarily associated with a significant change in expression (Fig. 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 (also see “Functional classification and phylogenetic analysis” in “Methods” and Additional File 2) 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 (Fig. 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 [70,71,72,73], 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 . 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 (Fig. 8b). Canonical RdDM is initiated through the activity of DNA-directed RNA-polymerase IV [75,76,77,78], 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 [79, 80] - triggers binding of Pol IV [81,82,83]. 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 [84,85,86], which can be cleaved into siRNAs by DICER-LIKE PROTEINs [87, 88], before being loaded to ARGONAUTEs (mainly AGO4 and 6 [89, 90]). At the target site, nascent transcripts produced by another core RdDM-polymerase, Pol V, can pair with appropriate AGO-coupled siRNAs . This finally recruits DOMAINS REARRANGED METHYLASE 2 (DRM2), an interactor of AGO4, to methylate cytosines at the target site [92,93,94].
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 Figs. 7 and 8, 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 [95,96,97]. In addition, expression of CLSY1 was significantly decreased in cold-treated sugar beets, as was the expression of three RDRs (RDR1, 3, 6; Fig. 8b) of which at least two participate in non-canonical RdDM at stages of alternative siRNA production in Arabidopsis . Furthermore, cold decreased expression of DRB4, which is thought to assist processing of dsRNA fragments to siRNA via DCL4 [51, 52], 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 [99, 100]. It is thought to normally promote DRM2-mediated DNA methylation at some target loci [101, 102].
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-CpG 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 [103,104,105]. 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 (Fig. 8b). And while AGO2 and Pol II – both of which were upregulated in cold-treated samples (Fig. 8b) - can be involved in providing siRNAs that eventually guide DNA methylation via DRM2 through non-canonical, partially Pol IV-independent mechanisms of RdDM , 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 (Fig. 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 Fig. 3d). DRM2 is able to methylate cytosines independent of their sequence context [92, 106]. 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 (Fig. 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 ). 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, 15, 108,109,110,111].
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 (Fig. 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 [110, 112]. 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 Fig. 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 . 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 (Fig. 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 [16,17,18,19]. 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 [18, 113]. Following further modification of the cleavage site, the gap is (filled and) repaired via DNA polymerase and ligase enzymes, e.g. LIG1 in Arabidopsis . Again, linking histone- to DNA modification, DME requires histone linker H1 and the histone remodeling complex FACT to demethylate DNA in heterochromatin during reproduction [114,115,116]. 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 [20, 117,118,119].
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 (Fig. 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.
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 .
In several plant species, ROS1 expression decreases when RdDM is inhibited, for example through knock-out of a core RdDM-component [121, 122]. 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 [121, 122], 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 upstream region comprising the promoter) was not significantly altered by cold.
Whereas upstream regions were overlapped by about equally many hyper- and hypomethylated DMRs in CpG context, about three quarters of all DMRs in CHH context coinciding with upstream regions were hypomethylated (Fig. 3e). Hypomethylation was also clearly favored regarding DMCs occurring within 5′ flanks, − particularly regarding non-CpG positions, with about twice as many upstream 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 upstream regions, 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 those areas. As recently shown, ROS1 demethylates preferentially promoters of otherwise repressed genes . 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 5′ regulatory regions. However, of the 33 genes carrying a hypomethylated DMC in non-CpG context in their upstream flanks, 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 , 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 . 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, Supplementary Fig. S4). 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 [19, 125]. 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).
A plant’s ability to adapt gene expression to changes in the environment can be crucial to its 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 , tomato , and poplar . 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.
Plant material and growth conditions
The sequencing data presented in this study was obtained from two hybrid sugar beet genotypes . Seed material of Beta vulgaris ssp. vulgaris GT1 (KWS-accession: 1NB0218) and GT2 (KWS-accession: 1NB0133) used in this study were provided by, and cultivation and collection of material was performed at KWS SAAT SE & Co. KGaA (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 24 h 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 (post-bisulfite adapter tagging, ) 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. . 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 (also see paragraph Availability of Data and Materials in section Declarations). All software and parameters described in the following sections are additionally summarized in Additional File 3.
Pre-processing and mapping
Raw sequence reads from both sequencing methods were inspected using fastQC  and multiQC . Reads were then filtered for N-reads, common contaminants, and known adapter sequences from both ends of the reads using BBDuk (version 38.69, ); with k-mer length between 11 and 23). Reads were additionally filtered based on read quality (Q20) and length (minlen = 35 bp). WGBS reads had an additional step after adapter trimming in which the first 18 nucleotides of each read were removed (ftl = 18) using BBDuk (v38.69, ) 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 (RefBeet-1.2.fna.gz, ; downloaded from ), using STAR (v2.5.0a, ). Cleaned WGBS reads were mapped to the same reference genome using BISMARK (v22.3, ). 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 (v1.19.0, ) within Bioconductor (v3.14, ) in R (v4.1.2 “Bird Hippie”, ), 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 (v1.0.0, ).
Analysis of differential gene expression
Mapped transcriptomic reads were quantified on gene-level for all predicted protein coding genes (BeetSet-2.unfiltered_genes.1408.gff3, [44, 50]) using featureCounts v1.5.0. provided with the Subread package . The output was used to analyze differential expression between CONTROL and COLD using DESeq2 (v1.39.0, ) 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 3 for main code snippets and non-default parameter settings).
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 10 kb pseudogaps inserted between individual sequence components (i.e. scaffolds). Based on the (cumulative) 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  in R  and finalized using Inkscape (v1.0.2–2, ).
Functional classification and phylogenetic analysis
A fasta file containing protein sequences for all annotated BeetSet-2 genes  was obtained from ‘The Beta vulgaris Resource’ . 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’ . 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 11 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 . 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 [140, 141], and subsequent visualization in iToL .
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 (submission ID: SUB9992759) under BioProject ID PRJNA748559, available via https://www.ncbi.nlm.nih.gov/bioproject/PRJNA748559. Refbeet1.2.2  and corresponding annotations (BeetSet-2, ) were downloaded from ‘The Beta vulgaris Resource’  and served as reference genome during mapping and any downstream analyses.
OECD/FAO. Sugar. In: OECD-FAO Agricultural Outlook 2021–2030. Paris: OECD Publishing; 2021. p. 150–62.
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.
Hoffmann CM, Kluge-Severin S. Growth analysis of autumn and spring sown sugar beet. Eur J Agron. 2011;34:1–9.
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.
Zhang X, Yazaki J, Sundaresan A, Cokus SJ, Chan SWL, Chen H, et al. Genome-wide high-resolution mapping and functional analysis of DNA methylation in Arabidopsis. Cell. 2006;126:1189–201.
Matzke MA, Mosher RA. RNA-directed DNA methylation: an epigenetic pathway of increasing complexity. Nat Rev Genet. 2014;15:394–408.
Kankel MW, Ramsey DE, Stokes TL, Flowers SK, Haag JR, Jeddeloh JA, et al. Arabidopsis MET1 cytosine methyltransferase mutants. Genetics. 2003;163:1109–22.
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.
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.
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.
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.
Stroud H, Greenberg MVC, Feng S, Bernatavichute YV, Jacobsen SE. Comprehensive analysis of silencing mutants reveals complex regulation of the Arabidopsis methylome. Cell. 2013;152:352–64.
Wendte JM, Schmitz RJ. Specifications of targeting heterochromatin modifications in plants. Mol Plant. 2018;11:381–7.
Choi J, Lyons DB, Zilberman D. Histone H1 prevents non-CG methylation-mediated small RNA biogenesis in Arabidopsis heterochromatin. Elife. 2021;10:e72676.
Stroud H, Do T, Du J, Zhong X, Feng S, Johnson LM, et al. Non-CG methylation patterns shape the epigenetic landscape in Arabidopsis. Nat Struct Mol Biol. 2014;21:64–72.
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.
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.
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.
Zhu J-K. Active DNA demethylation mediated by DNA glycosylases. Annu Rev Genet. 2009;43:143–66.
Liu R, Lang Z. The mechanism and function of active DNA demethylation in plants. J Integr Plant Biol. 2020;62:148–59.
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.
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.
Martin L, Fei Z, Giovannoni JJ, Rose JKC. Catalyzing plant science research with RNA-seq. Front Plant Sci. 2013;4:66.
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.
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.
Boyko A, Kovalchuk I. Epigenetic control of plant stress response. Environ Mol Mutagen. 2008;49:61–72.
Chinnusamy V, Zhu J-K. Epigenetic regulation of stress responses in plants. Curr Opin Plant Biol. 2009;12:133–9.
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.
Annacondia ML, Magerøy MH, Martinez G. Stress response regulation by epigenetic mechanisms: changing of the guards. Physiol Plant. 2018;162:239–50.
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:246.
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.
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.
Eichten SR, Springer NM. Minimal evidence for consistent changes in maize DNA methylation patterns following environmental stress. Front Plant Sci. 2015;6:308.
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.
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.
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.
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.
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.
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;107(3):727–39.
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.
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.
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.
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.
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.
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.
Vanyushin BF, Ashapkin VV. DNA methylation in higher plants: past, present and future. Biochim Biophys Acta Gene Regul Mech. 2011;1809:360–8.
Srivastava A, Karpievitch YV, 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.
Teng C-S, Wu B-H, Yen M-R, Chen PY. MethGET: web-based bioinformatics software for correlating genome-wide DNA methylation and gene expression. BMC Genomics. 2020;21:375.
Love MI, Huber W, Anders S. Moderated estimation of fold change and dispersion for RNA-seq data with DESeq2. Genome Biol. 2014;15:550.
Himmelbauer H, Dohm JC, Minoche AE. The Beta vulgaris Resource - sugar beet genetics and genomics https://bvseq.boku.ac.at/. Accessed 12 Feb 2021.
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.
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.
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.
Cavrak VV, Lettner N, Jamge S, Kosarewicz A, Bayer LM, Mittelsten SO. How a retrotransposon exploits the plant’s heat stress response for its activation. PLoS Genet. 2014;10:e1004115.
Liang Z, Anderson SN, Noshay JM, Crisp PA, Enders TA, Springer NM. Genetic and epigenetic variation in transposable element expression responses to abiotic stress in maize. Plant Physiol. 2021;186:420–33.
Galewski P, McGrath JM. Genetic diversity among cultivated beets (Beta vulgaris) assessed via population-based whole genome sequences. BMC Genomics. 2020;21(1). https://doi.org/10.1186/s12864-020-6451-1.
Niederhuth CE, Bewick AJ, Ji L, Alabady MS, Do KK, Li Q, et al. Widespread natural variation of DNA methylation within angiosperms. Genome Biol. 2016;17:194.
Sudan J, Raina M, Singh R. Plant epigenetic mechanisms: role in abiotic stress and their generational heritability. 3 Biotech. 2018;8:1–12.
Gayacharan C, Joel AJ. Epigenetic responses to drought stress in rice (Oryza sativa L.). Physiol Mol Biol Plants. 2013;19:379–87.
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.
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.
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:608540.
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.
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.
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.
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.
Zemach A, McDaniel IE, Silva P, Zilberman D. Genome-wide evolutionary analysis of eukaryotic DNA methylation. Science. 2010;328:916–9.
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.
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.
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:1627.
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.
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.
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.
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.
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.
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.
Herr AJ, Jensen MB, Dalmay T, Baulcombe DC. RNA polymerase IV directs silencing of endogenous DNA. Science. 2005;308:118–20.
Onodera Y, Haag JR, Ream TS, Nunes PC, Pontes O, Pikaard CS. Plant nuclear RNA polymerase IV mediates siRNA and DNA methylation- dependent heterochromatin formation. Cell. 2005;120:613–22.
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.
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.
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.
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.
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.
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.
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.
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.
Liu Q, Feng Y, Zhu Z. Dicer-like (DCL) proteins in plants. Funct Integr Genomics. 2009;9:277–86.
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.
Zilberman D, Cao X, Jacobsen SE. ARGONAUTE4 control of locus-specific siRNA accumulation and DNA and histone methylation. Science. 2003;299:716–9.
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.
Carbonell A. Plant ARGONAUTEs: features, functions, and unknowns. Methods Mol Biol. 2017;1640:1–21.
Cao X, Jacobsen SE. Role of the Arabidopsis DRM methyltransferases in de novo DNA methylation and gene silencing. Curr Biol. 2002;12:1138–44.
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.
Zhang H, Lang Z, Zhu J-K. Dynamics and function of DNA methylation in plants. Nat Rev Mol Cell Biol. 2018;19:489–506.
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.
Ream TS, Haag JR, Wierzbicki AT, Nicora CD, Norbeck AD, Zhu J-K, 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.
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.
Cuerda-Gil D, Slotkin RK. Non-canonical RNA-directed DNA methylation. Nat Plants. 2016;2:1–8.
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.
Finke A, Kuhlmann M, Mette MF. IDN2 has a role downstream of siRNA formation in RNA-directed DNA methylation. Epigenetics. 2012;7:950–60.
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.
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.
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.
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.
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.
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–74.
Gouil Q, Baulcombe DC. DNA methylation signatures of the plant chromomethyltransferases. PLoS Genet. 2016;12:e1006526.
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.
Jackson JP, Johnson LM, 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.
Johnson LM, Bostick M, Zhang X, Kraft E, Henderson IR, Callis J, et al. The SRA methyl-cytosine-binding domain links DNA and histone methylation. Curr Biol. 2007;17:379–84.
Law JA, Jacobsen SE. Establishing, maintaining and modifying DNA methylation patterns in plants and animals. Nat Rev Genet. 2010;11:204–20.
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–802.
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.
Frost JM, Kim MY, Park GT, Hsieh P-H, Nakamura M, Lin SJH, et al. FACT complex is required for DNA demethylation at heterochromatin during reproduction in Arabidopsis. Proc Natl Acad Sci. 2018;115:E4720–9.
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:e42530.
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.
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.
Li Y, Kumar S, Qian W. Active DNA demethylation: mechanism and role in plant development. Plant Cell Rep. 2018;37:77–85.
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.
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.
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.
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.
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.
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.
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.
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.
Miura F, Enomoto Y, Dairiki R, Ito T. Amplification-free whole-genome bisulfite sequencing by post-bisulfite adaptor tagging. Nucleic Acids Res. 2012;40:e136.
Andrews S. FASTQC. A quality control tool for high throughput sequence data. 2010. Available online at http://www.bioinformatics.babraham.ac.uk/projects/fastqc/.
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.
Bushnell B. BBMap: a fast, accurate, splice-aware aligner; 2014.
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.
Krueger F, Andrews SR. Bismark: a flexible aligner and methylation caller for bisulfite-Seq applications. Bioinformatics. 2011;27:1571–2.
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.
R Core Team. R: a language and environment for statistical computing; 2021. p. 2.
Liao Y, Smyth GK, Shi W. featureCounts: an efficient general-purpose read summarization program. Bioinformatics. 2014;30:923–30.
Inkscape v0.91. http://inkscape.org/en/download. Accessed 16 Sept 2015.
Inkscape Project. Inkscape. 2021.
Schwacke R, Ponce-Soto GY, Krause K, Bolger AM, 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.
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–6.
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.
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.
Letunic I, Bork P. Interactive tree of life (iTOL) v4: recent updates and new developments. Nucleic Acids Res. 2019;47:W256–9.
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.
Open Access funding enabled and organized by Projekt DEAL. 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.
Ethics approval and consent to participate
Sugar beet hybrid lines GT1 (KWS-accession: 1NB0218) and GT2 (KWS-accession: 1NB0133) 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
OC, KF-W, WK, and FL are employed by KWS SAAT SE & Co. KGaA (Einbeck, Germany). 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.
Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.
Additional file 1: Supplementary Figure S1.
Addition to Fig. 2: Methylation levels in the genome and in gene components under control conditions or after exposure to cold. Supplementary Figure S2. Addition to Fig. 2: Methylation levels of transposable and repetitive elements under control conditions or after exposure to cold. Supplementary Figure S3. Addition to Fig. 3: Overlaps of differential methylation with transposable elements. Supplementary Figure S4. Sequence composition in the neighborhood of fully methylated, non-methylated or differentially methylated cytosines. Supplementary Figure S5. Addition to Fig. 4: Methylation profiles along genes grouped by their expression level.
Additional file 2: Supplementary Table S1.
Merged annotation and Arabidopsis homologs of all DEGs. Supplementary Table S2. Databases used for annotation of DEGs and assignment of Arabidopsis homologs.
Additional file 3.
Tools and Implementation.
Rights and permissions
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.
About this article
Cite this article
Gutschker, S., Corral, J.M., Schmiedl, A. et al. Multi-omics data integration reveals link between epigenetic modifications and gene expression in sugar beet (Beta vulgaris subsp. vulgaris) in response to cold. BMC Genomics 23, 144 (2022). https://doi.org/10.1186/s12864-022-08312-2
- Beta vulgaris subsp. vulgaris
- Abiotic stress
- Cold response
- DNA methylation