Effects of 5-azaC on Iridoid Glycoside Accumulation and DNA Methylation in Rehmannia glutinosa

Iridoid glycoside is the important secondary metabolite and the main active component in Rehmannia glutinosa. However, the mechanisms that underlie the regulation of iridoid glycoside biosynthesis remain poorly understood in R. glutinosa. Herein, the analysis of RNA-seq data revealed that 3,394 unigenes related to the biosynthesis of secondary metabolites were identified in R. glutinosa. A total of 357 unigenes were involved in iridoid glycoside synthesis, in which the highly conservative genes, such as DXS, DXR, GPPS, G10H, and 10HGO, in organisms were overexpressed. The analysis of the above genes confirmed that the co-occurrence ratio of DXS, DXR, and GPPS was high in plants. Further, our results showed that under normal and 5-azacytidine (5-azaC) treatment, the expression levels of DXS, DXR, GPPS, G10H, and 10HGO were consistent with the iridoid glycoside accumulation in R. glutinosa, in which the application of the different concentrations of 5-azaC, especially 50 μM 5-azaC, could significantly upregulate the expression of five genes above and iridoid glycoside content. In addition, the changes in the spatiotemporal specificity of degree and levels of DNA methylation were observed in R. glutinosa, in which the hemi-methylation was the main reason for the change in DNA methylation levels. Similar to the changes in 5-methyl cytosine (5mC) content, the DNA demethylation could be induced by 5-azaC and responded in a dose-dependent manner to 15, 50, and 100 μM 5-azaC. Taken together, the expression of iridoid glycoside synthesis gene was upregulated by the demethylation in R. glutinosa, followed by triggering the iridoid glycoside accumulation. These findings not only identify the key genes of iridoid glycoside synthesis from R. glutinosa, but also expand our current knowledge of the function of methylation in iridoid glycoside accumulation.

It is well known that plant secondary metabolism is a quite complex physiological process and is dramatically affected at transcription level (Borges et al., 2017;Kurepin et al., 2017;Sanchita and Sharma, 2018). As a common epigenetic modification, DNA methylation is closely related to gene expression, transposon silencing, and heterochromatin formation (Qian et al., 2014;Bharti et al., 2015;Wang et al., 2021). It is a crucial means to regulate growth and development or responses to environments (Yamamuro et al., 2014;Duan et al., 2018). The level and status of DNA methylation are associated with plant growth and development, environmental changes, and plant species (Candaele et al., 2014;Duan et al., 2018). The abnormal state of DNA methylation can lead to the aberrant growth and development of plant (Yamamuro et al., 2014). Many studies have found that DNA methylation is involved in plant secondary metabolism and affects SMs accumulation (Tan et al., 2017;Wang et al., 2019;Badad et al., 2021). Under the treatment of DNA methylation inhibitor 5-azaC, the content of polysaccharides and alkaloids was increased significantly in Dendrobium (Ni et al., 2014), and the biosynthesis of tanshinones was enhanced in Salvia miltiorrhiza hairy roots (Jiang et al., 2018;Li et al., 2018). The Abbreviations: 10HGO, 8-hydroxygeraniol dehydrogenase; 5-azaC, 5azacytidine; 5mC, 5-methyl cytosine; DMAPP, dimethylallyl diphosphate; DXR, 1-deoxy-D-xylulose-5-phosphate reductase; DXS, 1-deoxyxylulose-5phosphate synthetase; G10H, geraniol 10-hydroxylase; GPP, geranyl diphosphate; GPPS, geranyl pyrophosphate synthase; HPLC, high-performance liquid chromatography; IPP, isoprene diphosphate; MEP, 2-C-methyl-D-erythritol-4phosphoric acid pathway; MSAP, methylation-sensitive amplified polymorphism; qRT-PCR, quantitative real-time PCR; SMs, secondary metabolites. accumulation of ganoderic acid and the upregulated expression of related synthetase enzyme genes were also induced by the applied 5-azaC in Ganoderma lucidum (Lan, 2016). Although the relationship between DNA methylation and SMs accumulation in plants has been explored, the mechanism of DNA methylation regulating plant secondary metabolism is still unclear.
Rehmannia glutinosa, a perennial herb with great medicinal and economic values, belongs to Rehmannia genus of Scrophulariaceae . Iridoid glycoside is one of the main active components in R. glutinosa . However, the regulatory mechanism of iridoid glycoside biosynthesis has not been reported in R. glutinosa. There are few studies on the relationship between DNA methylation and iridoid glycoside accumulation in plants. To reveal the role of DNA methylation modification in iridoid glycoside synthesis in R. glutinosa, the iridoid glycoside synthetase genes were cloned and their expression, iridoid glycoside accumulation, and DNA methylation in R. glutinosa were analyzed. The results not only provide a scientific basis for development and utilization of R. glutinosa, but also establish a foundation for the studies on epigenetic mechanisms underlying plant secondary metabolism.
Roots and leaves of R. glutinosa were collected, respectively, at Stage E (the elongate stage: the root is fleshy and cylindrical in early June), Stage I (the intumescent stage: the root displays expansion in the middle of September), and Stage M (the mature stage: the root is spindle-shaped in early December) and then stored at −80 • C. In addition, 100 collected root tubers of R. glutinosa in each stage were used for the treatment of 5-azacytidine (5-azaC) with three replicates. In the pilot experiment, R. glutinosa seedlings were treated with 15, 50, 100, 150, and 250 µM 5-azaC. According to the effects of 5-azaC on the growth of R. glutinosa seedlings, 15, 50, and 100 µM 5-azaC were selected further experiments. After sterilized for 10 min by 0.1% HgCl 2 , root tubers of R. glutinosa were washed for 3-5 min with sterile water, and then, these root tubers were treated with 15, 50, and 100 µM 5-azaC, respectively.

