Comparative transcriptome analyses revealed differential strategies of roots and leaves from methyl jasmonate treatment Baphicacanthus cusia (Nees) Bremek and differentially expressed genes involved in tryptophan biosynthesis

Baphicacanthus cusia (Nees) Bremek (B. cusia) is an effective herb for the treatment of acute promyelocytic leukemia and psoriasis in traditional Chinese medicine. Methyl jasmonate (MeJA) is a well-known signaling phytohormone that triggers gene expression in secondary metabolism. Currently, MeJA-mediated biosynthesis of indigo and indirubin in B. cusia is not well understood. In this study, we analyzed the content of indigo and indirubin in leaf and root tissues of B. cusia with high-performance liquid chromatography and measured photosynthetic characteristics of leaves treated by MeJA using FluorCam6 Fluorometer and chlorophyll fluorescence using the portable photosynthesis system CIRAS-2. We performed de novo RNA-seq of B. cusia leaf and root transcriptional profiles to investigate differentially expressed genes (DEGs) in response to exogenous MeJA application. The amount of indigo in MeJA-treated leaves were higher than that in controled leaves (p = 0.004), and the amounts of indigo in treated roots was higher than that in controlled roots (p = 0.048); Chlorophyll fluorescence of leaves treated with MeJA were significantly decreased. Leaves treated with MeJA showed lower photosynthetic rate compared to the control in the absence of MeJA. Functional annotation of DEGs showed the DEGs related to growth and development processes were down-regulated in the treated leaves, while most of the unigenes involved in the defense response were up-regulated in treated roots. This coincided with the effects of MeJA on photosynthetic characteristics and chlorophyll fluorescence. The qRT-PCR results showed that MeJA appears to down-regulate the gene expression of tryptophan synthase β-subunits (trpA-β) in leaves but increased the gene expression of anthranilate synthase (trp 3) in roots responsible for increased indigo content. The results showed that MeJA suppressed leaf photosynthesis for B. cusia and this growth-defense trade-off may contribute to the improved adaptability of B. cusia in changing environments.

Baphicacanthus cusia (Nees) Bremek (B. cusia) is generally distributed in southern China, Bangladesh, northeast India, Myanmar, Himalayan and the mid-south Peninsula [24]. The root and aerial parts of B. cusia are used as medicinal materials in Nan-Ban-Lan-Gen [25] and Indigo Naturalis [26], respectively. These were widely used as traditional Chinese medicine to remove heat from blood and eliminate toxicity in the human body [27]. Pharmacological studies have shown that Nan-Ban-Lan-Gen has many biological activities, such as antibacterial [28], antiviral [29,30], immunomodulatory [31,32] and anti-inflammatory activities [33]. Previous clinical studies indicated that Indigo Naturalis is good for the treatment of acute promyelocytic leukemia [34,35], ulcerative colitis [36,37], and psoriatic lesions [38]. And that the secondary metabolites, such as indirubin, indigo and tryptanthrin were the active components [39]. The molecular mechanism of the production of the active components B. cusia in response to biotic or abiotic stresses has not been reported.
In our previous study [40], tryptophan synthase was confirmed to be the candidate gene involved in biosynthesis of indican, which was one of the genes in the tryptophan biosynthesis pathway. Hence, we speculate that the key genes affecting the biosynthesis of indigo and indirubin are the genes involved in the tryptophan biosynthesis pathway, the upstream pathway for the biosynthesis of indican.
In this study, to obtain in-depth knowledge of indican biosynthesis upstream gene expression changes in MeJA-treated leaves and roots, we performed de novo high-throughput sequencing of B. cusia leaves and roots before and after MeJA treatment. The assembled unigenes were annotated by five databases: nr, SwissProt, GO, COG and KEGG. We focused on the differentially expressed genes (DEGs) in the MeJA-treated B. cusia leaves and roots. Furthermore, we identified several candidate genes associated with indican biosynthesis via the upstream tryptophan pathway by qRT-PCR. Meanwhile, we determined the content of indigo and indirubin in leaf and root tissues of B. cusia and measured the photosynthetic characteristics and chlorophyll fluorescence of leaves treated by MeJA. This is the first report on the transcriptional response of B. cusia leaves and roots treated by MeJA. The molecular mechanisms underlying MeJA treatment will promote research on the biological mechanisms involved in molecular breeding and secondary metabolite regulation of B. cusia. The transcriptome may help to clarify differentiated strategies of roots and leaves in response to exogenous application of methyl jasmonate in B. cusia.

