In this experiment, R. glutinosa Jinjiu and R. glutinosa 85-5 were used. R. glutinosa Jinjiu is a new breed recently bred, its main advantages are high yield, stable yield, high quality, and multiple resistances. R. glutinosa 85-5 also has high yield and high quality, but is not resistant to water stain and the degradation of some characteristics are more serious. The tuberous roots of R. glutinosa were kindly provided by Yongkang Liu and Cuihong Lu in Agricultural Research Institute of Wenxian County, Henan, China.

As described by Duan et al. (2018), tuberous roots of R. glutinosa were grown and cultured in test field of Henan Normal University, Xinxiang City, Henan, China. The type of planting soil was loam, planting density of R. glutinosa was 30×30 cm, and R. glutinosa was managed by conventional field management. Along with the growth and development of R. glutinosa, plant experiment materials were taken and studied, especially in the following three growth stages of R. glutinosa: the elongate stage (the root is fleshy and cylindrical), expansion stage (the root displays expansion), and the mature stage (the root is spindle-shaped). Roots, stems, and leaves of R. glutinosa were taken and stored at −80°C after quick-freezing with liquid nitrogen.

According to the extraction method of phenylic acid (Isaacson et al., 2006) and some improvements, proteins in roots of R. glutinosa Jinjiu were extracted. The concentration of protein extracts was determined by BCA method (Smith et al., 1985), protein extracts were detected by 12% SDS-PAGE electrophoresis, and the SDS-PAGE gel was visualized by CBB stain according to Candiano's protocol (Giovanni et al., 2004).

The shotgun comparative proteomic analysis of protein extracts from roots of R. glutinosa Jinjiu was studied by the iTRAQ 8-plex experiment. In this study, protein reduce, cysteine block, and digest were performed by FASP method (Jacek et al., 2009). iTRAQ labeling of protein sample was made according to the instruction of iTRAQ 8-plex kit (AB Sciex). After 50 μl tryptic digest (100 μg peptide) and iTRAQ reagent, the extracts were incubated for 2 h at room temperature, the labeled samples were pooled and collected, and then were dried in the vacuum freeze dryer for iTRAQ analysis.

The separation of peptides was performed by RPLC using the following method. The dry sample was re-suspended with 100 μl buffer A and the RPLC was employed on the Agilent 1200 HPLC System. The parameter of HPLC column (Agilent) was: Analytical Guard Column 4.6×12.5 mm 5-Micron, Narrow-Bore 2.1×150 mm 5 μm with 215 and 280 nm UV detection. The separation was performed at 0.3 ml/min. Elution of peptides was made with a mobile phase B of 2–90% over 8 min, each segment with 1min interval for 8-52 min was collected into one tube, four or five tubes were mixed with a total of 10 segments, and every segment was dried in a vacuum freeze dryer. Then the freeze-dried polypeptide samples were re-suspended with Nano-RPLC buffer A for RPLC-MS/MS analysis.

The online Nano-RPLC was employed on the Eksigent nanoLC-Ultra™ 2D System (AB SCIEX). The dissolved samples were loaded on C18 nanoLC trap column (100 μm×3 cm, C₁₈, 3μm, 150Å) and washed by Nano-RPLC Buffer A (0.1% FA, 2% ACN) at 2 μl/min for 10 min. An elution gradient of 5-35% acetonitrile (0.1% formic acid) in 70 min gradient was used on an analytical ChromXP C18 column (75 μm × 15 cm, C18, 3 μm 120 Å, ChromXP Eksigent) with spray tip. Data acquisition was performed with a Triple TOF 5600 System (AB SCIEX, USA) fitted with a Nanospray III source (AB SCIEX, USA) and a pulled quartz tip as the emitter (New Objectives, USA). For information dependent acquisition, survey scans were acquired in 250 ms and as many as 35 product ion scans were collected if they exceeded a threshold of 150 counts per second (counts/s) with a 2⁺-5⁺ charge-state. The total cycle time was fixed to 2.5 s, a rolling collision energy setting was applied to all precursor ions for CID, and dynamic exclusion was set for ½ of peak width (18 s).

