Abstract

   Jasmonates, including jasmonic acid (JA) and its derivatives such as
   methyl jasmonate (MeJA) or jasmonly isoleucine (JA-Ile), regulate plant
   responses to various biotic and abiotic stresses. In this study, we
   applied exogenous MeJA onto Senna tora leaves subjected to wounding and
   conducted a transcriptome deep sequencing analysis at 1 (T1), 3 (T3), 6
   (T6), and 24 (T24) h after MeJA induction, along with the pretreatment
   control at 0 h (T0). Out of 18,883 mapped genes, we identified 10,048
   differentially expressed genes (DEGs) between the T0 time point and at
   least one of the four treatment times. We detected the most DEGs at T3,
   followed by T6, T1, and T24. We observed the upregulation of genes
   related to JA biosynthesis upon exogenous MeJA application. Similarly,
   transcript levels of genes related to flavonoid biosynthesis increased
   after MeJA application and tended to reach their maximum at T6. In
   agreement, the flavonols kaempferol and quercetin reached their highest
   accumulation at T24, whereas the levels of the anthraquinones
   aloe-emodin, emodin, and citreorosein remained constant until T24. This
   study highlights an increase in flavonoid biosynthesis following both
   MeJA application and mechanical wounding, whereas no significant
   influence is observed on anthraquinone biosynthesis. These results
   provide insights into the distinct regulatory pathways of flavonoid and
   anthraquinone biosynthesis in response to MeJA and mechanical wounding.

   Keywords: Senna tora, transcriptomics, flavonoid, anthraquinone,
   chalcone synthase, wounding, methyl jasmonate

1. Introduction

   Plants are constantly faced with biotic and abiotic stresses and
   developed defense systems to respond to stresses [[36]1]. Secondary
   metabolites are one of the means by which plants survive environmental
   stresses [[37]2,[38]3]. Plant hormones regulate the biosynthesis of
   numerous secondary metabolites [[39]4,[40]5,[41]6]. Plant hormones are
   divided into several sub-classes according to their chemical structures
   and comprise auxins, gibberellins, ethylene, cytokinins, abscisic acid,
   jasmonic acid (JA), salicylic acid, brassinosteroids, and
   strigolactones [[42]7].

   Jasmonates are phytohormones whose accumulation responds to many biotic
   and abiotic stresses. They play an essential role in plant development,
   stress responses, and secondary metabolite regulation as essential
   signaling molecules. Jasmonates include all intermediates and products
   of the JA biosynthesis pathway and its metabolites. In addition to JA,
   the most studied jasmonates are methyl jasmonate (MeJA) and
   jasmonyl-isoleucine (JA-Ile), which are formed by jasmonic acid
   O-methyl transferase (JMT) and by JASMONATE RESISTANT 1 (JAR1),
   respectively, in the cytoplasm. MeJA, a methyl ester of JA, is volatile
   and can easily diffuse into the air and mediates plant-to-plant
   communication in response to biotic and abiotic stresses without having
   to rely on transmission through the vasculature [[43]8,[44]9]. JA-Ile
   is a biologically active isoleucine conjugate that regulates many
   stress responses and developmental responses. JA and MeJA are
   particularly important in defense responses to necrotrophic pathogens
   and wounding by insects, as well as diverse stages of plant development
   [[45]10,[46]11,[47]12,[48]13,[49]14,[50]15]. Collectively, jasmonates,
   whether produced within plants or applied exogenously, regulate the
   expression of many genes involved in plants development and responses
   to various stresses.

   Senna tora (L.) Roxb. is an annual shrub belonging to the family
   Fabaceae (subfamily Caesalpinioideae). It is also known as Cassia tora
   and generally distributed throughout India, Pakistan, Sri Lanka, and
   China [[51]16,[52]17,[53]18] and is used as an herbal medicine in these
   and nearby countries. This plant offers medicinal benefits on account
   of its anti-diabetes [[54]19], anti-dermatitis [[55]20,[56]21],
   anti-constipation [[57]22], and anti-bronchitis properties. These
   various effects are due to numerous active secondary metabolites that
   are anti-inflammatory, antibacterial, or antioxidant
   [[58]17,[59]23,[60]24,[61]25]. Such diverse medicinal ingredients are
   distributed in various parts of this plant, including the leaves,
   roots, and seeds [[62]26,[63]27,[64]28]. Among various metabolites,
   flavonoids tend to be found in relatively larger amounts in leaves and
   anthraquinones in seeds [[65]28,[66]29]. Both the flavonoids and
   anthraquinones are polyphenolic compounds, and their medicinal effects
   are mostly attributed to their antioxidant property [[67]30,[68]31].
   They modulate reactive oxygen species (ROS) levels as radical
   scavengers and protect plants from harmful UV irradiation as a UV
   protectant [[69]32]. In particular, kaempferol and quercetin, two
   flavonols, attracted much research attention because they are
   medicinally useful [[70]33,[71]34,[72]35].

   Previously, we showed that mechanical wounding accelerated the
   kaempferol and quercetin accumulation in S. tora leaves [[73]36]. In
   this study, we comprehensively analyzed the global expression patterns
   of genes involved in metabolism and the specific secondary metabolite
   production (flavonoid and anthraquinones) under the combined effects of
   MeJA and wounding. This study provides valuable insights into the
   regulatory pathways of flavonoid and anthraquinone biosynthesis in
   response to both biotic and abiotic stresses.

