Transcriptomic analyses reveal putative genes associated with bibenzyl biosynthesis in the traditional Chinese herb Dendrobium officinale

Background: Dendrobium plants are well known for their uses in traditional Chinese herbal medicine. Bibenzyl compounds are the main active compounds in Dendrobium officinale. However, the physiological and molecular basis for the biosynthesis of bibenzyl compounds in Dendrobium plants remains underexplored. Results: In this study, the accumulation of erianin and gigantol were studied as representative compounds of bibenzyl. Their presence in plant tissues were investigated. Our results show that root tissues contained the highest content of bibenzyl (erianin and gigantol). Based on the pre-experimental result that exogenous application of Methyl-Jasmonate promotes the biosynthesis of bibenzyl compounds in D. officinale root tissues, comparative transcriptomic analyses were conducted between the bibenzyl-accumulated root tissues and a control. In total, we identified 1,342 differentially expressed genes (DEGs) with 912 up-regulated and 430 down-regulated genes. Most of the identified DEGs are functionally involved in the JA signaling pathway and the biosynthesis of secondary metabolites. In particular, we identified 11 enzymatic genes functionally involved in bibenzyl biosynthesis. Conclusions: Our study provide insights on the identification of putative genes associated with bibenzyl biosynthesis and accumulation in Dendrobium plants, and also paves the way for future research on dissecting the physiological and molecular mechanisms of bibenzyl synthesis in plants as well as on how to best utilize genetic engineering and molecular modification techniques to genetically improve Dendrobium varieties by increasing the content of bibenzyl for drug production and industrialization.


