Skip to main content

Transcriptomic data reveals the dynamics of terpenoids biosynthetic pathway of fenugreek

Abstract

Medicinal plants are rich sources for treating various diseases due their bioactive secondary metabolites. Fenugreek (Trigonella foenum-graecum) is one of the medicinal plants traditionally used in human nutrition and medicine which contains an active substance, called diosgenin, with anticancer properties. Biosynthesis of this important anticancer compound in fenugreek can be enhanced using eliciting agents which involves in manipulation of metabolite and biochemical pathways stimulating defense responses. Methyl jasmonate elicitor was used to increase diosgenin biosynthesis in fenugreek plants. However, the molecular mechanism and gene expression profiles underlying diosgening accumulation remain unexplored. In the current study we performed an extensive analysis of publicly available RNA-sequencing datasets to elucidate the biosynthesis and expression profile of fenugreek plants treated with methyl jasmonate. For this purpose, seven read datasets of methyl jasmonate treated plants were obtained that were covering several post-treatment time points (6–120 h). Transcriptomics analysis revealed upregulation of several key genes involved in diosgenein biosynthetic pathway including Squalene synthase (SQS) as the first committed step in diosgenin biosynthesis as well as Squalene Epoxidase (SEP) and Cycloartenol Synthase (CAS) upon methyl jasmonate application. Bioinformatics analysis, including gene ontology enrichment and pathway analysis, further supported the involvement of these genes in diosgenin biosynthesis. The bioinformatics analysis led to a comprehensive validation, with expression profiling across three different fenugreek populations treated with the same methyl jasmonate application. Initially, key genes like SQS, SEP, and CAS showed upregulation, followed by later upregulation of Δ24, suggesting dynamic pathway regulation. Real-time PCR confirmed consistent upregulation of SQS and SEP, peaking at 72 h. Additionally, candidate genes Δ24 and SMT1 highlighted roles in directing metabolic flux towards diosgenin biosynthesis. This integrated approach validates the bioinformatics findings and elucidates fenugreek’s molecular response to methyl jasmonate elicitation, offering insights for enhancing diosgenin yield. The assembled transcripts and gene expression profiles are deposited in the Zenodo open repository at https://doi.org/10.5281/zenodo.8155183.

Peer Review reports

Introduction

Fenugreek (Trigonella foenum-graecum) is one of India’s oldest medicinal plants from the Fabaceae family and is generally used in Ayurvedic medicine. This versatile plant serves both edible and non-edible purposes, encompassing applications in spices, fodder, dietary supplements, pharmaceuticals, and therapeutic practices [1].

Plants produce a large and diverse group of organic compounds called secondary metabolites, which can have medicinal properties. Thus, these metabolites have a very high economic value, although, they have a complex structure and their chemical synthesis is usually complex and expensive, hindering their development and investigation. Nonetheless, given their significant value, particularly in medicinal industries, understanding their chemistry, activity, and biosynthesis becomes imperative [2].

Fenugreek is a source of a biogenic steroid which is important in the pharmaceutical industry. It can be a suitable alternative to sweet potatoes, due to their high content of a steroid saponin called diosgenin [3]. Diosgenin has many applications in the pharmaceutical industry. For example, diosgenin is used for the production of steroid-type drugs. Despite their immense importance, products derived from secondary metabolites are naturally produced in limited quantities. For this reason, researchers seek methods to increase the production of these compounds in plants.

The production of secondary metabolites in plants can be stimulated through applying elicitors [4]. Elicitors are (a)biotic factors that activate the defense response of the plants against e.g., herbivores. One of the well-known elicitors is Methyl jasmonate, which is an important organic compound that plays a variety of regulatory roles in plant growth and development including axis elongation during embryogenesis, flower development, leaf senescence, root formation, and stomatal opening [5, 6]. It has been shown that treatment with 0.01% methyl jasmonate increased diosgenin levels by 10.5 times in fenugreek seedlings [7]. In addition, methyl jasmonate upregulated the expression level of two pivotal genes of the mevalonate pathway, the metabolic route leading to diosgenin: 3-hydroxy-3-methylglutaryl-CoA reductase and sterol-3-β-glucosyl transferase [78]. The application of exogenous jasmonic acid and methyl jasmonate is responsible for the induction of reactive oxygen species and subsequent defense mechanisms in cultured cells and organs [9]. It is also responsible for the induction of signal transduction, the expression of many defense genes followed by the accumulation of secondary metabolites. In other words, the accumulation of secondary metabolites in plant, cell, and organ cultures often occurs when cultures are subjected to varied kinds of stresses including elicitors or signal molecules.

