Transcriptome Analysis of Secondary Metabolism Pathway, Transcription Factors, and Transporters in Response to Methyl Jasmonate in Lycoris aurea

Lycoris aurea, a medicinal species of the Amaryllidaceae family, is used in the practice of traditional Chinese medicine (TCM) because of its broad pharmacological activities of Amaryllidaceae alkaloids. Despite the officinal and economic importance of Lycoris species, the secondary mechanism for this species is relatively deficient. In this study, we attempted to characterize the transcriptome profiling of L. aurea seedlings with the methyl jasmonate (MeJA) treatment to uncover the molecular mechanisms regulating plant secondary metabolite pathway. By using short reads sequencing technology (Illumina), two sequencing cDNA libraries prepared from control (Con) and 100 μM MeJA-treated (MJ100) samples were sequenced. A total of 26,809,842 and 25,874,478 clean reads in the Con and MJ100 libraries, respectively, were obtained and assembled into 59,643 unigenes. Among them, 41,585 (69.72%) unigenes were annotated by basic local alignment search tool similarity searches against public sequence databases. These included 55 Gene Ontology (GO) terms, 128 Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways, and 25 Clusters of Orthologous Groups (COG) families. Additionally, 4,175 differentially expressed genes (DEGs; false discovery rate ≤ 0.001 and |log2 Ratio| ≥ 1) with 2,291 up-regulated and 1,884 down-regulated, were found to be affected significantly under MeJA treatment. Subsequently, the DEGs encoding key enzymes involving in the secondary metabolite biosynthetic pathways, transcription factors, and transporter proteins were also analyzed and summarized. Meanwhile, we confirmed the altered expression levels of the unigenes that encode transporters and transcription factors using quantitative real-time PCR (qRT-PCR). With this transcriptome sequencing, future genetic and genomics studies related to the molecular mechanisms associated with the chemical composition of L. aurea may be improved. Additionally, the genes involved in the enrichment of secondary metabolite biosynthesis-related pathways could enhance the potential applications of L. aurea.


INTRODUCTION
Plants produce the crucial chemicals for growth and development by their own, including the primary and secondary metabolites. The primary metabolites include carbohydrates, acids, amino acids, fat and so on. The secondary metabolites, also named as natural products or phytochemicals, are specific to some taxonomic groups, playing pivotal roles in interactions between the plant and environment, helping plants to defense against pathogens or herbivores (Dixon, 2001;Kennedy and Wightman, 2011;Kliebenstein and Osbourn, 2012). As the sources of drugs, insecticides, and flavor, most of the plant secondary metabolites are very important for human's daily life (Goossens et al., 2003;Hussain et al., 2012;Yang et al., 2012).
Lycoris aurea, belonging to Amaryllidaceae family, is a medicinally and ornamentally important species. It is rich of the secondary metabolites such as Amaryllidaceae-type alkaloids, and is widely used for traditional Chinese medicine (TCM). The Amaryllidaceae-type alkaloids which are isolated from Lycoris genus have been reported to exhibit immunostimulatory, antimalarial, tumor, and viral activities (Jin, 2009). For example, as one of the typical Amaryllidaceae alkaloids, galanthamine is a kind of reversible inhibitor of cholinesterase to increase acetylcholine sensitivity, and it has also been clinically used in the treatment of Alzheimer's disease (Harvey, 1995;Bores et al., 1996). Despite of the officinal, economic and cultural importance of Lycoris species, the secondary mechanism for this species are relatively limited.
The experimental approach based on sequencing the functional genomics was reported to facilitate gene discovery in plant secondary metabolism (Dixon, 2001;Goossens et al., 2003). Besides, RNA-sequencing (RNA-Seq) technology was used to obtain full-scale transcriptomic information from different plant species such as tea plant, Polygala tenuifolia, Chlorophytum borivilianum, and Atractylodes lancea, and provide a better insight into transcriptional and post-transcriptional regulation of the essential genes in the secondary metabolite biosynthetic pathways (Kalra et al., 2013;Li et al., 2015;Devi et al., 2016). Regarding to Amaryllidaceae-type alkaloids biosynthesis pathway, Kilgore et al. (2014) defined a 4 ′ -O-methyltransferase which is involved in the biosynthesis of galanthamine by sequencing transcriptome of Narcissus sp. aff. Pseudonarcissus. Previously, the de novo transcriptome was sequenced to produce the EST (comprehensive expressed sequence tag) dataset for Lycoris aurea, which provides one perspective of the regulatory and synthesized molecular mechanisms of Amaryllidaceaetype alkaloids . On the other hand, the inspection of the comparative transcriptome files between two databases provides another method to investigate the Abbreviations: DEG, differentially expressed gene; FDR, false discovery rate; RNA-Seq, high throughput sequencing of cDNA libraries; cDNA, complementary DNA synthesized from RNA; qRT-PCR, quantitative real-time polymerase chain reaction; GO, gene ontology; KEGG, Kyoto encyclopedia of genes and genomes; Nr, non-redundant database; COG, cluster of orthologous groups databases; CYP1, cyclophilin 1. NCBI, National Center for Biotechnology Information; MS, Murashige and Skoog; FPKM, Fragments per kilobase of exon per million fragments mapped reads; ABC, ATP-Binding Cassette. secondary metabolites biosynthesis, or/and find the objective genes involved in Salvia miltiorrhiza, Catharanthus roseus, sweet cherry, Atractylodes lancea, and so on (Guo et al., 2013;Yu and Luca, 2013;Wei et al., 2015;Huang et al., 2016).
A great diversity in the types of microbial, physical, or chemical factors, known as elicitors, could influence the levels of secondary metabolites in plants (Zhao et al., 2005;Vasconsuelo and Boland, 2007;Jimenez-Garcia et al., 2013;Verma and Shukla, 2015). For example, the elicitors such as nutrient supply, temperature, metal ions, light conditions and atmospheric CO 2 concentrations affect the production of secondary metabolites in plants (Nascimento and Fett-Neto, 2010;Gandhi et al., 2015). The elicitors derived from bacteria, fungi, viral pathogens and even plants also contribute to the variability in plant secondary metabolism (Berenbaum, 1995;Verma and Shukla, 2015). Previously, we showed that exogenous methyl jasmonate (MeJA) application accelerated the Amaryllidaceae alkaloids accumulation in Lycoris chinensis seedlings (Mu et al., 2009). In this study, by using elicitor MeJA treatment, the global expression patterns of genes involved in metabolism, particularly secondary metabolism, transcription factors, and transporter proteins were identified. Therefore, this transcriptome sequencing may help improve future genetics and genomics studies on molecular mechanisms associated with the secondary metabolites of L. aurea.

