Roles of DNA Methylation in Cold Priming in Tartary Buckwheat

Plants experience a wide array of environmental stimuli, some of which are frequent occurrences of cold weather, which have priming effects on agricultural production and agronomic traits. DNA methylation may act as an epigenetic regulator for the cold response of Tartary buckwheat (Fagopyrum tataricum). Combined with long-term field observation and laboratory experiments, comparative phenome, methylome, and transcriptome analyses were performed to investigate the potential epigenetic contributions for the cold priming of Tartary buckwheat variety Dingku1. Tartary buckwheat cv. Dingku1 exhibited low-temperature resistance. Single-base resolution maps of the DNA methylome were generated, and a global loss of DNA methylation was observed during cold responding in Dingku1. These sites with differential methylation levels were predominant in the intergenic regions. Several hundred genes had different DNA methylation patterns and expressions in different cold treatments (cold memory and cold shock), such as CuAO, RPB1, and DHE1. The application of a DNA methylation inhibitor caused a change of the free lysine content, suggesting that DNA methylation can affect metabolite accumulation for Tartary buckwheat cold responses. The results of the present study suggest important roles of DNA methylation in regulating cold response and forming agronomic traits in Tartary buckwheat.


INTRODUCTION
Plants can precisely perceive hypothermia through epigenetic regulation with short-term cold stress responses (Jung et al., 2013) or prolonged cold temperature changes (Luo and He, 2020). DNA methylation can function as an epigenetic regulator to potentially provide flexible genomic parameters for plants to respond to various cold stresses (Chinnusamy and Zhu, 2009;Huff and Zilberman, 2012;Sahu et al., 2013;Thiebaut et al., 2019). The stressful experiences of plants can change how they subsequently respond so that they have stronger stress tolerance when encountering sudden environmental changes in the future. In higher plants, this is known as "stress memory" or "stress imprinting" (Bruce et al., 2007), and the expression of stress memory genes is modulated by epigenetic mechanisms (D'Urso and Brickner, 2014;Chang et al., 2020). There have been many studies of stress memory, but the original report conducted by Zuther et al. (2019) proposed that cold memory can improve plant freezing tolerance by changes in gene expression and lipid and metabolite composition and defined the memory of cold acclimation as cold priming. Subsequently, there have been several reports of the involvement of epigenetics in plant stress memory (Ding et al., 2012;Yang et al., 2020b). However, the detailed mechanism of how epigenetic memory is involved in frozen memory has not yet been described.
Tartary buckwheat (Fagopyrum tataricum) is a pseudocereal that belongs to the genus Fagopyrum within the Polygonaceae family. Tartary buckwheat is strongly adapted to growth in adverse environments (such as harsh climates and nutrientpoor soils; Zhang et al., 2017Zhang et al., , 2018Zhou et al., 2017b). Additionally, Tartary buckwheat is a short-generation and a diploid, with highly enriched flavonoid products, facilitating its use as a potential model species to study low-temperature adaptability in plants. Recent studies in Tartary buckwheat have focused on the functional analysis of individual genes and developmental traits (Liu et al., 2018a(Liu et al., ,b, 2019bHe et al., 2019;Yao et al., 2019), but there have been no reports of global epigenetic regulation.
Low temperature or repeated diurnal temperature difference promotes the accumulation of flavonoids, mainly because low temperature greatly increases the activities of enzymes in the flavonoid synthesis pathway (Christie et al., 1994;Leyva et al., 1995;Caldwell et al., 2005), suggesting that cold-induced transcriptional events can lead to desirable agronomic traits in plants. Metabolites can also protect plant cell viability during adversity. Epigenetic information in "cis-memory" is stored as a state of local chromatin (e.g., by DNA methylation or histone modification) and "trans-memory" exists as movable factors (e.g., transcriptional repressors; Dean, 2017). Some studies have confirmed that epigenetic memory formed under stress may allow for the quick adaptation of plants to ambient temperature changes. The phenotypic variation that can be induced by epigenetic memory -but not hereditary variation -is important to cope with rapid changes in the environment (Latzel et al., 2016). Few details are known of how epigenetic memory allows for a response to a changing environment, and Tartary buckwheat is an ideal plant to investigate the related mechanisms.
In this study, long-term field testing and laboratory experiments revealed that Tartary buckwheat cv. Dingku1 presents frost resistance and other agronomic traits, such as higher germination rate, higher water content of seeds, and higher flavonoid content. Experiments were designed as different cold treatments, including cold memory (cold priming) and cold shock, and phenome, methylome, and transcriptome analyses were performed to investigate the frost resistance characteristics of Dingku1. The single-base resolution DNA methylomes of Dingku1 were characterized under different cold treatments and revealed the global loss of DNA methylation, with some locally hypermethylated sites. The Kyoto Encyclopedia of Genes and Genomes (KEGG) enrichment analysis revealed that DNA methylation significantly impacted the pathways of lysine degradation; pyrimidine metabolism; and the synthesis of isoquinoline alkaloids, metabolism-related genes (i.e.,FtDHE1,FtRPB,and FtCuAO), and metabolites (i.e., lysine) to precisely regulate frozen memory. Treatment with a DNA methylation inhibitor interfered with the lysine level during the cold response, indicating that DNA methylation is critical for proper stress responses in Tartary buckwheat. These findings provide comprehensive insights into the development of cold priming and suggest guidelines for future breeding efforts in Tartary buckwheat.

Field Experimental Site
Field experiments were conducted at the China Oat and Buckwheat Industrial Technology System Haidong Comprehensive Test Station. Tartary buckwheat varieties that were recently bred to be frost tolerant were used.
The core test base of the Haidong Comprehensive Test Station is located in Huangzhong, Qinghai Province; this site has a flat terrain and a cold and humid climate (Figure 1). The soil organic matter content was 98 g/kg, the content of the available N was 24.2 mg/kg, the available P (P 2 O 5 ) was 13.43 mg/kg, the available potassium (K 2 O) was 110.32 mg/kg, and the pH was 8.4 (2016-2019 data). Meteorological data were provided by the Huangzhong Meteorological Bureau. Agronomic parameters [average yield per mu, plant height, length of tillering, thousand kernel weight (TKW), and grain yield] were measured as described . Tartary buckwheat seeds were collected after maturity, and total flavonoids, total protein, total starch, and crude fat were determined according to national standards (GB5009.5-2016) and technical documents of the Standardization of Shenzhen City (SZDB/Z 349-2019).

Plant Materials and Cold Treatments in Laboratory Experiments
After 4-8 h of soaking the seeds in ddH 2 O, the seeds were disinfected using a 15% NaClO solution and then placed in a culture dish with two layers of gauze. The culture dish was moved to a greenhouse and cultured until germination. Threeweek-old seedlings were treated with 5-aza-2'-deoxycytidine (aza-dC, Sigma) that was added to the liquid nutrient medium (7 g/ml). Seedlings were divided into two groups in cold stress experiments. The first group was to test the changes in plant morphology, plant water content, and plant temperature before and after cold stress. The experiment was set at 4°C for 24 h, followed by recovery at room temperature (21°C) for 1 day. The second group was designed to simulate the field climate and was used to assess the character of different varieties. The experiment included plants subjected to cold memory (4°C for 6 h, then at room temperature for 18 h, repeated four times, and then placed at 0°C for 6 h), cold stock (not acclimated, directly exposed to 0°C for 6 h), and control groups with normal growth conditions.

High Throughput Phenotypic Tests in Laboratory Experiments
Indoor high-throughput phenotypic observations were performed using image processing based on visible light (morphological test), infrared light (relative temperature measurement), and near-infrared light (relative water content test). A commercial phenotyping system (Scanalyzer3D, LemnaTec GmbH, Würselen, Germany) was used for image acquisition. Images of each plant were taken from the top and side views. Images were acquired and data preprocessing was organized into LemnaBase, a central database interface for the phenotyping system (Guo et al., 2017). Imaging, data acquisition, and data analysis of seeds and plants before and after cold treatment were carried out. Indicators of area, circumference, expansion degree, density, area of external polygons, circumferential length of external polygons, vertical length, horizontal length, and minimum diameter were performed to determine phenotypic differences in detail Zhou et al., 2017a).
The chlorophyll values of plants were determined (Polyphenol-Chlorophyll Meter-Dualex Scientific+, Force-A, France) using the second and third fully expanded leaf near the plant center. Three indexes were measured: chlorophyll absorptivity (chlorophyll, Chl), anthocyanidin, and NBI (Cerovic et al., 2012). The one-way variance was calculated for multiple comparisons using SPSS 19.0, GraphPad Prim7 software and Microsoft Excel 2016.

The Anti-freezing Physiological Index
Electrolyte leakage tests were performed as previously described . Leaves were transferred to 15 ml tubes and placed in a freezer (XT5201-D31-R40C, XuTemp, China). The plants were exposed to freezing temperatures ranging from 10 to −6°C, and leaves without damage were then immersed in 10 ml ultrapure water (Milli-Q Advantage) and placed on a shaker at 4°C for 2 h. Electrolyte leakage was calculated as the ratio before and after leaves were boiled via a conductivity meter (DDSJ-308, Leici, China; Rohde et al., 2004).
The free proline was determined according to the following method. Briefly, 0.25 g of the plant material was weighed and combined with 1.75 ml 3% sulfosalicylic acid in a test tube and incubated in a boiling water bath for 10 min before centrifugation at 5,000 rpm for 10 min. Next, 0.5 ml of water, 0.5 ml of glacial acetic acid, and 1 ml of 2.5% acid trione were added to 0.5 ml of the supernatant, and the solution was developed for 30 min in boiling water. The absorbance at the wavelength at 520 nm was measured after cooling. According to the standard curve, the average content of free proline in each gram dry weight sample was calculated (Levy, 1980;Torrecillas et al., 1984).
After freezing treatment, the seedlings were transferred to 4°C for 12 h in dark conditions, put into normal condition for recovery for 3-5 days, and the number of seedlings that generated new leaves was counted as the survival rate (Ding et al., 2019). Membrane oxidation in cold stress was assessed by measuring the MDA level .

Whole-Genome Bisulfite Sequencing and Analysis
Total DNA was extracted using the QIAamp Fast DNA Tissue Kit (Qiagen, Dusseldorf, Germany) following the procedure given by the manufacturer. DNA samples were fragmented using sonication and subjected to bisulfite conversion and second-generation sequencing.
The Accel-NGS Methyl-Seq DNA Library Kit (Swift, MI, United States) was utilized to attach adapters to single-stranded DNA fragments for library construction. Paired-end, 2 × 150 bp sequencing was performed at the Hangzhou Lianchuan Biotechnology Center using an Illumina Hiseq 4,000 platform. Library construction, sequencing, and bioinformatics analysis are described in the SI Appendix, Supplementary Information Text 1.

RNA Analysis
Total RNA was extracted with a Trizol reagent (Invitrogen, CA, United States) following the procedure provided by the manufacturer. Next, 1 μg of RNA and oligo dT primers were used to synthesize cDNA in a 20 μl reaction to create the final cDNA library using an mRNA sequence sample preparation kit (Illumina, San Diego, CA, United States). The average insert size for the paired-end libraries was 300 bp (±50 bp). Paired-end sequencing was done on an Illumina Hiseq4000 (LC Sciences, United States) following the vendor's recommended protocol. Read mapping, transcript abundance estimation, and differential expression quantitation were performed as described in the SI Appendix, Supplementary Information Text 1.

Measurement of Rutin and Lysine Content
Fresh seedlings (1 g) of Dingku1 were frozen and ground in liquid nitrogen. The rutin content was analyzed by highperformance liquid chromatography (HPLC) from triplicate independent extractions as described previously . Briefly, dried seedlings (100 mg) were ground, then mixed with 500 μl Na-S™ buffer (2% sodium citrate, 1% HCl, and 0.1% benzoic acid; Beckman, United States) for 30 min in a mixer, and extracted for 10 min via ultrasonication. The free lysine content was analyzed using an A200 Amino Acid Analyzer (Aminosys, Germany), from triplicate independent extractions as described previously (Yang et al., 2016).

Evaluation of Traits of Tartary Buckwheat Varieties by Field-Scale Experiments
To screen for high-quality Tartary buckwheat resources adapted to the high altitude and cold climate of the Tibetan Plateau, dozens of Tartary buckwheat varieties were subjected to a longterm adaptive assessment at the Haidong comprehensive test station (Qinghai Academy of Animal Science and Veterinary Medicine) using advanced field real-time detection equipment and unmanned aerial vehicle (UAV) detection managements. Through long-term observation, three varieties were focused on for further field experimental tracking and indoor experiments: Heifeng1, Diku2, and Dingku1. The experimental site was established in Huangzhong, Qinghai Province (approximately 36°28'N, 101°37'E; Figure 1A), with an average altitude of 2,620 m, an annual average temperature of 3.7°C, and a growing season average temperature of 14°C ( Figure 1B). Plants were observed extensively, such as flower branch and seeds ( Figure 1C). Thereinto, Tartary buckwheat cv. Dingku1 exhibited traits of late flowering, higher germination rate, and higher content of total flavonoids and starch in seeds (Supplementary Figures S1, S2).

Evaluation of Frost Resistance in Tartary Buckwheat Varieties by Laboratory Bench-Scale Experiments
To confirm the findings of the field data, a large-scale and accurate phenotypic analysis of three plant varieties -Heifeng1, Diku2, and Dingku1 -was carried out under laboratory control conditions ( Figure 1B; Supplementary Figure S2). The morphological characteristics of seeds were observed within the 3 days before and after germination. The results showed that there was no significant difference in the morphology of the seeds from the three varieties during germination, but the seeds from Dingku1 were of slightly smaller volume than those of the other two varieties before and after germination, and the compact density was relatively larger (Figure 1B;  Supplementary Figures S2a, S3). During germination, the seed temperature and seed water content of Dingku1 were higher compared with seeds from the other two varieties (Supplementary Figures S1b,c, S4). An experiment of cold treatment (4°C, 24 h) followed by 1 day of recovery was carried out, and the frost resistance was assessed by observing, in detail, the morphology of seedlings with three leaves by optical imaging and image analysis technology. Compared with the other two varieties, the seedlings of Dingku1 were shorter, narrower in width, and more compact (Figure 2A; Supplementary Figure S5). Dingku1, however, did not differ significantly between groups and within groups in anthocyanin and chlorophyll test, and differences were observed in the control group only in nitrogen balance index test ( Figure 2B). The temperature and water content of seedlings were measured before and after cold treatment in three varieties. In contrast, Dingku1 had higher water content and plant temperature after cold treatment (Supplementary Figure S6).
Next, experiments of cold acclimation (memory) and freezing stress (shock) were designed to simulate field weather. Through physiological and biochemical indexes, it was determined that Dingku1 exhibited characteristics for stronger frost resistance (Figure 3). Compared to the other varieties, Dingku1 exhibited a lower ion leakage rate (Figure 3A), higher free proline content (Figure 3B), higher survival rate (Figure 3C), and lower malondialdehyde (MDA) content ( Figure 3D) during freezing stress. Therefore, Dingku1 was selected for additional study on the cold environmental adaptation mechanism.
In addition, the memory group (cold priming) with cold acclimation exhibited better performance (Figure 3), higher free proline content (Figure 3B), and higher survival rate ( Figure 3C) after being subjected to the second freezing stress, compared to the control plants grown at 16-21°C or the plants that experienced freezing shock without prior cold acclimation.

Tartary Buckwheat cv. Dingku1 DNA Methylomes During Cold Stress
DNA methylation is closely related to plant frost resistance (Song et al., 2015). To characterize the Tartary buckwheat Dingku1 methylome during cold stress, a whole-genome bisulfite sequencing was performed and single-base resolution maps of DNA methylation were generated for three test groups: the control group (normal temperature), cold memory group (priming: 4°C for 6 h, followed by 21°C for 18 h, repeated for 4 days, then 0°C for 6 h), and cold shock (0°C for 6 h directly, without priming).
Tartary buckwheat is the only Polygonaceae plant whose methylome has been reported to date. This analysis revealed the details of methylated base sites (Figure 4; Supplementary Table S2). The results show different 5-methylcytosine distributions in different regions for the control and the other two cold treatment groups, with most of the regions with decreasing in CpG methylation mapping to the VIII chromosome (Figures 4A-C). There are several functional genes with hypomethyl-modifications, such as the ubiquitin-activating enzyme E1, helicases, and hydrolases, with potential roles in the cold shock response that are regulated by DNA methylation (Supplementary Table S4). There was different methylation of CHH for all three treatments, with slightly higher methylation densities of CHH in the cold memory group compared to that in the other two treatments, as well as significantly lower methylation density in the cold shock group than that in the other treatments (Figures 4A-C).

A B
FIGURE 2 | Phenotypic observation and physiological test in cold experiment of 3-week-old seedlings. (A) Morphological observation before and after 4°C treatment. DW, Low temperature treatment groups from the top surface testing; DW-side, Low temperature treatment groups from the side face testing; NC, nonspecific control. (B) Physiological test between the control and cold treatments. The every mean value was from more than 30 independent plant measurements, and error bars indicated ± SD. Analysis was performed with one-way ANOVA followed by Tukey-Kramer post hoc analysis. * p < 0.05.
Frontiers in Plant Science | www.frontiersin.org 6 December 2020 | Volume 11 | Article 608540 An analysis of the distribution of DNA methylation in Tartary buckwheat showed the enrichment of DNA methylation in the proximal promoter region (5' terminal of TSS) and downstream region (3' terminal of TES; Figure 4D), indicating increased DNA methylation in the 5' and 3' flanking regions of genes. In addition, genes were characterized by a high enrichment of CG methylation and the moderate enrichment of CHG methylation and CHH methylation (Figure 4D), suggesting that CG was the primary type of cytosine methylation in Tartary buckwheat gene-transcription domains. Little difference was observed in the distribution of methylation among the three test groups (Supplementary Figure S7), and similar trends in methylation levels were observed for the three groups. In the control group, the total DNA methylation level was 24.62%, with 78% CG methylation, 40.28% CHG methylation, and 11.1% CHH methylation (Supplementary Table S3).

DNA Methylation in Dingku1 Decreases During Cold Stress
As shown in Supplementary Figure S8, the principal component analysis (PCA) revealed significant differences among the three groups. The average DNA methylation levels were calculated and a decrease in the global cytosine methylation level from the control group to the cold treatment group (from 23 to 24.5%) was found. Differentially methylated regions (DMRs) were identified and it was found that the change of DNA methylation was mostly due to a decrease in intergenic methylation ( Figure 5A).
To characterize the change of DNA methylation for different treatments, a method based on Fisher's exact test was used to identify DMRs between methylomes (Ausin et al., 2012).
The methylome of the cold memory group or the cold shock group (M/C and S/C, respectively) was compared with that of the control group. As shown in Figure 5, both hyper-DMRs and hypo-DMRs were detected in similar proportion in M/C and S/C (Figure 5B), but there were more hypo-DMRs in all comparisons, with 526,103 intergenic hypo-DMRs identified in the M/C comparison and 664,615 in the S/C comparison, whereas only 267,055 hyper-DMRs were observed for M/C and 272,784 for S/C. These results suggest a global decrease in DNA methylation in cold treatments, especially cold shock. Hypo-methylation was mainly observed in intergenic contexts, although decreased methylation was also observed in the CGI, promoter, intron, and exon regions to a lesser extent, suggesting that they have an important regulatory function of intergenic regions in the Tartary buckwheat genome.

Correlation Between DNA Methylation and Gene Expression in Cold Stress of Dingku1
Transcriptome profiles for Dingku1 were generated by RNA-seq analysis -with three biological replicates for seedlings under the same treatments used to assay methylation -to investigate whether the observed decrease in DNA methylation during a cold response was associated with changes in gene expression. A total of 31,391 genes and 34,067 transcripts were included in the analysis, with more upregulated genes and transcripts than downregulated genes and transcripts in the two treatment groups relative to the control group (Supplementary Figure S9). Furthermore, 24,931 differentially expressed genes (DEGs) were identified, of which 1,893 and 1748 were significantly differentially expressed in the cold memory group and the cold shock group, respectively, relative to the control group. Considering the association between DNA methylation and gene expression, hyper-and hypo-DMR-associated genes were analyzed during cold treatments. There were 377,514 and 382,440 hyper-DMR in the cold memory and the cold shock groups, respectively, which were significantly lower than the numbers of hypo-DMRs for the two groups. There were 742,756 and 889,778 hypo-DMRs in cold memory and cold shock, respectively.
Statistical analysis revealed the potential coupling of the changes in DNA methylation and gene expression. A total of 21,335 genes were identified as DMR-associated genes in the buckwheat genome (FPKM >1). Among these genes, 336 upregulated and 217 downregulated DEGs were identified as hypo-DMR-associated genes, and 143 downregulated and 189 upregulated DEGs were identified as hyper-DMR-associated genes in the cold memory compared to the control. Additionally, 341 upregulated and 205 downregulated DEGs were found to be hypo-DMR-associated genes, and 125 downregulated and 198 upregulated DEGs were determined to be hyper-DMRassociated genes in the cold shock vs. the control (Supplementary Table S5). As was found previously for tomato and orange fruit ripening (Lang et al., 2017;Huang et al., 2019a), DNA hypermethylation is associated with gene activation in Tartary buckwheat cold response. The results of the present study suggest that DNA methylation may play a positive role in regulating the Tartary buckwheat response when exposed to external stresses. In the cold treatments relative to the control (M/C and S/C), there were four clusters: up-up (upregulated DEGs with hyper-DMR), up-down (downregulated DEGs with hyper-DMR), down-up (upregulated DEGs with hypo-DMR), and down-down (downregulated DEGs with hypo-DMR).
Gene Ontology (GO) analysis was performed to understand the potential role of DNA methylation in Tartary buckwheat cold tolerance. With an overall trend of hypomethylation during cold treatment, the analysis was focused on down-up genes (DNA methylation decreased and the expression increased) in M/C and S/C, which were annotated as members of 203 and 199 terms, and the significantly enriched GO terms were 29 and 38 (p < 0.05), respectively. The most significantly enriched genes in M/C are related to integral components of membranes, suggesting that membrane structure is related to the formation of frozen memory, but these genes were not significantly expressed during cold shock (Supplementary Figure S10a). GO analysis revealed 30 ATP binding-related genes enriched in the S/C down-up cluster, suggesting that the activation of these genes occurs by cold shock induction of DNA hypomethylation (Supplementary Figure S10b). ATP binding is an important physiological activity in plants. For example, ATP binding to Cryptochrome2 (cry2) and some other plant cryptochromes promotes the activation of light receptors and increases stress resistance (Eckel et al., 2018). The same enrichment of GO terms was observed in different freezing treatments (richness factor > 0.8), including genes related to protein phosphorylation, DNA repair, endosome organization, endoplasmic reticulum organization, sulfuric ester hydrolase activity, D-alanine ligase activity, and mitotic cell cycle, suggesting that these tissue and cell activities were involved in the cold response and tolerance of plants. There were also obvious differences in the enrichment of GO terms for the different freezing treatments. For example, in the upregulated DEGs with hypo-DMR (down-up) for cold memory treatment relative to the control (M/C), dephospho-CoA kinase activity, voltage-gated calcium channel activity, and cell adhesion were significantly enriched (Supplementary Figure S10a).
For cold shock treatment relative to the control (S/C), glycolipid biosynthetic process, regulation of cell shape, and cell division were specifically enriched (Supplementary Figure S10b). The enriched GO terms were compared for the two different freezing treatments (M/S). The results showed that the number of enriched genes was more different (>80%) for the metabolic process, oxidation-reduction process, and ribonucleoprotein complex, as well as the ribosome, nucleolus, intracellular, and protein binding in M/S (Supplementary Figure S10c). Differential enrichment was observed for some enzymes (including racemase, epimerase, and catechol oxidase) and some metabolic processes (including glycerol metabolism, cellular carbohydrate metabolism, and malate metabolism; Supplementary Figure S10d), indicating differences in DNA replication, RNA transcription, and the tricarboxylic acid cycle due to different freezing treatments inducing different response mechanisms by demethylation regulation.
The effect of cold memory on the adaptability of Tartary buckwheat to a cold environment wanted to be better understood, and the KEGG analysis revealed the enrichment of genes in plant metabolic pathways in comparison of different treatment groups (Supplementary Figures S11a-c). Glycosylphosphatidylinositol (GPI)-anchor biosynthesis and D-Alanine metabolism showed a significant increase in the cold memory group (Supplementary Figure S11a). Plant GPI-anchored protein controls the cellulose content of the cell wall and guides the orientation for cell expansion (Cao et al., 2012). Frozen memory appears to activate the expression of this kind of protein to improve the freezing resistance of plants. Genes related to D-glutamine and D-glutamate metabolism, pyruvate metabolism, methane metabolism, biotin metabolism, and lipoic acid metabolism were significantly increased in the M/S cluster (Supplementary Figure S11c); however, little is known about the role of these pathways in countering freeze stress. There was an obvious difference between plants subjected to the slow accumulation of freezing signals and the rapid response when subjected to freezing. A two-dimensional scatter diagram of correlation analysis was constructed to illustrate the overall relationship between differentially expressed genes and DNA methylation (Figure 6A; Supplementary Table S7). As shown in the Venn diagram, there is an overlap of the DEGs identified in the different treatments ( Figure 6B) Table S6). The results showed shared and distinct aspects of the two cold coping strategies. The significantly different KEGG pathways were combined with metabolomics analysis (Figure 6C). The important genes in the pathways were demonstrated (Figure 7). There was a drop in methylation levels and an increase in the transcriptional level of copper amine oxidase (CuAO), which is a hypo-DMR-associated gene in cold treatment. This gene is involved in multiple pathways, such as the isoquinoline alkaloid biosynthesis pathway in M/C, as well as in the tropane, piperidine, and pyridine alkaloid biosynthesis pathway, and tyrosine metabolism pathway in S/C. RNA polymerase II's largest subunit (RPB1) and dehydrogenase E1 component (DHE1) showed differences in DNA methylation enrichment and transcriptional level in different cold treatments (M/S). The analysis of the present study suggests that DNA methylation can mediate the Tartary buckwheat cold response through the regulation of cold-induced genes.

