Transcriptome analysis provides insights into key gene(s) involved in steroidal alkaloid biosynthesis in the medicinally important herb Fritillaria

Fritillariae Cirrhosae Bulbus, which is a traditional antitussive and expectorant medicine in China, is derived from six original plants. Among these six species, all except for Fritillaria taipaiensis and Fritillaria unibracteata var. wabuensis have been recorded in the National Protected Plants III. Their major bioactive compounds, including imperialine and peimisine, are steroidal alkaloids, derived from the steroid synthesis pathway. However, research on the F. taipaiensis genome, transcriptomes and steroid-synthesis-related genes has been lacking for a long time. In this study, the first transcriptome sequencing of F. taipaiensis was performed. A total of 94,396,694 (~11.4 Gb) high-quality reads obtained using paired-end Illumina sequencing were de novo assembled into 190,350 transcripts. Functional annotation with multiple public databases identified an array of genes involved in steroidal alkaloid biosynthesis and other secondary metabolite pathways, including steroid and terpenoid backbone biosynthesis, and important oxidoreductases. The large number of differentially expressed transcripts together with CYPs suggests the involvement of these candidates in tissue-specific expression. The expression analysis revealed that bulbs may be the main site of the upstream steroidal alkaloid biosynthesis pathway, and we speculated that the downstream reactions, including the oxidation-reduction reaction catalyzed by CYPs, might also occur primarily in bulbs.

multiple public databases identified an array of genes involved in steroidal alkaloid biosynthesis and other secondary metabolite pathways, including steroid and terpenoid backbone biosynthesis, and important oxidoreductases. The large number of differentially expressed transcripts together with CYPs suggests the involvement of these candidates in tissue-specific expression. The expression analysis revealed that bulbs may be the main site of the upstream steroidal alkaloid biosynthesis pathway, and we speculated that the downstream reactions, including the oxidation-reduction reaction catalyzed by CYPs, might also occur primarily in bulbs.

Conclusion
In conclusion, the comprehensive transcriptome dataset created in this study will serve as a resource for the identification of potential candidates for the genetic manipulation of targeted bioactive metabolites and will also contribute to the development of functionally 3 relevant molecular marker resources to expedite molecular breeding and conservation efforts for F. taipaiensis.