There are many research studies focused on the effect of elicitors on secondary metabolites biosynthesis in medicinal plants, however, the underlying molecular mechanisms involved are not well studied. To understand the role of biosynthetic pathway genes and the effect of elicitors in increased biosynthesis and accumulation of diosgenin and other secondary metabolites, the expression patterns of genes involved in the biosynthesis pathway should be investigated.

To engineer the biosynthetic pathway, it is necessary to identify the biosynthetic pathway genes. Despite the fact the diosgenin biosynthetic pathway is partially identified (Fig. 1) [10], the data required for analyzing the molecular behavior of diosgenin biosynthetic pathway genes in fenugreek (e.g., RNA-sequencing or microarray analysis) is limited especially in response to elicitors [11]. It has been shown that meta-analysis can provide candidate genes for future research to improve our understanding of functions of secondary metabolites [12].

Fig. 1
figure 1

Biosynthetic pathway of diosgenin in fenugreek. Acetyl-coenzyme A (acetyl‐CoA) converts to cycloartenol via a series of multiple reactions. Cycloartenol is converted to either cycloartenol, leading to cholesterol formation (blue box), or 24 methylene cycloartenol, resulting in stigmasterol (green box). Cholesterol is then converted to 3β‐cholesta‐5‐en‐3,22‐diol, which ultimately results in the diosgenin biosynthesis (purple box). (Adapted from [10]). The heatmap shows the impact of treatment on expression changes of six genes on the biosynthetic pathway of diosgenin in fenugreek (see Result section for details)

The RNA-seq approach has been extensively used in the literature to characterize the expression profile of different samples and find differentially expressed genes (DEGs) [6, 13, 14]. To better understand the effect of methyl jasmonate on secondary metabolites of fenugreek and molecular mechanisms in secondary metabolite accumulation, in this study, we benefit from the raw RNA sequencing reads and performed de novo assembly to infer the transcriptome. DEGs during six different treatment stages of fenugreek compared to the control sample were found. Besides, an extensive Gene Ontology (GO) enrichment was conducted to identify important biological processes and molecular functions involved.

Results

De novo assembly of fenugreek transcriptome

The raw RNA-sequencing data were obtained from the National Center for Biotechnology Information (NCBI) database with accession number PRJNA508420, comprising fenugreek plant samples treated with methyl jasmonate harvested after 6, 12, 24, 48, 72, and 120 h of treatment together with the control treatment [13] (Supplementary Table 1). All the sequencing reads passed the quality control step (Supplementary Fig. 1). Since there was no reference transcript for fenugreek plants, we pursued the de novo assembly approach. We assembled each sample independently using the Trinity software [15], and deposited the transcripts in the public repository of Zenodo. For the rest of the analysis, we used the assembled transcripts based on the RNA-seq data of the control sample since it has the highest BUSCO completeness score (Supplementary Tables 1&2). Thus, we mapped the sequencing reads of other samples on this assembled transcript set (see Method section for details).

Differentially expressed genes of bioinformatics analysis show the significance of methyl jasmonate treatment in fenugreek

We analyze the expression profile of samples of six time points, namely, 6, 12, 24, 48, 72, and 120 h after methyl jasmonate treatment. To demonstrate the changes in the expression level of transcripts regulated by the methyl jasmonate application, we calculated the fold change (FC) in gene expression levels of different time points compared to the expression levels of the control sample (Fig. 2A and Supplementary Fig. 2). The volcano plot for the sample harvested after 120 h of treatment (Fig. 2A) shows the distribution of transcripts with increased or decreased expression levels. One can see that after 120 h the effect of methyl jasmonate on the plant is reduced resulting in a decrease in the fold change of expression (Supplementary Fig. 2). Besides, there are a few transcripts that experienced increased expression much more than the rest in the same timespan located on the top right side of the Volcano plot (Fig. 2A; Table 1). Of note, after 12 h of treatment, there are 2 transcripts with extreme increase in terms of expression level according to Supplementary Fig. 2 and Table 1. One of the transcripts is associated with the PAA1 protein belonging to the ATPase family, a group of enzymes that catalyze the breakdown of ATP to ADP or the reverse reaction [43].

Table 1 Transcripts with the highest increase in expression after the application of methyl jasmonate treatment associated to the point on the top right side of the Volcano diagram in Fig. 2A (FC: fold change)
Fig. 2
figure 2