RNA-Seq Data Analysis
The total RNA in Jinjiu root was extracted using a miRNA Isolation Kit (Ambion). mRNA of each sample was isolated from the total RNA by using beads with oligo (dT) and then was added with a fragmentation buffer to cleave the mRNA into short fragments. The mRNA above was employed as templates for the synthesis of first-strand cDNA using random hexamer primers. Sequencing of cDNA library was performed by Illumina NovaSeq technology. After obtaining the raw sequencing data, the adapter contaminations and trimming nucleotides with low-quality score were removed. After getting high-quality sequencing data by Trinity, reads are fragmented into smaller pieces, known as Kmer (Grabherr et al., 2011). These K-mers are then used as seeds to be extended into contigs and then component basing on contig overlapping. Finally, De Bruijn was applied to recognize transcripts in the components. Then, the randomness of mRNA fragments and the integrity of mRNA were checked by examining the distribution of inserts on unigenes, and the dispersion of inserts was evaluated by counting the length of inserts. The sufficiency of mapped data was checked by generating saturation curve of total gene number and total reads number.
The functions of the unigenes were annotated as NCBI nonredundant protein (NR), Swiss-Prot, COG, KOG, eggNOG 4.5, and KEGG (Buchfink et al., 2015). The KEGG Orthology of unigenes was analyzed by KOBAS (Xie et al., 2011), and GO Orthology of unigenes was analyzed by InterProScan (Jones et al., 2014). After predicting the amino acid sequence of unigenes, the domain was annotated by HMMER and Pfam database.

The Extraction of Nucleic Acid
Genomic DNA was extracted from root and leaf in R. glutinosa (Duan et al., 2018). The precipitated DNA was washed twice with 70% ethanol solution, and the dried DNA was dissolved in double-distilled water; then, 1% RNase was added to DNA solution. The integrity of genomic DNA was detected by 0.8% agarose gel electrophoresis, and the yield and purity of genomic DNA were determined by spectrophotometry.
Total RNA extraction from root and leaf in R. glutinosa was performed as described by Chen et al. (2021). DNase treatment and phenol-chloroform extraction were used to remove DNA. The integrity of total RNA was detected by 1.0% agarose gel. The yield and purity of total RNA were estimated. In addition, the first-strand cDNA was synthesized according to the instruction of TaKaRa kit (TaKaRa, Japan).

Cloning and Analysis of Target Genes
The synthetase genes of iridoid glycoside of R. glutinosa were cloned, such as DXS, DXR, GPPS, G10H, and 10HGO. The fulllength cDNA sequences were obtained by electronic cloning technology as follows: The partial sequences of target genes from transcriptome sequencing were used as probes to search SRA database by BLASTN. SRA sequence of R. glutinosa that matched the probe sequence was compared, spliced, and extended by DNAMAN 7.0, which was continually carried out until no more SRA sequence can be obtained.
According to cDNA sequence obtained by electronic cloning, the corresponding primers were designed. PCR product was detected by agarose gel electrophoresis. Next, the target fragments were cut and sequenced.

