De novo leaf and root transcriptome analysis to explore biosynthetic pathway of Celangulin V in Celastrus angulatus maxim

Celastrus angulatus Maxim is a kind of crucial and traditional insecticidal plant widely distributed in the mountains of southwest China. Celangulin V is the efficient insecticidal sesquiterpenoid of C. angulatus and widely used in pest control in China, but the low yield and discontinuous supply impeded its further popularization and application. Fortunately, the development of synthetic biology provided an opportunity for sustainable supply of Celangulin V, for which understanding its biosynthetic pathway is indispensable. In this study, six cDNA libraries were prepared from leaf and root of C. angulatus before global transcriptome analyses using the BGISEQ-500 platform. A total of 104,950 unigenes were finally obtained with an average length of 1200 bp in six transcriptome databases of C. angulatus, in which 51,817 unigenes classified into 25 KOG classifications, 39,866 unigenes categorized into 55 GO functional groups, and 48,810 unigenes assigned to 135 KEGG pathways, 145 of which were putative biosynthetic genes of sesquiterpenoid and triterpenoid. 16 unigenes were speculated to be related to Celangulin V biosynthesis. De novo assembled sequences were verified by Quantitative Real-Time PCR (qRT-PCR) analysis. This study is the first report on transcriptome analysis of C. angulatus, and 16 unigenes probably involved in the biosynthesis of Celangulin V were finally collected. The transcriptome data will make great contributions to research for this specific insecticidal plant and the further gene mining for biosynthesis of Celangulin V and other sesquiterpene polyol esters.


Background
With the increase of global population, natural insecticides become a trend of pesticide development for their environmental friendliness. Celastrus angulatus Maxim, extensively distributed in the mountains of southwest China, have been exploited as natural insecticide resource and folk medicines attributed to their active ingredients, including sesquiterpenes (β-agarofurans), alkaloids, and flavonoids [1].
Celangulin V, the major insecticidal ingredient of C. angulatus,, exhibits excellent activities against several insect species like Pieris rapae (Lepidoptera) and Locusta migratoria manilensis (Orthoptera) [1,2] targeting on Na + /K + -ATPase (sodium/potassium pump) of these insects. As precisely identified, Celangulin V is a kind of sesquiterpene polyol esters with β-dihydroagarofuran skeleton, holding seven hydroxyl groups [1]. In the past 20 years, a series of Celangulin V derivates have been isolated by using bioassay-directed fractionation. A total of 44 bioactive compounds were firstly reported, displaying antifeedant, narcotic or toxic effects against some important agricultural pests [3,4].
Celangulin V can be extracted from roots and leaves of C. angulatus, and the content in roots is higher than that in leaves [5][6][7]. The current supply is dependent on filed-grown plant with unstable yield. To address such problems, attempts have been made to produce Celangulin V in more eco-friendly means of biosynthesis. In general, the biosynthetic pathway of sesquiterpenes could be divided into three parts. Part I is isopentenyl diphosphate (IPP) and dimethylallyl pyrophosphate (DMAPP) biosynthesis [8,9]. Part II is carbo-cyclization of sesquiterpene biosynthesis. Part III is the modification pathway of complex functional groups in the sesquiterpene skeleton [10][11][12][13]. The whole reaction is catalyzed by a large variety of enzymes with substrate specificities (Fig. 1).
So far, the researches of C. angulatus mainly focused on its phytochemical properties, including the extraction of active ingredients, the identification of active substances and the evaluation of insecticidal activity [14]. As the most important bioactive component of C. angulatus, Celangulin V was also well studied in structures and physicochemical properties [15,16]. However, its biosynthetic and metabolic pathway remained poorly investigated, mainly owing to the absence of genomic or transcriptomic resources for this non-model species. Although no gene sequence information is available in public databases, we completed the transcriptome sequencing and analysis referring to the genome information of plants closely related to C. angulatus. Thus, comprehensive genomic information can be further analyzed for gene expression, molecular mechanisms and biological pathways [17,18]. This strategy has been efficiently used in traditional Chinese medicinal species such as Gardenia jasminoides [19], Pseudostellaria heterophylla [20], Panax notoginseng [21] and Cephalotaxus hainanensis [22]. In this study, functional genes involved in sesquiterpene biosynthesis in C. angulatus were screened using de novo transcriptome sequencing. Four unigenes encoding sesquiterpene synthases, eight unigenes encoding cytochrome P450s (CYP450s) and four unigenes encoding acyltransferases were identified from Leaf and Root transcription databases of C. angulatus through BGISEQ-500. The candidate transcript sequences were listed in the Additional file 1. Differentially expressed genes (DEGs) were also analyzed. Such databases could be used as an important resource to investigate the biosynthetic pathway of Celangulin V in C. angulatus. Furthermore, this database will supply important clues to explore biological characteristics genetically in other plants with close relationship to C. angulatus.