Bioinformatics analysis of assembled transcripts. (A). Volcano plot for up- and down-regulated transcripts after 120 h treatment with methyl jasmonate. Data points in green (close to zero) represent unchanged transcription level of a specific transcript compared to the control treatment. The blue points and red points represent the up- and down-regulated transcripts within that specific time point, respectively. (B). The stacked bar chart represents the number of up-regulated (cyan) and down-regulated (red) transcripts expressed after 6, 12, 24, 48, 72, and 120 h of methyl jasmonate treatment. (C). Intersection of up-regulated transcripts. (D). Intersection of down-regulated transcripts. (E). Distribution of transcripts with increased expression 6 (UP-reg6h), 12 (UP-reg12h), 24 (UP-reg24h), 48 (UP-reg48h), 72 (UP-reg72h) and 120 h (UP-reg120h) after methyl jasmonate treatment illustrated with a Venn diagram. (F). Distribution of transcripts with decreased expression

To better understand the functional aspect of methyl jasmonate treatments, we analyzed GO terms based on three categories including molecular functions, biological processes, and cellular components [1617]. The number of transcripts that are involved in a specific function is shown in Table 2 (Full list in Supplementary Fig. 3). Of note, GO functions including metal ion binding, protein phosphorylation, proteolysis, and translation were found to be associated with the highest number of transcripts.

Table 2 The number of transcripts involved in each function of GO term in three categories in different time points

Several of up- or down-expressed transcripts are in common at different time points after treatment

We calculated the number of up- and down-regulated genes in each time point sample compared to the control (Fig. 2B). The highest number of transcripts with increased expression occurred for the sample of 120 h after treatment with methyl jasmonate (i.e., 5491 transcripts, which is 2.1 higher than the number of up-regulated transcripts at 6 h after treatment). The highest number of down-regulated transcripts corresponds to the sample at 72 h after treatment with 14,239 transcripts. Besides, in the 24 h methyl jasmonate treated plants, the lowest number of transcripts were upregulated (i.e., 1501 transcripts). Cumulatively, the highest number of transcripts with modified expression levels (regardless of being up- or down-regulated) was in 120 h (19,409 transcripts), while the lowest in 24 h (Fig. 2B). Besides, Fig. 2C(D) shows the extent of shared transcript among different time points that were up (down) regulated with 469 (6594) genes in common indicating a high variation in the number of transcripts. Venn diagrams in Fig. 2E-F show the numbers of transcripts in common between different time points with increased and decreased expression, respectively. We can see that 469 (6549) transcripts with increased (decreased) expression levels were in common at all time points (center of Fig. 2E&F), showing the diversity of changes in expression levels.

Transcription level of genes involved in the diosgenin biosynthesis is affected by methyl jasmonate treatment

Among all transcripts assembled, we identify candidate genes involved in the diosgenin biosynthesis [18] including SQS, SEP, CAS, SMT1, Δ24, and BGL (Fig. 1). Our results showed an increase of transcription level over time for most of the candidate genes. However, the only gene that showed a decreasing pattern over time was BGL (Fig. 3F). Besides, Fig. 3D shows an increase in the transcription level of SMT1 followed by a decreasing pattern after 6 h. The lower levels of BGL might be due to the fact the plant might prefer to accumulate the glycosylated form of diosgenin rather than its free form. It has been also previously shown that fenugreek plants are able to accumulate diosgenin in its glycosylated form; in the biosynthesis pathway of diosgenin, Glycoside saponins such as dioscin convert to diosgenin during the hydrolysis process [19].

Fig. 3
figure 3

Dynamics of expression for six genes after the treatment. (A) Squalene synthase (SQS), (B) Squalene Epoxidase (SEP), (C) Cycloartenol Synthase (CAS), (D) Sterol Methyltransferase 1 (SMT1), (E) Delta24-sterol reductase (Δ24), (F) Beta-glucosidase (BGL)

Squalene synthase (SQS) is present in the first step of the biosynthesis pathway of diosgenin (Fig. 1). After applying methyl jasmonate, an increase in its expression can be seen with some fluctuations (Fig. 3 and Supplementary Fig. 4). At 6 h, its expression is two times higher than that of the control sample. The final decrease might be due to the fact that SQS is placed upstream of the whole pathway and the effect of the applied elicitor gradually diminished. Next, Squalene Epoxidase (SEP) oxygenates the product of the previous step, i.e., squalene. The application of methyl jasmonate elicitor increases the intermediate compound 2,3oxidosqualene metabolite by affecting the squalene metabolite. SEP expression witnessed an increase first, followed by a slight decrease staying at a level higher than the control (Fig. 3).