2. Results and Discussion

2.1. Identification of DEGs in Response to MeJA Treatment and Wounding

   The sequencing libraries derived from total RNA extracted from S. tora
   leaves before wounding (T0) and after 1 (T1), 3 (T3), 6 (T6), or 24
   (T24) h of MeJA application produced 33,502,278–46,846,336 trimmed
   reads—after processing raw reads. The clean reads encompassed between
   3,358,299,940 and 4,690,843,228 bases with a percentage of reads with a
   minimum quality score of Q30 of 96.1% to 97.1%. We mapped the clean
   reads to the S. tora reference genome [[74]37], returning expression
   information for 18,883 genes (mapping ratio of 95.5–98.0%,
   [75]Supplementary Table S1) from the three replicates for each time
   point. Thus, we obtained expression data for about 41.7% of all
   protein-coding genes annotated in the genome (out of 45,268 genes).

   To identify differentially expressed genes (DEGs) relative to the T0
   time point, we calculated the fold change in fragments per kilobase of
   transcript per million (FPKM) values for each sample over that of T0.
   Genes fulfilling the criteria of |Log[2](fold change FPKM)| > 1 and
   adjusted p-value < 0.05 were defined as DEGs. We obtained 10,048 DEGs
   out of the 18,883 mapped genes across all pairwise comparisons.
   Specifically, the numbers of DEGs identified at the T1, T3, T6, and T24
   time points were 4477, 7728, 6191, and 2460, respectively ([76]Figure
   1A). We detected the highest number of DEGs at T3, followed by T6, T1,
   and T24. This pattern indicates that the influence of MeJA application
   and wounding reaches its maximum at 3 h after the onset of MeJA
   treatment and wounding. We observed the same pattern for upregulated
   and downregulated DEGs.

Figure 1.

   [77]Figure 1
   [78]Open in a new tab

   Differential gene expression analysis among time points after methyl
   jasmonate and wounding treatments in S. tora. (A) Number of upregulated
   and downregulated DEGs at each time point relative to T0 (before
   wounding and spraying) samples. Venn diagrams showing the extent of
   overlap among total number of DEGs (B), upregulated DEGs (C), and
   downregulated DEGs (D) at each time point.

   We further analyzed the DEGs of four time points in a Venn diagram
   ([79]Figure 1B) that allowed us to find and exclude genes counted as
   DEGs at different time points. We detected the greatest number of DEGs
   at the area of T3 and T6, with 1838 genes, followed by DEGs detected at
   T3 alone (1748). Of the 10,048 DEGs, 1035 were differentially expressed
   across all four time points compared to T0. We plotted the
   corresponding Venn diagrams for upregulated and downregulated DEGs
   across the four time points ([80]Figure 1C,D), which revealed a similar
   trend. There is also a case where it can be called ‘mixed pattern
   DEGs’. This indicates DEGs that were upregulated at one time point but
   downregulated at another time point. For example, if there is a gene
   that is upregulated at T1, but downregulated at T3, this gene is found
   in the T1-T3 co-zone (848, 8.4%) in [81]Figure 1B. This gene is also
   counted as one of the genes belonging to T1-only (527, 11.0%) and
   T3-only (1098, 9.7%) regions in [82]Figure 1C,D, respectively. There
   were 350 mixed-pattern DEGs discovered in this study, which appears to
   be because responses to stimulus occur through complex processes. The
   regulation of expression levels of these genes over a time period
   appears to be worthy of an independent research theme.