DEGs in response to MeJA
In response to MeJA treatment, 8,355 DEGs were found to be significantly differentially expressed in 33,317 annotated unigenes, among which 2,664 DEGs were up-regulated and 3,335 DEGs were down-regulated in MeJA-treated roots and 761 DEGs were up-regulated and 1,595 DEGs were down-regulated in MeJA-treated leaves (Fig 1). These results suggest that leaves of B. cusia demonstrate greater suppression of unigenes than activation in response to MeJA treatment; in contrast, roots of B. cusia demonstrate greater activation than suppression of unigenes in response to MeJA treament.

KEGG pathway classification of unigenes
To identify the specific biological pathways of unigenes assembled above, the KEGG pathway database was employed to characterize the functional classification and pathway mapping according to sequence homology. Overall, 744 of 11,418 unigenes were classified into five main categories and 107 KEGG pathways, which included cellular processes, environmental information processing, genetic information processing, metabolism and organismal systems, in leaves and 1,528 of 11,980 unigenes were classified into the same five main categories and 121 KEGG pathways in roots. As shown in Table 1, in CL-VS-TL (S3 Table), the carbohydrate metabolism has the highest number of unigenes (150), followed by amino acid metabolism (95), biosynthesis of other secondary metabolites (79), lipid metabolism (66), and global and overview (63). In CR-VS-TR (S4 Table), the carbohydrate metabolism category also has the highest number of unigenes (254), but is followed by translation (161) of genetic information processing, global and overview (150), amino acid metabolism (124), lipid metabolism (121), and energy metabolism (108).

Validation of RNA-seq analysis by qRT-PCR
To verify the relative expression levels of DEGs involved in the biosynthesis of phenylalanine, tyrosine and tryptophan obtained by RNA sequencing, we carried out qRT-PCR on eight Comparative transcriptome analyses of B. cusia treated by MeJA and DEGs involved in tryptophan biosynthesis relative unigenes (aroF, aroK, aroA, aroC, TRP3, trpD, trpA-α and trpA-β) involved in the biosynthesis of tryptophan.
The results of qRT-PCR analysis revealed that the relative expression of candidate unigenes, which were down-regulated in CL-VS-TL, were consistent with the data of RNA-seq RPKM; in CR-VS-TR, the relative expression of all candidate unigenes was consistent with data of the RNA-seq RPKM except for trpA-β and aroC, the relative expression of trpA-β and aroC in CR-VS-TR were the opposite of the RNA-seq data (Fig 4), the relative expression of trpA-α have not detected.

Quantitative analysis of indigo and indirubin
The calibration curves for indigo and indirubin were prepared using five different concentration mix reference materials. The regression equations and correlation coefficients (r 2 ) were: Y = 7.1499X + 3.6795 (linear range from 13.92 to 125.25μg�L -1 ), r 2 = 0.9996 for indigo; Y = 16.342X + 2.7178 (linear range from 9.53 to 85.79μg�L -1 ), r 2 = 0.9997 for indirubin. The relative amounts (μg/g) of indigo and indirubin in the leaf and root tissues were calculated using the above equations. The amounts of indigo in TL were higher than those in CL (p = 0.004, Fig 5), and the amounts of indigo in TR were higher than those in CR (p = 0.048); there were no differences in the amounts of indirubin between CL and TL (p = 0.273), as well as between CR and TR (p = 0.904).

Influence of MeJA on appearance and traits of the leaves and roots
In  Table), the data were downloaded from http:// rp5.ru/archive.php?wmo_id=58847&lang=cn. The leaf area is calculated by the leaf length multiplied by 0.7 times the leaf width. The row data of length and width of B. cusia leaves was shown in supplementary file (S7 Table). The results indicated that the leaves showed different degrees of shrinkage after treated with methyl jasmonate (Fig 6, S6 Fig)

Influence of MeJA on photosynthetic parameters and chlorophyll fluorescence
MeJA stress severely affected gas exchange parameters and chlorophyll fluorescence compared to the control (Fig 7).

