BACKGROUND
ginseng belongs to species within the genus panax  that comprises approximately  <dig> species of slow-growing perennial plants with fleshy roots  <cit> .
north american ginseng  is native to eastern north america  <cit>  and panax species have been used for thousands of years to improve the overall health of human beings, as a remedy to promote vitality, assist the body functions, improve the immune system and protect against stress.
it also has been recommended for years as a traditional medicine for a multitude of diseases such as cancer, inflammation, diabetes, cardiovascular ills and obesity  <cit>  as well as being a recent source of natural extracts mass marketed as a cold prevention treatment  <cit> .
ginsenosides are considered to be the major bioactive compounds behind many of the claims of ginseng’s health benefits; they are triterpenoid saponins found nearly exclusively in ginseng and have been the target of considerable research effort  <cit> .
to date, more than  <dig> naturally occurring ginsenosides have been isolated from panax species and most of them can be classified into two groups based on the skeleton of their aglycones, namely dammarane-type and oleanane-type  <cit> .
the dammarane-type consists mainly of three varieties, classified according to their genuine aglycone moieties: 20s-protopanaxadiol , 20s-protopanaxatriol , and ocotillol.
rb <dig>  rb <dig>  rc, rd  and re and rg <dig>  are the most abundant six ginsenosides found in north american ginseng.
over 90% of total ginsenoside content from north american ginseng belongs to these two groups  <cit> .
north american ginseng contains high levels of rb <dig>  rd and re ginsenosides—higher than those of panax ginseng <cit> .
ro is the only saponin of the oleanane-type ginsenoside, found as a minor component in north american ginseng  <cit> .
ginsenosides are biosynthesized via the mevalonate pathway  <cit> .
using expressed sequence tag analysis, it has been possible to identify several candidate genes encoding for the enzymes farnesyl diphosphate synthase  and squalene synthase , involved in the various biosynthetic steps from isopentenyl pyrophosphate and dimethylallyl pyrophosphate to squalene  <cit> .
the cyclization of oxidosqualene is the branch point for the biosynthesis of ginsenosides and plant sterols.
the common steps from acetyle-coa to  <dig>  3-oxidosqualene have been widely studied  <cit> .
the  <dig>  3-oxidosqualene cyclases  that synthesize β-amyrin and dammarenediol-ii  <cit> , as well as the cyt p <dig> enzyme cyp716a <dig> that catalyses the formation of protopanaxadiol from dammarenediol-ii during ginsenoside biosynthesis have been found in panax ginseng <cit> .
however, the rest of the downstream pathways of ginsenoside biosynthesis remain largely unexplored.
over the past several years, next-generation sequencing technologies have revolutionized the analysis of genomic information  <cit> .
as applied to the transcriptome with rnaseq, it has been successfully used for transcript profiling, as well as snp discovery in a number of plant species and has dramatically improved the efficiency and speed of gene discovery  <cit> .
the application of  <dig> next generation sequencing technology has seen a rapid improvement in throughput, read length and accuracy in the past few years, with the gs flx titanium server used in this study able to generate one million reads with an average length of  <dig> bases with  <dig> % accuracy  <cit> .
meanwhile, the analysis of transcriptomic data often relies on aligning reads to a reference which is often not feasible for non-model plant species in which little genomic research has been performed.
we applied the transcriptome assembly program trinity to the assembly of a high quality reference transcriptome for north american ginseng.
trinity has been shown to recover most expressed transcripts as full-length sequences, and is also able to resolve alternative isoforms and duplicated genes, outperforming other de novo assembly tools  <cit> .
our application to ginseng resulted in  <dig>  ginseng root transcripts.
we fully annotated 84% of these transcripts using sequence similarity searches and protein domain scanning with publicly available databases.
in our results, we were able to identify predicted representatives for all of the known enzymes involved in the ginsenoside backbone biosynthesis and also profile their expression levels across seasonal development.
RESULTS
commercial production of ginseng usually results in harvests after  <dig> to 5 years.
in this study, three-year-old roots were collected, washed, sorted for uniformity and overwintered and grown under simulated growing season conditions in the biotron facility to minimize variation in environmental factors and soil pathogen infection.
root samples were collected over the fourth-year full growing season at seven development stages  <cit> : 1) budding - leaf emergence; 2) leaves - distinctive separation of leaf and stem; 3) flowering - plant in full flower; 4) green fruit set; 5) ripe fruit - fruit coloration fully red; 6) fruit drop – including early signs of leaf senescence  and 7) senescence - complete senescence of leaf and stem  .
sequencing and de novo assembly
tissue from the seven developing stages was collected from 4-year-old north american ginseng grown under controlled conditions in a growth chamber programmed to simulate southern ontario, canada commercial growing conditions .
ginsenoside content levels for each development stage in american ginseng are listed in table  <dig>  total rna was isolated from the roots of each stage, and messenger rna purified by ambion poly  purist™ mrna purification kits.
a half plate of  <dig> sequencing on the gs flx titanium platform was applied to each sample .
this generated from  <dig> to  <dig> mbp of sequence for each stage with an average sequence length of  <dig> base pairs.
all sequencing reads from the seven development stages were deposited in the ncbi  and can be accessed at the sequencing read archive  with the accession numbers srx <dig>  srx <dig>  srx <dig>  srx <dig>  srx <dig>  srx <dig> and srx <dig> for stages one through seven respectively.
a: samples were tested in duplicate and number of root samples tested in each stage.
filtering of the sequence data before assembly included steps that removed plastid contaminants and adaptor sequences, as well as trimmed the base pair bias present in the first 15 bp of the 5’ end and low quality bases  at the 3’ end.
quality score distributions for reads in each stage before and after filtering are provided in additional file  <dig>  after filtering,  <dig> , <dig> sequence reads remained for assembly with an average length of  <dig> base pairs .
unfortunately, our poly-a purification step failed to effectively filter the abundant rrna present and a large percentage of reads had to be filtered from each stage as contaminants arising from ribosomal or plastid rna sequences.
interestingly, after quality trimming the reads, a disproportionate amount of sequences with lengths of approximately  <dig> and  <dig> base pairs were found .
this is presumably an artefact of trends in sequence quality drop off at specific points during sequencing.
high quality reads from all stages were combined and provided to the transcriptome assembly program trinity.
the resulting assembled  <dig>  mb of transcriptome contained  <dig>  genes or  <dig>  transcripts when including the different gene isoforms trinity is capable of returning.
this number represents more than  <dig> times the number of north american ginseng sequences presently deposited in genbank .
transcript lengths ranged from  <dig> to  <dig>  base pairs with an average length of 896 bp and the majority of transcripts ranging between 500 bp and 2kb in size .
almost half of all genes assembled  possessed at least one isoform, with a total of  <dig>  splice variants identified by trinity and 11% of genes possessing  <dig> or more splice variants .
one gene  possessed  <dig> different isoforms, although we felt this may have been an artefact of the assembly process.
in a similarity comparison to  <dig>  panax quinquefolius ests in genbank,  <dig> % were present in our assembly with strong significance .
when genbank ests specifically derived from the panax quinquefolius rhizome were considered, this number increased to  <dig> %, suggesting a high quality, comprehensive sampling of the root developmental transcriptome.
to simplify identification and enable easy reference, all sequences in the assembly were assigned a unique identifier derived from the trinity graph component and appended with a splice number that followed the form of “pqx.y”, where “pq” stands for panax quinquefolias, “x” is the trinity component number and “y” is the splice variant number.
transcript annotation with public databases
to facilitate as complete an annotation as possible for the assembly, sequence similarity searches  <cit>  were performed against a collection  <dig>  ginseng ests from genbank, the arabidopsis genome , the uniprot plant protein annotation program  database and genbank’s non-redundant  protein database.
in addition, protein domain scanning using hidden markov models  from pfam were applied as well as the assignment of metabolic pathway information from the kyoto encyclopedia of genes and genomes  database.
overall, these efforts annotated  <dig> % of transcripts in the assembly.
comparison against the genbank nr database yielded  <dig>  hits , accounting for  <dig> % of transcriptome sequences.
a full  <dig> % of these hits showed strong homology .
as the nr database contains few ginseng proteins, we examined the different species with which homology was found.
the majority of hits  were found to be against grape , followed by castor oil plant  , black cottonwood poplar  , and soybean  .
similar results were found with searches against the plant protein annotation program  database from uniprot and the tair <dig> release of the arabidopsis genome yielding  <dig>   and  <dig>   transcripts with significant hits  .
as arabidopsis is the most thoroughly annotated plant, sequence homology to arabidopsis was also used to characterize the transcriptome based on gene ontology  information.
go annotation provides descriptions of gene products in terms of their associated molecular functions, cellular components, and biological processes.
using sequence homology to tair <dig>   <dig>  go terms were assigned to  <dig>  sequences categorized into  <dig> functional groups.
go assignments were most frequently related to biological processes , followed by cellular components  and molecular function  .
the assembly was scanned with protein domain hmm models from the pfam database in order to catalogue any significant matches  to known protein domains.
overall,  <dig>  hmms were scanned against the assembly resulting in annotation for  <dig>  transcripts possessing  <dig>  different protein domains.
this added annotation information to an additional  <dig> transcripts that had no hits in the prior homology searches.
the most abundant domain found was the protein kinase domain, present in  <dig>  transcripts.
this is a similar number to the  <dig>  kinases present in the arabidopsis genome and not surprising, as protein kinases play a role in a multitude of cellular processes, including division, proliferation, apoptosis, and differentiation.
finally, in order to assign metabolic information to our transcripts, the kaas tool  <cit>  was used to assign pathway information from the kegg database.
this resulted in a kegg orthology  number for  <dig>  transcripts that possessed homology with metabolic enzymes in the kegg database.
the pathways most strongly represented in the results were protein processing in the endoplasmic reticulum, ubiquitin mediated proteolysis, nucleotide excision repair, arginine and proline metabolism, peroxisome, amino sugar and nucleotide sugar metabolism, proteasome, starch and sucrose metabolism and glycolysis/gluconeogenesis.
all annotation information and significance scores were summarized and concatenated with transcript identifiers into a single line of annotation for each sequence in a fasta formatted transcriptome file .
expression profiles across development
to evaluate the contribution of each stage of development to the reference transcriptome and profile the relative expression levels of transcripts throughout development, we mapped all high quality reads from each stage back to the assembled transcriptome using the bwa tool  <cit> .
while this is useful in determining relative transcript abundance, it did have the unfortunate side effect of reporting unmapped reads for  <dig>  transcripts.
investigating further, it was found that nearly all unmapped transcripts were involved in alternative splicing and of the shortest transcripts in the assembly .
we presume that this is a consequence of the k-mer  assembly process used by trinity, perhaps creating isoform transcripts that are shorter than the long  <dig> reads from which they originally derive or possibly from read trimming by the alignment program.
nevertheless, expression information was obtained for over 85% of the transcriptome which comprised the longer contigs.
counts for mapped reads were then used to generate rpkm  values for each putative transcript in each developmental stage.
of the most abundantly expressed predicted proteins, storage proteins topped the list, with a number of panax ginseng annotated proteins: rnase-like major storage protein, specific abundant protein, tonoplast intrinsic protein, major latex-like protein and dehydrin  <dig> .
we also examined the level of overlap in transcript expression occurring between major developmental stages of the plant.
to limit the complexity in the number of comparisons needed, the seven developmental stages sampled were grouped into the broader developmental categories of budding , flowering , fruiting  and senescence .
overall,  <dig> ,  <dig> ,  <dig>  and  <dig>  predicted transcripts showed positive expression in each category respectively.
each category possessed around two thousand predicted transcripts uniquely expressed during that stage of root development, with the exception of flowering, which had  <dig> specific transcripts.
however, this is likely influenced by the fact that flowering is the only single stage group and thus represents a snapshot from a shorter period of developmental time compared to the other stages.
interestingly, the senescence stage had the largest number of predicted transcripts overlapping with other developmental stages.
altogether, a total of  <dig>  transcripts were found to be expressed throughout all stages of development .
as ginsenosides are the ultimate compound of commercial interest in ginseng harvest, we sought to inventory all putative transcripts in the assembly that could be implicated in the synthesis of ginsenosides and examine their expression profiles across seasonal development.
we therefore took all kegg orthology numbers associated with enzymes in the mevalonate pathway of terpenoid backbone biosynthesis along with those from the sesquiterpenoid and triterpenoid biosynthesis pathways leading to chair-chair-chair-boat conformation triterpenoids  and pulled all transcripts annotated with corresponding ko numbers in the assembly.
this returned  <dig> predicted genes  annotated and identified by kegg orthology.
in addition, as there was no appropriate assigned kegg ec for β-amyrin synthase, five β-amyrin synthase annotated genes  were pulled for analysis using homology annotation in the assembly .
* this gene possessed  <dig> splice variants and may be an artefact of assembly.
squalene synthase plays an important role as the precursor in backbone biosynthesis of the dammarenediol-type ginsenosides  <cit> .
there are three previously reported squalene synthases in panax ginseng <cit>  including squalene synthase , squalene synthase  <dig>  and squalene synthase  <dig> .
sqs <dig> and sqs <dig> were also previously found in panax quinquefolius <cit> .
although, we found  <dig> putative genes  with a significant squalene/phytoene synthase domain , only one  was annotated with kegg orthology to sqs.
this gene showed strong identity  with a. thaliana the sqs <dig> gene.
in north american ginseng, the majority of ginsenosides are known to be of the dammarene-type ginsenosides produced from protopanaxdiol and protopanaxtriol triterpenes.
dammarenediol-ii produces protopanaxdiol and protopanaxtriol, and ginsenosides are thought to be synthesized from subsequent hydroxylation of these products by cytochrome p <dig> enzymes and glycosylation by glycosyltransferases   <cit> .
our assembly contained  <dig> predicted transcripts annotated as cytochrome p450s; with  <dig> of these possessing high similarity to p <dig> sequences from panax ginseng and panax notoginseng est collections.
similarly, the assembly contained  <dig> predicted transcripts annotated as glycosyltransferases with  <dig> of these previously identified in panax ginseng, panax notoginseng and panax quinquefolias.
in order to identify potential candidates from these gene families that may be functioning in the latter stages of ginsenoside biosynthesis, we conducted a co-expression analysis with the transcript profiles for our putative dammarenediol synthase  and squalene expoxidase  found immediately upstream in the ginsenoside biosynthesis pathway.
the expression profiles for the two strongest sqe  and ds  annotated transcripts in our assembly showed very high co-expression .
this is not unexpected given that ds follows sqe in the biosynthesis pathway.
we reasoned that candidate downstream p <dig> and glycosyltransferase genes may be similarly co-expressed.
co-expression analysis between all putative p450s and our predicted ds transcript identified  <dig> candidate p450s highly co-expressed with ds  across the  <dig> developmental stages sampled .
in the case of pq <dig> , co-expression with ds was extremely high .
similarly,  <dig> glycosyltransferases annotated transcripts were highly co-expressed with our predicted ds  .
as before, one transcript  showed extremely high co-expression .
similar results were found in co-expression with the upstream sqe .
it should be noted that two splice variants from a single gene were also returned in the glycosyltransferase co-expression analysis.
however, these transcripts derived from a gene  for which  <dig> splice variants were identified by trinity.
as the annotation information for these  <dig> isoforms showed no strong consistency across isoforms, we concluded that this may be an artefact of the assembly, perhaps from the misinterpretation of a gene family, and were therefore dropped to simplify the analysis.
in order to provide biological validation to the  <dig> predicted transcripts showing coexpression with our putative ds gene, we conducted real-time pcr against the ds gene itself, the  <dig> putative glycosyltransferase genes and  <dig> putative cytochrome p <dig> genes across the final  <dig> stages of development.
we were able to successfully amplify all  <dig> genes using rt-pcr analysis and confirmed strong coexpression between our predicted ds,  <dig> of the glycosyltransferase transcripts and  <dig> of the p <dig> transcripts .
additional file  <dig> lists the primers used in this analysis.
before examining the relative expression profiles of our identified biosynthesis genes we sought to first look at the expression profiles for the entire transcriptome across seasonal development.
rpkm values for all predicted transcripts were thus hierarchically clustered, and displayed in a heat map to create a transcriptome-wide display of developmental expression .
many transcripts showed maximum expression within one or two specific stages of root development relative to other stages, with all stages of development possessing distinct clusters of genes showing dominant expression within that stage.
moreover, compared to other stages of development, senescence possessed a disproportionately large cluster of transcripts exhibiting a senescence-specific maximum in expression.
expression profiles for all predicted genes and their isoforms implicated in ginsenoside biosynthesis  as identified above were also hierarchically clustered and plotted in a heat map .
for representatives of downstream candidates, the most highly co-expressed p <dig> and glycosyltransferase transcripts identified in the co-expression analysis were also included in the collection.
as seen with the entire transcriptome, many enzymes displayed abundant expression, specific to one or a few developmental stages relative to the other stages.
to our surprise, stage  <dig>  and stage  <dig>  both possessed obvious clusters that encompassed putative representatives for almost all the enzymes in the biosynthesis pathway, suggesting that these two developmental stages are important points of ginsenoside biosynthesis in the developing plant.
in the case of stage  <dig>  this makes intuitive sense considering the traditional time of ginseng harvest after the fruit drop stage and onset of leaf senescence.
interestingly, when different isoforms of genes were present, they tended to exhibit expression maxima in different stages of development.
this was most apparent with the  <dig> different isoforms of sqe <dig> , for which at least one instance of an isoform with stage specific expression was found in all developmental stages except senescence stage  <dig> .
of course, due to the computational limitations of assembly and mapping, both the isoforms identified and the expression profiles should be interpreted with caution without further in planta validation.
discussion
ginsenosides are found exclusively in the plant genus panax and are viewed as the primary active compounds behind the claims of ginseng's health efficacy.
the ginsenosides, also known as triterpenoids belong to the class of terpenoid saponines and are synthesized via the mevalonate pathway, which shares a common pathway with sterol.
oxidosqualene is a precursor common to the biosynthesis of both steroid and triterpenoids in higher plants  <cit> .
the cyclization of  <dig>  3-oxidosqualene to oleanane and dammarane triterpenes is the first step in ginsenosides synthesis.
the primary enzymes of triterpenoid saponions synthases are β-amyrin synthase , and lupeol or dammarenediol-ii synthase .
putative orthologs for ds and as were both identified in our assembly.
β-amyrin synthase synthesizes oleanane-type triterpene and oleanane-type ginsenosides; however, they are found in only minor amounts in north american ginseng.
on the contrary, over 90% of the total ginsenosides of north american ginseng belong to the dammarane-type.
interestingly, in our results, predicted as transcripts appear most abundantly within developmental stages in which predicted ds enzymes are not abundant .
moreover, they are not associated with the distinctive co-expression clusters of abundant putative ginsenoside biosynthesis genes as seen with ds in stages  <dig> and  <dig>  this potentially emphasizes the minor role of as derived ginsenosides in north american ginseng.
at the same time, the clusters of predicted biosynthesis enzymes during flowering and fruit drop highlight a potential significance of damamarane-type synthesis and suggest that flowering and fruit drop could be primary points of ginsenoside biosynthesis during seasonal development, with most emphasis on the latter.
while we observed distinct stages of specific up-regulation for most predicted enzymes in the assembly, the presence of different isoforms with expression abundance specific to different stages of development is intriguing.
in particular, our putative sqe <dig> possessed  <dig> different isoforms that showed stage specific expression abundance in all stages of development save stage  <dig> .
curiously, an alternate predicted sqe <dig>  with only a single isoform was strongly expressed in this stage.
this seems to suggest that variation in alternative splicing may be a mechanism for producing varied forms of ginsenosides during seasonal development.
generally speaking, it also implies that alternative splicing may function as a means for directing variation in secondary metabolite production throughout the course of plant development.
however, due to the inherent computational limitations involved in assembly and mapping, further analysis in the form of qpcr and associated metabolic assays is needed to prove or disprove any such hypotheses.
of course, the assembly process will not be perfect with respect to isoform prediction and the transcripts themselves.
there is a strong potential for misassembly in the form of merged gene families, close paralogs, or even alleles of the same gene being misreported as isoforms.
while the comparison with ginseng ests in genbank is reassuring of the assembly quality, all predictions should be treated with caution in the absence of biological validation.
similarly, the mapping of reads to the assembly is limited by the presence of isoforms, as the actual point of origin for the read is confounded by the presence of potential multiple sources.
this introduces a level of stochastic noise to the expression analysis that is mostly confined to genes with several isoforms.
that said, real-time pcr was able to validate the presence of a number of transcripts within expected developmental stages, as well as confirm their coexpression and upregulation within the fruit drop stage of development.
transcripts for a predicted ds gene, six putative p450s and six putative glycosyltransferases were all confirmed as present in planta and expression levels for four of the p450s and five of the glycosyltransferases confirm a tight coregulation with the predicted ds gene across the last stages of development as seen in our expression data.
these predicted enzymes are thus strong candidates for controlling ginsenoside biosynthesis in the late stages of plant development.
while this study, as any next-generation sequencing study, would have benefited from an even larger amount of sequence data, additional sequencing over so many stages of development was unfortunately cost prohibitive.
nevertheless, the overall structure of the assembly with regard to number of genes, isoform frequency, length of transcripts and level of homology with existing est libraries and annotation patterns identified among the transcripts is all very supportive and indicative of a strong representation of the biology.
overall, we believe the analysis benefited considerably from the use of the much longer reads that  <dig> sequencing is capable of generating.
this added information translates into more reliable, longer, and complete transcripts, as well as more information for improved accuracy in the calling of alternative splicing among sequenced transcripts.
recent work by sun et al., found similar benefits to the use of  <dig> sequencing in the assembly of a transcriptome for north american ginseng  <cit> .
our work has expanded on these initial efforts through the generation of over six times the amount of high quality  <dig> reads sampled over seven stages of seasonal root development.
this provides an unprecedented snapshot of north american ginseng root development and a means to comprehensively capture root expressed transcripts and their relative abundance through developmental time.
in addition, we have also taken advantage of the recent trinity assembly algorithm, which is designed specifically for assembling transcriptomes.
altogether, this provides the benefit of a high confidence, reliable assembly of putative transcripts that includes potentially valuable splice variant information.
in the end, the combination of a fully annotated assembly along with corresponding expression data as sampled across developmental time has provided for an insightful look into the transcriptome of north american ginseng.
CONCLUSIONS
ginseng is a valuable medicinal herb that is widely used in traditional medicine and has considerable commercial value.
this study has focused on the ginseng root as it is the tissue almost exclusively used for medicinal applications and contains the majority of the commercially and medicinally valuable ginsenosides.
using extensive sampling of root tissue across seasonal development and next-generation sequencing we have generated a reliable north american ginseng transcriptome with extensive annotation and isoform information.
expression profiling of all ginsenoside backbone biosynthesis genes suggests that ginsenoside biosynthesis is tightly linked with the developmental stages of flowering and fruit drop.
we anticipate the assembly and associated expression data should prove invaluable to those conducting transcriptomic and metabolomic research against panax quinquefolius.

 ginseng belongs to species within the genus panax that comprises approximately species of slow-growing perennial plants with fleshy roots. north american ginseng is native to eastern north america and panax species have been used for thousands of years to improve the overall health of human beings, as a remedy to promote vitality, assist the body functions, improve the immune system and protect against stress. it also has been recommended for years as a traditional medicine for a multitude of diseases such as cancer, inflammation, diabetes, cardiovascular ills and obesity as well as being a recent source of natural extracts mass marketed as a cold prevention treatment. ginsenosides are considered to be the major bioactive compounds behind many of the claims of ginseng’s health benefits; they are triterpenoid saponins found nearly exclusively in ginseng and have been the target of considerable research effort. to date, more than naturally occurring ginsenosides have been isolated from panax species and most of them can be classified into two groups based on the skeleton of their aglycones, namely dammarane-type and oleanane-type. the dammarane-type consists mainly of three varieties, classified according to their genuine aglycone moieties: 20s-protopanaxadiol, 20s-protopanaxatriol, and ocotillol. rb rb rc, rd and re and rg are the most abundant six ginsenosides found in north american ginseng. over 90% of total ginsenoside content from north american ginseng belongs to these two groups. north american ginseng contains high levels of rb rd and re ginsenosides—higher than those of panax ginseng. ro is the only saponin of the oleanane-type ginsenoside, found as a minor component in north american ginseng. ginsenosides are biosynthesized via the mevalonate pathway. using expressed sequence tag analysis, it has been possible to identify several candidate genes encoding for the enzymes farnesyl diphosphate synthase and squalene synthase, involved in the various biosynthetic steps from isopentenyl pyrophosphate and dimethylallyl pyrophosphate to squalene. the cyclization of oxidosqualene is the branch point for the biosynthesis of ginsenosides and plant sterols. the common steps from acetyle-coa to 3-oxidosqualene have been widely studied. the 3-oxidosqualene cyclases that synthesize β-amyrin and dammarenediol-ii, as well as the cyt p enzyme cyp716a that catalyses the formation of protopanaxadiol from dammarenediol-ii during ginsenoside biosynthesis have been found in panax ginseng. however, the rest of the downstream pathways of ginsenoside biosynthesis remain largely unexplored. over the past several years, next-generation sequencing technologies have revolutionized the analysis of genomic information. as applied to the transcriptome with rnaseq, it has been successfully used for transcript profiling, as well as snp discovery in a number of plant species and has dramatically improved the efficiency and speed of gene discovery. the application of next generation sequencing technology has seen a rapid improvement in throughput, read length and accuracy in the past few years, with the gs flx titanium server used in this study able to generate one million reads with an average length of bases with % accuracy. meanwhile, the analysis of transcriptomic data often relies on aligning reads to a reference which is often not feasible for non-model plant species in which little genomic research has been performed. we applied the transcriptome assembly program trinity to the assembly of a high quality reference transcriptome for north american ginseng. trinity has been shown to recover most expressed transcripts as full-length sequences, and is also able to resolve alternative isoforms and duplicated genes, outperforming other de novo assembly tools. our application to ginseng resulted in ginseng root transcripts. we fully annotated 84% of these transcripts using sequence similarity searches and protein domain scanning with publicly available databases. in our results, we were able to identify predicted representatives for all of the known enzymes involved in the ginsenoside backbone biosynthesis and also profile their expression levels across seasonal development.  commercial production of ginseng usually results in harvests after to 5 years. in this study, three-year-old roots were collected, washed, sorted for uniformity and overwintered and grown under simulated growing season conditions in the biotron facility to minimize variation in environmental factors and soil pathogen infection. root samples were collected over the fourth-year full growing season at seven development stages: 1) budding - leaf emergence; 2) leaves - distinctive separation of leaf and stem; 3) flowering - plant in full flower; 4) green fruit set; 5) ripe fruit - fruit coloration fully red; 6) fruit drop – including early signs of leaf senescence and 7) senescence - complete senescence of leaf and stem. sequencing and de novo assembly tissue from the seven developing stages was collected from 4-year-old north american ginseng grown under controlled conditions in a growth chamber programmed to simulate southern ontario, canada commercial growing conditions. ginsenoside content levels for each development stage in american ginseng are listed in table  total rna was isolated from the roots of each stage, and messenger rna purified by ambion poly purist™ mrna purification kits. a half plate of sequencing on the gs flx titanium platform was applied to each sample. this generated from to mbp of sequence for each stage with an average sequence length of base pairs. all sequencing reads from the seven development stages were deposited in the ncbi and can be accessed at the sequencing read archive with the accession numbers srx srx srx srx srx srx and srx for stages one through seven respectively. a: samples were tested in duplicate and number of root samples tested in each stage. filtering of the sequence data before assembly included steps that removed plastid contaminants and adaptor sequences, as well as trimmed the base pair bias present in the first 15 bp of the 5’ end and low quality bases at the 3’ end. quality score distributions for reads in each stage before and after filtering are provided in additional file after filtering,, sequence reads remained for assembly with an average length of base pairs. unfortunately, our poly-a purification step failed to effectively filter the abundant rrna present and a large percentage of reads had to be filtered from each stage as contaminants arising from ribosomal or plastid rna sequences. interestingly, after quality trimming the reads, a disproportionate amount of sequences with lengths of approximately and base pairs were found. this is presumably an artefact of trends in sequence quality drop off at specific points during sequencing. high quality reads from all stages were combined and provided to the transcriptome assembly program trinity. the resulting assembled  mb of transcriptome contained genes or transcripts when including the different gene isoforms trinity is capable of returning. this number represents more than times the number of north american ginseng sequences presently deposited in genbank. transcript lengths ranged from to base pairs with an average length of 896 bp and the majority of transcripts ranging between 500 bp and 2kb in size. almost half of all genes assembled possessed at least one isoform, with a total of splice variants identified by trinity and 11% of genes possessing or more splice variants. one gene possessed different isoforms, although we felt this may have been an artefact of the assembly process. in a similarity comparison to panax quinquefolius ests in genbank, % were present in our assembly with strong significance. when genbank ests specifically derived from the panax quinquefolius rhizome were considered, this number increased to %, suggesting a high quality, comprehensive sampling of the root developmental transcriptome. to simplify identification and enable easy reference, all sequences in the assembly were assigned a unique identifier derived from the trinity graph component and appended with a splice number that followed the form of “pqx.y”, where “pq” stands for panax quinquefolias, “x” is the trinity component number and “y” is the splice variant number. transcript annotation with public databases to facilitate as complete an annotation as possible for the assembly, sequence similarity searches were performed against a collection ginseng ests from genbank, the arabidopsis genome, the uniprot plant protein annotation program database and genbank’s non-redundant protein database. in addition, protein domain scanning using hidden markov models from pfam were applied as well as the assignment of metabolic pathway information from the kyoto encyclopedia of genes and genomes database. overall, these efforts annotated % of transcripts in the assembly. comparison against the genbank nr database yielded hits, accounting for % of transcriptome sequences. a full % of these hits showed strong homology. as the nr database contains few ginseng proteins, we examined the different species with which homology was found. the majority of hits were found to be against grape, followed by castor oil plant, black cottonwood poplar, and soybean. similar results were found with searches against the plant protein annotation program database from uniprot and the tair release of the arabidopsis genome yielding and transcripts with significant hits. as arabidopsis is the most thoroughly annotated plant, sequence homology to arabidopsis was also used to characterize the transcriptome based on gene ontology information. go annotation provides descriptions of gene products in terms of their associated molecular functions, cellular components, and biological processes. using sequence homology to tair go terms were assigned to sequences categorized into functional groups. go assignments were most frequently related to biological processes, followed by cellular components and molecular function. the assembly was scanned with protein domain hmm models from the pfam database in order to catalogue any significant matches to known protein domains. overall, hmms were scanned against the assembly resulting in annotation for transcripts possessing different protein domains. this added annotation information to an additional transcripts that had no hits in the prior homology searches. the most abundant domain found was the protein kinase domain, present in transcripts. this is a similar number to the kinases present in the arabidopsis genome and not surprising, as protein kinases play a role in a multitude of cellular processes, including division, proliferation, apoptosis, and differentiation. finally, in order to assign metabolic information to our transcripts, the kaas tool was used to assign pathway information from the kegg database. this resulted in a kegg orthology number for transcripts that possessed homology with metabolic enzymes in the kegg database. the pathways most strongly represented in the results were protein processing in the endoplasmic reticulum, ubiquitin mediated proteolysis, nucleotide excision repair, arginine and proline metabolism, peroxisome, amino sugar and nucleotide sugar metabolism, proteasome, starch and sucrose metabolism and glycolysis/gluconeogenesis. all annotation information and significance scores were summarized and concatenated with transcript identifiers into a single line of annotation for each sequence in a fasta formatted transcriptome file. expression profiles across development to evaluate the contribution of each stage of development to the reference transcriptome and profile the relative expression levels of transcripts throughout development, we mapped all high quality reads from each stage back to the assembled transcriptome using the bwa tool. while this is useful in determining relative transcript abundance, it did have the unfortunate side effect of reporting unmapped reads for transcripts. investigating further, it was found that nearly all unmapped transcripts were involved in alternative splicing and of the shortest transcripts in the assembly. we presume that this is a consequence of the k-mer assembly process used by trinity, perhaps creating isoform transcripts that are shorter than the long reads from which they originally derive or possibly from read trimming by the alignment program. nevertheless, expression information was obtained for over 85% of the transcriptome which comprised the longer contigs. counts for mapped reads were then used to generate rpkm values for each putative transcript in each developmental stage. of the most abundantly expressed predicted proteins, storage proteins topped the list, with a number of panax ginseng annotated proteins: rnase-like major storage protein, specific abundant protein, tonoplast intrinsic protein, major latex-like protein and dehydrin. we also examined the level of overlap in transcript expression occurring between major developmental stages of the plant. to limit the complexity in the number of comparisons needed, the seven developmental stages sampled were grouped into the broader developmental categories of budding, flowering, fruiting and senescence. overall,,, and predicted transcripts showed positive expression in each category respectively. each category possessed around two thousand predicted transcripts uniquely expressed during that stage of root development, with the exception of flowering, which had specific transcripts. however, this is likely influenced by the fact that flowering is the only single stage group and thus represents a snapshot from a shorter period of developmental time compared to the other stages. interestingly, the senescence stage had the largest number of predicted transcripts overlapping with other developmental stages. altogether, a total of transcripts were found to be expressed throughout all stages of development. as ginsenosides are the ultimate compound of commercial interest in ginseng harvest, we sought to inventory all putative transcripts in the assembly that could be implicated in the synthesis of ginsenosides and examine their expression profiles across seasonal development. we therefore took all kegg orthology numbers associated with enzymes in the mevalonate pathway of terpenoid backbone biosynthesis along with those from the sesquiterpenoid and triterpenoid biosynthesis pathways leading to chair-chair-chair-boat conformation triterpenoids and pulled all transcripts annotated with corresponding ko numbers in the assembly. this returned predicted genes annotated and identified by kegg orthology. in addition, as there was no appropriate assigned kegg ec for β-amyrin synthase, five β-amyrin synthase annotated genes were pulled for analysis using homology annotation in the assembly. * this gene possessed splice variants and may be an artefact of assembly. squalene synthase plays an important role as the precursor in backbone biosynthesis of the dammarenediol-type ginsenosides. there are three previously reported squalene synthases in panax ginseng including squalene synthase, squalene synthase and squalene synthase. sqs and sqs were also previously found in panax quinquefolius. although, we found putative genes with a significant squalene/phytoene synthase domain, only one was annotated with kegg orthology to sqs. this gene showed strong identity with a. thaliana the sqs gene. in north american ginseng, the majority of ginsenosides are known to be of the dammarene-type ginsenosides produced from protopanaxdiol and protopanaxtriol triterpenes. dammarenediol-ii produces protopanaxdiol and protopanaxtriol, and ginsenosides are thought to be synthesized from subsequent hydroxylation of these products by cytochrome p enzymes and glycosylation by glycosyltransferases. our assembly contained predicted transcripts annotated as cytochrome p450s; with of these possessing high similarity to p sequences from panax ginseng and panax notoginseng est collections. similarly, the assembly contained predicted transcripts annotated as glycosyltransferases with of these previously identified in panax ginseng, panax notoginseng and panax quinquefolias. in order to identify potential candidates from these gene families that may be functioning in the latter stages of ginsenoside biosynthesis, we conducted a co-expression analysis with the transcript profiles for our putative dammarenediol synthase and squalene expoxidase found immediately upstream in the ginsenoside biosynthesis pathway. the expression profiles for the two strongest sqe and ds annotated transcripts in our assembly showed very high co-expression. this is not unexpected given that ds follows sqe in the biosynthesis pathway. we reasoned that candidate downstream p and glycosyltransferase genes may be similarly co-expressed. co-expression analysis between all putative p450s and our predicted ds transcript identified candidate p450s highly co-expressed with ds across the developmental stages sampled. in the case of pq, co-expression with ds was extremely high. similarly, glycosyltransferases annotated transcripts were highly co-expressed with our predicted ds. as before, one transcript showed extremely high co-expression. similar results were found in co-expression with the upstream sqe. it should be noted that two splice variants from a single gene were also returned in the glycosyltransferase co-expression analysis. however, these transcripts derived from a gene for which splice variants were identified by trinity. as the annotation information for these isoforms showed no strong consistency across isoforms, we concluded that this may be an artefact of the assembly, perhaps from the misinterpretation of a gene family, and were therefore dropped to simplify the analysis. in order to provide biological validation to the predicted transcripts showing coexpression with our putative ds gene, we conducted real-time pcr against the ds gene itself, the putative glycosyltransferase genes and putative cytochrome p genes across the final stages of development. we were able to successfully amplify all genes using rt-pcr analysis and confirmed strong coexpression between our predicted ds, of the glycosyltransferase transcripts and of the p transcripts. additional file lists the primers used in this analysis. before examining the relative expression profiles of our identified biosynthesis genes we sought to first look at the expression profiles for the entire transcriptome across seasonal development. rpkm values for all predicted transcripts were thus hierarchically clustered, and displayed in a heat map to create a transcriptome-wide display of developmental expression. many transcripts showed maximum expression within one or two specific stages of root development relative to other stages, with all stages of development possessing distinct clusters of genes showing dominant expression within that stage. moreover, compared to other stages of development, senescence possessed a disproportionately large cluster of transcripts exhibiting a senescence-specific maximum in expression. expression profiles for all predicted genes and their isoforms implicated in ginsenoside biosynthesis as identified above were also hierarchically clustered and plotted in a heat map. for representatives of downstream candidates, the most highly co-expressed p and glycosyltransferase transcripts identified in the co-expression analysis were also included in the collection. as seen with the entire transcriptome, many enzymes displayed abundant expression, specific to one or a few developmental stages relative to the other stages. to our surprise, stage and stage both possessed obvious clusters that encompassed putative representatives for almost all the enzymes in the biosynthesis pathway, suggesting that these two developmental stages are important points of ginsenoside biosynthesis in the developing plant. in the case of stage this makes intuitive sense considering the traditional time of ginseng harvest after the fruit drop stage and onset of leaf senescence. interestingly, when different isoforms of genes were present, they tended to exhibit expression maxima in different stages of development. this was most apparent with the different isoforms of sqe, for which at least one instance of an isoform with stage specific expression was found in all developmental stages except senescence stage. of course, due to the computational limitations of assembly and mapping, both the isoforms identified and the expression profiles should be interpreted with caution without further in planta validation. discussion ginsenosides are found exclusively in the plant genus panax and are viewed as the primary active compounds behind the claims of ginseng's health efficacy. the ginsenosides, also known as triterpenoids belong to the class of terpenoid saponines and are synthesized via the mevalonate pathway, which shares a common pathway with sterol. oxidosqualene is a precursor common to the biosynthesis of both steroid and triterpenoids in higher plants. the cyclization of 3-oxidosqualene to oleanane and dammarane triterpenes is the first step in ginsenosides synthesis. the primary enzymes of triterpenoid saponions synthases are β-amyrin synthase, and lupeol or dammarenediol-ii synthase. putative orthologs for ds and as were both identified in our assembly. β-amyrin synthase synthesizes oleanane-type triterpene and oleanane-type ginsenosides; however, they are found in only minor amounts in north american ginseng. on the contrary, over 90% of the total ginsenosides of north american ginseng belong to the dammarane-type. interestingly, in our results, predicted as transcripts appear most abundantly within developmental stages in which predicted ds enzymes are not abundant. moreover, they are not associated with the distinctive co-expression clusters of abundant putative ginsenoside biosynthesis genes as seen with ds in stages and this potentially emphasizes the minor role of as derived ginsenosides in north american ginseng. at the same time, the clusters of predicted biosynthesis enzymes during flowering and fruit drop highlight a potential significance of damamarane-type synthesis and suggest that flowering and fruit drop could be primary points of ginsenoside biosynthesis during seasonal development, with most emphasis on the latter. while we observed distinct stages of specific up-regulation for most predicted enzymes in the assembly, the presence of different isoforms with expression abundance specific to different stages of development is intriguing. in particular, our putative sqe possessed different isoforms that showed stage specific expression abundance in all stages of development save stage. curiously, an alternate predicted sqe with only a single isoform was strongly expressed in this stage. this seems to suggest that variation in alternative splicing may be a mechanism for producing varied forms of ginsenosides during seasonal development. generally speaking, it also implies that alternative splicing may function as a means for directing variation in secondary metabolite production throughout the course of plant development. however, due to the inherent computational limitations involved in assembly and mapping, further analysis in the form of qpcr and associated metabolic assays is needed to prove or disprove any such hypotheses. of course, the assembly process will not be perfect with respect to isoform prediction and the transcripts themselves. there is a strong potential for misassembly in the form of merged gene families, close paralogs, or even alleles of the same gene being misreported as isoforms. while the comparison with ginseng ests in genbank is reassuring of the assembly quality, all predictions should be treated with caution in the absence of biological validation. similarly, the mapping of reads to the assembly is limited by the presence of isoforms, as the actual point of origin for the read is confounded by the presence of potential multiple sources. this introduces a level of stochastic noise to the expression analysis that is mostly confined to genes with several isoforms. that said, real-time pcr was able to validate the presence of a number of transcripts within expected developmental stages, as well as confirm their coexpression and upregulation within the fruit drop stage of development. transcripts for a predicted ds gene, six putative p450s and six putative glycosyltransferases were all confirmed as present in planta and expression levels for four of the p450s and five of the glycosyltransferases confirm a tight coregulation with the predicted ds gene across the last stages of development as seen in our expression data. these predicted enzymes are thus strong candidates for controlling ginsenoside biosynthesis in the late stages of plant development. while this study, as any next-generation sequencing study, would have benefited from an even larger amount of sequence data, additional sequencing over so many stages of development was unfortunately cost prohibitive. nevertheless, the overall structure of the assembly with regard to number of genes, isoform frequency, length of transcripts and level of homology with existing est libraries and annotation patterns identified among the transcripts is all very supportive and indicative of a strong representation of the biology. overall, we believe the analysis benefited considerably from the use of the much longer reads that sequencing is capable of generating. this added information translates into more reliable, longer, and complete transcripts, as well as more information for improved accuracy in the calling of alternative splicing among sequenced transcripts. recent work by sun et al., found similar benefits to the use of sequencing in the assembly of a transcriptome for north american ginseng. our work has expanded on these initial efforts through the generation of over six times the amount of high quality reads sampled over seven stages of seasonal root development. this provides an unprecedented snapshot of north american ginseng root development and a means to comprehensively capture root expressed transcripts and their relative abundance through developmental time. in addition, we have also taken advantage of the recent trinity assembly algorithm, which is designed specifically for assembling transcriptomes. altogether, this provides the benefit of a high confidence, reliable assembly of putative transcripts that includes potentially valuable splice variant information. in the end, the combination of a fully annotated assembly along with corresponding expression data as sampled across developmental time has provided for an insightful look into the transcriptome of north american ginseng.  ginseng is a valuable medicinal herb that is widely used in traditional medicine and has considerable commercial value. this study has focused on the ginseng root as it is the tissue almost exclusively used for medicinal applications and contains the majority of the commercially and medicinally valuable ginsenosides. using extensive sampling of root tissue across seasonal development and next-generation sequencing we have generated a reliable north american ginseng transcriptome with extensive annotation and isoform information. expression profiling of all ginsenoside backbone biosynthesis genes suggests that ginsenoside biosynthesis is tightly linked with the developmental stages of flowering and fruit drop. we anticipate the assembly and associated expression data should prove invaluable to those conducting transcriptomic and metabolomic research against panax quinquefolius. 
