The present study focused on ‘FZ’, the earliest flowering P. mume cultivar (Yuan et al., 2020) with a relatively long flowering time and pleasant smell, which was planted in Jiufeng International Mei Garden (40°03′53″N, 116°05′49″E, 132 m a.s.l.), Haidian, Beijing. To reflect the bud developmental state of the whole plant, all buds on one to two twigs of each branch from three different sides of each experimental tree were collected as one biological replicate. Three biological replicates (1, 2, and 3) constituted one sample. Samples were frozen in liquid nitrogen immediately after being obtained and stored at −80°C until sequencing. We selected four samples in total, and took the S2/S4/S6/S8 bud development stage as the dominant stage, from January 9^th to March 13^th in 2019 (Figure 1). The dominant stages were determined by the following reasons: some of the key floral components of P.mume began to evaporate/emit in the early stages of flowering and peak when P. mume is in full bloom (S8). Slightly earlier than that, the biosynthesis of eugenol can be extracted before blooming (S6) (Zhang et al., 2019b) and the expression of key enzyme-encoding genes closely related to floral scent biosynthesis (like PmCFATs and PmBEATs) in the budding phase (S4 or S5) (Bao et al., 2019; Zhang et al., 2019b). The sample (S2) on January 9^th 2019 was also selected for sequencing, for considering that the buds of ‘FZ’ were already finished endo-dormancy on December 23^th 2018, and the methylation modification may be earlier than the onset of gene expression (Yuan et al., 2020).

Genomic DNA of all four sample groups was extracted using Plant Genomic DNA Kit DP305 (Tiangen Biotech, Beijing, China) and then used for library construction using Illumina’s standard DNA methylation analysis protocol and the Accel-NGS®Methyl-Seq DNA Library Kit. The sequencing process was performed on the Illumina Hiseq 6000 platform to generate raw reads, which were preprocessed with Trimmomatic v0.36 software and FastQC analysis (Langmead and Salzberg, 2012). Finally, the clean reads were obtained for subsequent analysis.

Total RNA was extracted with the DP422 extraction kit (Tiangen Biotech). RNA quality and quantity were determined using 1% agarose gels and a NanoPhotometer spectrophotometer (Implen, Calabasas, CA, United States), respectively. An RNA Nano 6000 Assay kit of the Bioanalyzer 2100 system (Agilent, Carpinteria, CA, United States) was used to assess RNA integrity. A total of 3 μg of RNA was used as the input material for sequencing library construction. The NEBNext Ultra RNA Library Prep Kit for Illumina (NEB, Ipswich, MA, United States) was used for reverse transcription. The libraries were also sequenced on the Illumina Hiseq 6000 platform to generate 150 bp paired-end reads (clean bases > 9G).