Against transcriptome data of R. glutinosa, raw data (.wiff, Sciex) was processed by Protein Pilot Software v. 5.0 (AB SCIEX, USA) using the Paragon algorithm (Shilov et al., 2007).

Furthermore, the experimental data from tandem MS was used to match the theory data to obtain results of protein identification, and protein identification was performed with the search option: emphasis on biological modifications.

In this study, the ID of reliable proteins was nonstandard protein ID; it was necessary to extract sequences from the transcriptome data of R. glutinosa (Sun et al., 2012). Based on sequence similarity and the sequence of retrieved transcriptome data of R. glutinosa, biological information functions of reliable proteins were analyzed using Arabidopsis thaliana as background population.

GO analysis was performed by the mainstream database David6.7 (http://david.abcc.ncifcrf.gov/) and QuickGO (http://www.ebi.ac.uk/QuickGO/) to describe GO classification annotation and enrichment analysis of reliable proteins. KEGG Pathway analysis and enrichment analysis results were obtained by mapping protein information to KEGG database. Furthermore, the analysis of PPI is based on string database (http://string.embl.de/) and cytoscape software (http://www.cytoscape.org/).

CTAB method was used to extract genomic DNA from roots and leaves of R. glutinosa and was performed as described by Duan et al. (2018). The yield and purity of genomic DNA were respectively determined with spectrophotometry at 260 nm, and the integrity of genomic DNA was detected by 0.8% agarose gel electrophoresis.

According to the instruction of RNAiso Plus (TaKaRa, Japan), total RNA was extracted from roots, stems, and leaves of R. glutinosa. In this experiment, DNase treatment and phenol-chloroform extraction were used to remove DNase, furthermore total RNA was dissolved in RNase-free dH₂O. As described in Duan et al. (2019), the integrity of total RNA was detected by 1.0% agarose gel electrophoresis and the yield and purity of total RNA was estimated according to the absorbance at 260 and 280 nm by UV spectrophotometer. Furthermore, the first-strand cDNA of total RNA from roots, stems, or leaves of R. glutinosa was synthesized by reverse transcription referring to the instruction of TaKaRa kit (TaKaRa, Japan).

Based on transcriptome sequencing results of R. glutinosa, the full-length cDNA sequences of target genes were obtained by electronic cloning. According to UTR regions of target genes, primers were designed, and their sequences were listed in Supplementary Table 1. In this study, cDNA and genomic DNA of R. glutinosa were respectively used as templates to obtain full-length cDNA or DNA sequence of target genes by PCR amplification. PCR amplification reaction was performed in a volume of 20.0 μl, which was composed of 2.0 μl template, 1.0 μl each primer (10 μM), 10.0 μl 2×Taq plus, and 6.0 μl ddH₂O. PCR amplification conditions were as follow: pre-denaturation for 5 min at 95°C, followed by 35 cycles (denaturation for 30 s at 94°C, annealing for 30 s and extension for 1 min at 72°C), and the final extension for 10 min at 72°C. Furthermore, the annealing temperatures of FPP/GGP synthase encoding genes were between 51 and 60°C in PCR amplification reaction.

Bioinformatics analysis of target genes were performed using the following method: ORF was analyzed with ORF Finder program (http://www.ncbi.nlm.nih.gov/gorf/gorf.html); amino acid sequence, MW, and other physicochemical properties of putative proteins were predicted by ExPASy proteomics server (http://www.expasy.ch/tools/protparam.html); and functional domain analysis of putative proteins were performed with NCBI-CDS (http://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). In addition, the signal peptide and transmembrane region in putative proteins were detected with Signal P 4.1 (http://www.cbs.dtu.dk/services/SignalP/) or TMHMM (http://www.cbs.dtu.dk/services/TMHMM/), similar sequences of target genes were analyzed by BLASTP (https://blast.ncbi.nlm.nih.gov/Blast.cgi), amino acid sequences of target genes were compared with DNAMAN, the phylogenetic tree of target genes was constructed by MEGA, and the conservative analysis of target genes was performed based on multiple alignments of homologs using string genome (https://string-db.org/).

In this experiment, gene expression of target genes in R. glutinosa was detected by qRT-PCR which was conducted by LightCycler 96 Real-time PCR reaction. The cDNA template of qRT-PCR was prepared from total RNA of R. glutinosa by PrimeScript RT reagent Kit with gDNA Eraser (Perfect Real Time: TaKaRa, Japan), and the internal reference gene was TIP41. Primers used in qRT-PCR were listed in Supplementary Table 1, and were synthesized by YingjieJi Trade Co., Ltd., (Shanghai, China).

The total volume of qRT-PCR reaction system was 20 μl, composed of 2.0 μl cDNA template, 1.0 μl each primer (10 μM), 10 μL SYBR^® Green Master Mix, and 6.0 μl ddH₂O. The reaction conditions of qRT-PCR were as follows: predenaturation for 5 min at 95°C, followed by 45 cycles of denaturation for 5 s at 95°C, and annealing for 20 s at 60°C. The dissolution conditions were 95°C for 15 s, 60°C for 1 min, and 65°C for 15 s, with continuous detection of signals. In addition, the relative expression levels of target genes were normalized and analyzed by the comparative Ct (2^−ΔΔct) method (Duan et al., 2019).

In this study, iTRAQ quantitative proteomic analysis was performed to study the relative abundance of proteins in roots of R. glutinosa. Peptides of six or more amino acids in length and with a maximum of two missed cleavages were exclusively considered for the analysis; the original data were obtained according to search results of protein mass spectrometry (ProteomeXchange identifier: PXD025914). After the database search, 8,929 proteins were detected, and 6,752 reliable proteins were obtained according to the screening criteria of reliable protein.

GO enrichment of all proteins were performed. The significant number of BP, CC, and MF was respectively 1704, 360, or 592 (P < 0.05) (Supplementary Figure 1). GO function and enrichment analysis were shown in Figure 1, most proteins were involved in metabolic process or cellular process (80% or so), about 42% proteins were in response to stimulus, about 22% proteins were involved in developmental process, and 4 and 3% proteins in root development or root morphogenesis, respectively (Figure 1A). Furthermore, there were some metabolic processes with higher significance (P < 0.01), such as small molecule metabolic process, nucleotide metabolic process, and carboxylic acid metabolic process, the biological process response to cadmium ion was also significant, and catabolic process had higher significance (Figure 1B).

In addition, almost every cellular compartment was enriched. Compared to cytosol, about 66% proteins were enriched in membrane-bounded organelle, such as nucleus, chloroplast, mitochondrion, vacuole, ribosome, endoplasmic and reticulum (Figure 1A), and cytosol, membrane-bounded organelle, and chloroplast had higher fold enrichment (P < 0.01) (Figure 1C). As shown in Figure 1A, about 58% proteins had catalytic activity, followed by hydrolase activity, transferase activity, or oxidoreductase activity. Furthermore, most proteins belong to ion binding, protein binding, or nucleotide binding, and protein binding or nucleotide binding had higher fold enrichment (P < 0.01) (Figure 1D).

Analysis and enrichment results of KEGG pathway in roots of R. glutinosa were listed in Supplementary File 1, there were 38 KEGG enrichments with significance (P < 0.05) (Supplementary Figure 1), such as metabolic pathways, carbon metabolism, biosynthesis of amino acids, and TCA cycle. The significant top ranking 10 entries (P < 0.01) were shown in Figure 2A, the P value was very high in carbon metabolism or biosynthesis of amino acids, but was very low in arginine biosynthesis. Furthermore, many proteins were enriched in carbon metabolism, biosynthesis of amino acids, metabolic pathway, spliceosome, and protein processing in endoplasmic reticulum, there were significantly more in the metabolic pathway with 795 proteins (Supplementary File 1). However, less proteins were enriched in proteasome, glycolysis, pyruvate metabolism, TCA cycle, and arginine biosynthesis; only 26 proteins were enriched in arginine biosynthesis (Supplementary File 1).

According to string database and Cytoscape software, PPI in roots of R. glutinosa was analyzed (Supplementary File 2), PPI network of proteins in the top ranking 10 KEGG pathways was shown in Figures 2B,C. There were more PPIs in carbon metabolism, metabolic pathway, and protein processing in endoplasmic reticulum, followed by biosynthesis of amino acids and spliceosome, yet PPIs were less in glycolysis, TCA cycle, or arginine biosynthesis, and were not found in proteasome (Figures 2B,C).

During the growth and development of R. glutinosa, some proteins were differentially expressed, and DEPs in KEGG enrichments of Figure 2 were further analyzed (Table 1). In some KEGG pathways, there were more DEPs between elongation stage and expansion stage, such as carbon metabolism, biosynthesis of amino acids, proteasome, pyruvate metabolism, and TCA cycle pathway, especially the pathway of proteasome and TCA cycle with 15.0% DEPs (Table 1). More DEPs were also found between other stages, such as proteasome between elongation stage and maturation stage or the pathway of pyruvate metabolism and arginine biosynthesis between expansion stage and maturation stage. However, as listed in Table 1, the percentage of DEPs was relatively lower in some KEGG pathways during the growth and development of R. glutinosa, especially spliceosome.

As shown in Supplementary File 1, about 457 proteins were enriched to Biosynthesis of secondary metabolites (ath01110). Some KEGG pathways about the active ingredients of R. glutinosa were also found, such as phenylpropanoid, terpene, stilbenoid, steroid, alkaloid, flavonoid, and so on (Table 2). In R. glutinosa, there were five KEGG regulatory pathways of terpene synthesis: terpenoid backbone biosynthesis pathway, ubiquinone and other terpenoid-quinone biosynthesis pathway, monoterpenoid biosynthesis pathway, diterpenoid biosynthesis pathway, and sesquiterpenoid and triterpenoid biosynthesis pathway (Table 2). Furthermore, 27 and 15 proteins, respectively, were enriched in terpenoid backbone biosynthesis or other terpenoid-quinone biosynthesis.

Based on string database, PPI analysis of proteins associated with terpenoid synthesis was performed (Supplementary Table 2). Except for AAT1 and TEGA10, most proteins could interact with each other, especially At5G27450, HMGS, ISPD, HST, IPP1, ISPF, and SPS2 (Figure 3A). Further analysis indicated that these terpenoid synthesis proteins could simultaneously interact with other proteins, in which most proteins interacted with CAC2, followed by GLU1 and MPA1, and there was also more interaction with ALDH10A8, PPC3, and ASP5 (Figure 3B).

During the growth and development of R. glutinosa, most proteins associated with terpenoid synthesis can be found differently expressed (P < 0.05) (Supplementary Table 3 and Figure 4A). Compared with that at elongation stage, some proteins were up-regulated at expansion stage and maturation stage, such as CHLPA, DXR, AAT1, ISPG, ISPE, PEN3, FPS2, and At3g20160, especially at expansion stag (P < 0.05) (Figure 4A). However, the expression of about 26% proteins at elongation stage was significantly higher than that at expansion stage (P < 0.05), in which some proteins were also up-regulated at maturation stage, such as 4CLL7, FTA, TGA10, COQ5, FACE1, At4g27270, and ISPH, and the expression of 4CLL7, FTA, and TGA10 at the maturation stage was significantly higher than that at elongation stage (P < 0.05) (Figure 4A).

As listed in Table 3, these proteins with prenyltransferase activity in R. glutinosa were divided into five families, namely FPP/GGPP synthase family, UPP synthase family, protein prenyltransferase family, phytoene/squalene synthase family and UbiA prenyltransferase family, and were enriched in some KEGG pathways, except for T27I1.12, RGTA1, and RGTB1 which belong to protein prenyltransferase subunit family. Further analysis found that most prenyltransferases in R. glutinosa belong to FPP/GGPP synthase family, such as FPS1, FPS2, At3g32040, At3g20160, GGPPS1, GGR, and SPS3, which were involved in terpenoid backbone biosynthesis and interacted with biotin carboxylase CAC2. GGPPS1 also interacted with Glucose-6-phosphate isomerase PGI1 (Supplementary Table 3). Furthermore, UPP synthase At5g58770 and protein prenyltransferase subunit alpha FTA were also involved in terpenoid backbone biosynthesis of R. glutinosa (Table 3).

In addition, protein expression of prenyltransferase in R. glutinosa (shown in Supplementary Table 3 and Figure 4B) showed that the expression of FPS2 and At3g20160 was higher at expansion stage (P < 0.05). Compared with that at elongation stage, RGTA1 was continuously up-regulated at expansion stage and maturation stage (P < 0.05). FTA, CHLG, and T27I1.12 all exhibited higher expression at maturation stage, but the expressions of GGR and HST continued to decrease at expansion stage and maturation stage, and GGPPS1, At5g58770, and PSY1 were down-regulated at maturation stage (P < 0.05).

In this study, some FPP/GGPP synthase encoding genes were cloned from R. glutinosa and their cDNA sequences were submitted to GenBank (Supplementary Figure 2), such as FPPS (MG770217), FPPS1 (MT680921), GGPS (MG770218), GGPS3 (MT680922), GGPS4 (MT680923), GGPS5 (MW298275), GPPS (MG770219), and GPPS2 (MW656184). Their full-length CDS sequences were, respectively, 1029, 1050, 1119, 987, 1119, 990, 1272, and 1263 bp (Table 4). Furthermore, intron sequence was found in some FPP/GGPP synthase encoding genes, such as FPPS, FPPS1, GGPS5, and GGPS2; their DNA sequences were shown in Supplementary Figure 3. The full-length DNA sequence of FPPS gene was 3572 bp and was composed of 11 exons and 10 introns, the longest intron was 1073 bp, the shortest intron was 86 bp (Supplementary Figure 3a). As shown in Supplementary Figure 3b, DNA sequence of FPPS1 gene was 4767 bp with 10 introns; the longest intron was 2015 bp, the shortest intron was 83 bp. DNA sequence of GPPS2 gene was 3193 bp and was composed of 6 exons and 5 introns, the longest intron was 625 bp and the shortest intron was 81 bp (Supplementary Figure 3d), yet only one intron with 1259 bp was found in GGPS5 gene (Supplementary Figure 3c).

CDD analysis results showed that FPP/GGPP synthase of R. glutinosa is all composed of three functional domains, namely polyprenyl-synt, IspA, and Trans-IPPS-HT, and belongs to Isoprenoid_biosyn_C1 superfamily which share the same isoprenoid synthase fold. The physical and chemical properties of FPP/GGPP synthase in R. glutinosa were listed in Table 5. PI was 5.5–6.1 and MW was 40 kD or so. MW of GGPS3 and GPPS5 was about 35 kD, but MW of GPPS and GPPS2 was 46 kD or so. FPP/GGPP synthases were all hydrophilic; only one transmembrane region was detected in FPPS (at 180th−220th) and the subcellular locations of FPP/GGPP synthases were different, for example, FPPS and FPPS1 were located in cytoplasm, GGPS4 and GPPS were located in mitochondrion, and others were located in chloroplast. The secondary structures of FPP/GGPP synthase were generally composed of four parts, namely alpha helix, extended strand, beta turn, and random coil. The proportion of alpha helix was the most (50–65%), followed by random coils (25–33%); the proportion of beta turn was the least (Table 5).

As can be seen in Figures 5A,B, the similarity was lower among FPP/GGPP synthases in R. glutinosa, but the identity between their domains was about 37.8%. Further analysis showed that the identity of GGPS and GGPS4 was as high as 98.4% and was also higher (81.6%) between GGPS3 and GGPS5, but the identity of GGPS with GGPS3 or GGPS5 was relatively low, at 38% or so. The similarity of FPPS and FPPS1 was also higher (81.0% identity), yet was 31.4% between GPPS and GPPS2. The identity among other FPP/GGPP synthases was also relatively low. In addition, the Phylogenetic tree of FPP/GGPP synthases in R. glutinosa was divided into two branches, FPPS and FPSS1 were clustered into one branch, other synthases were clustered together, in which GGPPS synthases were clustered into one sub-branch, such as GPPS, GPPS3, GPPS4, and GPPS5, while the other sub-branch was composed of GPPS and GPPS2 (Figure 5C).

Based on string genome, gene co-occurrence on FPP/GPP synthase of R. glutinosa was analyzed. It was found that FPP/GPP synthase was conservative in organisms and occurred in Archaea, Eukaryota, and Bacteria, in which FPPS, FPPS1, GGPS4, GGPS5, GPPS, or GPPS2 of R. glutinosa was more conservative, especially in viridiplantae (Figure 6A). Many genes in viridiplantae appeared to match FPP/GPP synthase of R. glutinosa (Table 6), gene co-occurrence for all FPP/GGPP synthases of R. glutinosa were all found in 39 species of plants, such as Arabidopsis thaliana, Brassica rapa, Erythranthe guttata, Glycine max, Gossypium raimondii, Nicotiana tomentosiformis, and Vitis vinifera. 5–7 FPP/GPP synthases were respectively matched in 14 species of plants, in which GGPS and GGPS4 had less co-occurrence. Only 1–2 synthases in 14 species of plants had co-occurrence, and mainly belong to GPPS. In addition, the co-occurrence of FPP/GPP synthase in plants was higher, especially GPPS and GPPS2 (85% co-occurrence), yet the co-occurrence of GGPS, GGPS4 in plants was 57% or so (Table 6).

Compared with other FPP/GPP synthases of R. glutinosa, GGPS and GGPS3 had higher similarity with their relative homologs in Bacteria, the similarity of FPPS or FPPS1 with the relative homologs was also higher in Eukaryota, but was lower in Bacteria and Archaea (Figure 6A). Further analysis showed that all FPP/GGPP synthases of R. glutinosa had higher similarity with their homologs in streptophyta, especially in Erythranthe guttata, followed by Nicotiana sylvestris, Nicotiana tomentosiformis, Vitis vinifera, Beta vulgaris, Solanum tuberosum, Beta vulgaris, and Arabidopsis thaliana (Figure 6B; Supplementary Table 4). In addition, FPPS or FPPS1 of R. glutinosa had higher similarity with their relative homologs in most plants compared to other synthases (Figure 6B).

As shown in Figure 7A, along with the growth and development of R. glutinosa 85-5, the expression of FPPS gene in root, stem, and leaf increased, and was higher at M stage (P < 0.05), especially in leaf. Similarly, the expression of FPPS gene in R. glutinosa Jinjiu also increased continuously and was higher in root and stem at M stage (P < 0.05). Furthermore, the expression profile of FPPS was obviously varied between R. glutinosa 85-5 and R. glutinosa Jinjiu, except for that in leaf at M stage, the expression of FPPS was generally higher in R. glutinosa Jinjiu, and was especially significant in root and stem (Figure 7A).

Along with the growth and development of R. glutinosa, the expression of GGPS gene in roots decreased continuously (P < 0.05), but increased initially and then decreased in stem and leaf, and was higher at I stage (Figure 7B). Compared with that in stem and leaf, the expression of GGPS was higher in root of R. glutinosa 85-5 at E stage, and was also higher in root of R. glutinosa Jinjiu at E–I stage (P < 0.05). Although the expression trend of GGPS gene was similar between R. glutinosa 85-5 and R. glutinosa Jinjiu and was both lower at M stage, the expression of GGPS was higher at E-I stage of R. glutinosa Jinjiu compared with that of R. glutinosa 85-5 (Figure 7B).

As can be seen from Figure 7C, the expression of GPPS gene in R. glutinosa 85-5 first increased and then decreased along with the growth and development of R. glutinosa, and was higher at I stage (P < 0.05). Compared with that in root, the expression of GPPS was significantly higher except in E stage (P < 0.05). Similarly, the expression of GPPS gene in R. glutinosa Jinjiu also first increased and then decreased, reached the highest peak at I stage, and the expression of GPPS gene in leaf was higher (P < 0.05) (Figure 7C). In addition, the expression of GGPS gene was much lower compared with that of FPPS and GPPS, especially at I–M stage of R. glutinosa (P < 0.05).