BACKGROUND
a detailed characterization of the transcriptome is not only a prerequisite for genome annotation, but is also crucial to relate gene expression and function to phenotypic variation, to estimate evolutionary rates or to identify genes under selection  <cit> .
while the first large-scale est  projects based on sanger sequencing had already revealed that transcriptomes are highly diverse and complex, this has over the last decade been corroborated on an unprecedented scale by massive parallel sequencing   <cit> .
until recently microarrays, based on either cdna, or in the case of model organisms, oligonucleotides, were the main tool to assess global patterns of gene expression.
these microarray studies provided a first characterization of temporal and spatial gene expression patterns in many different organisms, but they also had obvious limitations since they were dependent on already available sequence information from the organism of interest.
furthermore, at least with cdna arrays, most transcripts were only represented by the most common splice variant making the gene expression estimates unspecific  <cit> .
the emergence of mps techniques has profoundly transformed the landscape of both genome and transcriptome sequencing  <cit> .
in the case of transcriptome sequencing it has become possible to obtain both sequence data and estimates of gene expression levels in a single experiment  <cit> .
mps also allowed for very deep sequencing of transcriptomes revealing in principle all expressed parts of the genome.
in mouse , for example, despite a fully sequenced genome and millions of est sequences in public databases, mps of expressed genes from three different tissues revealed several thousand new features of known genes and close to  <dig> novel gene transcripts  <cit> .
similar results, with a very large number of previously unknown expressed parts of the genome, were recently obtained in drosophila melanogaster <cit> .
in plants, mps has been used to investigate gene expression and alternative splicing in model organisms.
as in animals, large parts of the genome are transcribed  and around 30-40% of the genes have more than one expressed splice variant  <cit> .
in addition, transcriptome sequencing has been employed in a large array of plant species to create basic genetic resources and to study gene expression dynamics .
most efforts have so far focused on angiosperms, less attention being paid to gymnosperms in spite of the availability of several very large est libraries in the spruce  and pine  genus.
by april  <dig> there were  <dig>  and  <dig>  est sequences deposited at ncbi from picea glauca and p. sitchensis, respectively.
efforts in the european picea abies  have been less intensive and only some  <dig>  ests had been deposited at ncbi by april  <dig>  based on sequence similarity, the data obtained in the north american species have been collapsed into  <dig>  putative unique transcripts  in p. glauca and  <dig>  unigenes in p. sitchensis .
in addition, more focused experimental and bioinformatic efforts have been made to obtain full-length cdnas  in both p. sitchensis and p. glauca and have so far lead to the identification of  <dig>   unique transcribed genes in p. glauca and around  <dig>  validated fl-cdna have been characterized in p. sitchensis <cit> .
these data sets have been very useful for creating microarray assays, detecting polymorphic markers to construct genetic maps and to estimate evolutionary rates  <cit> .
gymnosperms separated from angiosperms some  <dig> million years ago  <cit>  and they differ in many important aspects.
for example, in conifers, the overall recombination seems to be lower than in angiosperms and limited to genes  <cit>  and, in contrast to angiosperms where genome duplication has been pervasive, genome duplication is rare in conifers  <cit> .
the generally very large size of conifer genomes is rather a result from massive invasion of transposable elements early in the history of the group  <cit> .
one could therefore expect different mode and tempo of evolution of the transcriptome in these two groups of plants.
the gymnosperm group does, compared to the angiosperm group, contain fewer species, but many gymnosperm species have very large distribution ranges and are of great ecological and economical importance  <cit> .
the largest taxon within gymnosperms is the coniferophyta that harbors around  <dig> species, including many tree species such as pines, spruces and yews  <cit> .
within these three groups the phylogenetic relationships are not fully known since molecular divergence tends to be very low.
for example, within the spruce genera morphological and molecular markers often give conflicting results and multi-locus nuclear dna data indicates that incomplete lineage sorting is a major challenge for inferring species trees  <cit> .
the spruce lineage separated from the pine lineage around 120- <dig> millions years ago whereas the taxus lineage is more distantly related and are believed to have a common ancestor to the pine and spruce lineages between  <dig> and  <dig> million years ago  <cit> .
in this study we present results from deep sequencing of the p. abies needle transcriptome.
the data was used to create a reference transcriptome, containing a extensive set of genes expressed in actively growing spruce needles.
combined with available est from publicly available databases, the size of the characterized part of the transcriptome in p. abies now reaches the same order of magnitude as data available from north american spruce species, with thousands of putative full-length genes.
we used these data to address the following questions: first, to what extent can this type of short sequence data be assembled de novo to full length transcripts?
second, can the short reads be used to identify polymorphic sites?
third, by analyzing mrna samples from two time points; during day  and during night  can we detect differentially expressed transcripts and if so does the inferred pattern differ compared to angiosperms?
fourth, since recent studies of divergence between conifer species have suggested that there might be major differences in evolutionary rates between conifers and angiosperms  <cit>  we used the assembly together with data from p. glauca, p. sitchensis, pinus taeda and taxus mairei to evaluate this statement with a larger and more diverse data set.
by comparing not only our, but also two other spruce transcriptome sequences created by traditional sequencing methods, to both pinus taeda and taxus mairei we could minimize biases that might exist in any single assembly.
furthermore, since inferring pattern of divergence between species can be biased when closely related species  or highly divergent species  are used, the use of species that diverged at two different time points should reduce these problems.
RESULTS
rna extracted from needles of a single adult p. abies tree growing outside uppsala, sweden  was used to produce  <dig> different adapter-ligated cdna libraries.
in total, close to  <dig> million  <dig> bp long fragments were sequenced from both ends and close to 50% of the reads originated from the samples collected in the dark .
based on overlap between pair-end reads, the obtained library sizes were smaller than the expected 150- <dig> bp and around 50% of the reads contained detectable traces of adaptor sequence, which means that the fragment size was smaller than the read length of  <dig> bp.
thus, reads were filtered, not only based on quality, but also on presence of adaptor sequence, reducing the data set to  <dig> , <dig> fragments that was used in de novo assembly .
for mapping purposes we created a second data set where reads containing adaptor sequence were trimmed rather than removed and retained if they were longer than  <dig> bp after trimming.
this created a data set of  <dig> , <dig> fragments  of length between  <dig> and  <dig> bp .
collected number of sequences from the  <dig> different sequence libraries.
the number of reads used for assembly and mapping divided as pair-end  and single-end  reads retained after cleaning and filtering as described in material and methods.
de novo assembly
two types of short read de novo assemblers were used.
first, both pair-end and single-end reads were assembled using velvet  <cit>  and second, all reads were treated as single-end reads and assembled using abyss  <cit> .
the resulting contigs from both assemblers were merged with  <dig>  puts  assembled from ests at plantgdb.org  <cit>  to create a more comprehensive set of puts from p. abies.
the final assembly contained  <dig>  puts with a n <dig> and a n <dig> size of  <dig> bp and  <dig> bp, respectively .
in all downstream analyses we excluded puts shorter than  <dig> bp, leaving  <dig>  puts with a mean length of  <dig>  bp.
the total assembly length of puts longer than  <dig> bp was just over  <dig> mbp and the mean sequence depth of these puts as measured by rpkm  was  <dig> .
55% of the assembled p. abies puts showed similarity  to plant proteins when translated to protein sequences.
of the  <dig>  puts  <dig>  show a significant hit against either p. glauca or p. sitchensis and  <dig>  of these were similar over more than 90% of the sequence length of the p. abies put.
of the  <dig> transcripts longer than  <dig>  bp only  <dig> puts did not show any similarity to arabidopsis or other plant proteins.
taken together this left only  <dig>  p. abies puts with no or limited similarity to any plant protein sequence or any est sequence from other spruce species.
for complete annotation of puts longer than  <dig> bp against a. thaliana, p. glauca, p. sitchensis, plant proteins, pfam databases and repeat databases see additional file  <dig>
potential reading frames were identified from all puts longer than  <dig> nucleotides.
this led to the identification of  <dig>  puts that have properties suggesting that they are full-length transcripts including both a start and stop codon and have an open reading frame  of at least  <dig> amino acids .
comparison of the length distribution of both untranslated regions  and the orf length to the validated full-length data set from both p. sitchensis <cit>  and p. glauca <cit>  suggests that a fraction of our putative full-length puts includes full-length transcripts.
however, the length distributions of both the 5’-utrs, 3’-utrs and orfs were enriched for short sequences compared to what was reported from the p. sitchensis and p. glauca .
a direct comparison between a set of potentially orthologous sequences in p. abies and p. glauca shows that there is a large number of puts that likely are full-length as they are equally long in both species, but it is also clear from this comparison that neither of these two data sets are solely composed of full-length transcripts as many transcripts annotated as full-length in either p. glauca or p. abies appear to be truncated when compared to each other .
approximately one fifth of the  <dig>  full length p. abies transcripts where orthologous sequences were found in p. glauca displayed less than 5% difference in sequence length and hence likely represents full length transcripts.
all puts containing coding frames were also annotated by comparing them to the pfam-a domains.
a comparison of the number of different detectable domains in the p. abies transcriptome to the p. glauca transcriptome is consistent with shorter length of assembled p. abies transcriptome.
we detected fewer domains of all types in the p. abies assembly, but the distributions from the two species are indeed very similar .
identification of variable sites
by mapping short reads to the p. abies put assembly we extracted potential snps using samtools and the mpileup function on the mapped reads.
the number of snps identified was, as expected, largely dependent on the criteria used in the algorithm and as the data is from non-normalized rna the commonly used settings for genome sequencing will likely not give satisfactory results.
with very stringent criteria, including a cut-off for the observed allele frequency of the minor allele at 25% , read direction for alternate alleles  and coverage  <dig>  snps were found.
less stringent criteria obviously led to a higher number of identified snps and if we only filter on quality  <dig>  potential snps were identified.
from published p. abies population genetic data set the observed transition to transversion ratio is around  <dig> , which is very similar to the ratio observed with the criteria that includes both quality and observed minor allele frequency above  <dig>  .
the effect of using different quality criteria for identifying snps.
sanger qual = sanger sequence quality, af = allele frequency, depth = sequence depth at the base dp4 = depth of the two variants according to the read direction.
for details see the samtools specifications.
comparative molecular evolution in conifers
previous genetic studies from both nuclear and organellar markers within the picea group have revealed a fairly shallow phylogeny with partly contrasting results between different marker types  <cit> .
all these studies have been restricted to either small data sets using at most  <dig> gene fragments or only organelle markers.
here we combine est data from three spruce species , a yew  and a pine  to estimate phylogenetic trees and study the rate of molecular evolution in conifers.
in all five species we identified putative coding frames and aligned orthologous sequences from all species.
this resulted in a data set of  <dig>  aligned orfs between the three picea species.
this data set was reduced to  <dig>  when adding pine and further reduced to  <dig> when adding yew.
phylogenetic inference on this latter data set confirms the very close relationship between all three spruce species and we cannot with any confidence resolve their phylogeny .
due to this lack of resolution and the large reduction in data sets as more species were added, we restricted the analysis of synonymous and nonsynonymous divergence to pairwise comparisons between one of the three spruce species, on the one hand, and either the pine or the yew, on the other hand.
this created  <dig> pairwise data sets with an average  <dig>  orfs aligned in each comparison .
mean synonymous divergences between the three spruce species and pine and yew were  <dig>   and  <dig>  , respectively.
assuming a divergence time between pine and spruce of  <dig> to  <dig> million years  <cit>  and a divergence between spruce and yew of  <dig> to  <dig> million years  <cit>  these values correspond to a yearly substitution rate of  <dig>  × 10− <dig> to  <dig>  × 10− <dig> and 1 × 10− <dig> to  <dig>  × 10− <dig>  respectively.
the ratio of non-synonymous to synonymous divergence between the spruce species and pine and yew show very few genes with ratios larger than one, irrespective if the analysis were performed on only full-length orfs or on all identified orfs .
differentially expressed genes
since sequence data were collected from three technical replicates sampled in the dark during night and three replicates sampled during daytime in normal daylight conditions, we set out to detect differentially expressed genes.
the correlations of estimated levels of gene expression based on just counting the number of pair-end hits against any given reference gene were in general very high especially for the dark samples .
<dig>  transcripts were detected as significantly differentially expressed between the two sampling points using bayseq and a cut-off posterior probability at  <dig> , but among these transcripts only few showed a large fold change .
we did not observe any significant correlation between either expression level or differential expression and selective constraint as measured by dn/ds.
discussion
in spite of the remaining challenges associated to the assembly of short reads to full-length transcripts  <cit>  rna-seq has over the last years evolved into the de-facto standard for transcriptome analysis even in non-model species  <cit> .
in the present study we evaluated the quality of our assembly using a set of different criteria largely relying on comparisons to other plant species.
first, we compared basic summary statistics to data from other conifer species, where the transcriptome was obtained through sanger sequencing.
second, we further characterized the assembly by assessing properties like transcriptome annotation and similarity at the protein level to available plant protein sequences.
p. glauca, which has one of the best characterized transcriptome among conifer species, has a transcriptome assembly length of  <dig>  mbp, i.e.
<dig> mbp longer than the the p. abies assembly presented here.
in addition, the length of individual puts is on average shorter in the p. abies assembly than in both p. glauca and p. taeda .
this indicates that the p. abies transcriptome assembly only comprises a fraction of the total transcriptome and that a large proportion of the assembled puts corresponds to partial transcripts.
there are a number of possible reasons for this difference in assembly length despite of the large amount of short reads collected.
first, the assembly of short reads from transcriptome sequence is still far from optimal  <cit> , something that likely is elevated in this data set as the mean fragment length in the sequence library was shorter than the pair-end distance.
second, non-normalized mrna was used, which means that the read depth of lowly expressed genes is likely insufficient for assembly.
third, the short read libraries used here obviously will capture only a restricted part of the total transcriptome as they are based on a single tissue and only two sampling time points.
all this is in contrast to the other mentioned conifer species est assemblies that are, at least partially, composed of normalized mrna libraries and also include mrna extracted from several tissues including root, which in general have divergent expression patterns  <cit> .
in total over  <dig>  transcripts showed significant hits when compared to other spruce species, either at the protein level to plant proteins or at the nucleotide level, leaving only  <dig>  puts without any annotation.
this set of puts with low or no similarity to other species could represent erroneously assembled puts, but since estimates from the p. glauca transcriptome suggest that as much as  <dig>  mbp of the transcriptome is still missing  <cit> , they could also represent previously undiscovered spruce transcripts.
in theory they could also be transcripts specific to p. abies, but this seems less likely as both this data and previous studies have shown that spruce species are very closely related and even share many polymorphic sites  <cit> .
still, genome size between different spruce species is variable , which could imply significant differences in gene content between them.
finally, and perhaps most likely, they may represent assemblies containing mainly non-coding rna and partial transcripts with mostly utr sequence that, in general, show lower degree of conservation between species.
in summary, even if the p. abies transcriptome contains a large fraction of partial transcripts and is far from complete, the assembly presented here is of a high enough quality and size to be informative on the general properties of the transcriptome.
in fact, orf predictions identified  <dig>  puts from p. abies with properties that suggest that they are full-length transcripts.
a direct comparison of these puts towards full-length transcripts from p. glauca indicate that neither ours nor the p. glauca put assembly contains a comprehensive set of full-length transcripts and that despite millions of ests from spruce species there is still room for efforts aiming at capturing and sequencing full-length transcripts, such as oligo-capping, something that is of special importance in relation to upcoming genome annotations  <cit> .
in many of the puts with high enough sequence coverage we could extract potentially variable sites, but since sequencing technologies are evolving at a rapid pace, the understanding of the error profiles for a given technique is still sparse making this process nontrivial and prone to false positives.
this process is, in the present case, further complicated by the fact that the majority of tools developed to extract snps assumes that genome data is available and that sequence depth is even for the underlying genotypes, meaning that one expects around half the reads from each allele at heterozygous sites.
in addition, many tools make use of uneven coverage as a tool to detect copy number variants or duplicated regions making standard methods to discover snps difficult to use as the uneven coverage and allele specific gene expression will cause deviations from these expectations.
by merging population genetics sequence data sets from both unpublished and published projects in p. abies, where the data have been collected with sanger sequencing  <cit>  we compared the ratio of transitions to transversions to the ratio obtained here under different cut-offs.
adding both quality and coverage of the alternate alleles as criteria yielded almost identical transition to transversion ratio as the sanger data, suggesting that this criteria might be suitable for this data set .
hence, from a single individual study we have identified close to  <dig> thousand snps, but because all of these stem from only one individual a fraction of them will be either singletons or low frequency variants within the species.
even so these results suggests that mrna sequencing can be of great use for the identification of population genetic markers.
however, recent discussions have identified several issues in identification of variable sites from cdna something that calls for caution and further validation before using snps from cdna as a source for identification of genetic markers  <cit> .
the observed mean synonymous divergence of  <dig>  between spruce and pine obtained here, is similar to the  <dig>  reported as a mean over a phylogenetic tree of four conifer species  <cit>  and the  <dig>  reported from a pairwise comparison of p. sitchensis and p. taeda <cit> .
in the comparisons to yew we obtained a mean silent divergence around  <dig> .
these divergence estimates translate into average synonymous substitution rates of  <dig>  × 10− <dig> and  <dig>  × 10− <dig>  values that are also in line with previous estimates of  <dig>  −  <dig>  × 10− <dig> <cit> .
as previously noted this yearly synonymous substitution rate is much lower than estimates of substitution rate from annual angiosperm species.
data from numerous studies have obtained similar patterns and there is now substantial evidence supporting a lower annual substitution rate in gymnosperms compared to many annual angiosperms  <cit> .
since lower synonymous substitution rates per year have also been observed in angiosperms trees/shrubs compared to herbs  <cit>  one of the main causes seem simply to be that gymnosperms are trees with a long generation times and if substitution rates are calculated per generation rather than per year, the estimates between perennial gymnosperms presented here are similar to estimates from annual angiosperms.
the mechanisms underlying the slower substitution rates in perennials than in annuals remain unclear  <cit>  and we will not speculate on its causes here.
interestingly though, available data also suggest that polymorphism level is lower in trees, and in particular in conifers, than in herbs .
this apparent effect of generation time on nucleotide diversity could reflect the fact that organisms with long generation time and low substitution rate will require longer period of time to recover from past decline in effective population sizes.
also, for a given number of generations they will span much longer periods of time than annuals and thereby will be more likely to experience environmental changes causing variation in population size.
the ratio of non-synonymous to synonymous divergence at full length orfs between spruce species and either pine or yew reported here,  <dig>   and  <dig>  , respectively, are higher than estimates from previous analyses in conifer species  <cit> , where branch specific estimates were used , but lower than the pairwise estimate between p. taeda and p. sitchensis )  <cit> .
the deviation from the former study is likely due to the fact that they included only sequences present in all four conifer species used in their study thereby limiting severely the number of sequences  and likely enriching the data set for conserved sequences.
the difference with the latter study is harder to explain, but stems mainly from differences in the way coding frames and orthology were defined.
in buschiazzo et al.
<cit>   <dig> comparisons revealed a dn/ds ratio higher than one, and most of these estimates are based on short alignments where both the number of synonymous sites and changes were small.
this is in stark contrast to our results with very few values larger than one and a weaker correlation between dn/ds ratios and alignment length .
buschiazzo et al.
<cit>  also estimated dn/ds between arabidopsis thaliana and populus trichocarpa and obtained a value of  <dig> .
they therefore concluded that dn/ds was higher in conifers than in angiosperms.
while our results for the conifer data are within the same order of magnitude, it is worth noting that the estimate they obtained between a. thaliana and p. trichocarpa tends to be much lower than previous estimates from angiosperms.
for example, in estimates comparing a. thaliana and a. lyrata the mean over more than  <dig> genes was  <dig>  and within the poplar lineage p. trichocarpa and p. tremula the mean over almost  <dig> genes was  <dig>   <cit> .
it is also worth noting that the observed synonymous divergence between pine and spruce are more similar to the values obtained within the arabidopsis and populus lineages and since synonymous divergence seem to have a strong impact in estimates of dn/ds  <cit>  comparing ratios over different groups of species should only be done when the observed synonymous divergence is similar.
the latter, together with other uncertainty factors discussed above, suggests that caution is warranted before concluding that there are biologically meaningful differences in selective constraints between angiosperms and gymnosperms.
even when the correct orthologous sequences have been aligned, individual estimates of dn/ds values should not be directly used to infer positive selection as recent work have highlighted several problems in using dn/ds larger than one as indicative of positive selection  <cit> .
for example, the vast majority of human genes with dn/ds ratio larger than one is due to low ds values rather than high dn values.
hence, the classic interpretation of high dn/ds ratios is questionable as relatively low ds values will have the same effect on the ratio as high dn values.
furthermore, species pairs that are too closely related will tend to overestimate dn/ds.
in model organisms large scale studies of gene expression have identified numerous properties of genes and proteins that significantly correlate with selective constraints.
here we investigated whether level of gene expression in needles correlated with selective constraints, but found no significant pattern.
in animals, fungi and plants a correlation between expression and selective constraint, such that highly expressed genes tend to have lower dn/ds, has been reported  <cit> .
in multicellular organisms the correlation between dn/ds and expression breadth  is often stronger than with the actual level of gene expression  <cit> .
the lack of correlation in spruce could hence be due to the fact that we only measured gene expression in needles and have no data on expression breadth.
also, since our transcriptome reference sequence is largely created from short reads from needle rna our reference is biased and contains a specific set of genes that do not fully represent the complete transcriptome; e.g our approach has limited information from genes lowly expressed in needles.
the lack of correlation reported here should therefore be taken with a grain of salt as it might simply stem from lack of data rather than from an actual lack of correlation.
analysis of level of gene expression in samples collected in the dark compared to samples collected in the light revealed only a small set of genes showing at least a 2-fold difference in gene expression.
angiosperms studies suggest that around 20% of the transcriptome is differentially expressed between light and dark treatments, even though the number of genes varies depending on species, tissue and actual treatment  <cit> .
the pattern observed in our data is different and might suggest that the diurnally expressed genes in gymnosperm trees could be fewer than in angiosperms.
this is consistent with earlier studies, reporting a lack of clear diurnal expression pattern at key photosynthesis genes in gymnosperm species  <cit> .
to be noted here is that all expression estimates put forward here stems from a single individual and is hence not suitable for making strong statements at the species level.
CONCLUSIONS
our study highlights the versatility of next generation sequencing technology in generation of full-length expressed genes, identification of polymorphic sites and estimation of gene expression levels.
the available data from norway spruce have with this study been leveraged with available data from north american spruce species.
this allowed the comparison of sequence evolution in this group of plants for thousands of full-length genes.
by comparing our p. abies assembly to other gymnosperm species our analysis suggests that evolutionary constraints might not, as some previous work have suggested, be very dissimilar between gymnosperm and angiosperm species.
more importantly, our analysis does, together with the analysis of buschiazzo et al.
<cit> , highlights some of the problems associated with inference of selection from dn/ds in systems without annotated genome data.
even with a reference genome it is often difficult to unambiguously identify correct orfs and orthologous sequences and one should consequently refrain from strong statements about patterns of divergence and selection in systems without high quality annotation data.
however, as the collection of transcriptome sequence data is getting more streamlined and standardized  it will soon be possible to sample both more individuals and species and thereby facilitates more accurate characterization of selective constraints.
the addition of more species will also make it feasible to pinpoint on which branches of evolutionary trees selection has been acting  <cit> .
furthermore, the recent in-depth analysis of factors contributing to protein evolution in model plants  <cit>  will be within reach even in species without a reference genome and will hopefully reveal if there indeed is a difference in selective constraints between angiosperms and gymnosperms.

 a detailed characterization of the transcriptome is not only a prerequisite for genome annotation, but is also crucial to relate gene expression and function to phenotypic variation, to estimate evolutionary rates or to identify genes under selection. while the first large-scale est projects based on sanger sequencing had already revealed that transcriptomes are highly diverse and complex, this has over the last decade been corroborated on an unprecedented scale by massive parallel sequencing. until recently microarrays, based on either cdna, or in the case of model organisms, oligonucleotides, were the main tool to assess global patterns of gene expression. these microarray studies provided a first characterization of temporal and spatial gene expression patterns in many different organisms, but they also had obvious limitations since they were dependent on already available sequence information from the organism of interest. furthermore, at least with cdna arrays, most transcripts were only represented by the most common splice variant making the gene expression estimates unspecific. the emergence of mps techniques has profoundly transformed the landscape of both genome and transcriptome sequencing. in the case of transcriptome sequencing it has become possible to obtain both sequence data and estimates of gene expression levels in a single experiment. mps also allowed for very deep sequencing of transcriptomes revealing in principle all expressed parts of the genome. in mouse, for example, despite a fully sequenced genome and millions of est sequences in public databases, mps of expressed genes from three different tissues revealed several thousand new features of known genes and close to novel gene transcripts. similar results, with a very large number of previously unknown expressed parts of the genome, were recently obtained in drosophila melanogaster. in plants, mps has been used to investigate gene expression and alternative splicing in model organisms. as in animals, large parts of the genome are transcribed and around 30-40% of the genes have more than one expressed splice variant. in addition, transcriptome sequencing has been employed in a large array of plant species to create basic genetic resources and to study gene expression dynamics. most efforts have so far focused on angiosperms, less attention being paid to gymnosperms in spite of the availability of several very large est libraries in the spruce and pine genus. by april there were and est sequences deposited at ncbi from picea glauca and p. sitchensis, respectively. efforts in the european picea abies have been less intensive and only some ests had been deposited at ncbi by april based on sequence similarity, the data obtained in the north american species have been collapsed into putative unique transcripts in p. glauca and unigenes in p. sitchensis. in addition, more focused experimental and bioinformatic efforts have been made to obtain full-length cdnas in both p. sitchensis and p. glauca and have so far lead to the identification of unique transcribed genes in p. glauca and around validated fl-cdna have been characterized in p. sitchensis. these data sets have been very useful for creating microarray assays, detecting polymorphic markers to construct genetic maps and to estimate evolutionary rates. gymnosperms separated from angiosperms some million years ago and they differ in many important aspects. for example, in conifers, the overall recombination seems to be lower than in angiosperms and limited to genes and, in contrast to angiosperms where genome duplication has been pervasive, genome duplication is rare in conifers. the generally very large size of conifer genomes is rather a result from massive invasion of transposable elements early in the history of the group. one could therefore expect different mode and tempo of evolution of the transcriptome in these two groups of plants. the gymnosperm group does, compared to the angiosperm group, contain fewer species, but many gymnosperm species have very large distribution ranges and are of great ecological and economical importance. the largest taxon within gymnosperms is the coniferophyta that harbors around species, including many tree species such as pines, spruces and yews. within these three groups the phylogenetic relationships are not fully known since molecular divergence tends to be very low. for example, within the spruce genera morphological and molecular markers often give conflicting results and multi-locus nuclear dna data indicates that incomplete lineage sorting is a major challenge for inferring species trees. the spruce lineage separated from the pine lineage around 120- millions years ago whereas the taxus lineage is more distantly related and are believed to have a common ancestor to the pine and spruce lineages between and million years ago. in this study we present results from deep sequencing of the p. abies needle transcriptome. the data was used to create a reference transcriptome, containing a extensive set of genes expressed in actively growing spruce needles. combined with available est from publicly available databases, the size of the characterized part of the transcriptome in p. abies now reaches the same order of magnitude as data available from north american spruce species, with thousands of putative full-length genes. we used these data to address the following questions: first, to what extent can this type of short sequence data be assembled de novo to full length transcripts? second, can the short reads be used to identify polymorphic sites? third, by analyzing mrna samples from two time points; during day and during night can we detect differentially expressed transcripts and if so does the inferred pattern differ compared to angiosperms? fourth, since recent studies of divergence between conifer species have suggested that there might be major differences in evolutionary rates between conifers and angiosperms we used the assembly together with data from p. glauca, p. sitchensis, pinus taeda and taxus mairei to evaluate this statement with a larger and more diverse data set. by comparing not only our, but also two other spruce transcriptome sequences created by traditional sequencing methods, to both pinus taeda and taxus mairei we could minimize biases that might exist in any single assembly. furthermore, since inferring pattern of divergence between species can be biased when closely related species or highly divergent species are used, the use of species that diverged at two different time points should reduce these problems.  rna extracted from needles of a single adult p. abies tree growing outside uppsala, sweden was used to produce different adapter-ligated cdna libraries. in total, close to million bp long fragments were sequenced from both ends and close to 50% of the reads originated from the samples collected in the dark. based on overlap between pair-end reads, the obtained library sizes were smaller than the expected 150- bp and around 50% of the reads contained detectable traces of adaptor sequence, which means that the fragment size was smaller than the read length of bp. thus, reads were filtered, not only based on quality, but also on presence of adaptor sequence, reducing the data set to, fragments that was used in de novo assembly. for mapping purposes we created a second data set where reads containing adaptor sequence were trimmed rather than removed and retained if they were longer than bp after trimming. this created a data set of, fragments of length between and bp. collected number of sequences from the different sequence libraries. the number of reads used for assembly and mapping divided as pair-end and single-end reads retained after cleaning and filtering as described in material and methods. de novo assembly two types of short read de novo assemblers were used. first, both pair-end and single-end reads were assembled using velvet and second, all reads were treated as single-end reads and assembled using abyss. the resulting contigs from both assemblers were merged with puts assembled from ests at plantgdb.org to create a more comprehensive set of puts from p. abies. the final assembly contained puts with a n and a n size of bp and bp, respectively. in all downstream analyses we excluded puts shorter than bp, leaving puts with a mean length of bp. the total assembly length of puts longer than bp was just over mbp and the mean sequence depth of these puts as measured by rpkm was. 55% of the assembled p. abies puts showed similarity to plant proteins when translated to protein sequences. of the puts show a significant hit against either p. glauca or p. sitchensis and of these were similar over more than 90% of the sequence length of the p. abies put. of the transcripts longer than bp only puts did not show any similarity to arabidopsis or other plant proteins. taken together this left only p. abies puts with no or limited similarity to any plant protein sequence or any est sequence from other spruce species. for complete annotation of puts longer than bp against a. thaliana, p. glauca, p. sitchensis, plant proteins, pfam databases and repeat databases see additional file potential reading frames were identified from all puts longer than nucleotides. this led to the identification of puts that have properties suggesting that they are full-length transcripts including both a start and stop codon and have an open reading frame of at least amino acids. comparison of the length distribution of both untranslated regions and the orf length to the validated full-length data set from both p. sitchensis and p. glauca suggests that a fraction of our putative full-length puts includes full-length transcripts. however, the length distributions of both the 5’-utrs, 3’-utrs and orfs were enriched for short sequences compared to what was reported from the p. sitchensis and p. glauca. a direct comparison between a set of potentially orthologous sequences in p. abies and p. glauca shows that there is a large number of puts that likely are full-length as they are equally long in both species, but it is also clear from this comparison that neither of these two data sets are solely composed of full-length transcripts as many transcripts annotated as full-length in either p. glauca or p. abies appear to be truncated when compared to each other. approximately one fifth of the full length p. abies transcripts where orthologous sequences were found in p. glauca displayed less than 5% difference in sequence length and hence likely represents full length transcripts. all puts containing coding frames were also annotated by comparing them to the pfam-a domains. a comparison of the number of different detectable domains in the p. abies transcriptome to the p. glauca transcriptome is consistent with shorter length of assembled p. abies transcriptome. we detected fewer domains of all types in the p. abies assembly, but the distributions from the two species are indeed very similar. identification of variable sites by mapping short reads to the p. abies put assembly we extracted potential snps using samtools and the mpileup function on the mapped reads. the number of snps identified was, as expected, largely dependent on the criteria used in the algorithm and as the data is from non-normalized rna the commonly used settings for genome sequencing will likely not give satisfactory results. with very stringent criteria, including a cut-off for the observed allele frequency of the minor allele at 25%, read direction for alternate alleles and coverage snps were found. less stringent criteria obviously led to a higher number of identified snps and if we only filter on quality potential snps were identified. from published p. abies population genetic data set the observed transition to transversion ratio is around, which is very similar to the ratio observed with the criteria that includes both quality and observed minor allele frequency above. the effect of using different quality criteria for identifying snps. sanger qual = sanger sequence quality, af = allele frequency, depth = sequence depth at the base dp4 = depth of the two variants according to the read direction. for details see the samtools specifications. comparative molecular evolution in conifers previous genetic studies from both nuclear and organellar markers within the picea group have revealed a fairly shallow phylogeny with partly contrasting results between different marker types. all these studies have been restricted to either small data sets using at most gene fragments or only organelle markers. here we combine est data from three spruce species, a yew and a pine to estimate phylogenetic trees and study the rate of molecular evolution in conifers. in all five species we identified putative coding frames and aligned orthologous sequences from all species. this resulted in a data set of aligned orfs between the three picea species. this data set was reduced to when adding pine and further reduced to when adding yew. phylogenetic inference on this latter data set confirms the very close relationship between all three spruce species and we cannot with any confidence resolve their phylogeny. due to this lack of resolution and the large reduction in data sets as more species were added, we restricted the analysis of synonymous and nonsynonymous divergence to pairwise comparisons between one of the three spruce species, on the one hand, and either the pine or the yew, on the other hand. this created pairwise data sets with an average orfs aligned in each comparison. mean synonymous divergences between the three spruce species and pine and yew were and, respectively. assuming a divergence time between pine and spruce of to million years and a divergence between spruce and yew of to million years these values correspond to a yearly substitution rate of  × 10− to  × 10− and 1 × 10− to  × 10− respectively. the ratio of non-synonymous to synonymous divergence between the spruce species and pine and yew show very few genes with ratios larger than one, irrespective if the analysis were performed on only full-length orfs or on all identified orfs. differentially expressed genes since sequence data were collected from three technical replicates sampled in the dark during night and three replicates sampled during daytime in normal daylight conditions, we set out to detect differentially expressed genes. the correlations of estimated levels of gene expression based on just counting the number of pair-end hits against any given reference gene were in general very high especially for the dark samples. transcripts were detected as significantly differentially expressed between the two sampling points using bayseq and a cut-off posterior probability at, but among these transcripts only few showed a large fold change. we did not observe any significant correlation between either expression level or differential expression and selective constraint as measured by dn/ds. discussion in spite of the remaining challenges associated to the assembly of short reads to full-length transcripts rna-seq has over the last years evolved into the de-facto standard for transcriptome analysis even in non-model species. in the present study we evaluated the quality of our assembly using a set of different criteria largely relying on comparisons to other plant species. first, we compared basic summary statistics to data from other conifer species, where the transcriptome was obtained through sanger sequencing. second, we further characterized the assembly by assessing properties like transcriptome annotation and similarity at the protein level to available plant protein sequences. p. glauca, which has one of the best characterized transcriptome among conifer species, has a transcriptome assembly length of mbp, i.e. mbp longer than the the p. abies assembly presented here. in addition, the length of individual puts is on average shorter in the p. abies assembly than in both p. glauca and p. taeda. this indicates that the p. abies transcriptome assembly only comprises a fraction of the total transcriptome and that a large proportion of the assembled puts corresponds to partial transcripts. there are a number of possible reasons for this difference in assembly length despite of the large amount of short reads collected. first, the assembly of short reads from transcriptome sequence is still far from optimal, something that likely is elevated in this data set as the mean fragment length in the sequence library was shorter than the pair-end distance. second, non-normalized mrna was used, which means that the read depth of lowly expressed genes is likely insufficient for assembly. third, the short read libraries used here obviously will capture only a restricted part of the total transcriptome as they are based on a single tissue and only two sampling time points. all this is in contrast to the other mentioned conifer species est assemblies that are, at least partially, composed of normalized mrna libraries and also include mrna extracted from several tissues including root, which in general have divergent expression patterns. in total over transcripts showed significant hits when compared to other spruce species, either at the protein level to plant proteins or at the nucleotide level, leaving only puts without any annotation. this set of puts with low or no similarity to other species could represent erroneously assembled puts, but since estimates from the p. glauca transcriptome suggest that as much as mbp of the transcriptome is still missing, they could also represent previously undiscovered spruce transcripts. in theory they could also be transcripts specific to p. abies, but this seems less likely as both this data and previous studies have shown that spruce species are very closely related and even share many polymorphic sites. still, genome size between different spruce species is variable, which could imply significant differences in gene content between them. finally, and perhaps most likely, they may represent assemblies containing mainly non-coding rna and partial transcripts with mostly utr sequence that, in general, show lower degree of conservation between species. in summary, even if the p. abies transcriptome contains a large fraction of partial transcripts and is far from complete, the assembly presented here is of a high enough quality and size to be informative on the general properties of the transcriptome. in fact, orf predictions identified puts from p. abies with properties that suggest that they are full-length transcripts. a direct comparison of these puts towards full-length transcripts from p. glauca indicate that neither ours nor the p. glauca put assembly contains a comprehensive set of full-length transcripts and that despite millions of ests from spruce species there is still room for efforts aiming at capturing and sequencing full-length transcripts, such as oligo-capping, something that is of special importance in relation to upcoming genome annotations. in many of the puts with high enough sequence coverage we could extract potentially variable sites, but since sequencing technologies are evolving at a rapid pace, the understanding of the error profiles for a given technique is still sparse making this process nontrivial and prone to false positives. this process is, in the present case, further complicated by the fact that the majority of tools developed to extract snps assumes that genome data is available and that sequence depth is even for the underlying genotypes, meaning that one expects around half the reads from each allele at heterozygous sites. in addition, many tools make use of uneven coverage as a tool to detect copy number variants or duplicated regions making standard methods to discover snps difficult to use as the uneven coverage and allele specific gene expression will cause deviations from these expectations. by merging population genetics sequence data sets from both unpublished and published projects in p. abies, where the data have been collected with sanger sequencing we compared the ratio of transitions to transversions to the ratio obtained here under different cut-offs. adding both quality and coverage of the alternate alleles as criteria yielded almost identical transition to transversion ratio as the sanger data, suggesting that this criteria might be suitable for this data set. hence, from a single individual study we have identified close to thousand snps, but because all of these stem from only one individual a fraction of them will be either singletons or low frequency variants within the species. even so these results suggests that mrna sequencing can be of great use for the identification of population genetic markers. however, recent discussions have identified several issues in identification of variable sites from cdna something that calls for caution and further validation before using snps from cdna as a source for identification of genetic markers. the observed mean synonymous divergence of between spruce and pine obtained here, is similar to the reported as a mean over a phylogenetic tree of four conifer species and the reported from a pairwise comparison of p. sitchensis and p. taeda. in the comparisons to yew we obtained a mean silent divergence around. these divergence estimates translate into average synonymous substitution rates of  × 10− and  × 10− values that are also in line with previous estimates of  −   × 10−. as previously noted this yearly synonymous substitution rate is much lower than estimates of substitution rate from annual angiosperm species. data from numerous studies have obtained similar patterns and there is now substantial evidence supporting a lower annual substitution rate in gymnosperms compared to many annual angiosperms. since lower synonymous substitution rates per year have also been observed in angiosperms trees/shrubs compared to herbs one of the main causes seem simply to be that gymnosperms are trees with a long generation times and if substitution rates are calculated per generation rather than per year, the estimates between perennial gymnosperms presented here are similar to estimates from annual angiosperms. the mechanisms underlying the slower substitution rates in perennials than in annuals remain unclear and we will not speculate on its causes here. interestingly though, available data also suggest that polymorphism level is lower in trees, and in particular in conifers, than in herbs. this apparent effect of generation time on nucleotide diversity could reflect the fact that organisms with long generation time and low substitution rate will require longer period of time to recover from past decline in effective population sizes. also, for a given number of generations they will span much longer periods of time than annuals and thereby will be more likely to experience environmental changes causing variation in population size. the ratio of non-synonymous to synonymous divergence at full length orfs between spruce species and either pine or yew reported here, and, respectively, are higher than estimates from previous analyses in conifer species, where branch specific estimates were used, but lower than the pairwise estimate between p. taeda and p. sitchensis ). the deviation from the former study is likely due to the fact that they included only sequences present in all four conifer species used in their study thereby limiting severely the number of sequences and likely enriching the data set for conserved sequences. the difference with the latter study is harder to explain, but stems mainly from differences in the way coding frames and orthology were defined. in buschiazzo et al. comparisons revealed a dn/ds ratio higher than one, and most of these estimates are based on short alignments where both the number of synonymous sites and changes were small. this is in stark contrast to our results with very few values larger than one and a weaker correlation between dn/ds ratios and alignment length. buschiazzo et al. also estimated dn/ds between arabidopsis thaliana and populus trichocarpa and obtained a value of. they therefore concluded that dn/ds was higher in conifers than in angiosperms. while our results for the conifer data are within the same order of magnitude, it is worth noting that the estimate they obtained between a. thaliana and p. trichocarpa tends to be much lower than previous estimates from angiosperms. for example, in estimates comparing a. thaliana and a. lyrata the mean over more than genes was and within the poplar lineage p. trichocarpa and p. tremula the mean over almost genes was. it is also worth noting that the observed synonymous divergence between pine and spruce are more similar to the values obtained within the arabidopsis and populus lineages and since synonymous divergence seem to have a strong impact in estimates of dn/ds comparing ratios over different groups of species should only be done when the observed synonymous divergence is similar. the latter, together with other uncertainty factors discussed above, suggests that caution is warranted before concluding that there are biologically meaningful differences in selective constraints between angiosperms and gymnosperms. even when the correct orthologous sequences have been aligned, individual estimates of dn/ds values should not be directly used to infer positive selection as recent work have highlighted several problems in using dn/ds larger than one as indicative of positive selection. for example, the vast majority of human genes with dn/ds ratio larger than one is due to low ds values rather than high dn values. hence, the classic interpretation of high dn/ds ratios is questionable as relatively low ds values will have the same effect on the ratio as high dn values. furthermore, species pairs that are too closely related will tend to overestimate dn/ds. in model organisms large scale studies of gene expression have identified numerous properties of genes and proteins that significantly correlate with selective constraints. here we investigated whether level of gene expression in needles correlated with selective constraints, but found no significant pattern. in animals, fungi and plants a correlation between expression and selective constraint, such that highly expressed genes tend to have lower dn/ds, has been reported. in multicellular organisms the correlation between dn/ds and expression breadth is often stronger than with the actual level of gene expression. the lack of correlation in spruce could hence be due to the fact that we only measured gene expression in needles and have no data on expression breadth. also, since our transcriptome reference sequence is largely created from short reads from needle rna our reference is biased and contains a specific set of genes that do not fully represent the complete transcriptome; e.g our approach has limited information from genes lowly expressed in needles. the lack of correlation reported here should therefore be taken with a grain of salt as it might simply stem from lack of data rather than from an actual lack of correlation. analysis of level of gene expression in samples collected in the dark compared to samples collected in the light revealed only a small set of genes showing at least a 2-fold difference in gene expression. angiosperms studies suggest that around 20% of the transcriptome is differentially expressed between light and dark treatments, even though the number of genes varies depending on species, tissue and actual treatment. the pattern observed in our data is different and might suggest that the diurnally expressed genes in gymnosperm trees could be fewer than in angiosperms. this is consistent with earlier studies, reporting a lack of clear diurnal expression pattern at key photosynthesis genes in gymnosperm species. to be noted here is that all expression estimates put forward here stems from a single individual and is hence not suitable for making strong statements at the species level.  our study highlights the versatility of next generation sequencing technology in generation of full-length expressed genes, identification of polymorphic sites and estimation of gene expression levels. the available data from norway spruce have with this study been leveraged with available data from north american spruce species. this allowed the comparison of sequence evolution in this group of plants for thousands of full-length genes. by comparing our p. abies assembly to other gymnosperm species our analysis suggests that evolutionary constraints might not, as some previous work have suggested, be very dissimilar between gymnosperm and angiosperm species. more importantly, our analysis does, together with the analysis of buschiazzo et al., highlights some of the problems associated with inference of selection from dn/ds in systems without annotated genome data. even with a reference genome it is often difficult to unambiguously identify correct orfs and orthologous sequences and one should consequently refrain from strong statements about patterns of divergence and selection in systems without high quality annotation data. however, as the collection of transcriptome sequence data is getting more streamlined and standardized it will soon be possible to sample both more individuals and species and thereby facilitates more accurate characterization of selective constraints. the addition of more species will also make it feasible to pinpoint on which branches of evolutionary trees selection has been acting. furthermore, the recent in-depth analysis of factors contributing to protein evolution in model plants will be within reach even in species without a reference genome and will hopefully reveal if there indeed is a difference in selective constraints between angiosperms and gymnosperms. 