Cycloartenol Synthase (CAS) converts 2,3oxidosqualene into cycloartenol. Its expression is first reduced and then increased at the later stages of the methyl jasmonate application (Fig. 3). It is implied that some intermediate pathway genes are affected later than the ones at the beginning of the pathway and hence are influenced later. Often the production of the previous compound in the pathway (in this case, 2,3oxidosqualene) might act as an enhancer of the CAS expression levels which needs time within such a multistep pathway. This gene is known as a rate-limiting enzyme, a key enzyme whose activity determines the overall rate of the metabolic pathway. A decrease in its expression probably causes methyl jasmonate to initially increase with having little effect on CAS, but then, methyl jasmonate can affect the CAS from upstream pathways justifying the later increase. Another important gene is sterol Methyltransferase 1 (SMT1) located in the pathway of diosgenin biosynthesis after CAS. This gene produces 24-methylene cycloartanol in a parallel pathway with Δ24 reductase. In 6 h after the treatment, there was an increase in SMT1 gene, followed by a decrease, reaching a lower level than the control.

After the effect of CAS on 2,3oxidosqualene and the production of cycloartenol metabolite, Δ24 facilitates cycloartanol formation (Fig. 1). RNA-seq data proved a continuous increase in Δ24 levels over time (Fig. 3). The increase in the expression levels of this gene reaches approximately 1.6 times the value at 72 h after the treatment. Next is Beta-glucosidase (BGL) which belongs to the glycosidase enzyme family found in many organisms. These are converted in various processes, including the metabolism of sugars, biomass, phytohormones activity, and cell wall ligninization with the involvement of plant chemical defense against pathogens [19]. Therefore, it is less likely to be present freely in the plant. Its expression is the lowest at 72 and 48 h after applying the treatment (Fig. 3).

Validation of bioinformatics analysis using real-time PCR

To validate the obtained results from the bioinformatics analysis we have grown three accessions of fenugreek plants in the greenhouse with the same conditions described by Zhou [7]. For the real-time PCR reaction, we used 4 genes including SQS, SEP, CAS, and SSR (Δ24). The results of the real-time PCR analyses confirmed that the overall trend of the expression of these genes is increasing. However, we are aware of differences in culture conditions (e.g., altitude, relative humidity, and genetic background of accessions) of the current study with the one where RNA-seq analysis was performed expecting minor effects on the results.

Relative gene expression of Squalene synthase (SQS) and fold changes compared to the control

SQS gene synthesizes one molecule of squalene by connecting two molecules of farnesyl pyrophosphate, which is the first committed step in the biosynthesis of sterols. We studied the relative expression of SQS in three fenugreek populations of Shiraz (Iran), Bushehr (Iran) and Egypt at different time points after treatment with methyl jasmonate (Table 3; Fig. 4A). The highest fold change corresponds to Bushehr’s sample at 72 h after treatment (FC = 155.03) and the lowest was related to Shiraz’s sample at 12 h after treatment (FC = 0.803). The general trend is an increase compared to the control, except at the time of 120 h after treatment. This result is consistent with the results of bioinformatics analysis, which shows an increase in gene expression in 72 h after treatment.

Table 3 Fold changes of the SQS, SEP, CAS and SSR genes at different time points after methyl jasmonate treatment
Fig. 4
figure 4

The expression profile of three populations after methyl jasmonate treatment (ns: not statistically significant, *: from 0.01 to 0.05, **: from 0.001 to 0.05, ***: from 0.0001 to 0.05, ****: less than 0.0001)

Relative gene expression of SEP, CAS and SSR

SEP gene performs the first oxygenation step in the biosynthesis of sterols by epoxidizing the squalene molecule. Because this gene is upstream of the diosgenin biosynthesis pathway, it can be one of the effective genes in increasing diosgenin metabolite. The results of the minimum significant difference analysis show a significant difference in the Bushehr mass at 72 h after applying the treatment compared to the control (at the significance level of 0.05). The highest multiplication changes occurred in the Bushehr mass at 72 h after treatment, 45 times higher than that of the control. The highest fold changes compared to the control of SEP is for the sample of Bushehr at 72 h after treatment (FC = 43.67) reported in Table 3; Fig. 4B. However, the lowest changes related to the sample Egypt at 12 h after the application of treatment (FC = 5.42).

CAS converts 2,3-oxido-squalane to cycloartenol, which is a necessary substrate to convert to phytosterols and other plant sterols. The mass of Bushehr and Egypt at 72 h after applying the treatment has increased about 45 times compared to the control. In Bushehr and Shiraz, all time intervals except 12 h, have a significant difference with the control at the significant level of 0.05 (Table 3; Fig. 4C). The highest (lowest) fold change was seen for the Egypt sample 72 (12) hours after treatment.

