CIR treatment changed the root growth parameters, lignification and saponin content
To investigate the influence of the CIR treatments on the market value of B. chinense DC. plants. The inflorescences in treated plants (defined as T) were removed as soon as inflorescence buds emerged since the day of the first inflorescence bud appeared. The treatment was continuously conducted for 75 days. Then roots were sampled for parameters measurement and chemical components analyses. Root biomass was an important index of market value of cultivar of B. chinense DC.. The degree of lignification was also important for its market value since the dried roots of B. chinense DC. were cut into decoction pieces in traditional processing and serious lignification means higher mechanical strength and processing cost. Thus, we examined the effects of CIR on the root biomass and growth parameters, the total xylem cross-sectional area (CSA), total lignin contents, as well as the contents of saponin a, c, d which are the main medicinal ingredients of B. chinense DC..
The results showed that the biomass and growth parameters including root fresh weight, root dry weight, the taproot length, root head diameter and lateral root number of B. chinense DC. were all significantly higher in treated plants than in control plants. Especially for the fresh weight, dry weight, and lateral root number which showed 1.71-, 1.76-, and 0.87-fold increase, respectively (Fig. 1a-e). For the lignification investigation, although there was little difference of xylem CSA between treated and control plants, the lignin contents in roots of treated plants were found to be a 33.73 percent decrease (Fig. 1f and g). Thus the CIR treatment reduced the level of lignifications. The microstructures of root were also observed with light microscope. Like control plants (Additional file 1: Fig. S1a-d), intact anatomical structures in cross-sections of roots could be also observed in treated plants (Additional file 1: Fig. S1e-h). Furthermore, HPLC analysis showed that the content of saponin a was significantly higher in treated plants with a 32.01 percent increase, while the content of saponin c and d were not significantly changed (Fig. 1h-j).
Numbers of differentially expressed genes increased along with the duration of CIR
The differentially expressed genes (DEGs) were identified by pair-wise comparisons of these libraries, T1 vs FB, T2 vs IG, T3 vs R, T4 vs H. A total of 172 DEGs were identified from transcriptome of T1 vs FB, with 80 up-regulated and 92 down-regulated genes, while 243 (151 up-regulated and 92 down-regulated), 1974 (1385 up-regulated and 589 down-regulated), 3024 (1788 up-regulated and 1236 down-regulated) genes from compares of T2 vs IG, T3 vs R and T4 vs H, respectively (Fig. 2a). Obviously, the DEGs were obviously increased in the later stages.
In order to identify common and unique DEGs in response to CIR at different time points, a Venn diagram was plotted (Fig. 2b). Among these DEGs, none common DEG was identified in four groups. Little common DEG genes were found for the first three time points. Four genes were shared between the “T1 vs FB” and “T2 vs IG” comparisons and 17 genes were shared between the “T2 vs IG” and “T3 vs R” comparisons. However, the “T3 vs R” and “T4 vs H” shared 738 DEGs.
To validate the accuracy of RNA-seq data, 16 DEGs were selected stochastically from transcriptome dataset and qRT-PCR analysis was carried out (Additional file 2: Table S1). Compared with the sequencing results, the qRT-PCR results showed the similar trends (Additional file 3: Fig. S2). Hence, the reliability of the expression level of RNA-seq data was confirmed.
Thus, the operation of CIR triggered or depressed the expression of various genes, and the DEG numbers of four groups displayed a markedly increasing tendency along with the duration of CIR operation.
GO and KEGG pathway enrichment analysis of DEGs
To further investigate the main biological functions of these DEGs from four groups, the DEGs were subjected to the GO and KEGG database. In general, 4813 DEGs from four groups were GO-annotated into 44 terms. 172 DEGs in “T1 vs FB” were GO-annotated into 31 terms. In “T2 vs IG”, 243 DEGs were GO-annotated into 30 terms. In “T3 vs R”, 1974 DEGs were GO-annotated into 43 terms, and 3024 DEGs in “T4 vs H” were GO-annotated into 41 terms (Fig. 3 and Additional file 4: Table S2). Although the DEG numbers were diverse, “cellular process”, “metabolic process”, and “single-organism process” were the most enriched terms in biological process (BP), and “catalytic activity” and “binding” were the most enriched terms in molecular function (MF) for all the four time points.
There were five, four, six and five KEGG categories based on pathway hierarchy1 (Additional file 5: Table S3). Among the pathway subgroups based on pathway hierarchy 2, the pathway of Carbohydrate metabolism contained the most genes, followed by Lipid metabolism (Fig. 4 and Additional file 5: Table S3). All of the KEGG-annotated DEGs were assigned to 23, 31, 86 and 97 KEGG pathways, respectively (Additional file 5: Table S3). For DEGs of “T1 vs FB” and “T2 vs IG”, “sesquiterpenoid and triterpenoid biosynthesis” and “Protein processing in endoplasmic reticulum” was the top enriched pathway, respectively. The top four enriched pathways in both “T3 vs R” and “T4 vs H” were “starch and sucrose metabolism”, “pentose and glucuronate interconversions”, “plant hormone signal transduction” and “phenylpropanoid biosynthesis” (Additional file 5: Table S3) .
Major DEGs profiling of sugar metabolic pathways
Based on the KEGG pathway enrichment analyses, DEGs involved in the carbohydrate metabolism were the most abundant, so we analyzed the “starch and sucrose metabolism”, “galactose metabolism” and “pentose and glucuronate interconversions” pathways.
38 and 48 DEGs from “T3 vs R” and “T4 vs H” were enriched in “starch and sucrose metabolism” pathway, respectively. Among them, 36 DEGs were up-regulated and two were down-regulated in “T3 vs R”, while 38 DEGs were up-regulated and 10 were down-regulated in “T4 vs H” (Fig. 5a and b and Additional file 6: Table S4). Four and five DEGs from “T3 vs R” and “T4 vs H” were enriched in “galactose metabolism” pathway, respectively (Additional file 6: Table S4 and Additional file 7: Fig. S3a and c). Both “T3 vs R” and “T4 vs H” had 31 DEGs enriched in “pentose and glucuronate interconversions” pathway (Additional file 6: Table S4 and Additional file 7: Fig. S3b and d). Most of these DEGs were up-regulated. In “T3 vs R”, 29 genes were up-regulated, while in “T4 vs H”, the number was 28. The annotation showed that some DEGs including xylosidase gene, xylan synthase gene, UDP-glucuronate decarboxylase gene (UXS), α-galactosidase gene, β-fructofuranosidase gene, endoglucanase gene and glucosidase gene were related to the xylose, galactose and glucose biosynthesis and metabolism (Additional file 6: Table S4 ), indicating the activation of the production of xylose, galactose and glucose related metabolites. These sugars are important in the biosynthesis of saponin [2]. Some DEGs including UDP-glucuronate 4-epimerase gene, alpha-1,4-galacturonosyltransferase gene (GAUT), pectinesterase gene (PE) , pectinesterase inhibitor gene (PEI), exopolygalacturonase gene and pectate lyase gene (PL) were related to the pectin metabolism (Additional file 7: Fig. S3), and pectin is a major component of primary plant cell wall. Notably, CIR treatments increased the abundance of PL, endoglucanase, PE, and PEI transcripts in these two groups. One PEI gene (Cluster-10318.20452) had a 4.66-fold change in “T4 vs H”, followed one exopolygalacturonase gene (Cluster-10318.10716) with a 4.36-fold change in “T3 vs R” (Fig. 5 and Additional file 6:Table S4).
Expression of DEGs involved in phenylpropanoid and terpenoid biosynthesis pathways
Lignin was generated through phenylpropanoid biosynthesis pathway which had been well elucidated [17]. Saponin was biosynthesized through terpenoid biosynthesis pathway [18]. In our study, the DEGs from phenylpropanoid and terpenoid biosynthesis pathways were detected. Three, 19 genes and 34 genes from “T2 vs IG”, “T3 vs R”, and “T4 vs H” were enriched in phenylpropanoid biosynthesis pathway, respectively. Among the DEGs, some upstream genes were up-regulated. For instance, one phenylalanine ammonia lyase gene (PAL) was up-regulated in “T4 vs H”, and three hydroxycinnamoyl-CoA:shikimate/quinate hydroxycinnamoyltransferase genes (HCT) were up-regulated in “T4 vs H”. One and six beta-glucosidase genes were up-regulated by CIR in “T2 vs IG” and “T3 vs R”, respectively. One caffeic acid O-methyltransferase gene (COMT; Cluster-10318.50166) was up-regulated in “T3 vs R”, and another gene in this family (Cluster-10318.33017) was up-regulated in “T4 vs H”. Interestingly, we observed the serious reduction of the expression level of seven genes encoding cinnamyl alcohol dehydrogenase (CAD) in “T3 vs R” or “T4 vs H” group. CAD is a key enzyme in monolignin biosynthesis, and down-regulation of CAD gene generally leads to reduction of lignin content [19]. The last step in lignin biosynthesis pathway is catalyzed by peroxidase (POD) which transforming cinnamyl alcohols to lignin. A number of POD genes were identified in this study. Both up-regulated genes and down-regulated genes were exhibited in this gene family. (Fig. 6 and Additional file 8: Table S5).
Likewise, CIR modulated the expression of most terpenoid biosynthetic genes in four groups, and most of these genes were up-regulated in the last three groups. Particularly, these DEGs including one acetyl-CoA acetyltransferase gene (AACT), one hydroxymethylglutaryl-CoA synthase gene (HMGS), four hydroxymethylglutaryl-CoA reductase genes (HMGR), two mevalonate diphsphate decarboxylase genes (MVD), two 1-deoxy-dxylulose-5-phosphate synthase genes (DXS), one 4-diphosphocytidl-2C-methyl-D-erythritol synthase gene (MCT), two diphosphomevalonate decarboxylase genes (DPMD), six farnesyl diphosphate synthase genes (FPS), nine germacrene D synthase genes (GERD), one sesquiterpene synthase gene (SQS), one squalene synthase gene (SS), one squalene epoxidase gene (SE), two lupeol synthase genes (LUP) and one β-amyrin synthase gene. All the above genes except one FDPS (Cluster-10318.59484) were up-regulated by CIR in “T2 vs IG”, “T3 vs R” and “T4 vs H” (Fig. 7 and Additional file 9: Table S6).
Transcription factors involved in response to CIR
In total, 295 DEGs encoding TFs were obtained in transcriptome data. Specifically, there were five DEGs (three up-regulated and two down-regulated) in “T1 vs FB”, seven DEGs (three up-regulated and four down-regulated) in “T2 vs IG”, 103 DEGs (55 up-regulated and 48 down-regulated) in “T3 vs R” and 180 DEGs (95 up-regulated and 85 down-regulated) in “T4 vs H”. These genes encoding TFs included members of ERF, MYB-related, bHLH, NAC, MYB, AUX/IAA, bZIP families etc. Among all the differentially expressed TFs in “T3 vs R”, MYB-related was the predominant family with eight DEGs, including one up-regulated gene (1/8), followed by the AUX/IAA (7/7) and MYB (4/6) families. In “T4 vs H”, ERF was the predominant family with 20 DEGs, including nine up-regulated genes (9/20), followed by the bHLH (13/14) and MYB (7/9) families (Additional file 10: Table S7).
Hormone signal transduction-related DEGs involved in response to CIR
The expression of genes involved in hormone signal transduction pathways was altered following the CIR treatment. According to the KEGG analysis of DEGs, 30 and 46 DEGs in this pathway were identified from “T3 vs R” and “T4 vs H”, respectively (Fig. 8). Among different hormone signal transduction pathways, auxin signal transduction pathway contained the largest number of identified DEGs for both the “T3 vs R” and “T4 vs H” groups, followed by the cytokinine and brassinosteroid signal transduction pathways.
Auxin, inducing three kinds of early auxin response genes: AUXIN/INDOLE-3-ACETIC ACID (AUX/IAA), GRETCHEN HAGEN3 (GH3), and SMALL AUXIN UP RNA (SAUR) which had been identified in diverse plant species [20], was supposed to be the primary hormones that regulate the lateral roots development and elongation [21]. We observed significant changes of 10 and 11 genes involved in auxin signal transduction pathways of “T3 vs R” and “T4 vs H”, respectively, including one gene encoding transport inhibitor response protein (TIR), twelve auxin-responsive protein genes (nine AUX/IAA genes and three SAUR genes), and five auxin response factor genes (ARF). TIR functioned as auxin receptor [22]. Most of AUX/IAA genes were up-regulated, while all ARF genes were down-regulated in “T3 vs R” and “T4 vs H”. In addition, the three genes encoding SAUR included two up-regulated genes and one down-regulated gene, and one gene (Cluster-10318.15908) had a 3.39-fold change in “T4 vs H” (Additional file 11: Table S8).
Among the cytokinine signal transduction-related genes, 100% of DEGs (6/6) were up-regulated in “T3 vs R”, and 66.7% of DEGs (6/9) were up-regulated in “T4 vs H”. Three and four genes encoding histidine kinases (HKs), the cytokinine-binding receptors, were up-regulated in “T3 vs R” and “T4 vs H”, respectively. And one HK gene was down-regulated in “T4 vs H”. A multistep HK→histidine-containing phosphotransfer protein (HP)→response regulator (RR) phosphorelay pathway was involved in response to cytokinine [23]. One gene encoding HP functioned as phosphor-histidine intermediates between the HK and RR was up-regulated in “T4 vs H”. Three RR genes were up-regulated in “T3 vs R”, and while one RR gene was up-regulated and two RR genes were down-regulated in “T4 vs H” (Additional file 11: Table S8).
In brassinosteroid (BR) signal transduction pathway, 80% of DEGs (4/5) were up-regulated in “T3 vs R”, and 77.8% of DEGs (7/9) were up-regulated in “T4 vs H”. Brassinosteroid insensitive (BRI) functioned as the brassinosteroid receptor and it is a member of the leucine-rich repeat receptor-like kinase family. One gene encoding BRI kinase inhibitor which is the inhibitor of brassinosteroid receptor (BKI) was down-regulated in “T3 vs R” and “T4 vs H”. Moreover, two brassinosteroid-signaling kinase genes (BSK, Cluster-10318.45934 and Cluster-10318.84972), the important positive regulator of BR signal, were up-regulated in “T3 vs R” and “T4 vs H”. In general, BR triggered active cell division associated with increased transcripts of cell cycle-related genes, especially that of cyclin D3 genes [24]. Notably, seven cyclin D3 genes were identified and all genes were up-regulated by CIR treatment (Additional file 11: Table S8).
In addition, DEGs enriched in GA (“T4 vs H”), ABA (“T3 vs R” and “T4 vs H”), Ethylene (“T4 vs H”), JA (“T3 vs R” and “T4 vs H”), SA (“T3 vs R” and “T4 vs H”) signal pathways were alsoidentified, and displayed a trend of down-regulation except the JA signal transduction-related genes (Fig. 8).