2.2. GO and KEGG Enrichment Analysis for Functional Classification of DEGs

   We performed a gene ontology (GO) term enrichment analysis on the DEGs
   to predict their functions following MeJA treatment and wounding. We
   identified 14 GO terms (level-2) containing more than 5% of DEGs in
   each of the three key categories (level-1, biological process; 6718,
   cellular component; 8033, and molecular function; 7689 genes). As shown
   in [83]Figure 2, the number of GO terms is proportional to the number
   of DEG. The number of genes associated with a significant GO term at
   the T3 time point was higher than at any other time point. Among the
   three key categories, we further analyzed the GO terms for the
   biological process category to predict gene functions in detail. To
   this end, we separately analyzed upregulated and downregulated DEGs
   (p-value cut off < 0.05), as shown in [84]Figure 3A and [85]Figure 3B,
   respectively. We selected the top 10 GO terms with the highest rich
   factor (%) for each time point.

Figure 2.

   [86]Figure 2
   [87]Open in a new tab

   Gene ontology term enrichment analysis of differentially expressed
   genes. The three major categories of biological process, cellular
   component, and molecular function are shown at level 2.

Figure 3.

   [88]Figure 3
   [89]Open in a new tab

   Enrichment analysis of differentially expressed genes in S. tora leaves
   after MeJA treatment and wounding. Gene ontology enrichment analysis of
   upregulated DEGs (A) and downregulated DEGs (B). Kyoto Encyclopedia of
   Genes and Genomes (KEGG) pathway enrichment analysis of upregulated
   DEGs (C) and downregulated DEGs (D). The results of enrichment analysis
   are shown for the top 10 gene rich factors at each time point. The rich
   factor is the ratio (%) of the number of DEGs assigned to each GO term
   (or KEGG pathway) relative to the number of genes assigned to that GO
   term (or KEGG pathway). The size of the dot indicates the gene rich
   factor. The color of the dot represents the false discovery rate (FDR)
   values. All p-values were below 0.05. The dot plots were generated
   using the ggplot function in the R package ggplot2.

   We determined that the upregulated DEGs at T1 ([90]Figure 3A) are
   related to defense responses, such as ‘systemic acquired resistance’
   (GO:0009627), ‘response to wounding’ (GO:0009611), ‘response to
   molecule of bacterial origin’ (GO:0002237), ‘response to chitin’
   (GO:0010200), ‘defense response to fungus, incompatible interaction’
   (GO:0009817), ‘cellular response to jasmonic acid stimulus’
   (GO:0071395), ‘cellular response to heat’ (GO:0034605), and
   ‘calcium-mediated signaling’ (GO:0019722). We speculate that this
   pattern reflects the rapid response of plants to externally applied
   MeJA or to wounding. At later time points (T3, T6, and T24), the GO
   term ‘jasmonic acid metabolic process’ (GO:0009694), which includes
   numerous JA biosynthetic genes, was highly enriched after external MeJA
   application. We also noticed the enrichment of DEGs at the T24 time
   points for the GO term ‘flavonoid biosynthetic process’ (GO:0009813). A
   similar GO analysis with downregulated DEGs ([91]Figure 3B) revealed
   that GO terms related to photosynthesis are highly represented at all
   time points: ‘response to red light’ (GO:0010114), ‘response to blue
   light’ (GO:0009637), ‘photosynthesis’ (GO:0015979),
   ‘protein-chromophore linkage’ (GO:0018298), ‘chlorophyll metabolic
   process’ (GO:0015994), ‘photomorphogenesis’ (GO:0009640), and ‘response
   to light intensity’ (GO:0009642). A decrease in photosynthetic activity
   following exposure to biotic and abiotic stresses was reported
   [[92]38,[93]39].

   We also performed a Kyoto Encyclopedia of Gene and Genomes (KEGG)
   pathway enrichment analysis to obtain additional information about the
   gene functions. This analysis indicated that upregulated DEGs are
   associated with the KEGG pathway ‘flavonoid biosynthesis’ (ko00941) at
   all time points ([94]Figure 3C). Moreover, downregulated DEGs were
   related to the KEGG pathway ‘photosynthesis’ (ko00195), consistent with
   the GO term analysis above ([95]Figure 3D). Thus, GO and KEGG
   enrichment analyses revealed that genes associated with flavonoid
   biosynthesis are significantly upregulated, while those related to
   photosynthesis are downregulated, in response to wounding and MeJA
   treatment.