Plant materials
Leaf and root tissues of C. angulatus were randomly picked from the exposition park of Northwest A&F University, followed by transcriptome analysis on May 11, 2017. Leaf and root samples were harvested from three plants for RNA extraction and three biological replicates (Leaf_1 and Root_1, Leaf_2 and Root_2, Leaf_3 and Root_3) were performed. Tissues were rinsed in water, cut into small pieces, frozen in liquid nitrogen immediately, and stored at − 80°C for further analyses. The C. angulatus was authenticated by Prof. Wenjun Wu and Prof. Zhaonong Hu.
cDNA library preparation and BGISEQ-500 sequencing for transcriptome analysis Total RNA was extracted using an RNA plant Plus Reagent (Tiangen, Beijing, China) according to the manufacturer's protocol. The extracted RNA was checked using a NanoDrop 2000 (Thermo, CA, USA), and the RNA concentration and integrity were assessed using the RNA Nano 6000 Assay Kit of the Agilent Bioanalyzer 2100 system (Agilent, CA, USA) to ensure that the RNA Integrity Number (RIN) values were above 7.0.
Oligo (dT) beads were used to isolate poly(A) + mRNA, which was fragmented to 250 bp. Fragmentation of the RNA and reverse transcription of double-strand cDNA (ds cDNA) by N6 random primer. The synthesized cDNA was subjected to end-repair and then was 3′ adenylated. Adaptors were ligated to the ends of these 3′ adenylated cDNA fragments; The ligation products were purified and many rounds of PCR amplification were performed to enrich the purified cDNA template using PCR primer; Denature the PCR product by heat and the single strand DNA is cyclized by splint oligo and DNA ligase; Each cDNA library was sequenced in a single lane of the BGISEQ-500 system with paired-end sequencing length of 100 bp according to the manufacturer's instructions at the Beijing Genomics Institute (BGI-Shenzhen, China). The amounts of reads generated per sample was 10-11 Gb to obtain deep coverage of transcripts for de novo assembly.
De novo assembly and functional annotation analysis of BGISEQ-500 sequencing To obtain high-quality clean read data for de novo assembly, the raw reads from BGISEQ-500 were filtered by discarding the reads with adaptor sequences, reads with ambiguous "N" bases larger than 5% and removing the low-quality reads in which more than 20% bases had a Q-value < 15. All the downstream analyses were based on the resulting clean reads. We used the software Trinity (v2.0.6) [23], which was efficient to form contigs in reconstructing full-length transcripts across a broad range of expression levels and sequencing depths, with default parameters and a minimum contig length of 150 bp for assembly generation. These contigs were then further processed with sequence clustering software TGICL (v2.0.6) [24] to remove the redundant Trinity generated contigs.
To determine the functional annotation of the unigenes, a BLASTx search was performed with an E-value of 10 − 5 against protein databases, including Nr (non-redundant) protein database, SwissPort, KOG (euKaryotic Orthologous Group database), and KEGG (Kyoto Encyclopedia of Genes and Genomes protein database). Besides, a BLASTn search was also performed against the Nt database. With Nr annotation, the Blast2GO [25] and InterProScan5 program was used to obtain the GO (Gene ontology) and InterPro annotation of unigenes, respectively. GO classification was then performed using WEGO software [26] to illustrate the distribution of gene functions including Biological Process, Cellular Component and Molecular Function.

Differentially expressed unigene analysis
After assembly, clean reads were mapped to unigenes using Bowtie2 (v2.2.5) [27], and then gene expression level was calculated with RSEM (v1.1.12) [28]. To compare the difference of gene expression among different samples, the FPKM (Fragments per kilobase per transcript per million mapped reads) method was used [29]. The DEseq2 (Fold Change > = 2.00 and Adjusted Pvalue <= 0.05) and PossionDis (Fold Change > = 2.00 and FDR < = 0.001) were proposed to identify DEGs, and the P-value and FDR (false discovery rate) for each gene were calculated. DEGs were required to have thresholds of "log2 ratio ≥ 1" and "FDR < 0.001" [30]. Next, GO and KEGG analysis were again performed on the DEGs.