Change in DNA Methylation Level Affects the Expression of Different Genes and Metabolites
The correlation between promoter methylation levels and transcript levels of genes related to lysine, pyrimidine, alkaloid, and flavonoid metabolism were further examined. Promoter methylation levels and gene expression levels were determined for individual genes for the two cold treatments. For most genes, the expression levels and promoter methylation levels  exhibited a negative correlation ( Figure 8A). For isoquinoline alkaloid biosynthesis, pyrimidine metabolism, and lysine degradation, the expression levels of FtCuAO, FtRPB1, and FtDHE1 showed increased trends, while their promoter cytosine methylation levels were decreased in cold treatments ( Figure 8A). The transcript abundances were measured for three selected unigenes by qRT-PCR to validate transcriptome data. Consistent cold responses were detected between the qRT-PCR analysis and RNA-Seq data ( Figure 8A). CuAO acts in the isoquinoline alkaloid biosynthesis pathway, relating to H 2 O 2 production with polyamine catabolism (Rea et al., 2004) and ABA-induced stomatal closure (Fraudentali et al., 2019). CuAO catabolizes polyamines and is associated with stress responses (Rea et al., 2004;Fraudentali et al., 2019). FtCuAO showed a significant expression increase in the cold shock group by qRT-PCR, which was consistent with the observed response detected by RNA-Seq. The FtRPB1 gene is involved in purine and pyrimidine metabolism and exhibited a significant increase by qRT-PCR in the cold shock group. FtDHE1 participates in lysine degradation and showed a significant decrease in the memory group but a significant increase in the cold shock group based on the qRT-PCR analysis. These genes also exhibited demethylation ( Figure 8A).
To further analyze the significance of DNA methylation for metabolic pathways, the DNA methylation inhibitor 5-aza-2'deoxycytidine (aza-dc) was added into the Dingku1 nutrient solution. Rutin and lysine levels were detected in plants grown in this solution. The results showed that aza-dc had little effect on rutin accumulation but had a great influence on lysine level, especially in the cold memory group (Figure 8B), and we further examined the expression of FtDHE1 under aza-dc treatment, the results also showed that FtDHE1 expression has been disturbed (Figure 8C), suggesting a role of DNA methylation in lysine synthesis and metabolism pathway in cold priming.

