BACKGROUND
crocus sativus is a triploid, sterile geophyte which has been cultivated and used as spice and medicinal plant since thousands of years  <cit> .
it is vegetatively propagated by corms which maintains genetic characteristics of the plant but limits its improvement.
crocus belongs to family iridaceae, members of which have relatively large but poorly characterized genomes  <cit> .
crocus genus consists of about  <dig> species and many of them are considered as economically valuable.
the main crocus producing countries are iran, greece, spain, italy and india .
c. sativus has 2n = 3x =  <dig> chromosomes.
it is thought to be sterile triploid form of c. cartwrightianus  <cit> .
the desiccated stigma of c. sativus forms saffron and is source of various carotenoids and unique compounds called apocarotenoids which are produced by oxidative tailoring of carotenoids  <cit> .
in fact crocus is the only plant which produces apocarotenoids like crocin, picrocrocin and safranal in significant amounts.
these compounds provide color, flavour and aroma to this crop making it world’s costliest spice  <cit> .
in addition, these compounds have a broad spectrum of pharmacological properties as well  <cit> .
since crocus produces important carotenoids and apocarotenoids, it becomes imperative to have a holistic approach for identification and isolation of genes involved in their metabolic pathways.
moreover, carotenoid biosynthesis and degradation is thought to be tightly regulated throughout the life cycle of a plant and dynamic changes appear to occur in their composition so as to stay in tune with prevailing developmental requirements and environmental conditions  <cit> .
also their regulation is relatively a complex phenomenon  <cit>  with cross talk and integration of various pathways at multiple levels so as to achieve metabolic flexibility and robustness in response to environmental signals.
therefore, transcriptome sequencing would pave way for elucidation of metabolic pathway in crocus and corresponding regulatory networks.
further, the genome organization of iridaceae family is not yet known and therefore transcriptome sequencing of crocus would be first step to provide a gene atlas for other members of this familyas well.
the characterization of saffron transcriptome is in fact a prerequisite to shed light on essential biological processes including the molecular basis of flavour and color biogenesis, genomic organization and flower development of c. sativus in particular and family iridaceae in general.
although, a few studies have been carried out on the c. sativus plant to understand the flower development and apocarotenoid biosynthesis  <cit> , molecular basis of these essential processes is still largely unknown.
this study presents first report on transcriptome sequencing and analysis of crocus stigma and rest of the flower  using high-throughput illumina sequencing.
the assembled transcripts were annotated to identify genes involved in various steps of the carotenoid/apocarotenoid pathway.
the transcriptome aids in significant addition to the number of genes which are potentially involved in apocarotenoid biosynthesis.
the transcriptome was further screened for the identification of the genes encoding transcription factors  so as to generate transcription factor database from crocus which may help in unravelling the regulatory mechanism of apocarotenoid biosynthesis.
finally expression specificity of the assembled transcripts provided an accurate estimation of the biological processes involved in different tissues .
the crocus transcriptome provides a platform for understanding the molecular basis of carotenoid/apocarotenoid pathway and various other biological processes pertaining to this crop.
methods
tissue sampling, cdna library construction and sequencing
crocus sativus l. tissue used in this study was collected from experimental farm at indian institute of integrative medicine, srinagar, jammu and kashmir, india .
the voucher specimen was deposited at janaki ammal herbarium , iiim, jammu.
the details of the specimen are: .
the most important crocus apocarotenoids which include crocin, picrocrocin and saffranal are synthesized and accumulated in stigmas and their amount increases from yellow, through orange to scarlet stage.
also the compounds increase in quantity from pre- anthesis stage to anthesis and then follow a decline post anthesis.
in view of this, we collected flowers at scarlet stage.
further, three flowers each were collected at pre-anthesis, anthesis and one day after anthesis.
stigmas were hand-picked from flowers.
the tissue samples were frozen in liquid nitrogen and stored at -80 °c for further use.
total rna was extracted from stigma and rest of the flower using trizol .
the quality and quantity of total rna was analysed using nanodrop spectrophotometer and their integrity was further evaluated using bioanalyzer .
high quality rna isolated from three independent tissue samples  was pooled for library preparation.
the cdna libraries of stigma and rest of the flower tissue were constructed using illumina truseq rna preparation kit v <dig>  following manufacturer’s instructions.
the libraries were quantified using qubit dsdna br assay kit .
the insert size of libraries was verified using bioanalyzer dna  <dig> chip.
further for generating clusters, 10pm of each library was loaded onto the flow cell using truseq pe cluster kit v <dig> on cluster station .
the flow cell containing clonally amplified clusters was loaded onto the genome analyzer iix  and paired-end   sequencing was performed.
the sequence reads which are obtained from sequencer often contain adapter sequences, low quality reads and very short length reads.
this data is therefore processed in order to remove such reads.
paired reads were quality filtered using ngs qc toolkit v  <dig>   <cit> .
the cutoff for quality score is >  <dig> q <dig> score and should have high quality bases > 70 % of read length.
high quality reads were used for de novo assembly using trinity software with k-mer of  <dig>  the assembly resulted in contigs and singletons which together form set of unigenes.
annotation of the assembled transcriptome
crocus transcriptome was annotated using blastx similarity search against ncbi-nr database.
homology search was also made against other databases like swiss-prot and uniprot trembl databases.
further go term and interpro domain annotation for the assembled transcripts was performed using the trapid annotation server   <cit> .
kegg orthology assignments for the transcripts were obtained using the kaas server   <cit> .
for identification of transcription factors, homology search was done against planttfdb  <cit> .
in all cases analysis was performed using the default parameters.
differential expression analysis
the assembled transcripts were filtered to omit sequences <  <dig> nucleotides.
further trimmed mean of m values  normalization was performed on the raw reads for each transcript using the noiseq package.
then noiseq count filter was used to remove transcripts with counts per million  <  <dig> in both samples which avoids noise from lowly expressed transcripts.
further, the filtered transcripts were analysed with the bioconductor noiseq package  <cit>   to identify those differentially expressed between flower and stigma.
the gene ontology enrichment analysis was performed on the go mapping done by trapid server step using a custom r script to select significantly enriched go classes in the differentially expressed transcripts compared to the total transcriptome .
finally custom perl and r scripts were used to associate differentially expressed transcripts with their kegg identifiers and convert the expression information into a format suitable to be visualized in ipath <dig>  <cit> .
quantitative real time pcr
for real time pcr, total rna was extracted using trizol reagent and used for cdna synthesis by reverse transcription kit  following manufacturer’s instructions.
qrt-pcr was performed in triplicates in abi stepone real time .
the reaction was carried out in a total volume of 20 μl which consisted of 10 μl of 2x sybr green master mix,  <dig>  μm gene specific forward and reverse primers and 100 ng of template cdna.
the cycling parameters were 95 °c for 20 s, followed by 40 cycles of 95 °c for 15 s and 60 °c for 1 min.
the sequence of the primers used in this study is given in additional file  <dig>  the relative quantification method  was used to evaluate quantitative variation between the replicates examined.
the amplification of 18s rrna gene was used as an endogenous control to normalize the data.
sequence analysis and phylogenetics
the deduced amino acid sequences of the ccd proteins were aligned by a multiple sequence comparison using the log-expectation  alignment tool  with the default parameters  <cit> .
the phylogenetic analysis was done using the neighbor-joining method and  <dig> bootstrap replicates were employed in each analysis to maximize the statistical significance  <cit> .
the phylogenetic trees were constructed and visualized by mega <dig>  software  <cit> .
the accession numbers of the genes used for phylogenetic analysis are given in additional file  <dig>  all the phylogenetic data including the sequence files, alignment and phylogenetic tree have been submitted to dryad  .
RESULTS
crocus transcriptome sequencing and de novo assembly
the sequencing of crocus cdna libraries generated  <dig> , <dig> raw reads from flower and  <dig> , <dig> from stigma.
approximately 14 % of raw reads were removed post filtering of adapter sequences, low quality and short reads.
further,  <dig> , <dig> flower and  <dig> , <dig> stigma reads were assembled de novo using trinity which resulted in  <dig>  transcripts.
these transcripts were classified into  <dig>  unigenes comprised of  <dig> clusters and  <dig>  singletons.
the statistical summary of data is outlined in table  <dig>  the average contig length was  <dig>  bp, gc content  <dig>  % and n <dig> was 753 bp.table  <dig> summary of crocus sativus transcriptome
functional annotation and classification of crocus transcriptome
for comprehensive annotation of crocus transcriptome, similarity search for the sequences was done using blastx against nr database in ncbi, with e- value cutoff of 10− <dig>  interpro, swissprot and uniprot trembl analysis were also used for annotation.
among the  <dig>  transcripts  <dig>  % could be annotated.
the putative functions assigned to crocus transcripts are available as additional file  <dig>  crocus belongs to family iridaceae and whole genome and/or transcriptome of none of the members of this family has been sequenced so far and that may be one of the reasons for relatively low homology results.
although est database of crocus was developed earlier but that reports only  <dig> unique transcripts  <cit> .
the sequences with unknown homology may represent genes involved in metabolic processes which are unique to this plant and whose intermediates and enzymes have not been identified so far.
moreover, many transcriptome studies of other plant species have also reported functional annotation of around half or even less percent of unigenes.
for example, in case of transcriptome analysis of cymbidium sinence, only  <dig>  % of the unigenes could be annotated  <cit> .
further, significant match could be assigned to only 55 % of the unigenes in bamboo transcriptome  <cit>  while in case of amaranthus tricolor only  <dig>  % of the unigenes showed significant similarity  <cit> .
the crocus unigenes were further classified according to gene ontology annotation into three categories viz biological process, cellular component and molecular function .
the top  <dig> classes from each category are shown in fig.
<dig>  in the category of biological process most of the genes belonged to metabolic process class  followed by cellular process , primary metabolic process  and cellular metabolic process .
in case of cellular component, the predominant categories were cell  and cell part .
as far as molecular function is concerned, binding  and catalytic activity  were the major classes.
the go term abundance results show similarity with previous transcriptome studies, for example, in gardenia  <cit> , safflower  <cit>  and chrysanthemum  <cit> .fig.
<dig> go classification of crocus unigenes.
bar chart represents percent genes belonging to top  <dig> biological process, cellular component and molecular function categories
in order to elucidate the biochemical pathways represented in crocus transciptome, the sequences were searched against kegg pathway database.
a total of  <dig> unigenes  were assigned to  <dig> kegg pathways .
the most abundant pathways were ribosome comprising of  <dig> genes followed by carbon metabolism , starch and sucrose metabolism  and biosynthesis of amino acids .
further, protein processing , spliceosome  and oxidative phosphorylation  also represented significantly higher number of genes.
we also identified classes like plant hormone signal transduction , plant pathogen interaction  and ubiquitin mediated proteolysis .
the representative top  <dig> classes are depicted in fig.
<dig>  in many other transcriptome studies also, the above classes represented the predominant categories.
for example, in horse gram the highest number of genes was mapped to ribosome biosynthesis  <cit>  while in chickpea the predominant classes were ribosome, spliceosome and biosynthesis of amino acids  <cit> .
these pathways and the genes involved thereof might be involved in growth and developmental pathways and also in plant response to various environmental cues.fig.
<dig> kegg mapping of crocus unigenes representing number of genes mapped to top  <dig> classes
identification of genes involved in carotenoid/apocarotenoid biosynthesis
carotenoid biosynthesis occurs through mep pathway  <cit> .
the genes encoding enzymes which are involved in this pathway have been isolated in many plant species  <cit> .
in case of crocus also, the pathway has been elucidated to a large extent  <cit> , however, there are still gaps which remain to be filled in terms of the enzymes and the intermediates involved in the pathway.
in this study we report identification of many carotenoid/apocarotenoid pathway genes which were known from other plants but not isolated from crocus.
this would provide a knowledgebase for understanding the biosynthesis of carotenoids and their subsequent degradation to apocarotenoids in this plant.
core carotenoid biosynthetic pathway is initiated by condensation of two molecules of geranylgeranyl diphosphate  to form phytoene .
this step is catalysed by phytoene synthase  enzyme.
this is the rate limiting step of this pathway  <cit> .
phytoene, is then desaturated into lycopene by two related enzymes of phytoene desaturase  and ζ-carotene desaturase  in plants  <cit> .
this lycopene is acted upon by two cis-trans isomerases of z-iso  <cit>  and crtiso  <cit>  and is converted from poly-cis-configured phytoene into the all-trans form lycopene.
while psy and pds were already known from crocus  <cit> , we have identified zds, z-iso and crtiso for the first time.
lycopene represents the branching point of this pathway and is cyclized either to form α-carotene by the action of lycopene ε-cyclase  and lycopene β-cyclase  or is converted to β-carotene by β-lcy alone  <cit> .
β-lcy has been identified from crocus earlier  <cit> , however, ε-lcy from crocus was not known and has been identified in this study.
futher, α- and β-carotene are hydroxylated to produce lutein and zeaxanthin, respectively by the action of β-ring carotene hydroxylase .
many isoforms of bch have been identified from crocus earlier  <cit> .
there is a class of enzymes called carotenoid cleavage dioxygenases  which cleave double bonds of carotenoids at different positions resulting in the formation of apocarotenoids.
zeaxanthin acts as substrate for ccd <dig> which cleaves it into crocetin dialdehyde and β-cyclocitral.
ccd <dig> was very recently identified to be the enzyme which has 7′ <dig> cleavage activity and catalyses the first step in crocin biosynthesis in crocus  <cit> .
crocetin dialdehyde is converted into crocetin by aldehyde dehydrogenase.
in this study we have identified many transcripts encoding aldehyde dehydrogenases which might include the one involved in converting crocetin dialdehyde into crocetin.
further, crocetin and β-cyclocitral are converted into crocin and picrocrocin respectively by udp-glucosyl transferases.
many transcripts coding for udp-glucosyl transferases have been identified in this study.
so far only a few udp-glucosyl transferases were known from crocus  <cit>  and we have added new isoforms to the existing database.
the picrocrocin is further converted into saffranal which is supposed to be a beta glucosidase action.
gene encoding this enzyme has not been identified in crocus, so far.
here we have identified many transcripts coding for this enzyme.
therefore, the present study would provide a platform for generating knowledge about the substrate specificities and activities of the enzymes identified so as to understand the apocarotenoid biosynthetic pathway.
further, there are many other ccd isoforms which act on a wide range of substrates, cleave them at different positions and produce a myriad of apocarotenoid products  <cit> .
analysing the enzyme activities and their substrate specificities would help to identify new compounds.
this knowledge can further be useful for designing metabolic engineering strategies for either enhanced production of known compounds or production of new metabolites.fig.
<dig> schematic representation of apocarotenoid pathway showing different enzymes and intermediates involved in the pathway.
genes marked in red are identified first time in this study
considering other branch points of this pathway, zeaxanthin is also converted into violaxanthin by zeaxanthin epoxidase  which is in turn converted into neoxanthin by neoxanthin synthase .
violaxanthin and neoxanthin are cleaved by 9-cis-epoxycarotenoid dioxygenase  to produce xanthoxin, the direct substrate for phytohormone abscisic acid  synthesis.
while nced from crocus is known  <cit> , the genes encoding the enzymes which act upstream  are identified for the first time in this study.
the schematic representation of the carotenoid/apocarotenoid pathway is given in fig.
<dig> and the list of pathway genes identified is given in table  <dig> table  <dig> genes involved in carotenoid/apocarotenoid pathway
agenes identified first time in the present study
identification of transcription factors
in order to identify transcription factor encoding genes in crocus transcriptome, homology search was performed against the plant transcription factor database .
around  <dig> unigenes  show similarity to transcription factors representing  <dig> families .
out of these, transcription factors belonging to myb family were most abundant  followed by c3h , far <dig> , mads box  .
in many other studies also, the above mentioned classes represented highly expressed transcription factor families.
for example, in catharanthus roseus, myb was the most abundant family  <cit> .
in horse gram, most highly represented transcription factors belonged to c3h, bhlh and ap <dig> domain families  <cit> .
transcription factors perform key roles in plant growth and development.
they are also involved in regulation of secondary metabolism and also coordinate plant’s response to various environmental cues.
several members of myb family have been shown to regulate different secondary metabolic pathways .
zinc finger proteins of different classes are also involved in regulation of plant secondary metabolism  <cit> .
mads box genes regulate a range of plant processes including flower development  <cit>  plant reproduction  <cit>  etc.
so far regulatory pathway controlling crocus apocarotenoid metabolism is not known.
although many transcription factors were identified from crocus  <cit> , their role in apocarotenoid regulation was not experimentally validated.
recently, a sand domain ultrapetala transcription factor was isolated from our laboratory and was shown to regulate apocarotenoid biosynthesis in crocus  <cit> .
no other information is available as far regulation of apocarotenoid biosynthesis is concerned.
this transcription factor database would therefore, be an important asset to characterize the regulatory pathways of crocus secondary metabolism.
this may also provide base for identification of regulators of crocus flower development, its tripartite stigma and would help in addressing other question related to biology of this plant.fig.
<dig> percent crocus transcripts representing top  <dig>  transcription factor families and  conserved domain families
domain analysis
we aimed to get insight about the molecular mechanism underlying the phenomenon of colour and flavour biogenesis in crocus.
go classification showed the major classes in molecular function comprised of binding, catalytic activity and nucleotide binding.
further, kegg mapping indicated that genes involved in protein processing, spliceosome and oxidative phosphorylation were significantly represented.
in order to extrapolate this information and build up knowledge about the mechanism of action of the genes involved in such activities we performed conserved domain analysis for crocus transcriptome.
the major domain classes are given in fig.
4b.
the results showed that most of the genes were those with protein kinase like domain  followed by protein kinase catalytic domain , serine/threonine protein kinase domain , serine/threonine-protein kinase-like domain  and tyrosine-protein kinase catalytic domain  .
these domains are involved in phosphorylation of proteins at various amino acid residues.
protein modification by phosphorylation and dephosphorylation is a crucial mechanism that controls activity of proteins  <cit>  and as a result regulates important cellular functions in eukaryotes such as cell division, differentiation, signal transduction etc.
predominant presence of genes with different kinase domains in crocus transcriptome is indicative of the fact that even though transcriptional regulation is the central regulatory mechanism for most of the biological processes, diverse post translational events also play important role.
these multiple layers of regulation help plant to synchronize properly with developmental stages and environmental cues.
differential gene expression, gene ontology and pathway enrichment
in crocus, the stigma part of the flower has commercial significance because its dried state forms saffron which is the site of biosynthesis and accumulation of apocarotenoids.
although attempts have been made to identify and isolate genes involved in apocarotenoid biosynthesis, little has been done so far.
in order to identify genes involved in biosynthesis and regulation of apocarotenoids in crocus stigma, we studied differential gene expression  in crocus stigma vs. rest of the flower.
a total of  <dig> transcripts  were differentially expressed out of which  <dig> transcripts  were upregulated in flower and  <dig> transcripts  in stigma .
in order to identify the major functional categories represented by differentially expressed genes, go enrichment analysis was carried out .
we observed that in flower upregulated genes, the top five classes in biological process were metabolic process, cellular process, primary metabolic process, macromolecule metabolic process and protein metabolic process .
in case of the genes upregulated in stigma, metabolic process, cellular process, primary metabolic process, cellular metabolic process and macromolecule metabolic process were the top five categories .
in case of the molecular function category, the top five classes represented by flower upregulated genes were catalytic activity, hydrolase activity, ion binding, cation binding and metal ion binding, while as in case of stigma top five classes were binding, catalytic activity, nucleotide binding, purine nucleotide binding and ribonucleotide binding.
further, genes related to transport were also highly represented in flower while in stigma genes related to nucleotide activity were more prevalent.
we also performed kegg analyses on differentially expressed genes in flower and stigma.
the results showed that in case of flower, genes involved in ribosome, starch and sucrose metabolism, plant-pathogen interaction, phenylpropanoid biosynthesis and pentose and glucuronate interconversions were enriched  while in case of stigma, genes related to protein processing in endoplasmic reticulum, photosynthesis, carotenoid biosynthesis, legionellosis represented the major classes  .
enrichment of carotenoid biosynthesis genes in stigma confirms the fact that biosynthesis of carotenoids and their subsequent degradation into apocarotenoids occurs mainly in stigma.fig.
<dig> differential expression analysis of crocus transciptome.
a, b bar chart showing go terms associated with genes differentially expressed in flower and stigma.
c, d pie chart representing kegg classes of flower and stigma upregulated genes
we also investigated the expression of genes involved in carotenoid/apocarotenoid biosynthesis and observed that most of them were more significantly enriched in stigma  which is in accordance with the fact that stigma is the actual site of biosynthesis of apocarotenoids in crocus.
for example, psy and pds which catalyse initial steps in carotenoid biosynthetic pathway are expressed more in stigma.
further, ziso and zcd which are involved in converting phytoene into lycopene are also upregulated in stigma.
it was quite interesting to see that lycopene β-lcy and bch were upregulated in stigma therefore increasing the metabolic flux towards production of zeaxanthin.
earlier reports have also shown increased expression of these genes in crocus stigma  <cit> .
recently it was shown that ccd <dig> cleaves zeaxanthin at 7′ <dig> double bond and results in the formation of crocetin and picrocrocin which are apocarotenoids responsible for color and flavor of saffron  <cit> .
dge analysis showed that ccd <dig> is also upregulated in stigma thereby confirming the earlier reports.
further, another isoform ccd4b is also expressed more in stigma.
this enzyme cleaves carotenoids like beta-carotene at 9′ <dig> double bond and forms apocarotenoids like β-cyclocitral and β-ionone.
earlier reports have also shown enhanced expression of ccd4b in stigma  <cit> .
thus our results are in agreement with earlier reports and also provide molecular proof for stigma being the actual site of biosynthesis of apocarotenoids which includes crocin, picrorocin, β-cyclocitral and β-ionone.fig.
<dig> heat-map showing expression patterns of  regulatory  carotenoid pathway genes differentially upregulated in stigma and flower
one of the aims of the present work was to identify transcription factors which regulate the biosynthesis and accumulation of apocarotenoids in crocus in tissue and developmental stage specific manner.
our data suggested that transcription factors belonging to various families like mads box, myb, zinc finger, wrky, phd etc.
were differentially expressed in both flower and stigma .
we observed that different transcription factor families were enriched in stigma and flower.
in this context we observed that transcription factors belonging to mads box were represented more in stigma while nac, bhlh and wrky families were enriched more in flower suggesting that they perform specialized roles in different tissues.
there were many other transcription factors exhibiting expression in both the tissue types.
for example, zinc finger proteins with an <dig> domain were more prevalent in stigma upregulated genes while those with zat domains were upregulated in flower.
this indicates that different members of the same family might play different roles in different tissues or under different conditions.
experimental validation of differential expressed genes by quantitative realtime pcr
qrt-pcr of ten selected genes was performed in order to validate the differential gene expression obtained by rna-seq.
the genes selected were ccd4b, ccd <dig>  bch, psy, pds, gt, zinc-finger, mads box, myb and, wrky.
the results indicated that expression pattern as obtained by qrt-pcr corroborated with that obtained by rna-seq for all the genes .
this confirmed the reliability of rna-seq data.
we also performed statistical analysis on the data obtained for these genes from rna-seq and qpcr and observed that there was very good correlation between the two  .fig.
<dig> qrt-pcr validation of differentially expressed genes.
a expression of  <dig> genes was validated using qrt-pcr and compared with their expression obtained from rna-seq.
qrt-pcr was performed and values were normalized relative to the expression levels of 18s rrna in the same cdna sample.
data are the means  of three replicates expression data are presented as expression values of genes in stigma sample relative to their expression in flower sample  correlation of gene expression results obtained from qpcr and rna-seq
phylogenetic analyses of ccd gene family
carotenoid cleavage dioxygenases  form a group of enzymes which are involved in cleavage of carotenoids leading to the production of apocarotenoids.
ccds are specific to the double bond they cleave but often exhibit substrate promiscuity which is responsible for the diversity of apocarotenoids found in nature.
apocarotenoids play role in many aspects of plant growth and development.
the ccd family is ancient and has its members present in bacteria, plants and animals.
members of the ccd family share several characteristics like, they require a fe2+ for catalytic activity; they have four conserved histidines which coordinate iron binding and they contain a conserved peptide sequence at their carboxyl terminus  <cit> .
in order to get an insight about the phylogenetic relationship between various ccds, a neighbour-joining phylogenetic tree was constructed with  <dig> ccd genes from  <dig> plant species.
the amino acid sequences from various ccds were grouped into five clusters, named ccd <dig>  ccd <dig>  ccd <dig>  ccd <dig> and nced .
in each of these clusters, ccds were present in two groups corresponding to monocotyledon and dicotyledon species.
further, genes or isoforms within the ccd sub-family and belonging to same species were grouped in the same branch, for example, crocus sativus ccd4a/b, ccd8a/b.
phylogenetic analysis also revealed that two major duplications had occurred in ccd subfamilies.
ccd duplication, which occurred in the moss lineages, ultimately led to the emergence of two lineages that evolved into ccd <dig> and ccd <dig>  this result suggests that, ccd7/ <dig> genes had similar evolutionary trends than ccd <dig> and ccd <dig> sub-family.
in general, ccds belonging to a particular cluster have similar cleavage activity.
recently a new ccd isoform  was identified from crocus sativus, which plays an important role in crocin biosynthesis.
csccd <dig> clusters with ccd <dig> subfamily but is distinct from ccd <dig> subfamily as far as its catalytic activity is concerned  <cit> .
this suggests that ccd <dig> might have evolved from ccd <dig> and developed different cleavage site specificity.
thus phylogenetic analysis may help in understandin g functional diversity of ccd gene family.fig.
<dig> phylogenetic tree for ccd proteins.
tree was constructed by aligning  <dig> ccd proteins from  <dig> plant species including crocus using the muscle program and subjected to phylogenetic ana lysis by the distance with neighbor-joining method using the mega <dig>  software.
the numbers on the nodes indicate the bootstrap values after  <dig> replicates
CONCLUSIONS
in the present study, crocus sativus transcriptome was sequenced, assembled and annotated.
the database led to the identification of many new candidate genes involved in carotenoid/apocarotenoid pathway.
identification of transcription factors provides a platform for unravelling the regulatory pathway of crocus flower development and apocarotenoid biosynthesis.
differential gene expression and pathway mapping confirmed enrichment of apocarotenoid pathway genes in stigma thus confirming that stigma is the site of apocarotenoid biosynthesis.
the transcript resource generated would therefore facilitate and enhance our understanding of biosynthetic pathway of carotenoids and their subsequent tailoring into apocarotenoids and the mechanism that regulates carotenoid-apocarotenoid metabolic flux.
availability of supporting data
the illumina sequence data have been submitted as bioproject  to ncbi sequence read archive under accession number .
all the other supporting data are included as additional files.
additional files
additional file 1:
sequence of the primers used in this study.
additional file 2:
accession numbers of the genes used for phylogenetic analysis.
additional file 3:
functional annotation of
crocus
transcriptome.
additional file 4:
go and kegg classes of
crocus
transcripts.
additional file 5:
details of transcription factor and conserved domain families.
additional file 6:
details of differentially expressed transcripts and flower and stigma upregulated genes.
additional file 7:
go and kegg classes in flower upregulated genes.
additional file 8:
go and kegg classes in stigma upregulated genes.
abbreviations
gogene ontology
keggkyoto enclyclopedia of genes and genomes
ccdcarotenoid cleavage dioxygenase
ggdpgeranylgeranyl diphosphate
psyphytoene synthase
pdsphytoene desaturase
shoib ahmad baba and tabasum mohiuddin share first authorship.
swaraj basu and mohit kumar swarnkar share second authorship.
shoib ahmad baba tabasum mohiuddin swaraj basu and mohit kumar swarnkar contributed equally to this work.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
na and aks: conceptualized and designed experiments and wrote manuscript; sab, tm, ahm, zaw and na performed experiments; sb and mks did analysis.
all authors read and approved the final manuscript.

 crocus sativus is a triploid, sterile geophyte which has been cultivated and used as spice and medicinal plant since thousands of years. it is vegetatively propagated by corms which maintains genetic characteristics of the plant but limits its improvement. crocus belongs to family iridaceae, members of which have relatively large but poorly characterized genomes. crocus genus consists of about species and many of them are considered as economically valuable. the main crocus producing countries are iran, greece, spain, italy and india. c. sativus has 2n = 3x =  chromosomes. it is thought to be sterile triploid form of c. cartwrightianus. the desiccated stigma of c. sativus forms saffron and is source of various carotenoids and unique compounds called apocarotenoids which are produced by oxidative tailoring of carotenoids. in fact crocus is the only plant which produces apocarotenoids like crocin, picrocrocin and safranal in significant amounts. these compounds provide color, flavour and aroma to this crop making it world’s costliest spice. in addition, these compounds have a broad spectrum of pharmacological properties as well. since crocus produces important carotenoids and apocarotenoids, it becomes imperative to have a holistic approach for identification and isolation of genes involved in their metabolic pathways. moreover, carotenoid biosynthesis and degradation is thought to be tightly regulated throughout the life cycle of a plant and dynamic changes appear to occur in their composition so as to stay in tune with prevailing developmental requirements and environmental conditions. also their regulation is relatively a complex phenomenon with cross talk and integration of various pathways at multiple levels so as to achieve metabolic flexibility and robustness in response to environmental signals. therefore, transcriptome sequencing would pave way for elucidation of metabolic pathway in crocus and corresponding regulatory networks. further, the genome organization of iridaceae family is not yet known and therefore transcriptome sequencing of crocus would be first step to provide a gene atlas for other members of this familyas well. the characterization of saffron transcriptome is in fact a prerequisite to shed light on essential biological processes including the molecular basis of flavour and color biogenesis, genomic organization and flower development of c. sativus in particular and family iridaceae in general. although, a few studies have been carried out on the c. sativus plant to understand the flower development and apocarotenoid biosynthesis, molecular basis of these essential processes is still largely unknown. this study presents first report on transcriptome sequencing and analysis of crocus stigma and rest of the flower using high-throughput illumina sequencing. the assembled transcripts were annotated to identify genes involved in various steps of the carotenoid/apocarotenoid pathway. the transcriptome aids in significant addition to the number of genes which are potentially involved in apocarotenoid biosynthesis. the transcriptome was further screened for the identification of the genes encoding transcription factors so as to generate transcription factor database from crocus which may help in unravelling the regulatory mechanism of apocarotenoid biosynthesis. finally expression specificity of the assembled transcripts provided an accurate estimation of the biological processes involved in different tissues. the crocus transcriptome provides a platform for understanding the molecular basis of carotenoid/apocarotenoid pathway and various other biological processes pertaining to this crop. methods tissue sampling, cdna library construction and sequencing crocus sativus l. tissue used in this study was collected from experimental farm at indian institute of integrative medicine, srinagar, jammu and kashmir, india. the voucher specimen was deposited at janaki ammal herbarium, iiim, jammu. the details of the specimen are:. the most important crocus apocarotenoids which include crocin, picrocrocin and saffranal are synthesized and accumulated in stigmas and their amount increases from yellow, through orange to scarlet stage. also the compounds increase in quantity from pre- anthesis stage to anthesis and then follow a decline post anthesis. in view of this, we collected flowers at scarlet stage. further, three flowers each were collected at pre-anthesis, anthesis and one day after anthesis. stigmas were hand-picked from flowers. the tissue samples were frozen in liquid nitrogen and stored at -80 °c for further use. total rna was extracted from stigma and rest of the flower using trizol. the quality and quantity of total rna was analysed using nanodrop spectrophotometer and their integrity was further evaluated using bioanalyzer. high quality rna isolated from three independent tissue samples was pooled for library preparation. the cdna libraries of stigma and rest of the flower tissue were constructed using illumina truseq rna preparation kit v following manufacturer’s instructions. the libraries were quantified using qubit dsdna br assay kit. the insert size of libraries was verified using bioanalyzer dna chip. further for generating clusters, 10pm of each library was loaded onto the flow cell using truseq pe cluster kit v on cluster station. the flow cell containing clonally amplified clusters was loaded onto the genome analyzer iix and paired-end sequencing was performed. the sequence reads which are obtained from sequencer often contain adapter sequences, low quality reads and very short length reads. this data is therefore processed in order to remove such reads. paired reads were quality filtered using ngs qc toolkit v. the cutoff for quality score is >  q score and should have high quality bases > 70 % of read length. high quality reads were used for de novo assembly using trinity software with k-mer of the assembly resulted in contigs and singletons which together form set of unigenes. annotation of the assembled transcriptome crocus transcriptome was annotated using blastx similarity search against ncbi-nr database. homology search was also made against other databases like swiss-prot and uniprot trembl databases. further go term and interpro domain annotation for the assembled transcripts was performed using the trapid annotation server. kegg orthology assignments for the transcripts were obtained using the kaas server. for identification of transcription factors, homology search was done against planttfdb. in all cases analysis was performed using the default parameters. differential expression analysis the assembled transcripts were filtered to omit sequences <  nucleotides. further trimmed mean of m values normalization was performed on the raw reads for each transcript using the noiseq package. then noiseq count filter was used to remove transcripts with counts per million  <  in both samples which avoids noise from lowly expressed transcripts. further, the filtered transcripts were analysed with the bioconductor noiseq package to identify those differentially expressed between flower and stigma. the gene ontology enrichment analysis was performed on the go mapping done by trapid server step using a custom r script to select significantly enriched go classes in the differentially expressed transcripts compared to the total transcriptome. finally custom perl and r scripts were used to associate differentially expressed transcripts with their kegg identifiers and convert the expression information into a format suitable to be visualized in ipath. quantitative real time pcr for real time pcr, total rna was extracted using trizol reagent and used for cdna synthesis by reverse transcription kit following manufacturer’s instructions. qrt-pcr was performed in triplicates in abi stepone real time. the reaction was carried out in a total volume of 20 μl which consisted of 10 μl of 2x sybr green master mix,  μm gene specific forward and reverse primers and 100 ng of template cdna. the cycling parameters were 95 °c for 20 s, followed by 40 cycles of 95 °c for 15 s and 60 °c for 1 min. the sequence of the primers used in this study is given in additional file the relative quantification method was used to evaluate quantitative variation between the replicates examined. the amplification of 18s rrna gene was used as an endogenous control to normalize the data. sequence analysis and phylogenetics the deduced amino acid sequences of the ccd proteins were aligned by a multiple sequence comparison using the log-expectation alignment tool with the default parameters. the phylogenetic analysis was done using the neighbor-joining method and bootstrap replicates were employed in each analysis to maximize the statistical significance. the phylogenetic trees were constructed and visualized by mega software. the accession numbers of the genes used for phylogenetic analysis are given in additional file all the phylogenetic data including the sequence files, alignment and phylogenetic tree have been submitted to dryad.  crocus transcriptome sequencing and de novo assembly the sequencing of crocus cdna libraries generated, raw reads from flower and, from stigma. approximately 14 % of raw reads were removed post filtering of adapter sequences, low quality and short reads. further,, flower and, stigma reads were assembled de novo using trinity which resulted in transcripts. these transcripts were classified into unigenes comprised of clusters and singletons. the statistical summary of data is outlined in table  the average contig length was  bp, gc content  % and n was 753 bp.table summary of crocus sativus transcriptome functional annotation and classification of crocus transcriptome for comprehensive annotation of crocus transcriptome, similarity search for the sequences was done using blastx against nr database in ncbi, with e- value cutoff of 10− interpro, swissprot and uniprot trembl analysis were also used for annotation. among the transcripts  % could be annotated. the putative functions assigned to crocus transcripts are available as additional file crocus belongs to family iridaceae and whole genome and/or transcriptome of none of the members of this family has been sequenced so far and that may be one of the reasons for relatively low homology results. although est database of crocus was developed earlier but that reports only unique transcripts. the sequences with unknown homology may represent genes involved in metabolic processes which are unique to this plant and whose intermediates and enzymes have not been identified so far. moreover, many transcriptome studies of other plant species have also reported functional annotation of around half or even less percent of unigenes. for example, in case of transcriptome analysis of cymbidium sinence, only  % of the unigenes could be annotated. further, significant match could be assigned to only 55 % of the unigenes in bamboo transcriptome while in case of amaranthus tricolor only  % of the unigenes showed significant similarity. the crocus unigenes were further classified according to gene ontology annotation into three categories viz biological process, cellular component and molecular function. the top classes from each category are shown in fig.  in the category of biological process most of the genes belonged to metabolic process class followed by cellular process, primary metabolic process and cellular metabolic process. in case of cellular component, the predominant categories were cell and cell part. as far as molecular function is concerned, binding and catalytic activity were the major classes. the go term abundance results show similarity with previous transcriptome studies, for example, in gardenia, safflower and chrysanthemum.fig. go classification of crocus unigenes. bar chart represents percent genes belonging to top biological process, cellular component and molecular function categories in order to elucidate the biochemical pathways represented in crocus transciptome, the sequences were searched against kegg pathway database. a total of unigenes were assigned to kegg pathways. the most abundant pathways were ribosome comprising of genes followed by carbon metabolism, starch and sucrose metabolism and biosynthesis of amino acids. further, protein processing, spliceosome and oxidative phosphorylation also represented significantly higher number of genes. we also identified classes like plant hormone signal transduction, plant pathogen interaction and ubiquitin mediated proteolysis. the representative top classes are depicted in fig.  in many other transcriptome studies also, the above classes represented the predominant categories. for example, in horse gram the highest number of genes was mapped to ribosome biosynthesis while in chickpea the predominant classes were ribosome, spliceosome and biosynthesis of amino acids. these pathways and the genes involved thereof might be involved in growth and developmental pathways and also in plant response to various environmental cues.fig. kegg mapping of crocus unigenes representing number of genes mapped to top classes identification of genes involved in carotenoid/apocarotenoid biosynthesis carotenoid biosynthesis occurs through mep pathway. the genes encoding enzymes which are involved in this pathway have been isolated in many plant species. in case of crocus also, the pathway has been elucidated to a large extent, however, there are still gaps which remain to be filled in terms of the enzymes and the intermediates involved in the pathway. in this study we report identification of many carotenoid/apocarotenoid pathway genes which were known from other plants but not isolated from crocus. this would provide a knowledgebase for understanding the biosynthesis of carotenoids and their subsequent degradation to apocarotenoids in this plant. core carotenoid biosynthetic pathway is initiated by condensation of two molecules of geranylgeranyl diphosphate to form phytoene. this step is catalysed by phytoene synthase enzyme. this is the rate limiting step of this pathway. phytoene, is then desaturated into lycopene by two related enzymes of phytoene desaturase and ζ-carotene desaturase in plants. this lycopene is acted upon by two cis-trans isomerases of z-iso and crtiso and is converted from poly-cis-configured phytoene into the all-trans form lycopene. while psy and pds were already known from crocus, we have identified zds, z-iso and crtiso for the first time. lycopene represents the branching point of this pathway and is cyclized either to form α-carotene by the action of lycopene ε-cyclase and lycopene β-cyclase or is converted to β-carotene by β-lcy alone. β-lcy has been identified from crocus earlier, however, ε-lcy from crocus was not known and has been identified in this study. futher, α- and β-carotene are hydroxylated to produce lutein and zeaxanthin, respectively by the action of β-ring carotene hydroxylase. many isoforms of bch have been identified from crocus earlier. there is a class of enzymes called carotenoid cleavage dioxygenases which cleave double bonds of carotenoids at different positions resulting in the formation of apocarotenoids. zeaxanthin acts as substrate for ccd which cleaves it into crocetin dialdehyde and β-cyclocitral. ccd was very recently identified to be the enzyme which has 7′ cleavage activity and catalyses the first step in crocin biosynthesis in crocus. crocetin dialdehyde is converted into crocetin by aldehyde dehydrogenase. in this study we have identified many transcripts encoding aldehyde dehydrogenases which might include the one involved in converting crocetin dialdehyde into crocetin. further, crocetin and β-cyclocitral are converted into crocin and picrocrocin respectively by udp-glucosyl transferases. many transcripts coding for udp-glucosyl transferases have been identified in this study. so far only a few udp-glucosyl transferases were known from crocus and we have added new isoforms to the existing database. the picrocrocin is further converted into saffranal which is supposed to be a beta glucosidase action. gene encoding this enzyme has not been identified in crocus, so far. here we have identified many transcripts coding for this enzyme. therefore, the present study would provide a platform for generating knowledge about the substrate specificities and activities of the enzymes identified so as to understand the apocarotenoid biosynthetic pathway. further, there are many other ccd isoforms which act on a wide range of substrates, cleave them at different positions and produce a myriad of apocarotenoid products. analysing the enzyme activities and their substrate specificities would help to identify new compounds. this knowledge can further be useful for designing metabolic engineering strategies for either enhanced production of known compounds or production of new metabolites.fig. schematic representation of apocarotenoid pathway showing different enzymes and intermediates involved in the pathway. genes marked in red are identified first time in this study considering other branch points of this pathway, zeaxanthin is also converted into violaxanthin by zeaxanthin epoxidase which is in turn converted into neoxanthin by neoxanthin synthase. violaxanthin and neoxanthin are cleaved by 9-cis-epoxycarotenoid dioxygenase to produce xanthoxin, the direct substrate for phytohormone abscisic acid synthesis. while nced from crocus is known, the genes encoding the enzymes which act upstream are identified for the first time in this study. the schematic representation of the carotenoid/apocarotenoid pathway is given in fig.  and the list of pathway genes identified is given in table  table genes involved in carotenoid/apocarotenoid pathway agenes identified first time in the present study identification of transcription factors in order to identify transcription factor encoding genes in crocus transcriptome, homology search was performed against the plant transcription factor database. around unigenes show similarity to transcription factors representing families. out of these, transcription factors belonging to myb family were most abundant followed by c3h, far, mads box. in many other studies also, the above mentioned classes represented highly expressed transcription factor families. for example, in catharanthus roseus, myb was the most abundant family. in horse gram, most highly represented transcription factors belonged to c3h, bhlh and ap domain families. transcription factors perform key roles in plant growth and development. they are also involved in regulation of secondary metabolism and also coordinate plant’s response to various environmental cues. several members of myb family have been shown to regulate different secondary metabolic pathways. zinc finger proteins of different classes are also involved in regulation of plant secondary metabolism. mads box genes regulate a range of plant processes including flower development plant reproduction etc. so far regulatory pathway controlling crocus apocarotenoid metabolism is not known. although many transcription factors were identified from crocus, their role in apocarotenoid regulation was not experimentally validated. recently, a sand domain ultrapetala transcription factor was isolated from our laboratory and was shown to regulate apocarotenoid biosynthesis in crocus. no other information is available as far regulation of apocarotenoid biosynthesis is concerned. this transcription factor database would therefore, be an important asset to characterize the regulatory pathways of crocus secondary metabolism. this may also provide base for identification of regulators of crocus flower development, its tripartite stigma and would help in addressing other question related to biology of this plant.fig. percent crocus transcripts representing top transcription factor families and conserved domain families domain analysis we aimed to get insight about the molecular mechanism underlying the phenomenon of colour and flavour biogenesis in crocus. go classification showed the major classes in molecular function comprised of binding, catalytic activity and nucleotide binding. further, kegg mapping indicated that genes involved in protein processing, spliceosome and oxidative phosphorylation were significantly represented. in order to extrapolate this information and build up knowledge about the mechanism of action of the genes involved in such activities we performed conserved domain analysis for crocus transcriptome. the major domain classes are given in fig. 4b. the results showed that most of the genes were those with protein kinase like domain followed by protein kinase catalytic domain, serine/threonine protein kinase domain, serine/threonine-protein kinase-like domain and tyrosine-protein kinase catalytic domain. these domains are involved in phosphorylation of proteins at various amino acid residues. protein modification by phosphorylation and dephosphorylation is a crucial mechanism that controls activity of proteins and as a result regulates important cellular functions in eukaryotes such as cell division, differentiation, signal transduction etc. predominant presence of genes with different kinase domains in crocus transcriptome is indicative of the fact that even though transcriptional regulation is the central regulatory mechanism for most of the biological processes, diverse post translational events also play important role. these multiple layers of regulation help plant to synchronize properly with developmental stages and environmental cues. differential gene expression, gene ontology and pathway enrichment in crocus, the stigma part of the flower has commercial significance because its dried state forms saffron which is the site of biosynthesis and accumulation of apocarotenoids. although attempts have been made to identify and isolate genes involved in apocarotenoid biosynthesis, little has been done so far. in order to identify genes involved in biosynthesis and regulation of apocarotenoids in crocus stigma, we studied differential gene expression in crocus stigma vs. rest of the flower. a total of transcripts were differentially expressed out of which transcripts were upregulated in flower and transcripts in stigma. in order to identify the major functional categories represented by differentially expressed genes, go enrichment analysis was carried out. we observed that in flower upregulated genes, the top five classes in biological process were metabolic process, cellular process, primary metabolic process, macromolecule metabolic process and protein metabolic process. in case of the genes upregulated in stigma, metabolic process, cellular process, primary metabolic process, cellular metabolic process and macromolecule metabolic process were the top five categories. in case of the molecular function category, the top five classes represented by flower upregulated genes were catalytic activity, hydrolase activity, ion binding, cation binding and metal ion binding, while as in case of stigma top five classes were binding, catalytic activity, nucleotide binding, purine nucleotide binding and ribonucleotide binding. further, genes related to transport were also highly represented in flower while in stigma genes related to nucleotide activity were more prevalent. we also performed kegg analyses on differentially expressed genes in flower and stigma. the results showed that in case of flower, genes involved in ribosome, starch and sucrose metabolism, plant-pathogen interaction, phenylpropanoid biosynthesis and pentose and glucuronate interconversions were enriched while in case of stigma, genes related to protein processing in endoplasmic reticulum, photosynthesis, carotenoid biosynthesis, legionellosis represented the major classes. enrichment of carotenoid biosynthesis genes in stigma confirms the fact that biosynthesis of carotenoids and their subsequent degradation into apocarotenoids occurs mainly in stigma.fig. differential expression analysis of crocus transciptome. a, b bar chart showing go terms associated with genes differentially expressed in flower and stigma. c, d pie chart representing kegg classes of flower and stigma upregulated genes we also investigated the expression of genes involved in carotenoid/apocarotenoid biosynthesis and observed that most of them were more significantly enriched in stigma which is in accordance with the fact that stigma is the actual site of biosynthesis of apocarotenoids in crocus. for example, psy and pds which catalyse initial steps in carotenoid biosynthetic pathway are expressed more in stigma. further, ziso and zcd which are involved in converting phytoene into lycopene are also upregulated in stigma. it was quite interesting to see that lycopene β-lcy and bch were upregulated in stigma therefore increasing the metabolic flux towards production of zeaxanthin. earlier reports have also shown increased expression of these genes in crocus stigma. recently it was shown that ccd cleaves zeaxanthin at 7′ double bond and results in the formation of crocetin and picrocrocin which are apocarotenoids responsible for color and flavor of saffron. dge analysis showed that ccd is also upregulated in stigma thereby confirming the earlier reports. further, another isoform ccd4b is also expressed more in stigma. this enzyme cleaves carotenoids like beta-carotene at 9′ double bond and forms apocarotenoids like β-cyclocitral and β-ionone. earlier reports have also shown enhanced expression of ccd4b in stigma. thus our results are in agreement with earlier reports and also provide molecular proof for stigma being the actual site of biosynthesis of apocarotenoids which includes crocin, picrorocin, β-cyclocitral and β-ionone.fig. heat-map showing expression patterns of regulatory carotenoid pathway genes differentially upregulated in stigma and flower one of the aims of the present work was to identify transcription factors which regulate the biosynthesis and accumulation of apocarotenoids in crocus in tissue and developmental stage specific manner. our data suggested that transcription factors belonging to various families like mads box, myb, zinc finger, wrky, phd etc. were differentially expressed in both flower and stigma. we observed that different transcription factor families were enriched in stigma and flower. in this context we observed that transcription factors belonging to mads box were represented more in stigma while nac, bhlh and wrky families were enriched more in flower suggesting that they perform specialized roles in different tissues. there were many other transcription factors exhibiting expression in both the tissue types. for example, zinc finger proteins with an domain were more prevalent in stigma upregulated genes while those with zat domains were upregulated in flower. this indicates that different members of the same family might play different roles in different tissues or under different conditions. experimental validation of differential expressed genes by quantitative realtime pcr qrt-pcr of ten selected genes was performed in order to validate the differential gene expression obtained by rna-seq. the genes selected were ccd4b, ccd bch, psy, pds, gt, zinc-finger, mads box, myb and, wrky. the results indicated that expression pattern as obtained by qrt-pcr corroborated with that obtained by rna-seq for all the genes. this confirmed the reliability of rna-seq data. we also performed statistical analysis on the data obtained for these genes from rna-seq and qpcr and observed that there was very good correlation between the two.fig. qrt-pcr validation of differentially expressed genes. a expression of genes was validated using qrt-pcr and compared with their expression obtained from rna-seq. qrt-pcr was performed and values were normalized relative to the expression levels of 18s rrna in the same cdna sample. data are the means of three replicates expression data are presented as expression values of genes in stigma sample relative to their expression in flower sample correlation of gene expression results obtained from qpcr and rna-seq phylogenetic analyses of ccd gene family carotenoid cleavage dioxygenases form a group of enzymes which are involved in cleavage of carotenoids leading to the production of apocarotenoids. ccds are specific to the double bond they cleave but often exhibit substrate promiscuity which is responsible for the diversity of apocarotenoids found in nature. apocarotenoids play role in many aspects of plant growth and development. the ccd family is ancient and has its members present in bacteria, plants and animals. members of the ccd family share several characteristics like, they require a fe2+ for catalytic activity; they have four conserved histidines which coordinate iron binding and they contain a conserved peptide sequence at their carboxyl terminus. in order to get an insight about the phylogenetic relationship between various ccds, a neighbour-joining phylogenetic tree was constructed with ccd genes from plant species. the amino acid sequences from various ccds were grouped into five clusters, named ccd ccd ccd ccd and nced. in each of these clusters, ccds were present in two groups corresponding to monocotyledon and dicotyledon species. further, genes or isoforms within the ccd sub-family and belonging to same species were grouped in the same branch, for example, crocus sativus ccd4a/b, ccd8a/b. phylogenetic analysis also revealed that two major duplications had occurred in ccd subfamilies. ccd duplication, which occurred in the moss lineages, ultimately led to the emergence of two lineages that evolved into ccd and ccd this result suggests that, ccd7/ genes had similar evolutionary trends than ccd and ccd sub-family. in general, ccds belonging to a particular cluster have similar cleavage activity. recently a new ccd isoform was identified from crocus sativus, which plays an important role in crocin biosynthesis. csccd clusters with ccd subfamily but is distinct from ccd subfamily as far as its catalytic activity is concerned. this suggests that ccd might have evolved from ccd and developed different cleavage site specificity. thus phylogenetic analysis may help in understandin g functional diversity of ccd gene family.fig. phylogenetic tree for ccd proteins. tree was constructed by aligning ccd proteins from plant species including crocus using the muscle program and subjected to phylogenetic ana lysis by the distance with neighbor-joining method using the mega software. the numbers on the nodes indicate the bootstrap values after replicates  in the present study, crocus sativus transcriptome was sequenced, assembled and annotated. the database led to the identification of many new candidate genes involved in carotenoid/apocarotenoid pathway. identification of transcription factors provides a platform for unravelling the regulatory pathway of crocus flower development and apocarotenoid biosynthesis. differential gene expression and pathway mapping confirmed enrichment of apocarotenoid pathway genes in stigma thus confirming that stigma is the site of apocarotenoid biosynthesis. the transcript resource generated would therefore facilitate and enhance our understanding of biosynthetic pathway of carotenoids and their subsequent tailoring into apocarotenoids and the mechanism that regulates carotenoid-apocarotenoid metabolic flux. availability of supporting data the illumina sequence data have been submitted as bioproject to ncbi sequence read archive under accession number. all the other supporting data are included as additional files. additional files additional file 1: sequence of the primers used in this study. additional file 2: accession numbers of the genes used for phylogenetic analysis. additional file 3: functional annotation of crocus transcriptome. additional file 4: go and kegg classes of crocus transcripts. additional file 5: details of transcription factor and conserved domain families. additional file 6: details of differentially expressed transcripts and flower and stigma upregulated genes. additional file 7: go and kegg classes in flower upregulated genes. additional file 8: go and kegg classes in stigma upregulated genes. abbreviations gogene ontology keggkyoto enclyclopedia of genes and genomes ccdcarotenoid cleavage dioxygenase ggdpgeranylgeranyl diphosphate psyphytoene synthase pdsphytoene desaturase shoib ahmad baba and tabasum mohiuddin share first authorship. swaraj basu and mohit kumar swarnkar share second authorship. shoib ahmad baba tabasum mohiuddin swaraj basu and mohit kumar swarnkar contributed equally to this work. competing interests the authors declare that they have no competing interests. authors’ contributions na and aks: conceptualized and designed experiments and wrote manuscript; sab, tm, ahm, zaw and na performed experiments; sb and mks did analysis. all authors read and approved the final manuscript. 
