Full-length transcriptome profiling of Gentiana straminea Maxim. provides new insights into iridoid biosynthesis pathway

Preparation and collection of samples for transcriptome sequencing and qPCR analysis

Samples of G. straminea individuals were collected in August 2023 during the flowering stage from Yushu, Maqin County, Qinghai Province, China (N34°38′380, E100°23′546, altitude 4,200 m). Fresh tissues from five types—root, stem, leaves, flower, and ovary—were collected, washed with sterilized water, wrapped in foil, and then preserved in liquid nitrogen. Tissues of these five types were further utilized for library construction and SMRT sequencing. Additionally, for the quantification of gene expression via qPCR, two additional tissue types—non-embryonic callus (NEC) and embryonic callus (EC)—were included. The generation of EC and NEC tissues was carried out following the protocol outlined by He, Yang & Zhao (2011), utilizing leaves as explants. Three biological replicates were collected for all samples.

RNA extraction and SMRT sequencing

Samples of G. straminea were used for total RNA extraction on ice, according to the manufacturer’s protocol, with TRIzol reagent (Life Technologies, Karlsbad, CA, USA). An Agilent 2100 Bioanalyzer and agarose gel electrophoresis were used to determine the integrity of the total RNA. A Nanodrop microspectrophotometer (Thermo Fisher Scientific, Waltham, MA, USA) was used to check the purity and concentration of the RNA. The Clontech SMARTer PCR cDNA Synthesis Kit was used to reverse transcribe the oligo (dT) magnetic bead-enriched mRNA to cDNA. PCR cycle optimization was employed to identify the ideal number of amplification cycles for subsequent large-scale PCRs. Double-stranded cDNAs were generated with the optimized cycle number. Additionally, size selection at >5 kb and equal mixing without size selection of cDNA were performed with the BluePippin^TM Size Selection System. The next step in the construction of the SMRTbell library was carried out by large-scale PCR. The sequencing primer was matched with the SMRTbell template by annealing, and then linked to the polymerase. Sequencing was conducted using the PacBio Sequel II platform at Gene Denovo Biotechnology Co.

The raw sequencing reads from the cDNA library were classified via the Pacific Biosciences Iso-Seq pipeline, with high-quality circular consensus sequence (CCS) first extracted. Transcript integration was assessed according to whether the CCS reads contained all 5′primers, the 3′primer and the poly-A sequences. Full length sequences (FLs) were those that contained all three sequences. After the removal of primers, barcodes and poly A tails, full-length nonchimeric (FLNC) reads were obtained. Reads less than 50 bp in length were discarded. The entire isoform was generated by clustering the FLNC reads. Minimap2 was used for similar FLNC reads, which were then clustered hierarchically to obtain a consistency sequence (unpolished consensus isoforms). The consistency sequence was then further corrected via the Quiver algorithm. The high-quality isoforms (prediction accuracy ≥ 0.99) were used for subsequent analysis.

Library construction and Illumina sequencing

Total RNA was enriched by Oligo(dT) beads to form mRNA, then was fragmented into short fragments. With random primers, fragments were transcribed into cDNA, then syntheisized the second-strand cDNA with DNA polymerase I, Rnase H, dNTP and buffer. the obtained cDNA was purified with QiaQuick PCR extraction kit (Qiagen, Venlo, The Netherlands), end repaired, poly(A) added, and ligated to Illumina sequencing adapters. The ligated products were screened by agarose gel electrophoresis, amplified by PCR, and sequenced by Gente Denovo Biotechnology Co. (Guangzhou, China) using Illumina HiSeq^TM 4000. High quality clean reads were obtained by fastp (Version 0.18.0), with removing adapters, containing more than 10% of unknown nucleotides (N) and low-quality reads.

Isoform expression and differential expression analysis

Using the full-length transcriptome as the reference, the clean and high-quality reads were mapped using RSEM (version 1.2.8) to determine the isoform expression in five different tissues of G. straminea. The expression levels of isoforms from each sample were calculated and normalized to FPKM. Differentially expressed genes (DEGs) were identified using the DESeq2 software with —log2 (Fold Change)— ≥ 2 and a false discovery rate (FDR) below 0.05 as the screening criteria.

Functional annotation, structure analysis