Phylogenetic analysis
The phylogeny was inferred using the Neighbor-Joining method [31]. The bootstrap consensus tree inferred from 500 replicates is taken to represent the evolutionary relationship of the taxa analyzed [32]. Branches corresponding to partitions reproduced in less than 50% bootstrap replicates are collapsed. The percentage of replicate trees in which the associated taxa clustered together in the bootstrap test (500 replicates) are shown next to the branches. The analysis involved 38 terpene synthase and 52 CYP450 nucleotide sequences, respectively. All positions containing gaps and missing data were eliminated. There was a total of 479 and 1094 positions in the final dataset, respectively. Evolutionary analyses were conducted in MEGA7 [33].

qRT-PCR analysis
Total RNA was extracted as indicated above. Each RNA sample was treated with RNase-free DNase (Tiangen, Beijing, China) following the manufacturer's protocol in an effort to remove any residual genomic DNA (gDNA). DNase-treated RNA (2 mg) was subjected to reverse transcriptase reactions using M-MLV reverse transcriptase (Tiangen, Beijing, China) according to the manufacturer's instructions. The sequences of the specific primer sets are listed in Additional file 2. The Actin (CL5382.Contig2_All) gene was used as an endogenous control, qRT-PCR was performed in 96-well plates in a Bio-Rad CFX96 real-time PCR system (Bio-Rad, CA, USA) with a SYBR Green-based PCR assay. The final volume for each reaction was 20 μL with the following components: 2 μL diluted cDNA template (1 mg/mL), 10 μL SYBR Green Mix (Bio-Rad, CA, USA), 0.4 μL forward primer (10 μM), 0.4 μL reverse primer (10 μM) and 7.2 μL ddH 2 O. The reaction was conducted under the following conditions: 95°C for 3 min, followed by 40 cycles of denaturation at 95°C for 10 s and 60°C for 30 s. The melting curve was obtained by heating the amplicon from 65°C to 95°C at increments of 0.5°C per 5 s. Each qRT-PCR analysis was performed with three biological replicates. The relative quantification of gene expression was computed using the 2 -ΔΔCt method.

Active ingredient content in different tissues of C. angulatus
Prof. Wenjun Wu have detected the active ingredient content in different tissues of C. angulatus. He extracted the adhesive paste of the root bark, bark, xylem, leaves and fruits of C. angulatus and determined the bioactivities of the extracts by feed poisoning method. As a result, the active ingredients of C. angulatus were mainly distributed in the root bark, secondly in the leaves. However, bark and xylem contained little active ingredients, and the extracts from the fruit stimulated the growth of larvae of Mythimna separate at low concentrations (Table 1) [34] .

C. angulatus transcriptome sequencing and unigene assembly
In order to explore gene expression profile of C. angulatus, six cDNA libraries prepared from three leaf tissues (Leaf_1, Leaf_2 and Leaf_3) and three root tissues (Root_1, Root_2 and Root_3) were sequenced using the BGISEQ-500 platform. A total of 110. 16 (Table 3).
After de novo assembly of two C. angulatus tissues, 104,950 unigenes were finally obtained with an average length of 1200 bp (Table 3). Among these, 48,715 unigenes have a length longer than 1 kb (> 1000 bp) and 16,934 unigenes have a length range between 1000 bp and 1500 bp. A detailed summary of the sequencing and assembly results is shown in Table 3 and the length distribution of all unigenes is shown in Fig. 2.

