Mining of candidate genes involved in the biosynthesis of dextrorotatory borneol in Cinnamomum burmannii by transcriptomic analysis on three chemotypes

Plant material, chemicals and reagents

Ten-year-old Cinnamomum burmannii of low-content borneol type (LBT, 0 <D-borneol <1 mg/g) and borneol-free type (FBT) were planted in the Medicinal Botanical Garden of Guangzhou University of Chinese Medicine. Samples of ten-year-old C. burmannii of high-content borneol type (HBT; D-borneol ≥ 1 mg/g) samples were donated by Guangdong Huaqingyuan Biotechnology Co. LTD. Authentic standard D-borneol (purity >98%) was purchased from Shanghai Yuanye Bio-Technology Co., Ltd.

Gas chromatograohy-mass spectrometry (GC/MS) analysis of chemical contents

Leaves collected from different individual trees of C. burmannii were ground into fine powder in liquid nitrogen. Powder samples (0.05–0.10 g) were weighed accurately, soaked in 2.0 mL petroleum ether using an ultrasonic cleaner for 30 min, and centrifuged at 16,000 rpm for 5 min in preparation for GC/MS analysis. Three technical replicates were performed for each of the three biological samples. An Agilent 7890B Gas Chromatograph with 5977A Inert Mass Selective Detector (Agilent, USA) was used to analyze the extraction. Helium was used as the carrier gas with a flow rate of 1 mL/min. A Cyclosil-B GC column (30 m × 0.25 mm × 0.25 µm) with initial temperature of 50 °C was used to separate samples. The GC column temperature program was as follows: oven temperature was increased from 50 °C to 180 °C at 2 °C/min, and then from 180 °C to 300 °C at 4 °C/min, then held at 300 °C for 4 min. Injector: split mode with split ratio of 20:1, injection volume of 1.0 µL and inlet temperature of 250 °C. D-Borneol was confirmed by comparing its retention time with that of the known standard, which was determined under the same conditions. The concentration of D-borneol in the plant sample was calculated from its respective standard curves.

RNA extraction, cDNA library preparation and sequencing

