Metabolism and transcriptome profiling provides insight into the genes and transcription factors involved in monoterpene biosynthesis of borneol chemotype of Cinnamomum camphora induced by mechanical damage

Plant material, chemicals, and reagents

We purchased the borneol chemotype of C. camphora from the Ji’an Yu Feng Natural Species Co., Ltd (Ji’an, China), and these plants were grown in an artificial climate box (Shanghai Yiheng Instrument Co. Ltd, Shanghai, China), at 25 °C, with a 12-h-light/12-h-dark photoperiod. The D-borneol compound was bought from the Shanghai Yuanye Bio-Technology Co., Ltd (Shanghai, China).

Mechanical damage (MD) treatment

Nine plants showing uniform growth and no disease symptoms were divided into three groups, each containing three biological samples. A stress treatment in the form of MD was applied to six of the nine BCC plants, by referring to two studies (Reymond et al., 2000; Ross, Kupper & Jacobs, 2006). Specifically, we used a medical hemostatic button to damage about 40% of each leaf’s area on a plant. These damaged leaves were then harvested at 2 h and 6 h (respectively, MD_2h and MD_6h) after wounding, leaving three undamaged BCC plants to serve as the control group (CK). We collected all their leaves and immediately froze them using liquid nitrogen and stored them at −80 °C until their use.

Extraction of chemical contents

After imposing the MD treatment, the leaves collected at different times were quickly placed in liquid nitrogen and ground into a powder. From each sample, 0.15 g of powder was accurately weighed and put into an Eppendorf tube with 2.0 mL of petroleum ether. After 30 min of sonication (using an ultrasonic cleaner), the samples were centrifuged at 12,000 rpm for 10 min prior to their Gas Chromatography-Mass Spectrometer (GC/MS) analysis. Three technical replicates were performed for each of the three biological replicate samples.

GC-MS analysis of chemical contents

To analyze the leaf extracts, we used an Agilent 7890B gas chromatograph equipped with a 5977A inert mass selective detector (Agilent, USA). Helium was the carrier gas, at a flow rate of one mL/min, while a Cyclosil-B GC column (30 m × 0.25 mm × 0.25 µm) at an initial temperature of 50 °C separated the analytes. The GC column’s temperature program consisted of an oven temperature rising from 50 °C to 180 °C at a rate of 2 °C/min, and then from 180 °C to 300 °C at 4 °C/min, then holding it at 300 °C for 4 min. For the injector, it was used in the split mode with a split ratio of 20:1, with an injection volume of 1.0 µL and an inlet temperature of 250 °C. To identify the ensuing components, their recorded mass spectra were compared to existing data stored in the NIST14/Wiley275 Mass Spectral Library. Further, we confirmed the presence of D-borneol by comparing its retention time to that of a known standard, this determined under the same GC/MS conditions above. The concentration of D-borneol per plant sample was calculated from its standard curve, while the contents of other components were quantified using D-borneol as an external standard. To determine the proportion of each component, its peak area/total peak area ratio was calculated.

RNA extraction, cDNA library preparation, and sequencing

