Transcriptome-wide identification of WRKY transcription factors and their expression profiles in response to methyl jasmonate in Platycodon grandiflorus

ABSTRACT

Platycodon grandiflorus, a perennial flowering plant widely distributed in China and South Korea, is an excellent resource for both food and medicine. The main active compounds of P. grandiflorus are triterpenoid saponins. WRKY transcription factors (TFs) are among the largest gene families in plants and play an important role in regulating plant terpenoid accumulation, physiological metabolism, and stress response. Numerous studies have been reported on other medicinal plants; however, little is known about WRKY genes in P. grandiflorus. In this study, 27 PgWRKYs were identified in the P. grandiflorus transcriptome. Phylogenetic analysis showed that PgWRKY genes were clustered into three main groups and five subgroups. Transcriptome analysis showed that the PgWRKY gene expression patterns in different tissues differed between those in Tongcheng City (Southern Anhui) and Taihe County (Northern Anhui). Gene expression analysis based on RNA sequencing and qRT-PCR analysis showed that most PgWRKY genes were expressed after induction with methyl jasmonate (MeJA). Co-expressing PgWRKY genes with triterpenoid biosynthesis pathway genes revealed four PgWRKY genes that may have functions in triterpenoid biosynthesis. Additionally, functional annotation and protein-protein interaction analysis of PgWRKY proteins were performed to predict their roles in potential regulatory networks. Thus, we systematically analyzed the structure, evolution, and expression patterns of PgWRKY genes to provide an important theoretical basis for further exploring the molecular basis and regulatory mechanism of WRKY TFs in triterpenoid biosynthesis.

KEYWORDS:

• Platycodon grandiflorus
• WRKY transcription factors
• gene family
• expression patterns
• methyl jasmonate stress

1 Introduction

The dried roots of Platycodon grandiflorus (Jacq.) A.DC is commonly used as a traditional medicine across northern and southern regions of China.¹ P. grandiflorus contains triterpenoid saponins, polysaccharides, flavonoids, phenols, sterols, polyalkynes, fatty acids, trace elements, and other chemical components, among which triterpenoid saponins are the main active compounds.^2,3 P. grandiflorus also has various pharmacological effects, including anti-inflammatory,^4,5 anti-tumor,^6,7 and hepatoprotective ^8,9 effects. Owing to changes in the environment and climate, the chemical content of P. grandiflorus varies in different regions. At the same time, the growth and development of plants are also regulated by a variety of hormones. Among them, methyl jasmonate (MeJA) is widely present in plants and activates the coordinated expression of a series of response genes through the transduction of JA signaling, thereby regulating plant growth and development, stress response and secondary metabolism.^10–13 For example, MeJA can increase the content of artemisinin by increasing the expression level of AaWRKY9.¹⁴ LrWRKY4 and LrWRKY12 were responded to MeJA and enhanced Lilium regale resistance to gray mold.¹⁵ PpWRKY46 and PpWRKY53 can interact to regulate MeJA-mediated energy metabolism in peach disease resistance.¹⁶

WRKY TFs are one of the largest TF families in plants. They contain a highly conserved WRKY motif that regulates various physiological processes, including stress, growth, and development, thus forming a complete and complex molecular signaling network.^17,18 These factors feature a DNA-binding domain, a highly conserved WRKY motif consisting of 60 amino acids (aa) with a highly-conserved N-terminal WRKYGQK polypeptide sequence, and a C-terminal zinc finger motif, which can be classified into C2H2 and C2HC types.^19,20 Based on the number of WRKY TF structural domains and the zinc finger structure type, WRKY TFs are classified into three categories:²¹ Class I contains two WRKY motifs with a C2H2-type zinc finger; Class II contains one WRKY motif with a C2H2-type zinc finger structure, which can be further divided into five subgroups named IIa, -b, -c, -d, and -e based on phylogenetic relationships; and Class III contains one WRKY motif with a C2HC-type zinc finger.^22–24

The WRKY transcription factor family was first identified as SPF1 in the cDNA of sweet potatoes (Lpomoea batatas) in 1994.²⁵ Over time, WRKY TFs have been discovered various plants, such as Arabidopsis thaliana,²⁶ Oryza sativa,²⁷ Phakopsora pachyrhizi,²⁸ and Panax ginseng.²⁹ WRKY TFs regulate plant growth and development, secondary metabolic pathways, and responses to biotic and abiotic stressors.³⁰ Recent studies have suggested that WRKY TFs in help the seed germination,³¹ the formation and development of plant roots, stems, and leaves,³² as well as growth and development, including reproduction ^33,34 and senescence.^35,36 WRKY TFs also function in secondary metabolic pathways of several medicinal plants. For example, the AaWRKY1 gene in Artemisia annua has been shown to regulate of artemisinin biosynthesis,³⁷ and overexpression of the PqWRKY1 gene activates the triterpene biosynthetic pathway in Panax quinquefolius.³⁴ WRKY are also indispensable to abiotic and biotic stresses, which enhance the resistance of Arabidopsis to pathogens ^38,39 and increase the tolerance of plants to drought,^40,41 high temperatures,^42,43 low temperatures,^44,45 and other stresses.

Our research has shown that the content of P. grandiflorus triterpenoids from the southern regions of China is higher than that of the northern regions.^46–48 We also identified the transcriptomes of P. grandiflorus in combination with PacBio and DNBSEQ sequencing platforms.⁴⁹ The WRKY gene family, a critical transcriptional factor in regulating secondary metabolism and stress response in plants, plays a vital role in the growth and development of P. grandiflorus. However, only a few studies have been performed on the WRKY gene family of P. grandiflorus. The expression trends of this gene family in different parts of P. grandiflorus and their response to abiotic stress remain unclear.