go analysis showed that the predominant classes in biological process category include metabolic process followed by cellular process and primary metabolic process. therefore, deep transcriptome sequencing of crocus stigma and rest of the flower tissue was done to identify the genes and transcriptional regulators involved in the biosynthesis of these compounds. transcriptome of stigma and rest of the flower tissue was sequenced using illumina genome analyzer iix platform which generated, flower and, stigma reads. transcripts showed similarity to transcription factors representing families out of which myb family was most abundant. the database would also help in investigating many questions related to saffron biology including flower development. differential expression analysis indicated that key carotenoid/apocarotenoid pathway genes including phytoene synthase, phytoene desaturase and carotenoid cleavage dioxygenase are enriched in stigma thereby providing molecular proof for stigma to be the site of apocarotenoid biosynthesis. kegg mapping analysis indicated that pathways involved in ribosome, carbon and starch and sucrose metabolism were highly represented. inspite of their importance, the biosynthetic pathway of crocus apocarotenoids is not fully elucidated. many genes involved in carotenoid/apocarotenoid pathway were identified for the first time in this study which includes zeta-carotene isomerase and desaturase, carotenoid isomerase and lycopene epsilon-cyclase.
 crocus sativus stigmas form rich source of apocarotenoids like crocin, picrocrocin and saffranal which besides imparting color, flavour and aroma to saffron spice also have tremendous pharmacological properties. inspite of their importance, the biosynthetic pathway of crocus apocarotenoids is not fully elucidated. moreover, the mechanism of their stigma specific accumulation remains unknown. therefore, deep transcriptome sequencing of crocus stigma and rest of the flower tissue was done to identify the genes and transcriptional regulators involved in the biosynthesis of these compounds.  transcriptome of stigma and rest of the flower tissue was sequenced using illumina genome analyzer iix platform which generated, flower and, stigma reads. sequences were assembled de novo using trinity resulting in transcripts which were classified into unigenes comprising of clusters and singletons. a comprehensive functional annotation and gene ontology analysis was carried out.  % of the transcripts showed similarity to sequences present in public databases while rest could be specific to crocus. transcripts showed similarity to transcription factors representing families out of which myb family was most abundant. many genes involved in carotenoid/apocarotenoid pathway were identified for the first time in this study which includes zeta-carotene isomerase and desaturase, carotenoid isomerase and lycopene epsilon-cyclase. go analysis showed that the predominant classes in biological process category include metabolic process followed by cellular process and primary metabolic process. kegg mapping analysis indicated that pathways involved in ribosome, carbon and starch and sucrose metabolism were highly represented. differential expression analysis indicated that key carotenoid/apocarotenoid pathway genes including phytoene synthase, phytoene desaturase and carotenoid cleavage dioxygenase are enriched in stigma thereby providing molecular proof for stigma to be the site of apocarotenoid biosynthesis.  this data would provide a rich source for understanding the carotenoid/apocarotenoid metabolism in crocus. the database would also help in investigating many questions related to saffron biology including flower development. electronic supplementary material the online version of this article contains supplementary material, which is available to authorized users. keywords crocussaffronapocarotenoidsilluminade novo assemblyissue-copyright-statement© the author 2015 