The Δ24-reductase gene (a.k.a., SSR) encodes the side branch of reductase by sterol enzyme and is responsible for converting cycloartenol to cycloartenol. It is considered the enzyme that determines the path of cholesterol synthesis [20]. Bushehr population has a significant difference with the control (at the significant level of 0.05) in all intervals except 12 h after treatment. The highest fold changes occurred in the Bushehr sample at 72 h after treatment (about 150 times) compared to the control. However, the lowest fold changes compared to the control was 3.47 for the sample of Egypt 12 h after treatment (Table 3; Fig. 4D). The general trend of fold changes for all three genes is an increase, except for 120 h after the treatment (compared to the control).

Discussion and conclusion

Fenugreek is known as a medicinal plant and a source of a biogenic steroid which is important in the pharmaceutical industry. The two main secondary metabolites in fenugreek are trigonelline and diosgenin. Diosgenin is a saponin found in fenugreek, in addition to its many therapeutic effects, known as the most active compound of fenugreek [21]. However, the amount of this metabolite in the fenugreek plant is low. Treatment with methyl jasmonate can result in increased expression of genes involved in the diosgenin biosynthesis which subsequently is reflected in an increase in the amount of this substance. To investigate the effect of methyl jasmonate on the increase or decrease of diosgenin, we investigated its effect on the genes involved in the diosgenin biosynthesis. For this purpose, the RNA-seq datasets deposited in the NCBI database were used with the Trinity software to assemble the transcriptome.

Differential gene expression analysis revelaed a decrease in the BGL gene expression pattern over time after methyl jasmonate application. also performed. This trend is likely due to plants tendency to preserve diosgenin in its pure form by preventing the formation of glucose molecules. For this purpose, the BGL enzyme should be less active resulting in its expression reduction. The increased count of up- and down-regulated transcripts at the 120-hour may result from the influence of methyl jasmonate on additional biosynthetic pathways and biological processes. These pathways may either exhibit lower sensitivity to the treatment or rely on intermediate metabolites generated throughout the process. This aligns with our anticipated effects of methyl jasmonate on gene expression. Upregulation of SQS and SEP in response to methyl jasmonate has been describedin other plants like Panax ginseng. Additionally, the application of methyl jasmonate is known to elevate diosgenin levels in T. foenum-graecum [22]. To validate the accuracy of the bioinformatics analysis, real-time PCR was conducted in the laboratory. As expected, an increase in SQS and SEP gene expression was observed. In the Bushehr population, both genes of SQS and SEP showed significant elevation 72 h after treatment ( expressed 155 and 43 times higher than the control sample) confirming the bioinformatic findings. The Egypt population showed the highest fold changes for the CAS gene expression at the same time point, consistent with previouslyy observed bioinformatic trends. Notably, the real-time PCR analysis showed a substantial150-fold increase in the SSR key gene, aligning with the bioinformatics results. Since this is a key gene in the diosgenin biosynthesis pathway, increasing its expression can have a direct effect in producing more diosgenin as the final product of the pathway.

Two pivotal genes, Δ24 and SMT1, which take part in the branch point, exhibited relatively contrasting transcriptional expression patterns. The gene expression result, alongside with diosgenin content, suggested that the Δ24 plays an essential role in directing the metabolic flux toward diosgenin production, further validating our findings. Conversely, our results clearly indicates the opposite impact of SMT1 in the diosgenin biosynthesis pathway. Previous studies on transgenic plants demonstrated that the over-expression of SMT1 results in decreased cholesterol production, while its down‐regulation causes higher biosynthesis of cholesterol.Real-time PCR results also revealed that a high expression of Δ24 resulted in a higher amount of diosgenin, whereas increased SMT1 expression is associated with reduced diosgenin production. This suggests that the branch point regulates diosgenin biosynthesis [1].. The obtained results show that SMT1 had a significant downward trend in expression after 12 h with minimum increases at120 hours, while Δ24 expression gene exhibited an increasing trend. Overall, our results indicate that the elicitor positively regulates the diosgenin biosynthetic pathway genes, particularly SMT1 and Δ24, potentially enhancing the final diosgenin yield in fenugreek plants.

Overall, the bioinformatic data were validated by real-time PCR experiments, indicating the precision of the RNA-sequencing data analysis. This valuable agreement indicates the accuracy of the RNA-seq analyses. It might not be always the case to have complete agreement between in silico analysis and wet lab results for all samples and all time points since the populations, cultivation, and greenhouse conditions and several other factors can affect the results. However, the general trend of gene expression in this study for the investigated genes was in agreement using both approaches.