Quantitative Real-Time PCR
The quantitative real-time PCR (qRT-PCR) was carried out in LightCycler 96 fluorescence quantitative PCR instrument. TIP41 was used as the internal reference gene, and all the primer sequences were presented in Supplementary Table 1. According to AceQ Qpcr SYBR Green MasterMix kit (Vazyme, Nanjing, China) instructions, the configuration of reaction system is performed. The expression level of target gene was calculated by 2 − Ct .

Determination of Iridoid Glycoside
The fresh roots and leaves from R. glutinosa were dehydrated in the oven at 50 • C to constant weight. The dried samples were then ground into powder. One gram powder was mixed with 10 mL 70% ethanol solution in 100 mL plugged conical bottle, and the mixture was filtered by ultrasonic technique. The filtrate was filtered and shaken in a 50 mL volumetric flask. Then, 2.5 mL filtrate was diluted to 50 mL.
In this experiment, catalpol was used as reference substance. The determination of iridoid glycoside was conducted as follows: 1 mL of test solution and 2 mL of 1 M HCL were mixed in 10 mL plugged test tube. The tube was heated in a water bath at 90 • C for 15 min and then placed at room temperature for 15 min. The mixture was added with 0.5 mL of dinitrophenylhydrazine ethanol solution and 30 mL 70% ethanol solution. Then, this reaction solution was placed at room temperature for 1 h and determined by spectrophotometer at 463 nm.

Determination of 5-Methyl Cytosine
The 5-methyl cytosine (5mC) content in R. glutinosa was detected by high-performance liquid chromatography (HPLC; Li J. Y. et al., 2020). After genomic DNA was hydrolyzed in order with DNase I, nuclease P1, and alkaline phosphatase, this digestion system was centrifuged for 5 min at 12,000 rpm. The supernatant was filtered with 0.45 µM organic microfiltration membrane, and 5mC content was determined by HPLC.
These chromatography conditions performed in this study were as follows: The mobile phase was composed of 50 mM KH 2 PO 4 and 8% methanol (92:8) with 0.5 ml/min flow velocity, and the analytical column was Agilent C18 Zorbax XDB column (4.6 × 150 mm, 5 µm particle size). Furthermore, nonmethylated cytosines (C) and 5mC in genomic DNA could be detected at 285 nm according to the retention time of C and 5mC.

MSAP Amplification
The methylation-sensitive amplified polymorphism (MSAP) amplification was performed as described by Duan et al. (2020). Genomic DNA was digested with EcoRI/Msp? (M) and EcoRI/HpaII (H) and then was ligated with EcoRI adapter and MspI-HpaII adapter for 15 h at 16 • C. Subsequently, the products above were prepared for MSAP amplification. After the MSAP pre-amplification, the product was used as template in MSAP selective amplification.
DNA methylation levels were quantified by MSAP binary data (Duan et al., 2018). The presence or absence of one clear band was, respectively, scored as "1" and "0." On the basis of the presence or absence in H and M, the banding patterns of DNA methylation were divided into three classes (Supplementary Figure 2), in which (1) the presence of band in H and M was considered to be no-methylation (class I), (2) the presence of band only in H was considered to be DNA hemimethylation (class II), and (3) the presence of band only in M was considered to be DNA full methylation (class III). In addition, DNA total methylation level (%) = (II + III)/(I + II + III) × 100, DNA hemi-methylation level (%) = II/(I + II + III) × 100, DNA full methylation level (%) = III/(I + II + III) × 100.

Statistical Analysis
All data were analyzed using the SPSS 17.0 software (SPSS Inc, Chicago, USA) and presented as the means ± SD from three independent experiments in triplicate for each. The statistical analysis of all data was performed by one-way analysis of variance followed by Duncan's multiple comparison test.