annotation efforts using a number of public databases resulted in detailed annotation information for transcripts. ginseng including north american ginseng is one of the most widely used medicinal plants. rna extracts from root samples from seven development stages of north american ginseng were subjected to sequencing, filtered for quality and used in the de novo assembly of a collective root reference transcriptome consisting of transcripts. the assembly generated provides a comprehensive annotated reference for future transcriptomic study of north american ginseng. as the root is the primary tissue harvested commercially for ginsenosides, next generation sequencing was applied to the characterization and assembly of the root transcriptome throughout seasonal development. transcript expression profiles across seasonal development suggest a primary dammarane-type ginsenoside biosynthesis occurs just prior to plant senescence, with secondary ginsenoside production occurring throughout development. its success is thought to be due to a diverse collection of ginsenosides that serve as its major bioactive compounds. among our results, we found all of the known enzymes involved in the ginsenoside backbone biosynthesis and used co-expression analysis to identify a number of candidate sequences involved in the latter stages ginsenoside biosynthesis pathway. transcript profiles suggest ginsenoside biosynthesis occurs at distinct stages of development. in addition, genes were assigned to metabolic pathways using the kyoto encyclopedia of genes and genomes. a collection of putative ginsenoside biosynthesis genes were identified and candidate genes predicted from the lesser understood downstream stages of biosynthesis.
 ginseng including north american ginseng is one of the most widely used medicinal plants. its success is thought to be due to a diverse collection of ginsenosides that serve as its major bioactive compounds. however, few genomic resources exist and the details concerning its various biosynthetic pathways remain poorly understood. as the root is the primary tissue harvested commercially for ginsenosides, next generation sequencing was applied to the characterization and assembly of the root transcriptome throughout seasonal development. transcripts showing homology to ginsenoside biosynthesis enzymes were profiled in greater detail.  rna extracts from root samples from seven development stages of north american ginseng were subjected to sequencing, filtered for quality and used in the de novo assembly of a collective root reference transcriptome consisting of transcripts. annotation efforts using a number of public databases resulted in detailed annotation information for transcripts. in addition, genes were assigned to metabolic pathways using the kyoto encyclopedia of genes and genomes. among our results, we found all of the known enzymes involved in the ginsenoside backbone biosynthesis and used co-expression analysis to identify a number of candidate sequences involved in the latter stages ginsenoside biosynthesis pathway. transcript profiles suggest ginsenoside biosynthesis occurs at distinct stages of development.  the assembly generated provides a comprehensive annotated reference for future transcriptomic study of north american ginseng. a collection of putative ginsenoside biosynthesis genes were identified and candidate genes predicted from the lesser understood downstream stages of biosynthesis. transcript expression profiles across seasonal development suggest a primary dammarane-type ginsenoside biosynthesis occurs just prior to plant senescence, with secondary ginsenoside production occurring throughout development. data from the study provide a valuable resource for conducting future ginsenoside biosynthesis research in this important medicinal plant. north american ginsengtranscriptomenext generation sequencingginsenoside 