Materials and methods

RNA sequencing data for fenugreek

In this study, RNA sequencing data for fenugreek were downloaded from the SRA of the NCBI database. In addition to the control sample of the fenugreek plant with accession code SRR8281656, there are six time points after treatment with methyl jasmonate including 6 h (SRR8281657), 12 h (SRR8281658), 24 h (SRR8281659), 48 h (SRR8281660), 72 h (SRR8281654) and 120 h (SRR8281655) which all come from high throughput sequencing machine, Illumina HiSeq 2000. All the reads are paired end with a length of 150 bp each. The SRA toolkit v3.0 was used to download the raw sequencing reads in SRA format with command line prefetch in the Linux operating system. This is followed by fastq-dump to convert them to FASTQ format. Then, FastQC v0.11.8 was used to check the quality of the downloaded sequencing reads.

The de novo transcriptome assembly

The Trinity software package v2.15 [15] was used to assemble the raw sequencing reads of fenugreek [13]. This took 73 CPU hours using 64GB of RAM. The TrinityStats.pl script from the Trinity package was used to report the assembly statistics including median contig length, GC percentage, contig N50, and total assembled bases. To evaluate the completeness of the assembled transcripts, BUSCO v4.12 (with argument -m transcriptome) was used with the eukaryota_odb10 dataset. We looked for the transcript level variation of the selected genes. For this purpose, we calculated FPKM (fragments per kilobase exon per million mapped fragments) at different time points after methyl jasmonate treatment.

BlAST analysis and GO enrichment

To compare assembled transcripts with proteins that are available in public databases, Basic Local Alignment Search Tool (BLAST) v2.14 was used against the universal protein resource knowledgebase (UniProtKB) including around 249 million sequences across the Tree of Life. In order to analyze the transcriptomes in terms of functional annotation, the Trinotate software v4.0 [23] was used, and Gene Ontology (GO) [1617] terms associated with the assembled transcripts were found.

Differential expression gene analysis

The raw sequencing reads of all samples were mapped to the assembled transcripts of the control samples using the bowtie2 software v1.3. Then, using the samtools package v.1.12, mapped sequences were extracted and compressed to produce the BAM file. Then, the Express package v1.5 was used to quantify the number of reads that are mapped to the transcripts resulting in the expression profile. Statistical analysis and visualization were done using R software v4.1 and Rstudio v1.4 based on in-house scripts. To estimate DEGs, the deseq2 package v1.38 in R was used, resulting in the fold changes and p-values [24]. Results were visualized using the ggplot2 package v3.4, ggVennDiagram v1.2, and UpsetR v1.4. In a Volcano diagram, the most upregulated (or downregulated) genes are toward the right (or left). The p-value threshold was 0.05 and the log fold change of > 1 was considered as the expression change threshold. Since the p-value has a negative transformation, there is an inverse relationship between the y-axis and the p-value. Besides, each point corresponds to each transcript and when a point is close to the center, it indicates a gene with an unchanged expression level. Genes with statistically significant increased expression are shown in blue. However, down-regulated transcripts are plotted in red.

We also benefited from the UpsetR plot to show the number of transcripts that were in common among all time points or were present in only a few time intervals [25, 26]. Each point in a row of the UpsetR plot shows the number of genes with increased expression level at one time point (or a few time points), but without expression changes in the rest of the time points.

Real-time PCR reaction

Real-time PCR reaction was performed in the laboratory to study the expression of genes including SQS, CAS, SEP, and SSR. Three plant populations including Shiraz (Iran), Bushehr (Iran), and Egypt were considered. Plants were treated with methyl jasmonate at a concentration of 0.015% and in time intervals of 12, 24, 48, 72, and 120 h after the treatment. They were removed and placed in a -80 °C freezer with liquid nitrogen. Then RNA was extracted from plant samples and the extracted RNA was converted into cDNA. Then, real-time PCR reaction was performed. First, it started with a 15-minute reaction at 95 °C followed by 35 cycles of 20 s at 95 °C and 35 cycles at 60 °C. Data were analyzed using the Livak analysis method (RQ2 = 2-ΔΔCT) [27]. The gene glyceraldehyde 3-phosphate dehydrogenase (GAP) was used as the reference gene [12, 14]. Because the number of samples was small, the Shapiro-Wilk normality test was used [8]. One-way ANOVA and non-parametric tests were used in the Prism software for variance analysis to compare the differences among groups of data. The results are reported as mean ± SE (standard error).

Data availability