Background
Fritillariae Cirrhosae Bulbus is a kind of traditional Chinese medicine and a famous food that has been used in China for thousands of years due to its antitussive, expectorant and anti-asthma activities [1]. Fritillariae Cirrhosae Bulbus is derived from the dried bulbs of six species from the Fritillaria genus, including Fritillaria cirrhosa D. Don, F. unibracteata Fritillariae Cirrhosae Bulbus has increased dramatically [2]. The various plants of this medicine are mainly distributed in the western region of China, including Tibet, Yunnan, Sichuan, Qinghai and Gansu provinces, at an altitudinal range of 1800 to 4700 m [3]. Due to the increasing demand for Fritillariae Cirrhosae Bulbus and the long growing period of its plant constituents, many populations remain at risk of being depleted [4]. Among the six species mentioned above, all except for F. taipaiensis have been recorded in the National Protected Plants Ⅲ. Artificial cultivation may be an important way to resolve the current imbalance between supply and demand. According to a previous study, F. taipaiensisis the most suitable species for cultivation at low altitudes [5].
Previous investigations on the pharmacology and phytochemistry suggested that the steroidal alkaloids, especially isosteroidal alkaloids, were responsible for the antitussive, expectorant, anti-inflammatory and antinociceptive effects of these herbs [10][11][12][13][14]. F. taipaiensis is no different from F. cirrhosa as a commodity due to their similar chemical constituents [15,16]. Steroidal alkaloids are mainly synthesized from squalene through isopentenyl diphosphate (IPP) and dimethylallyl diphosphate (DMAPP) via cytosolic mevalonate (MVA) and plastidial 2-C-methyl-D-erythritol 4-phosphate (MEP/DOXP) pathways, respectively. Furthermore, several hydroxylation steps catalyzed by CYP450s contribute to the structural and functional diversification of various steroidal alkaloids [17,18]. However, due to the nonavailability of genomic resources, this medicine,F. taipaiensis, is less studied at the molecular level. Until now, the steroidal alkaloid biosynthesis pathway in this medicinal plant has been unclear.
RNA sequencing (RNA-Seq) is a common approach to transcriptome profiling that uses next-generation sequencing (NGS) technologies. Compared with other methods, this approach has several significant advantages, such as being able to reveal the precise location of transcription boundaries at a single-base resolution and avoiding background signals with the detection of digital signals, while requiring less sample RNA and not requiring the design of specific probes with existing genomic sequences. Thus, RNA-Seq has been shown to be an efficient tool for genome-wide transcriptome profiling and elucidation of important candidates involved in complex biosynthetic pathways, irrespective of genome complexity, even in the case of nonmodel plant species [19]. RNA-Seq has been used to discover key enzyme genes in various species. Meanwhile, many tools have been developed to increase the data credibility, such as trinity and cd-hit.
In this study, a comprehensive transcriptome of F. taipaiensis was sequenced using the RNA-Seq approach to identify key genes involved in the complex steroidal alkaloid pathway, and qRT-PCR was used to verify gene expression levels. Current data provide the first genome-wide transcriptional insights into the functions of steroidal alkaloid genes 5 and their differential expression in bulbs or leaves of F. taipaiensis of different ages. In the future, the outcome of available data will serve as a resource to expedite cutting edge research for the upscaling of targeted secondary metabolite production through genetic engineering of F. taipaiensis. Efforts were also made to identify potential transcription factors and CYPs that play key roles in the regulation and diversification of secondary metabolites.
Three genotypes for different year-old tissues located at a distance of 10 m from each other were randomly considered. Bulb or leaf tissues were collected from each genotype, and the details are presented in Table 1. The collected samples of F. taipaiensis were first washed under running tap water. Then, the surface of samples was sterilized with 75% ethanol for RNA-Seq. All samples were stored in RNA Later liquid (Beijing HT-biotech Co., Ltd., China) and frozen immediately in dry ice and stored at -80°C until RNA isolation.
Samples for HPLC-ELSD were dried at 45℃ in a drying oven (101FXB-2, Shanghai, China) until a constant weight was achieved. Subsequently, the dried samples were powdered to a homogeneous size, sieved through a No. 4 mesh and then stored at room temperature. RNA isolation, sequencing and de novo assembly Total RNA was extracted from individual samples using a Plant Universal Total RNA Kit (BioTeke Biotech Beijing Co., China) according to the manufacturer's instructions. The 6 concentration of RNA was determined using an Agilent 2100 Bioanalyzer system (Agilent Technologies Inc., USA), and integrity was checked on a denaturing agarose gel. Equimolar concentrations of RNA from three genotypes for each tissue were pooled together for RNA-Seq library preparation to remove biological bias.
RNA-Seq libraries were prepared using the NEBNext Ultra RNA Library Prep Kit for Illumina (New England Biolabs, Inc., UK) according to the manufacturer's instructions. The libraries were purified using Beckman AMPure XP beads (Beckman Coulter, Inc., USA), and an Agilent 2100 Bioanalyzer (Agilent Technologies, USA) was used for library size estimation.
An ABI 7500 real-time PCR system (Applied Biosystems Co., USA) with KAPA SYBR green fast qPCR Kits (Kapa Biosystems, USA) was used to quantify the libraries. Furthermore, for cluster generation, 10 pM of these libraries were loaded onto the flow cell of a cluster station using a TruSeq PE Cluster Kit (Illumina Inc., USA). Clonally amplified clusters were used for paired-end (PE) (2 × 125) sequencing using HiSeq 2500 with a TruSeq SBS Kit v3-HS (Illumina). Trimmomatic (v 0.30) was used to filter raw reads, and reads with 99% probability of no error (minimum phred score 20 for each base) were utilized for assembly [20]. . Trinity (v r20140717) was used for de novo transcriptome assembly with default parameters and a minimum transcript length of 200 base pairs [21]. Cd-hit was used to cluster biological sequences to reduce sequence redundancy and obtain unigenes [22].
TransDecoder (v r20130225), which comes included in the Trinity software distribution, was used to perform open reading frame annotation [21].

