De novo transcriptome analysis using 454 pyrosequencing of the Himalayan Mayapple, Podophyllum hexandrum

454 sequencing of the Mayapple cell culture transcriptome

Comparison between default and optimized parameters of Newbler

Here we present a simple workflow for 454 sequencing, assembly, annotation and other analyses (Figure 1). Newbler is frequently used in de novo pyrosequencing projects[23]. The comparative denovo assembly was carried out using Newbler with default and optimized parameters[24]. The optimized parameter generated 40,380 assembled sequences, comprising 12,940 contigs and 27,440 singlets with an N50 of 463 and 240 for contigs and singlets respectively. Newbler with optimized parameters gave the best results in terms of the numbers of assembled contigs and singlet, N50, mean contig/singlets length and total bases of contigs /singlets (bp) (Table 2). Additional file1 (A and B) shows the distribution of contig lengths generated by Newbler using default and optimized parameters respectively. Further analysis of the singlet generated by Newbler default assembly were excluded because their mean length was below 200 bp.

Functional annotation of assembled transcripts and determination of FPKM values

The annotation of transcripts was carried out using green plants of non-redundant (nr) protein database at NCBI by BLASTX[25]. BLASTX resulted in the annotation of 3,249 contigs out of 3,372 assembled contigs obtained using Newbler default parameters (Additional file2) whereas 40,380 transcripts from among 12,940 contigs and 27,440 singlet generated using Newbler optimized parameters (Additional file3). Using default parameters, transcripts showed significant similarity with P.trichocarpa followed by Oryza sativa, Ricinus communis and so forth, (Figure 2A) while, using optimized parameters, significant similarity was achieved with M. truncatula, followed by Glycine Max, Sorghum bicolor, P. trichocarpa and others (Figure 2B). GO assignments[26] were used to classify the functions of the predicted transcript contigs to determine the E-value distribution, sequence similarity distribution, evidence code distribution of sequences, evidence code distribution of BLAST hits, annotation score distribution, annotation distribution and GO-level distribution of transcripts generated by default (Additional file4) and optimized Newbler assembly ( Additional file5). The GO annotation of transcripts from Newbler default and optimized assemblies are represented in Additional file6 and Additional file7 respectively. GO-level sequence distribution for Biological processes, Molecular functions and Cellular components of the transcript contigs and singlet generated by Newbler default and optimized assembly are shown in Figure 3 (A, B respectively). Digital expression profiling by FPKM of each transcript generated from Newbler using default and optimized parameters assembly were also determined (Additional file8 and Additional file9).

To identify the biological pathways that are active in P. hexandrum cell culture, we mapped the annotated sequences to the reference canonical pathways in KEGG[27]. Among the annotated sequences generated by Newbler using default and optimized parameters, 321 and 1069 unique non-redundant sequences were involved in a particular KEGG pathway of which 32 and 100 unique sequences could be assigned to secondary metabolism respectively (Additional file10 and Additional file11).

Protein domains encoded by the P. hexandrum transcriptome that may represent genes involved in podophyllotoxin biosynthesis

We were interested in probable podophyllotoxin pathway genes that could be identified from the transcriptome, therefore Conserved Domains Database (CDD), Pfam, and Tigr databases were searched for domains encoding CADs, monooxygenases, peroxidases (POD), pinoresinol reductases, DPOs, SDGs, and methyl transferases. Our search identified transcripts coding for domains of CAD, SDG, monooxygenase, POD, methyl transferase, NADB Rossmann superfamily, Flavin utilizing monooxygenases superfamily, Uroporphyrinogen decarboxylase methyltransferase (URO-D CIMS) like superfamily, Isoprene-C2-like reductase (ISOPREN C2) like superfamily, Cytochrome oxidase (CypX) superfamily and Oxidoreductase q1 (Oxidored q1) superfamily (Additional file12 and Additional file13). According to the hypothetical scheme of podophyllotoxin pathway[18], matairesinol is converted to podophyllotoxin by two consecutive methyl group additions forming a compound like yatein. We were also interested in finding methyl transferases that can transfer two methyl groups to the same substrate at the same time. A Uroporphyrinogen IIIC methyl transferase from P. hexandrum was identified in our previous studies, which is known to function in transferring two methyl groups from S-Adenosyl-L-methionine (SAM) to its substrate[28, 29]. Therefore, in addition to finding SAM dependent methyl transferases, we also identified transcripts encoding URO-D CIMS domains.

Combining KAAS-KEGG pathway analysis with domain searching for phenylpropanoid and probable downstream podophyllotoxin pathway genes