of the puts longer than bp in this reference assembly, % show similarity to ests from other spruce species and of the remaining puts, show similarity to protein sequences from other plant species, leaving puts with limited similarity to currently available plant proteins. by predicting coding frames and comparing not only the norway spruce puts, but also puts from the close relatives picea glauca and picea sitchensis to both pinus taeda and taxus mairei, we obtained estimates of synonymous and non-synonymous divergence among conifer species. a detailed knowledge about spatial and temporal gene expression is important for understanding both the function of genes and their evolution. with the development of next generation sequencing, a single experiment can now simultaneously identify the transcribed part of a species genome and estimate levels of gene expression. in addition, we detected close to snps of high quality and estimated gene expression differences between samples collected under dark and light conditions. our study yielded a large number of single nucleotide polymorphisms as well as estimates of gene expression on transcriptome scale. a de novo assembly of these reads, together with publicly available expressed sequence tag data from norway spruce, was used to create a reference transcriptome. the estimates of the dn/ds ratio reported here are in general much lower than and only a few genes showed a ratio larger than in agreement with a recent study we find that the synonymous substitution rate per year is an order of magnitude smaller than values reported for angiosperm herbs. mrna from actively growing needles of norway spruce was sequenced using next generation sequencing technology. for the vast majority of species, transcriptomes are still largely uncharacterized and even in those where substantial information is available it is often in the form of partially sequenced transcriptomes. however, if one takes generation time into account, most of this difference disappears.
 a detailed knowledge about spatial and temporal gene expression is important for understanding both the function of genes and their evolution. for the vast majority of species, transcriptomes are still largely uncharacterized and even in those where substantial information is available it is often in the form of partially sequenced transcriptomes. with the development of next generation sequencing, a single experiment can now simultaneously identify the transcribed part of a species genome and estimate levels of gene expression.  mrna from actively growing needles of norway spruce was sequenced using next generation sequencing technology. in total, close to million fragments with a length of bp were sequenced resulting in gbp of raw data. a de novo assembly of these reads, together with publicly available expressed sequence tag data from norway spruce, was used to create a reference transcriptome. of the puts longer than bp in this reference assembly, % show similarity to ests from other spruce species and of the remaining puts, show similarity to protein sequences from other plant species, leaving puts with limited similarity to currently available plant proteins. by predicting coding frames and comparing not only the norway spruce puts, but also puts from the close relatives picea glauca and picea sitchensis to both pinus taeda and taxus mairei, we obtained estimates of synonymous and non-synonymous divergence among conifer species. in addition, we detected close to snps of high quality and estimated gene expression differences between samples collected under dark and light conditions.  our study yielded a large number of single nucleotide polymorphisms as well as estimates of gene expression on transcriptome scale. in agreement with a recent study we find that the synonymous substitution rate per year is an order of magnitude smaller than values reported for angiosperm herbs. however, if one takes generation time into account, most of this difference disappears. the estimates of the dn/ds ratio reported here are in general much lower than and only a few genes showed a ratio larger than 