We conducted a comprehensive investigation of the WRKY gene family in P. grandiflorus. We identified 27 P. grandiflorus WRKY genes and analyzed sequence alignments, conserved motifs, phylogenetic trees, GO annotations, and protein interactions. RNA-sequencing (RNA-seq) data were used to study the expression profile of PgWRKY genes in different tissues and their response trends under MeJA treatment. The expression of selected PgWRKY genes under MeJA treatment was verified by quantitative real-time PCR (qRT-PCR). Based on the tissue-specific expression of the PgWRKY gene and co-expression of genes involved in the triterpenoid biosynthesis pathway of P. grandiflorus, we proposed the potential regulatory TFs involved in the triterpenoid biosynthesis in P. grandiflorus, which laid a foundation for further study of the WRKY gene family in P. grandiflorus.

2 Materials and methods

2.1 Plant materials and sample collection

Biennial Platycodon grandiflorus were collected from Tongcheng City (Southern Anhui) and Taihe County (Northern Anhui) respectively, and were identified as P. grandiflorus by Professor Peng Huasheng. It was separated into roots, stems, and leaves with three biological replicates per organ. After a quick rinse with sterile water, samples were snap frozen in liquid nitrogen and stored at −80°C for subsequent analysis.

The material used for hormone treatment was Tongcheng P. grandiflorus tissue culture seedlings provided by College of Pharmacy, Anhui University of Chinese Medicine, Hefei, China, and treated with 50 μM MeJA for 0, 3, 6, 9, 12, 24 and 48 h. Root tissues of plants were harvested, all treatments were three biological replicates. Freeze immediately in liquid nitrogen and store at −80°C.

2.2 Analysis of saponin contents in P. grandiflorus

The root, stem and leaf tissues were dried under 55°C, ground into powder and passed through a 50-mesh sieve. 10 mg powder samples were weighed and extracted with 8 mL of 70% methanol by ultrasonic extraction for 1 h (100 W, 40 Hz). The extracted samples were passed through 0.45 μm filter membrane and injected into the column.⁴⁶ The contents of eight Platycodon saponins were determined by ultra-high performance liquid chromatography-Orbitrap-tandem mass spectrometry (UHPLC-Orbitrap-MS/MS). All analyses were performed on an Ultimate HPLC system (Thermo Fisher Scientific, Waltham, MA, United States). The samples were separated on an Agilent Eclipse XDB-C18 column (4.6 mm × 250 mm, 5 μm). The mobile phase was 0.075% acetic acid water (A)-methanol (B). The gradient curves were as follows: 0 ~ 3 min, 90 ~ 80% A; 3 ~ 11 min, 80 ~ 77% A; 11 ~ 20 min, 77 ~ 5% A; 20 ~ 25 min, 5 ~ 90% A. The mobile phase flow rate was 0.3 mL/min, and the injection volume was 2 μL.

2.3 Identification of WRKY Gene family in P. grandiflorus