Plant Growth and Treatment
L. aurea seeds were collected from Institute of Botany, Jiangsu Province and Chinese Academy of Sciences, Nanjing, China. The seeds were surface sterilized with 75% alcohol (v/v), and germinated on half-strength Murashige and Skoog (MS) medium (pH 5.8) in the dark at room temperature for 10 days. Afterwards, the seedlings were transferred into plastic pots containing a mixture of soil and vermiculite (3:1, v/v) and cultured in a growth chamber under 14 h light (25 • C)/10 h dark (22 • C). After 12 months growth, the seedlings were treated with 100 µmol L −1 methyl jasmonate (MJ100) for 6 h. MeJA was dissolved with 1% DMSO (v/v) to prepare the stock solution. Seedlings grown in MeJA-free solution (1% DMSO) were used as control (Con). The seedlings were harvested and immediately frozen in liquid nitrogen and stored at −80 • C.

RNA Isolation, cDNA Library Construction and Illumina Sequencing
Total RNA of the samples were extracted using RNAiso Plus reagent (Takara Bio, Dalian, China) following the manufacturer's instruction. RNA samples were examined with a spectrophotometer (Thermo Fisher Scientific, Inc. Waltham, MA, USA) and electrophoresed on a 1% agarose gel. The construction of the cDNA libraries and the RNA-Seq assay were performed by the OE Biotech Company (Shanghai, China).
Poly (A) mRNA was enriched referring to the previous method  by using NEBNext R Poly(A) mRNA Magnetic Isolation Module (New England Biolabs, Ipswich, MA, USA), and fragmented to short pieces. These short fragments were as applied as the templates for cDNA. The cDNAs were then subjected to end-repair using T 4 DNA polymerase and phosphorylation using Klenow DNA polymerase. Then, a base "A" was added to the 3 ′ ends of the repaired cDNA fragments. All of the short fragments were linked to sequencing adapter. The resulting fragments were selected as the PCR templates after electrophoresis. Finally, the four libraries were sequenced using Illumina HiSeq TM 2,000.