Identification and Annotation of Unigenes Related to Iridoid Glycoside Biosynthesis
To identify the genes of iridoid glycoside biosynthesis, RNAsequencing (RNA-seq) technology was applied to analyze the transcriptomic data from R. glutinosa (Jinjiu). For the RNAseq analysis, the adaptor and low quality of raw data were removed. All clean reads were assembled into expressed sequence tag clusters (contigs), and the above contigs were de novo assembled into transcripts using the Trinity in pairedend method, in which yielded a total of 75,698 unigenes. The analysis of KEGG pathway showed that 3,394, 328, and 357 unigenes were involved in the secondary metabolite biosynthesis, terpenoid biosynthesis, and iridoid glycoside biosynthesis, respectively (Supplementary Tables 2-4). A total of 357 unigenes of iridoid glycoside biosynthesis were further annotated by using KOG database, in which 37.33, 31.33, and 7.65% unigenes were clustered to the energy production and conversion, the carbohydrate transport and metabolism, and the secondary metabolites biosynthesis, transport, and catabolism, respectively ( Figure 1A). The cellular component enrichment from GO database confirmed that the above 357 unigenes were located in the integral component of membrane, intracellular membrane-bounded organelle, cytoplasm, obsolete cytoplasmic part, chloroplast, endoplasmic reticulum, plastid, endoplasmic reticulum membrane, etc. (Figure 1B).
Subsequently, topGO lineage was plotted, and the significance of each node was analyzed. A total of 82 GO terms were classified into 13 gradations (Figure 2). A total of seven unigenes were significantly enriched in dimethylallyl diphosphate biosynthetic process; 10 unigenes were obviously gathered in isopentenyl diphosphate biosynthetic process and methylerythritol 4-phosphate pathway; and seven unigenes were involved in L-proline biosynthetic process.

Cloning and Analysis of Iridoid Glycoside Synthetase Genes
To further explore the key genes of iridoid glycoside biosynthesis, five genes, such as DXS, DXR, GPPS, G10H, and 10HGO, were cloned from R. glutinosa (Supplementary Figure 3). This gene structure was analyzed (Supplementary Figures 4, 5,  Supplementary Table 5). The physicochemical properties of amino acids were predicted. DXS, GPPS, and G10H were hydrophilic protein, and DXR and 10HGO were hydrophobic protein. 10HGO was located in cytoplasm, while others were located in chloroplastb (Supplementary Table 6).
FIGURE 2 | TopGO directed acyclic graph of unigenes related to iridoid glycoside synthesis. Boxes represent top10 GO term with the most significant enrichment. Each box (or ellipse) contains a description and a P-value. Arrow represents containment relationship.
at E, I, and M stage, but the opposite trends were observed in G10H. Unlike the I stage, our results revealed that at E and M stage, the expression level of DXS in root was higher than those in leaf. Although the GPPS transcript level in root was higher than those in leaf at E and I stage, the GPPS expression in root was inhibited at M stage.

Effects of 5-azaC Treatment on Accumulation of Iridoid Glycoside
To explore the effects of DNA methylation on the gene expression of DXS, DXR, GPPS, G10H, and 10HGO, the 15, 50, and 100 µM 5-azaC, a DNA methylation inhibitor, were applicated in R. glutinosa at E, I, and M stage, respectively. Unexpectedly, our results showed that at E, I, and M stage, the expressions of DXS, DXR, GPPS, G10H, and 10HGO in root and leaf were upregulated under 15 or 50 µM 5-azaC treatment, compared with the control groups (Figures 5A,B). However, the expression of DXS, DXR, GPPS, G10H, and 10HGO in root was downregulated at M stage under 100 µM 5-azaC treatment (Figure 5A), but the similar trends were not observed in leaf. In addition, the expression level of these genes in both root and leaf reached the highest at M stage under 50 µM 5-azaC treatment (Figures 5A,B). The classification analysis showed that in root, the expression trend between DXS and DXR was similar, and the expression of GPPS and G10H revealed the similar trend ( Figure 5A). Unlike, 10HGO was divided into one branch alone in root. In leaf, DXS showed the expression trend similar to those of GPPS, and the expression trend of 10HGO was similar to those of G10H ( Figure 5B). However, DXR was divided into one branch ( Figure 5B).
Subsequently, the iridoid glycoside accumulation in root and leaf was further determined under normal and 5-azaC treatment. In control group, the iridoid glycoside accumulation in R. glutinosa was increased gradually from E to M stage in root (P < 0.05; Figure 6A). Although there were no significant differences of iridoid glycoside content in leaf between I and M stage, iridoid glycoside content of at both I and M stage was higher than those at E stage ( Figure 6B). From E to M stage, the iridoid glycoside content in root ( Figure 6A) was significantly higher than those in leaf (Figure 6B), especially M stage (nearly 2-fold). Compared with the control group, the iridoid glycoside content in root was significantly higher under 50 µM 5-azaC treatment at E, I, and M stage, especially E stage ( Figure 6A). Unlike the I and M stage, both 15 and 100 µM 5-azaC treatment induced the iridoid glycoside accumulation in root at E stage, compared to control groups ( Figure 6A). Compared with control, the application of 15 µM 5-azaC could induce the accumulation of iridoid glycoside in leaf at I stage ( Figure 6B). Under 50 µM 5-azaC treatment, the iridoid glycoside content in root was higher than those in leaf (nearly 2-fold) at M stage ( Figure 6B). Under 100 µM 5-azaC treatment, the iridoid glycoside content in root was above 3-fold than those in leaf at M stage ( Figure 6B).