The sequences of the isoforms were checked against the non-redundant protein (Nr) database of the NCBI (http://www.ncbi.nlm.nih.gov), the Clusters of Orthologous Genes/EukaryoticOrthologous Group (COG/KOG) database (http://www.ncbi.nlm.nih.gov/COG), the Kyoto Encyclopedia of Genes and Genomes (KEGG) database (http://www.genome.jp/kegg), and the Swiss-Prot protein database (http://www.expasy.ch/sprot) via the BLASTx program (http://www.ncbi.nlm.nih.gov/BLAST/), with an E value threshold of 1e⁻⁵, to assess the similarity of the sequences to those of genes from other species. Gene Ontology (GO) annotation was analyzed using isoforms from the Nr annotation results by Blast2GO software. The top 20 scoring isoforms and no fewer than 33 high-scoring segment pair hits (HSPs) were selected for the Blast2GO analysis. Isoforms were functionally classified using WEGO software. Transcription factors (TFs) were predicted via hmmscan by aligning the protein coding sequences to the Plant TFdb (https://planttfdb.gao-lab.org/).

Identification of the putative iriodoid biosynthesis pathway

Based on previous literature and known information from public databases, candidate isoforms for iridoid biosynthesis were screened from the sequence annotation files derived from third-generation sequencing. These candidate isoforms were subsequently annotated against the Nr, Swiss-Prot, KEGG, and KOG databases. Isoforms annotated with Enzyme Commission numbers (EC numbers) and keywords associated with terpenoid backbone biosynthesis, cytochrome P450, and iridoids were selected for further analysis. The screened isoforms were mapped to reference pathways, such as terpenoid backbone biosynthesis (map000900) and monoterpenoid biosynthesis (map00902), within the KEGG database. With reference to published core iridoid biosynthesis pathways (Oudin et al., 2007; Miettinen et al., 2014; Ni et al., 2019; Liu et al., 2017; Salim et al., 2014; Rai et al., 2013; Xu et al., 2022). These putative enzyme isoforms were assigned to corresponding biochemical reaction steps, thereby establishing a putative pathway for iridoid biosynthesis in G. straminea.

Protein-protein interaction network analysis

Functional enrichment analysis of protein-protein interaction (PPI) networks for genes involved in the iridoid biosynthesis pathway was conducted using the String database (https://cn.string-db.org/) (Szklarczyk et al., 2014). The protein sequences corresponding to key enzymatic components of the iridoid biosynthetic machinery in G. straminea were queried against the string database, with A. thaliana designated as the reference model organism for orthologous mapping. After sequence homology analysis, the highest-confidence protein isoforms (with open reading frames, ORFs) were selected based on maximal sequence identity from the mapping interface. Network construction parameters were optimized as follows: Minimum required interaction score: High confidence (score ≥ 0.700); max number of interactors to show: first shell limited to query proteins only; visualization output: dynamic Scalable Vector Graphics (SVG) network representation.

Identification and bioinformatic analysis of G(G)PPSs in G. straminea

For identification of GsG(G)PPSs, local BLAST search was performed using GGPPSs from Arabidopsis thaliana (Beck et al., 2013) or Chimonanthus praecox (Kamran et al., 2020) as queries. A threshold of e-value <10⁻¹⁰ was applied for preliminary screening. After searching, sequences of putative GsG(G)PPSs were further subjected to CDD (https://www.ncbi.nlm.nih.gov/cdd/) and InterPro (https://www.ebi.ac.uk/interpro/result/InterProScan/) for domain confirmation (Paysan-Lafosse et al., 2023). Prediction and analysis of the physicochemical properties of the GsG(G)PPS amino acid sequences were performed via ExPASy (https://web.expasy.org/protparam/) (Artimo et al., 2012). Sequences were submitted to SignalP4.1 server for prediction of the signal peptide (https://services.healthtech.dtu.dk/services/SignalP-4.1/) (Petersen et al., 2011). Subcellular localization was determined via WoLF PSORT (https://wolfpsort.hgc.jp/), transmembrane structure was predicted by HMHMM2.0 (https://services.healthtech.dtu.dk/services/TMHMM-2.0/). In addition, the annotated sequence information was submitted to the MEME website (https://meme-suite.org/meme/doc/meme.html), using 6–100 residues as the optimal motif size to search for 10 conserved motifs and predicted the conserved protein motifs in the sequence (Bailey et al., 2015). Similar to the GsGGPPS SSU, GsGGPPS, and GsGPPS amino acid sequences were downloaded from NCBI BLAST (https://blast.ncbi.nlm.nih.gov/Blast.cgi). Protein sequence alignment was performed via DNAMAN. Protein structure prediction was performed via SWISS-MODEL (https://swissmodel.expasy.org/).

Phylogentic analysis of G(G)PPS gene family

The GsG(G)PPSs obtained, and G(G)PPSs from other species were incorporated for phylogenetic analysis. G(G)PPSs in Arabidopsis thaliana and Nicotiana tabacum genomes were identified through BLAST by using the ensemble database (https://asia.ensembl.org/Multi/Tools/Blast). G(G)PPS homologues from other species available in the NCBI database were included for phylogenetic analyses of the G(G)PPS family. Details on all the G(G)PPSs used for phylogenetic analysis were listed in Table S1. These species protein sequences were performed multiple sequence alignment using the ClustalW program integrated in MEGA11.0. Phylogenetic inference of G(G)PPSs was conducted using the neighbor-joining method in MEGA 11.0 software, with a bootstrap test of 1,000 replicates (Tamura, Stecher & Kumar, 2021) and the best-fit substitution model JTT+G+I. The refinement of the evolutionary tree was completed using the online software Evoview (https://www.evolgenius.info/evolview/#/).

Expression analysis of key enzymes by real—time quantitative PCR

First-strand cDNA synthesis was performed using a cDNA reverse transcription kit (PrimeScriptTMII 1st Strand cDNA Synthesis Kit), following the protocol provided. Primers for reverse transcription-quantitative polymerase chain reaction (RT-qPCR) were designed using the OligoArchitect online sever and synthesized by Sangon Biotech (Shanghai) Co., Ltd. The primers sequence shown in Table S2, qPCR was performed using TB Green Premix Ex Taq with a 20 μL reaction system, which included 10 μL of TB Green Premix, 0.8 μL each of forward and reverse primers (10 μM), two μL of cDNA, 6.4 μL of ddH₂O. The reaction procedure consisted of the following steps: pre-denaturation at 94 °C for 5 min, denaturation at 94 °C for 30 s, annealing at 53 °C for 30 s, extension at 72 °C for 30 s, followed by 40 cycles. The GAPDH gene was utilized as the internal reference for relative expression analysis. The quantification of gene expressions was conducted using three biological replicates. Relative expression was calculated using the Ct (2^−ΔΔCt) method, following the approach described by Livak & Schmittgen (2001), with root expression serving as the control. The significance analysis of difference tissue was conducted with means of gene expression by the Duncan test at 5%.

Transcriptome sequencing of G. straminea

Both SMRT and Illumina sequencing were performed on the root, stem, leaf, flower, and ovary tissues of G. straminea. The average amount of raw data generated was 6.5 GB for second-generation sequencing per sample (Table 1). For PacBio sequencing, a total of 62.47 GB of raw data was obtained. The per-base coverage depth reached approximately 45× (circular consensus sequences corrected depth). A total of 23,318,162 subreads were generated from third-generation sequencing. After self-correction and merging, 499,496 circular consensus sequences (CCS) were generated, with an average CCS read length of 2,789 bp. The distributions for the number of passes and read lengths of these CCS reads are shown in Figs. S1A–S1B. The full-length nonchimeric sequences with high-precision CCS reads were identified, and similar FLNC reads were clustered hierarchically to obtain consensus sequences (Fig. S1C). A total of 41,785 high-quality isoforms (HQs) and 140 low-quality isoforms (LQs) were obtained after further correction. After removing redundant sequences, the total length of the isoforms was 84,861,577 bp. A total of 32,776 isoforms were obtained, and the lengths ranged from 165 to 10,169 bp, with an average length of 2,589.14 bp, an N50 of 2,767 bp, and a guanine-cytosine (GC) content of 41.43%. The number and length distribution of the isoforms are shown in Fig. S1D.

Functional annotation of the full-length transcriptome of G. straminea

The HQ unigenes were annotated against four functional annotation databases Nr, Swiss-Prot, KEGG, and KOG. A total of 31,434 (95.9%) unigenes were successfully annotated, while 1,342 were unannotated. The highest number of unigenes (31,235; 97.47%) was annotated to Nr database, followed by the KEGG database and the Swiss-Prot database, in which 30,990 (94.55%) and 27,622 (84.27%) unigenes, respectively, were annotated. The lowest number of unigenes was annotated in the KOG database (22,753; 69.42%). Summary, 21,742 common unigenes (66.34%) were annotated in all four databases (Fig. 1A). These findings were compared with those for 414 species annotated in the Nr database (top ten shown in Fig. 1B). The species with the most annotated sequence information was Coffea arabica, with 8,701 (27.86%) unigenes, followed by Coffea eugenioides, Coffea canephora, and Olea europaea, with 5,318 (17.03%), 3,512 (11.24%), and 1,051 (3.36%) unigenes, respectively.

The KOG analysis identified 22,753 unigenes, which could be classified into 25 categories (Fig. 2). The largest number of annotated genes was associated with general function prediction only 4,733 genes (20.80%), followed by 4,146 genes (18.22%) annotated to signal transduction mechanisms; 2,859 genes (12.57%) annotated to posttranslational modifications, protein turnover, and chaperones; 1,617 genes (7.11%) annotated to carbohydrate transport and metabolism; and 1,533 genes (6.74%) annotated to RNA processing and modification. The lowest number was observed for cell motility (37; 0.16%). In addition, 1,210 (5.30%) genes with unknown functions were identified.

The unigenes annotated by GO function analysis were associated with 51 GO terms, which were grouped into three categories: cellular component, molecular function, and biological process (Fig. S2). The top three GO enriched terms in the biological process category were cellular process, metabolic process, and response to stimulus, with 21,052, 18,417 and 7,559 genes, respectively. The top three enriched GO terms in the molecular function category were binding (18,864), catalytic activity (16,549), and transporter activity (3,175). In the cellular component category, cellular anatomical entity (16,552) and protein-containing complex (6,946) terms were highly enriched.

In the KEGG database, 30,990 unigenes of G. straminea were annotated and classified into five main categories and 19 subclasses (shown in Table S3). Of these annotated unigenes, 9,485 were mapped to specific KEGG pathways. The greatest number of genes (4,647; 48.99%) was annotated in metabolism pathways, followed by secondary metabolite biosynthesis (2,494; 26.29%), carbon metabolism (826; 8.71%), and biosynthesis of amino acids (635; 6.69%) (Table S4). Genes annotated to secondary metabolite biosynthesis pathways may be related to the synthesis of the medicinal components of G. straminea. In addition to carbon metabolism, the biosynthesis of amino acids and other metabolic pathways may be related to cellular osmotic regulation and the oxidative stress response. These annotated genes provide important sequence information for investigating the biosynthetic mechanism of the metabolites of G. straminea.

In this study, 708 annotated genes were found to participate in 20 standard KEGG secondary metabolism pathways in the transcriptome of G. straminea (Table S5); among these genes, 121 genes were annotated to the terpenoid backbone biosynthesis pathway and 67 genes were enriched in terpenoids (monoterpenoid, diterpenoid, sesquiterpenoid and triterpenoid biosynthesis) pathways. Furthermore, there were 84 genes involved in phenylpropanoid biosynthesis, while 41 genes were involved in flavonoids, isoflavonoid, flavone and flavonol biosynthesis. Additionally, 89 genes related to the synthesis of various alkaloids (indole, isoquinoline, tropane, piperidine, and pyridine alkaloid biosynthesis) were annotated, as shown in Table S5.

Predicting TFs

According to the assembly results, 1,151 genes were annotated as TFs, distributed in 51 TF families. The largest number of genes belonged to the GRAS family, with 128 genes (11.12%), followed by the ARF, C3H, bHLH, and WRKY families, with 92, 70, 66 and 66 genes, respectively. The least common families were the NF-YB (1), M-type (1), Whirly (1), AP2 (1), and YABBY (1) families. The ten TF families with the greatest number of genes in G. straminea are shown in Fig. S3.

DEGs analysis

A total of 32,470 genes were detected, and venn diagram analysis revealed that 31,330 genes were commonly found in the five tissues (Fig. 3B). In the comparisons of root and stem, root and leaf, root and flower, root and ovary, a total of 9,809, 10,503, 13,195, 9,699 DEGs were identified, respectively. Among these, 6,594, 5,762, 6,572, 5,727 DEGs were up-regulated and 3,260, 4,741, 6,623, 3,972 DEGs were down-regulated, respectively. In the contrast between leaf and stem, leaf and flower, leaf and ovary, a total of 6,980, 10,475, 10,006 DEGs were separately detected. Of these, 4,030, 4,707, 5,002 DEGs were up-reaulated and 2,950, 5,768, 5,004 DEGs were down-regulated. Furthermore, 8,855, 7,021 DEGs were identified in the stem vs flower, stem vs ovary comparison groups, respectively. Comparing with ovary, in the tissue of flower, 3,456 DEGs were up-regulated, and 2,218 DEGs were down-regulated, as shown in Fig. 3A.

In the comparisons of five different tissues, DEGs annotated were mainly enriched in the metabolic pathways, biosynthesis of secondary metabolites. Additionally, DEGs genes were enriched in the pentose and glucuronate interconversions in the root-vs-flower and root-vs-ovaries; carbon metabolism in leaves-vs-stem and flowers-vs-ovaries; amino sugar and nucleotide sugar metabolism in roots-vs-stem (Fig. S4).

Analysis of iridoid biosynthesis genes in G. straminea

Iridoid compound, which are common secondary metabolite components found in various medicinal plants, are the main components of G. straminea and have significant biological activity. By combining these results with previous research results (Ni et al., 2019; Liu et al., 2017), we identified a putative pathway for iridoid biosynthesis and the isoforms involved (Fig. 4). Our results revealed that 117 isoforms encoded 19 key enzymes (Table S6). The expression levels of these key enzyme isoforms in seven tissues are shown with a heatmap (Fig. 4), of which, GCPE, STR, ISPE, DXR, ISPH,7-DLS showed relatively high expression in leaves, other genes showed different expression patterns in different tissues. To further elucidate the interactions among these genes. The STRING database was employed to construct a protein-protein interaction (PPI) network for proteins annotated in the iridoid biosynthesis pathway of G. straminea. A. thaliana was used as the reference organism, the annotated protein sequences were uploaded to STRING database. Initial analysis revealed only 18 associated proteins, with subsequent network formation showing just 15 functionally connected proteins (Fig. 5A). The resulting PPI network comprised 15 nodes and 82 edges (Fig. 5A), exhibiting an average node degree of 10.9 and an average local clustering coefficient of 0.863. The PPI enrichment p-value was highly significant (<1.0e−16). STRING database annotations indicated experimentally validated interactions between: GGPPS SSU and GGPPS (Purple edges), AACT and HMGCS (light blue edges). The other interacting pairs were derived from curated databases and predicted interactions. The protein interaction network-associated genes displayed distinct tissue-specific expression profiles (Fig. 5B).

Bioinformatics analysis of G(G)PPS

Our results revealed that ten isoforms had GPPS/GGPPS annotations, three of which had open reading frames (ORFs). Two of these were annotated as GGPPS, while one isoform was annotated as GPPS, among the two GGPPSs, one was categorized as GGPPS small subunits (GGPPS SSU), and one was classified as a typical GGPPS. The prediction results revealed that the amino acid length of G(G)PPS (SSU) ranged from 342 to 424 aa, with corresponding molecular weights were 37.9 kDa and pI values were 5.81 (Table 1). Two of Gs GGPPS (SSU) possessed nagative GRAVY values, ranging from −0.050 to −0.187, indicating that these proteins have hydrophilicity. Gs GPPS had a positive GRAVY value (0.049), suggesting hydrophobicity of them. Three of Gs G(G)PPS (SSU) were identified no signal peptide. Two of Gs GGPPS (SSU) were localized in the chloroplast, Gs GPPS were predicted to be mitochondrion. Transmembrane structures were not predicted in any of the G(G)PPS proteins by TMHMM2.0 predictions (Table 2). Pfam protein structural domain prediction revealed a distinctive polyprenyl-synt domain shared by all the G(G)PPS proteins (Figs. S5A and S5B).

G(G)PPS usually contains two highly conserved aspartic acid-rich regions-FRAM and SARM with the sequences of DD(XX)₁₋₂D (D is aspartic acid, and X refers to any amino acid). The first conserved region FRAM (DDXXXXD) is consistent with the binding site of the substrate dimethylallyl diphosphate (DMAPP), and the second conserved region SARM (DDXXD) corresponds to the binding site of the substrate isopentenyl diphosphate (IPP). which affects the catalytic activity of G(G)PPS. Some G(G)GPPS proteins also have the characteristic sequence CXXXC (C is cysteine, and X refers to any hydrophobic amino acid) of their structural domain, which is essential for proteins-proteins interactions (Beck et al., 2013). Sequence alignment results revealed that the Gs GGPPS SSU sequences were similar to those of Pj GGPPS SSU, Ae GGPPS SSU, Ca GGPPS SSU, and Si GGPPS SSU2, with identity values of 81.74%, 81.55%, 81.49% and 81.61%, respectively, according to DNAMAN (Table S7). The identities of the Gs GGPPS sequences were similar to those of Cr GGPPS, Ca GPPS, Ce GGPPS and Gj GGPPS, with values of 74.06%, 72.29%, 72.04% and 71.28%, respectively (Table S7). The Gs GPPS sequences were similar to those of Ce SPPS, Ca SPPS, Cr GPPS1, Cr GPPS2, Gsy FPPS, Si SPPS and Na SPPS, with identities of 91.84%, 91.76%, 92%, 91.84%, 91.53%, 90.68%, and 90.75%, respectively (Table S7). Gs GGPPS SSU were enriched with one FARM (DD(XX)₂D) and two CXXXC regions. The Gs GGPPS subunit underwent a change in the second aspartic acid enrichment motif, from D to E, i.e, DDXXE (Fig. S5A). Gs GGPPS was enriched with one FARM region (DD(XX)₂D), one SARM region each (DDXXD), and one CXXXC region (Fig. S5A). Gs GPPS was enriched with two SARM regions (DDXXD) (Fig. S5B).

The analysis of the conserved motifs in Gs G(G)PPS revealed that both Gs GGPPS SSU and Gs GGPPS contained motif 1, 2, 4, 6, 9. Gs GGPPS SSU and Gs GPPS shared motif 3, 5, 7. Gs GGPPS and Gs GPPS shared common motif 7 and 8 (Fig. 6). The difference in motif composition may influence the function of Gs G(G)PPSs, leading to changes in catalytic activity, protein subcellular localization, and other aspects. Protein structure prediction revealed that Gs GGPPS SSU, Gs GGPPS, and Gs GPPS exhibited high structural similarity to their homologous counterparts from Mucuna pruriens (GMQE = 0.85), Handroanthus impetiginosus (GMQE = 0.81), and Catharanthus roseus (GMQE = 0.79), respectively. Gs G(G)PPS mainly contained α-helices and random coils in the tertiary structure (Fig. S6).

Phylogentic analysis of G(G)PPSs

Phylogenetic analysis revealed that the G(G)PPSs identified can be categorized into three distinct branches. Among them, Gs GGPPS, together with large subunits of GGPPS (GGPPS LSUs) and GGPPSs from other species, clustered into group 1. Gs GGPPS SSU along with the small subunits of GGPPS (GGPPS SSUs) from other species, was grouped into the second branch (group 2). Gs GPPSs formed the third branch, together with GPPS, SPPS, and FPPS from various other species (group 3) (Fig. 7). Gs G(G)PPSs were categorized into three distinct groups based on their sequence and functional divergence.

Expression analysis by real-time quantitative PCR

To further analyze the expression patterns of genes annotated to iridoid biosynthesis pathway in different tissues, RT-qPCR was performed on key genes of the iridoid biosynthesis pathway. We selected the upstream and downstream initiators of the MVA (AACT, MVD) and MEP (DXS, ISPH, GCPE) pathways, and these genes (IDI, GPPS) serving as key genes involved in intermediate production and skeletal formation steps. As shown in Fig. 8, the qPCR analysis revealed distinct tissue-specific expression patterns of these genes; AACT, IDI, ISPH, and GCPE showed their highest expression levels in leaves, DXS and GPPS were most abundantly expressed in stems, while MVD exhibited peak expression in non-embryogenic callus. Conversely, the lowest expression levels were observed for AACT, DXS, and ISPH in root tissues, and for IDI, MVD, GCPE, and GPPS in flowers. Significant differential expression across tissues was observed for all genes, with particularly pronounced variations for ISPH, GCPE, and GPPS. The RNA-seq results demonstrated some similar but also distinct patterns: AACT and GPPS showed highest expression in embryogenic callus, ISPH and GCPE in leaves, while MVD was uniquely most highly expressed in roots. For most genes (DXS, IDI, ISPH, GCPE, GPPS), the expression trends across different tissues were generally consistent between qPCR and RNA-seq data. However, discrepancies were noted for AACT and MVD in roots and leaves, as well as for DXS and MVD in the ovary.