The R. glutinosa transcript database obtained by RNA-Seq collected over 87,665 unigenes from R. glutinosa leaves and tuberous roots. The coding sequence (CDS) of WRKY genes from A. thaliana was used to determine homology in the R. glutinosa transcript database by basic local alignment (BLASTn). To remove redundancy, sequences were assembled using the SeqMan function of the DNAStar software package and manually adjusted. Only sequences that shared >95% of the matches are considered redundant. The open reading frames (ORFs) were predicted using ORF Finder¹ and translated into amino acid sequences. Finally, to confirm that the obtained sequences were WRKY members, all primary identified non-redundant amino acid sequences of the WRKY members were submitted to the website http://pfam.sanger.ac.uk to prediction of the WRKY structural domain. Only the sequences that shared the WRKY domains were confirmed to be WRKY members. The R. glutinosa WRKY sequences reported in this work have been submitted to GenBank under accession numbers MZ285925-MZ285961.

To identify potential protein motifs in R. glutinosa WRKY, we used the MEME version 4.9.1 tool (http://meme-suite.org/tools/meme; Bailey et al., 2009) with the following parameters: The distribution of motifs was 0 or 1 per sequence; the maximum number of motifs was 10 motifs, the minimum motif width was 6, and the maximum width was 50. In addition, only motifs with e-value≤1e⁻¹⁰ were retained for further analysis. Subsequently, the detected motifs were searched in protein databases using the SMART (http://smart.embl.de/; Letunic et al., 2015) program. The amino acid sequences of Arabidopsis and rice WRKY proteins were downloaded from the NCBI database, and phylogenetic trees were constructed using the neighbor-joining method and bootstrap analysis (1,000 replicates) of MEGA 6.06 (www.megasoftware.net; Tamura et al., 2013).

The total RNA was extracted from hairy roots with elicitor treatment and 12 different tissues of R. glutinosa using TRIzol reagent (Invitrogen) according to the instruction and then treated with RNase-free DNaseI (Invitrogen). The RNA quality and concentration were determined using a Nanodrop™ 2000 spectrophotometer (Thermo Fisher Scientific, United States) and a Bioanalyzer 2100 (Agilent, United States). RNA-seq analysis of SA- and MeJA-treated hairy roots of R. glutinosa were conducted using Illumina HiSeq™ 2000 platform (Project ID: F15FTSCCKF2309). The H₂O₂-treated hairy roots (Project ID: F18FTSCCKF1733) and 12 R. glutinosa tissues (Project ID: F19FTSCCKF2067) were subjected to RNA-Seq analysis using BGISEQ-500 from BGI-Tech (Shenzhen, China). The RNA-Seq data (SRR5438036, SRR5438037, and SRR5438042) from SA-treated R. glutinosa hairy roots have been deposited in NCBI under project number PRJNA382479. RNA-seq data (CRA004677 CRA004688 and CRA004689) from MeJA- and H₂O₂-treated R. glutinosa hairy roots and 12 R. glutinosa tissues have been deposited in Genome Sequence Archive² under project number PRJCA006052, PRJCA006054, and PRJCA006055.

Transcripts from tuberous roots and leaves transcriptomes of R. glutinosa (Project ID: F13FTSCCKF1467) were used as references for read mapping and gene annotation (Wang et al., 2017a). The clean reads from each sample were individually aligned to the reference transcripts of R. glutinosa using the Bowtie2 software (Langmead et al., 2009), and the abundance of gene transcripts was estimated using the RSEM method (Li and Dewey, 2011) and measured as fragments per kilobase of transcript per million fragments sequenced (FPKM; Trapnell et al., 2010). Genes with the FPKM fold change absolute value≥2 and controlling for false discovery rate adjusted value of p<0.001 were designated as differentially expressed genes. The expression profiles of DGEs from different samples were analyzed by hierarchical clustering, and a heat map of expression values was generated using the T-MeV 4.9.0 software (Howe et al., 2011).

Total RNA was extracted using TRIzol reagent (Invitrogen) according to the manufacturer’s instructions and treated with RNase-free DNase I (Invitrogen). cDNA was synthesized from 1μg of total RNA using the PrimeScript™ II 1st Strand cDNA Synthesis Kit (TaKaRa Bio, Dalian). Relative expression levels of genes were analyzed by qRT-PCR using SYBR® Premix Ex Taq™ II (Tli RNaseH Plus; Takara Bio, Dalian) on a Bio-Rad iQ5 Real-Time PCR System (Bio-Rad, United States), as described by Wang et al. (2017a). The RgTIP41 gene (GenBank accession number KT306007) was used as an endogenous control to normalization relative expression levels based on the 2^-ΔΔCt method (Schmittgen and Livak, 2008). Data were analyzed using ANOVA and Student’s t test (p<0.05). The specific primer sequences used in this study are listed in Supplementary Table S1.

The content of acteoside in the hairy roots and different tissues of R. glutinosa was determined according to a high-performance liquid chromatographic (HPLC) method established by the group (Wang et al., 2017a). Total PhG content in the samples was determined using a UV spectrophotometer, referring to the method established in a previous study (Yi et al., 2017).

For overexpression of RgWRKY37 vectors, primers of the complete RgWRKY37 ORF sequences were designed using Primer-BLAST online, software,³ and the sequences are shown in Supplementary Table S1. The RgWRKY37 sequence was amplified using PrimeSTAR® HS DNA Polymerase (Takara, Dalian), and the amplification product was purified. The purified product was subsequently cloned into a pBI121 to construct the RgWRKY37 overexpression vector, which is driven by the CaMV35S promoter with a fragment of the neomycin phosphotransferase II (NPTII) gene, conferring resistance to kanamycin. The recombinant vector was transformed into the Agrobacterium rhizogenes strain MSU440 using the freeze-thaw method. Specific primers (Supplementary Table S1) were designed to amplify the complete RgWRKY37 ORF sequence without the stop codon. The amplified fragment was inserted into the pBWA(V)HS-GLosgfp vector with a CaMV35S promoter and fused to the N-terminal region of the green fluorescent protein (GFP) gene, to generate the CaMV35S:RgWRKY37-GFP recombinant plasmid, and transformed it into the Agrobacterium tumefaciens LBA4404 strain using the freeze-thaw method.

The A. rhizogenes MSU440 with RgWRKY37-OE vector was initiated from glycerol stock and grown overnight in LB liquid medium at 28°C with shaking (180rpm) until OD₆₀₀ was 0.6. Leaves from 25-day-old seedings of R. glutinosa were cut into 0.5×0.5 leaf disks and infected with A. rhizogenes MSU440 liquid inoculation medium for 5min. The infested leaf disks were placed on MS solid medium containing 100μmol/l acetosyringone (AS) and incubated in the dark at 26°C for 2days. After coculture, the leaf explants were transferred to MS solid medium containing 250mg/L timentin and 100μmol/L AS for hairy roots induction. Numerous hairy roots were observed emerging from the wound sites after approximately 2–4weeks. When the hairy roots had grown to a length of 2–3cm, they were separated from the explants and cultured in the dark at 26°C on MS medium with a gradual decrease in timentin at a 20-day interval to obtain bacteria-free culture. Well-grown hairy root was inoculated in flasks containing 50ml MS liquid medium and cultured on a gyratory shaker at 120rpm in the dark (26°C). Total DNA was isolated from R. glutinosa hairy roots and non-transformed (NT) hairy roots cultured in suspension for 50days using a modified CTAB method (Cuc et al., 2008). The DNA was used as template for transgenic identification by PCR with specific primers for pBI121-D and rolB-D (Supplementary Table S1). The amplification system and procedure were identical to that of Wang et al. (2017a).

To construct gene-metabolite regulatory network of RgWRKY37, the RgWRKY37 overexpressed hairy root and wild-type (WT) hairy root cultured for 40days were used to RNA-seq analysis, respectively. The RNA-Seq data (CRA004786) from RgWRKY37 overexpressed R. glutinosa hairy roots have been deposited in Genome Sequence Archive under project number PRJCA006249. Correlations between the WRKY genes, enzyme genes, and acteoside were calculated using the Pearson correlation coefficient using DPS v7.05 software (Tang and Zhang, 2013) based on the co-occurrence principle between mRNA and metabolite levels. The acteoside content and the FPKM values of RgWRKY37 and key enzyme genes in acteoside biosynthesis were used to construct coexpression networks. The expression correlation matrix was generated with Cytoscape software (Shannon et al., 2003) to measure the similarity of expression between pairwise genes. Gene pairs with r>0.60 (positive coexpression) or r<−0.60 (negative coexpression) were considered significantly coexpressed.

Protoplasts were prepared from Nicotiana benthamian mesophyll cell and transformed, according to previously described method (Yoo et al., 2007). The protoplasts were transformed with highly purified RgWRKY37-GFP-containing plasmid DNA. After 24h of transformation, the protoplasts were observed under a laser confocal scanning microscope (Olympus FV10-ASW, Tokyo, Japan), and the protoplasts expressing GFP alone were used as control.

In order to predict the WRKY gene of R. glutinosa, 71 WRKY gene cDNA sequences of Arabidopsis were downloaded from GenBank⁴ as search sequences. A total of 115 unigenes were retrieved in the R. glutinosa transcriptome database by BLASTn. There were 47 sequences with bases between 200 and 500bp and 68 sequences larger than 500bp. Online prediction of the 68 sequences larger than 500bp was performed using the online tool ORF Finder.⁵ The results showed that 44 cDNA sequences had complete ORFs, and 37 WRKY gene cDNA sequences were eventually obtained after removing four redundant sequences (Table 1). It can be seen that the length of the 37 identified WRKY genes of R. glutinosa ranges from 835 to 2,541bp, and the predicted amino acid sequence length ranges from 192 to 712 aa. Analysis using Batch CD (CDD) search software and SMART software confirmed that all 37 sequences with complete ORFs were all cDNA of the WRKY gene.

The predicted amino acid sequence of WRKY in R. glutinosa showed a homologous alignment on NCBI, and it was found that 36 sequences had the highest consistency with the WRKY sequence of Sesamum indicum, a plant of the family Pedaliaceae. RgWRKY17 has the highest sequence identity with the homologous protein of S. miltiorrhiza homologous protein from Labiatae. RgWRKY7 and RgWRKY37 have the highest sequence identity with the homologous protein of Erythranthe guttata, which also belongs to Scrophulariaceae. Pedaliaceae, Labiatae, and Scrophulariaceae all belong to the Sympetalae subclass of Tubiflorae. RgWRKY23 and RgWRKY32 have the highest sequence identity with Pyrus bretschneideri from Rosaceae and Handroanthus impetiginosus from Bignoniaceae, respectively.

To further examine the evolutionary relationships among R. glutinosa, Arabidopsis, and rice WRKY proteins, multiple sequence alignments of the complete amino acid sequences of all WRKY proteins were performed using the MAFFT program. An unrooted phylogenetic tree (Figure 1) was then constructed using the neighbor-joining method of MEGA6.06 (Tamura et al., 2013) with the multiple sequence alignment files. Based on the number of WRKY domains (WDs) and the features of the specific zinc finger motifs, all 37 RgWRKY genes were classified into three main groups, with five subgroups in group II. Twelve R. glutinosa WRKY genes (RgWRKY1 to RgWRKY12) with two WRKY domains belong to group I and have a zinc finger motif of C–X₄–C–X_22-23–H–X₁–H. The other 20 R. glutinosa WRKY proteins with the zinc finger structure of C–X_4-5–C–X₂₃–H–X₁–H were assigned to group II, which comprised 54% of the total number of RgWRKY genes. The 20 RgWRKY genes of group II are unevenly distributed among the five subgroups: group IIa (two: RgWRKY34, RgWRKY35), group IIb (one: RgWRKY30), group IIc (nine: RgWRKY13, RgWRKY14, RgWRKY15, RgWRKY16, RgWRKY17, RgWRKY18, RgWRKY19, RgWRKY23, and RgWRKY27), group IId (six: RgWRKY20, RgWRKY21, RgWRKY22, RgWRKY24, RgWRKY25, and RgWRKY26), and group IIe (two: RgWRKY31, RgWRKY36). In contrast to group I, group II genes have only one WRKY domain. Instead of the C₂H₂ pattern, group III genes contain a C₂HC zinc finger motif (C–X₇–C–X₂₃–H–X₁–C) and five of the 37 RgWRKY genes (RgWRKY28, RgWRKY29, RgWRKY32, RgWRKY33, and RgWRKY37) belong to this group. Detailed information about the classification of the genes and the WRKY domains and the profile of zinc finger motifs can be found in Table 1.

To better understand the conservation and diversification of the R. glutinosa WRKY proteins, the putative motifs of all WRKY proteins were predicted by MEME motif analysis, and 10 distinct motifs were identified. The locations of the WRKY domains, zinc finger binding motifs, and any conserved motifs are shown in Figure 2. Motif 1 plus Motif 2 was found in all 37 WRKY proteins, and Motif 3 plus Motif 8 was found only in the N-terminal of 12 WRKY proteins. SMART analysis revealed that Motif 1 plus Motif 2 and Motif 3 plus Motif 8 are conserved WRKY domains and zinc finger motifs.

The most prominent structural feature of WRKY proteins is the WRKY domain, which has been shown to interact with the W-box (C/T)TGAC(T/C) to activate a large number of defense-related genes (Eulgem et al., 2000). The WRKY domain consists of a highly conserved WRKYGQK heptapeptide stretch at the N-terminus, followed by a zinc finger motif (Eulgem et al., 2000). A multiple sequence alignment of the core WRKY domain, spanning approximately 60 amino acids of all 37 RgWRKY proteins, is shown in Figure 3. A total of 36 RgWRKY proteins were found to have highly conserved WRKYGQK sequences, while RgWRKY27 varies by a single amino acid and has the variant WRKYGKK sequence, which is consistent with two Lilium regale Group IIc members (LrWRKY6 and LrWRKY7; Cui et al., 2018). It is hypothesized that variation in this WRKY domain may alter the binding specificity in the DNA targets, but this remains to be demonstrated. As previously described (Eulgem et al., 2000), the metal-chelating zinc finger motif (C–X_4-5–X_22-23–H–X₁–H or C–X₇–C–X₂₃–H–X₁–C) is another important structural characteristic of WRKY proteins. Interestingly, all zinc finger motifs of WRKR group II are C–X₇–C–X₂₃–H–X₁–C: In contrast to group III WRKY in rice (Wu et al., 2005) and barley (Mangelsen et al., 2008), there are no RgWRKY in group III proteins containing a C–X₇–C–X₂₄–H–X₁–C zinc finger motif, which have the same domain characteristic as the grape WRKR group III proteins, perhaps suggesting that this may be a feature of monocotyledonous species.

Salicylic acid and MeJA play important regulatory roles as signaling molecules in plant response to adversity stress and the accumulation of secondary metabolites (Pieterse et al., 2009; Danaee et al., 2015). Our previous research found that SA and MeJA could promote the accumulation of acteoside in the R. glutinosa hairy roots at appropriate concentrations (Wang et al., 2017a). Relative to untreated hairy roots, SA induced significant increase in acteoside content at a wide range of concentration (10 – 40μmol/L), but MeJA only slightly increased accumulation of acteoside at lower concentration (5μmol/L). To identify WRKY genes in response to SA and MeJA treatments, the expression levels of these 37 RgWRKY genes were determined by comparing 12 and 24h samples with the control sample suing RNA-seq analysis. The results indicate that these RgWRKY genes showed differential expression patterns under SA and MeJA treatments in hairy roots of R. glutinosa (Figure 4A; Supplementary Table S2). Twelve and 10 RgWRKY genes were found to be differentially expressed at atleast one time point after MeJA and SA treatments, respectively. Four genes, RgWRKY13, RgWRKY15, RgWRKY35, and RgWRKY37, were differentially expressed and upregulated at 12 and 24h after SA treatment, and the log2 (0/12h) and log2 (0/24h) for RgWRKY35 and RgWRKY37 were both greater than 2.3. There were six genes differentially expressed at 12 and 24h after MeJA treatment. Among them, RgWRKY9, RgWRKY31, and RgWRKY34 were all upregulated, and their log2 (0/12h) and log2 (0/24h) were greater than 2.0, while RgWRKY15, RgWRKY18, and RgWRKY37 were all downregulated. Interestingly, both RgWRKY15 and RgWRKY37 were differentially expressed after SA and MeJA treatments, showing upregulated expression after SA treatment but inhibited expression after MeJA treatment. In addition, RgWRKY34 significantly upregulated differential expression at 12h after SA treatment and two time points after MeJA treatment, and RgWRKY23 was significantly differentially expressed at 12h after SA and MeJA treatment.

Based on their expression patterns of the 37 RgWRKY genes in R. glutinosa hairy roots under SA and MeJA treatments, eight RgWRKY genes were selected for qRT-PCR analysis to detect their relative expression levels at 3, 9, 12, and 24h after SA and MeJA treatment. Similar to the RNA-seq analysis, qRT-PCR results showed that eight RgWRKY genes were significantly upregulated at 12 and 24h after SA treatment (Figure 4B). Among them, RgWRKY5, RgWRKY7, RgWRKY22, RgWRKY23, and RgWRKY34 genes had the highest expression levels at 3 or 9h after SA treatment, and RgWRKY13 had higher expression levels at 12 and 24h after SA treatment. The expression levels of RgWRKY35 and RgWRKY37 genes showed a highly significant increase after SA treatment, while the expression levels of RgWRKY37 increased by 10.45-fold–19.03-fold, and RgWRKY35 increased by 23.92-fold–69.27-fold. After MeJA treatment, the expression levels of RgWRKY23 and RgWRKY34 were significantly upregulated at 4h, and the other six RgWRKY genes were all insensitive to the induction or downregulated. RgWRKY23 and RgWRKY34 are two genes that simultaneously respond to SA and MeJA treatments. In addition, the expression of RgWRKY27 gene was detected by qRT-PCR, the log2 (0/12h) and log2 (0/24h) of which were greater than 6.0, but the statistical value was not significant after SA treatment. The results of qRT-PCR analysis showed that the expression level of RgWRKY27 in R. glutinosa hairy roots was significantly higher than that of the control at 3, 12, and 24h after SA treatment (Supplementary Figure S2). The results of the above expression patterns suggest that different WRKY genes may participate in the biosynthesis of acteoside by responding to different elicitors.

Hydrogen peroxide (H₂O₂), as a phytohormone, plays important roles in regulating plant growth and development (Černý et al., 2018). In previous studies, exogenous application of H₂O₂ induced SA biosynthesis in Nicotiana tabacum leaves (Leon et al., 1995). In addition, SA treatment promotes the production of H₂O₂ in A. thaliana (Rao et al., 1997). The content of acteoside in the hairy roots of R. glutinosa was detected at 12, 24, and 36h after H₂O₂ treatment, and it was found that the acteoside content at 24 and 36h after H₂O₂ treatment was 1.52 and 1.84 times that of the control, respectively (Supplementary Figure S3). RNA-seq analysis showed that 115, 77, and 123 DGEs in Control-vs-H12, Control-vs-H24, and Control-vs-H36 were enriched in the KEGG pathway of phenylpropanoid biosynthesis, respectively (Supplementary Figure S4, Supplementary Table S3). Overall, 37 WRKY genes had different response patterns to H₂O₂ induction (Figure 5A, Supplementary Table S4). There were four, three, and six WRKY genes that were differentially expressed at 12, 24, and 36h after H₂O₂ treatment, respectively. Only the RgWRKY33 gene was differentially expressed at the three time points after H₂O₂ treatment. RgWRKY9 and RgWRKY35 were significantly upregulated at least two time points after H₂O₂ treatment. Furthermore, it was found that RgWRKY13, RgWRKY22, RgWRKY28, RgWRKY34, and RgWRKY37 were significantly upregulated at least two time points after H₂O₂ treatment. qRT-PCR analysis revealed that RgWRKY34, RgWRKY35, and RgWRKY37 showed similar expression profiles under H₂O₂ treatment and had the highest expression levels at 24h after H₂O₂ treatment (Figure 5B).

In previous studies, we found that the content of medicinal ingredients (i.e., catalpol and acteoside) presented typical tissue characteristics (Wang et al., 2017b; Zhi et al., 2018). We used HPLC method to detect the acteoside content in different R. glutinosa tissues, including tender leaf (L1), fully expanded leaf (L3), old leaf (L5), top of stem (S1), middle piece of stem (S2), lower stem (S3), SS, HTR, MTR, YB, MaB, and MF. The determined results indicated that the acteoside contents were higher in leaves and floral organs than in stems, seed stock, and tuberous roots (Figure 6).

To investigate the tissue-specific expression of R. glutinosa WRKY genes, RNA-seq analysis was used to determine the expression patterns of 37 RgWRKY genes in these tissues. The results showed that the expression of all 37 RgWRKY genes was detected in at least one of the 12 examined tissues with FPKM values greater than 1 (Supplementary Table S5). Among them, five genes, RgWRKY1, RgWRKY8, RgWRKY11, RgWRKY20, and RgWRKY26, showed high expression (|log2(FPKM)|>1) in all 12 tissues and organs (Figure 7, Supplementary Table S5). The expression pattern of RgWRKY genes indicated that many of them showed a relatively high expression in the old leaf (Figure 7A). Interestingly, the expression of key enzyme genes, including RgC3H, RgC4H, RgCuAO, RgHCT, RgPAL, RgPPO, RgTyDC, and RgUGT (CL592.Contig1), in the acteoside biosynthesis pathway, was higher in the old leaf (Supplementary Table S6). The results of qRT-PCR analysis showed that RgWRKY23 and RgWRKY34 were predominantly expressed in old leaf (Figure 7B). RgWRKY7, RgWRKY35, and RgWRKY37 were more highly expressed in old leaf (Figure 7B), which is similar to the acteoside content in various tissues of R. glutinosa. Correlation analysis of acteoside content and gene expression levels of the 12 tissues of R. glutinosa revealed that the FPKM values of 10 RgWRKY genes were positively correlated to acteoside content (Supplementary Table S7). However, RgWRKY13 and RgWRKY22 showed high expression levels in seed stocks and tuberous roots, especially in the MTRs (Figure 7B), which was obviously different from the accumulation pattern of the acteoside. In conclusion, these results suggest that RgWRKY7, RgWRKY23, RgWRKY34, RgWRKY35, and RgWRKY37 might be involved in the biosynthesis of acteoside.

Based on the correlation analysis of expression levels of the 37 RgWRKY genes and acteoside content in hairy roots of R. glutinosa under SA and H₂O₂ treatments (Supplementary Tables S8, S9) and their correlation coefficients in 12 tissues of R. glutinosa plants (Supplementary Table S7), we found the expression pattern of RgWRKY37 was highly associated with acteoside content, so we selected RgWRKY37 for further studies to verify the functions of WRKY genes in acteoside biosynthesis. We constructed an overexpressing RgWRKY37 recombinant plasmid driven by the CaMV 35S promoter and transformed it into hairy root via Agrobacterium tumefaciens-mediated transformation. In total, 19 transgenic hairy root lines of R. glutinosa were obtained. Among them, six transgenic R. glutinosa hairy root lines grew well in MS medium containing kanamycin (Figure 8A) and were confirmed as positive lines by PCR analysis (Figure 8B). It was observed that the hairy root lines overexpressing the RgWRKY37 gene were slight darker than wild-type hairy root lines (Supplementary Figure S5A), and it suggest that the overexpressed RgWRKY37 influenced the development and secondary metabolism of R. glutinosa hairy roots. We also found that the dry weight of the 35S-RgWRKY37 hairy root lines was varied (Supplementary Figure S5B), and it was probably that growth rate of the different hairy root lines was inconsistent with each other; this phenotype may be irrelevant to RgWRKY37. To further analysis the function of RgWRKY37 in acteoside biosynthesis, the content of acteoside in transgenic and wild-type hairy root lines was determined by HPLC. The result showed that the acteoside content of the six tested transgenic hairy roots increased by 74.44–147.11% compared to the wild type (Figure 8C). Furthermore, the total PhGs content was significantly increased in the R. glutinosa hairy roots of RgWRKY37 overexpressed (Figure 8D). The results showed that overexpression of RgWRKY37 could enhance the accumulation of acteoside and total PhGs in R. glutinosa hairy roots.

To further demonstrate that RgWRKY37 was positively regulated acteoside biosynthesis, two positive transgenic hairy root lines, lines 13 and 19, were randomly selected for qRT-PCR analysis. The relative expression levels of RgWRKY37 were found to be significantly increased in lines 13 and 19 compared to the control (Figure 9A). Furthermore, we performed qRT-PCR analysis of the relative expression levels of enzyme genes involved in the acteoside biosynthesis. The results showed that the relative expression of genes encoding UDP-glucose glucosyltransferase (UGT, CL592.Contig1), polyphenol oxidase (PPO), copper-containing amine oxidase (CuAO), 4-coumarate-CoA ligase (4CL), shikimate O-hydroxycinnamoyltransferase (HCT), and tyrosine decarboxylase (TyDC) was significantly upregulated in the two transgenic hairy root lines compared to the wild type (Figures 9B,C).

To determine the relationship among the RgWRKY37, enzyme genes, and metabolites, gene expression levels of three transgenic hairy root lines (line 10, 13, and 19) and three wild-type hairy root lines (line 1, 2, and 3) were detected by RNA-seq analysis. Similar to the results of qRT-PCR detections, the expression levels of RgWRKY37 and eight structural genes (RgUGT, RgC4H, RgTyDC, Rg4CL, RgPPO, RgALDH, RgHCT, and RgCuAO) were generally higher in transgenic hairy lines than in wild-type lines (Figure 10A; Supplementary Table S10). Pearson’s correlation analysis was conducted between the FPKM values of RgWRKY37 and enzyme genes and between FPKM values of enzyme genes and acteoside content in transgenic hairy roots and wild-type lines. Gene-metabolite correlation network was constructed for anchoring target genes of RgWRKY37. The results showed that RgTyDC, RgC4H, Rg4CL, RgCuAO, and RgHCT were positively correlated with acteoside, while RgPPO and RgUGT were weakly correlated with it (Figure 10B). RgPPO and Rg4CL were coexpressed with RgWRKY37, and RgHCT, RgC4H, RgCuAO, RgUGT, and RgTyDC were weakly correlated with RgWRKY37 (Figure 10B).

To further investigate the relationship of RgWRKY37 and acteoside biosynthetic genes in elicitor treatments and in different tissues, three gene-metabolite correlation networks were constructed with acteoside content, FPKM values of RgWRKY37 and the 10 enzyme genes in SA- and H₂O₂-treated R. glutinosa hairy roots and 12 tissues, respectively. Three enzyme genes including RgHCT, RgUGT, and RgPPO were positively correlated with acteoside content and RgWRKY37 levels in SA-treated hairy roots, while RgTyDC and RgC4H were negatively correlated with them (Figure 10C). The coexpression of RgHCT, RgUGT, and RgTyDC with acteoside content and RgWRKY37 level in H₂O₂-treated hairy roots showed weak positive correlation, while RgPPO and RgC4H showed weak negative correlation with them (Figure 10D). It is presumed that the mechanism of acteoside accumulation in R. glutinosa hairy roots in response to H₂O₂ treatment may be different from SA treatment, and it is still unclear how RgWRKY37 regulate acteoside biosynthesis under SA and H₂O₂ treatments. In addition, there was no enzyme gene significantly correlated with RgWRKY37 in different tissues of R. glutinosa plant, and it could be the lower FPKM values of RgWRKY37 in most of the tissues (Supplementary Figure S6). Based on the above analysis, we concluded that RgUGT, RgHCT, RgPPO, and RgTyDC may be the target genes of RgWRKY37 in acteoside biosynthesis pathway.

The cis-acting element W-box in the promoter of WRKY regulated genes is predominant binding motif for plant WRKYs (Chen et al., 2018). Promoter sequences are necessary to verify whether RgWRKY37 protein activates their transcription of the candidate enzyme genes. However, R. glutinosa is lack of a high-quality reference genome, so only partial promoter sequences of RgUGT, RgPPO, and RgTyDC were obtained. Even so, six and two putative W-box sequences were identified in the promoter of RgHCT and RgUGT, respectively (Supplementary Figure S7), suggesting that RgWRKY37 might directly induce RgHCT and RgUGT expression.

All these data suggest that RgWRKY37 is indeed involved in regulating the biosynthesis of acteoside in R. glutinosa hairy roots.

The PSORT program was used to predict the subcellular localization of RgWRKY37, and predictions showed that the RgWRKY37 protein contains two putative nuclear localization signals (⁹³PAPKDRR⁹⁹ and ¹⁰³KRRK¹⁰⁶; Supplementary Figure S8). To confirm this prediction, the RgWRKY37 protein was fused to GFP driven by the CaMV 35S promoter to construct a plant expression vector (Figure 11A). Transient transformation of tobacco leaves with the recombinant vector showed that the RgWRKY37-GFP construct was predominantly localized to the nucleus in tobacco protoplasts, while the fluorescence of the GFP construct was presented throughout the cytoplasm and nucleus (Figure 11B). This suggests that RgWRKY37 is located in the nucleus, consistent with its function as a transcription factor.