2.3. Expression Analysis of the Flavonoid Biosynthetic Genes and
Quantification of the Flavonoids Kaempferol and Quercetin

   We investigated the transcript levels of genes encoding enzymes
   involved in flavonoid biosynthesis pathway (ko00941). We determined
   that 25 genes were the upregulated DEGs in the pathway. In the
   flavonoid biosynthesis pathway [[96]40,[97]41,[98]42], illustrated in
   [99]Figure 4A, 4-hydroxycinnamoyl CoA is generated through the
   phenylpropanoid biosynthetic pathway using phenylalanine as a
   substrate, and then three manlonyl CoAs are successively added to form
   a chalcone that forms the basic skeleton of flavonoids. This is the
   first committed step and is catalyzed by CHS. Among the 12 CHS genes of
   the S. tora genome [[100]37], 11 CHSs, except CHS10, were upregulated
   DEGs after MeJA treatment and wounding. In particular, CHS5, CHS6, and
   CHS7 showed significantly increased expression levels and reached
   maximum at T6. The flavonoid biosynthesis genes showed their maximum
   FPKM values mostly at T6 or T24 ([101]Figure 4B and [102]Supplementary
   Table S2). In addition, we quantified the levels of kaempferol and
   quercetin in the same samples; we observed a marked increase at T24
   ([103]Figure 4C). Together, these results suggest that kaempferol and
   quercetin are synthesized in response to the abundant expression of
   genes encoding flavonoid biosynthesis enzymes at T6 and/or T24,
   reaching a maximum accumulation at T24. Our results from this
   experiment are consistent with those of many previously published
   studies reporting that wounding [[104]43,[105]44] or MeJA treatment
   [[106]45,[107]46,[108]47,[109]48] increased the content of flavonoids.

Figure 4.

   [110]Figure 4
   [111]Open in a new tab

   Expression analysis of genes involved in flavonoid biosynthesis induced
   after MeJA treatment and wounding in S. tora leaves. (A) Flavonoid
   biosynthesis pathway in plants. PAL, phenylalanine ammonia-lyase; C4H,
   cinnamate−4−hydroxylase; 4CL, 4−coumarate−CoA ligase; CHS, chalcone
   synthase; CHI, chalcone isomerase; F3′H, flavonoid 3′-monooxygenase;
   F3H, flavanone 3−hydroxylase; and FLS, flavanol synthase. (B) Heatmap
   representation of the expression of genes involved in flavonoid
   biosynthesis. The FPKM (fragments per kilobase of transcript per
   million mapped reads) values for each gene were Z-score normalized. The
   heatmap was plotted using the R package pheatmap. (C) Kaempferol and
   quercetin contents in S. tora leaves at the indicated time points.
   Asterisk indicates significant differences from T0, with p-value < 0.05
   (Student’s t-test).

2.4. Expression Analysis of the Jasmonic Acid Biosynthetic Genes

   JA biosynthesis starts in the chloroplast and is completed in the
   peroxisome; it is then modified into MeJA or JA-Ile in the cytoplasm
   ([112]Figure 5A) [[113]8,[114]49]. Most JA biosynthetic genes were
   upregulated at T1 and some remained upregulated at T3 ([115]Figure 5B,C
   and [116]Supplementary Table S3). Application of exogenous MeJA and
   wounding dramatically increased the expression of two JMT genes at T1.
   This observation indicates that externally added MeJA and/or wounding
   influences the biosynthesis of endogenous MeJA. The JAR1 gene encoding
   the enzyme that conjugates isoleucine to JA [[117]50] was also strongly
   expressed at T1. This finding may reflect the need for JA-Ile to
   propagate signals for internal transmission [[118]51] as well as MeJA
   to be released outside the plant.

