*2.6. Correlation between Gene Expression Levels and DNA Methylation*

The DNA methylation of genes is correlated with changes in their expressions. To reveal the potential roles of expression level differences in WT versus RT regulated by DNA methylation, we analyzed the distribution of gene methylation and expression levels within each chromosome and gene functional region. As shown in Figure 3a, highly expressed genes were distributed in regions of low methylcytosine density. This negative correlation between gene expression levels and DNA methylation in general was apparent on scatter plots and heat maps (Figure S15).

To analyze the relationship between levels of expression and methylation within gene bodies (including 2-kb upstream and downstream domains), we divided genes into high, medium, low, and no expression groups [high expression (FPKM ≥ FPKM\_75%); medium expression (FPKM\_25% ≤ FPKM < FPKM\_75%); low expression (1 ≤ FPKM < FPKM\_25%); and no expression

(Figure 6).

(FPKM < 1). FPKM stands for fragments per kilobase of exon per million fragments mapped, and FPKM\_25% and FPKM\_75% refer to values at the boundary of the 25th and 75th percentiles of expression levels, respectively]. And their methylation levels were plotted as a function of location [48]. As shown in Figure 6, the non-expressed genes displayed high methylation levels at CG sites within domains that are 2-kb upstream or downstream of the gene body and at CHG sites in all of the gene regions. Within the gene body and 2-kb downstream regions, non-expressed genes showed high methylation levels at CHH sites, but gene expression was positively correlated with CG methylation within the gene body. Similarly, a positive correlation was detected between gene expression and CHH methylation levels within 2-kb upstream of the TSS. Interestingly, two peaks in CG and CHG methylation levels were observed in the 2-kb upstream and downstream regions of highly expressed genes. *Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 9 of 27 methylation levels were observed in the 2-kb upstream and downstream regions of highly expressed genes.

**Figure 5.** Analysis of differentially expressed genes (DEGs). (**a**) Detection and distribution of genes differentially expressed between white petal tissues (WT) and red petal tissues (RT). Green and red dots indicate upregulated and downregulated DEGs, respectively. (**b**) Cluster dendrogram of DEGs. (**c**) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of upregulated and downregulated DEGs. **Figure 5.** Analysis of differentially expressed genes (DEGs). (**a**) Detection and distribution of genes differentially expressed between white petal tissues (WT) and red petal tissues (RT). Green and red dots indicate upregulated and downregulated DEGs, respectively. (**b**) Cluster dendrogram of DEGs. (**c**) Kyoto Encyclopedia of Genes and Genomes (KEGG) pathway enrichment of upregulated and downregulated DEGs.

Another statistical approach, projecting methylation levels onto a coordinate system of expression levels (*x*-axis) and frequency (*y*-axis), was used to further explore gene expression and DNA methylation (Figure 7). We divided the methylation levels of gene bodies and promoters into five groups according to Xu et al. [48], and found that DNA methylation and gene expression levels were negatively correlated in most genes. In contrast, genes whose promoters lacked mCHH sequences were only weakly expressed (Figure 7); this situation may have been responsible for the positive correlation between gene expression and methylation levels within 2-kb upstream regions

FPKM\_75% = 45.62].

Another statistical approach, projecting methylation levels onto a coordinate system of expression levels (*x*-axis) and frequency (*y*-axis), was used to further explore gene expression and DNA methylation (Figure 7). We divided the methylation levels of gene bodies and promoters into five groups according to Xu et al. [48], and found that DNA methylation and gene expression levels were negatively correlated in most genes. In contrast, genes whose promoters lacked mCHH sequences were only weakly expressed (Figure 7); this situation may have been responsible for the positive correlation between gene expression and methylation levels within 2-kb upstream regions (Figure 6). *Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 10 of 27