Transcriptome Assembly and Functional Annotation
The raw data of the RNA sequencing were purified by trimming adapters, removing reads containing poly-N, and rejecting the low-quality data (quality value ≤ 10 or unknown nucleotides larger than 5%) to get the clean reads. Meanwhile, the proportion of nucleotides with quality values greater than 20 (Q20) and GC content of the clean data were calculated. Then, all of the clean reads were assembled by using Trinity program (Grabherr et al., 2011). Firstly, for each library, the certain short reads with overlap regions were assembled into longer contiguous sequences (contigs). Then, based on the paired-end information, the distance of different contigs was recognized by mapping the clean reads, to obtain the sequence of the transcripts. Finally, performing the sequence of potential transcript to TGI Clustering tool, the unigenes were obtained (Pertea et al., 2003). For gene functional annotation, all of the assembled unigenes were aligned to the public databases, including National Center for Biotechnology Information (NCBI) non-redundant protein (Nr) and nucleotide (Nt) database, the Swiss-Prot protein database, Gene Ontology (GO, http://wego.genomics.org.cn/cgi-bin/wego/ index.pl) database, Clusters of Orthologous Groups (COGs) database, and the Kyoto Encyclopedia of Genes and Genomes (KEGG, http://www.genome.jp/kegg/kegg2.html) database, with a cut-off of E ≤ 10 −5 .

Differentially Expressed Genes (DEGs) Analysis
The expression level of unigene was calculated following fragments per kilobase of exon per million fragments mapped reads (FPKM) method (Mortazavi et al., 2008). To identify DEGs between two groups, the ratio of the FPKM values (using 0.001 instead of 0 if the FPKM was 0) were taken as the fold-changes in the expression of each gene. In order to compute the significance of the difference in transcript abundance, the false discovery rate (FDR) control method was used to identify the threshold of the P-value in multiple tests (Reiner et al., 2003). In this result, only fold change with|log2 (MJ100_FPKM/Con_FPKM)|≥ 1, and an FDR ≤ 0.001 were taken as the threshold for significantly differential expression. The log 2 -transformed FPKM value for DEGs was applied to generate heat map by MeV 4.7 (Howe et al., 2011). Meanwhile, the DEGs were annotated with GO and KEGG databases.

Validation of DEGs with Quantitative Real-Time PCR (qRT-PCR)
qRT-PCR was used to confirm the expression of MeJA-responsive genes in L. aurea. Total RNA was extracted from tissue sample using the method described above. The first-strand cDNA was synthesized using a PrimeScript TM II First Strand cDNA synthesis kit (Takara Bio, Dalian, China) according to the manufacturer's protocol. The primers were designed using the software Primer Premier 5.0, and listed in Table S1. The quantified expression levels of the tested genes were normalized against the housekeeping genes Cyclophilin 1 (CYP1) (Ma et al., 2016). qRT-PCR was performed using SYBR Premix Ex Taq TM II kit (Takara) and run on qTOWER2.2 Real-Time PCR System (Analytik Jena AG, Jena, Germany). Conditions for quantitative analysis were as follows: 94 • C for 2 min; 35 cycles of 94 • C for 15 s, 60 • C for 20 s, 72 • C for 10 min. Data for each sample were calculated using 2 − CT method (Livak and Schmittgen, 2001).

Transcriptome Sequencing Profile of L. aura under MeJA Treatment
To develop a comprehensive overview of the L. aurea transcriptome under MeJA treatment, two Solexa/Illumina libraries Con and MJ100, were designed for RNA-Seq. These two libraries (Con and MJ100) produced 4.82 Gbyte and 4.65 Gbyte of clean data, respectively. In addition, paired-end reads of each library are with GC percentage and Q20 percentage of 97.52 and 45.24%, 97.42 and 44.83%, respectively. Subsequently, short reads with average lengths of 244 and 273 bp of the two libraries were collected into 181,881 and 156,621 contigs, respectively ( Table 1). Taking the distance of paired-end reads into account, these contigs were assembled into non-redundant unigenes (Table S2). In total, 59,643 unigenes in the range of 250-14,908 bp (with an average length of 740 bp) were obtained (Figure 1).

Sequence Annotation and Classification
All of the unigenes were annotated by BLAST search in the public databases ( number of unigenes was annotated in each GO term shown in Figure 2 and Table S3. Based on the similarity to sequences with known functions, 115,216 sequences were assigned to the biological process (BP) category, 95,500 to the cellular component (CC) category, 34,965 to the molecular function (MF) category. Because several unigenes were assigned to more than one category, the total assigned number was bigger than the total number of the unigenes. In MF category, genes assigned to "catalytic activity" (15,547) and "binding" (14,379) constituted the largest proportion, accounting for 85.59% of the total. Within the CC category, "cell" (23,719), "cell part" (23,717) and "organelle" (19,872) were highly represented. Moreover, "cellular process" (19,233) and "metabolic processes" (18,289) were the main groups under BP category (Figure 2; Table S3). For functional prediction and classification, all unigenes were subjected to BLAST search against the COG database. There were 15,370 unigenes assigned to COG functional classification and divided into 25 specific categories based on the homology (Table 2; Figure 3). The "General functional prediction only" (4,997) was the largest category, followed by "Transcription" (3,056), "Replication, recombination and repair" (2,562), "Post-translational modification, protein turnover, chaperones" (2,362) and "Signal transduction mechanisms" (2,000). "Nuclear structure" (7) and "Extracellular structures" (6) were the smallest COG categories. KEGG pathway mapping was also carried out. The results showed that 24,651 unigenes were mapped to 128 predicted metabolic pathways (Figure 4; Table S4). The largest category was "Metabolism" including "Global map" (8,736), "Carbohydrate metabolism" (3,464), "Lipid metabolism" (2,955), "Amino acid metabolism" (1,390), "Biosynthesis of other secondary metabolites" (1,105), "Nucleotide metabolism" (1,033), "Energy metabolism" (936), "Metabolism of terpenoids and polyketidesgly" (818), "Glycan biosynthesis and metabolism" (619), "Metabolism of cofactors and vitamins" (586), and "Metabolism of other amino acids" (521) (Figure 4, A class; Table S4).

Differential Expression Analysis of Assembled L. aurea Transcripts under MeJA Treatment
To identify transcriptional responses of unigenes under MeJA treatment, reads from Con and MJ100 samples, were mapped to the obtained non-redundant unigenes. According to the FPKM values, the mappable reads were used to estimate the transcription levels. More than 97.0% of unigenes had FPKM values in the range of 1-100 ( Figure 5A). Then the expression levels of unigenes in both samples were calculated. The unigenes which had at least a two-fold change with FDR ≤ 0.001 ( Figure 5B) were screened and taken as differentially expressed genes (DEGs). Totally, we identified 4,165 DEGs between Con and MJ100 samples. Among them, 2,281 DEGs were found up-regulated and 1,884 down-regulated ( Figure 5C).

Functional Classification of the DEGs
To further elaborate the functions of DEGs, we performed GO enrichment analysis, using Fisher's exact test with an FDR adjusted P ≤ 0.01 as the cutoff. Of the 4,165 DEGs, 1,806 were assigned GO annotations (Table S5). For example, In the BP category, "cellular process, " "metabolic process, " "singleorganism process, " and "response to stimulus" were the topfour DEGs group (56.14%). In the CC category, the DEGs were annotated to "cell, " "cell part, " and "organelle" comprised the largest proportion (70.28%). Moreover, in the MF category, the genes which were associated with "catalytic activity" and "binding" took the biggest part (84.72%) of the DEGs (Table S6). Additionally, these DEGs were similarly enriched in the BP      categories, such as "metabolic process, " "response to stimulus, " and "cellular process" (Table S6). KEGG pathway enrichment analysis for DEGs was also performed. Of the 4,165 DEGs, 1,547 genes were assigned a KEGG ID and categorized into 121 pathways (Table S5). Of these, 29 pathways were significantly overrepresented under MeJA treatment, containing "Biosynthesis of secondary metabolites, " "Glycerophospholipid metabolism, " "Metabolic pathway, " "Plant hormone signal transduction, " and "Plant-pathogen interaction" (Table S7).

Verification of RNA-Seq Data by qRT-PCR
To test the reliability of the RNA-Seq data, the expression level of 15 transporter genes were selected for qRT-PCR assays (Figure 9). These candidates included 10 ABC transporters and 5 drug transmembrane transporters (Table S1; Figure 8). The RNA-Seq data were compared with the transcript abundance patterns of the MeJA treatment and control. Our results showed that almost all of the expression comparisons of qRT-PCR assay were in fairly good match with the RNA-Seq data, even if the fold-change of some genes in their expression level detected by sequencing and qRT-PCR did not match perfectly. These data confirmed the reliability of the RNA-Seq results (Figure 9). Only one ABC transporter gene (CL8796.Contig3), which was analyzed by qRT-PCR, revealed significant difference in expression level comparing with the RNA-Seq data (Figure 9). In a word, the expression patterns of all unigenes consistent with the RNA-Seq data, indicating that our experimental results were reliable.

DISCUSSION
The RNA-Seq should lead to the transcriptional profiling analysis. By using Illumina HiSeq TM 2,000 sequencing, we got two libraries (Con and MJ100), producing 4.82 Gbyte and 4.65 Gbyte of clean data, respectively, which are larger than the previous 454 GS platform database . The transcriptome assembly was carried out by using Trinity software, which has good performing in none reference genome assembling (Grabherr et al., 2011). It showed that the short reads with average lengths of 244 and 273 bp of the two libraries, were collected into 181,881 and 156,621 contigs, respectively (Figure 1). Sequencing profile is an effective tool to obtain functional genes . For example, inspection of the leaf epidermis enriched transcription database, an ABC transporter CrTPT2 that mediates the transporter of anticancer drug components in Catharanthus roseus was identified (Yu and Luca, 2013). Even more, combined the danshen (Salvia miltiorrhiza) EST database with the next-generation sequencing of mRNA from induced hairy root, Guo et al. (2013) identified CYP76AH1 which catalyzes turnover of multiradiene in tanshinones biosynthesis, leading to a successful heterologous production of ferruginol in yeast cell. In this study, combining with the previous EST database of L aurea , the entire transcriptome information obtained would be helpful for the future functional genomic research in L. aurea.
It has become clear that as an elicitor, MeJA is the main signal of secondary metabolite production across the plant kingdom, from angiosperms to gymnosperm (De Geyter et al., 2012). MeJA treatment triggers the majority of secondary metabolites biosynthesis (i.e., terpenoids, phenylpropanoids, and alkaloids) through an extensive transcriptional reprogramming (Zhao et al., 2005;Pauwels et al., 2009;De Geyter et al., 2012;Misra et al., 2014). Previously, we also showed that exogenous MeJA application accelerated the Amaryllidaceae alkaloids accumulation in L. chinensis seedlings (Mu et al., 2009). Here, by using trancriptome sequencing, we investigated the MeJA-responsive transcriptional changes, and identified 4,165 DEGs (including 2,281 up-regulated and 1,884 downregulated) between Con and MJ100 samples in L. aurea ( Figure 5C). They were categorized into 121 KEGG pathways (Table S5), and involved in 12 biosynthesis pathways of secondary metabolites (Figure 6). In general, Amaryllidaceae are regarded as derivatives of the common precursor 4 ′ -O-methylnorbelladine via intramolecular oxidative phenolcoupling (Eichhorn et al., 1998;Park, 2014), which belongs to the isoquinoline alkaloid biosynthesis pathway. Recently, a norbelladine 4 ′ -O-methyltrasferase (NpN4OMT) gene involved FIGURE 9 | Validation and comparative relative expression of 15 selected unigenes between the control and MeJA libraries in L. aurea. Differential expression between control (black column) and MeJA-treated sample (gray column) was compared. The normalized expression level (FPKM) of RNA-Seq is indicated on the x-axis. The relative expression levels were determined by quantitative real-time PCR using Cyclophilin 1 as an internal reference; each PCR was repeated three times and the error bars represent standard deviation (SD).
in the biosynthesis of galanthamine in Narcissus ap. aff. pseudonarcissus has been characterized (Kilgore et al., 2014). When aligned against transcriptome database, the potential homologous gene of NpN4OMT in L. aurea was also found, and its expression level showed more than two times higher under MeJA treatment, when compared with the control sample (Table S7). It suggested that DEGs database reserving a large number of genetic information which would be helpful for characterizing the key enzymes associated with the Amaryllidaceae alkaloids biosynthesis.
The expression of plant secondary pathway is under the tight control of a large number of TFs at different levels . For example, TFs involved in Jasmonates (JAs) signaling cascades usually regulate the transcription of multiple genes in a biosynthesis pathway, so as to improve the production of secondary metabolites (Zhou and Memelink, 2016). Several types of TFs shown as regulators of secondary metabolite biosynthesis in plants have been identified, belonging to the families AP2/ERF, bHLH, MYB, and WRKY (Naoumkina et al., 2008;Shoji et al., 2010;Todd et al., 2010;Yang et al., 2012;Yu et al., 2012;Misra et al., 2014). In this study, at least 1,591 unigenes encoding transcription factors (TFs) were found. Among them, 147 DEGs of TFs were simultaneously detected (Tables S5, S8). They were divided into 17 subfamilies, and most of them were extensively up-regulated in response to MeJA treatment (Figure 7; Table 3), suggesting their possible involvement in the regulation of secondary metabolite biosynthesis in L. aurea. The WRKY TF family is unique to plants, and is characterized by a conserved WRKY domain which specifically binds to the W-box sequence (Rushton et al., 2010;Zhou and Memelink, 2016). Previous reports indicated that many WRKY TFs may be regulated by wound and JA signal, and are reported as regulators of genes involved in various secondary metabolite pathways Phukan et al., 2016). As shown in this study, among the DEGs, WRKY cluster is the largest family and most of them were up-regulated under MeJA treatment (Table 3; Figure 7). This result suggested the potentially important role of WRKY TFs in the regulation of gene expression in L. aurea. In addition, the bHLH TF MYC2 has been reported as a master regulator in the JAs signaling network (Kazan and Manners, 2013). Transcription levels of the key enzymes of the primary and secondary metabolism were also regulated by MYC2 (Dombrecht et al., 2007). In this study, we also noticed that most of the bHLH TF genes (including 3 MYC genes), were up-regulated under the MeJA treatment in L. aurea (Figure 7). The expression changes of these transcription factors may reveal their key functions.
The verity of secondary metabolites has a great deal of differences among the species; the specific compounds are often restricted in a few species, or even within a few varieties within a species (Smetanska, 2008). Moreover, the biosynthesis and storage of these compounds are usually tissueand developmental stage-specific. In plants, the secondary metabolites always have to be stored in the special tissues or subcellular compartment which is distinct from the biosynthesis part. The mechanism of the biosynthesis and accumulation of plant secondary metabolites in such appropriate pattern attracted more attention (Yazaki et al., 2008;Nour-Eldin and Halkier, 2013). In most cases, plant secondary metabolites are transported intercellularly, intracellularly, and in an intratissue fashion by specific transporters (Yazaki, 2005;Yazaki et al., 2008;Nour-Eldin and Halkier, 2013). Meanwhile, the membrane transporter is comparatively specific and highly controlled for each secondary metabolite (Yazaki, 2005;Nour-Eldin and Halkier, 2013;Lv et al., 2016). ABC transporter family, based on hydrolysis of ATP, is proved to take a main part of transporting the secondary metabolites (Yazaki, 2005(Yazaki, , 2006. Since many plant secondary metabolites are medicinally used, the ABC transporters are associated with the drug resistant (Fletcher et al., 2010). In addition, the ABC transporters that presumably take part in secondary metabolite transport were also characterized among the MeJA up-regulated transcripts. For example, the ABCG transporter unigenes related to artemisinin yield in Artemisia annua, which are responsive to methyl jasmonate treatment, have been identified . In Panax ginseng, a novel PDR-type ABC transporter gene PgPDR3 induced by MeJA was also characterized (Zhang et al., 2013). It has been suggested that ABC transporters are often associated with the special compound transport, including alkaloids, terpenoids, and polyphenols, etc. (Sakai et al., 2002;Yazaki, 2005Yazaki, , 2006Yazaki et al., 2008;Shoji, 2014;Lv et al., 2016). Among the DEGs of L. aurea treated with MeJA, a large proportion of ABC transporters were identified (Figure 8; Table S8). Thus, the results will be helpful to identify the potential ABC transporters for translocation of secondary metabolites in L. aurea.
In this study, a large-scale unigene investigation of L. aurea under MeJA treatment was performed by Illumina sequencing. In total, we found 4,175 DEGs, and many transcripts were encoded by putative genes including transporter proteins, transcription factors and enzymes involved in secondary metabolism pathway. The data we obtained provided comprehensive information on gene discovery, transcriptome profiling, and transcriptional regulation of L. aurea. Additionally, our findings highlight the significance of JA signaling and Amaryllidaceae alkaloids synthesis in L. aurea, and provide a foundation for subsequent genomic research in the future.

AUTHOR CONTRIBUTIONS
ReW designed the project, performed the experiments and wrote the manuscript. RoW, SX, and NW analyzed the data and wrote the manuscript. BX and YJ provided the plants and revised the manuscript. All authors read and approved the final manuscript.