Based on the existing transcriptome sequencing data of the research group (SRA accession number: PRJNA688328), all presumed WRKY proteins of P. grandiflorus were retrieved by searching against the transcriptome annotation sequence. The sequences with complete WRKY domains were selected, and the homology of amino acid sequence above 97 % were deleted. The amino acid sequences of WRKY in A. thaliana published on the Arabidopsis Information Resource website (http://www.arabidopsis.org/) were downloaded and aligned to remove sequences with incomplete WD via the MEGA-X. Subsequently, all validated candidate genes with probable WRKY domains were corroborated using the SMART program (http://smart.embl-heidelberg.de) to obtain sequences with the WRKY domain only. The putative WRKY genes from P. grandiflorus were named PgWRKY1-PgWRKY27 based on their relationship with the A. thaliana sequence. The isoelectric points (pIs) and theoretical molecular weights (MWs) of each PgWRKY protein were predicted online using the Expasy portal (https://web.expasy.org/protparam/). Cell localization was predicted using the WoLF PSORT tool (https://wolfpsort.hgc.jp/).

2.4 Multiple sequence alignment and protein structure analysis

Accurately grouping gene families is an integral part of the functional analysis of transcription factors. Based on the AtWRKYs classification in A. thaliana, 72 AtWRKYs were downloaded from the Arabidopsis Information Resource (https://www.arabidopsis.org/), and all identified PgWRKYs were divided into different groups according to the classification of AtWRKYs. Multiple sequence alignment of domains in PgWRKYs was executed using the ClustalW program with default parameters,⁵⁰ and the results were colored using the GeneDoc tool. The conserved motifs of PgWRKYs were analyzed using multiple expectation maximization for motif elicitation (MEME: http://meme-suite.org/tools/meme) with the following parameters: repetitive time was “any,” maximum motif number was 10, and motif width was between 6 and 50 residues.^51,52 The MEME results were displayed using TBtool software.⁵³

2.5 Phylogenetic analysis of PgWRKY

The amino acid sequences of WRKY in P. grandiflorus and A. thaliana were aligned using ClustalW,⁵⁰ and then imported to construct an evolutionary relationship tree using the maximum likelihood method.⁵⁴ The substitution model, equal input model, rates among sites, uniform rates, gaps, missing data treatment, and partial deletion parameters were used in phylogenetic tree construction. The trees were visualized and optimized using the iTOL tool (https://itol.embl.de/).

2.6 Expression profile analysis of PgWRKY genes

We analyzed the expression of WRKY genes in different tissues of P. grandiflorus and their response to MeJA treatment. Our group completed transcriptome sequencing of different tissues of P. grandiflorus (SRA accession number: PRJNA688328), and 12 RNA-Seq datasets for MeJA processing were obtained from NCBI (SRX5506849, SRX5506850, SRX5506841, SRX5506842, SRX5506839, SRX5506840, SRX5506845, SRX5506846, SRX5506843, SRX56847, SRX550656847, SRX55065684RX). Clean reads obtained from RNA-seq were mapped to the reference full-length transcripts using the Bowtie2 tool.⁵⁵ The expression levels were then calculated for each sample using RNA-seq by expectation-maximization (RSEM) software, and the transcript per kilobase fragment mapping (FPKM) was used to normalize the read count.⁵⁶ To analyze the expression of genes related to triterpenoid synthesis pathway in P. grandiflorus, we obtained KEGG annotations from transcriptome data (SRA accession number: PRJNA688328). We screened genes related to synthesis of platycodins according to annotation information of map00900 (terpenoid backbone biosynthesis) and map00909 (sesquiterpenoid and triterpenoid biosynthesis) in KEGG database.⁴⁹ The differentially expressed genes (DEGs) in different tissues of roots, stems and leaves were further screened according to the their FPKM values. A total of 7 DEGs encoding 6 enzymes related to the MVA pathway, 8 DEGs encoding 6 enzymes involved in the MEP pathway, and 9 DEGs of 5 enzymes of the downstream of triterpenoid saponin synthesis pathway were identified. Pearson correlations between genes from the P. grandiflorus triterpenoid biosynthetic pathway and PgWRKYs of Tongcheng P. grandiflorus were calculated using the SPSS software. Heat maps were plotted using the TBtools software.⁵³

2.7 Gene expression analysis by quantitative real-time PCR (qRT-PCR)

Total RNA was extracted using RNAiso Plus (TaKaRa, Japan), according to the manufacturer’s instructions. RNA (1 μg) was reverse-transcribed to cDNA using a reverse transcription kit (PrimeScript™ II 1st Strand cDNA Synthesis Kit, Waryong, China). qRT-PCR was performed to determine the relative expression levels of PgWRKY genes under MeJA treatment, and β-actin was used as an internal control. Primer sequences are listed in Supplementary Table 1. The reaction system was prepared according to the instructions of the TB Green® Premix Ex Taq™ II kit. The PCR conditions were as follows: 95°C for 30s, 95°C for 15s, 60°C for 30s, 72°C for 15s, and 40 cycles. Three biological replicates and three technical replicates were set up, melting curves were plotted, and the relative gene expression levels were analyzed using the 2^−ΔΔCt method.⁵⁷

2.8 Gene ontology annotation and protein-protein interaction analysis

Gene ontology (GO) annotation was obtained from Beijing Genomics Institute (BGI; https://report.bgi.com/ps/login/login.html) and visualized using OmicShare (https://www.omicshare.com/tools/). The STRING 11.5 (http://string-db.org/) database was used to predict Arabidopsis WRKY protein-protein interaction networks.⁵⁸ Fourteen AtWRKY proteins were identified, with a confidence parameter of 0.4. AtWRKY proteins with interactions were mapped to PgWRKY proteins by homology.

3 Results

3.1 Identification of WRKY Gene family in P. grandiflorus

Twenty-seven candidate PgWRKY genes were identified in the full-length transcriptome of P. grandiflorus ⁴⁹ using AtWRKY genes in A. thaliana as queries. We examined the sequences of all candidate proteins using the conserved domain database (https://www.ncbi.nlm.nih.gov/Structure/cdd/wrpsb.cgi). We also verified that all 27 PgWRKY proteins contained the WRKY DNA-binding domain. As shown in Table 1, all the identified PgWRKY proteins contained at least one highly conserved heptapeptide WRKY domain; however, a difference was found in the WRKYGQK motif and the zinc-finger-like motif. Further analysis indicated that 26 PgWRKY proteins contain the conserved WRKYGQK heptapeptide domain, and the PgWRKY10 protein contains a WRKYGKK sequence. The ORF length for the PgWRKY genes ranged from 576 bp (PgWRKY10) to 2196 bp (PgWRKY2), encoding 191–731 aa. The deduced 27 PgWRKY proteins have MWs from 21.29 kDa (PgWRKY10) to 78.99 kDa (PgWRKY2) and pIs from 4.91 (PgWRKY16) to 9.66 (PgWRKY24). The subcellular localization prediction showed that most members (except PgWRKY17) were present in the nucleus. Supplementary Table 2 contains specific information on the WRKY genes of P. grandiflorus.

Table 1. Features of PgWRKY genes identified in P. grandiflorus.

Name	isorfrom	orf	AA	WRKYdomin	Zinc-finger type	Domain number	Group
PgWRKY1	isoform_28040	1548	515	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY2	isoform_67094	2196	731	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY3	isoform_124836	1560	518	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY4	isoform_1243	1532	509	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY5	isoform_3239	964	320	WRKYGQK	C2HC	1	III
PgWRKY6	isoform_56097	1108	368	WRKYGQK	C2H2	1	II b
PgWRKY7	isoform_41099	1100	365	WRKYGQK	C2H2	1	II d
PgWRKY8	isoform_40330	896	297	WRKYGQK	C2H2	1	II c
PgWRKY9	isoform_31604	1744	580	WRKYGQK	C2H2	1	II b
PgWRKY10	isoform_50511	576	191	WRKYGKK	C2H2	1	II c
PgWRKY11	isoform_22766	956	317	WRKYGQK	C2HC	1	III
PgWRKY12	isoform_19381	672	223	WRKYGQK	C2H2	1	II c
PgWRKY13	isoform_44539	1076	357	WRKYGQK	C2HC	1	III
PgWRKY14	isoform_125380	1232	409	WRKYGQK	C2H2	1	II e
PgWRKY15	isoform_107654	1036	344	WRKYGQK	C2HC	1	III
PgWRKY16	isoform_10015	796	264	WRKYGQK	C2H2	1	II e
PgWRKY17	isoform_114563	588	218	WRKYGQK	C2H2	1	II c
PgWRKY18	isoform_2554	1004	333	WRKYGQK	C2H2	1	II a
PgWRKY19	isoform_130708	1352	449	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY20	isoform_122078	1820	605	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY21	isoform_44238	1044	346	WRKYGQK	C2H2	1	II d
PgWRKY22	isoform_57380	1048	348	WRKYGQK	C2H2	1	II e
PgWRKY23	isoform_23542	1184	393	WRKYGQK	C2H2	1	II c
PgWRKY24	isoform_144642	1108	368	WRKYGQK	C2	1	II c
PgWRKY25	isoform_92941	1584	527	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY26	isoform_72652	1212	403	WRKYGQK/ WRKYGQK	C2H2	2	I
PgWRKY27	isoform_38617	640	212	WRKYGQK	C2HC	1	III

3.2 Multiple sequence alignment and structure analysis

Seven Arabidopsis proteins from different groups were randomly selected as representative sequences (Group I: AtWRKY58, Group IIa: AtWRKY40, Group IIb: AtWRKY6, Group IIc: AtWRKY56, Group IId: AtWRKY21, Group IIe: AtWRKY35, and Group III: AtWRKY46). The detailed structures of the WRKY domain and zinc-finger type are displayed in Figure 1. The WRKYGQK heptapeptide is a signature of WRKY proteins. Twenty-six PgWRKY proteins contained the highly conserved WRKYGQK sequence, while PgWRKY10 had a single amino acid substitution, K, for Q (Figure 1 and Table 1). The second motif is a zinc-finger structure containing two types of zinc finger motifs: C-X4-5-C-X22-23-H-X-H (C2H2) and C-X7-C-X23-H-X-C (C2HC). Twenty-one PgWRKY proteins included C2H2 type zinc finger motifs, five PgWRKY proteins displayed C2HC type zinc finger motifs, and a partial absence of the zinc-finger motif sequence was present in PgWRKY24. A comparison of the WRKY TFs identified in the P. grandiflorus transcriptome is presented in Supplementary Figure 1.

Figure 1. Alignment of 27 P. grandiflorus (PgWRKY) and 8 A. thaliana (AtWRKY) WRKY domain sequences. For group I WRKY proteins, N-terminal and C-terminal WRKY domains are represented by N and C, respectively. In the red solid box are WRKY domains. Zinc finger structures are marked by a red dotted box.

To gain insight into the functional regions of PgWRKY proteins, the MEME program was employed to reveal conserved motifs among 27 PgWRKY proteins. The MEME results showed that the number of motifs in the PgWRKYs ranged from 2–8 aa, and the width of the 10 identified motifs ranged from 6–41 aa. Ten motifs were identified in the structure of the PgWRKYs, and their details are shown in Supplementary Figure 2. As shown in Figure 2, motif 1 encoded the conserved WRKY domain, while motifs 2 and 3 encoded the conserved zinc-finger structure. Notably, similar motif compositions were observed in the same group of PgWRKY proteins. Most members within the same clade shared a unique motif composition compared to that in the other clades. For example, motif 1 and motif 2 were found in all 27 PgWRKY proteins, while motif 3 was not found in Group III members and was found in all other groups. Motif 7 exists only in two members of Group I and Group IIc (PgWRKY24 and PgWRKY26). Only members of Group I and Group II have motif 8. Motif 9 was found only in Group III. Generally, the closely-related PgWRKYs in the phylogenetic tree shared similar gene and motif compositions, suggesting that PgWRKYs in the same group may play similar functional roles.

Figure 2. Conserved motif analysis of PgWRKYs. The motif composition of PgWRKY was analyzed by the MEME tool. The detailed information of the ten motifs is in Supplementary Figure 2.

3.3 Phylogenetic analysis of PgWRKY

To investigate the evolution of PgWRKY family members, a phylogenetic tree was constructed using the maximum-likelihood estimate method using MEGA-X software and based on multiple alignments between full-length protein sequences of 72 AtWRKYs and 27 PgWRKYs. Phylogenetic analysis indicated that the classification of PgWRKY proteins was the same as that of A. thaliana, thereby confirming the classification accuracy of PgWRKY proteins. According to the constructed phylogenetic tree containing PgWRKYs and AtWRKYs (Figure 3), PgWRKYs were classified into three primary groups (Groups I, II, and III). Among the 27 PgWRKY protein sequences, eight were assigned to Group I, fourteen to Group II, and five to Group III. Group II had the largest number of WRKY proteins and was divided into five major subgroups: one PgWRKY protein belonged to IIa (PgWRKY18), two to IIb (PgWRKY6 and PgWRKY9), six to IIc (PgWRKY8, PgWRKY10, PgWRKY12, PgWRKY17, PgWRKY23, and PgWRKY24), two to IId (PgWRKY7 and PgWRKY21) and, three to IIe (PgWRKY14, PgWRKY16, and PgWRKY22). Subgroups IIa and IIb were two subgroups in the same branch, whereas subgroups IId and IIe were derived from one clade. Notably, among these groups or subgroups, most members of the PgWRKYs were clustered in Group I. Overall, the classification of PgWRKYs confirms their diversification, which suggests that different family members may have varied functions.

Figure 3. Phylogenetic tree of the total WRKY sequences from P. grandiflorus and A. thaliana. The WRKY sequences were used for phylogenetic analysis using the MEGA-X software. The arcs with different colors represent seven subgroups of WRKY proteins. The solid black star and hollow circle represent WRKY sequences from P. grandiflorus and A. thaliana, respectively.

3.4 PgWRKYs expression profiles in different tissues of P. grandiflorus

In this study, RNA-seq analysis was used to study the potential physiological functions of PgWRKY TFs (Figure 4 and Supplementary Table 3). The expression levels of 27 genes in the roots, stems, and leaves of P. grandiflorus collected from Tongcheng southern Anhui (TC), and Taihe northern Anhui (TH), were investigated. Five genes were highly expressed in roots, and PgWRKY8 and PgWRKY17 in Group IIc were highly expressed in the roots of TC. PgWRKY1, PgWRKY2, and PgWRKY4 were highly expressed in the roots of TH, as were all members of Group I. Furthermore, PgWRKY6 was highly expressed in leaves, whereas eight genes (PgWRKY9, PgWRKY13, PgWRKY15, PgWRKY18, PgWRKY22, PgWRKY24, PgWRKY25, and PgWRKY27) were highly expressed in the stems. Based on the expression abundance of different groups, most members of Group III were highly expressed in all tissues (FPKM>5), and the genes in Groups I and II showed tissue-specific expression patterns, which might play an important role in the growth and development of P. grandiflorus. In addition, Figure 4B shows similar trends in the expression of roots from the two origins, differing in that the expression of PgWRKY genes was more concentrated in the roots from Tongcheng City. The expression trend of the PgWRKY gene in the stem of TC is higher than that in TH. In contrast, the PgWRKY gene showed an opposite trend with higher expression in the TH leaf compared to the TC. This indicates a difference in the expression level of PgWRKY genes in P. grandiflorus from northern and southern Anhui province.

Figure 4. Expression of PgWRKY genes in P. grandiflorus tissues. (A) The heatmap of all PgWRKY genes expression in different tissues. (B) The violin plot of all PgWRKY genes expression in different tissues.

3.5 Determination of triterpenoid saponins

We used UHPLC-Orbitrap-MS/MS method to determine root, stem and leaf samples from Tongcheng City and Taihe County (3 replicates per organ). The contents of 8 saponins were expressed as mean and standard deviation (Figure 5). In general, the content of triterpenoid saponins in roots was higher than that in stems and leaves. In different origins, the content of P. grandiflorus saponins also showed certain differences. The contents of Platycodin D, Platycodin D3, Deapioplatycodin D3 and Platycoside D in the roots and stems of TC were higher than those of TH. But the Platycodin D2 and Platicodigenin in TH Root is more than twice that of TC. Although the contents of Deapioplatycodin D and Deapioplatycodin D3 were lower in the two kinds of P. grandiflorus, the contents in TC stems and leaves were higher than those in TH. In addition, Platycodin D and Platycodin D2 were higher in TH leaves. Surprisingly, Platycoside E and Platicodigenin were not detected in the leaves from both origins. The total saponins content in roots and stems of TC was higher than that in TH, but the situation in leaves was opposite.

Figure 5. Saponin contents in different tissues of P. grandiflorus collected from Tongcheng southern Anhui (TC) and Taihe northern Anhui (TH) (mean ± SD, n = 3). ***Data compared with TC, P < 0.05.

3.6 Expression analysis of PgWRKY genes under MeJA treatment

The sequence and heat map of this experiment were based on the MeJA-treated transcriptome data (Supplementary Table 4). As shown in Figure 6, the expression levels of most PgWRKYs treated with MeJA varied over time. The WRKY expression level was upregulated after 12 h, whereas it was downregulated at 24 and 48 h. Specifically, the transcription levels of nine genes (PgWRKY1, PgWRKY6, PgWRKY9, PgWRKY14, PgWRKY17, PgWRKY19, PgWRKY23, PgWRKY24 and PgWRKY26) increased rapidly under MeJA treatment, while those of the other 7 genes (PgWRKY2, PgWRKY3, PgWRKY8, PgWRKY10, PgWRKY11, PgWRKY20, and PgWRKY21) were lower under MeJA treatment than that in the control. In addition, the expression performance of some genes changed significantly in specific time frames. For example, the expression levels of PgWRKY2 and PgWRKY8 at 24 h were lower than those in the control but were significantly upregulated at 48 h.

Figure 6. Expression of PgWRKY genes in response to MeJA treatment. (A) The expression of PgWRKY genes for MeJA treatments. (B) Violin plot of PgWRKY genes before and after MeJA treatment.

Fourteen PgWRKY genes were initially selected based on their expression levels in the RNA-seq data, and their performance was further validated by qRT-PCR analysis (Figure 7). The expression trends of PgWRKY8, PgWRKY11, PgWRKY16, and PgWRKY26 were consistent with the transcriptome data. There were certain differences between the actual observations and transcriptome data, which can be explained by the individual differences between the material transcriptome sequencing and qRT-PCR analysis. Notably, the expression levels of PgWRKY17, PgWRKY18, and PgWRKY26 increased rapidly at 3 h. After 9 h of treatment, PgWRKY11, PgWRKY14, PgWRKY15, and PgWRKY16 also increased rapidly. From these data, it can be deduced that MeJA regulates the expression of these genes.

Figure 7. qRT-PCR analysis of PgWRKY genes under MeJA treatment. The Y- and X-axis represent the relative expression level and the time course of MeJA treatment, respectively. Roots were sampled at 0, 3, 6, 9, 12, 24, and 48 h after MeJA treatments. Data represent the mean ± SD of three technical repeats.

3.7 Co-expression analysis of candidate triterpenes saponin biosynthesis and PgWRKY genes

WRKY TFs are involved in various plant physiological activities. To elucidate the relationship between PgWRKYs and triterpenoid biosynthesis in P. grandiflorus, we constructed a co-expression network of PgWRKYs and genes related to the upstream 2-C-methyl-D-erythritol-4-phosphate (MEP) and mevalonate (MVA) pathways and downstream triterpene skeleton formation pathway (Figure 8 and Supplementary Table 5). Co-expression analysis suggested that the correlation between PgWRKYs and genes related to the biosynthesis of triterpenoids in P. grandiflorus could be divided into three main clusters. In cluster I, several PgWRKY genes showed strong correlations with upstream genes in the MEP pathway, including PgDXS1, PgDXS2, PgDXS3, and PgMCT. The PgWRKYs in cluster I also showed strong correlations with MVA and some downstream genes, such as PgAACT1, PgMK, PgMDC, PgGPPS2, and PgSQS. In contrast, PgWRKYs in cluster II showed a low correlation with genes in the P. grandiflorus triterpenoid biosynthetic pathway and only correlated strongly with upstream PgHDS and PgHMGR. Some PgWRKYs in cluster III showed a strong positive correlation with pathway genes. PgDXS1 and PgCMK in the MEP pathway and PgIDI, PgHMGS, and PgMDC in the MVA pathway were strongly positively correlated with PgWRKYs in cluster III. Downstream, PgGPPS1, PgSQE1, Pgβ-AS1, and Pgβ-AS3 also showed a strong correlation with PgWRKYs in cluster III. PgWRKYs in cluster III may be involved in and regulate the biosynthesis of triterpenoids in P. grandiflorus. Co-expression analysis showed that MEP pathway genes were relevant to PgWRKYs in cluster I, and the MVA pathway and downstream pathway genes were relevant to PgWRKYs in cluster III. Four genes, PgWRKY2, PgWRKY9, PgWRKY10, and PgWRKY24, were selected for their modulatory potential in triterpene biosynthesis.

Figure 8. The Pearson’s correlation coefficients of PgWRKYs with triterpenoid biosynthesis pathway in P. grandiflorus. 1-Deoxy-d-xylulose-5-phosphate synthase (PgDXS); 2-C-methyl-d-erythritol 4-phosphate cytidylyltransferase (PgMCT); 4-diphosphocytidyl-2-C-methyl-D-erythritol kinase (PgCMK); 2-C-methyl-d-erythritol 2,4-cyclodiphosphate synthase (PgMCS); (E)-4-Hydroxy-3-methylbut-2-enyl-diphosphate synthase (PgHDS); 4-Hydroxy-3-methylbut-2-en-1-yl diphosphate reductase (PgHDR); Isopentenyl diphosphate isomerase (PgIDI); Acetyl-CoA C-acetyltransferase (PgAACT); Hydroxymethylglutaryl-CoA synthase (PgHMGS); Hydroxymethylglutaryl-CoA reductase (PgHMGR); Mevalonate kinase (PgMK); Diphosphomevalonate decarboxylase (PgMDC); Geranyl pyrophosphate synthase (PgGPPS); Farnesyl diphosphate synthase (PgFPPS); Squalene synthase (PgSQS); Squalene monooxygenase (PgSQE); β-Amyrin synthase (Pgβ-AS).

3.8 GO annotation and interaction analysis of specific PgWRKY proteins

To further explore the biological and molecular functions of PgWRKY proteins, GO annotation and protein-protein interaction (PPI) network analysis were performed on the 27 PgWRKY proteins in this study. The GO annotation results consisted of biological processes (BP), cellular components (CC), and molecular functions (MF) (Figure 9). Most PgWRKYs regulate cellular, metabolic, developmental, reproductive, and multicellular organismal processes, as well as responses to stimuli. The molecular functions of PgWRKYs are related to transcriptional regulatory activity and binding. Cellular components of this protein family include organelles and cells.

Figure 9. GO analysis of PgWRKY proteins. The results are grouped into three main categories: biological process, cellular component, and molecular function. The X-axis indicates the number of genes.

The STRING 11.5 online software was used to construct a WRKY protein-protein interaction network to discover the presence of PPIs. We identified five high confidence interacting proteins in the Arabidopsis WRKY family: VQ proteins (such as Meckel syndrome, type 1 [MKS1] and sigma factor binding protein 1 [SIB1]) involved in regulating plant defense responses, mitogen-activated protein kinase (MPK) 3 and MPK4 proteins involved in plant responses to pathogens and stresses, and genomes uncoupled 5 (GUN5) proteins involved in abiotic stress responses. As shown in Figure 10, the PgWRKY33 protein is highly homologous to Arabidopsis WRKY33, suggesting that it may have potentially stronger interactions with most plant defense proteins (MPK3, MPK4, MKS1, and SIB1). Similarly, the PgWRKY18 protein is highly homologous to Arabidopsis WRKY40 and is presumed to have a strong interaction with the internal members of the WRKY family.

Figure 10. Protein-protein interaction network of specific PgWRKY proteins. Black and red color characters represent A. thaliana and P. grandiflorus, and the thick lines represent the strength of interaction.

4 Discussion

The WRKY gene family is one of the largest TF families that regulate plant growth and secondary metabolism.^18,59 Highly advanced high-throughput sequencing technology and available public databases have facilitated the comprehensive study of several gene families.⁶⁰ The WRKY family has been studied in several model plants, cash crops, and medicinal plants such as Arabidopsis thaliana,⁶¹ Sorghum bicolor,⁶² Zea mays,⁶³ Panax ginseng,²⁹ Isatis tinctoria,⁶⁴ and Andrographis paniculata.⁶⁵ Our previous studies combined second- and third-generation sequencing technologies to systematically analyze the transcriptomes of different organs of P. grandiflorus.⁴⁹ Based on transcriptomic data in P. grandiflorus, 27 PgWRKYs were identified of P. grandiflorus for the first time. This study provides a reference for a further functional understanding of the WRKY gene family.

Most PgWRKY proteins contain a highly conserved heptapeptide at the N-terminus (WRKYGQK). However, at heptapeptide variant WRKYGKK was found in PgWRKY10. These differences may affect the ability of WRKY TFs to bind to the W-box. For example, soybean WRKY TFs with WRKYGKK do not properly bind to the W-box.⁶⁶ NtWRKY12 in Nicotiana tabacum has the WRKYGKK motif, and these transcription factors bind to the TTTTCCAC sequence but not to the W-box (TTGACT/C).⁶⁷ Therefore, we hypothesize that the heptapeptide variants may alter the binding properties of TFs to W-box and affect the transcriptional levels of PgWRKY transcription factors-targeted genes.

Based on sequence alignment and phylogenetic tree analysis (Figure 1 and Figure 3), the WRKY family members were grouped similarly to Arabidopsis: Group I, Group II, and Group III, while Group II was further subdivided into subgroups IIa, IIb, IIc, IId, and IIe.¹⁹ In this study, Group I had eight members, Group II had fourteen members, accounting for 51.85% of all PgWRKY genes, Group III had five members. These results are consistent with WRKY group sizes in Arabidopsis, maize, and ginseng.^19,29,63 Among these subgroups, IIc had the most members with six PgWRKY genes, accounting for 42.86% of the genes in Group II, which is similar to the results for soybean, Arabidopsis, and Andrographis.^19,28,65 Rinerson et al. ⁶⁸ proposed four main lineages of WRKY transcription factor families in flowering plants: Group I + IIc, Group IIa + IIb, Group IId + IIe, and Group III, reflecting the evolution pathway of WRKY family members. This was also observed in the P. grandiflorus WRKY gene family. For example, members of Groups IIa and IIb were in the same evolutionary branch, and members from Groups IId and IIe were grouped into the same branch, suggesting that these subgroups arose indirectly from a common ancestor.

These motifs are associated with transcriptional activity and protein-protein interactions of target genes.⁶⁹ Therefore, the characterization and functions of TFs can be identified by analyzing conserved motif information. All PgWRKY genes had more than two conserved motifs, consistent with reports from other plant species.^70,71 Each group member had the same or similar motif composition, suggesting that WRKY genes from the same group have similar protein structures and biological functions.

RNA-seq is an essential tool for gene function and structure studies at the global level to reveal the molecular mechanisms of specific biological processes. It has been widely used in basic research and drug development.^72,73 We used transcriptome data form different tissues of P. grandiflorus from two origins to determine the WRKY family gene expression (Figure 4). In the present study, two PgWRKYs of P. grandiflorus were two origins expressed differentially in roots, stems, and leaves. This indicated that PgWRKYs function in the growth and development of P. grandiflorus. All PgWRKY genes were expressed in at least one tissue (FPKM > 1), and their expression varied significantly among the different tissues. As shown in Figure 4B, PgWRKYs in roots and stems were higher than those in leaves, and the median FPKM values of root and stem tissues of Tongcheng City (south Anhui province) P. grandiflorus were higher than those of Taihe County (north Anhui province). This study also analyzed the differences in the content of 8 kinds of Platycodon saponins in different tissues of P. grandiflorus from the two origins mentioned above (Figure 5). The highest saponin content was found in root tissues, followed by stems and leaves. The content of saponins in the roots and stems of Tongcheng P. grandiflorus is higher than that of Taihe. These are consistent with the expression trend of PgWRKY genes. Therefore, PgWRKYs may regulate the metabolic process of triterpenoids.

Notably, five genes showed the highest expression levels in roots. One gene was highly expressed in the leaves, whereas eight genes were highly expressed in the stems. These tissue-specific expression patterns suggest that PgWRKYs may be involved in tissue-specific developmental and signal transduction processes. For example, roots interact with the external environment, and genes highly expressed in the roots are more likely to participate in various stress response activities. AtWRKY39, AtWRKY48, and AtWRKY57 regulate plant heat tolerance, drought resistance, and immunity against pathogens.^74–76 PgWRKY8, PgWRKY17, and their orthologs in P. grandiflorus were specifically expressed in roots. They likely regulate plant growth and development by helping the underground plant parts defend against stress. Similarly, PgWRKY8, a highly expressed gene in the roots, and its homologous AaWRKY40 in A. annua are involved in the metabolism of terpenoids.⁷⁷ This indicates that PgWRKY8 may also regulate the metabolism of terpenoids in the roots. In addition, the wheat WRKY transcription factor gene TaWRKY71-1 was specifically expressed in leaves, and plants overexpressing TaWRKY71-1 had larger leaves and better growth ability than its parent common wheat JN177.⁷⁸ It is suggested that PgWRKY6 highly expressed in leaves has similar functions. Studying tissue-specific expression of PgWRKYs can provide critical information for further understanding the function of the P. grandiflorus WRKY gene family.

MeJA is involved in several physiological processes, such as plant growth and stress response. It regulates secondary metabolism in both angiosperms and gymnosperms.^79–81 Triterpene biosynthesis is regulated by MeJA signaling, and many MeJA-responsive TFs can regulate the MeJA signaling pathway.^14,82 Therefore, WRKY genes in response to MeJA in P. grandiflorus may impact triterpenoid biosynthesis. Previous studies have shown that the accumulation of triterpenoid saponins in the hairy roots of P. grandiflorus increased after MeJA induction, and the transcription levels of triterpenoid biosynthesis pathway genes (PgHMGS, PgHMGR, PgMK, and PgMVD) increased rapidly at 3–6 h of MeJA treatment.⁸³ This study used publicly-available transcriptome data to study the PgWRKY genes expression under MeJA induction. A few genes (e.g., PgWRKY8, PgWRKY11, PgWRKY16, and PgWRKY26) showed the same expression trend as the RNA-seq analysis. But we found that the expression levels of some genes (PgWRKY15, PgWRKY17, PgWRKY18, and PgWRKY26) increased significantly at 3–6 h. This was consistent with the gene expression trend of triterpenoid biosynthesis pathway after MeJA treatment. It is speculated that PgWRKY15, PgWRKY17, PgWRKY18, and PgWRKY26 may indirectly regulate the biosynthesis of triterpenoid saponins by interacting with terpenoid synthase genes.

WRKY TFs are involved in plant development and regulate the biosynthesis of secondary metabolites. Medicinal plants with Terpenes are the main medicinal components of P. ginseng, A. apiacea, P. notoginseng, and P. quinquefolium. The terpenoid biosynthesis strategy in plants includes the MVA and MEP pathways.^84,85 They are governed by multi-level network regulation, which can be accomplished by various structural genes in the biosynthetic pathway, and secondary regulation is accomplished by TFs.⁸⁶ P. ginseng transcription factors PgWRKY4X binds to the squalene epoxidase (PgSE) promoter to upregulate ginsenoside biosynthesis-related genes and highly improve the accumulation of ginsenosides.⁸⁷ Overexpression of PqWRKY1 activates the triterpenoid biosynthesis pathway in transformed Arabidopsis.⁸⁸ In this study, we found that four genes (PgWRKY2, PgWRKY9, PgWRKY10, and PgWRKY24) were positively correlated with triterpene biosynthesis pathway genes. These four genes were distributed in phylogenetic Groups I and II, which is in accordance with the related research results for P. ginseng.²⁹ The above co-expression results reflect the potential regulatory function of WRKYs in triterpenoid biosynthesis, and overexpression or knockout analyses of these PgWRKY genes will help clarify their functions.

GO annotation of gene families and the study of potential protein-protein interaction networks help understand the functions of these genes in this family.⁸⁹ GO analysis showed that the WRKY gene family of P. grandiflorus functions in three aspects: biological process, cellular component and molecular function. In A. thaliana, most studies on protein-protein interactions of WRKY TFs have focused on the AtWRKY33 protein involved in plant defense responses.⁹⁰ For example, Arabidopsis AtWRKY33 and AtWRKY25 interact with MKS1, a VQ protein that uses AtMPK4 as a substrate. Activated MPK4 phosphorylates MKS1 and releases AtWRKY33 after infection with the bacterial pathogen Pseudomonas syringae or flagellin. The released AtWRKY33 acts on the promoter of PAD3, which encodes a biosynthetic enzyme involved in producing phytoalexins.⁹¹ In addition, dual-targeted SIB1 activates AtWRKY33 and functions in plant defenses against necrotrophic pathogens.⁹² The protein-protein interaction network of P. grandiflorus showed that PgWRKY26 was highly homologous to AtWRKY33. These results suggest that PgWRKY26 may mediate interactions between different signaling pathways and hinge on the regulatory network of plant defense responses.

5 Conclusion

In total, 27 P. grandiflorus WRKY genes were identified in this study. The classification, evolutionary relationships, gene structure, and conserved motifs of the P. grandiflorus WRKY gene family were studied. The expression patterns of PgWRKY genes in different tissues of P. grandiflorus and their response to MeJA indicated that these genes might have essential functions in the growth, development, and secondary metabolism of P. grandiflorus. Co-expression analysis revealed a potential regulatory role of PgWRKYs in triterpenoid biosynthetic pathway genes. GO annotation and PPI analysis of PgWRKY proteins were performed to explore the functions and regulatory mechanisms of PgWRKY genes in P. grandiflorus. In conclusion, this study provides valuable information for further research on the regulatory mechanisms of WRKY TFs in plant growth, secondary metabolism, and resistance to various stressors.

6 Abbreviations

TFs, Transcription factors; P. grandifloras, Platycodon grandiflorus; MeJA, methyl jasmonate; MW, molecular weight; pI, isoelectric point; AA, amino acid; ORF, open reading frame; h, hours; qRT-PCR, quantitative real-time polymerase chain reaction; CK, control check; RNA-seq, RNA sequencing; MEME, multiple expectation maximizations for motif elicitation; RSEM, RNA-seq by expectation-maximization; GO, gene ontology; PPI, protein-protein interaction; BP, biological processes; CC, cellular components; MF, molecular functions; MKS1, Meckel syndrome, type 1; SIB1, sigma factor binding protein 1; MPK, mitogen-activated protein kinase; GUN5, genomes uncoupled 5; MVA, mevalonate; MEP, 2-C-methyl-D-erythritol-4-phosphate UHPLC-Orbitrap-MS/MS, ultra-high performance liquid chromatography-Orbitrap-tandem mass spectrometry.

Data availability

The raw data supporting the conclusions of this article will be made available by the authors without undue reservation to any qualified researcher.

Author contributions

LZ and SG conceived the study. BC and ND prepared and cultivated the plants. ML and JL prepared the samples, extracted RNA, and performed the qRT-PCR experiments. XC, JW, and HY assisted with data analysis and interpretation. JL and HY wrote the manuscript, and LZ revised the manuscript. All authors have reviewed, read, and approved the final manuscript.

Acknowledgments

We would like to thank BGI Genomics for their assistance with the experiments.

Disclosure statement

No potential conflict of interest was reported by the author(s).

Supplementary material

Supplemental data for this article can be accessed online at https://doi.org/10.1080/15592324.2022.2089473

Additional information

Funding

This research was funded by the National Natural Science Foundation of China (No. U21A20406) and the Exploration Project of Anhui University of Chinese Medicine (2021zxts52).