**Figure 6.** Distributions of methylation levels within gene bodies and their 2-kb upstream and downstream regions. Genes were classified into four groups according to their expression levels: no expression (pink; FPKM <1); low expression (green; 1≤ FPKM < FPKM\_25%); medium expression (blue; FPKM\_25% ≤ FPKM < FPKM\_75%); and high expression (red; FPKM ≥ FPKM\_75%). FPKM stands for fragments per kilobase of exon per million fragments mapped, and FPKM\_25% and FPKM\_75% refer to values at the boundary of the 25th and 75th percentiles of expression levels, respectively. Next, the different gene regions (gene body and 2-kb upstream and downstream) were divided into 50 bins, and the methylation levels of each were averaged. WT and RT indicate white petal tissues and red petal tissues, respectively. The *x* and *y*-axes represent gene body regions and DNA methylation levels by sequence context (mCG, mCHG, and mCHH), respectively. (**a**,**b**) Distributions of methylation levels within gene bodies and 2-kb upstream and downstream regions in samples WT [(**a**); FPKM\_25% = 4.60, FPKM\_75% = 45.31] and RT [(**b**); FPKM\_25% = 4.91, **Figure 6.** Distributions of methylation levels within gene bodies and their 2-kb upstream and downstream regions. Genes were classified into four groups according to their expression levels: no expression (pink; FPKM <1); low expression (green; 1≤ FPKM < FPKM\_25%); medium expression (blue; FPKM\_25% ≤ FPKM < FPKM\_75%); and high expression (red; FPKM ≥ FPKM\_75%). FPKM stands for fragments per kilobase of exon per million fragments mapped, and FPKM\_25% and FPKM\_75% refer to values at the boundary of the 25th and 75th percentiles of expression levels, respectively. Next, the different gene regions (gene body and 2-kb upstream and downstream) were divided into 50 bins, and the methylation levels of each were averaged. WT and RT indicate white petal tissues and red petal tissues, respectively. The *x* and *y*-axes represent gene body regions and DNA methylation levels by sequence context (mCG, mCHG, and mCHH), respectively. (**a**,**b**) Distributions of methylation levels within gene bodies and 2-kb upstream and downstream regions in samples WT [(**a**); FPKM\_25% = 4.60, FPKM\_75% = 45.31] and RT [(**b**); FPKM\_25% = 4.91, FPKM\_75% = 45.62].

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 11 of 27

**Figure 7.** Frequencies of expressed genes anchored by methylation sites within gene bodies and promoters. For each methylation sequence context (mCG, mCHG and mCHH), expressed genes were divided into five groups according to methylation level: group 1 (red; methylation level < level\_20%); group 2 (yellow-green; level\_20% ≤ methylation level < level\_40%); group 3 (green; level\_40% ≤ methy\_level < level\_60%); group 4 (blue; level\_60% ≤ methy\_level < level\_80%); group 5 (pink; methy\_level ≥ level\_80%). Level\_20%, \_40%, \_60%, and \_80% represent values at the boundaries of the 20th, 40th, 60th, and 80th percentiles of methylation levels. WT and RT indicate white petal tissues and red petal tissues, respectively. The *x* and *y*-axes represent gene expression levels and gene frequencies, respectively. **Figure 7.** Frequencies of expressed genes anchored by methylation sites within gene bodies and promoters. For each methylation sequence context (mCG, mCHG and mCHH), expressed genes were divided into five groups according to methylation level: group 1 (red; methylation level < level\_20%); group 2 (yellow-green; level\_20% ≤ methylation level < level\_40%); group 3 (green; level\_40% ≤ methy\_level < level\_60%); group 4 (blue; level\_60% ≤ methy\_level < level\_80%); group 5 (pink; methy\_level ≥ level\_80%). Level\_20%, \_40%, \_60%, and \_80% represent values at the boundaries of the 20th, 40th, 60th, and 80th percentiles of methylation levels. WT and RT indicate white petal tissues and red petal tissues, respectively. The *x* and *y*-axes represent gene expression levels and gene frequencies, respectively.

#### *2.7. DEGs with Methylation Modification and DMR-Related Gene Expression 2.7. DEGs with Methylation Modification and DMR-Related Gene Expression*

Analysis of DEG DNA methylation levels revealed large differences between WT and RT (Figure 8). In WT samples, high methylation levels were observed within DEGs, including downregulated and upregulated ones (i.e., involving promoters), at CG, CHG, and CHH sites. In contrast, DEGs in RT samples exhibited low methylation levels, with especially low levels observed at CHG and CHH sites within gene body domains. We also observed that DEG expression level differences between WT and RT samples were positively correlated with ratios of mCG/CG (*r* = 0.037, *p*< 0.01), mCHG/CHG (*r* = 0.046, *p*< 0.01), and mCHH/CHH (*r* = 0.037, *p*< 0.01) in gene body regions. However, in 2-kb downstream regions, a negative correlation (*r* = −0.038, *p*< 0.01) was observed between CHG methylation and expression level differences between WT and RT (Figure Analysis of DEG DNA methylation levels revealed large differences between WT and RT (Figure 8). In WT samples, high methylation levels were observed within DEGs, including downregulated and upregulated ones (i.e., involving promoters), at CG, CHG, and CHH sites. In contrast, DEGs in RT samples exhibited low methylation levels, with especially low levels observed at CHG and CHH sites within gene body domains. We also observed that DEG expression level differences between WT and RT samples were positively correlated with ratios of mCG/CG (*r* = 0.037, *p*< 0.01), mCHG/CHG (*r* = 0.046, *p*< 0.01), and mCHH/CHH (*r* = 0.037, *p*< 0.01) in gene body regions. However, in 2-kb downstream regions, a negative correlation (*r* = −0.038, *p*< 0.01) was observed between CHG methylation and expression level differences between WT and RT (Figure S16).

S16). Next, an expression analysis of DMR-related genes (located in both gene body and promoter regions) identified 1154 mCG–DMR-related genes, 852 mCHG–DMR-related genes, and 4282 mCHH–DMR-related genes. Hypomethylated DMR-related genes were mostly associated with highly expressed genes, whereas the hypermethylation of DMR-related genes was associated with decreased expression levels (Figure 9). Methylation levels of DMR-related genes were negatively correlated with their expression levels with one exception: hypomethylation at CG sites within gene body regions was positively correlated with gene expression (Figure 10).

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 12 of 27

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 12 of 27

**Figure 8.** Methylation levels of differentially expressed genes (DEGs) by sequence context (mCG, mCHG, and mCHH) and gene region in white petal tissues (WT) and red petal tissues (RT). The terms "all", "up", and "down" indicate all, upregulated, and downregulated DEGs, respectively. DNA methylation levels within gene bodies and promoters of genes differentially expressed between WT and RT are shown ("WT vs. RT. Genebody" and "WT vs. RT. Promoter", respectively). **Figure 8.** Methylation levels of differentially expressed genes (DEGs) by sequence context (mCG, mCHG, and mCHH) and gene region in white petal tissues (WT) and red petal tissues (RT). The terms "all", "up", and "down" indicate all, upregulated, and downregulated DEGs, respectively. DNA methylation levels within gene bodies and promoters of genes differentially expressed between WT and RT are shown ("WT vs. RT. Genebody" and "WT vs. RT. Promoter", respectively). mCHH–DMR-related genes. Hypomethylated DMR-related genes were mostly associated with highly expressed genes, whereas the hypermethylation of DMR-related genes was associated with decreased expression levels (Figure 9). Methylation levels of DMR-related genes were negatively correlated with their expression levels with one exception: hypomethylation at CG sites within gene body regions was positively correlated with gene expression (Figure 10).

regions) identified 1154 mCG–DMR-related genes, 852 mCHG–DMR-related genes, and 4282