BLASTX analysis and KAAS-KEGG pathway mapping of transcripts from Newbler default and optimized parameter identified cDNAs encoding Phenylalanine ammonia lyase (PAL), hydroxyl cinnamoyl transferase (HCT), cinnamate-4-hydroxylase (C4H), 4-Coumarate Ligase (4CL), 4-Coumarate CoA Ligase (4CL), cinnamoyl reductase (CCR), CAD, sinapyl alcohol dehydrogenase (SAD), β-Glucosidase (β-GLUC) and POD as being involved in the phenylpropanoid pathway. Podophyllotoxin pathway is hypothesized to start with CAD converting coniferaldehyde to coniferyl alcohol, which then undergoes dirigent-mediated coupling to form pinoresinol. Specific reductases, dehydrogenases and methyl transferases are then believed to convert pinoresinol to podophyllotoxin. We surveyed the CDD results for cDNAs with domains that may represent genes of this pathway and identified transcripts containing Phenylcoumaran benzylic ether reductase (PCBER), SDGs, monooxygenases, SAM dependent methyl transferases and URO-D CIMS-like domains. A scheme has been presented combining the BLASTX annotation, KAAS-KEGG mapping and domain search for phenylpropanoid pathway transcripts and transcripts with domains that may be part of podophyllotoxin biosynthetic pathway ( Figure 4).

Transcription factors related to secondary metabolism

Controlled transcription of biosynthetic genes is an important mechanism for regulating secondary metabolite production in plant cells[30]. Certain TFs are known to be involved in the regulation of secondary metabolism, such as R2R3-MYB, basic helix-loop-helix (bHLH) proteins like CrMYC2, AP2/ERF family proteins, WRKY, NAC, DOF, bZIP, HD-ZIP, and TFIIIA zinc finger TFs. We identified 96 transcripts from Newbler default assembly (Additional file14) and 961 transcripts from optimized parameter assembly (Additional file15) that may encode TFs (Figure 5A, B respectively). Amongst them, notable transcripts were AP2-EREBP, NAC, bHLH, MYB or MYB related, bZIP, mTERF, WRKY, C2C2-CO-like and C2C2-Dof. A number of plant MYB TFs regulating the phenylpropanoid biosynthetic pathway, identified from many species, including Arabidopsis, apple, grape, maize, petunia and snapdragon, most of which are R2R3-MYB TFs[31] can be correlated with our study as 48 transcripts coding for MYB or MYB related TFs have been identified from the optimized Newbler assembly. Again, R2R3-MYB transcription factor MYB12 in A. thaliana has been shown to function as a flavonol specific activator of flavonoid biosynthesis[32]. Transcriptional regulation of flavonoid biosynthesis, a major branch of phenylpropanoid pathway, controlled by a set of R2R3 MYB transcription factors, have been reported in several plants such as Prunus persica, Epimedium sagittatum as well[33, 34]. Other than this TF, 18 transcripts coded for bHLH TFs have been identified here. The bHLH-domain of the maize R-gene is reported to participate in anthocyanin formation and serve as a link between flavonoid formation and histone modification[35]. Amongst the diverse functions, bHLH transcription factors also regulate the biosynthetic pathway of flavonoids, in several plant species[31]. 1 DOF family TF has been identified in our analyses. AtDOF4 is known to influence phenylpropanoid metabolism in an environmental and tissue-specific manner by positively regulating the production of hydroxycinnamic acids in the hypocotyl and flower buds, and negatively regulating flavonoid biosynthesis in pollen grains[36]. Together, TFs identified here and related to the phenylpropanoid pathway can be explored further in the regulation of podophyllotoxin biosynthesis.

In silico SSR marker identification

SSRs can be divided into genomic SSRs and EST-SSRs. EST-SSRs are more evolutionary conserved than noncoding sequences and therefore have a relatively high transferability[33, 37]. Next-generation sequencing has identified EST-SSRs in many plant species[38–41]. However, there have been no reports of EST-SSRs in P. hexandrum to date.

