In vitro grown B. rapa ssp. pekinensis (Chinese cabbage) variety Seoul (Kyoungshin seed, Korea) seedlings were transplanted into pots (25 × 25 cm) containing a sterile mixture of soil, sand, and Biopeat (Seminis Asia Ltd., Korea). The plants were then raised under controlled glasshouse conditions until they reached maturity. Plant samples were collected from seven different developmental stages, such as five-day-old seedling (SDL), one-month-old root (R), shoot (S), and leaf (L), two-month-old roots (2R), inner leaves (IL), and outer leaves (OL). All samples were collected and stored at –80 °C, immediately after sampling, for further analysis. For stress treatments, surface sterilized B. rapa seeds were grown on half strength Murashige and Skoog (MS Duchefa, Germany) medium with 0.8% agar and 0.5% sucrose for 10 days at 21 °C, under a 16/8-h photoperiod, and then the plantlets were shifted to a liquid MS medium supplemented with 0.5% sucrose. Different stress treatments, such as 100 μM of indole-3-acetic acid (IAA), gibberellic acid (GA), methyl jasmonate (MeJA), salycilic acid (SA), abscicic acid (ABA), and 3% glucose (MS basal medium devoid of sucrose) were given to seedlings and incubated for 24 h. MeJA and SA treatments were given to both shoots and roots simultaneously, whereas other stresses were applied through the root system. Low- and high-temperature treatments were given by incubating the plantlets at 4 °C for 24 h and 37 °C for 2 h. Wound stress was induced by cutting the leaves with a scalpel, and the samples were collected after five minutes of wounding. Root and shoot tissue samples were collected from all stress treatments, with the exception of high-temperature treatments, in which the whole plant was taken for quantitative real-time PCR (qRT-PCR) gene expression analysis.