**Figure 9.** Distributions of differential expression of differentially methylated region (DMR)-related genes between white petal tissues (WT) and red petal tissues (RT). CG, CHG, and CHH refer to CG-DMR, CHG-DMR, and CHH-DMR-related genes, respectively. The *x*-axis indicates methylation-level differences, and the *y*-axis indicates gene differential expression levels. Red dots **Figure 9.** Distributions of differential expression of differentially methylated region (DMR)-related genes between white petal tissues (WT) and red petal tissues (RT). CG, CHG, and CHH refer to CG-DMR, CHG-DMR, and CHH-DMR-related genes, respectively. The *x*-axis indicates methylation-level differences, and the *y*-axis indicates gene differential expression levels. Red dots represent hypermethylated (hypomethylated) DMR-related genes associated with downregulated (upregulated) expression. Blue dots represent hypomethylated (hypermethylated) DMR-related genes associated with upregulated (downregulated) expression.

**Figure 9.** Distributions of differential expression of differentially methylated region (DMR)-related genes between white petal tissues (WT) and red petal tissues (RT). CG, CHG, and CHH refer to CG-DMR, CHG-DMR, and CHH-DMR-related genes, respectively. The *x*-axis indicates methylation-level differences, and the *y*-axis indicates gene differential expression levels. Red dots