To extract the total RNA of the samples collected at each time point after the MD treatment, we used the EASYspin Plus Plant RNA kit (Aidlab, Beijing, China). To check the integrity and quantity of the RNA samples, 1% agarose gel electrophoresis and a NanoDrop 2000C Spectrophotometer (Thermo Scientific, USA) were respectively used. To qualify for the analyses, the RNA samples had to meet these criteria: an OD 260/280 = 1.8∼2.2, OD 260 / 230 ≥ 2.0, an RIN ≥ 6.5, and a 28S:18S ≥ 1.0; these were sent to Majorbio (Shanghai, China) to build the cDNA library after which transcriptome sequencing was carried out on an Illumina HiSeq 4000 platform (Illumina Inc., San Diego, CA, USA). This generated paired-end (PE) reads whose sequence quality was initially assessed using fastx_toolkit_0.0.14 software (http://hannonlab.cshl.edu/fastx_toolkit/). To clean the raw reads, we first removed any adaptor sequences, empty reads, and low quality sequences (i.e., those reads with unknown sequences ‘N’ or less than 25 bp in length), using SeqPrep (https://github.com/jstjohn/SeqPrep) and Sickle software (https://github.com/najoshi/sickle) set to their default parameters.

De novo transcriptome assembly and annotation

Next, this high quality-sequenced clean data set underwent de novo assembly in TRINITY v2.5.0 software (https://github.com/trinityrnaseq/trinityrnaseq) set to default parameters, which integrates three independent software modules (Inchworm, Chrysalis and Butterfly) to process and splice large amounts of RNA-Seq data in turn. First, the Inchworm tool assembles the clean data into linear contigs; next, Chrysalis clusters the overlapping contigs into sets of connected components; finally, the Butterfly tool constructs the transcripts. The assembly results were then further filtered and optimized by TransRate (http://hibberdlab.com/transrate/), as well as clustered, to obtain the non-redundant unigenes via CD-HIT (http://weizhongli-lab.org/cd-hit/). For this assembly, we used BUSCO (Benchmarking Universal Single-Copy Orthologs, http://busco.ezlab.org) (Simao et al., 2015), set to its default parameters, to estimate the completeness of a conserved set of genes.

To annotate the transcripts, we used the BLAST program with an E-value cut-off of 1E⁻⁵, running against four public databases: NCBI non-redundant protein sequences (Nr, https://blast.ncbi.nlm.nih.gov/), Protein family (Pfam, http://pfam.sanger.ac.uk/), Swiss-Prot (http://web.expasy.org/docs/swiss-protguideline.html), and Clusters of Orthologous Groups of proteins (COG, http://www.ncbi.nlm.nih.gov/COG/). Functional annotation according to the Gene Ontology (GO) terms was implemented in Blast2GO v2.5.0 (http://www.geneontology.org). We also performed a Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways’ analysis, in KOBAS v2.1.1 (http://www.genome.jp/kegg/), to systematically explore the functions of the genes. Their transcripts were also blasted with the PlantTFDB (http://planttfdb.cbi.pku.edu.cn/), using the HMMER method, to obtain the homologous TFs and their family information.

Differential gene expression analysis

To estimate the relative expression levels, the FPKM (fragments per kilo bases per million reads) of transcripts were determined in RSEM v1.2.15 software (http://deweylab.github.io/RSEM/) set to its default parameters. The raw counts were statistically analyzed directly, in DESeq2 v1.10.1 software, using a negative binomial distribution. Those transcripts having a significantly different expression between pairs of treatment groups were obtained by using threshold parameters of a P-adjusted <0.05 and —log2FC— ≥ 1. Next, these designated differentially expressed genes (DEGs) underwent a KEGG enrichment analysis in KOBAS v2.1.1 (http://www.genome.jp/kegg/).

Screening of candidate TPS genes and TFs related to monoterpene biosynthesis

To screen out candidate TPS genes, we first selected those unigenes directly annotated as TPSs by KEGG. Then, a local BLASTP analysis for candidate TPSs sequences was carried out against the unigenes, using previously reported related sequences from other plants. The selected unigenes had an identity >40%, bit score >100, and a sequence length >800 bp. Furthermore, to locally search for unigenes in the Pfam-A database, we used hmmsearch in HMMER 3.0 for the two Pfam domains (PF01397 and PF03936), with an e-value cut-off of 1 × 10⁻⁵. To not overlook any TPS candidates, the TERZYME program was relied upon to identify and classify TPSs (Priya et al., 2018). Finally, we integrated results of the KEGG annotation, BLASTP, HMM and TERZYME searches, and all candidate TPS genes were further examined to confirm the presence of both PF01397 and PF03936 domains using the Pfam database and the NCBI’s conserved domains database (CDD). We subjected the candidate CcTPSs to a phylogenetic analysis—this done in MEGA6.0, using the maximum likelihood (ML) method with n = 1000 bootstrap replicates and best models of JTT+G. The resulting tree’s presentation was augmented in EvolView v3 (Subramanian et al., 2019). The conserved motifs of CcTPSs were discovered using MEME (http://meme-suite.org/tools/meme). Further, to obtain the Spearman correlation coefficients between terpene synthases and volatile terpenoids, we used IBM SPSS 19.0 (Armonk, NY, United States). Finally, using Cytoscape (Shannon et al., 2003), a gene–metabolite network map was constructed according to those Spearman results (using only correlation coefficients >0.6, P-values <0.05).

To uncover those TFs associated with monoterpnene biosynthesis, we tested the differentially expressed AP2/ERF-, WRKY-, bZIP-, BHLH-, MYB- family genes for their correlation with the candidate monoterpene synthase. This Spearman correlation analysis was done with expression data for TF genes and candidate monoterpene synthase (also in IBM SPSS 19.0). For those Spearman correlation coefficients >0.8 with a P-value <0.05, a graphic of the putative TPSs–TFs networks were constructed in Cytoscape. The sequences of those TFs positively correlated with monoterpene synthases were submitted online to the website ORF finder (https://www.ncbi.nlm.nih.gov/orffinder/), to determine the presence of full-length sequences. Then, the deduced amino sequences were submitted to the online website of PfamScan (https://www.ebi.ac.uk/Tools/pfa/pfamscan/), to determine the presence of conserved domains.

Gene expression analysis by quantitative real-time PCR (qRT-PCR)

To verify the reliability of the obtained RNA sequencing data, we relied on qRT-PCR, as previously described (Yang et al., 2020). The expression levels of 23 DEGs—9 related to the biosynthesis of monoterpenes, 6 TFs, and 8 related to the jasmonic acid (JA) signaling pathway—were evaluated using the 2^−△△Ct method, with CcActin serving as the reference gene. The measurements from three biological replicates were performed in triplicates. All primers of the 23 unigenes, designed with Primer Premier 5.0 software, can be found in Table S1.

Mechanical damage (MD) increased the production of terpenes

Higher contents of terpenes in MD-treated plants relative to the CK group were evidenced except for one monoterpene and one sesquiterpene. Specifically, from BCC leaves we detected a total of 21 terpenes induced by MD treatment, including 14 monoterpenes (accounting for 88%, of which 63% was D-borneol) and 7 sesquiterpenes (accounting for 11%), using GC/MS (Fig. 1A, 1B, Table S2). Seventeen terpenes (11 monoterpenes, 6 sesquiterpenes) had the highest content 6 h after the MD treatment, and for 2 monoterpenes their content peaked after 2 h and then declined by 6 h post-MD. Furthermore, the carene, pseudolimonene, camphene, sabinene, eucalyptol, D-borneol, terpineol, elemene, and germacrene B contents were all prominently enhanced in the MD_6h samples compared with the CK group, while their levels in the MD_2h samples also tended to surpass those of CK, albeit not significantly (Fig. 1C).

Transcriptome sequencing and annotation, expression analysis of unigenes

Based on the ability of MD to increase D-borneol and other monoterpenes’ production, deep transcriptome sequencing of the induced BCC collected at different hours was performed, as well as the CK group. In total, 67.66 Gb of high-quality base pair data (at least 7.51G per sample and 22.55 G for each group), with a GC content of 48.62%, were obtained (Table 1). Our PCA analysis of the resulting transcription data revealed a clear separation of the MD_6 h group from both MD_2 h and CK groups, suggesting a high repeatability of the sequencing data (Fig. 2A).

A total of 86,452 unigenes were assembled de novo, these having a mean length of 837 bp and an N₅₀ of 1,300 bp (Table 2). Most of the unigenes (65,483; 75.74%) were shorter than 1,000 bp, while only a few genes (7,308; 8.4%) were longer than 2,000 bp (Fig. 2B). We then evaluated the quality of the assembled sequences by examining the mapped ratio of the clean data to the assembled unigenes of each sample using RSEM. The result showed that 81.05% of the clean reads could be perfectly mapped to the reference unigenes, indicating the throughputs and sequencing quality were sufficient to ensure the reliability of further analyses (Table 1). Moreover, the ‘TransRate’ assembly score, which aims to assess the accuracy and completeness of transcriptome assembly, reached 31.45% for C. campora, a value higher than that of many of other assemblies (Smith-Unna et al., 2016). We also assessed the assembly completeness using BUSCO, this lent further support to the high quality of the assemblies, in which 80.70% of single-copy orthologs were successfully captured (Table 2). All RNA-Seq raw data was submitted to CNGBdb (accession number CNP0000956) and NCBI (accession number PRJNA679978). And the the de novo transcript assemblies was submitted to CNGBdb (accession number CNA0020891, https://db.cngb.org/search/project/CNP0000956/).

Using the BLAST algorithm, sequence annotations were performed against the NR, Swiss-Prot, Pfam, COG, GO, and KEGG databases. From this, 27,727 unigenes (32.07%) were annotated by the KEGG database, 25,714 unigenes (29.74%) were identified in the GO database, 13,356 unigenes (15.45%) were annotated against the COG database, and 34,072 unigenes (39.41%) had significant matches to the NR database. In addition, 41,035 (47.47%) and 41,466 (47.96%) of the unigenes displayed significant similarities to known proteins in the Pfam and Swiss-Prot database, respectively (Fig. 2C). When summed, we found 15,578, 22,906, and 26,732 of the unigenes expressed in the CK, MD_2h, and MD_6h groups, respectively (Fig. 2D). Among those, 1,582 unigenes were expressed only in the CK group, 5,441 specifically in the MD_2h group, and 9,472 unigenes uniquely so in the MD_6h group (Fig. 2E), indicating that the MD treatment was able to significantly induce the expression of genes in BCC. These results provided a general overview of the assembled transcriptome that we then used in our further analysis of DEGs.

Comparing the DEGs induced by MD treatment to the enrichment analysis

To fully explore potential DEGs induced by the MD treatment, we compared three groups of C. camphora plants. The DEGs were designated according to an expression level of —log 2 (fold-change) —>1 and FDR <0.05 in each pairwise comparison. This revealed 12,180 (9,643 upregulated, 2,537 downregulated) and 13,466 (9,222 upregulated, 4,244 downregulated) unigenes that were significantly differentially expressed relative to the control samples in the MD-treated samples at 2 h and 6 h post-treatment, respectively. Between MD_2h and MD_6h, we identified 10,477 unigenes significantly expressed (4,570 upregulated, 5,907 downregulated) (Figs. S1A–S1C). Among the three pairwise comparisons, 2,539 genes were commonly regulated. Finally, a total of 19,149 DEGs were screened among the three comparison groups (Fig. S1D).

To exhaustively explore the presumed functions of the above DEGs, we then performed KEGG and GO enrichment analyses to uncover overrepresented pathways and biological functions in the BCC leaves as induced by their mechanical damaging. Based on the GO analysis, within the main molecular function (MF), cellular component (CC), and biological process (BP) categories the most enriched GO terms of the upregulated DEGs in the CK vs. MD_2 h comparison were protein serine/threonine phosphatase activity, recognition of pollen, and intrinsic component of membrane (Table S3); in the CK vs. MD_6 h comparison, the over-presented GO terms contained genes were related to transcription, the DNA-templated extracellular region, and heme binding activity (Table S4). For the MD_2 h vs. MD_6 h comparison, the largest percentages of unigenes identified within the CC, BP, and MF categories were extracellular region, reactive oxygen species metabolic process, and monooxygenase activity (Table S5).

Additionally, KEGG pathway enrichment analyses showed that the upregulated DEGs from the CK vs. MD_2h comparison could be assigned to 128 pathways, while the top 3 enriched pathways were “plant-pathogen interaction”, “plant hormone signal transduction”, and “phenylpropanoid biosynthesis” (Fig. S1E, Table S6). In comparing the CK with the MD_6 h group, the most overrepresented pathways for the upregulated DEGs were “phenylpropanoid biosynthesis”, “alpha-linolenic acid metabolism”, “plant-pathogen interaction”, “phenylalanine, tyrosine, and tryptophan biosynthesis”, and “selenocompound metabolism” (Fig. S1F, Table S7). Regarding the MD_2h vs. MD_6h comparison, the KEGG pathways of “ribosome”, “phenylpropanoid biosynthesis”, “plant-pathogen interaction”, “MAPK signaling pathway-plant”, and “flavonoid biosynthesis” were the five most enriched (Fig. 3S1G, Table S8).

Identifying the DEGs related to terpene biosynthesis induced by MD

To further assess the effect of MD induction on volatile terpene biosynthesis, we next focused on the expression pattern of genes involved in the biosynthesis of terpenes. A total of 20 genes participating in the terpene backbone synthesis pathway had their activities significantly altered by MD. These 20 DEGs consisted of 3 homologues of HMGR, 3 homologues of HMGS, 5 homologues of DXS, 2 homologues of IDI, 1 homologues of GPPS, 2 homologues of FPPSs, and 4 homologues of GGPPSs. The expression of each DEGs was substantially increased in MD_2h or MD_6h when compared with the CK group, except for one putative DXS (CcDXS1) and one GGPPS (CcGGPPS2) (Fig. 3, Table S9).

We spotted 13 TPS-encoding unigenes (Table 3). According to our phylogenetic analysis based on the deduced amino acid sequences in combination with well-characterized terpene synthase genes, the 13 BCC TPSs could be divided into five subfamilies (Fig. 4A). The CcTPS1 (TRINITY_DN14566_c0_g1), CcTPS2 (TRINITY_DN50688_c1_g1), CcTPS3 (TRINITY_DN43620_c0_g4), and CcTPS5 (TRINITY_DN39759_c0_g1) all belong to the TPS-b subfamily, which are responsible for the synthesis of cyclic monoterpenes and contain the RR (X) 8W motif in their N-terminal region. Both CcTPS4 (TRINITY_DN45508_c0_g2) and CcTPS6 (TRINITY_DN49732_c0_g3) could be assigned to the TPS-g subfamily, which produces acyclic monterpenes and lacks the RRX8W motif in their N-terminal region. Another three, CcTPS7 (TRINITY_DN34440_c0_g1), CcTPS8 (TRINITY_DN40799_c0_g1), and CcTPS9 (TRINITY_DN46911_c0_g1), are members of the TPS-a subfamily clade, which may encode a sesquiterpene synthase. Next, CcTPS10 (TRINITY_DN47476_c0_g1), CcTPS11 (TRINITY_DN47476_c0_g2), and CcTPS13 (TRINITY_DN51000_c1_g1) belong to the TPS-e/f subfamily, while only CcTPS12 (TRINITY_DN49635_c3_g3) belongs to the TPS-c subfamily that is responsible for the synthesis of diterpenes. Using the MEME program, the conserved motif analysis revealed that those TPS members clustering in the same subfamily showed similar motif characteristics (Figs. 4B–4C). Of these 14 TPSs, 10 were defined as differentially expressed TPS genes, whose expression pattern were shown Fig. 3. We used these 10 putative TPSs, comprising 5 putative monoterpenoid synthase, 1 sesquiterpenoid synthase, and 4 diterpenoid synthase genes, in the further analyses.

The FPKM values of the 10 candidate TPS genes and the content of terpenes were integrated and a correlation analysis was performed to clarify the relationship between them. According to the ensuing correlation networks (Fig. 4D, Table S10), CcTPS1 was positively associated with six monoterpenes (carene, pseudolimonene, camphene, sabinene, L-limonene, D-limonene) as was CcTPS3 and CcTPS4 with four monoterpenes (borneol, bornyl acetate, camphor, terpineol). Presumably, CcTPS1 might directly influence the biosynthesis of carene, pseudolimonene, camphene, sabinene, L-limonene, and D-limonene in BCC leaf tissue, while CcTPS3 may partake in the biosynthesis of borneol, camphor, and terpineol, with CcTPS4 perhaps participating in the biosynthesis of borneol, terpineol, and bornyl acetate.

Identifying candidate transcription factor in response to MD

Since transcriptional reprogramming is primarily controlled by TFs, we further identified the differentially expressed TFs among the three treatment comparison groups. In this study, 673 differential expressed putative TFs belonging to 30 families were identified in BCC among the three group comparisons. There are about 58 TF families in higher plants, for which the major TF families, such as AP2/ERF, bHLH, MYB, WRKY, and bZIP, are linked with responses to biotic and abiotic stresses (Karanja et al., 2019; Lin et al., 2013; Liu & Liu, 2020; Pandey et al., 2017; Zhao et al., 2011). In line with other published findings, we found many TFs of the MYB family (127 genes), along with those of the AP2/ERF family (81 genes), BHLH family (55 genes), WRKY family (37 genes), and bZIP family (35 genes) (Table S13). The expression patterns of these five TF families were depicted in Figs. 5A–5E and conveyed in Table S14, which indicated that most of the unigenes in these five TF families (82/127 for MYB, 70/81 for AP2/ERF, 38/55 for BHLH, 31/37 for WRKY, and 29/37 for bZIP) were upregulated in the MD-treated samples at either 2 h or 6 h vis-à-vis the control samples. These screened TFs offered a pool of candidates that can be harnessed to potentially enhance plant resistance to MD.

Mining of candidate transcription factors related to monoterpenoid biosynthesis

Recently, the MYB-, bHLH-, ERF-, WRKY-, and bZIP-type TFs were shown capable of regulating the expression of plant monoterpene synthase genes (Sarker, Adal & Mahmoud, 2019). To select candidate TFs able to regulate the expression of monoterpene synthase in BCC, the abovementioned 335 DEGs belonging to these five TF families were submitted to a gene co-expression network analysis with the three selected candidate monoterpene synthase genes (CcTPS1, CcTPS3, CcTPS4). The results showed that three homologues of the BHLH family, six homologues of the bZIP family, five homologues of the WRKY family, 12 homologues of the ERF family, and 16 homologues of the MYB family all had strong, positive correlation coefficient values (R²-values > 0.8) with CcTPS1 and CcTPS4. In stark contrast, no TFs were found positively related to CcPTS3. Based on these results, interaction networks between the two monoterpene synthase and 42 transcription factors were organized (Fig. 6A, Table S15). Building on this, we carried out a more in-depth analysis of the sequences of these candidate TFs. First, by using the ORF finder, 37 of the TFs corresponding to 5 WRKY genes, 15 MYB genes, 10 ERF/AP2 genes, 5 bZIP genes, and 2 BHLH genes were discovered to contain a complete ORF. Next, the deduced amino acid sequences of each of these 37 TFs were submitted to the PfamScan, to determine whether a given sequence contained the conserved domains of the corresponding gene family. As a result, the 5 WRKY all contained complete WRKY domains (PF03106.15) and likewise, the 15 MYB genes all harbored a MYB DNA-binding domain with a Pfam ID of PF00249.31. The 10 ERF/AP2 genes, five bZIP genes, and two BHLH genes were confirmed to belong to their corresponding gene families since their respective AP2 (PF00847.20), bZIP_1 (PF00170.21), and HLH (PF00010.26) domains were identified. Taken together, these 37 TFs may play important roles in regulating monoterpene synthase expression to affect monoterpene biosynthesis in BCC. Therefore, this promising pool of candidate TFs is the focus of our future research.

qRT-PCR validation of DEGs from the RNA-Seq analysis

To validate the reliability of the RNA-Seq data, we chose 15 genes (nine monoterpene biosynthesis related genes, six transcription factors) to conduct the qRT-PCR analysis. These results basically agreed with the above transcriptome results; hence, our original analysis of the C. camphora transcriptome is robust (Fig. 6B).