Determination of Genomic DNA Methylation
The above results have found that under 5-azaC treatment, the iridoid glycoside accumulation in root was higher than those in leaf at M stage. In order to explore the variation of genome methylation, the methylation ratio was detected by MSAP amplification in R. glutinosa at M stage under 5-azaC treatment ( Figure 7A). The total and full methylation ratios in 15 µM and 50 µM 5-azaC treatment were higher than those in other groups. The hemi-methylation ratio reached the highest under 15 µM 5-azaC treatment, and the difference among other groups is very little. In addition, the total and full methylation ratios were decreased significantly under 100 µM 5-azaC treatment, compared with control group (Figure 7A). To elucidate the changes in genomic DNA methylation, the full and hemi-methylation levels in both root and leaf were determined from E to M stage. According to full and hemimethylation levels, the total methylation levels were calculated. Compared with leaf, the total and full methylation levels of root were lower at I and M stage (Figures 7B,C). As shown in Figure 7B, the total methylation and hemi-methylation levels presented increasing trend and reached the highest at I stage in root, while full methylation level was almost no change. Hemi-methylation levels in root were significantly higher than full methylation levels (Figure 7B), indicating that most CCGG sequences are keeping the state of hemi-methylation. Total methylation levels were also increased in leaf, especially M stage (P < 0.05) (Figure 7C).
Compared with the control, the 5mC accumulation in root and leaf was inhibited by the application of 5-azaC at E, I, and M stage (Figures 7D,E). Under the 5-azaC treatment, the 5mC contents in root and leaf were decreased gradually in dosedependent manner (Figures 7D,E).

Demethylation Occurs Upon the 5-azaC Treatment
To explore the dynamic changes in methylation sites under 5-azaC treatment, the methylation status of genomic DNA could be divided into two categories (B and C), in which the eight band types (B1-B4 and C1-C4) were identified, respectively ( Table 2). The band patterns for the increase in DNA methylation were presented in type B, and it included the four categories: I → III, I → IV, II → III, and II → IV. By contrast, the band patterns for DNA demethylation were revealed by type C, in which four type C were defined as the IV → II, IV → I, III → II, and III → I. These results of Table 2 confirmed that under 15, 50, and 100 µM 5-azaC treatment, the band numbers of methylation in root were raised by 3, 7, and 5, respectively, compared to the control. As listed in Table 2, our results showed that compared to the control, the band numbers of demethylation in root under 15, 50, and 100 µM 5-azaC treatment were increased as 7, 20, and 25, respectively. Similar to the changes in methylation levels in root, our results further revealed that compared to control, the band numbers of methylation in leaf under 15, 50, and 100 µM 5-azaC treatment were increased as 4, 9, and 5, respectively ( Table 2). Meanwhile, the band numbers of demethylation were increased as 9, 21, and 29, respectively (Table 2). However, the total band numbers of methylation and demethylation in leaf under 5-azaC treatment were more than those in root, especially under 100 µM 5-azaC treatment, suggesting that compared with root, the dynamic transformation of methylation was sensitive to 5-azaC treatment in leaf.
Collectively, the band numbers of demethylation in root and leaf under 5-azaC treatment were more than those in  methylation, suggesting that DNA demethylation was induced significantly by the application of 5-azaC. Moreover, the demethylation degree was more significant under 5-azaC treatment in dose-dependent manner. In addition, it was also shown that in root and leaf, the band patterns for the increase in methylation were almost II → IV, and the band patterns for the increase in demethylation were mainly IV → II.