Functional annotation of C. angulatus unigenes
The contigs in six transcriptome sequencing databases were integrated and assembled into a total of 104,950 unigenes. After the assembly, annotations for the assembled unigenes were carried out via BLAST in seven public databases (NR, NT, Swissprot, KEGG, KOG, InterPro and GO), a total of 66,193,53,315,44,141,48,810,51,817,56,377 and 39,866 were aligned, respectively. There were 71,479 unigenes (68.11% of all unigenes) annotated within at least one functional database. (Additional file 3).
Functional classification of C. angulatus unigenes by KOG, GO and KEGG KOG classification was used to further evaluate the completeness and effectiveness of the C. angulatus annotation. 51,817 unigenes were classified into 25 functional classifications (Fig. 3). The most dominant term was "General function prediction only" and 12,607 unigenes (24.3%) matched it, followed by "Signal transduction mechanisms" (6864, 13.3%), "Posttranslational modification, protein turnover, chaperones" and so on. Notably, 1768 unigenes in the "secondary metabolites biosynthesis transport and catabolism" category may play important roles in the biosynthesis of Celangulin V and other secondary metablites. GO analysis is an international standard system of gene function classification, and the mainly terms including "biological process", "cellular component" and "molecular function". Among the 104,950 unigenes, 39,866 unigenes were assigned to 55 GO terms. 24 groups were involved in biological process, in which "metabolic process" and "cellular process" occupied the mainly categories. 17 groups were mainly existed in "cellular component" and the high percentage of unigenes was associated with "cell", "cell part" and "membrane". 14 groups existing in "molecular function", "catalytic activity" and "binding functions" play the most assignments in this category (Fig. 4).
The sesquiterpene polyol ester Celangulin V is the major insecticidal active component in C. angulatus. 145 unigenes involved in sesquiterpenoid and triterpenoid biosynthesis were found in the KEGG pathways, and the detailed metabolic pathway is shown in Additional file 4. The pathway will be useful for further studies on the biosynthesis of Celangulin V.

DEGs in the leaf vs root in C. angulatus
DEGs of the six transcriptome libraries were used to discover the unigenes with significant differences in expression. The expression of unigenes was calculated by FPKM. The different analysis methods were as follows: Leaf_1 Vs Root_1, Leaf_2 Vs Root_2, and Leaf_3 Vs Root_3 libraries were analyzed, and the DEGs were used to comment on all three replicates for GO classification and KEGG pathway analysis. A total of 19,440 DEGs were obtained, including 9024 up-regulated and 10,416 down-regulated unigenes in Leaf Vs Root. Furthermore, 5667, 1994 and 2120 unigenes expressed uniquely in Leaf_1, Leaf_2, and Leaf_3, respectively, and 75,738 unigenes were expressed in all three libraries, but at different levels; 3612, 2621 and 4541 unigenes expressed uniquely in Root_1, Root_2 and Root_3, respectively, and 70,006 unigenes were expressed in all three libraries, but at different levels (Fig. 5).
A total of 11,948 DEGs could be annotated in KEGG based on sequence homologies, which were annotated into 135 pathways and most genes were related to metabolic pathways (2854 DEGs) and biosynthesis of secondary metabolites (1764 DEGs). For the GO classification analysis, about 6366 7082 and 7890 unigenes could be annotated in "biological process", "cellular component", "molecular function" of GO terms respectively. Furthermore, the number of up-regulated genes was more than the number of down-regulated genes in the three pathways.

Biosynthetic genes of terpenoid backbone in C. angulatus
Based on the KEGG annotation, a total of 310 contigs/unigenes are correlated with the terpene   Table 5).

Validation and expression analysis of key genes
To confirm the accuracy of the BGISEQ-500 sequencing and FPKM calculated results, we selected 15 unigenes and used qRT-PCR to determine their relative expression level in the leaf and root tissues of C. angulatus. All 15 unigenes were putative sesquiterpene biosynthetic genes containing five higher-expressed unigenes (CL1103.Contig4_All, CL5925.Contig1_All, Unigene27813_All, CL1173.Conti-g1_All, and Unigene225446_All), five lower-expressed

Discussion
BGISEQ-500 sequencing and sequence annotation C. angulatus is a kind of excellent traditional insecticidal plant due to their ample active ingredients of sesquiterpenes (β-agarofurans), alkaloids and flavonoids. Even though Celangulin V is considered as the most vital active constituents of C. angulatus, little is known about the biosynthetic and metabolic mechanisms of this sesquiterpene polyol ester. The aims of this study were to generate a large amount of cDNA sequence data that would facilitate more detailed studies on C. angulatus, and to identify the genes related to sesquiterpene polyol esters biosynthesis and metabolism. The availability of transcriptome data for C. angulatus will meet the initial information needs for functional studies of this species and its relatives. Six RNA-seqs were performed using BGISEQ-500 sequencing, which generated a total of 104,950 unigenes. 71,479 (68.11%) unigenes provided a significant BLAST result. This is the first report about transcriptome study of C. angulatus, providing adequate references to study on other plants with close relationship to C. angulatus.