Functional annotation and classification
To find the putative functions of assembled transcripts of F. taipaiensis, a similarity search using BLASTx [23] was performed against publicly available protein databases, including NCBI nonredundant protein database (Nr) and Swiss-Prot with an e-value cut-off of ≤1e − 5 . F. taipaiensis transcripts were classified into three major categories: biological process, cellular component and molecular function according to Gene Ontology (GO) terms using WEGO software (http://wego.genomics.org.cn/). Transcripts were further functionally categorized into different classes using the COG database

Identification of differentially expressed genes
To measure the expression pattern of transcripts in each tissue, high-quality reads were mapped onto the final de novo assembled transcripts of F. taipaiensis using RSEM v 1.2.4 [24]. The expression level for each transcript was measured in terms of fragments per kilobase per million reads (FPKM) by normalizing read counts. Furthermore, the edgeR package was used to evaluate the differential gene expression using read counts in pairwise comparisons of different year-old bulbs [25]. Based on statistical analysis, genes having a p-value < 0.05 and log2-fold change ≥ 2 were considered differentially expressed genes. A heatmap representing the tissue-specific gene expression pattern (log2-fold change) for different pathways was generated using Multiple Experiment Viewer (MEV v4.9.0).

Quantitative real-time PCR (qRT-PCR) analysis
gDNA Eraser (Takara Bio Co., USA) treatment was performed on total RNA to remove DNA contamination. First-strand cDNA was synthesized from 1 μg of the total RNA using the PrimeScript RT reagent kit (Takara, Japan) according to the manufacturer's instructions and was diluted by 10X. Gene-specific primers ( Table 2) were designed using Primer3 Input (http://primer3.ut.ee) and synthesized by SinoGenoMax Co., Ltd, China. The relative expression of steroidal alkaloid pathway genes in seven samples was analyzed on a CFX96 Real-time System (Bio-Rad Lab., USA) using SYBR ® Premix Ex Taq™ (Tli RNaseH Plus) (Takara, Japan). qRT-PCR was performed with three technical replicates based on which standard error was calculated. β-actin was used as a reference gene for establishing equal amounts of cDNA in each reaction. The relative gene expression and fold change was calculated using the 2 − ΔΔCT method with rhizome as a control tissue [26]. Other reagents were of HPLC grade or the highest grade commercially available.

Solution preparation
Stock solutions (500 μg/ml) of the alkaloid standards peimisine, sipeimine, peimine and peiminine were prepared. Total alkaloids of samples were prepared as follows: the dried bulb powders were soaked in 4 mL ammonium hydroxide for 2 h and then reflux extracted with 60 mL of a mixture of chloroform:methanol (80:20, v/v) at 80℃ for 3 h. After extraction, the total extracted solution was filtered through qualitative filter paper.
Subsequently, the filtrate was concentrated to dryness in an evaporating dish under a fume hood. The residue was dissolved to exactly 2 mL with methanol using a volumetric 9 flask. The solution above was filtered through a 0.22 μm membrane filter prior to injection into the HPLC system. Instrumentation and HPLC conditions for determining alkaloids in F. GO has been widely used for functional analysis and inferring the biological significance of genomic and transcriptomic datasets, in which there are three ontologies: molecular functions, biological processes and cellular components. A total of 51,212 transcripts were assigned 4954 GO terms, of which 1929 terms were categorized into molecular functions, 2529 were categorized into biological processes and 496 were categorized into cellular components. Among the molecular functions, the most abundant of the second-level GO terms was related to catalytic activity (69.35%), followed by binding (12.80%) and transporter activity (11.05%). Among the cellular components, organelle parts (26.80%) and cell parts (26.30%) were the most represented in the second level, followed by protein-containing complexes (23.64%) and membrane parts (11.00%). In biological processes, cellular processes (22.37%) and metabolic processes (20.66%) were abundant, followed by single-organism processes (19.28%) and biological regulation (16.04%) ( Figure   2). were clearly the main categories, followed by signal transduction mechanisms (133 DEGs, 7.33%) ( Figure 5). A total of 1,453 DEGs were assigned 1438 GO terms, of which 546 terms were categorized into molecular functions, 172 were categorized into biological processes and 720 were categorized into cellular components. Interestingly, some CYP 450 genes were found to be differentially expressed, including CYP71 and CYP734 families.
Furthermore, a total of 520, 891 and 985 unique DEGs were obtained in one-year-old bulbs 13 vs two-year-old bulbs, one-year-old bulbs vs four-year-old bulbs and two-year-old bulbs vs four-year-old bulbs, respectively.

Secondary metabolic pathway analysis
Metabolic pathway analysis enables us to understand the interactions of genes in particular pathways and their related biological functions. A total of 19 pathways (174 entries) related to secondary metabolite biosynthesis were identified from the KEGG database. The identification of these pathways helped us to analyze secondary metabolite biosynthesis in F. taipaiensis. Of these, five major pathways, carotenoid, steroid, ubiquinone and other terpenoid-quinone, terpenoid backbone, and phenylpropanoid biosynthesis pathways, were well represented in our data. Out of 515 transcripts involved in these pathways, 33 recorded specific differential expression in bulb tissues from different growth years. Furthermore, most of the genes involved in these pathways were highly expressed in one-year-old bulbs.

Steroidal alkaloid pathway gene expression and content analysis
In plants, steroidal alkaloids are mainly synthesized from lanosterol and cycloartenol via cholesterol and sitosterol, respectively [28,29]. . The genes related to steroidal alkaloids via cholesterol have not been characterized so far; therefore, genes related to steroidal alkaloid biosynthesis via sitosterol were considered in this study, which comprises two parts: terpenoid backbone and steroid biosynthesis, according to KEGG classification. In total, 14 genes corresponding to 22 transcripts involved in the steroid biosynthesis pathway were identified in our study. In this biosynthesis pathway, farnesyl diphosphate farnesyltransferase (FDFT1), which is also called squalene synthase (SQS), catalyzes the first two-step reaction in which two identical molecules of farnesyl pyrophosphate (FPP) are converted into squalene. Then, squalene can be converted into 2,3-oxidosqualene based on oxidation by squalene epoxidase (SQLE), which is the branching point between triterpenoid and steroid [31] . During steroid biosynthesis, the chair-boat-chair-boat conformational change of 2,3-oxidosqualene results in the formation of cycloartenol, which ultimately leads to the synthesis of steroidal alkaloids through various modifications catalyzed by CYPs, isomerase, methyltransferases, reductase, etc.
To validate the RNA-Seq data and study the specific expression among different tissues and growth years, 19 genes involved in the steroidal alkaloid biosynthesis pathway were selected for quantitative real-time PCR (qRT-PCR) analysis, and the results are shown in Figure 6. In terpenoid backbone biosynthesis, the expression levels of all MEP/DOXP pathway genes, including DXS, DXR, ispD, ispE, ispF, ispG and ispH, were found to be clearly higher in leaves than in bulbs, whereas six genes involved in the MVA pathway, including ACAT, HMGR, HMGS, MVK, PMVK and MVD, showed diametrically opposite trends.
However, the IDI gene in the next step, which can convert DMAPP and IPP to each other, showed no clear trend in bulbs or leaves. Interestingly, GGPS and FDPS were described as having the same reaction to FPP, but they showed completely different trends: GGPS was highly expressed in leaves, whereas FDPS was clearly more highly expressed in bulbs than in leaves. In the steroid biosynthesis pathway, FDFT1 and SQLE were both highly expressed in leaves, while the expression level of CAS was similar in leaves or bulbs, and it was relatively high in the three-year-old bulb.
The HPLC-ELSD method was used in this study to quantify 4 alkaloids from 11 samples of F. taipaiensis in 5 growth stages; the analyte contents are shown in Table 3. Among the four alkaloid standards, only peimisine was detected in bulbs of F. taipaiensis, and the content ranged from 0.0230% to 0.0660%. The quantitative results showed that the peimisine content in F. taipaiensis increased with age from one to three years. However, when the plants were cultivated for over 4 years, the content decreased with age. Therefore, the peimisine content of the three-year-old bulb was the highest in either the 2016 or 2017 samples. Meanwhile, some samples could be regarded as the same sample taken in different years due to the continuous sampling within two years, such as three- year-old bulbs in 2016 and four-year-old bulb in 2017s. A longitudinal comparison was carried out within these samples (Figure 7). The results showed that the peimisine content in the three-year-old bulb (2016) was significantly higher than that of the four-year-old bulb (2017). Meanwhile, alkaloids in four-year-old leaves were extracted and detected, but no peak was obtained, which means that the steroidal alkaloid content in leaves was below the detection limit.
Cytochrome P450s and some other oxidoreductase genes Plant cytochrome P450 (CYP450s) is a B-type cytochrome protease gene superfamily with heme as a prosthetic group. CYP450s catalyze a wide variety of monooxygenation/hydroxylation reactions in metabolic reaction pathways, including those for cinnamates, fatty acids, hormone homeostasis, lignin, terpenoids, sterols, alkaloids, and saponins [32] . A total of 205 CYP genes corresponding to 620 transcripts classified under 58 families were identified in F. taipaiensis. Among the various CYPs, the CYP71 family was the most abundant. In total, 42 CYP genes related to 17 families were found to be differentially expressed in at least one pairwise comparison. The tissue-specific expression of these CYP genes revealed that 21, 13, and 8 genes were highly expressed in one-year-old, two-year-old, and four-year-old bulbs, respectively. Furthermore, some modification genes that may be involved in the steroidal alkaloid biosynthesis pathway were analyzed by qRT-PCR ( Figure 8). Sterol 14-demethylase (CYP51G1), steroid 22-alpha-hydroxylase (CYP90B1, DWF4), sterol delta-14 reductase (FK), delta (24)-sterol reductase (DIM, DWF1) and 7-dehydrocholesterol reductase (DWF5) were highly expressed in bulbs. Especially for CYP90B1 and DWF1, both genes were barely expressed in leaves but had high expression levels in the bulbs, and the expression levels were fifty to one thousand times that of the leaves. Interestingly, similar to the CAS gene, CYP90B1 showed the highest expression in three-year-old bulbs in both 2016 and 2017 samples. Moreover, the five genes above shared one common trait: the expression level was lowest in the perennial bulbs. In addition, acyl-CoA synthetase (ACS12) and CYP734A6 were more highly expressed in leaves than bulbs, but 3-epi-6-deoxocathasterone 23monooxygenase (CYP90D1) had no clear trends of expression levels in bulbs or leaves.
The qRT-PCR result matches the RNA-Seq sequencing result by 55.56%.

Transcription factors
TFs are major regulatory elements that play significant roles in gene expression, plant secondary metabolism and response to environmental stress by binding to specific cisregulatory elements of the promoter regions. TF families, including bHLH, bZIP, MYB, NAC, and WRKY, were reported to be involved in the regulation of secondary metabolites and abiotic and biotic stress responses in many plant species [33,34] . In our study, a total of 1,362 transcripts were assigned to 62 TF families. Among these, WRKY (165) The RNA-Seq approach based on NGS technology has a wide variety of applications, 18 including in plant, animal and human cells [36][37][38][39][40] . Moreover, data analysis software has been developed for de novo assembly, which overcame the shortcomings of short reads, making assembly of transcriptome results from nonreference genomic species reliable [21] . In this case, the time was ripe for transcriptome sequencing of F. taipaiensis, and we present the first Illumina-based RNA-seq derived transcriptome for this species. The pairwise reads (33 million to 37 million) obtained from this study are sufficient for reliable de novo transcriptome characterization and accurate quantification of gene expression patterns.
In plants, a variety of metabolites are derived formally from isoprenoid units, including terpenoids, steroidal alkaloids, carotenes, phytol and gibberellin. Isosteroidal alkaloids, belonging to steroidal alkaloids, are the main pharmacodynamic components of F. taipaiensis. In higher plants, IPP, as the common precursor of all isoprenoid compounds, is synthesized via two different pathways: MVA in the cytoplasm and MEP/DOXP in plastids [41]. From the leaf and bulb samples, it was evident that chloroplasts are common to leaves, but few are found in bulbs; this difference may explain why the genes in the MEP pathway were upregulated in leaves, whereas genes in the MVA pathway showed higher expression levels in bulbs.
FDPS can catalyze the synthesis of geranyl diphosphate (GPP), FPP and geranylgeranyl pyrophosphate (GGPP) from IPP and DMAPP, which can be further synthesized into various metabolites. It is generally accepted that the MVA pathway mainly synthesizes sesquiterpenoids, triterpenoids and steroids [42] ), while the MEP pathway is mainly related to the synthesis of monoterpenes, diterpenes, phytol, gibberellin, and carotenoid.
In this study, the results of the HPLC-ELSD experiment demonstrated that isosteroidal alkaloids mainly accumulate in bulbs of F. taipaiensis. Furthermore, previous studies have found that genes in the MEP pathway play a major role in coordinating the synthesis of chloroplast-localized isoprenoid-derived compounds, such as carotenoid and chlorophyll [43] . Thus, it is speculated that the IPP required by the isosteroidal alkaloids mainly comes from the MVA pathway in bulbs of F. taipaiensis.
The CAS gene, which catalyzes the formation of the sterol cycloartenol by cyclization of 2,3-epoxysqualene, is a key enzyme gene in the steroid biosynthesis pathway [44,45] .
This catalytic reaction could provide precursors for the biosynthesis of sterols and triterpenoids and is considered to be a key site for the biosynthesis branch of sterols and triterpenoids. More interestingly, the relative peimisine content was highest in the threeyear-old bulbs, which matched the expression trend of the CAS gene in this study. Thus, it could be speculated that the high expression of enzymes related to the steroid biosynthesis pathway leads to a high accumulation of steroidal alkaloids in the third year in bulbs.
The precursors, such as cycloartenol, need to undergo a series of redox reactions to form the secondary metabolites. Thus, some oxidoreductases, such as the cytochrome P450 enzymes (CYP450), play an essential catalytic role in redox reactions. The CYP450 superfamily is a large and diverse group of membrane-bound hemoproteins belonging to monooxygenase, which is the first enzyme class to be classified as a superfamily [46,47]; it is also an ancient superfamily, including more than 1,000 families and 2,500 subfamilies [48] . In plants, CYP450s have multiple subfamily members, including CYP51, CYP71, CYP99, CYP701, CYP736, etc. CYP90B1/DWF4, which belongs to family 90, subfamily B, has been confirmed to catalyze the C22-alpha-hydroxylation step in brassinosteroid biosynthesis by converting campestanol to 6-deoxocathasterone and 6-oxocampestanol to cathasterone. Fujita et al. confirmed that CYP90B1 could catalyze the C-22 hydroxylation of various C 27 , C 28 and C 29 sterols [49] . Interestingly, this enzyme gene was discovered in this study, and qRT-PCR results demonstrated that it was only highly expressed in bulbs, and the expression level corresponded with the trend of peimisine accumulation. Thus, we hypothesize that this gene may catalyze C-22 hydroxylation in the biosynthesis of F. taipaiensis alkaloids. CYP51G1, a type of multifunctional oxidase known to be involved in sterol and steroid biosynthesis, was also identified in this study [50,51] .
CYP51G1/CYP51A2, known as sterol 14-demethylase, is involved in sterol biosynthesis, belongs to family 51, subfamily A and catalyzes the reaction removing methyl from the C-14 of certain phytosterols. Kim et al. noted that this enzyme plays a crucial role in controlling plant growth and development through a sterol-specific pathway [52].
In addition to the CYP450 enzyme genes, some other oxidoreductase genes have also been verified by qRT-PCR, and DWF1/DIM has been the most prominent among all genes. DWF1/DIM is a bifunctional protein that is necessary for both the isomerization and reduction in 24-methylenecholesterol in the brassinosteroid biosynthesis pathway [6] .
Sumeer et al. noted that increasing the transcript levels of DWF1 could enhance the accumulation of withanolides-a group of steroidal lactones, possibly by affecting downstream phytosterol biosynthesis [53] . According to the significant difference in the expression levels of this gene between different tissues, we suspect that DWF1 also has a certain function in the steroidal alkaloid synthesis pathway of F. taipaiensis. However, no studies have shown whether these genes will play a role in the biosynthesis of peimisine and other alkaloids. Thus, we will continue to focus on whether this occurs and what roles the above enzyme genes may play in the synthesis of sterol alkaloids in subsequent studies of F. taipaiensis.
As shown before, relative steroid accumulation is highest in the third year and decreases in the fourth year. However, in actual production of medicinal materials, we must consider not only the content of active ingredients per gram but also the final yield of the medicine. Thus, biomass must be an important indicator. As shown in Figure 9, the four-year-old bulbs of F. taipaiensis are at least twice the size of three-year-old bulbs. Thus, the reason why the alkaloid content in the four-year-old bulbs was relatively low may be that the accumulation of secondary metabolites cannot keep up with the increasing bulb mass. For production, the larger mass and relatively high secondary metabolite accumulation are more suitable for harvesting. Based on the above results and discussion, we recommend harvesting F. taipaiensis in the fourth year.

Conclusions
Overall, the lack of a reference genome for F. taipaiensis has made it difficult to determine the exact number of genes involved in the metabolic pathways of this plant species. Here, we conducted Illumina sequencing to profile the transcriptome of F. taipaiensis. To the best of our knowledge, this is the first attempt at using Illumina PE sequencing technology to produce bulb and leaf transcriptomes of Fritillariae Cirrhosae Bulbus through de novo sequencing and assembly without a reference genome. A large number of candidate transcripts were matched with known enzymes in public databases, providing a substantial gene resource for further research in F. taipaiensis. Our study also provides reference sequences for evolutionary analyses of metabolomes among the members of the Fritillaria family, which is rich in ethnomedicinal plants.

Declarations
Ethics approval and consent to participate Not applicable Consent for publication Not applicable Availability of data and material The datasets generated during the current study are available in the NCBI Sequence Read

22
Archive database (The data is being uploaded).

Competing interests
The authors declare that they have no competing interests.    Size comparison of bulbs of different ages