DISCUSSION
R. glutinosa is a perennial herb and is one of the traditional Chinese medicines . It is well known that iridoid glycoside, which belongs terpenoids, has a wide range of pharmacological effects and important economic value (Keeling and Bohlmann, 2006;Bohlmann and Zerbe, 2012). The previous study confirmed that most prenyltransferases involved in terpenoid backbone biosynthesis were the member of FPP/GGPP synthase families, and all prenyltransferases were interacted with biotin carboxylase CAC2 . A reference genome of R. glutinosa was reported using Nanopore technology, Illumina, and Hi-C sequencing, in which the expansion of the UDP-dependent glycosyltransferases and the terpene synthase gene families was demonstrated (Ma et al., 2021). Nevertheless, the regulation of genes involved in iridoid glycoside biosynthesis was rarely reported in R. glutinosa.
In our study, 357 unigenes involved in iridoid glycoside biosynthesis were identified and annotated (Supplementary Table 4). GO enrichment indicated that these unigenes above were located in the ample organelle ( Figure 1B), suggesting that they were widely distributed in cell. In current, the study of secondary metabolism in plants is mainly focused on the key synthetic enzymes, such as DXS, DXR, GPPS, G10H, and 10HGO. Herein, these genes were cloned from R. glutinosa, and the further analysis showed that the intron was only existed in DXS (Supplementary Figure 5). DXS, DXR, GPPS, G10H, and 10HGO were predicted exist in cytoplasm and chloroplast (Supplementary Table 6). The co-occurrence ratios of DXS, DXR, or GPPS were higher than those in G10H and 10HGO in plant, and they were conservative in Viridiplantae (Figure 3, Table 1). Subsequently, the expression level of DXS, DXR, GPPS, G10H, and 10HGO was detected by qRT-PCR in root and leaf, respectively (Figure 4). DXS, DXR, GPPS, and 10HGO in root were higher than those in leaf, and the expression levels of DXS at I stage and GPPS at M stage in leaf were higher than those in root. Analogously, the expression level of MtDXS2 is abundant in root, but low in tissues above-ground in M. truncatula (Zhou W. et al., 2016). Unexpectedly, the expression level of G10H in leaf was above 3-fold change compared with root (Figures 4A,B). Similarly, SmG10H transcripts were much more abundant in the leaf than in either the root or the stem of Swertia mussotii .
DNA methylation could affect plant growth and development (Candaele et al., 2014;Yamamuro et al., 2014;Zhou W. et al., 2016), SMs accumulation, and the expression of synthetase gene associated with SMs (Qi et al., 2016;Wang B. et al., 2018;Wang et al., 2019). Studies have shown that the gene expression of the secondary metabolic biosynthesis pathway was upregulated under 5-azaC (a DNA methylation inhibitor)  treatment (Ni et al., 2014;Jiang et al., 2018;Li et al., 2018;Bacova et al., 2019;Luo et al., 2021). In our research, under 15 and 50 µM 5-azaC treatment, the expression of DXS, DXR, GPPS, G10H, and 10HGO was upregulated in root and leaf, especially 50 µM (Figure 5). In addition, the accumulation of secondary metabolites was also increased under 5-azaC treatment (Jiang et al., 2018;Li et al., 2018;Bacova et al., 2019;Luo et al., 2021). Ample studies revealed that the expression level of DXS, DXR, GPPS, G10H, and 10HGO could affect the accumulation of the corresponding secondary metabolites. For example, the increase in carotenoid and vitamin E content was observed, when DXS was overexpressed (Estevez et al., 2001). By contrast, the content of artemisinin decreased significantly in Artemisia annua with DXR no-expression (Wang C. et al., 2018). Overexpression of geranyl diphosphate synthase small subunit 1 (LcGPPS.SSU1) could enhance the monoterpene production in L. cubeba and tobacco (Zhao et al., 2019). G10H overexpression increased significantly the accumulation of strictosidine, vindoline, and catharanthine (Pan et al., 2012). Taken together, the iridoid glycoside accumulation was detected under 5-azaC treatment. It was found that the iridoid glycoside accumulation was also increased obviously under 15 and 50 µM 5-azaC treatment in root, especially 50 µM 5-azaC ( Figure 6A). Iridoid glycoside content in root was higher than those in leaf under control and 5-azaC treatment group (Figures 6A,B). It was worth noting that under 50 and 100 µM 5-azaC treatment, the iridoid glycoside content in root was above 2-fold than those in leaf at M stage (Figures 6A,B). Similarly, ganoderic acid (GAs) accumulation was promoted by 5-azaC in Ganoderma lucidum, and the key enzyme genes (hmgr, sqs, se, and ls) of GAs biosynthesis expression level were upregulated (Lan, 2016). The SmDMLs (DEMETER-like DNA glycosylases) expression was downregulated by 5-azaC treatment in Salvia miltiorrhiza, in which the tanshinone content was increased significantly (Jiang et al., 2018;Li et al., 2018), suggesting that the tanshinone accumulation was pronounced by 5-azaC induced DNA demethylation. The effect of 5-azaC on secondary metabolites accumulation and genes expression was also confirmed in Chlamydomonas reinhardtii (Bacova et al., 2019), Dendrobium (Ni et al., 2014), and Pogostemon cablin (Luo et al., 2021). In summary, we speculated that the upregulation of DXS, DXR, GPPS, G10H, and 10HGO induced by 5-azaC treatment was involved in the iridoid glycoside accumulation in R. glutinosa. Previous studies confirmed that the gene expression level was increased by applicating 5-azaC, in which the DNA demethylation was induced (Jiang et al., 2018;Li et al., 2018;Zhu, 2020). Our study showed that at M stage of R. glutinosa, the total and full methylation ratios under 15 and 50 µM 5-azaC treatment were higher compared with control, but decreased significantly under 100 µM 5-azaC ( Figure 7A). The hemi-methylation ratio was higher than other groups under 15 µM 5-azaC treatment ( Figure 7A). Compared with leaf, the total and full methylation levels of root were lower at I and M stage (Figures 7B,C). Genome DNA demethylation levels could be presented by 5mC content. In potato, the 5mC levels were reduced by 5-azaC treatment (Law and Suttle, 2005). 5mC content and total methyltransferase activity were decreased, when treated by 5-azaC in peony (Zhang et al., 2020). Our results showed that 5mC contents in root and leaf were decreased significantly under 5-azaC treatment in dose-dependent manner (Figures 7D,E). In addition, the internal reason for the methylation changed was also revealed in this study. The methylation increasing band patterns were almost II → IV, and demethylation was mainly IV → II (Table 2). Moreover, the methylation changes in the leaf were more active than those in root. In addition, 5mC content was decreased under 15 and 50 µM 5-azaC, but full methylation level showed an increasing trend in root (Figures 7A,D). We speculated that CCGG methylation site could be detected by MSAP, but m C m CGG/GGC m C m sites could not be detected. Thus, the increase or decrease in DNA methylation could not be determined due to the limitation of MSAP, such as I → II, II → I, III → IV, or IV → III. These methylation sites might transform into full methylation (band type III) or hemi-methylation sites (band type II), followed by the methylation ratio which was higher than those of control.
At present, the study showed that the change trend of DNA methylation was dynamically related to isoflavone synthesis in Robinia pseudoacacia, and the isoflavone synthesis was blocked by 5-azaC treatment (Xia, 2006). Overexpressing demethylase gene ROSI in transgenic tobacco could promote significantly the gene expression of flavonoid metabolism (Bharti et al., 2015). Moreover, the upregulation of the expression of genes involved in flavonoid synthesis was induced, when the demethylase gene was overexpressed in transgenic hybrid poplars, followed by increasing the flavonoid accumulation in apical meristem, scale bud, and apical bud (Daniel et al., 2017). Ample studies indicated that DNA demethylation was induced by 5-azaC treatment, which increased the expression of key enzyme genes in secondary metabolic pathway, resulting in the secondary metabolite accumulation.
Overall, our results suggested that under 5-azaC treatment, the demethylation of iridoid glycoside synthetase genes in R. glutinosa might activate the expression of these genes, followed by inducing the iridoid glycoside accumulation. However, the terpenoid synthesis pathway is very complex, and the mechanisms of regulation of DNA methylation on iridoid glycoside synthesis were unclear. Thus, the roles of epigenetic regulation in iridoid glycoside synthesis were still investigated in near future.

DATA AVAILABILITY STATEMENT
The datasets used and analyzed in the current study are available from the corresponding author on reasonable request.

AUTHOR CONTRIBUTIONS
TD and HD conceptualized the study. JS, TD, and YL designed the methodology. TD, SS, YW, RY, and XD were involved in formal analysis. TD, SS, and YW wrote-original draft. JS, PC, and HD were involved in writing-review and editing. HD was involved in funding acquisition. All authors have read and agreed to the published version of the manuscript.