SSRs were identified with MISA search tool (http://pgrc.ipk-gatersleben.de/misa/), which is based on the criteria that a dinucleotide or a trinucleotide pattern should appear at least six times, and tetra, penta and hexa nucleotide patterns should appear five times each (Additional file16 for Newbler default, Additional file17 for Newbler optimized). SSR distribution and SSR mining of transcripts identified a total of 1,011 SSRs from 40,380 transcripts, with 94 transcripts containing more than one SSR. The most abundant repeat type was dinucleotides (68.6%) and the dominant tandem repeat motifs were (AT)6 and (AT)7 representing 19.4% and 25.7% respectively.

Transcriptome wide survey of miRNA targets in P. hexandrum cell cultures

MiRNAs are known to regulate many developmental and effector genes at the posttranscriptional level[38, 42]. Using oligonucleotide arrays, miRNAs have been shown to be differentially expressed between tissues and during the maturation of the fruit in the grapevine[43]. Wong et al.[44], predicted three wood related genes, flavonol synthase-like, xyloglucan fucosyltransferase and glucan synthase-like genes to be the targets of miR170, miR172 and miR319, respectively, and suggested that these miRNAs might be directly involved in regulation of the phenylpropanoid pathway and hemicellulose biosynthesis pathway. Downregulation of Flavonol synthase by miR170 would redirect the precursor 4-coumaroyl CoA to lignin biosynthesis.

We identified precursor and mature miRNAs in the P. hexandrum cell culture transcriptome, by searching the contigs and singlet in the public miRNA database (Additional file18 and Additional file19)[41, 45]. Transcripts targeted by miRNAs that are possibly related to phenylpropanoid and podophyllotoxin biosynthesis include cytochrome b6, cytochrome p450 like, cell wall associated hydrolase, cell wall associated protein, laccase, and cytochrome f. Deoxypodophyllotoxin 6-hydroxylase is a cytochrome p450 dependent monooxygenase, that catalyzes the introduction of a hydroxyl group in position 6 of deoxypodophyllotoxin[46]. A cytochrome p450 protein is known to catalyze the biosynthesis of a lignan, (+)- sesamin, by forming two methylenedioxy bridges[47]. Laccases have auxiliary roles in stereoselective coupling to 8-8′ linked lignans[48].

Comparative qRT-PCR of selected phenylpropanoid pathway genes in cell culture, callus and rhizome and podophyllotoxin accumulation in P. hexandrum cell culture

Podophyllotoxin content in P. hexandrum rhizome is known to vary from 4% to 10% in reference to age, altitude, net photosynthetic rate, stomatal conductance, carbon uptake and number of leaves[49–52]. Green calli were found to have 40–50 μg/g podophyllotoxin which were used to generate cell suspension culture. Enhanced accumulation of podophyllotoxin was observed from 3 day old cell culture and increased up to 12 day (Figure 6A) as noted by HPLC analysis (Figure 6B, C). However, no significantly increased accumulation was observed till 18 days. Hence, we chose to compare relative gene expression profile of selected phenylpropanoid pathway transcripts amongst the calli, 12 day old cell culture, and the rhizome. Transcripts of CAD identified here, CAD1, CAD5, and CAD8 share sequence similarity to Arabidopsis CAD1, CAD5, and CAD8 respectively as identified by BLASTX analysis. Primers were designed from transcript contigs for qRT-PCR analysis (Table 3). PAL in 12 days cell culture shows 12.12 fold upregulation with respect to the callus, while upregulation of PAL expression levels in rhizome is insignificant (Figure 6D). C4H is upregulated in the cell culture and the rhizome by 8 and 3.4 folds respectively in comparison to that of callus. The high FPKM value of C4H in the 12 day cell culture samples may correlate with this observation (Table 3). Furthermore, upregulation of CAD 1,5, and 8 can be correlated as well with the FPKM value on the higher side.

Sample preparation and 454 pyrosequencing

Calli were induced from mature leaves of P. hexandrum in MS medium[53] supplemented with 2.68 μM Napthyl acetic acid (NAA) and 8.88 μM Benzylaminopurine (BAP)[54]. Cell suspension cultures were initiated from freshly subcultured green calli of P. hexandrum in modified liquid MS medium[55] containing 60mM total N content, 1.25 mM potassium dihydrogen phosphate, 6% glucose and 11.41 μM Indole acetic acid (IAA). An inoculum of 5 g cells was used in 50 ml of cell suspension culture medium. Six flasks containing cell suspension cultures were shaken in the dark at 110 rpm for 12 days. The cells were collected by centrifugation at 1000 × g for 5 min, and immediately put in liquid nitrogen and used for RNA isolation. Total RNA was isolated from 12 days old cell suspension cultures using Purelink miRNA isolation kit (Invitrogen, CA, USA). Total RNA was quantified by NanoDrop technology, checked on a 1% denaturing agarose gel and on a bioanalyzer 2100 (Agilent Technologies, Palo Alto, CA). Removal of rRNA from total RNA was performed using RiboMinus Plant kit for RNA seq (Invitrogen cat. no. A10838-08) using the standard procedure and then concentrated by RiboMinus concentration module (cat no. K155005), according to the manufacturer’s instructions. Library preparation performed using a cDNA Rapid Library Preparation Method Manual-GS FLX Titanium Series, according to the manufacturer’s instructions. For transcriptome sequencing, 1 μg of Ribo-minus total RNA from each sample was used for fragmentation using ZnCl₂ solution, followed by ds-cDNA synthesis using a standard cDNA synthesis kit (Roche). This ds-cDNA was then subjected to fragment end-repair followed by adaptor ligation using Rapid Library Prep kit (Roche). emPCR amplification of the cDNA library was performed according to the manufacturer’s instructions (Roche). Clonally amplified cDNA library beads obtained from the emPCR amplification reaction were deposited on a PTP for sequencing using pyrosequencing chemistry. The next-generation sequencing run for whole transcriptome analysis was performed on a Roche 454 GS FLX.

De novo assembly

Raw reads obtained from 454 pyrosequencing were preprocessed by removing low quality reads, and adapter/ primer sequences using PRINSEQ. The high quality reads were uniqed and mapped to non-coding RNA database Rfam (v11.0) using gsMapper. Reads that mapped to ncRNAs sequences were excluded and remaining reads were used for further analysis.The preprocessed reads were then assembled using Newbler with default parameters and optimized parameters. Optimized parameters were set by checking “Use duplicate reads”, “Extend low depth overlaps”, “Reads limited to one contig” and “Single Ace file options”. The sequence data generated in this study have been deposited at NCBI in the Short Read Archive database under the accession (SRX180870 and SRX180386).

Functional annotation, GO mapping, pathway analysis, FPKM value determination and EST-SSR identification

Annotation of the transcripts was carried out using green plants of non-redundant (nr) protein database NCBI using BLASTX. GO mapping was carried out with BLAST 2GO (GO; http://www.geneontology.org). KEGG maps and an enzyme classification number (EC number) were built for pathway analysis. FPKM values for the transcripts were determined using the formula, FPKM = (Number of reads Mapped × 10⁹) / (Length of Transcript × Total Number of Reads). Here number of reads mapped were calculated by mapping reads on assembled transcripts using CLC Genomics Work bench with a mismatch, insertion, deletion cost of 2, 3 and 3 respectively. Potent EST-SSR markers were identified by MISA, a customized Perl script tool freely available for prediction of SSRs[56].

Protein domains and transcription factor identification in P. hexandrum

Transcripts were searched against a conserved domain database (CDD v3.07) with an E-value cut-off of 0.01 for different domains. For the identification of transcription factor families represented in the P. hexandrum cell culture transcriptome, the transcript contigs were searched against all the transcription factor protein sequences at the plant transcription factor database (http://plntfdb.bio.uni-potsdam.de) using BLASTX with an E-value cutoff 1E-06.

MiRNA target identification in P. hexandrum cell cultures

Conserved miRNAs and their target cDNAs, were found by aligning transcripts against the mature and precursor sequences of known plant miRNAs deposited in miRBase version 19 http://www.mirbase.org/ using CLC Genomic Work bench with a mismatch, insertion, deletion cost of 2, 3 and 3 respectively.

Lignan extraction and high performance liquid chromatography (HPLC) analysis

Lignans were extracted from P. hexandrum cells[57]. In brief, 100 mg of cells were extracted with 2 ml ethanol for 20 min at 60°C in microtubes and sonicated for 15 min. The supernatant was collected after centrifugation and evaporated to dryness under a vacuum. Extracts were dissolved in methanol and used for HPLC analysis. Podophyllotoxin (Sigma-Aldrich, Bangalore, India) was used as a standard. Podophyllotoxin extractions were performed with three biological replicates.

For HPLC, a Waters 2998 photodiode-array detector was set at 290 nm, and separation was carried out using an XTerra RP18, 5 μm, (4.6 × 250 mm i.d.) column. Data analyses were performed with Empower 2 software. Chromatographic conditions were essentially as previously described[7] and standardized in our laboratory[29].

Quantitative RT-PCR (qRT-PCR)

Total RNA was extracted with a Purelink miRNA isolation kit (Invitrogen) using three biological replicates. RNA (2 μg) was reverse transcribed with oligo(dT) primers using RevertAid H Minus First Strand cDNA Synthesis kit (Fermentas, USA). PCR amplification was performed with Power SYBR Green PCR Master Mix (Applied Biosystems, Japan) on an Applied Biosystems 7500 Real-Time PCR System (Applied Biosystems, USA). Relative expression levels were calculated using the ∆-∆Ct method. All primers for qRT-PCR of selected phenylpropanoid pathway genes have been designed from the sequences obtained from optimized Newbler assembly (Additional file20). Actin primers were designed as reported (Acc. No. FL640971.1), Forward primer: 5′-CTCGGGAGGTGCCACCACC-3′ and Reverse primer 5′-GATGGAAGCTGCTGGGTATTCA-3′.