Rong Wang
Sheng Xu
Ning Wang
Bing Xia
Yumei Jiang
Ren Wang- 1Institute of Botany, Jiangsu Province and Chinese Academy of Sciences, Nanjing, China
- 2The Jiangsu Provincial Platform for Conservation and Utilization of Agricultural Germplasm, Nanjing, China
- 3Key Laboratory of Biology and Genetic Improvement of Soybean, National Center for Soybean Improvement, Ministry of Agriculture, Nanjing Agricultural University, Nanjing, China
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, anti-malarial, 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 Amaryllidaceae-type alkaloids (Wang et al., 2013). On the other hand, the inspection of the comparative transcriptome files between two databases provides another method to investigate the 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 CO2 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.
Materials and Methods
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 (Yu et al., 2016) by using NEBNext® 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 T4 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™ 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 log2-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™ 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™ 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).
Results
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).
Table 1. Summary of Illumina HiSeq™ 2,000 assembly and analysis of L. aurea transcriptome sequences.
Figure 1. Length distribution of assembled unigenes.
Sequence Annotation and Classification
All of the unigenes were annotated by BLAST search in the public databases (Table 2). The results revealed that 40,636 unigenes (68.13%) had significantly matched in the Nr database, 30,216 (50.66%) in the Nt database, and 26,345 (44.17%) in the Swiss-Prot database. Taken the entire public databases together, there were 41,585 unigenes (69.72%) successfully annotated (Table 2). For GO analysis, there were 31,157 unigenes divided into three ontologies, the percentage and summary 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).
Table 2. Functional annotation of non-redundant unigenes against the public databases.
Figure 2. Gene ontology (GO) categorization of assembled unigenes. The unigenes were categorized based on gene ontology annotation, and the proportion of each category is displayed based on the ontologies biological process, cellular component, and molecular function.
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).
Figure 3. Clusters of orthologous groups (COG) classifications of putative proteins. All putative proteins were aligned to the COG database and can be classified functionally into at least 25 molecular families.
Figure 4. KEGG metabolism pathway categories of assembled unigenes. Main functional categories are the Metabolism (A), Genetic Information Processing (B), Environmental Information Processing (C), Cellular Processes (D), and Organismal Systems (E). Bars represent the numbers of L. aurea assignments of unigenes with BLASTX matches to each KEGG term.
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).
Figure 5. Transcriptomes of L. aurea under MeJA treatment. (A) Number of unigenes expressed in each sample. (B) Scatter-plot graphs of the differential gene expression patterns between Con and MJ100 libraries. DEGs were determined using a threshold of FDR ≤ 0.001 and |log2 Ratio| ≥ 1. Red spots represent up-regulated DEGs and green spots indicate down-regulated DEGs. Those shown in blue are unigenes that did not change significantly under MeJA treatment. (C) Number of differentially expressed genes (DEGs) showing up- (red) or down- (green) regulation between the samples.
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,” “single-organism process,” and “response to stimulus” were the top-four 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).
Identification of MeJA-Responsive Genes Expressed during MeJA Treatment
We analyzed 209 DEGs involving in 12 biosynthesis pathways of secondary metabolites (Figure 6). The results showed that the largest subcategory was phenylpropanoid biosynthesis (58), followed by flavonoid biosynthesis (36), flavone and flavonol biosynthesis (35), and stilbenoid, diarylheptanoid and gingerol biosynthesis (35). To identify MeJA-related genes in L. aurea, we used unigene sequences in BLAST searches of the public databases, and found 1,591 and 2,111 unigenes encoding transcription factors (TFs) and transporter proteins (TPs). Among them, 147 and 138 DEGs of TFs and TPs were detected respectively (Tables S5, S8). In total, the DEGs of TFs were divided into 17 subfamilies including WRKY, APETALA2/Ethylene-Response Factors (AP2/ERF), basic Helix-Loop-Helix (MYBs), bZIPs, and so on. Most of them were extensively up-regulated in response to MeJA treatment in L. aurea (Figure 7; Table 3). For example, 26 out of 32 WRKY TFs and all the 14 MYB TFs are up-regulated under MeJA treatment (Table 3). Meanwhile, numerous genes encoding transporters were also included in the sets of DEGs we detected (Table S8). The results showed that, under MeJA treatment, the expression level of 138 transporter genes varied (Figure 8). Among them, zinc transporter is the largest cluster family (21, 15.22%), and most of them were down-regulated under MeJA treatment. Notably, the expression level of genes encoding ATP-Binding Cassette (ABC) transporters (20, 14.49%), ammonium transporters (3, 2.17%), amino acid/peptide/protein transporters (23, 16.67%), drug transporters (11, 7.97%), electron transporters (6, 4.35%), iron ion transporters (5, 3.62%), magnesium transporters (3, 2.17%), nucleobase-ascorbate transporters (3, 2.17%), nucleoside transporters (2, 1.45%), proton transporters (5, 3.62%), sugar transporters (10, 7.25%), sulfate transporters (4, 2.90%), and uncharacterized transporters (7, 5.07%) changed under MeJA treatment (Figure 8).
Figure 6. Expression profiling of DEGs in secondary metabolism pathway. (A) List of DEGs involved in different secondary metabolism pathway under MeJA treatment. (B) Distribution of DEGs involved in different secondary metabolism pathway under MeJA treatment.
Figure 7. Expression profiling of DEGs annotated as transcription factors.
Table 3. Regulated transcription factors (TFs) under the MeJA treatment.
Figure 8. Expression profiling of DEGs annotated as transporters.
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.
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).
Discussion
The RNA-Seq should lead to the transcriptional profiling analysis. By using Illumina HiSeq™ 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 (Wang et al., 2013). 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 (Wang et al., 2009). 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 (Wang et al., 2013), 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 down-regulated) 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 phenol-coupling (Eichhorn et al., 1998; Park, 2014), which belongs to the isoquinoline alkaloid biosynthesis pathway. Recently, a norbelladine 4′-O-methyltrasferase (NpN4OMT) gene involved 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 (Yang et al., 2012). 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 (Chen et al., 2012; 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 tissue- and 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, 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 (Zhang et al., 2012). 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, 2005, 2006; Yazaki 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.
Conflict of Interest Statement
The authors declare that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest.
The reviewers ZM, JZ and handling Editor declared their shared affiliation, and the handling Editor states that the process nevertheless met the standards of a fair and objective review.
Acknowledgments
This work was supported by the National Natural Science Foundation of China (31270339, 31301798, 31572151), the Natural Science Foundation of the Jiangsu Province (BK20161380), the “333 project” of Jiangsu Province (BRA2015432), and the Jiangsu Provincial Public Institutions Program for Research Conditions and Building Capacity (BM2015019).
Supplementary Material
The Supplementary Material for this article can be found online at: http://journal.frontiersin.org/article/10.3389/fpls.2016.01971/full#supplementary-material
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.
References
Berenbaum, M. (1995). Phototoxicity of plant secondary metabolites: insect andmammalian perspectives. Arch. Insect. Biochem. Physiol. 29, 119–134. doi: 10.1002/arch.940290204
Bores, G. M., Huger, F. P., Petko, W., Mutlib, A. E., Camacho, F., Rush, D. K., et al. (1996). Pharmacological evaluation of novel Alzheimer's disease therapeutics: acetylcholinesterase inhibitors related to galanthamine. J. Pharmacol. Exp. Ther. 277, 728–738.
Chen, L., Song, Y., Li, S., Zhang, L., Zou, C., and Yu, D. (2012). The role of WRKY transcription factors in plant abiotic stresses. Biochim. Biophys. Acta 1819, 120–128. doi: 10.1016/j.bbagrm.2011.09.002
De Geyter, N., Gholami, A., Goormachtig, S., and Goossens, A. (2012). Transcriptional machineries in jasmonate-elicited plant secondary metabolism. Trends Plant Sci. 17, 349–359. doi: 10.1016/j.tplants.2012.03.001
Devi, K., Mishra, S. K., Sahu, J., Panda, D., Modi, M. K., and Sen, P. (2016). Genome wide transcriptome profiling reveals differential gene expression in secondary metabolite pathway of Cymbopogon winterianus. Sci. Rep. 6:21026. doi: 10.1038/srep21026
Dixon, R. A. (2001). Natural products and plant disease resistance. Nature 411, 843–847. doi: 10.1038/35081178
Dombrecht, B., Xue, G. P., Sprague, S. J., Kirkegaard, J. A., Ross, J. J., Reid, J. B., et al. (2007). MYC2 differentially modulates diverse jasmonate-dependent functions in Arabidopsis. Plant Cell 19, 2225–2245. doi: 10.1105/tpc.106.048017
Eichhorn, J., Takada, T., Kita, Y., and Zenk, M. H. (1998). Biosynthesis of the Amaryllidaceae alkaloid galanthamine. Phytochemistry 49, 1037–1047. doi: 10.1016/S0031-9422(97)01024-8
Fletcher, J. I., Haber, M., Henderson, M. J., and Norris, M. D. (2010). ABC transporters in cancer: more than just drug efflux pumps. Nat. Rev. Cancer 10, 147–156. doi: 10.1038/nrc2789
Gandhi, S. G., Mahajan, V., and Bedi, Y. S. (2015). Changing trends in biotechnology of secondary metabolism in medicinal and aromatic plants. Planta 241, 303–317. doi: 10.1007/s00425-014-2232-x
Goossens, A., Häkkinen, S. T., Laakso, I., Seppänen-Laakso, T., Biondi, S., De, S. V., et al. (2003). A functional genomics approach toward the understanding of secondary metabolism in plant cells. Proc. Natl. Acad. Sci. U.S.A 100, 8595–8600. doi: 10.1073/pnas.1032967100
Grabherr, M. G., Haas, B. J., Yassour, M., Levin, J. Z., Thompson, D. A., Amit, I., et al. (2011). Full-length transcriptome assembly from RNA-Seq data without a reference genome. Nat. Biotechnol. 29, 644–652. doi: 10.1038/nbt.1883
Guo, J., Zhou, Y. J., Hillwig, M. L., Shen, Y., Yang, L., Wang, Y. J., et al. (2013). CYP76AH1 catalyzes turnover of miltiradiene in tanshinones biosynthesis and enables heterologous production of ferruginol in yeasts. Proc. Natl. Acad. Sci. U.S.A 110, 12108–12113. doi: 10.1073/pnas.1218061110
Harvey, A. L. (1995). The pharmacology of galanthamine and its analogues. Pharmacol. Ther. 68, 113–128. doi: 10.1016/0163-7258(95)02002-0
Howe, E. A., Sinha, R., Schlauch, D., and Quackenbush, J. (2011). RNA-Seq analysis in MeV. Bioinformatics 27, 3209–3210. doi: 10.1093/bioinformatics/btr490
Huang, Q., Huang, X., Deng, J., Liu, H., Liu, Y., Yu, K., et al. (2016). Differential gene expression between leaf and rhizome in Atractylodes lancea: a comparative transcriptome analysis. Front. Plant Sci. 7:348. doi: 10.3389/fpls.2016.00348
Hussain, M., Sarfaraj, F. S., Ansari, S., Rahman, M., Akhlaquer, A. I. Z., and Saeed, M. (2012). Current approaches toward production of secondary plant metabolites. J. Pharm. Bioallied. Sci. 4, 10–20. doi: 10.4103/0975-7406.92725
Jimenez-Garcia, S. N., Vazquez-Cruz, M. A., Guevara-Gonzalez, R. G., Torres-Pacheco, I., Cruz-Hernandez, A., and Feregrino-Perez, A. A. (2013). Current approaches for enhanced expression of secondary metabolites as bioactive compounds in plants for agronomic and human health purposes–a review. Pol. J. Food Nutr. Sci. 63, 67–78. doi: 10.2478/v10222-012-0072-6
Jin, Z. (2009). Amaryllidaceae and Sceletium alkaloids. Nat. Prod. Rep. 26, 363–381. doi: 10.1039/b718044f
Kalra, S., Lal Puniya, B., Kulshreshtha, D., Kumar, S., Kaur, J., Ramachandran, S., et al. (2013). De novo transcriptome sequencing reveals important molecular networks and metabolic pathways of the plant, Chlorophytum borivilianum. PLoS ONE 8:e83336. doi: 10.1371/journal.pone.0083336
Kazan, K., and Manners, J. M. (2013). MYC2: the master in action. Mol. Plant 6, 686–703. doi: 10.1093/mp/sss128
Kennedy, D. O., and Wightman, E. L. (2011). Herbal extracts and phytochemicals: plant secondary metabolites and the enhancement of human brain function. Adv. Nutr. 2, 32–50. doi: 10.3945/an.110.000117
Kilgore, M. B., Augustin, M. M., Starks, C. M., O'Neil-Johnson, M., May, G. D., Crow, G. A., et al. (2014). Cloning and characterization of a Norbelladine 4′-O-methyltransgerase involved in the biosynthesis of the Alzheimer's drug galanthamine in Narcossis sp. aff. Pseudonarcissus. PLoS ONE 9:e103223. doi: 10.1371/journal.pone.0103223
Kliebenstein, D. J., and Osbourn, A. (2012). Making new molecules-evolution of pathways for novel metabolites in plants. Curr. Opin. Plant Biol. 15, 415–423. doi: 10.1016/j.pbi.2012.05.005
Li, Z., Cheng, Y., Cui, J., Zhang, P., Zhao, H., and Hu, S. (2015). Comparative transcriptome analysis reveals carbohydrate and lipid metabolism blocks in Brassica napus L. male sterility induced by the chemical hybridization agent monosulfuron ester sodium. BMC Genomics 16:206. doi: 10.1186/s12864-015-1388-5
Livak, L. J., and Schmittgen, T. D. (2001). Analysis of relative gene expression data using real-time quantitative PCR and the 2-ΔΔCT. Methods 25, 402–408. doi: 10.1006/meth.2001.1262
Lv, H. J., Li, J. J., Wu, Y. Y., Garyali, S., and Wang, Y. (2016). Transporter and its engineering for secondary metabolites. Appl. Microbiol. Biotechnol. 100, 6119–6130. doi: 10.1007/s00253-016-7605-6
Ma, R., Xu, S., Zhao, Y. C., Xia, B., and Wang, R. (2016). Selection and validation of appropriate reference genes for quantitative real-time PCR analysis of gene Expression in Lycoris aurea. Front. Plant Sci. 7:536. doi: 10.3389/fpls.2016.00536
Misra, R. C., Maiti, P., Chanotiya, C. S., Shanker, K., and Ghosh, S. (2014). Methyl jasmonate-elicited transcriptional responses and pentacyclic triterpene biosynthesis in sweet basil. Plant Physiol. 164, 1028–1044. doi: 10.1104/pp.113.232884
Mortazavi, A., Williams, B. A., McCue, K., Schaeffer, L., and Wold, B. (2008). Mapping and quantifying mammalian transcriptomes by RNA-Seq. Nat. Methods 5, 621–628. doi: 10.1038/nmeth.1226
Mu, H. M., Wang, R., Li, X. D., Jiang, Y. M., Wang, C. Y., Quan, J. P., et al. (2009). Effect of abiotic and biotic elicitors on growth and alkaloid accumulation of Lycoris chinensis seedlings. Z. Naturforsch. C 64, 541–550. doi: 10.1515/znc-2009-7-813
Naoumkina, M. A., He, X. Z., and Dixon, R. A. (2008). Elicitor-induced transcription factors for metabolic reprogramming of secondary metabolism in Medicago truncatula. BMC Plant Biol. 8:132. doi: 10.1186/1471-2229-8-132
Nascimento, N. C., and Fett-Neto, A. G. (2010). Plant secondary metabolism and challenges in modifying its operation: an overview. Methods Mol. Biol. 643, 1–13. doi: 10.1007/978-1-60761-723-5_1
Nour-Eldin, H. H., and Halkier, B. A. (2013). The emerging field of transport engineering of plant specialized metabolites. Curr. Opin. Biotechnol. 24, 263–270. doi: 10.1016/j.copbio.2012.09.006
Park, J. B. (2014). Synthesis and characterization of norbelladine, a precursor of Amaryllidaceae alkaloid, as an anti-inflammatory/anti-COX compound. BMC Lett. 24, 5381–5384. doi: 10.1016/j.bmcl.2014.10.051
Pauwels, L., Inzé, D., and Goossens, A. (2009). Jasmonate-inducible gene: what does it mean? Trends Plant Sci. 14, 87–91. doi: 10.1016/j.tplants.2008.11.005
Pertea, G., Huang, X., Liang, F., Antonescu, V., Sultana, R., Karamycheva, S., et al. (2003). TIGR Gene Indices clustering tools (TGICL): a software system for fast clustering of large EST datasets. Bioinformatics 19, 651–652. doi: 10.1093/bioinformatics/btg034
Phukan, U. J., Jeena, G. S., and Shukla, R. K. (2016). WRKY transcription factors: molecular regulation and stress responses in plants. Front. Plant Sci. 7:760. doi: 10.3389/fpls.2016.00760
Reiner, A., Yekutieli, D., and Benjamini, Y. (2003). Identifying differentially expressed genes using false discovery rate controlling procedures. Bioinformatics 19, 368–375. doi: 10.1093/bioinformatics/btf877
Rushton, P. J., Somssich, I. E., Ringler, P., and Shen, Q. J. (2010). WRKY transcription factors. Trends Plant Sci. 15, 247–258. doi: 10.1016/j.tplants.2010.02.006
Sakai, K., Shitan, N., Sato, F., Ueda, K., and Yazaki, K. (2002). Characterization of berberine transport into Coptis japonica cells and the involvement of ABC protein. J. Exp. Bot. 53, 1879–1886. doi: 10.1093/jxb/erf052
Shoji, T. (2014). ATP-binding cassette and multidrug and toxic compound extrusion transporters in plants: a common theme among diverse detoxification mechanisms. Int. Rev. Cell Mol. Biol. 309, 303–346. doi: 10.1016/B978-0-12-800255-1.00006-5
Shoji, T., Kajikawa, M., and Hashimoto, T. (2010). Clustered transcription factor genes regulate nicotine biosynthesis in tobacco. Plant Cell 22, 3390–3409. doi: 10.1105/tpc.110.078543
Smetanska, I. (2008). Production of secondary metabolites using plant cell cultures. Adv. Biochem. Eng. Biotechnol. 111, 187–228. doi: 10.1007/10_2008_103
Todd, A. T., Liu, E., Polvi, S. L., Pammett, R. T., and Page, J. E. (2010). A functional genomics screen identifies diverse transcription factors that regulate alkaloid biosynthesis in Nicotiana benthamiana. Plant J. 62, 589–600. doi: 10.1111/j.1365-313X.2010.04186.x
Vasconsuelo, A., and Boland, R. (2007). Molecular aspects of the early stages of elicitation of secondary metabolites in plants. Plant Sci. 172, 861–875. doi: 10.1016/j.plantsci.2007.01.006
Verma, N., and Shukla, S. (2015). Impact of various factors responsible for fluctuation in plant secondary metabolites. J. Appl. Res. Med. Arom. Plant 2, 105–113. doi: 10.1016/j.jarmap.2015.09.002
Wang, R., Xu, S., Jiang, Y., Jiang, J., Li, X., Liang, L., et al. (2013). De novo sequence assembly and characterization of Lycoris aurea transcriptome using GS FLX titanium platform of 454 pyrosequencing. PLoS ONE 8:e60449. doi: 10.1371/journal.pone.0060449
Wang, Z., Gerstein, M., and Snyder, M. (2009). RNA-Seq: a revolutionary tool for transcriptomics. Nat. Rev. Genet. 10, 57–63. doi: 10.1038/nrg2484
Wei, H., Chen, X., Zong, X., Shu, X., Gao, D., and Liu, Q. (2015). Comparative transcriptome analysis of genes involved in anthocyanin biosynthesis in the red and yellow fruits of sweet cherry (Prunus avium L.). PLoS ONE 10:e0121164. doi: 10.1371/journal.pone.0121164
Yang, C. Q., Fang, X., Wu, X. M., Mao, Y. B., Wang, L. J., and Chen, X. Y. (2012). Transcriptional regulation of plant secondary metabolism. J. Integr. Plant Biol. 54, 703–712. doi: 10.1111/j.1744-7909.2012.01161.x
Yazaki, K. (2005). Transporters of secondary metabolites. Curr. Opin. Plant Biol. 8, 301–307. doi: 10.1016/j.pbi.2005.03.011
Yazaki, K. (2006). ABC transporters involved in the transport of plant secondary metabolites. FEBS Lett. 580, 1183–1191. doi: 10.1016/j.febslet.2005.12.009
Yazaki, K., Sugiyama, A., Morita, M., and Shitan, N. (2008). Secondary transport as an efficient membrane transport mechanism for plant secondary metabolites. Phytochem. Rev. 7, 513–524. doi: 10.1007/s11101-007-9079-8
Yu, C. G., Xu, S., and Yin, Y. L. (2016). Transcriptome analysis of the Taxodium ‘Zhongshanshan 405’ roots in response to salinity stress. Phytochem. Rev. 100, 439–444. doi: 10.1016/j.plaphy.2016.01.009
Yu, F., and Luca, V. D. (2013). ATP-binding cassette transporter controls leaf surface secretion of anticancer drug components in Catharanthus roseus. Proc. Natl. Acad. Sci. U.S.A. 110, 15830–15835. doi: 10.1073/pnas.1307504110
Yu, Z. X., Li, J. X., Yang, C. Q., Hu, W. L., Wang, L. J., and Chen, X. Y. (2012). The jasmonate-responsive AP2/ERF transcription factors AaERF1 and AaERF2 positively regulate artemisinin biosynthesis in Artemisia annua L. Mol. Plant 5, 353–365. doi: 10.1093/mp/ssr087
Zhang, L., Lu, X., Shen, Q., Chen, Y., Wang, T., Zhang, F., et al. (2012). Identification of Putative Artemisia annua ABCG transporter unigenes related to artemisinin yield following expression analysis in different plant tissues and in response to methyl jasmonate and abscisic acid treatments. Plant Mol. Biol. Rep. 30, 838–847. doi: 10.1007/s11105-011-0400-8
Zhang, R., Huang, J., Zhu, J., Xie, X., Tang, Q., Chen, X., et al. (2013). Isolation and characterization of a novel PDR-type ABC transporter gene PgPDR3 from Panax ginseng C.A. Meyer induced by methyl jasmonate. Mol. Biol. Rep. 40, 6195–6204. doi: 10.1007/s11033-013-2731-z
Zhao, J., Davis, L. C., and Verpoorte, R. (2005). Elicitor signal transduction leading to production of plant secondary metabolites. Biotechnol. Adv. 23, 283–333. doi: 10.1016/j.biotechadv.2005.01.003
Keywords: Amaryllidaceae alkaloids, Lycoris aurea, methyl jasmonate, transcriptome sequencing, secondary metabolite
Citation: Wang R, Xu S, Wang N, Xia B, Jiang Y and Wang R (2017) Transcriptome Analysis of Secondary Metabolism Pathway, Transcription Factors, and Transporters in Response to Methyl Jasmonate in Lycoris aurea. Front. Plant Sci. 7:1971. doi: 10.3389/fpls.2016.01971
Received: 01 October 2016; Accepted: 12 December 2016;
Published: 05 January 2017.
Edited by:
Ping Ma, University of Georgia, USACopyright © 2017 Wang, Xu, Wang, Xia, Jiang and Wang. This is an open-access article distributed under the terms of the Creative Commons Attribution License (CC BY). The use, distribution or reproduction in other forums is permitted, provided the original author(s) or licensor are credited and that the original publication in this journal is cited, in accordance with accepted academic practice. No use, distribution or reproduction is permitted which does not comply with these terms.
*Correspondence: Ren Wang, jswangren@aliyun.com
†These authors have contributed equally to this work.