The datasets analyzed during the current study are available in the Zenodo open repository at https://doi.org/10.5281/zenodo.8155183 and the NCBI database at https://www.ncbi.nlm.nih.gov/bioproject/PRJNA508420.

References

  1. Naika N MB, Sathyanarayanan, Sajeevan RS, Bhattacharyya T, Ghosh P, Iyer MS, Jarjapu M, Joshi AG, Harini K, Shafi KM, Kalmankar N. Exploring the medicinally important secondary metabolites landscape through the lens of transcriptome data in fenugreek (Trigonella Foenum Graecum L). Sci Rep. 2022;12(1):13534.

  2. Karuppusamy SA. A review on trends in production of secondary metabolites from higher plants by in vitro tissue, organ and cell cultures. J Med Plants Res 2009, 3(13), 1222–39.

  3. Vendl O, Wawrosch C, Noe C, Molina C, Kahl G, Kopp B. Diosgenin contents and DNA fingerprint screening of various yam (Dioscorea sp.) genotypes. Z für Naturforschung C. 2006;61(11–12):847–55.

    Article  CAS  Google Scholar 

  4. Klimek-Szczykutowicz M, Dziurka M, Blažević I, Đulović A, Apola A, Ekiert H, Szopa A, Livak KJ, Schmittgen TD.: Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) Method. Methods 2001, 25(4), 402–408.

  5. Cheong JJ, Do Choi Y. Methyl jasmonate as a vital substance in plants. TRENDS in Genetics 2003, 19(7), 409–13.

  6. Zhang Y, Kashkooli AB, Blom S, Zhao T, Bouwmeester HJ, Kappers IF. Capsicum Terpenoid Biosynthetic Module is Affected spider-mite Herbivory Plant Mol Biology. 2023;113(4):303–21.

    CAS  Google Scholar 

  7. Chaudhary S, Chikara SK, Sharma MC, Chaudhary A, Alam Syed B, Chaudhary PS, Mehta A, Patel M, Ghosh A, Iriti M. Elicitation Diosgenin Prod Trigonella foenum-graecum (Fenugreek) Seedlings Methyl Jasmonate Int J Mol Sci. 2015;16(12):29889–99.

    CAS  Google Scholar 

  8. Ciura J, Szeliga M, Grzesik M, Tyrka M. Changes Fenugreek Transcriptome Induc Methyl Jasmonate Steroid Precursors Revealed RNA-Seq Genomics. 2018;110(4):267–76.

    CAS  Google Scholar 

  9. Ho TT, Murthy HN, Park SY. Methyl Jasmonate Induc Oxidative Stress Accumulation Secondary Metabolites Plant cell Organ Cultures Int J Mol Sci. 2020;21(3):716.

    CAS  Google Scholar 

  10. Mohamadi M, Ebrahimi A, Amerian M. The expression enhancement of some genes involved in th diosgenin biosynthesis pathway in fenugreek treated with different levels of melatonin under salinity stress. Iran J Field Crop Sci. 2021;52(4):235–47.

    Google Scholar 

  11. Maroufi A, Lotfi M, Esmaeli A, Dastan D. Relative expression of the key genes of diosegnin biosynthesis in fenugreek (Trigonella foenum-graesum) in response to salicylic acid and methyl jasmonat. Cell Mol Res (Iranian J Biology). 2021;34(3):440–54.

    Google Scholar 

  12. Tahmasebi A, Ebrahimie E, Pakniyat H, Ebrahimi M, Mohammadi-Dehcheshmeh M. Tissue-specific transcriptional biomarkers in medicinal plants: Application of large-scale meta-analysis and computational systems biology. Gene 2019, 691, 114–24.

  13. Zhou C, Li X, Zhou Z, Li C, Zhang Y. Comparative transcriptome analysis identifies genes involved in diosgenin biosynthesis in Trigonella foenum-graecum L. Molecules 2019, 24(1), 140.

  14. Valipour E, Nooshabadi VT, Mahdipour S, Shabani S, Farhady-Tooli L, Majidian S, Noroozi Z, Mansouri K, Motevaseli E, Modarressi MH. Anti-angiogenic Eff testis-specific gene Antigen 10 Prim Endothelial Cells Gene. 2020;754:144856.

    CAS  Google Scholar 

  15. Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z. Trinity: Reconstructing full-length Transcriptome without Genome RNA-Seq data Nat Biotechnol. 2011;29(7):644.

    CAS  PubMed  Google Scholar 

  16. Central GO, Aleksander SA, Balhoff J, Carbon S, Cherry JM, Drabkin HJ, Ebert D, Feuermann M, Gaudet P, Harris NL, Hill DP. The Gene Ontology knowledgebase in 2023. Genetics, 2023 224(1).

  17. Nicheperovich A, Altenhoff AM, Dessimoz C, Majidian S. OMAMO: orthology-based Altern Model Organism Selection Bioinf. 2022;38(10):2965–6.

    CAS  Google Scholar 

  18. CiuraJ,SzeligaM,GrzesikM,TyrkaM:Next-generation sequencing of representational difference analysis products for identification of genes involved in diosgenin biosynthesis in fenugreek (Trigonella foenum-graecum).Planta 2017,245,977–991.

  19. Aminkara S, Shojaeiyan A, Monfaredb S, Ayyari M. Special effect of ionic liquids on extraction of diosgenin from fenugreek (Trigonella foenum-graecum l.) by ultrasonic assistance. Nat Prod Chem Res. 2018;6(322):2.

    Google Scholar 

  20. Sonawane PD, Pollier J, Panda S, Szymanski J, Massalha H, Yona M, Unger T, Malitsky S, Arendt P, Pauwels L, Almekias-Siegl E.: Plant cholesterol biosynthetic pathway overlaps with phytosterol metabolism. Nature plants 2017, 3(1): 16205.

  21. Bitarafan Z, Asghari HR, Hasanloo T, Gholami A, Moradi F, Khakimov B, Liu F, Andreasen C.:The effect of charcoal on medicinal compounds of seeds of fenugreek (Trigonella foenum-graecum L.) exposed to drought stress.Ind Crops Prod 2019,131,323–329.

  22. Cao H, Nuruzzaman M, Xiu H, Huang J, Wu K, Chen X, Li J, Wang L, Jeong JH, Park SJ, Yang F. Transcriptome analysis of methyl jasmonate-elicited Panax ginseng adventitious roots to discover putative ginsenoside biosynthesis and transport genes. International journal of molecular sciences 2015, 16(2), 3035–57.

  23. Griffith M, Walker JR, Spies NC, Ainscough BJ, Griffith OL. Inf RNA Sequencing: Web Resource Anal Cloud PLoS Comput Biology. 2015.;11(8):e1004393.

  24. McDermaid A, Monier B, Zhao J, Liu B, Ma Q. Interpretation of differential gene expression results of RNA-seq data: review and integration, Briefings in Bioinformatics 2019, Volume 20, Issue 6, Pages 2044–54.

  25. Conway JR, Lex A, Gehlenborg N. UpSetR: an R package for the visualization of intersecting sets and their properties. Bioinformatics. 2017;33(18):2938–40.

    Article  CAS  PubMed  PubMed Central  Google Scholar 

  26. Ezazi R, Ahmadzadeh M, Majidian S, Stefani E, Pindo M, Donati C. Responses of cucumber (Cucumis sativus L.) rhizosphere microbial community to some agronomic management practices. FEMS Microbiol Ecol. 2021;97(8):fiab107.

    Article  CAS  PubMed  Google Scholar 

  27. Livak KJ, Schmittgen TD. Analysis of relative gene expression data using real-time quantitative PCR and the 2(-Delta Delta C(T)) method. Methods. 2001;25(4):402–8.

