In this experiment, Rehmannia glutinosa “Jinjiu” was used. In early April, root tubers were planted in test field of Henan Normal University, Xinxiang City, Henan, China (N35°18′13.71″, E113°55′15.05″).

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₂, 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.

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 K-mer (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 non-redundant 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.

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).

The synthetase genes of iridoid glycoside of R. glutinosa were cloned, such as DXS, DXR, GPPS, G10H, and 10HGO. The full-length 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.

The following methods were adopted to analyze target genes: ORF of genes was predicted by ORF finder program (http://www.ncbi.nlm.nih.gov/gorf/gorf.html); the amino acid sequence, molecular weight (MW), isoelectric point (PI), and other physicochemical properties were predicted by ExPASy proteomics server (http://www.expasy.ch/tools/protparam.html); the subcellular localization was predicted by Plant-mPLoc (http://www.csbio.sjtu.edu.cn/bioinf/plant-multi/). In addition, the conservative analysis of DXS, DXR, GPPS, G10H, or 10HGO was performed by NCBI CD-search and NCBI BLASTP (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Amino acid sequences were compared, and phylogenetic trees were constructed by MEGA 7.0. Based on STRING, the co-occurrence of DXS, DXR, GPPS, G10H, and 10HGO was analyzed (https://string-db.org/).

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 (Wang et al., 2021).

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.

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₂PO₄ 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, non-methylated cytosines (C) and 5mC in genomic DNA could be detected at 285 nm according to the retention time of C and 5mC.

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 hemi-methylation (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.

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.

To identify the genes of iridoid glycoside biosynthesis, RNA-sequencing (RNA-seq) technology was applied to analyze the transcriptomic data from R. glutinosa (Jinjiu). For the RNA-seq 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 paired-end 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.

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). The phylogenetic trees showed that DXS, GPPS, G10H, and 10HGO from R. glutinosa could be clustered into one branch with Scutellaria barbata, Sesamum indicum, Handroanthus impetiginosus, and Striga asiatica, respectively (Supplementary Figure 6). The homologous genes of DXS, DXR, GPPS, G10H, and 10HGO in plants were searched by BLAST (Supplementary Figure 7, Supplementary Table 7). The homology of DXS, GPPS, G10H, and 10HGO was high with those in Sesamum indicum, and DXR showed high homology with Handroanthus impetiginosus or Osmanthus fragrans.

Based on STRING database, the gene co-occurrence on DXS, DXR, GPPS, G10H, and 10HGO of R. glutinosa was analyzed (Figure 3), suggesting that they were very conservative in evolution and exist in Bacteria, Eukaryota, and Archaea (Figure 3A). Further analysis confirmed that DXS, DXR, GPPS, G10H, and 10HGO of R. glutinosa showed high similarity with their homologs in Viridiplantae (Figure 3B), especially Erythranthe guttata, Aquilegia coerulea, Citrus clementina, Citrus sinensis, Gossypium raimondii, Eucalyptus grandis, Malus domestica, Morus notabilis, and Nelumbo nucifera (Supplementary Table 8). As shown in Table 1, DXS, DXR, GPPS, G10H, and 10HGO were matched, respectively, by many genes in Viridiplantae, whose gene co-occurrence was found in 19 species in plant kingdom, such as Aquilegia coerulea, Citrus clementine, Citrus sinensis, Erythranthe guttata, Eucalyptus grandis, etc. DXS, DXR, GPPS, and 10HGO were matched in seven species in plants, respectively, in which less co-occurrence was observed in G10H. DXS, DXR, and GPPS of R. glutinosa were matched, respectively, in 33 species, such as Aegilops tauschii, Amborella trichopoda, Boechera stricta, Brassica rapa, Camelina sativa, and Zea mays. In addition, the co-occurrence of DXS, DXR, or GPPS in plants was higher (90% or so), and the co-occurrence of G10H and 10HGO in plants was only 30.6 and 38.9%, respectively (Table 1).

Next, the expression of DXS, DXR, GPPS, G10H, and 10HGO was gain insighted in R. glutinosa. As shown in Figures 4A,B, the expression level of five genes above in root and leaf was increased gradually from E to M stage. The expression level of DXR and 10HGO, especially 10HGO, in root was higher than those in leaf 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.

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).

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 hemi-methylation 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 dose-dependent manner (Figures 7D,E).

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.