Discussion
From Fig 3, we can see that Erythrose 4-phosphate is a precursor in tryptophan biosynthesis, which is an intermediate in the pentose phosphate pathway and the Calvin cycle. The Calvin cycle is a series of chemical reactions occurring in the chloroplast during photosynthesis. Photosynthesis affects the production of Erythrose 4-phosphate, which in turn affects the biosynthesis of tryptophan. Therefore, this study measured photosynthesis to demonstrate that exogenous methyl jasmonate may affect the synthesis of tryptophan through photosynthesis. In this study, chlorophyll fluorescence measurements show that leaves' photosynthesis was reduced in the leaves, and that photosynthesis was suppressed after treatment with 22.29 μM MeJA, and DEGs of CL-VS-TL involved in "photosystem", "photosynthetic membrane", "photosynthesis-antenna proteins" were down-regulated, which further verified that MeJA induce leaf senescence in B. cusia. The appearance and traits of the leaves and roots treated with different concentrations of MeJA showed that the leaf areas decreased and the numbers of hair roots increased. These results were consistent with previous reports about the effects of MeJA on leaf senescence [41][42][43] and lateral roots formation [44][45][46][47]. The effects of different concentrations of methyl jasmonate on Baphicacanthus cusia (Nees) Bremek in different ecological environments were further studied in the future. RNA-seq has been applied to gene expression analysis on the genome-wide level in a number of plants [48][49][50]. In a previous study [51], qRT-PCR assays were employed to analyze the expression stability of ten candidate reference genes and the expression levels of two genes involved in synthesis of terpenoid indole alkaloids of B. cusia after 24 h MeJA treatment; However, the qRT-PCR assay is limited by the number of sequencing fragments. The RNA-seq method allows for comprehensive and accurate quantitative information on gene expression to acquire better gene regulatory profiles and to identify more DEGs. In this study, we obtained 33,317 annotated unigenes and analyzed 8,355 DEG expression levels, especially genes related to tryptophan synthesis. The results revealed that the relative expression of all the DEGs involved in tryptophan synthesis were down-regulated in leaves treated by MeJA after seven days. However, the relative expression levels of aroF, aroK, trpD and trp3 were up-regulated in roots, which is inconsistent with the results of the previous study [51].
After treatment with methyl jasmonate, the indigo content of the leaves and roots from treated B. cusia increased. This was consistent with the results of the previous study [51]. However, expression levels of the eight related genes involved in phenylalanine, tyrosine and tryptophan biosynthesis were down-regulated in leaves and gene expression levels of aroA were down-regulated, while the other seven genes were up-regulated in roots of B. cusia treated by MeJA. This indicates that exogenous MeJA regulates the tryptophan biosynthetic pathway, which further affects the expression of downstream genes involved in indican biosynthesis. Genes involve in the indican metabolism pathway, which is upstream of indigo and indirubin biosynthesis, may also be induced by MeJA treatment; such genes include cytochrome P450, UDP-glycosyltransferase, glucosidase and tryptophan synthase [40]. Our data provide a valuable resource for discovering candidate genes related to indigo, indirubin and indican biosynthesis in response to MeJA, especially TrpA. TrpA is a heterodimeric enzyme with two α and β subunits. The α-subunit catalyzes indole production and the β-subunit catalyzes tryptophan yield [52,53]. In this study, expression levels of TrpA β-subunits were down-regulated in leaves, which would increase indole biosynthesis. Hence, DEGs of trpA β-subunit may be the key upstream genes leading to the production of indigo in leaves of B. cusia. However, unlike leaves, DEGs of trp3 may be the key upstream genes leading to the production of indigo in roots of B. cusia.
The optimal defense theory (ODT) and the growth-differentiation balance hypothesis (GDBH) were two excellent theories that attempted to explain the expression patterns of chemical defense in plants. The ODT proposed that high fitness value plant parts would be highly defended, but the GDBH speculated that slow-growing plant parts would be highly defended [54]. Tryptophan acts as a biochemical precursor of auxin, which plays a key role in the plant life cycle and development [55]. In this study, we found that DEGs related to growth and development processes were down-regulated in the treated leaves, such as "photosystem", "photosynthesis-antenna proteins", catalytic activity and "negative regulation of biological process", while most of the unigenes involved in the defense response were up-regulated in treated roots, such as "catalytic activity", "carbohydrate metabolic process", "sesquiterpenoid and triterpenoid biosynthesis", and oxidoreductase activity. Similar findings from GO and KEGG functional annotation support the view of JAs playing a key role in regulating resource distribution between the defense and growth competition processes [56]. These results support established theories of ODT [57,58]. This growth-defense trade-off may help B. cusia improve its adaptability by shifting the energy conservation of down-regulated photosynthesis to a relevant defense response in a changing environment.