All nucleotide and protein sequences of MYB and other genes used in this study were extracted from the Brassica database (BRAD) (http://brassicadb.org/brad), the National Centre for Biotechnology Information (NCBI) sequence database (http://www.ncbi.nlm.nih.gov/), and the Arabidopsis Information Resource (TAIR) Arabidopsis genome database (http://www.arabidopsis.org/) [11,12]. The putative orthologous MYB protein sequences from other crops were identified through the NCBI BLASTP program, using B. rapa ssp. pekinensis MYB28 and 29 amino acid sequences [11,12] as a query sequence, with a cut-off E-value of 1^–10. The top 5–10 orthologous MYB proteins selected from each homology search were used for multiple sequence alignments and phylogenetic analysis. The sequence identity between the MYB proteins was analyzed with the ClustalW2 (http://www.ebi.ac.uk/Tools/msa/clustalw2) program, by utilizing the default parameter [13]. The phylogenetic tree was built by the neighbour-joining method, using the Molecular Evolution Genetics Analysis (MEGA 4.02) program [14] to analyze the evolutionary relationships among the MYB proteins. The reliability of the tree was achieved with 1000 bootstrap replicates.

The putative structural domain identification in the B. rapa MYB28 and 29 paralogous proteins and Arabidopsis MYB28 and 29 proteins were analyzed by using the PROSITE (http://www.expasy.org/tools/scanprosite/) program [15]. For putative cis regulatory motif analysis, the promoter regions (–1 to –1.5 kb) from gene sequences of all the BrMYB28 and BrMYB29 paralogs and their DSGs, namely, methylthioalkylmalate synthase 1 (MAM1), MAM3, cytochrome P450 belonging to the CYP79 family gene (CYP79F1), CYP83 family gene (CYP83A1), C-S lyase superroot-1 (SUR1), sulfotransferase 5C (ST5C), flavin monooxygenase glucosinolate S-oxygenase 2 (FMO-GSOX2), and flavin monooxygenase glucosinolate S-oxygenase 5 (FMO-GSOX5) in the AGSL biosynthetic pathway, were extracted from the BRAD database. The distributions of putative cis regulatory motifs were identified with the help of the PLACE (http://www.dna.affrc.go.jp/PLACE/signalscan.html) and PlantPAN (http://plantpan.mbc.nctu.edu.tw/) software programs [16,17]. Apart from the common letters, such as A, G, C or T designated for nucleotides, the other letters used in the PLACE database were: K: G or T; N: A, C, G or T; W: A or T; Y: C or T. The flow chart of the work presented in the materials and method section was represented in Fig. 1.

Total RNA was isolated from seven different developmental stages (SDL, R, S, L, 2R, IL, and OL) and from various stress-treated B. rapa plants, using TRI reagent (Sigma–Aldrich, Saint-Louis, USA), followed by treatment with DNase I (Promega, MA, USA). cDNA synthesis was performed with 3 μg of total RNA, oligodT primer, and Super Script-II RT, according to the manufacturer's instructions (Invitrogen, Carlsbad, CA, USA). The B. rapa actin (BrAct – FJ969844) gene was used as an internal control for relative quantification of BrMYBs and their DSGs transcripts. The nucleotide multiple sequence similarities between the BrMYB amplicons were analyzed by using ClustalW2. qRT-PCR was performed in the CFX 96 Real-Time PCR system (Bio-Rad, Hercules CA, USA) using the SYBR PCR kit (Bio-Rad, USA). The PCR amplification reactions were performed as follows: 95 °C for 5 min, followed by 35 cycles at 94 °C for 30 s, annealing temperature for 30 s and an extension at 72 °C for 2 min using the primers listed in Table S1. The specificity of amplification was confirmed by the melt curve analysis and agarose gel electrophoresis. The relative gene expression was quantified by 2^−ΔΔCt method [18]. Comparative threshold (Ct) values were normalized to actin control and compared to quantify the relative transcript levels.

Each experiment was carried out independently, at least thrice, with three replicates. The resulting data from different experiments were analyzed using Origin 8.1 (Origin lab corporation, Northampton MA USA) software program. Statistical analyses were done using one-way analysis of variance (ANOVA) followed by Tukey's multiple comparison tests, and the standard error (SE) was calculated using the n values of each experiment.

Identification of putative structural domains in BrMYBs and Arabidopsis MYB28 and MYB29 proteins, using the PROSITE tool (http://www.expasy.org/tools/scanprosite/), indicated the presence of a specific DNA-binding domain (HTH) in all the BrMYBs, except in BrMYB29.3 (Fig. 2). In addition to the DNA-binding domain (DBD), the BrMYB28.3 paralog had a putative transmembrane domain between the amino acid region 332 and 354. BrMYB29.2 possessed only a single HTH motif along with an MYB-like domain and did not contain the regulatory protein–protein interaction domain (Fig. 2).

The multiple sequence alignment results revealed that B. rapa MYB28 paralogs, such as BrMYB28.1, 28.2, and 28.3 showed 83, 80, and 73% sequence identity, respectively, with AtMYB28 (AAS10113), at the protein level. Similarly, Brassica MYB29 paralogs, such as BrMYB29.1, 29.2, and 29.3 amino acid sequences showed 73, 98, and 72% sequence identity, respectively, with AtMYB29 (AAS10087) (Fig. S1). Furthermore, in the phylogenetic analysis, all the IGSL and AGSL regulators had formed distinct clusters (Fig. S2).

The analysis of the putative regulatory motifs from the promoter region of the DSGs from the key steps of the AGSL biosynthetic pathway, such as MAM1, MAM3 (chain elongation), CYP79F1, CYP83A1, SUR1, ST5C (core structure formation), FMO-GSOX2, and FMO-GSOX5 (secondary side chain modification) was carried out to analyze the motifs related to the binding sites of the MYB TF genes, and also the stress-specific and tissue-specific gene expressions. Similarly, putative regulatory motifs from the promoter regions of Brassica MYB TF genes (BrMYB28 and BrMYB29 paralogs) and Arabidopsis MYB28, 29 genes were also analyzed, to find the motifs related to the stress-specific and tissue-specific gene expressions. The proximal promoters (–101 to –1000) function as a docking site for TFs, and therefore, we selected the promoter region covering –1 to –1500 bp from the translational start site of the DSGs and all paralogs of the BrMYB28, 29 genes as well as the Arabidopsis MYB28, 29 genes for in silico cis motif analysis. The results showed the presence of highly conserved regulatory elements, such as, CCAAT, CGCG, TATA boxes, sulfur-response motifs, binding sites for MYB, as well as motifs related to abiotic and biotic stresses, hormones, glucose signalling, circadian, and organ-specific. Duplication frequency of putative cis regulatory motifs and their functions in the UPS region of BrMYB28, 29 paralogs and their DSGs as well as the Arabidopsis MYB28, 29 genes are presented in the supplementary table (Tables S1 and S2). The respective positions of the important putative cis motif clusters in the 5′ UPS region of BrMYB28, 29 paralogs and their DSGs are shown in the figure (Fig. S3a and S3b).

In order to monitor the quantitative gene expression variation of BrMYB28 and BrMYB29 paralogous genes, qRT-PCR was performed at different developmental stages of B. rapa. The sequence similarities between the amplicons of BrMYBs were shown in Fig. S4. The gene expression data showed a profound variation in tissues of different developmental stages, such as seedlings, one-month-old root, shoot, and leaf, and two-month-old roots, inner leaves and outer leaves. BrMYB28.3 gene expression was ubiquitous and exceptionally high in all the stages compared to other closely related BrMYB genes, such as BrMYB28.1 and BrMYB28.2 genes (Fig. 3a). Among the different developmental stages, gene expression of all three paralogs of BrMYB28 was exclusively higher in the two-month-old outer leaf tissues. Similarly, among the BrMYB29 paralogs (BrMYB29.1 and BrMYB29.3) studied, BrMYB29.1 showed high transcript abundance in all the developmental stages compared to BrMYB29.3 (Fig. 3b). Both BrMYB29 paralogs were preferentially expressed in two-month-old root samples as well as in two-month-old outer leaf tissues.

The gene expression variations of B. rapa AGSL biosynthetic pathway genes (paralogs that are closely related to the Arabidopsis AGSL biosynthetic genes based on sequence homology) was analyzed at different developmental stages of B. rapa. The gene expression profile for the following AGSL biosynthetic genes, such as MAM1, MAM3, CYP79F1, CYP83A1, SUR1, ST5C, FMO-GSOX2, and FMO-GSOX5 was analyzed in all the selected developmental stages (SDL, 1R, S, L, 2R, IL and OL) through qRT-PCR (Fig. 4). Among the genes studied, ST5C, CYP83A1, and MAM3 showed a ubiquitous and strong expression, irrespective of the developmental stages. However, the expressions of MAM1, CYP79F1, SUR1, FMO-GSOX2 and FMO-GSOX5 genes revealed to be of low level (Fig. 4). Most of the AGSL biosynthetic pathway genes were also found to be highly expressed in two-month-old outer leaf tissues as well as in the seedling stages.

The stress-specific gene expression variations of BrMYBs were analyzed under various stress treatments (Fig. 5a and b). Among the various stresses given to plants, a strong induction in most of the BrMYBs was observed in the wounding, IAA, and low- and high-temperature treatments. To investigate the possible involvement of MYB TF genes in phytohormone-mediated signaling in B. rapa ssp. pekinensis, we examined the expression of BrMYB genes in response to IAA, GA, ABA, MeJA and SA treatments. The phytohormone IAA strongly induced the upregulation of all BrMYBs in the aerial parts, whereas the expression of all genes was found to be the least in the root tissues. In contrast, GA triggered the expression of BrMYB29.1 in root tissues and BrMyb28.3 in shoot tissues (Fig. 5a). In the aerial parts of Brassica, MeJA induced the upregulation of BrMYB29.3 and BrMYB28.1 genes, whereas in the roots it induced the upregulation of BrMYB28.3 and BrMYB29.1 genes. In response to SA treatment, the expression of most of the BrMYBs increased significantly in roots except BrMYB29.3 gene. However, there were some differences in the ABA-induced expression levels of BrMYBs. The ABA-induced expression of BrMYB28.2 and BrMYB29.1 was observed in shoot and root tissues, respectively, and the expression of BrMYB28.1 was least expressed among the BrMYBs, in response to the ABA treatment (Fig. 5a).

Similarly, the other abiotic stresses (glucose, high- and low-temperature treatments) and biotic stress (wound stress) specific gene expressions of BrMYBs in B. rapa were analyzed by qRT-PCR. Glucose treatment slightly induced upregulation of the BrMYB28.2 paralog in both root and shoot tissues, while it slightly induced the BrMYB29.3 paralog in the shoot tissues. Low-temperature treatment induced the expression of most of the BrMYBs in the root tissues, whereas the expression of BrMYB28.1 was the lowest in this treatment. Among the BrMYB paralogs, BrMYB29.1 and BrMYB28.3 genes showed the strongest expression, whereas the BrMYB28.2 gene revealed the least expression in plants treated at high temperature (Fig. 5b). BrMYB paralogs, such as BrMYB28.2, 28.3, and BrMYB29.1 genes were predominantly expressed in wound stressed plants. In contrast, BrMYB28.1 showed least expression in shoots as well as in root tissues of wound stressed plants (Fig. 5b).