Download references

Acknowledgements

Authors would like to thank Dr. Amin Ebrahimi for his support during q-PCR analysis experiments and Dr. Irene Consuelo Julca Chavez for helpful discussions.

Funding

The project was partly funded by Tarbiat Modares University.

Author information

Authors and Affiliations

Authors

Contributions

SL, ABK, AS and SM hypothesized the research. SL and SM performed the bioinformatics analysis. SL performed the molecular experiments. ABK, SL and SM drafted the manuscript. SL and SM prepared the graphs. All authors reviewed the final version of the manuscript.

Corresponding author

Correspondence to Arman Beyraghdar Kashkooli.

Ethics declarations

Ethics approval and consent to participate

Not applicable.

Consent for publication

Not applicable.

Competing interests

The authors declare no competing interests.

Additional information

Publisher’s Note

Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations.

Electronic supplementary material

Below is the link to the electronic supplementary material.

Supplementary Material 1

Rights and permissions

Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. The Creative Commons Public Domain Dedication waiver (http://creativecommons.org/publicdomain/zero/1.0/) applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Javan, S.L., Kashkooli, A.B., Shojaeiyan, A. et al. Transcriptomic data reveals the dynamics of terpenoids biosynthetic pathway of fenugreek. BMC Genomics 25, 390 (2024). https://doi.org/10.1186/s12864-024-10253-x

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: https://doi.org/10.1186/s12864-024-10253-x

Keywords