Using the HISAT2 (https://daehwankimlab.github.io/hisat2/) and Bismark (https://github.com/FelixKrueger/Bismark.git) software to aligned the clean reads in transcriptome and methylome (Kim et al., 2015; Krueger and Andrews, 2011; Langmead and Salzberg, 2012; Wang et al., 2014) to the reference genome of P. mume (http://prunusmumegenome.bjfu.edu.cn/) (Zhang et al., 2012), respectivly. The sequencing depth and coverage were estimated according to duplicates that aligned to a unique genomic region (Krueger and Andrews, 2011; Wang et al., 2014). The non-conversion rate (r) of bisulfite treatment was defined as the rate of the number of sequenced cytosines at all cytosine reference positions divided by the number in the lambda genome.

mC sites were detected using Bismark (Krueger and Andrews, 2011) and defined using a binomial test, with thresholds of sequence depth ≥5 and q-value ≤ 0.01 (Lister et al., 2009; Gifford et al., 2013; Habibi et al., 2013). ML of C site was calculated as follows: ML = mC/(mC + umC) (C here indicate one single C site). The “mCG,” “mCHG,” and “mCHH” contexts (methylcytosine occurring at CG, CHG, and CHH regions, respectively) and their densities, as well as their distributions in each chromosome, were analyzed using previous methods (Chan et al., 2005; Krzywinski et al., 2009; Lister et al., 2009; Zhong et al., 2013). Between samples, the global ML and distributions were also tested (Song et al., 2013). DMRs and differentially methylated loci were identified using DSS software (Feng et al., 2014; Wu et al., 2015; Park and Wu, 2016). Among chomosomes, every 200bp were grouped in one bin. In different functional region of genes, each region is evenly divided into 50 bins. In genesbody and its up- and downstream, each region is evenly divided into 20 bins. The average ML (ML_avg) and methylation density (MD) of each bin was calculated as follows: MD = mCX/CX (coverage > 5X, C here indicate all C site within each bin, CX means CG, CHG or CHH) (Wang et al., 2015). M L a v g = ∑ i = 1 n M L i n

Differentially methylated genes were defined as genes with DMRs within 2 kb of their promoter or gene body (from TSS to TES) region. To clarify the possible regulatory effects of methylation in the natural flowering process of P. mume, the differentially expressed genes (DEGs) between two samples were determined using the following criteria: false discovery rate (FDR) < 0.05 and |fold change| ≥0.

The GOseq R soft package was used for Gene Ontology (GO) enrichment analysis of DMR-related genes that were aligned with GO (p-value < 0.05) (Young et al., 2010). Kyoto Encyclopedia of Genes and Genomes (KEGG) analysis was used to explore pathways of DMR-associated genes (Kanehisa et al., 2008). We also used KOBAS software to test DMR-related genes’ statistical enrichment in the KEGG pathways (Mao et al., 2005).

To show the involvement of methylation in the biosynthesis of P. mume’s flower scent, DMGs enriched in eight pathways that are related to floral scent biosynthesis were listed. We then summarized and drew the key pathway of floral scent biosynthesis, phenylpropanoid biosynthesis, and marked the processes involving DMGs with red asterisks. For further information about the methylation involved in regulating the expression of floral scent-related genes, the FPKM (Fragments Per Kilobase of exon model per Million mapped fragments) values of 145 summarized genes in previous reaserch that encoding the key enzymes in the biosynthesis of the dominant floral volatile of P. mume were used to make heat maps using the heatmap option and the normalized FPKM values of genes (row) in TBtools software. The DMGs were marked and classified according to region (promoter, gene body), context (CG, CHG, CHH), and ML change (hyper, hypo) properties on heatmap.

The qRT-PCR was performed on the CFX96 TouchTM Real-Time PCR Detection System (Bio-Rad, Hercules, California, United States) using the following parameters: 95°C for 5 min, 30 cycles of 95°C for 5 s, and 60°C for 30 s, concluding with a melting-curve stage for 10 s at 95°C, 5s at 65°C, and 5 s at 95°C. Each reaction consisted of 0.8 μL 1^st strand cDNA, 10 μL SYBR Premix Ex Taq (Takara, Dalian, Japan), 0.8 μL each of 10 mM primer pairs, and 7.6 μL H₂O. Each sample was assessed in three biological replicates for each sample and normalized using PmPP2A as an internal control. The transcription levels were determined using the 2^-△△Ct method. The correlation coefficients between RNA-seq and qRT-PCR were calculated with the CORREL function in an Excel spreadsheet. The specific primers were designed online by Primer 3 plus (http://www.primer3plus.com/cgi-bin/dev/primer3plus.cgi) and are listed in Supplementary Table S9.

To investigate the modification of the biosynthesis of floral scent by DNA methylation during the natural flowering process of mei, we performed WGBS of P. mume buds at four different flowering stages: S2, S4, S6, and S8 (hereafter referred to as Fm1–Fm4) (Figure 1). Each stage was sequenced with three biological replicates. For each sample, at least 32 million clean reads were produced (Supplementary Table S2). Approximately 60% of the clean reads were mapped to the reference genome using Bismark (Krueger and Andrews, 2011) (Supplementary Table S2). The average genome coverage was 21.67 × (Supplementary Table S3). All of our sequenced methylated regions had ∼21 × average coverage per chromosome and ∼10 × average coverage per cytosine site, covering an average of 83.8% of the genome (Supplementary Figure S1; Supplementary Table S4).

Among the 82,368,087 C (cytosine) sites in the P. mume genome, the CHH context was the most common and the CG context was the least (Figure 2A). In each sample, ∼8% (7.32–8.46%) of the total C sites were methylated (Table 1). The mC sites mainly featured the mCG context (42.44–46.94%), with the mCHG (28.40–30.19%) and mCHH contexts (22.86–29.14%) comprising the rest (Figure 2B). Among the three contexts of C sites, 33.90%–35.72% of CG regions, 15.60%–16.97% of CHG regions, and 2.21–3.26% of CHH regions were methylated (Table 1). Meanwhile, it was shown that chromosome Pm3 was highly methylated (Supplementary Figure S2). In terms of ML, most mCG and mCHG regions had higher ML (80–100%), while most mCHH regions showed lower ML (10–40%) (Figure 2C). The analysis of methylation within chromosomes showed that both high ML_avg and MD areas were enriched in pericentromeric regions, where the gene density is low (Supplementary Figure S3). The results are consistent with findings from the petals of P. mume cv. ‘FT’ (Jiang et al., 2020), apples, and Arabidopsis.

The average values of ML in the mCG, mCHG, and mCHH contexts were 38.13% (36.29–39.32%), 17.25% (16.47–17.96%), and 1.85% (1.69–2.15%), respectively, while the average global ML in flower buds of P. mume was 8.35% during the flowering process (Table 2). Among different gene functional regions, all CX contexts had the highest ML in the TE region (Figure 3A). Among the gene body and 2 kb regions up- and downstream of it, the CG context had higher ML in the gene body, while the CHG and CHH contexts had lower ML in the gene body than in other regions (Figure 3B).

In the search for DMR-related genes, we detected 5492 (3461 hyper- and hypomethylated), 4744 (2796 hyper- and 1948 hypomethylated), and 7677 (6320 hyper- and 1357 hypomethylated) DMRs within three methylome comparisons between F2v1, F3v2, and F4v3, respectively. During the flowering process, the number of hyper-DMRs was greater than that of hypo-DMRs (Figure 4A). As shown in this Figure 4B, the number of DMRs was mainly from the CHH context, which occurs in promoter regions, and the number of hyper-DMRs was larger than that of hypo-DMRs. The result is consistent with the gradual increase of ML in promoter regions in the CHH context during flowering (Figures 3A,B). In contrast, the number of hypo-DMRs in the final stage of flowering was mainly contributed to by hypomethylated CG-DMRs in promoter regions (Figure 4B).

Next, we detected 4005 (promoter, 3178; gene body, 1006), 3515 (promoter, 2765; gene body, 881), and 5306 (promoter, 4123; gene body, 1533) DMGs, which were associated with DMRs, in comparisons of Fm2v1, Fm3v2 and Fm4v3 (“m” here indicate methylome). In brief, the results showed that the genes related to CHH context DMR within promoter region constituted the largest proportion of DMGs among all methylome comparisons (Figure 5), which is consistent with previous results (Supplementary Table S5). We revealed the functions of those DMGs by KEGG enrichment analysis. Remarkably, some crucial DMGs were enriched in eight KEGG pathways related to floral scent biosynthesis. Overall, 97 (66 promoter, 44 gene body) and 17 DMGs were enriched in phenylpropanoid biosynthesis (ko pmum00940) and its upstream pathway phenylalanine, tyrosine, and tryptophan biosynthesis (ko pmum00400), respectively. Moreover, 20, 11, 14, 4 and 45 DMGs were enriched in the pathways of terpenoid backbone biosynthesis (ko pmum00900), monoterpenoid biosynthesis (ko pmum00902), diterpenoid biosynthesis (ko pmum00904), sesquiterpenoid and triterpenoid biosynthesis (ko pmum00909), and phenylalanine metabolism (ko pmum00360). Furthermore, 371 DMGs were enriched in the biosynthesis of secondary metabolites (ko pmum01110) (Supplementary Table S5).

Using the same samples of mei buds, we also performed RNA-seq to investigate the DEGs involved in floral scent biosynthesis during the natural flowering process. For each sample, at least 50 million (average 62 million, 9.32G over 30× coverage) clean reads were produced (Supplementary Table S6). Over 91% of the clean reads were mapped to the reference genome using Bismark (Supplementary Table S6) (Krueger and Andrews, 2011).

At the chromosome level, low ML and high gene expression levels were distributed in the same region on the chromosome. Higher gene density was also observed in this region (Figure 6A). All genes were divided into four groups according to their expression levels from low to high, and their average ML was determined in different gene regions (gene body and 2 kb up- and downstream of it, TSS, and TES). Among all methylated genes (MGs) at each flowering stage, a higher expression level of the gene was associated with a lower methylation level of the TSS region. Furthermore, higher average ML was detected in the promoter regions of the CG and CHG contexts of unexpressed genes (Figure 6B).

The MGs were divided into five groups (groups 1–5) according to their ML from low to high, and the gene distribution frequencies (y-axis) at different expression levels (x-axis) were represented graphically for each group. In the promoter regions of the CG- and CHG-context MGs, a higher ML was associated with more genes being distributed in the low-expression region log2(FPKM+1) <1, and fewer genes being distributed in 2 < log2(FPKM+1) < 8. The CHH-context MGs also had similar distributions in frequency, although genes with the highest ML (group 5) appeared to have more genes in region 2 < log2(FPKM+1) < 8 than group 4 (Figure 6C). The results indicated that ML in the promoter regions was negatively correlated with expression level. In the gene body regions, the MG group with the highest ML was dominated by a low expression level in all CX contexts. From groups 1 to 3, a higher ML was associated with a higher frequency of genes in the expression level of region 2 < log2(FPKM+1) < 8. The results indicated that ML in gene body regions is positively correlated with expression level (Figure 6C).

The comparison of F2v1 in Figure 6D is here taken as an example. In the gene body region, there were 32 upregulated and 43 downregulated genes in hyper-DMG, and 28 upregulated and 41 downregulated genes in hypo-DMG. Meanwhile, in the promoter region, there were 253 upregulated and 268 downregulated genes in hyper-DMG, and 130 upregulated and 107 downregulated genes in hypo-DMG. Similar results were also found in F3v2 and F4v3. Some of the DMGs have both hyper- and hypomethylated regions (Figure 6D). These results showed that, when limiting the analysis to certain genes, the expression levels of many genes do not always conform to having a negative correlation with the ML in the promoter and a positive correlation with the ML in the gene body.

We detected 7706 (upregulated 3848; downregulated, 3858), 5016 (upregulated 2259; downregulated, 2457), and 12,451 (upregulated 5502; downregulated, 6949) DEGs in compairsion Fr2v1, Fr3v2 and Fr4v3, respectively (“r” here indicate transcriptome) (Figure 6D). The results of KEGG enrichment analysis showed that a total of 214 DEGs were involved in seven floral scent biosynthesis pathways. As the key pathway for biosynthesis of the characteristic aroma of P. mume, phenylpropanoid biosynthesis (ko pmum00940) was enriched the most for DEGs (96 DEGs). Meanwhile, its upstream pathway, phenylalanine, tyrosine, and tryptophan biosynthesis (ko pmum00400), featured 32 DEGs. Moreover, 21, 39, 13, 14, and 10 DEGs were involved in phenylalanine metabolism (ko pmum00360), terpenoid backbone biosynthesis (ko pmum00900), monoterpenoid biosynthesis (ko pmum00902), diterpenoid biosynthesis (ko pmum00904), and sesquiterpenoid and triterpenoid biosynthesis (ko pmum00909), respectively (Supplementary Table S7).

Over 90% of the volatile components of P. mume fragrance were shown to be synthesized by the phenylpropane biosynthetic pathway. We also integrated some genes that previously characterized in the biosynthesis of other characteristic fragrance components in P. mume into the phenylpropane biosynthetic pathway, and then marked the processes involving DMGs with red asterisks (Figure 7). The results clearly showed that DMGs were widely involved in almost all steps in the phenylpropanoid biosynthesis process, including not only the dominant and key components, but also the low (FPKM value < lower quartile) or none (FPKM value < 1) part in P. mume fragrance.

We further summarized 145 genes that were related to the key enzyme in the phenylpropane biosynthetic pathway, as identified in our published results of P. mume floral fragrance research (Bao et al., 2019; Bao et al., 2020; Zhang et al., 2019b; Zhao et al., 2017), including 2 PmPALs, 6 PmPAASs, 37 PmCNLs, 11 PmBARs, 81 PmBEATs, 4 PmCFATs, and 4 PmEGSs. According to the heatmap results, among 47 candidate genes that exhibited an expression pattern matching the pattern of floral scent emission, 22 were differentially methylated during the flowering process (Figure 7, Supplementary Table S10). Some of these genes were proven to encode key enzymes involved in the characteristic fragrance of P. mume in previous studies, such as PmCFAT1a/1c (Zhang et al., 2019b), PmBEAT36/37 (Bao et al., 2019), and PmBAR8/9/10 (Bao et al., 2020); or to have the same expression pattern as in other P. mume cultivars, such as PmPAL2, PmPAAS3, and PmCNL1/3/5/6/14/17/20 (Zhao et al., 2017; Bao et al., 2020) (Supplementary Table S8).

Among them, PmCFAT1a, which have been proved to produced eugenol through an enzyme assay (Zhang et al., 2019b), had four hyper-DMR in its genebody region during F1 to F2, involving all three CX contextsthe with ML changes from: 26.19–44.65% (CG), 20.05–39.34% (CHG), 2.43–5.64% and 1.06–3.45% (CHG). The ML of all hyper-DMR exceeded the average ML of the sample (Table 2; Figure 7). Pm011009 and Pm011010 (PmBEAT36 and PmBEAT37) have been reported to significantly positive affect the synthesis of benzyl acetate in the petal cells of P. mume when overexpressed (Bao et al., 2019). Pm011009 had CG context hypo-DMR (90.49–81.02%) in the promoter region while the gene expressing were significantly increased in F3 to F4. Pm011010, which have been reported to inhibited the expression of Pm011009 in protoplasts, had CHH context hyper-DMR (22.41–39.87%) in its promoter region and CG context hypo-DMR (41.82–17.63%) in its genebody region from F3 to F4, although the expressing were also increased (Figure 7). We randomly selected 6 genes to do qRT-PCR to verify the transcriptome results. The qRT-PCR results showed a similar trend to the FPKM value (Supplementary Figure S4), proving that the transcriptome results are reliable.