S17).

represent hypermethylated (hypomethylated) DMR-related genes associated with downregulated

genes associated with upregulated (downregulated) expression.

**Figure 10.** Combined maps of correlations between differentially methylated region (DMR)-related genes and expression levels. Each of the 12 subfigures is divided into four sections: upper left, box plot of DMR-related gene expression; upper right, scatter plot of DMR-related gene expression vs. methylation level; bottom left, comparison and correlation statistics; bottom right, box plot of DMR-related gene methylation level. WT and RT indicate white petal tissues (gray) and red petal tissues (blue), respectively. **Figure 10.** Combined maps of correlations between differentially methylated region (DMR)-related genes and expression levels. Each of the 12 subfigures is divided into four sections: upper left, box plot of DMR-related gene expression; upper right, scatter plot of DMR-related gene expression vs. methylation level; bottom left, comparison and correlation statistics; bottom right, box plot of DMR-related gene methylation level. WT and RT indicate white petal tissues (gray) and red petal tissues (blue), respectively.

#### *2.8. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment of DMR-Related Genes Associated with DEGs 2.8. Gene Ontology (GO) and Kyoto Encyclopedia of Genes and Genomes (KEGG) Enrichment of DMR-Related Genes Associated with DEGs*

To investigate the effect of DNA methylation on genes with unique functions, we performed GO and KEGG enrichment analyses of gene sets generated from an association analysis between DMR-related genes and 506 DEGs. We first identified 285 DMR-related genes associated with gene-body domains of downregulated or upregulated DEGs. Of these 285, 126, and 125 hypermethylated DMR-related genes were associated with upregulated and downregulated DEGs, respectively, whereas 27 and 18 hypomethylated DMR-related genes were associated with upregulated and downregulated DEGs, respectively. In addition, six upregulated and five downregulated DEGs overlapped with both hypermethylated and hypomethylated DMR-related genes (Figure 11). We also discovered that 282 associated DEGs, 61 of which also overlapped with DMR-related genes within gene body regions, were hypermethylated (213) or hypomethylated (83) within promoter domains. A total of 43 and 40 DMR-related genes that were hypomethylated within promoter regions were associated with downregulated and upregulated DEGs, respectively; these numbers were larger than the number of genes hypomethylated within promoter domains (Figure To investigate the effect of DNA methylation on genes with unique functions, we performed GO and KEGG enrichment analyses of gene sets generated from an association analysis between DMR-related genes and 506 DEGs. We first identified 285 DMR-related genes associated with gene-body domains of downregulated or upregulated DEGs. Of these 285, 126, and 125 hypermethylated DMR-related genes were associated with upregulated and downregulated DEGs, respectively, whereas 27 and 18 hypomethylated DMR-related genes were associated with upregulated and downregulated DEGs, respectively. In addition, six upregulated and five downregulated DEGs overlapped with both hypermethylated and hypomethylated DMR-related genes (Figure 11). We also discovered that 282 associated DEGs, 61 of which also overlapped with DMR-related genes within gene body regions, were hypermethylated (213) or hypomethylated (83) within promoter domains. A total of 43 and 40 DMR-related genes that were hypomethylated within promoter regions were associated with downregulated and upregulated DEGs, respectively; these numbers were larger than the number of genes hypomethylated within promoter domains (Figure S17).

*Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 14 of 27

**Figure 11.** Venn diagram and KEGG pathway enrichment of differentially methylated region (DMR)-related differentially expressed genes (DEGs). DEGs are anchored by DMRs within gene body regions. DMR\_Hypergenes and DMR\_Hypogenes indicate hypermethylated and hypomethylated DMR-related genes, respectively. DEG\_downgenes and DEG\_upgenes indicate downregulated and upregulated DEGs, respectively. CG.DMRgenes, CHG.DMRgenes, and CHH.DMRgenes refer to CG-DMR, CHG-DMR, and CHH-DMR-related genes, respectively. **Figure 11.** Venn diagram and KEGG pathway enrichment of differentially methylated region (DMR)-related differentially expressed genes (DEGs). DEGs are anchored by DMRs within gene body regions. DMR\_Hypergenes and DMR\_Hypogenes indicate hypermethylated and hypomethylated DMR-related genes, respectively. DEG\_downgenes and DEG\_upgenes indicate downregulated and upregulated DEGs, respectively. CG.DMRgenes, CHG.DMRgenes, and CHH.DMRgenes refer to CG-DMR, CHG-DMR, and CHH-DMR-related genes, respectively.

GO enrichment analysis was then carried out to functionally classify the DEGs associated with DMR-related genes into biological process, cellular component, and molecular function categories (Figure S18). We also performed a KEGG analysis on the associated DEGs, which revealed the enrichment of the following pathways that play important roles in flower color formation: flavonoid biosynthesis (ko pmum00941), the biosynthesis of secondary metabolites (ko pmum01110), phenylpropanoid biosynthesis (ko pmum00940), and plant hormone signal transduction (ko pmum04075) (Figure 11, Figures S17, S19, and S20). Examination of methylation level modifications of key structural and transcription factor genes differentially expressed between WT and RT indicated that *Pm031359* (*UGT79B6*) was hypermethylated at CHH sites in exon regions, and *Pm027422* (*YUCCA2*) and *Pm008425* (*MYB305*) were hypermethylated at CHH sites in the promoters of downregulated genes. In upregulated DEGs, *Pm008680* (*UFGT*) and *Pm019063* (*CHI*) were hypomethylated at CG and CHH sites within exons and CG sites within promoters. Hypermethylation was also detected in promoter and intron regions of the *Pm013782* (*DFRA*) gene (Supplementary Data 1). GO enrichment analysis was then carried out to functionally classify the DEGs associated with DMR-related genes into biological process, cellular component, and molecular function categories (Figure S18). We also performed a KEGG analysis on the associated DEGs, which revealed the enrichment of the following pathways that play important roles in flower color formation: flavonoid biosynthesis (ko pmum00941), the biosynthesis of secondary metabolites (ko pmum01110), phenylpropanoid biosynthesis (ko pmum00940), and plant hormone signal transduction (ko pmum04075) (Figure 11, Figures S17, S19 and S20). Examination of methylation level modifications of key structural and transcription factor genes differentially expressed between WT and RT indicated that *Pm031359* (*UGT79B6*) was hypermethylated at CHH sites in exon regions, and *Pm027422* (*YUCCA2*) and *Pm008425* (*MYB305*) were hypermethylated at CHH sites in the promoters of downregulated genes. In upregulated DEGs, *Pm008680* (*UFGT*) and *Pm019063* (*CHI*) were hypomethylated at CG and CHH sites within exons and CG sites within promoters. Hypermethylation was also detected in promoter and intron regions of the *Pm013782* (*DFRA*) gene (Supplementary Data 1).

## *2.9. Detection of Methylated TEs Regulating Candidate Gene Expression 2.9. Detection of Methylated TEs Regulating Candidate Gene Expression*

At the global genome level, TE density was heaviest in genomic regions with a high methylation density and lightest in areas of low gene density and low gene expression (Figure 3a). Most long terminal repeat (LTR) TEs were located in 2-kb upstream and downstream regions of At the global genome level, TE density was heaviest in genomic regions with a high methylation density and lightest in areas of low gene density and low gene expression (Figure 3a). Most long

gene body domains, while *Helitron* TEs were mainly scattered throughout the gene body and long interspersed nuclear elements (LINE) were approximately evenly distributed (Figure 12a). Class-II

**3. Discussion** 

terminal repeat (LTR) TEs were located in 2-kb upstream and downstream regions of gene body domains, while *Helitron* TEs were mainly scattered throughout the gene body and long interspersed nuclear elements (LINE) were approximately evenly distributed (Figure 12a). Class-II TEs in intragenic regions were the most heavily methylated of the intergenic and intragenic class-I and class-II TEs. In addition, more methylated sites were found within TEs in the WT sample than in the RT sample. An analysis of the methylation sequence context of TEs revealed that mCG and mCHG accounted for a large proportion of methylated sites, especially within intragenic areas (Figure 12b). *Int. J. Mol. Sci.* **2018**, *19*, x FOR PEER REVIEW 15 of 27 TEs in intragenic regions were the most heavily methylated of the intergenic and intragenic class-I and class-II TEs. In addition, more methylated sites were found within TEs in the WT sample than in the RT sample. An analysis of the methylation sequence context of TEs revealed that mCG and mCHG accounted for a large proportion of methylated sites, especially within intragenic areas (Figure 12b).

**Figure 12.** Methylation of transposable elements (TEs) in the genomes of *Prunus mume*. (**a**) Distribution of TEs density within gene body and 2-kb upstream and downstream domains. (**b**) Distribution of methylated TEs according to sequence context (mCG, mCHG, and mCHH) in intergenic (red) and intragenic (blue) domains. WT and RT indicate white petal tissues and red petal tissues, respectively. Class I and class II refer to class-I and class-II TEs, respectively. (**c**) Genome-wide length percentages of differentially methylated regions (DMRs) within TEs and other different functional regions (leftmost bar) and the methylome-wide percentage distributions of DMRs with different methylation types in different functional regions. **Figure 12.** Methylation of transposable elements (TEs) in the genomes of *Prunus mume*. (**a**) Distribution of TEs density within gene body and 2-kb upstream and downstream domains. (**b**) Distribution of methylated TEs according to sequence context (mCG, mCHG, and mCHH) in intergenic (red) and intragenic (blue) domains. WT and RT indicate white petal tissues and red petal tissues, respectively. Class I and class II refer to class-I and class-II TEs, respectively. (**c**) Genome-wide length percentages of differentially methylated regions (DMRs) within TEs and other different functional regions (leftmost bar) and the methylome-wide percentage distributions of DMRs with different methylation types in different functional regions.

A total of 1833 DMRs (1437 hypermethylated and 396 hypomethylated) comprising 158 mCG, 302 mCHG, and 1373 mCHH sequences were detected within TEs between WT and RT samples. Hypermethylated mCG, hypomethylated mCG, hypermethylated mCHG, hypomethylated mCHG, hypermethylated mCHH, and hypomethylated mCHH within TEs accounted for 8.3%, 4.8%, 11.6%, 9.2%, 12.6% and 10.7%, respectively, of DMRs in genomic functional regions (Figure 12c). Among these DMRs, 197 were located in 196 TEs inserted into 106 DEGs associated with DMR-related genes (Supplementary Data 3). The gene *Pm008680* (*UGFGT*), with three hypomethylated-mCHH *Copia* TE insertions, was upregulated in WT relative to RT, as was *Pm013782* (*DFRA*) carrying three inserted hypermethylated TEs (two hypermethylated-mCHG *L1* and one hypermethylated-mCHG *Copia*  elements). In contrast, the *Pm031359* (*UGT79B6*) gene, harboring a hypermethylated-mCHH *Copia*  TE insertion in its 2-kb upstream region, and the *Pm031359* (*UGT79B6*) gene, with a hypermethylated-mCHG *Helitron* TE insertion in its 2-kb downstream region, were both downregulated (Supplementary Data 1). A total of 1833 DMRs (1437 hypermethylated and 396 hypomethylated) comprising 158 mCG, 302 mCHG, and 1373 mCHH sequences were detected within TEs between WT and RT samples. Hypermethylated mCG, hypomethylated mCG, hypermethylated mCHG, hypomethylated mCHG, hypermethylated mCHH, and hypomethylated mCHH within TEs accounted for 8.3%, 4.8%, 11.6%, 9.2%, 12.6% and 10.7%, respectively, of DMRs in genomic functional regions (Figure 12c). Among these DMRs, 197 were located in 196 TEs inserted into 106 DEGs associated with DMR-related genes (Supplementary Data 3). The gene *Pm008680* (*UGFGT*), with three hypomethylated-mCHH *Copia* TE insertions, was upregulated in WT relative to RT, as was *Pm013782* (*DFRA*) carrying three inserted hypermethylated TEs (two hypermethylated-mCHG *L1* and one hypermethylated-mCHG *Copia* elements). In contrast, the *Pm031359* (*UGT79B6*) gene, harboring a hypermethylated-mCHH *Copia* TE insertion in its 2-kb upstream region, and the *Pm031359* (*UGT79B6*) gene, with a hypermethylated-mCHG *Helitron* TE insertion in its 2-kb downstream region, were both downregulated (Supplementary Data 1).

Since variation in petal color is highly prized in ornamental plant species, *P. mume* trees
