RNA-, miRNA-, and degradome-sequencing datasets
For RNA-seq, 66.69 Gb clean reads were obtained from two groups through three independent biorepeats; quality control assessment showed Q20 values of 96.66% (Table 1). Then, 74,603 unigenes with an N50 of 1,464 bp were assembled, and the genes were annotated using the NR, Swiss-prot, COG, GO, and KEGG databases (Additional file 1, Additional file 2, Table 2). Here, 11,956 unigenes were differentially expressed between the callus cells and tissues (Additional file 1d, Additional file 2, Table 2). Most differentially expressed genes (DEGs) were mainly located at the plasma membrane and the integral module of membranes (Additional file 3a, Additional file 4) and significantly involved in “plant–pathogen interaction” and “plant hormone signal transduction” (Figure 1a, Figure 1b, Additional file 5).
Through miRNA-seq, 493 miRNAs were detected, 161 of which were newly identified in T. media (Additional file 6). A total of 95 miRNAs, including 35 novel miRNAs, were considered to be differentially expressed between the two groups with a p-value (Student T-test) of less than 0.05 (Figure 1a, Additional file 7, Additional file 8).
Degradome sequencing revealed that 1,829 unigenes were degraded by 347 miRNAs, leading to 2,432 degradation targets; of these, 323 unigenes were degraded by more than one miRNA (Additional file 9). Among the 347 identified miRNAs, cme-MIR166e-p5_2ss9CT19GC, mtr-MIR171c-p5_2ss1TC17GC, and ath-miR5021_R-1_1ss1TA degraded the most targets, specifically, 163, 93, and 84 unigenes, respectively. Among the degraded targets, two SPL-like TFs were degraded by the largest number of miRNAs, up to 11 (Additional file 9). SPLs (SQUAMOSA promoter-binding protein-like) form a plant-specific transcription factor family and participate in key activities; IPA1, for example, participates in the formation of plant architecture [16]. SPLs have recently been found to be tightly regulated by miRNAs in many plants, which indicates that they are crucial targets of miRNA and important nodes in the regulatory networks of plants [16-19].
Degraded differentially expressed (DE) targets were detected in approximately all DE pathways and significantly involved in 19 pathways; this result suggests that miRNAs are regulators involved in the transcriptional reprogramming of callus cells (Figure 2a, Additional file 10).
Callus cells are highly active in taxol biosynthesis but indirectly regulated by targeting miRNAs
The taxol content of newly induced callus cells was 1.34 mg/gDW (dry weight), which is 2.32 times higher than that in their parent tissues, and most biosynthesis genes were upregulated (Figure 3). Additionally, the amount of 10-deacetylbaccatin III, an intermediate precursor of taxol, in callus cells was 1.02 mg/gDW; such content is 3.15-fold higher than the content of parent tissues. The contents of two other taxanes in callus cells, namely, 10-deacetyl taxol and baccatin III, were also high (Figure 3b).
Taxol biosynthesis genes were found to be active, and 7 out of 12 known taxol biosynthesis genes were significantly upregulated in callus cells (Figures 3a and 3c). The expression of the rate-limiting gene, 10-deacetylbaccatin III-10-O-acetyl transferase (DBAT), which was barely expressed in tissues, increased by 70.9 times in callus cells. The expression of Taxadiene synthase (TASY), which is involved in the first step of taxol synthesis, increased by over 13.96-fold in callus cells; this gene also showed a high expression level in parent tissues. Five other genes, namely, phenylpropanoyltransferase (BAPT), taxadiene 5-alpha hydroxylase (T5H), taxane 2-alpha-O-benzoyltransferase (DBBT), 5-alpha-taxadienol-10-beta-hydroxylase (T10H), and taxadienol acetyl transferase (TAT), were significantly upregulated in callus cells (Figure 3c). These results suggest that TASY, DBAT, BAPT, T5H, and DBBT are critically important for taxol biosynthesis.
Not all biosynthesis genes were upregulated in callus cells. Phenylalanine ammonia-lyase (PAM) and taxane 2-alpha hydroxylase (T2H) were downregulated; in particular, the former was barely detectable in callus cells (Figure 3c). The functions of these two genes requires further elucidation.
Seven homologue genes of T7H, DBBT, BAPT, TAT, T10H, T13H, PAM, and DBTNBT were targeted by 10 miRNAs, none of which were differentially expressed (Figure 3a). T7H, DBBT, and BAPT are targeted by several miRNAs [7, 20-22]. Here, for the first time, BAPT was found to be targeted by miRNAs; T10H was targeted by miR5248 and miR397a, whereas BAPT was targeted by gma-miR6300 and the Taxus-specific miRNA PC-5p-97202_13 (Figure 4c). Moreover, PC-5p-97202_13 was identified in Taxus spp. for the first time and found to degrade 34 targets (Figure 4c, Additional file 11).
Degradome sequencing revealed that no taxol biosynthesis genes were degraded by any miRNA; however, a homologue of T10H was degraded by mtr-miR5248_2ss6AT21AT (Additional file 11). While T10H-like was significantly upregulated by 4.08-fold in callus cells, the expression of mtr-miR5248_2ss6AT21AT did not differ, thus suggesting that taxol biosynthesis is not directly regulated by miRNAs in callus cells.
Taken together, these 10 miRNAs which targeted taxol biosynthesis genes degraded 226 genes that are comprehensively involved in various primary metabolic processes and common defense activities, including “plant-hormone signaling transduction,” and “plant–pathogen interaction pathways” (Figures 4a and 4b, Additional file 11).
Biosynthesis of most secondary metabolites is upregulated but barely regulated by miRNAs
Most genes involved in the biosynthesis of secondary metabolites, especially flavonoids, phenylpropanoids, lignin, and lignans, were remarkably active in callus cells (Figures 1c and 5d). In particular, the methylerythritol phosphate (MEP) (non- mavalonic acid (MVA)) pathway, which produces terpenoid precursors, was significantly upregulated. By contrast, the MVA pathway, which is another means to produce terpenoid precursors, was downregulated (Figures 3a and 5d). Previous reports have indicated that MEP is a highly effective and efficient means to produce terpenoid precursors and likely a positive factor for callus cells to produce additional taxanes [23].
Only 17 DE genes of secondary metabolite biosynthesis were degraded by 15 miRNAs (Table 3). However, only five genes, namely, PER25, GT4, CYP86A22, UGT85A24, and SNL6, were upregulated because the four miRNAs that could degrade them were repressed (Table 3). Indeed, 598 DEGs involved in the biosynthesis of secondary metabolites were targeted by oppositely expressed miRNAs. These results indicate that miRNAs are capable of directly regulating secondary metabolism but they do not preferentially target metabolites in T. media callus cells, thus suggesting that a more cost-effective regulatory system for secondary metabolism should be available.
Callus cells strengthens the metabolism of biological materials and weakens the common defense system mainly mediated by miRNAs
Among the 43 significant DE pathways found in callus cells, the pathways related to metabolism of biological substances were all strengthened (Figure 1c, Additional file 5a, Additional file 12). For example, the metabolism of most amino acids, such as tyrosine and phenylalanine, was significantly enhanced. In addition, the DEGs of “Glycolysis” and “Citrate cycle” were remarkably upregulated; both of these pathways are key processes for decomposing sugar for living energy (Figure 1c, Additional file 5a, Additional file 12). Interestingly, “Photosynthesis” (35 downregulated vs. 3 upregulated) and “Photosynthesis-antenna proteins” (13 downregulated) were extremely weakened (Figure 1c, Additional file 5a, Additional file 12). This finding indicates a transformation in living behavior from autotrophic tissues to heterotrophic cells.
Degradome sequencing confirmed that degraded DEGs are significantly enriched in the metabolism of several biological materials, such as “Starch and sucrose metabolism,” “Ascorbate and aldarate metabolism,”, and “Aminoacyl-tRNA biosynthesis” (Figure 2b). In addition, most of the degraded DEGs involved in these pathways were regulated by oppositely expressed miRNAs, thus suggesting that miRNAs effectively regulate these pathways (Figure 2c, Additional file 13).
Biotic and abiotic stresses are the main threats to living plants; to address these stresses, plants execute a number of response activities, such as “Plant–pathogen interaction,” “Plant hormone signal transduction,” “Phagosome,” and “Endocytosis” [24-30]. Callus cells showed more downregulated genes involved in abiotic and biotic stress responses compared with tissues but showed upregulated heat/cold/light responses (Figure 5a). For “Plant–pathogen interaction” and “Plant hormone signal transduction”, 70% (434) and 62% (209) DEGs were downregulated (Figures 5b and 5c, Additional file 5a). For example, jasmonic acid-amino synthetase 1 (JAR1) and nonexpresser of pathogenesis-related 1 (NPR1), which positively function in JA- and SA-signaling transduction, were significantly downregulated (Additional file 14, Additional file 15). These results suggest that callus cells are fragile in terms of stress response and disease resistance.
MiRNAs are the main factors regulating common defense activities and significantly enriched in “Plant–pathogen interaction” “Plant hormone signal transduction”, and “Endocytosis” [31]. Among the pathways observed, “Plant–pathogen interaction” and “Plant hormone signal transduction” were the most enriched pathways regulated by miRNAs; indeed, 77 and 59 DE degraded targets of pathogen interaction and plant hormone signal transduction were respectively detected in callus cells (Figure 2b, Additional file 10). In addition, degraded DE targets were detected in “Phagosome”, “Spliceosome” and “Base excision repair”, which are other pathways related to defense. Such findings indicate that miRNA is an important regulator of these primary metabolism pathways and common defense activities. In particular, pathways, such as “plant–pathogen interaction” and “plant hormone signal transduction”, are the prior regulation targets of miRNAs.
TFs are important targets of miRNAs
Among 1,830 degraded targets, 635 (34.7%) were TFs and various transcriptional regulators and degraded by 236 miRNAs, thus constituting 894 pairs of miRNA–TF modules (Additional file 9). Among the 635 degraded TFs, PHD, C3H, bHLH, WRKY, MYB, and NAC were mostly regulated by miRNAs (Figure 6a). A total of 156 miRNAs degraded more than one TF, and cme-MIR166e-p5_2ss9CT19GC degraded 71 TFs (Figure 6b, Additional file 16). A total of 426 miRNA-TF modules showed contrasting expression patterns, and 292 pairs (68.6%) were constituted by downregulated miRNA and upregulated TF, thus indicating that callus cells repress the expression of miRNAs to regulate bioactivities (Additional file 16). Previous reports also concluded that miRNAs repress the expression of most genes until these genes are needed [32].
Among the TFs, ERF, bHLH, NAC, and MYB had the highest number of degraded fragments because one TF were degraded by several miRNAs, leading to a higher number of degraded fragments than DEGs (Figure 6b). ERF, SBP, NZZ/SPL, and NF-YA were the most enriched targets of degradation, thus suggesting their importance in transcriptional reprograming from tissues to callus cells (Figure 6b).
The roles of these TFs were analyzed here on the basis of studies on several known Taxus TFs [33-36]. TcERF15, TcMYC2a, TcJAMYC1/2/4, and TcWRKY1/8/20/26/47/52 were highly correlated with taxol biosynthesis genes with Pearson coefficients greater than 0.9 (Figure 7b, Additional file 17). Interestingly, 4 homologues, TcMYC2a and TcJAMYC1/2/4, showed contrasting expression patterns in callus cells. These results are highly consistent with results from functional studies. TcMYC2a from Taxus chinensis and TcJAMYC1/2/4 from Taxus cuspidata function as positive and negative regulators in taxol biosynthesis, respectively [34].
Candidate TFs and miRNAs involved in the regulation of taxol biosynthesis
To determine candidate TF regulators, Pearson coexpression correlations were analyzed in combination with the previous dataset of long-term subcultured T. chinensis cells [7]. A total of 451 TFs were similarly or oppositely expressed with taxol biosynthesis genes (Figure 7a), and 346 of 451 TFs were highly correlated with taxol biosynthesis genes with coefficient values greater than 0.9 (Figure 7b, Additional file 17). T13H, T2H, and PAM were correlated with 161, 145, and 141 TFs, respectively. Although DBAT and T7H were correlated with 107 TFs, they were simultaneously correlated with only 38 TFs (Figure 7b, Additional file 17). TF families, NAC, WRKY, bHLH, and ERF were the most coexpressed candidate regulators.
Among the 346 miRNAs mediating TF degradation, the expression patterns of only 28 miRNAs were similar to those of taxol biosynthesis genes. Subsequent coexpression analysis confirmed that 21 miRNAs were closely related to taxol biosynthesis with high coefficient values (>0.9 or <−0.9). Nine taxol biosynthesis genes, namely, T5H, TAT, T10H, DBBT, DBTNBT, T7H, T13H, PAM, and T2H, were related to miRNAs. PAM and T13H were mainly related to 10 miRNAs; this relationship indicates that PAM and T13H are crucial regulatory targets. Hbr-MIR6173-p5_1ss9TG was coexpressed with 4 taxol biosynthesis genes. In addition, a novel miRNA first identified in Taxus, PC-5p-97202_13, simultaneously targeted DBBT and DBTNBT; this characteristic suggests its important regulatory roles in taxol biosynthesis (Figures 4c, 4d, and 7c). Although these coexpressed miRNAs do not degrade taxol biosynthesis genes, they may regulate taxol biosynthesis through other means.