Total RNA of samples collected from those three groups was extracted using TRizol reagent according to the manufacturer’s instructions and then was digested with RNase-free DNase to eliminate the contamination of genomic DNA. The integrity and quantity of the RNA samples were analyzed using 1% agarose gel electrophoresis and with a NanoDrop 2000 Spectrophotometer (Thermo Scientific, USA). Qualified RNA samples (OD₂₆₀/₂₈₀ = 1.8∼2.2, OD₂₆₀/₂₃₀ ≥ 2.0, RIN ≥ 6.5, 28S:18S ≥ 1.0) were immediately frozen at −80 °C until use. To analyze the transcriptome, 2 µg qualified RNA samples from three biological replicates were mixed and sent to Majorbio (Shanghai, China) to construct a cDNA library on an Illumina Hiseq X Ten (Illumina Inc., San Diego, CA, USA). Paired-end (PE) reads were generated and checked by the software fastx_toolkit_0.0.14 (http://hannonlab.cshl.edu/fastx_toolkit/) to assess the quality of sequences.

De novo transcriptome assembly and annotation

High-quality sequencing data was de novo assembled using TRINITY v2.5.0 (https://github.com/trinityrnaseq/trinityrnaseq) with default parameters, which integrates three independent software modules (Inchworm, Chrysalis and Butterfly) to process and splice a large amount of RNA-Seq data in turn. First, the Inchworm software was used to assemble the clean data into linear contigs. Second, the overlapping contigs were clustered into sets of connected components using the Chrysalis software. Finally, transcripts were constructed using the Butterfly software. The assembly results were further filtered and optimized using TransRate (http://hibberdlab.com/transrate/) and clustered to obtain non-redundant unigenes using CD-HIT (http://weizhongli-lab.org/cd-hit/). For gene identification and expression analysis, the reads from different chemotypes were co-assembled, and for gene sequence analysis, reads from different chemotypes were assembled separately.

Annotation of unigenes was performed using the BLAST program with an E-value cut-off of 1E-5 against the following three databases: NCBI non-redundant protein sequences (Nr, https://blast.ncbi.nlm.nih.gov/), Protein family (Pfam, https://pfam.sanger.ac.uk/) and Swiss-Prot (https://web.expasy.org/docs/swiss-protguideline. html). Functional annotation using gene ontology (GO) terms was analyzed using Blast2GO version 2.5.0 (http://www.geneontology.org/). Kyoto encyclopedia of genes and genomes (KEGG) pathways analysis was performed using KOBAS version 2.1.1 (https://www.genome.jp/kegg/) for systematic analysis of gene function.

Differential gene expression analysis

Fragments Per Kilobases per Million-reads (FPKM) was used to estimate relative expression levels of transcripts using the software RSEM v1.2.15 (http://deweylab.github.io/RSEM/) with default parameters. Raw counts were directly analyzed statistically using DESeq2 version 1.10.1 software based on the negative binomial distribution. Unigenes with expression differences among the groups were obtained with the parameters p-adjust <0.05 and —log2FC— ≥ 1. KEGG enrichment analysis of differentially expressed genes (DEGs) was performed using KOBAS version 2.1.1 (http://www.genome.jp/kegg/).

Screening of candidate TPS genes and transcription factors related to D-borneol biosynthesis

Candidate TPS genes were first identified by their KEGG annotation and selected by their length and RPKM expression value. Up-regulated genes were then selected for further analysis. The coding sequences of the candidate TPS genes were translated into amino acid sequences, and a phylogenetic tree was reconstructed from the aligned sequences using the neighbor-joining algorithm in MEGA6.0 (Tamura et al., 2013) with 1,000 bootstrap replicates. Transcription factors were analyzed using the HMMER method (Thiriet-Rupert et al., 2016) and their sequences aligned using the plantTFDB (http://planttfdb.cbi.pku.edu.cn/); they were divided into different families according to their conserved domains.

Validation of DEGs by qRT-PCR analysis

Quantitative reverse-transcription polymerase chain reaction (qRT-PCR) was conducted using TransStar Tip Green qPCR SuperMix (TransGen Biotech, Beijing, China) and the CFX96 Touch Deep Well platform (Bio-Rad, USA) with a total reaction volume of 20 µL, comprising 1 µL cDNA, 10 µL Tip Green qPCR Supermix, 0.4 µL 10 µM each primer and 8.2 µL nuclease-free water. Nine differentially expressed unigenes related to monoterpenenoid biosynthesis and six candidate transcription factors related to D-borneol accumulation were selected, including CbTPS2 (c122670.graph_c0) and CbTPS3 (c129067.graph_c0), CbGPPS_5 (c117157.graph_c0), CbDXR (c84883.graph_c0), CbDXS_1 (c126058. graph_c0), CbDXS_3 (c134634.graph_c0), CbHDS_1 (c134509.graph_c0), CbHDR (c121089.graph_c0), CbHMGR_1 (c133256.graph_c0), CbAP18/ERF-1 (c130684.graph_c1), CbbHLH-1 (c117041.graph_c0), CbbHLH-2 (c127743.graph_c0), CbMYB-1 (c120253.graph_c0), CbMYB-2 (c123637.graph_c0), CbWRKY-1 (c122191.graph_c0). The β-tubulin gene (c121493.graph_c0) was used as a reference. Primer premier 5 software (Premier Biosoft Intl., CA USA) was used to design primers for the selected unigenes; coding sequences and primers are listed in Data S1–S2. The qRT-PCR procedure was as follows: 94 °C for 5 min, 40–45 cycles of 94 °C for 5 s, 60 °C for 15 s and 72 °C for 10 s, followed by a dissociation stage. Expression levels were evaluated using the 2^−△△Ct method, and the Ct values for all genes were normalized to the Ct value of β-tubulin.

D-borneol quantitative determination in three chemotypes of C. burmanni

In previous work, we collected three different chemotypes (HBT, LBT, BFT) of C. burmannii exhibiting similar morphology (Figs. 2A–2C). These chemotypes possessed high, low and absent levels of D-borneol, respectively. Here, we used a GC/MS method to analyze the composition of volatile oil in these C. burmannii samples (Figs. 2D–2J). A total of 28 compounds were detected in the three chemotype samples, of which nine were common to all three. D-borneol was the most abundant component of the volatile terpenoids in both HBT and LBT groups, accounting for 42.23% and 57.26%, respectively. In the BFT group, the main component was cymene, accounting for 27.44% of the volatile terpenoids (Figs. 2K–2N, Data S3–S4).

Transcriptome analysis

To identify genes related to the biosynthesis of D-borneol, we performed deep transcriptome sequencing of the three C. burmanni chemotypes using an Illumina HiSeq 2000 platform. In total, 21.21 Gb of clean sequence data with GC content of 47.55% was obtained from all samples after filtering and removing adapter sequences from the raw data: 8.40 Gb from the HBT group, 6.39 Gb from the LBT group and 6.42 Gb from the BFT group. Raw data and quality control data statistics are summarized in Table 1. All RNA-Seq raw data has been submitted to CNGBdb with accession number CNP0000810.

For further analysis, high-quality clean reads were assembled using the TRINITY program, resulting in 142,673 transcripts and 100,218 unigenes with N50 of 1,625 bp and 1,128 bp and average length of 935 bp and 713 bp, respectively. The majority of the reads for both transcripts and unigenes were shorter than 1,000 bp in length, accounting for 92% and 88%, respectively. Only 2% of the transcripts and 3% of the unigenes were longer than 2,000 bp in length. All assembled data are summarized in Table 2 and Fig. 3A. Moreover, 80.80% of the clean reads were perfectly mapped to the assembled unigenes by RSEM, indicating that these mapped genes were of high quality and could be used to conduct the subsequent analysis (Table 1). Sequence of the 100,218 unigenes were shown in Data S5.

All resulting 100,218 unigenes were used for BLAST searches against the GO, KEGG, NR, Pfam and Swiss-Prot databases for annotation, resulting in 28,798, 15,376, 44,495, 28,107 and 24,238 annotated unigenes for the above database, respectively (Fig. 3B, Data S6). In total, 45,485 unigenes were annotated to at least one database, and 9,730 unigenes shared annotation in all five databases (Fig. 3C). The species distribution of the annotated unigenes is shown in Fig. 3D.

Identification of DEGs

A total of 8,860 DEGs were identified among the comparison groups. The expression patterns of these DEGs could be classified into four clusters (Figs. 3E–3F, Data S7). Cluster 1 consisted of 2,600 genes significantly down-regulated in the LBT group in comparison with both the HBT and BFT groups. Cluster 2 exhibited the opposite expression pattern to that of cluster 1, containing 3,721 unigenes that were particularly up-regulated in the LBT group in comparison to both the HBT and BFT groups. Cluster 3 contained 2,539 unigenes significantly up-regulated in the BFT group in comparison to both the HBT and LBT groups. Ten genes significantly down-regulated in the BFT group in comparison to both the HBT and LBT groups formed cluster 4. We used a Wayne diagram (Fig. 3G) to show the number of DEGs between each group. There were 3,140 up-regulated unigenes and 3,643 down-regulated unigenes in the HBT vs LBT comparison. In the HBT vs BFT comparison, we identified 1,303 up-regulated unigenes and 924 down-regulated unigenes. In the LBT vs BFT comparison, 3,902 unigenes were up-regulated while 2,357 unigenes were down-regulated (Fig. 3H).

To gain a deep understanding of the biological functions of these 8,860 DEGs, we performed KEGG pathway enrichment analysis. We used dot plots to show the top 20 most enriched KEGG pathway categories with p value <0.05, which indicated that the top three most enriched KEGG pathways were “cutin, suberin and wax biosynthesis,” “vancomycin resistance” and “fatty acid degradation” (Fig. 4A, Data S8). The top 10 most enriched GO terms under the three general sections Biological Process (BP), Cellular Component (CC) and Molecular Function (MF) are summarized in Fig. 4B and Data S9. The most enriched GO term in the MF category was protein serine/threonine kinase activity. Within the CC category, the most abundant GO term was chloroplast. In the BP category, the most associated term was protein phosphorylation. In addition, terpene synthase activity was also significantly enriched.

Mining of candidate genes related to D-borneol biosynthesis

Monoterpenes derived from plants are biosynthesized in the plastid through the MEP pathway. We detected 30 unigenes involved in the MEP pathway, but only three of them exhibited significantly different expression levels among the three comparison groups. Most of the unigenes (11/30), including one DEG, were annotated as DXS genes, which encode the first rate-limiting enzyme in the MEP pathway. Relatively higher expression of the DXS gene in the HBT group compared to that in the other two groups might contribute to the differences in D-borneol accumulation in C. burmannii. Six unigenes were annotated as encoding GPPS, responsible for GPP formation. It is well known that GPPSs exist as homomeric and heteromeric structures in plants (Burke, Klettke & Wise, 2004; Rai et al., 2013; Wang & Dixon, 2009; Yin et al., 2017) and that plants contain multiple GGPPS or GGPPS-related enzymes, some functioning as large or small subunits (LSU/SSU) of heteromeric GPPS (Tholl et al., 2004; Wang & Dixon, 2009). For example, there are 12 GGPPS paralogs in the Arabidopsis genome (Beck et al., 2013; Lange & Ghassemian, 2003); however, AtGPPS11 and AtGGPPS12 were characterized as LSU and SSU, with their interaction involved in monoterpene biosynthesis in Arabidopsis flowers (Chen, Fan & Wang, 2015). To confirm whether the C. burmannii GPPSs formed homomeric or heteromeric structures, we conducted a phylogenetic analysis of seven KEGG-annotated CbGPPS/GGPPS proteins from C. burmannii and the amino acid sequences of other characterized GPPS/GGPPS proteins. We identified six GPPSs, including two putative homomeric GPPSs, three putative GPPS.SSUs and a putative GPPS.LSU. Only one unigene might encode a GGPPS (Fig. S1). It has been reported that monoterpene biosynthesis in gymnosperms is mainly regulated by homomeric GPPS and SSUs of heteromeric GPPS (Schmidt et al., 2010; Tholl et al., 2004). The 2 putative homomeric GPPSs, 3 putative GPPS.SSUs identified in this study may provide a possible explanation for the monoterpene richness of C. burmanni.

There were 40 unigenes annotated as TPS by KEGG, among which 11 were associated with the sesquiterpenoid and triterpenoid biosynthesis pathway (ko00904), eight were monoterpenoid synthase unigenes (ko00902) and 21 encoded diterpenoid synthase (ko00909). Since D-borneol is a monoterpene, the eight putative monoterpenoid synthase genes were selected for further analysis. The expression patterns of these eight unigenes encoding monoterpenoid synthase, along with those of genes related to the MEP metabolic pathway (Ko00900), were normalized and shown in Figs. 5A–5B and Data S10. Considering the relatively high accumulation of D-borneol in the HBT group compared to the other two groups, we rationalized that metabolism-related genes might exhibit more activity in the HBT group than in the other two groups. Inspiringly, two unigenes encoding monoterpenoid synthase (designated CbTPS2 and CbTPS3) were significantly up-regulated in the HBT group compared with the LBT and BFT groups (Data S11). Furthermore, both unigenes had complete open reading frames. Sequence comparison using BLAST indicated that CbTPS2 (c122670.graph_c0) exhibited highest similarity [53.27% amino acid (aa) sequence identity] to an (E)-beta-ocimene/myrcene synthase (ADR74206.1) from Vitis vinifera, while CbTPS3 (c129067.graph_c0) presented similarity in aa sequence with a S-(+)-linalool synthase (RVW86370.1) from Vitis vinifera (Martin et al., 2010).

These two candidate TPS genes were subjected to a phylogenetic analysis, along with other TPS genes that have been functionally characterized. Similar to other well-characterized BPPS genes, CbTPS2 fell within the TPS-b subfamily (Fig. S2). Furthermore, sequence alignment showed that CbTPS2 had conserved mono-TPS domains, including RRX8W, RXR, DDXXD and (N,D) D (L,I,V) X (S,T) XXE (NSE/DTE) (Fig. S3). To find out whether the candidate unigenes differed between each chemotype, the transcriptomes of the three groups were assembled separately, and the sequences were compared. The assembled sequences of each chemotypes were shown in Data S6–S14. The amino acid sequences of CbTPS2 in the three different C. burmannii chemotypes differed only by two amino acids, and these sites were not located within the conserved domain, suggesting that these differences may not affect the function of these proteins (Fig. S4). However, it is interesting that the complete sequence of CbTPS3 could only be found in the HBT and LBT groups, whereas only a small fragment was found in the BFT group (Fig. S5). This may be due to poor assembly of the BFT transcriptome, or perhaps CbTPS3 is a pseudogene in BFT. We speculated that CbTPS2 and CbTPS3 are probably candidate genes responsible for the biosynthesis of D-borneol in C. burmannii.

Mining of candidate transcription factors related to D-borneol accumulation

Transcription factors play an important role in regulating the biosynthesis and accumulation of secondary metabolites (Hong et al., 2012). In our transcriptome analysis, we identified 1,108 transcription factors, which could be assigned to 47 families (Data S15). Several TF families, including MYB, basic helix-loop-helix (bHLH), WRKY, basic region/leucine zipper motif (bZIP) and APETALA2/ethylene-response factor (AP2/ERF) families, function in regulating the structural genes involved in monoterpenoid biosynthesis (Su et al., 2019). These five TF families cover most of the TFs identified in this study (Fig. 5C). We identified 116 MYB, 65 AP2/ERF, 63 BHLH, 60 bZIP and 47 WRKY candidate genes, among which 35 MYB, 27 AP2/ERF, 24 BHLH, 10 bZIP and 9 WRKY genes were classified as DEGs (Figs. 5D–5H, Data S16). Notably, one WRKY (c122191.graph.c0), two BHLH (c117041.graph.c0, c127743.graph.c0), one AP2/ERF (c130684.graph.c1) and three MYB candidate genes (c68268.graph.c0, c120253. graph.c0, c123637.graph.c0) were significantly up-regulated in both the HBT and LBT groups compared to the BFT group. The different transcription factor expression patterns might be responsible for regulating monoterpenoid biosynthesis. Moreover, the up-regulated transcription factors in the HBT and LBT groups compared to the BFT group provide us with a candidate TF pool related to the biosynthesis and accumulation of D-borneol.

Validation of DEGs using Quantitative real-time PCR (qRT-PCR)

To verify the accuracy of RNA-Seq data, we performed expression level analysis on nine unigenes related to the terpenoid biosynthesis pathway and six candidate transcription factors related to D-borneol accumulation using the qRT-PCR method. As shown in Fig. 6, their expression patterns were similar to those revealed from RNA-Seq data. Thus, these results help confirm the reliability of the transcriptomic analysis.