Figure 5.

   [119]Figure 5
   [120]Open in a new tab

   Expression analysis of genes involved in jasmonate biosynthesis
   following MeJA treatment and wounding in S. tora leaves. (A) Jasmonate
   biosynthesis pathway in plants. LOX, lipoxygenase; AOS, allene oxide
   synthase; AOC, allene oxide cyclase; OPR3, OPDA reductase, ACX,
   acyl−coenzyme A oxidase; AIM, ABNORMAL INFLORESCENCE MERISTEM (fatty
   acid methyltransferase); KAT, 3−ketoacyl−CoA thiolase; JMT, jasmonic
   acid carboxyl methyltransferase; JAR1, jasmonate resistant 1. The FPKM
   values for each gene were Z−score normalized. (B) Heatmap
   representation of the expression of genes involved in jasmonate
   biosynthesis. (C) Expression levels of JA biosynthesis genes.

2.5. Expression Analysis of the Ethylene, Salicylic Acid, and Abscisic Acid
Biosynthetic Genes

   Ethylene [[121]52,[122]53,[123]54], salicylic acid (SA) [[124]55], and
   abscisic acid (ABA) [[125]56,[126]57] are also important phytohormones
   that regulate growth, development, and response to biotic and abiotic
   stresses. Hence, we investigated the expression tendencies of the
   biosynthesis genes of these hormones over time periods. Unlike JA
   biosynthesis genes that showed immediate increases in expression level
   and maximum at T1, the genes of ethylene ([127]Figure 6A–C and
   [128]Supplementary Table S3) and salicylic acid ([129]Figure 6D–F and
   [130]Supplementary Table S3) biosynthesis pathways showed maximum
   expression levels at T6. Ethylene, SA, and jasmonates, including MeJA,
   form the basis of plant defense systems [[131]58]. It was reported that
   MeJA induces the ethylene biosynthesis organ specifically [[132]59],
   and that SA biosynthesis has an inverse correlation with MeJA
   [[133]60]. RNA-seq analysis in this study also suggests that the
   biosyntheses of ethylene and SA occur after endogenous MeJA
   biosynthesis induced by exogenous MeJA treatment and/or wounding.
   Although research on cross-talk between them is in progress, much
   effort is still needed to fully understand it.

Figure 6.

   [134]Figure 6
   [135]Open in a new tab

   Expression analysis of genes involved in ethylene, salicylic acid (SA),
   and abscisic acid (ABA) biosynthesis following MeJA application and
   wounding. (A) Ethylene biosynthesis pathway. (B) Heatmap representation
   of the expression of genes involved in ethylene biosynthesis. (C)
   Expression level in FPKM of ethylene biosynthesis genes. (D) Salicylic
   acid biosynthesis pathway. (E) Heatmap representation of the expression
   of genes involved in salicylic acid biosynthesis. (F) Expression level
   in FPKM of salicylic acid biosynthesis genes. (G) Abscisic acid
   biosynthesis pathway. (H) Heatmap representation of the expression of
   genes involved in abscisic acid biosynthesis. SAM, S−adenosylmethionine
   synthase; ACS, 1−aminocyclopropane−1−carboxylate synthase; ACO,
   1−aminocyclopropane−1−carboxylate oxidase; CM, chorismate mutase; PAT,
   prephenate aminotransferase; ADT, arogenate dehydratase; PAL,
   phenylalanine ammonia-lyase, BCH, beta-carotene hydroxylase; ZEP,
   zeaxanthin epoxidase; NSY, neoxanthin synthase; NCED,
   9−cis−epoxycarotenoid dioxygenase; ABA2, ABA DEFICIENT 2 (xanthoxin
   dehydrogenase); and AAO, abscisic aldehyde oxidase.

   Unlike ethylene and SA, the expression pattern for ABA biosynthetic
   genes was complicated and showed two opposite expression patterns after
   wounding and MeJA treatment. We identified 14 ABA biosynthesis DEGs of
   which 10 DEGs were downregulated at T1, and the remaining four were
   upregulated (maximum at T3), as shown in [136]Figure 6H and
   [137]Supplementary Table S3. It is clear that wounding or MeJA
   treatment affected the expression of these genes, as the expression
   levels at T1 or T3 were significantly increased or decreased compared
   to T0. However, we were unable to elucidate their correlation in this
   study. Previous references reported that ABA biosynthesis was not the