Transcriptional reprogramming strategies and miRNA-mediated regulation networks of Taxus media induced into callus cells from tissues

Background Taxus cells are a potential sustainable and environment-friendly source of taxol, but they have low survival ratios and slow grow rates. Despite these limitations, Taxus callus cells induced through 6 months of culture contain more taxol than their parent tissues. In this work, we utilized 6-month-old Taxus media calli to investigate their regulatory mechanisms of taxol biosynthesis by applying multiomics technologies. Our results provide insights into the adaptation strategies of T. media by transcriptional reprogramming when induced into calli from parent tissues. Results Seven out of 12 known taxol, most of flavonoid and phenylpropanoid biosynthesis genes were significantly upregulated in callus cells relative to that in the parent tissue, thus indicating that secondary metabolism is significantly strengthened. The expression of genes involved in pathways metabolizing biological materials, such as amino acids and sugars, also dramatically increased because all nutrients are supplied from the medium. The expression level of 94.1% genes involved in photosynthesis significantly decreased. These results reveal that callus cells undergo transcriptional reprogramming and transition into heterotrophs. Interestingly, common defense and immune activities, such as “plant–pathogen interaction” and salicylic acid- and jasmonic acid-signaling transduction, were repressed in calli. Thus, it’s an intelligent adaption strategy to use secondary metabolites as a cost-effective defense system. MiRNA- and degradome-sequencing results showed the involvement of a precise regulatory network in the miRNA-mediated transcriptional reprogramming of calli. MiRNAs act as direct regulators to enhance the metabolism of biological substances and repress defense activities. Given that only 17 genes of secondary metabolite biosynthesis were effectively regulated, miRNAs are likely to play intermediate roles in the biosynthesis of secondary metabolites by regulating transcriptional factors (TFs), such as ERF, WRKY, and SPL. Conclusion Our results suggest that increasing the biosynthesis of taxol and other secondary metabolites is an active regulatory measure of calli to adapt to heterotrophic culture, and this alteration mainly involved direct and indirect miRNA-induced transcriptional reprogramming. These results expand our understanding of the relationships among the metabolism of biological substances, the biosynthesis of secondary metabolites, and defense systems. They also provide a series of candidate miRNAs and transcription factors for taxol biosynthesis.