Plant materials and MeJA treatment
The Baphicacanthus cusia (Nees) Bremek (B. cusia) samples were collected from our experimental field at the Fujian Agriculture and Forestry University (26.0822, 119.2398). The B. cusia samples were propagated from cuttings and planted in our experimental field. When shoots were rooted and well-established, seedlings were selected from the experimental field and planted into shade plots. All leaves of treated groups (TL: Treated Leaf and TR: Treated Root) were sprayed with a solution of 0.01% (v/v) Tween 20 containing 22.29 μM MeJA (Sigma-Aldrich), and the leaves of control groups (CL: Controlled Leaf and CR: Controlled Root) were sprayed with 0.01% (v/v) Tween 20 without MeJA to the point of runoff. After treatment, the second to fourth leaves from the top of plants were harvested for physicochemical and molecular analysis. The collected leaves and roots were frozen in liquid nitrogen, and then stored at -80˚C for future analyses. Three biological replicates from the independent control and treated B. cusia were prepared for RNA sequencing. Each biological replicate used three root or leaf plants. qRT-PCR was performed using three biological and three technical replicates.

RNA isolation, sequencing, assembly and bioinformatics analysis
Total RNA was isolated with the EASYspin plant RNA kit (Aidlab Ltd, Beijing, China). RNA purity was validated using a Biodrop spectrophotometer (Biochrom Ltd). RNA integrity was analyzed using Agilent Bioanalyzer 2100 (Agilent Technologies). The construction of the RNA library and RNA sequencing of the libraries were performed by commercial service providers Gene denovo Biotechnology Co. (Guangzhou, China) under the Illumina HiSeq 4000. The RNA-seq data of the treated-by-MeJA and control B. cusia were deposited in the NCBI SRA repository: SRA628524 (SRP124081: PRJNA415260). The raw reads were processed and analyzed using the method previously reported [40]. The assembled unigenes were BLASTX searched and annotated against the databases of NCBI non-redundant protein (NR), Swis-sProt, euKaryotic Orthologous Groups (KOG), Gene Ontology (GO), and Kyoto Encyclopedia of Genes and Genomes (KEGG) with an E-value threshold of 1E-5. The gene abundances were calculated and normalized to reads per kb per million reads (RPKM). The significant DEGs were identified with the threshold of fold change >2 and false discovery rate (FDR) <0.05. The calculated p-value was obtained by FDR correction, and the q-value is the multi-hypothesis test corrected p-value.

Validation of qRT-PCR
Expression levels of the eight candidate genes, aroF, aroK, aroA, aroC, trp3, trpD, trpA-α and trpA-β were evaluated; these genes are involved in phenylalanine, tyrosine and trytophan biosynthesis. Primers of the eight candidate genes and glyceraldehyde-3-phosphate dehydrogenase (GAPDH) used for the reference gene were designed to amplify short regions using primer 3 web (version 4.1.0, http://bioinfo.ut.ee/primer3/, S5 Table). RNA was reverse transcribed using the Fast Quantity RT Kit (TIANGEN, China) according to the manufacturer's specifications. qRT-PCR was carried out by using the SYBR Premix Ex Taq kit (TaKaRa Bio Inc). The amplification was executed with the following PCR program: 90 s at 95˚C, 40 cycles of 5 s at 95˚C for denaturation, 15 s at 60˚C for annealing, 20 s at 72˚C for elongation, and 65˚C~90˚C for melting curve analysis. qRT-PCR was performed on the ABI PRISM 7500 Real-Time PCR System (Applied Biosystems, US). The relative expression ratios of each candidate gene were calculated based on the comparative CT (2 −ΔΔCt ) method.

Analysis with high-performance liquid chromatography
The HPLC analysis was carried out using an LC-20AT HPLC system (Shimadzu, Japan). Methanol and water, with a volume ratio of 70: 30 was used as the mobile phase for both indigo and indirubin, and the flow rate was 1.0 ml/min. To prepare the solutions, samples of the leaves and roots were extracted by a Soxhlet extractor (frequency of 40 KHz, power of 500 W). The filtrates were combined and swept in a rotary evaporator. The mixed reference solutions of indigo and indirubin were dissolved in N, N-dimethyl formamide in triplicate. After filtrated with a 0.45 μm filter, 10 μl of sample solutions were injected into HPLC in duplicate.

Measurement of photosynthetic characteristics and chlorophyll fluorescence
The chlorophyll fluorescence of minimal fluorescence (F0), maximal fluorescence (Fm), variable fluorescence (Fv), Fv/F0 and Fv/Fm were measured in fully expanded leaves of B. cusia using FluorCam6 Fluorometer (Photon Systems Instruments, Czech Republic). The gas exchange parameters of stomatal conductance (GS), intercellular CO 2 concentration (CI), transpiration rate (EVPA), net photosynthetic rate (PN), and photosynthetically active radiation (PAR) were measured using the portable photosynthesis system CIRAS-2 (Hansatech, UK). All determinations were repeated in triplicates. The measurements were carried out at 10 am in sunny weather.
Supporting information S1