Background
Plants produce vital secondary metabolites for growth and development and also in response to environmental stresses. These secondary metabolites (such as alkaloids, phenolics, flavonoid, and terpenoids) often accumulate within a specific group of plants or given tissues, which play crucial roles in helping plants in defense against various biotic and abiotic stresses [1][2][3]. In particular, these secondary metabolites provide essential resources for new drug innovations, insecticides, and flavors [4][5][6].
The Dendrobium plants (an herb in the family Orchidaceae) are known as traditional Chinese medicinal herbs (referred to as shihu in Mandarin). Dendrobium is widely distributed across Asia and the Pacific Islands [7]. Previous studies documented the health benefits (antipyretic, ophthalmic and regulative of immune system) of Dendrobium plants and their contribution to Chinese medicines [8]. In particular, owing to its wealth of active compounds with antitumor and antioxidants functions, D. officinale has received tremendous interest in Asian countries. In recent years it has been cultivated in many regions to provide feedstock for making cosmetic and medicinal products.
The active medicinal ingredients in D. officinale include: polysaccharides, alkaloids, phenols, terpenes, flavonoids and bibenzyl [9][10]. In Dendrobium plants, the content of sesquiterpene alkaloid is the main measure of its quality and medicinal efficacy [11][12]. Notably, previous studies have found that bibenzyl compounds (belonging to sesquiterpene alkaloids) might be the only bioactive ingredients in D.
The biosynthesis pathway of secondary metabolites involves various physiological factors and regulatory modifications in different plants that have developed their strategies in response to environmental changes or stresses. Usually, the accumulation of secondary metabolites (such as bibenzyl compounds) is low among tissues [5]. However, market demands require a higher content of bibenzyl compounds in D. officinale tissues to meet the threshold for drug-making. There is, therefore, an immediate need to dissect the physiological and molecular mechanisms underlying how bibenzyl compounds are biosynthesized in D. officinale tissues. The identification of rate-limiting enzymes or regulatory factors, which are responsible for the biosynthesis of bibenzyl compounds, is thus also area in need of further exploration. It is exceedingly helpful to use genetic engineering and 6 molecular modification techniques to create improved varieties for commercial purposes. Investigating the physiological and molecular basis of the accumulation of bibenzyl compounds is an essential prerequisite to understanding molecular and genetic factors that regulate the biosynthesis of bibenzyl compounds in D. officinale tissues. Moreover, studies have found that jasmonate (JA), a plant-specific signaling molecule, is widely involved in the biosynthesis of diverse secondary metabolites.
The exogenous application of methyl-jasmonate (MeJA) often results in the strong activation of secondary metabolites biosynthesis and has frequently been applied to induce the biosynthesis of secondary metabolites in plants [27][28][29].
With the rapid advancement in RNA-seq technology, transcriptomic data offers both a great opportunity and powerful tool for the discovery of crucial rate-limiting enzymes or regulators, which control the production of some secondary metabolites in plants under different conditions [5,[30][31]. For instance, based on transcriptomic analysis, several putative rate-limiting enzyme genes responsible for the biosynthesis of terpenoids in Eugenia uniflora [32], lignin in Apium graveolens [33], and flavonoid in Phyllanthus emblica, Dracaena cambodiana, and Solanum viarum [34][35][36] have been identified. For Dendrobium plants, several studies have reported on flavonoid biosynthetic pathway analysis and the gene mining of key enzymes [37] as well as alkaloid biosynthetic pathway analysis and the identification of key enzyme genes [12,[38][39]. However, the bibenzyl biosynthetic pathway and the potential genes responsible for regulating bibenzyl biosynthesis in Dendrobium plants remain underexplored.
In this study, we investigated the accumulation of bibenzyl in various tissues of D.
officinale including in the leaf, root, basal stem and upper stem. We conducted comparative transcriptomic analyses to unravel the putative genes involved in the 7 biosynthesis of bibenzyl in D. officinale. To our knowledge, this study is the first report that focus on the identification of putative genes associated with the bibenzyls biosynthesis in Dendrobium plants. This study will aid our understanding of unique genes involved in the synthesis of bibenzyl in D. officinale as well as provide new insights for future research into the molecular mechanisms of the genes involved in bibenzyl biosynthesis.

Investigation of Bibenzyl Accumulation Among Tissues
According to previous studies, bibenzyl compounds (mainly including erianin and gigantol) are the main bioactive ingredients in D. officinale for medicine production [15,18,40]. To investigate tissues with the highest accumulation of bibenzyl contents, we initially assessed the content of bibenzyl compounds in four different tissues (leaf, root, basal stem and upper stem tissues) of a three-year-old D.
officinale plant. As shown in Table 1, the root tissues had the highest content of erianin and gigantol, which were 2.63±0.69 and 37.01±2.16 µg/g, respectively, followed by the basal stem (0.61±0.01 and 22.67±0.15 µg/g). Within the upper stem, erianin was not detectable while the content of gigantol was lower compared to the root and the basal stem tissues. However, we could detect neither erianin nor gigantol in the leaf tissues. These results indicate that bibenzyl biosynthesis and accumulation mostly occur in the root tissues. officinale with a methyl-jasmonate (MeJA) solution at different concentrations (0.2 8 mM, 0.5 mM, and 1.5 mM). We observed that the contents of bibenzyls (erianin and gigantol) in root tissues as showed in Fig 1a and b were significantly higher (13.01fold and 8.43-fold increase, respectively) after being treated for 36 hours. However, the 0.5 mM concentration seemed to be the optimal accumulation for bibenzyl induction in root tissue. These results clearly show that the accumulation of bibenzyl was significantly induced after 24 hours by exogenous MeJA in root tissues.

Transcriptome Sequencing Datasets
To dissect the bibenzyl biosynthetic pathway and putative essential rate-limiting genes associated with bibenzyl biosynthesis in D. officinale, we carried out transcriptomic analysis using high-throughput RNA-seq technology.

Identification of Differentially Expressed Genes
To identify putative genes associated with bibenzyl biosynthesis in D. officinale, we analyzed the differentially expressed genes (DEGs) between the two sequenced datasets. Parameters of the False Discovery Rate (FDR) were set at < 0.05 and |log2 Fold change| >2 for identifying DEGs. In total, 1,324 DEGs, consisting of 912 upregulated and 430 down-regulated DEGs, were identified compared to the CK (

Identification of Candidate Genes Involved in Bibenzyl Biosynthesis
Our principal objective was to identify candidate genes involved in bibenzyl biosynthesis in this study. Bibenzyl, an alkaloid belonging to the group of sesquiterpene [14,41], is a downstream product of mevalonate (MVA) and methylerythritol 4-phosphate (MEP) pathway in plants [12,42]. In our dataset, most of the critical enzymes involved in these pathways such as hydroxymethylglutaryl- anticancer and immunomodulatory activities [18,19,45,16,46,47]. Several studies have been conducted on the identification of putative genes involved in polysaccharide biosynthesis [9,41,48]. However, little is known regarding the physiological and molecular bases of bibenzyl biosynthesis in planta. To our knowledge, this study is the first investigation into the biosynthesis of bibenzyl at the physiological and molecular levels. Bibenzyl compounds are composed of various formulas, with a pair of benzyl radicals. Of them, the erianin and gigantol are structurally similar, and they are represented as bibenzyl compounds [46,47,40,49]. As a result, in this study, the contents of erianin and gigantol were measured as representatives of bibenzyl compounds.
A recent study detected the total alkaloid content of D. officinale in the leaf and found a significant increase after exogenous MeJA treatment [50]; however, the contents of erianin and gigantol were not detectable in leaf tissues in our study. In our study, we found that bibenzyl compounds mainly accumulate in the roots and the basal parts of the stem tissues. This means that the roots of Dendrobium plants may be more critical in the extraction of bioactive ingredients of antioxidant and anticancer compounds than other tissues. Many studies have found that the 13 biosynthesis of secondary metabolites (such as terpenoids, phenylpropanoids, flavonoids, and alkaloids) could be induced by hormone JA signaling [50][51][52]. As expected, the accumulation of bibenzyl was induced by the exogenous hormone MeJA in this study. The content of bibenzyl significantly increased between 24 to 36 hours after MeJA treatment, suggesting the rapid accumulation of bibenzyl during this timeframe. We reasonably assumed that most genes involved in bibenzyl biosynthesis were actively expressed. Thus, we compared global gene expressions between the root tissues at the highest accumulation of bibenzyl period and that of the control (untreated).
Although several generated transcriptomic data have been available for Dendrobium [9,48,53,50], in our study, based on the root tissues 23,131 unigenes were annotated, and this number was smaller when compared to previous data, but it is comparable to the identified unigenes using transcriptomic data of root tissues generated by Chen et al., (2017). In total, 1,324 DEGs, including 912 up-regulated and 430 down-regulated DEGs, were identified through the comparison of expression changes between our two libraries. Compared to a recent investigation of transcriptomic analyses using exogenous MeJA treatment in D. officinale leaves [50], we identified fewer DEGs. It is likely that the fewer identified unigenes and DEGs in our study resulted from the difference of tested tissues between these studies (only root tissues were investigated in this study, whereas other studies investigated either a leaf, stem or mixed tissues). Most of the identified DEGs were mainly enriched in the JA signaling pathway and biosynthesis of secondary metabolites in this study. The induced DEGs were significantly enriched in the GO terms related to the diverse processes of secondary metabolites (including sesquiterpene and flavonoid biosynthesis) and in response to JA induction in the signaling pathway. As expected, many DEGs involved in the JA signal pathway were identified, such as PLD, DAD1, LOX, and AOS. In particular, JAZ and MYC are wellknown responses to exogenous MeJA treatment in plants [54][55][56][57][58][59].
In D. officinale, several transcriptomic analyses have been performed, examining topics such as prediction of putative metabolic pathways [53], identification of functional genes in various metabolic pathways [9,48,50,60] and uncovering environmental responses [61]. In this study, our main objective was to identify potential candidate genes involved in the biosynthesis of bibenzyl (a group of sesquiterpenes) in D. officinale. Generally, sesquiterpenes are derived from farnesyl diphospate (FPP), which is provided by the MVA and MEP pathways in the initial stage of sesquiterpenes biosynthesis in plants [42]. Recently, Chen et al. (2019) have identified several genes involved in FPP biosynthesis. These identified genes usually function at the initial stages of sesquiterpenes biosynthesis, and most of these identified genes such as HMGS, MK, PMK, DXS, DXR, CMK, MDS, and HDR appeared in our data, suggesting that these genes participated in regulating the initial biosynthesis of sesquiterpenes in D. officinale and also participated in the early biosynthesis of bibenzyl. In particular, we identified 11 enzyme genes (including two PALs, two C4Hs, two 4CLs, two P450s, and three BBSs) functionally involved in bibenzyl biosynthetic processes. Although it is well known that the enzyme genes CYP450, PAL, C4H and 4CLs with multi-members are functionally involved in various metabolic processes in plants [62][63][64][65][66], this study identified that these 11 enzyme genes may be especially active in regulating bibenzyl biosynthesis in D. officinale. It is possible that these genes (or some of them) may be ratelimiting for bibenzyl biosynthesis, and the expression levels of these genes may

Illumina Sequencing Data Analysis
To improve the sequence quality, reads with poly-N and low-quality fragments were removed to obtain high-quality reads. The clean reads were mapped onto the reference genome of the species using HISAT2 [68]. The fragments per kilobase per million map reads (FPKM) was calculated by Cufflinks to estimate the level of gene expression. DEGSeq was used to detect differently expressed genes (DEGs) between CK and MeJA treated samples. Genes were identified significantly differently expressed with a False Discovery Rate (FDR) of <0.05 and |log2 (fold change) | >2.
Functional GO and KEGG enrichment analyses of DEGs were respectively performed. P-value and FDR correction were used to determine the significance of enriched pathways. Multi-Experiment Viewer (MeV) (version 4.9.0) was used to generate heatmaps of DEGs.

Analysis of Differentially Expressed Genes (DEGs)
Differential expression analysis was designed to identify genes with differential expressions between different samples, according to the negative binomial distribution test in DESeq [69] software (http://bioconductor.org/packages/release/bioc/html/DESeq.html). The amount of differential expression of the genes was calculated while the NB (negative binomial distribution test) was used to test the difference in the number of reads. The gene expression was estimated by using the baseMean value. To compare and analyze whether the same gene in two samples is differentially expressed, two criteria were selected: one is FoldChange, which is the fold change of the same gene expression level in the two samples; the other is the p-value or FDR (adjusted p-value). The default screening difference condition was p-value<0.05 and the differential multiple was greater than 2. To identify genes related to bibenzyl biosynthesis, unigenes were searched against the flavonoid biosynthetic pathway.
Validation of genes related to bibenzyl biosynthesis by qRT-PCR  Table 1   Table 1 Comparison of bibenzyls contents (erianin and gigantol) in different tissues of D.
officinale (Mean ± Standard error). ND denotes "not detectable".  Identification of differentially expressed genes involved in the JA signal pathway (a) and TFs