Terpenoid backbone biosynthetic genes and their differential expression patterns in C. angulatus
In previous studies on C. angulatus, Celangulin V were shown to be distributed in many organs, but mainly accumulated in root bark, the contents of Celangulin V showed higher in root than in leaf, and changed in different developmental stages [5][6][7]. We speculated that the genes expressed higher in root might encode some enzymes responsible for the biosynthesis of Celangulin V. Therefore, choosing the root and the leaf for comparative transcriptome analysis will greatly facilitate dissection of the genes involved in the organ-specific biosynthesis of Celangulin V. This approach is widely used for mining and identifying novel genes in biosynthesis of secondary metabolites in plants [35][36][37]. IPP and its isomer DMAPP, the universal biological precursors of all isoprenoids (basic C5 isoprene unit), can be obtained either through MVA/MEP pathways [38] (Fig. 1) in all eukaryotic cells and cytoplasm and mitochondria of plants or MEP pathway in bacteria, other prokaryotes and plastids in plants [39].
Totally, 22 up-regulated unigenes related to MVA pathway and 38 up-regulated unigenes related to MEP pathway possibly involved in terpenoid backbone biosynthesis are listed in Table 5. Numerous studies have been conducted with regard to engineer microorganisms for the production of different isoprenoids. The yeast Saccharomyces cerevisiae was always chosen as a microbial host to enhance the production of sesquiterpene by deregulating the MVA pathway, and the Escherichia coli was usually chosen as the microbial host for deregulating the MEP pathway. For example, the production of artemisinic acid engineered by MVA pathway achieved nearly 500 times production as previously reported [40]. Other natural products, such as patchoulol, farnesol, limonene and lycopene have also got higher titers through MEP or MVA metabolic pathway engineering methods [38,41,42]. In this study, the 22 up-regulated genes involved in MVA pathway and 38 up-regulated genes involved in MEP pathway were the candidate genes that can enhance the availability of the precursor (FPP) for biosynthesis of sesquiterpenes, and have the potential to improve the production of Celangulin V.  Sesquiterpene biosynthetic genes and their differential expression patterns in C. angulatus Terpenes, the largest class of natural products, are formed by terpene cyclases from linear oligoprenyl diphosphate precursors [43]. The phylogentic analyses of plant TPS protein sequences recognized seven major clades or subfamilies (TPS-a, TPS-b, TPS-c, TPS-d, TPS-e/f, TPS-g and TPS-h) [44,45]. C. angulatus belongs to angiosperms. TPSs involved in angiosperm sesquiterpene synthesis are mainly distributed in TPS-a subfamily, with a small part dispersed in TPS-g/h subfamily. However, those involved in angiosperm monoterpene synthesis are mainly distributed in TPS-b subfamily [46]. CL7773.Contig1_All, CL8776.Conti-g2_All and CL8776.Contig1_All fell into TPS-a subfamily according to the phylogenetic analysis. Therefore, these unigenes maybe the sesquiterpene synthases and are predicted to be germacrene D synthase and germacradienol synthase [47,48]. The CL7773.Contig8_All is predicted to be valencene/7-epi-alpha-selinene synthase, and the recombinant enzyme from Vitis vinifera gave 49.5% (+)-valencene and 35.5% (−)-7-epi-α-selinene [47]. These four unigenes maybe involved in the Celangulin V biosynthesis (Table 5). Many studies have proposed the cyclization mechanisms of sesquiterpene synthases [49,50]. The different cyclization processes can form diverse structures [51]. FPP cyclization reactions can proceed through C1-C6, C1-C7, C1-C10 and C1-C11 bond formation depending on which carbon-carbon double bond reacts with the initially formed allylic carbocation [52]. Celangulin V is a sesquiterpene polyol ester with β-dihydroagarofuran skeleton. Based on the chemical structural characteristics of β-dihydroagarofuran, we speculated that these four unigenes encoding sesquiterpene synthases could convert FPP to β-dihydroagarofuran through C1-C10 bond formation. In the future research, we will identify the function of candidate genes for sesquiterpene synthases belonging to the Celangulin V biosynthetic pathway by experiment.
Celangulin V biosynthesis modified genes and their differential expression patterns in C. angulatus Plant CYP450s [53][54][55] and acyltransferases [56] modify the sesquiterpene carbon ring and are important in biotechnology due to their ability to biosynthesize diverse secondary metabolites with various biological properties. The CYP71 clan is the largest P450 family and functionally diverse in all plant species [54], including a total of 110 members divided into 18 subfamilies and showing large clusters of duplicated genes and taxa-specific subfamily blooms [57]. The CYP71 clan is also called the cradle of monoterpenoid and sesquiterpenoid diversity. P450s involved in sesquiterpenoids metabolism have spawned exclusively in five families of the CYP71 clan: CYP71AV, CYP71BA, CYP71BL, CYP71D, CYP706B [58]. So far, the functions of many CYP450s have been identified, participating in the biosynthesis of terpenoid-based natural products in various plants [54]. At present, six bioactive sesquiterpenes have been successfully produced in heterologous systems by using CYP450 coding sequences [53]. In the capsidiol biosynthesis, 5-epi-aristolochene 1,3-dihydroxylase (CYP71D20) played a decisive role and the preferred reaction order of hydroxylation is at C1 followed by C3 [59]; CYP706B1 (part of the CYP71 clan) catalyzed hydroxylation of δ-cadinene at C8 in gossypol's biosynthesis [60]; CYP71BA1 from Zingiber zerumbet (Zingiberaceae) catalyzed the conversion of a-humulene to 8-hydroxy-a-humulene at C8 in the zerumbone biosynthesis [61]; CYP71BL2 from Lactuca sativa [62] and CYP71BL3 from Cichorium intybus [63] specifically catalyzed 6a-hydroxylation by mediate hydroxylation of germacrene A acid at positions adjacent to the carboxy group. We obtained eight unigenes that were predicted to be the CYP71D9, CYP71D10, CYP71D11 may involve in Celangulin V biosynthesis (Table 5). It is not well documented that acyltransferases are recruited as catalysts in the biosynthesis of plant specialized sesquiterpenoid compounds. Take artemisinin [40] and capsidiol [64] for examples, the acylated groups are not needed in their biosynthesis. However, some plant terpenoids need an acylation post-modification, which is depend on BAHD superfamily acyltransferases [65]. The BAHD acyltransferases are named after the first four biochemically characterized enzymes of the group, which use acyl-CoA thioesters as donor molecules [66]. For example, two BAHD alcohol acetyltransferases were identified from L. x intermedia glandular trichome database, which convert a variety of monoterpenes to monoterpene esters using coenzyme A as a cofactor [67]. Several taxoid-O-acetyl transferases (like TAX 9 and TAX 14), obtained from a previously isolated family of Taxus acyl/ aroyl transferase cDNA clones, could catalyze some acylation steps of Taxol biosynthesis [68,69]. We obtained four unigenes that were predicted to be the BAHD acyltransferase from transcriptome data, probably participating in biosynthesis of Celangulin V, according to acylated groups on C5, C6, C8, C9 and C14 positions of β-dihydroagarofuran skeleton.
Accordingly, we proposed the biosynthetic pathway of Celangulin V: β-dihydroagarofuran is the parent skeleton, highly oxygenated at C1, C3, C5, C6, C8, C9, C14 sites by one or more CYP450s to produce dihydroagarofuran polyols; Then, C5, C6 hydroxyl groups are acetylated by acetyltransferases; C9, C14 hydroxyl groups are butyrylated by butyryltransferase; C8 hydroxyl group is benzoylated by benzoyl transferase; C1, C3 are free. Based on our transcriptome sequencing, we annotated eight up-regulated CYP450 unigenes and four up-regulated BAHD acyltransferases in the sesquiterpene biosynthesis. In the future research, we will identify the function of candidate CYP450 and acyltransferase genes.

Conclusions
Celangulin V is an efficient insecticidal sesquiterpenoid and widely used in pest control in China. In this study, leaf and root transcriptome were generated by next -generation sequencing, obtaining 104,950 unigenes with Fig. 7 qRT-PCR validation of selected 15 DEGs at Leaf and Root. The relative expression level of each selected gene was determined by 2-ΔΔCT. Each bar represents the mean ± STD of triplicate assays. Values with different letters indicate significant differences at P < 0.05 according to Duncan's multiple range tests. This figure was generated with OriginPro 8.5.1SR2 average length of 1200 bp. The transcriptome analysis revealed 16 unigenes probably involved in Celangulin V biosynthesis, in which four unigenes encoded sesquiterpene synthases (germacrene D synthase, germacradienol synthase, valencene synthase and 7-epi-alpha-selinene synthase), eight unigenes encoded CYP450s (CYP71D9, CYP71D10, CYP71D11) and four unigenes encoded BAHD acyltransferases. We further proposed the Celangulin V biosynthetic pathway. The transcriptome data not only provides valuable information on sesquiterpene polyol ester with β-dihydroagarofuran skeleton but also assists in mining genes involved in the production of other plants with close relationship to C. angulatus.