DISCUSSION
Temperature is a key factor affecting growth and development in plants. However, the temperature fluctuates under natural conditions, both daily and seasonally. Different plants have evolved a variety of mechanisms to sense complex and variable temperature signals and to regulate their growth, development, or behavior to adapt to changes in environmental temperature. Plants use complex mechanisms to grow under natural conditions and perceive changes of ambient temperature and store temperature memories to better adapt and form desired agronomic traits. In this study, Dingku1 (a Tartary buckwheat variety) was identified as having a high flavonoid content and good temperature tolerance (Figures 2, 3). Repeated cold simulation (cold priming) was used together with high-throughput measurement to find patterns of DNA methylation and relationships between genes and metabolites. The results showed a significant involvement of DNA methylation with enrichment changes of genes and metabolites related to isoquinoline alkaloid biosynthesis, lysine degradation, and pyrimidine metabolism (Figures 6-8).
Environmental temperatures affect virtually all aspects of plant growth and developmental processes. Plant cells can create various chromatin states for stress-responsive gene expression that are required for the adaptation to harsh environmental conditions (Zhao et al., 2020). The epigenetic marks are deemed as environment-dependent patterns through calculations using large populations (Alonso et al., 2019). Environmental factors modulate the epigenomic landscape and regulate adaptive responses in plants. Cytosine methylation (5-mC) is an epigenetic mark associated with developmental programs and stress responses and maintains genome stability by preventing mobilization of transposable elements (TE; Law and Jacobsen, 2010). Environmental cold temperatures can induce the transcription of transposons and lead to desirable agronomic traits in plants (Butelli et al., 2012). The 5-mC profile of a locus reflects the balance between pathways allowing accurate maintenance or a change in DNA methylation, with dynamics affected by abiotic and biotic stresses as well as developmental programs to result in phenotypic changes. This is seen by the overall change of the DNA methylation level during ripening of orange (Huang et al., 2019a), tomato (Zhong et al., 2013;Gallusci et al., 2016;Lang et al., 2017), and strawberry fruit (Cheng et al., 2018). DNA methylation may also be involved in the heterotic traits in broccoli . These observations suggest the potential to alter the methylome to alter adaptation in plants and promote breeding (Molinier, 2020).
The results of the present study show a global loss of DNA methylation during cold treatments, with significant changes in chromosome VIII in cold shock ( Figure 4C). No similar results were reported, and this result will continue to be validated in the future. A detailed analysis of regions in S/C showing differential expression was conducted and DEGs were classified by GO and KEGG analysis. Terms involved in the biological process were enriched, such as the "carbohydrate metabolic process, " "transport, " "proteolysis, " "oxidation-reduction process and lipid metabolic process, " "metal ion binding, " and "amino acid metabolism, " suggesting DEGs regulated by methylation are closely related to the metabolic activity (Supplementary Figure S13).
Prior studies have noted the relationship between epigenetic mechanisms and environmental changes, such as sulfur homeostasis by DNA and histone methylation ; salicylic acid metabolism in heterosis regulated by decreased DNA methylation 1 (DDM1; Zhang et al., 2016); indole-3-acetic acid (IAA) metabolism and transport regulated by epigenetic factors (Mateo-Bonmati et al., 2019); and the relative amounts of different forms of acetyl-CoA, which can be altered by environmental and metabolic factors. Together, these indicate that histone acetylation dynamics integrate metabolic activity to regulate plant responses to stress (Hu et al., 2019). Differential KEGG pathways were identified through joint bisulfite sequencing and transcriptome analysis ( Figure 6C). There was an altered expression of key genes regulated by DNA methylation in these metabolic pathways (Figure 7). For example, genes involved in the isoquinoline alkaloid biosynthesis, the purine metabolism pathway, and the tropane, piperidine, and tyrosine metabolism pathway were upregulated. These results suggest that carbohydrate metabolism, carbon-nitrogen budget, and secondary metabolism are constitutively promoted in cold treatment. Genes related to amino acids and pyrimidine accumulation were upregulated in the cold shock treatment, potentially due to the differential accumulation of lysine and thymidine in both cold treatments.
To the knowledge of the authors, this is the first report on the regulation of lysine metabolic pathways by DNA methylation in plants, and the result of the present study shows different regulatory mechanisms for lysine degradation under different cold treatments (cold memory group and cold shock group). Lysine is an essential amino acid, and lysine level represents crop quality. Lysine is the classic target site of epigenetic modifications (histone methylation and acetylation; Zhang, 2008), and high accumulation of free lysine in endosperm induces multiple plant stress responses . Lysine metabolism is involved in the plant stress response in various ways (Yang et al., 2020a), and the results showed that lysine degradation is involved in the plant freezing response, with different lysine accumulation levels under various freezing treatments affected by DNA methylation (Figure 8B). Rutin is a flavonoid substance that is produced from phenylalanine as a precursor. Freezing treatment affects phenylalanine biosynthesis (Supplementary Table S2), with changes in methylation, but there was little difference between the two cold treatments (Figure 8B). To test the significance of DNA methylation for lysine and rutin enrichment, the DNA methylation inhibitor 5-aza-2'-deoxycytidine was applied during cold treatment, and a significant effect of DNA methylation on lysine content was observed. Most of the differentially expressed sequences were mapped to intergenic regions, which was likely caused by incomplete genome annotation. De novo sequencing should be carried out in the future to deeply analyze the genome of Tartary buckwheat varieties with excellent traits.
Plant memory involves multiple physiological, proteomic, transcriptional, and epigenetic changes, with the important role of epigenetic modification in plant memory confirmed by numerous studies (Iwasaki and Paszkowski, 2014;Kinoshita and Seki, 2014;Dean, 2017;Lamke and Baurle, 2017;He and Li, 2018;Turgut-Kara et al., 2020). In this study, the relationship between the effects of repeated environmental low-temperature stimulation (cold priming) and DNA methylation patterns, as well as changes in gene expression and metabolite enrichment after cold memory generation, was investigated, which leads to changes in some agronomic traits in crops (e.g., free lysine content). Future efforts to improve crops should utilize molecular module theory based on multiomics to improve important agronomic traits and increase tolerance to extreme ambient temperatures (Zhang et al., 2019). Here, differentially expressed genes, metabolites, and the possible roles of DNA methylation modification in cold priming in Tartary buckwheat were investigated. These findings provided comprehensive insights for the role of DNA methylation in cold priming (i.e., cold memory) and facilitated the breeding of ideal agronomic traits in Tartary buckwheat varieties.

DATA AVAILABILITY STATEMENT
The sequencing data generated in this study have been deposited in National Center for Biotechnology Information's Gene Expression Omnibus and are accessible through the GEO Series accession nos. GSE138547, GSE138497 (BS-Seq), and GSE138546 (RNA-seq).

AUTHOR CONTRIBUTIONS
YS, ZJ, and LA contributed to conceive, design, and coordinate the experiments. YH, XM, LL, and XJ performed experiments. YS and ZJ analyzed the data. LA participated to the manuscript revision. YS wrote the manuscript. All authors contributed to the article and approved the submitted version.