Background
Taxol (generic name: paclitaxel), was first identified in the bark of Taxus brevifolia (Pacific yew) and is commonly used as a clinical drug for several types of cancer, such as lung cancer, and other solid tumors [1]. The substance is mainly extracted from the leaves and stems of Taxus spp., such as Taxus media. T. media is the natural hybrid of the maternal Taxus cuspidata and paternal Taxus baccata and contains considerable amounts of taxol in nearly all of its tissues [2]. However, although T. media has been extensively cultivated, current taxol supplies cannot meet clinical requirements due to limited land resources and the low content of the desired products in the collected plant tissues.
Induction of callus cells of Taxus spp. is considered a promising means to produce taxol and can prevent the severe misuse of conventional plant species [3]. However, most Taxus calli newly induced from tissues turn brown and cease to grow within 6-12 months, and this tissue often requires approximately 1-2 years or more to grow for industrial use [4]. Cells cultured through long-term suspension subculture show drastically reduced taxol yields and exhibit numerous other problems, such as heterogeneity and cellular ploidy [5][6][7]. These issues inhibit the extensive industrial use of Taxus cells for taxol production. Therefore, clarifying the regulatory mechanism of taxol biosynthesis is beneficial to address these challenges.
In a previous study, we investigated the long-term (10 years) subculture of Taxus cells and found that longterm subculture changes the living behavior of Taxus cells, resulting in reductions in taxol and other secondary metabolites and strengthening of their primary metabolism. These results indicate a vague regulatory mode in which long-term subcultured cells undergo decrease taxol biosynthesis [7].
Taxus suspension cells require approximately 1-2 years or more of subculture before they can be used for industrial culture. During the induction and formation of suspension cells, newly induced Taxus calli, which develop within 6 months of initiation of induction, contain more taxol than their parent tissues [4]. To clarify the regulatory mechanism of taxol biosynthesis in Taxus cells, we applied multiomics technology to compare callus cells with parent tissues.
Multiomics technologies have been applied to clarify complicated regulatory networks; these technologies can reveal the details of regulatory mechanisms and validate the results of high-throughput sequencing analysis [8][9][10]. MiRNAs have previously been recognized to be key regulators of taxol biosynthesis and transcriptional reprogramming in Taxus cells [7]. Moreover, several miRNAs regulate TFs and epigenetic factors in tissue dedifferentiation in many plants [11][12][13]. For example, Wu et al. found that miRNAs direct DNA methylation at loci in which they are produced, as well as in trans to their target genes, and play roles in gene regulation [14]. Shen et al. confirmed that miRNAs participate in dedifferentiation by regulating key functional genes enriched in the pathways of plant hormone signal transduction [15]. Therefore, miRNAs may be crucial factors in the regulatory mechanisms of taxol biosynthesis in callus cells.
In this work, multiomics analyses were conducted to detect the reprogrammed transcriptional profiles of newly induced callus cells, reveal the key pathways and factors influencing taxol biosynthesis and transcriptional alterations, and verify and screen key miRNAs and targets.

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 1464 bp were assembled, and the genes were annotated using the NR, Swiss-prot, COG, GO, and KEGG databases (Additional files 1 and 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" (Fig. 1a, b, 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 (Fig. 1a, Additional files 7 and 8).
Degradome sequencing revealed that 1829 unigenes were degraded by 347 miRNAs, leading to 2432 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][17][18][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 (Fig. 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 (Fig. 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, 10deacetyl taxol and baccatin III, were also high (Fig. 3b).
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 (Fig. 3c). The functions of these two genes requires further elucidation.
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 "planthormone signaling transduction," and "plant-pathogen interaction pathways" ( Fig. 4a and b, 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 (Figs. 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 (Figs. 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 costeffective 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 (Fig. 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 Fig. 1 Transcriptional alterations in callus cells Total six samples from two groups were high-throughput sequenced. Six samples were validated by Pearson correlation analysis (a), and they were obviously separated into two groups, callus cells and tissues. All differentially expressed genes were annotated with KEGG database, and the significant DE pathways were showed in (b). These DEGs were analyzed by Mapman3.6.0, and a metabolism overview were showed in (c) upregulated; both of these pathways are key processes for decomposing sugar for living energy (Fig. 1c, Additional file 5a, Additional file 12). Interestingly, "Photosynthesis" (35 downregulated vs. 3 upregulated) and "Photosynthesis-antenna proteins" (13 downregulated) were extremely weakened (Fig. 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" (Fig. 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 (Fig. 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][25][26][27][28][29][30]. Callus cells showed more downregulated genes involved in abiotic and biotic stress responses compared with tissues but showed upregulated heat/cold/light responses (Fig. 5a). For "Plant-pathogen interaction" and "Plant hormone signal transduction", 70% (434) and 62% (209) DEGs were downregulated ( Fig. 5b and c, Additional file 5a). For example, jasmonic acid-amino synthetase 1 (JAR1) Fig. 2 Function annotation of miRNA and degraded targets. Degraded DE genes were annotated to be involved in nearly all differentially expressed pathways (a). Moreover, the DE degraded genes significantly enriched in 19 pathways (b). Among the opposite expressed miRNAs and targets, there were six miRNAs and six degraded targets found to have the most targets and be targeted respectively (c). The DE miRNAs and their DE targets mostly enriched in 5 pathways, including "Plant-Pathogen Interaction", "Plant Hormone signaling transduction", "Ascorbate and aldarate metabolism", "Starch and sucrose metabolism" and "Aminoacyl-tRNA biosynthesis". FC was short for foldchange and nonexpresser of pathogenesis-related 1 (NPR1), which positively function in JA-and SA-signaling transduction, were significantly downregulated (Additional files 14 and 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 "Plantpathogen 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 (Fig. 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 1830 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 (Fig. 6a). A total of 156 miRNAs degraded more than one TF, and cme-MIR166e-p5_2ss9CT19GC degraded 71 TFs (Fig. 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 (Fig. 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 (Fig. 6b).

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 (Fig. 7a), and 346 of 451 TFs were highly correlated with taxol biosynthesis genes with coefficient values greater than 0.9 (Fig. 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 (Fig. 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 Fig. 4 Annotation of degraded targets of 10 miRNAs. These 10 miRNAs, which targeted to taxol biosynthesis genes, degraded 226 genes totally. And these degraded targets were functional annotated with GO (a) and KEGG (b), they were significantly involved in the pathways that callus cells mainly transcriptional reprogrammed. The pre-miRNA of Pc-5p-97202_13, which targeted to taxol biosynthesis genes was predicted to form a stable hairpin in structure (c) biosynthesis (Figs. 4c, d, and 7c). Although these coexpressed miRNAs do not degrade taxol biosynthesis genes, they may regulate taxol biosynthesis through other means.
"Plant-pathogen interaction" and the JA/SA-signaling pathways help plants defend themselves against pathogens and external stimuli [26,40]. Our results showed that callus cells dramatically downregulate the genes of these pathways, such as JAR1 and NPR1, which are positive regulators of the JA-and SA-signaling pathways, respectively (Additional files 5 and 15). Moreover, although several pathways did not show a significant change in callus cells, they tended to be downregulated. A total of 150 (81.7%) According to these results of increasing of biosynthesis of secondary metabolites and reducing of defense pathways, callus cells appeared to use secondary metabolites for defense instead of common defense pathways. Given that callus cells transition to heterotrophic feeding on substances in the medium, this alteration may be related to cost-effectiveness. Specifically, "Photosynthesis" was completely downregulated in callus cells, and the metabolism of several sugars and that of almost all amino acids contained in medium were enhanced. Rojas et al. indicated that the downregulation of some primary metabolic processes, such as photosynthesis, could save energy and promote defense responses; they also revealed that the upregulation of amino acid metabolism contributes to plant defenses [37]. Moreover, secondary metabolites could be used for immediate defense materials, whereas common defense pathways, such as "Plant-pathogen Interactions" and "JA/SA signaling transduction", require time to respond. Thus, accumulation of secondary metabolites is a cost-effective strategy for survival. Table 3 Degraded DE genes involving in biosynthesis of secondary metabolites FC is short for fold change (Callus/Tissues), inf and -inf indicated they were not detected in Tissues and Callus cells respectively. AO1 (abscisic aldehyde oxidase 1), DLO2 (DMR6-like oxygenase 2), FLS (Flavonol synthase), PLR1 (Pyridoxal reductase 1), SNL6 (Cinnamoyl-CoA reductase-like SNL6), PER25 (Peroxidase 25), CBP1 (Serine carboxypeptidase 1), SDR2b (Short-chain dehydrogenase/reductase 2b) and UGT85A24 (7- This inference is supported by a previous study on T. chinensis cells subjected to 10 years of subculture. In long-term subcultured cells, the metabolism of biological substances is enhanced, whereas common defense activities are attenuated. Remarkable reductions in taxol and other secondary metabolites suggest that long-term subculture cells abandon unnecessary defense activities [7]. These findings indicate that taxol accumulation is a survival strategy that incurs extra costs for wild plants. Thus, callus cells would prefer to choose a defense system with increased cost effectiveness.
MiRNAs are the tools used by callus cells to reprogram their metabolism. RNAs play decisive roles in changes in primary metabolism and defense systems, including biological substance metabolism, "plant-pathogen interactions", and "Endocytosis." Although miRNAs did not regulate the secondary metabolism of callus cells, they degraded a large number of TFs, such as MYB, ERF, bHLH, and other factors related to the biosynthesis of secondary metabolites. Moreover, a growing number of studies have revealed that regulation by miRNA-TFs is a main function mode in plants [11,41,42].

Conclusion
Newly induced Taxus callus cells exhibited reprogramed transcriptional profiles to survive and adopt to novel situations (Fig. 8). Increasing accumulation of taxol and other secondary metabolites appeared to be an active mode to compensate for weakened immune and defense systems because callus cells repress pathways such as "Plant-pathogen interaction" and "JA/SA-signaling transduction", which are crucial for plants to defend themselves against stresses. In addition, callus cells repressed activities such as "nucleotide/protein mismatch repair" and "spliceosome" but reinforced the metabolism and biosynthesis of biological substances. These alterations suggest that callus cells prefer a cost-effective survival strategy. Subsequent analysis revealed that miRNAs are key tools for the transcriptional reprogramming of callus cells. The reduction of defense and immune systems and increased metabolism and biosynthesis of essential substances mainly depended on miRNAs. Although a small proportion of the genes involved in the biosynthesis of secondary metabolites were degraded by miRNAs, miR-NAs regulated secondary metabolism by degrading a series of TFs. Several candidate TFs were highly correlated and coexpressed with taxol biosynthesis genes, and some TFs were regulated by miRNAs. Our results reveal the candidate regulators of TFs and miRNAs for taxol biosynthesis, as well as the significance of taxol biosynthesis for plant regulation networks.

Plant materials and callus induction
The 4-years-old Taxus media Rehder, obtained from Hubei green tide biotechnology company, were cultivated from cuttings at nursery in Huazhong University of Science and Technology. The tissues were fresh stems Fig. 6 Transcription factors were important targets of miRNAs. More than half of TF families were significantly degraded by miRNAs (a), and several TFs were degraded by more than one miRNAs so that degraded fragments were much more than the DE TFs, suggesting these TFs played crucial roles in transcriptional reprogramming in callus cells (b) and callus cells were induced with MS medium containing 0.3 mg/L 2,4-D, 0.4 mg/L NAA from the fresh stems from the same plant, these callus cells were maintained for 6 months with 62# (modified Gamborg's B5, the hormones concentration were 0.5 mg/L NAA, 0.5 mg/L 6-BA and 0.2 mg/L 2.4-D) medium at 25°C without daylight [43]. Six months culturing made callus cells form a stability growth state so that it was beneficial for us to screen the key regulate factors.
High-throughput mRNA sequencing, assembly and differential expression genes screening Total 6 samples belong to two groups (named as Tissues and Callus) were used to detect their transcriptional profiles. Total RNA isolation and library construction were conducted by LC-bio Tech (Hangzhou, China). The average insert size for the paired-end libraries was 300 bp (±50 bp), then the paired-end sequencing was performed on an Illumina Hiseq 4000 (1 × 150 bp for each end) platform following the vendor's recommended protocol. Firstly, Cutadapter V1.9 and fqtrim V0.94 were used to remove the reads that contained adaptor contamination, low quality bases and undetermined bases [44]. Then FastQC V0.10.1 was used to verify sequence quality (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/) including the Q20, Q30 and GC-content of the clean data. All downstream analyses were based on clean data of high quality, which did not contain any adaptor contamination (if the read was smaller than 100 bp after removing the adapter, the read were deleted), low quality bases (if the average quality of every 6 bp of single read is lower than 20, the entire 3′-flanking sequence was deleted) and undetermined bases (if the unknown base is more than 5% in a read, the read will be deleted). De novo assembly of the transcriptome was performed with Trinity 2.4.0. All assembled Unigenes were aligned against the non-redundant (Nr) protein database (http://www.ncbi.nlm.nih.gov/), Gene ontology (GO) (http://www.geneontology.org), SwissProt (http:// www.expasy.ch/sprot/), Kyoto Encyclopedia of Genes and Genomes (KEGG) (http://www.genome.jp/kegg/) and eggNOG (http://eggnogdb.embl.de/) databases using Fig. 7 Candidate TFs and miRNAs related to regulation of taxol biosynthesis. Previously, transcriptional profiles during long-term subculture process were detected. Thus, integrative analysis of these expression profiles with Callus&Tissues were helpful to identify the key transcription factors. All DE TFs were analyzed their expression patterns comparing with 12 known taxol biosynthesis genes (a), and the co-expression relation values of candidate TFs were showed in (b). The candidate miRNAs and co-expression values were showed in (c) DIAMOND [45] with a threshold of E-value<1e-5. Salmon was used to perform expression level for Unigenes by calculating TPM [46], and R package edgeR was applied for calculating statistical significance (p-value), then FDR were used to detected the false rate [47]. Since there were three independent bio-repeats of each individual group, the differentially expressed unigenes were selected not only with |log2 (fold change)| > 1, but also with FDR < 0.05. The GO and KEGG enrichment analysis of DEGs were validated by Fisher's Exact Test with a cutoff p-value< 0.05.

MiRNA analysis and targets prediction
Part of total RNAs from the same 6 samples of transcriptome were also used to conduct library for miRNAseq. The miRNA library was constructed with TruSeq Small RNA Sample Prep Kits (Illumina, San Diego, USA), and they were single-end sequenced by Illumina Hiseq2000/2500 with 1 × 50 bp. ACGT101-miR (LC Sciences, Houston, Texas, USA) was used to process raw reads by removing adapter dimers, junk (unknown nucleotide bases more than 10%), low complexity, common RNA families (rRNA, tRNA, snRNA, snoRNA) and repeats as described by Hu et al., 2018 [48]. Subsequently, unique sequences with length in 18~25 nucleotide were mapped to specific species precursors in miRBase 21.0 by BLAST search to identify known miRNAs that allowed length variation at both 3′ and 5′ ends and one mismatch inside of the sequence. For novel miRNAs, the unmapped sequences were BLASTed against the shallow sequencing Taxus baccata genomes [49], and the hairpin RNA structures containing sequences were predicted from the flank 120 nt sequences using RNAfold software (http://rna.tbi.univie.ac. at/cgi-bin/RNAfold.cgi).
Firstly, expression levels of miRNAs were based on normalized deep-sequencing counts, the normalization procedures were according to Zhao et al., 2017 [50]. Differential expression of miRNAs was analyzed by using Student T-test, when the p-value was less than 0.05, the miRNAs were considered as differentially expressed miRNAs [50] (Additional file 18).
To predict the genes targeted by most abundant miR-NAs, Target Finder ver. 50 were used to identify miRNA binding sites using default settings [51].

Degradome sequencing and data processing
Approximately 20 μg total RNA that mixed from 6 samples were used to prepare degradome library followed as with some modification [52]. Then, the libraries were sequenced using the 5'adapter only, resulting in the sequencing of the first 36 nucleotides of the inserts that represented the 5′ ends of the original RNAs. And the single-end sequencing (36 bp) were performed on an Illumina Hiseq2500.
After removing adaptors and low quality reads, the extracted sequencing reads were then used to identify  Figure were significantly differentially expressed, and the thickness of the arrows indicated the regulation strength by miRNAs. Some pathways also significantly changed but they were not regulated by miRNA, such as isoflavonoid biosynthesis. The pathways in green indicated most DEGs of these pathways were downregulated, red indicated most DEGs were upregulated, black means the upregulated genes of these pathways were as much as the downregulated ones. Most pathways of secondary metabolism were barely regulated by miRNAs directly, probably the regulation model was a "miRNA-TF-enzyme genes" mode potentially cleaved targets by the CleaveLand pipeline [53]. The degradome reads were mapped to the transcripts obtained above. Only the perfect matching alignment(s) for the given read would be kept for degradation analysis. All resulting reads (t-signature) were reverse complemented and aligned to the miRNA identified in our study. All the identified targets were subjected to BlastX analysis to search for similarity, and then to GO analysis to uncover the miRNA-gene regulatory network on the basis of biological process, cellular module and molecular function.

Taxanes extraction and quantitative analysis
The tissues and callus cells were firstly grinding with liquid Nitrogen, then taxanes were extracted using the method described previously [54]. And the quantitative analysis was conducted using HPLC as the previous reports [36].