BACKGROUND
douglas-fir  franco) is an ecologically highly variable species that occurs in two main varieties in north america.
the natural range of the coastal or green douglas-fir  extends over  <dig>  km from the pacific northwestern fog belt and the adjacent summer-dry coastal range and cascade mountains to the drier coastland of southern california.
the interior or blue douglas-fir  is distributed over more than  <dig>  km along the dry continental climates of the montane to the subalpine rocky mountains from alberta to colorado with isolated populations reaching into mexico.
douglas-fir grows from sea level on vancouver island up to  <dig>  m altitude in the southern rocky mountains
<cit> .
within its natural range, douglas-fir has evolved into a variety of genetically diverse populations adapted to contrasting ecozones .
douglas-fir populations differ in their response to frost
<cit> , drought
<cit> , and along environmental gradients
<cit> .
like most conifer species, it is able to cope with limitations in soil water availability within its natural range
<cit> .
there is a negative relationship between shoot water potential and the photosynthesis rate
<cit> , which decreased by about 70% in water-stressed trees with a pre-dawn shoot water potential of about - <dig>  mpa.
in conifers such as douglas-fir or pine, the recovery of photosynthesis upon rainfall and rewatering occurs within days together with the rapid recovery of predawn shoot water potential from stressed , or mildly stressed  to values higher than - <dig>  mpa
<cit> .
this high ecological, genetical and physiological diversity provides an excellent system to study the adaptation of conifer trees to contrasting environments.
due to its rapid growth and favorable wood quality, douglas-fir is an economically relevant species for timber production.
in europe, the area of douglas-fir production is rising
<cit> .
forest practitioners appreciate the resistance of douglas-fir against many european pathogens
<cit> .
it is also expected that douglas-fir is better adapted to future climate conditions in central europe than e.g.
norway spruce
<cit> .
the intergovernmental panel on climate change  expects increasing summer temperatures and decreasing precipitation in central europe in the coming years
<cit> .
a trend towards warmer summers and more frequent summer droughts was reported in recent studies and simulations
<cit> .
for this reason, it is important for forest managers to select suitable tree species or provenances that are adapted to the anticipated future climate.
currently, coastal douglas-fir provenances are more frequently planted in central europe due to their superior growth performance compared to interior douglas-firs
<cit> .
the identification and characterization of differentially adapted provenances of coastal and interior douglas-fir varieties has therefore high practical value.
because of the large genome size of douglas-firs , transcriptome analysis is a cost-effective and suitable approach for the identification of candidate genes for adaptive traits and molecular markers that are linked to phenotypic variation.
transcriptomes of many species have been analyzed by next-generation sequencing technologies
<cit> , and numerous coding single-nucleotide polymorphisms  were identified in conifer species such as pinus contorta, picea glauca and pinus taeda <cit> .
douglas-firs, like other forest trees, have a high level of genetic diversity
<cit> .
for example, one study identified  <dig> snps in  <dig> candidate genes for cold-hardiness  in coastal douglas-firs
<cit> .
for this reason, transcriptome sequencing of different provenances will lead to candidate genes for differential adaptation and to many new genetic markers for the characterization of different populations.
the purpose of this study was to establish a catalogue of douglas-fir putative unique transcripts  enriched for drought stressed genes and to identify genetic polymorphisms as resource for further analysis such as resequencing projects, association studies, and gene expression profiling.
RESULTS
sequencing and assembly
the sequencing of twelve cdna libraries resulted in  <dig> , <dig> reads with an average length of  <dig> bp.
after preprocessing, the number of reads decreased to  <dig> , <dig>  read numbers were not equally distributed among libraries .
the dinm, dins and diwc libraries consisted of less than  <dig>  reads each and the diwm library of less than  <dig>  reads .
the average length of the reads decreased to  <dig> bp after pre-processing .
more than 99% of reads in each library were used for the construction of the assembly after quality trimming, with the exception of the diwm library .
a total of  <dig> , <dig>  reads were assembled into  <dig>  isotigs  of at least  <dig> bp length.
additional file
<dig> contains the origin and the number of assembled reads.
all isotigs were clustered in  <dig>  isogroups.
the mean isotig length was  <dig>  bp , the mean coverage per base was  <dig>  reads , and the mean number of reads per isotig was  <dig>  .
for  <dig>  isotigs longer than  <dig> bp, the mean coverage increased to  <dig>   reads per base.
furthermore, the mean number of reads per isotig reached  <dig>  .
length of the isotig was positively correlated with the number of reads .
d = douglas-fir, c/i = coastal/interior, n/w = needle/wood tissue, c/m/s = no/mild/severe drought stress.
based on the results of the assembly, we constructed a set of puts as outlined in the methods section.
<dig>  of  <dig>  reads with a length >  <dig> bp initially labelled as singletons were mapped to isotigs and were considered as false positive singletons.
therefore, the final put set consisted of  <dig>  sequences  with an average sequence length of  <dig>  bp .
as no reference sequence of douglas-fir was available, we used the put set as reference for the following analysis, including functional annotation and snp detection.
functional annotation of the puts
for functional annotation, we compared all puts to the ncbi nr database using blastx with an e-value cutoff of e− <dig>  at least one blast hit was obtained with  <dig>  transcripts.
if only the best hit of each transcript is considered, a total of  <dig>  different sequences  were hit in the nr database.
the largest number of hits was against picea sitchensis sequences, followed by vitis vinifera .
in the subsequent analysis, blast2go assigned at least one go term to  <dig>  transcripts.
for the three main go categories,  <dig>  transcripts were assigned a go term from the molecular function category,  <dig>  from the biological process, and  <dig>  from the cellular component category.
to compare the distribution of go terms of the douglas-fir transcriptome with the distribution of go terms of transcriptomes from other species, we also applied blast2go to the arabidopsis thaliana and picea sitchensis sequences downloaded from tair and ncbi, respectively.
we chose these two species for comparison because a. thaliana is a well studied model species with a well studied transcriptome and p. sitchensis is the species with most top blastx hits in our douglas-fir put set.
figure
<dig> and additional file
<dig> show that the distributions of go terms at go level  <dig> to  <dig> for each of the three ontology classes are highly similar for all three species.
identification of treatment-specific puts
the isotigs  were divided and clustered according to the origin of their reads .
about one third  of the isotigs contained reads from all three treatments and were therefore grouped in the cms group.
the average length of isotigs in the cms group was  <dig>  bp .
each of the cm, cs, and ms groups contained 6-8% of the isotigs with an average isotig length of  <dig> to  <dig> bp.
the c, m, and s groups contained 14-15% of the isotigs in each case.
the average lengths of those isotigs were between  <dig> and  <dig> bp.
the search for specific keywords in the blastx results revealed that  <dig>  different isotigs coming from  <dig> isogroups had a blastx hit containing one of the keywords related to stress .
<dig> of those isotigs coming from  <dig> isogroups were part of the m, s, or ms groups and will serve as top candidate genes in future studies.
we expected that cms group sequences are more conserved than sequences assigned to the remaining groups because drought stress specific sequences may evolve faster or are of a more recent evolutionary origin than common or widely expressed genes.
percentages of isotigs  with a blastx hit against the nr, ara, and picea databases within the groups .
c = control, m = mild stress, s = severe stress, cm = control and mild stress, cs = control and severe stress, ms = mild and severe stress, cms = control, mild and severe stress.
number of isotigs  with a blastx hit containing a keyword for each group .
c = control, m = mild stress, s = severe stress, cm = control and mild stress, cs = control and severe stress, ms = mild and severe stress, cms = control, mild and severe stress.
1aba = abscisic acid.
2snf = sucrose non-fermenting.
3dreb = dehydration-responsive element-binding.
4erd = early responsive to dehydration.
5cipk = cbl-interacting protein kinase.
6cdpk = calcium-dependent protein kinase.
7cbl = calcineurin b-like protein.
8pks = phytochrome kinase substrate.
to test this hypothesis, we determined the proportion of significant blastx hits within each group of isotigs against the nr, the ara, and the picea databases .
most hits were observed in the cms group  and the least number of hits in the m group .
however, there is a highly significant correlation between the average length of isotigs and percent blast hits , and also between the total sequence length of each isotig group with the proportion of blast hits .
hence, the differences in the proportion of blast hits among classes of isotigs are not a result of differential evolutionary conservation, but of the amount of sequence data in each class.
results of similarity searches with blastx and functional annotation using blast2go subdivided by transcript length in bp.
snp identification
snp detection was performed with three different programs, gsmapper, ssahasnp, and bwa/samtools, to minimize the number of false positives.
puts obtained from the assembly served as reference.
the programs detected  <dig>  ,  <dig>  , and  <dig>   snps, resulting in a total number of  <dig>  different snps.
however, only  <dig>  snps were detected by all three tools .
these snps were selected for further analysis because we consider them as most reliable true positive polymorphisms.
these snps were distributed over  <dig>  different puts of  <dig>  different isogroups.
most transcripts harbored only a single snp and  <dig>  transcripts contained more than three snps.
a total of  <dig> snps were detected in the most polymorphic put.
in the  <dig>  transcripts with at least one snp and a significant match against the nr database,  <dig>  snps were classified as synonymous and  <dig>  as non-synonymous.
in addition, we estimated the polymorphism level of the transcriptome by dividing the number of snps with the total number of different nucleotides in puts .
if only the most reliable snps are considered, the estimated nucleotide diversity  is very low.
using all snps identified by newbler, bwa/samtools and ssahasnp separately, resulted in estimated polymorphism levels of  <dig> % ,  <dig> % , and  <dig> % , respectively.
to investigate differences in the level of genetic diversity between coastal and interior douglas-firs, we divided the snps into groups depending on whether their alleles segregated in coastal or interior provenances, or in both .
the majority of snps are polymorphic in both coastal and interior provenances , but coastal provenances have a higher number of provenance-specific alleles, as seen in the comparison of ci/c  versus ci/i  groups.
partitioning of snps into groups depending on the origin  of sequence reads.
c: only reads of coastal libraries; i: only reads of interior libraries, ci: reads of coastal as well as interior libraries; ci/ci: both possible nucleotides were confirmed by reads of coastal and interior libraries; c/i: one of the possible nucleotides at the snp position was confirmed only by reads of coastal libraries, the other nucleotide only by reads of interior libraries; etc.
number of snps with a specific composition of reads.
ci/ci: variant and reference nucleotide appeared in reads from coastal and interior libraries; ci/c, ci/i: variant or reference nucleotide appeared only in reads of the coastal or interior libraries, the other one appeared in reads of both kind of libraries; c/i: variant or reference nucleotide appeared only in reads of the coastal libraries, the other one only in reads of the interior libraries; c/c, i/i: variant and reference nucleotides appeared only in reads of coastal or interior libraries.
discussion
sequencing and assembly
next generation sequencing  has now a major impact on the genome-wide analysis of transcriptomes in non-model species
<cit> .
to achieve a comprehensive characterization of the protein-coding genome of douglas-fir, we exposed young seedlings from different provenances to drought stress treatments and generated normalized cdna libraries to enrich for rare transcripts or genes not constitutively expressed.
all libraries were assembled into a single assembly to maximize the number of reads per transcript and to improve the quality of assembly and snp annotation.
there is a strong relationship between the number of reads and the length of a transcript, confirming the observation that longer transcripts consist of more reads than shorter transcripts
<cit> .
the number and average length of the reads of four libraries  were below expectation
<cit> , probably because of problems during the sequencing process.
however, we did not exclude these libraries, because they contributed the same proportion of reads to the assembly as the other libraries .
<dig> % of all reads were assembled into isotigs during assembly, which is a high proportion compared to similar  <dig> transcriptome assemblies.
for example, 88% of reads were assembled in melitaea cinxia <cit>  and eucalyptus grandis <cit> , 78% in pandinus imperator <cit> , and 48% in pinus contorta <cit> .
one cause for the high proportion in our data is the stringent preprocessing of reads, which excluded most uninformative reads prior to the assembly.
the number of puts in the assembly  exceeds the number of expected genes in conifer genomes, which ranges from  <dig>  to  <dig>  genes
<cit> .
nevertheless, the put set is smaller than the one obtained with pinus contorta
<cit> , but larger than in other  <dig> transcriptome sequencing projects .
it should be noted that it is difficult to compare numbers of transcripts in different projects, because they are influenced by the genome and transcriptome sizes of the sequenced organism, the assembly method used
<cit> , and the number of reads used for the assembly.
overall, the high number of transcripts compared to the expected number of genes is likely a result of incompletely assembled genes.
the average length of transcripts is  <dig>  bp , which is approximately half of the expected average gene length in eukaryotes
<cit> .
functional annotation of puts
we annotated the puts by using stringent blastx searches against the nr database from ncbi with a cutoff e-value of e− <dig>  assuming that each unique blast hit represents a different transcribed gene, we tagged  <dig>  genes of the douglas-fir genome.
this number is similar to other projects in pinus contorta with  <dig>  tagged genes
<cit>  and is lower than the estimated total number of about  <dig>  transcribed genes in white spruce picea glauca <cit> .
if we further assume that the number of unique blast hits equals the number of transcribed genes and that the transcriptome sizes of douglas-fir and white spruce are comparable, it seems that the put set generated in this study tags about two thirds of the protein-coding genes of douglas-fir.
the missing third likely comprise  cdnas that were excluded from assembly because of low quality;  genes that are expressed at different developmental stages, growth conditions, or tissues and were thus not represented in the twelve cdna libraries despite the normalization process; and  non-conserved genes, which are either lineage-specific and not yet present in the nr database or rapidly evolving genes with e-values > e− <dig> in the blastx comparisons.
since about 75% of the transcripts do not result in significant blastx hits, one may speculate that a large proportion represents non-conserved genes.
gene ontologies  provide a standardized set of terms to describe genes and gene products consistently in different species and databases
<cit> .
go terms are widely used for annotation and for comparisons of gene products of different species .
the similarity of the go annotation distributions of douglas-fir puts to the well-characterized transcriptome of a. thaliana and the protein sequences of p. sitchensis  suggests that the put set broadly represents the douglas-fir transcriptome and can be viewed as being representative for further applications and investigations.
identification of drought stress related genes
dividing puts consisting of multiple reads  by the origin of their reads is a simple, but useful method to identify potential treatment-specific sequences.
about 50% of isotigs consist of reads from the single treatment groups c, m, or s. on average, they are shorter than isotigs containing reads from at least two different treatments.
the isotigs of the m and s groups, but also the ms group, were considered to be top candidates for drought stress tolerance or resistance.
however, since most treatment-specific isotigs consist of only two or three reads that originated from a total of twelve cdna libraries, we reasoned that the power of a statistical model to detect significant differences is low.
therefore, we compared the extent of evolutionary conservation between groups as judged by the proportion of significant blast hits.
this analysis is based on the notion that widely expressed genes are under stronger selective constraint than treatment-specific genes
<cit> .
under the assumption that constitutively expressed genes are expressed in all different treatments, we expected that cms group isotigs are more conserved than isotigs from the c, m, and s groups.
since the libraries were normalized and cdna levels do not represent true expression levels, we restricted our analysis to presence-absence patterns.
the differences in the fractions of blastx hits in single treatment groups and the cms group suggested that genes expressed in all three treatments are more conserved.
however, if groups are corrected for the total amount of sequence data, cms group isotigs are not more conserved than treatment-specific isotigs, because the main determinant for a blast hit is isotig length which is longer in cms isotigs .
this pattern was also observed in white spruce
<cit> .
in addition to testing the general hypothesis that treatment-specific genes are less conserved than widely-expressed genes, we also parsed blastx results for drought stress related keywords to find potential drought stress related puts.
we expect that the  <dig>  transcripts with a blastx hit containing one of the keywords are probably involved in the douglas-firs response to drought .
more than  <dig>  candidate puts are part of the cms group and only  <dig> candidates are part of the m, s, and ms groups.
this reflects that the response to drought seems to be mainly facilitated through changes in gene expression levels via up- or down-regulation.
the small set of  <dig> puts exclusively induced by drought stress appears to contribute to a specific drought response, but this needs to be further verified because their expression pattern may reflect a sampling artifact.
even though the function of those puts may not be conserved across large evolutionary distances, the identified puts serve as top candidates for further analysis of sequence and expression variation in comparisons of differentially adapted  douglas-fir provenances.
analysis of genetic variation
the construction of the cdna libraries representing genetically different provenances allowed the detection of snps for later analysis of patterns of genetic variation.
the two most important results are the high proportion of shared polymorphisms and the strong influence of the snp calling algorithm on the number of detected snps.
by taking a conservative approach and considering only those snps that were called by all three programs, only  <dig>   snps were obtained, which is only about half of the number detected with gsmapper, which identified the lowest number of snps .
since the numbers of called snps differed highly between snp detection tools, our results indicate that snp calling from next generation sequencing data depend to a high degree on the software tools used.
therefore, results should be interpreted with caution, if relying on a single snp detection approach only.
to our knowledge there are no systematic studies yet that compared the accuracy of different snp callers with next-generation sequencing data.
a comparison of the snp density  of the most reliable snps with published data shows that the former is an underestimate of the true level of sequence variation in douglas-fir.
the snp density is  <dig> snp for every  <dig>  bp, whereas other studies estimated an average snp density from  <dig> snp per  <dig> bp to  <dig> snp per  <dig> bp using sanger sequencing protocols
<cit> .
the reasons for the large difference to the reported snp density are probably the stringency criteria used and the better quality of base-calling using sanger sequencing.
if we take only the snps identified by bwa/samtools or ssahasnp in account, the snp density increases to  <dig> snp per  <dig> bp and  <dig> snp per  <dig> bp, respectively.
nevertheless, our sequence data make a significant contribution to the number of douglas-fir snps available for further applications.
until now, only around  <dig>  snps have been identified in douglas-fir
<cit> .
if only the most reliable snps are considered, a key result is the large number of snps whose alleles are segregating in both the coastal and interior provenances .
in only 5% of snps  the two alleles are specific to coastal and interior provenances, respectively.
this high proportion of shared polymorphisms indicates either a high level of shared ancestral polymorphisms between the two main douglas-fir varieties, or recent, possibly pollen-mediated gene flow.
the comparison of snps, in which only one of the two alleles is shared between coastal and interior provenances suggest a higher level of genetic diversity in coastal provenances because three times as many snps are polymorphic for both alleles in the coastal  than in the interior accessions .
this difference is also observed for snps which were called only in either the interior or coastal provenances because no reads were available from the other provenance, respectively .
although these results are consistent with earlier studies on the genetic diversity of douglas-fir varieties
<cit> , they are also certainly influenced by the different numbers of reads originating from coastal and interior cdna libraries .
since there are 70% more reads from the coastal provenances, the probability of finding a polymorphism in these provenances is increased and needs to be accounted for in further conclusions.
different numbers of reads can be accounted for by using methods for population genetic inference developed for next-generation sequencing that account for differences in read numbers from individuals or pools of individuals in estimating allele frequencies and population parameters
<cit> .
however, such an approach does not work in the present study because allele frequencies depend on the sampled individuals in a library, the gene expression level and the effect of normalization on read numbers.
unbiased population genetic estimators like tajima’s Π can be calculated from  <dig> data
<cit> , but as the coverage at most snp positions is much smaller than the total number of individuals, the results are not reliable.
the development of genotyping and resequencing arrays using the present set of puts to estimate snp allele frequencies and population genetic inference will allow accurate and unbiased estimates of nucleotide diversity.
CONCLUSIONS
in this study we established a catalogue of douglas-fir putative unique transcripts  enriched for drought stress induced genes.
although the real magnitude of the transcriptome is yet unknown, we estimate that the majority of the transcriptome has been tagged by the put set presented here.
this is based on the results of the functional annotation and the comparison of the go term distributions with those of arabidopsis thaliana and picea sitchensis.
by analyzing sequence variation in the puts we detected thousands of new snps.
furthermore, we identified drought stress specific candidate sequences.
taken together these data represent a useful resource for the next steps in the characterization of the douglas-fir genome and transcriptome and the association of genetic variation with phenotypic traits such as adaptation to different ecogeographic environments.

 douglas-fir franco) is an ecologically highly variable species that occurs in two main varieties in north america. the natural range of the coastal or green douglas-fir extends over km from the pacific northwestern fog belt and the adjacent summer-dry coastal range and cascade mountains to the drier coastland of southern california. the interior or blue douglas-fir is distributed over more than km along the dry continental climates of the montane to the subalpine rocky mountains from alberta to colorado with isolated populations reaching into mexico. douglas-fir grows from sea level on vancouver island up to m altitude in the southern rocky mountains. within its natural range, douglas-fir has evolved into a variety of genetically diverse populations adapted to contrasting ecozones. douglas-fir populations differ in their response to frost, drought, and along environmental gradients. like most conifer species, it is able to cope with limitations in soil water availability within its natural range. there is a negative relationship between shoot water potential and the photosynthesis rate, which decreased by about 70% in water-stressed trees with a pre-dawn shoot water potential of about - mpa. in conifers such as douglas-fir or pine, the recovery of photosynthesis upon rainfall and rewatering occurs within days together with the rapid recovery of predawn shoot water potential from stressed, or mildly stressed to values higher than - mpa. this high ecological, genetical and physiological diversity provides an excellent system to study the adaptation of conifer trees to contrasting environments. due to its rapid growth and favorable wood quality, douglas-fir is an economically relevant species for timber production. in europe, the area of douglas-fir production is rising. forest practitioners appreciate the resistance of douglas-fir against many european pathogens. it is also expected that douglas-fir is better adapted to future climate conditions in central europe than e.g. norway spruce. the intergovernmental panel on climate change expects increasing summer temperatures and decreasing precipitation in central europe in the coming years. a trend towards warmer summers and more frequent summer droughts was reported in recent studies and simulations. for this reason, it is important for forest managers to select suitable tree species or provenances that are adapted to the anticipated future climate. currently, coastal douglas-fir provenances are more frequently planted in central europe due to their superior growth performance compared to interior douglas-firs. the identification and characterization of differentially adapted provenances of coastal and interior douglas-fir varieties has therefore high practical value. because of the large genome size of douglas-firs, transcriptome analysis is a cost-effective and suitable approach for the identification of candidate genes for adaptive traits and molecular markers that are linked to phenotypic variation. transcriptomes of many species have been analyzed by next-generation sequencing technologies, and numerous coding single-nucleotide polymorphisms were identified in conifer species such as pinus contorta, picea glauca and pinus taeda. douglas-firs, like other forest trees, have a high level of genetic diversity. for example, one study identified snps in candidate genes for cold-hardiness in coastal douglas-firs. for this reason, transcriptome sequencing of different provenances will lead to candidate genes for differential adaptation and to many new genetic markers for the characterization of different populations. the purpose of this study was to establish a catalogue of douglas-fir putative unique transcripts enriched for drought stressed genes and to identify genetic polymorphisms as resource for further analysis such as resequencing projects, association studies, and gene expression profiling.  sequencing and assembly the sequencing of twelve cdna libraries resulted in, reads with an average length of bp. after preprocessing, the number of reads decreased to, read numbers were not equally distributed among libraries. the dinm, dins and diwc libraries consisted of less than reads each and the diwm library of less than reads. the average length of the reads decreased to bp after pre-processing. more than 99% of reads in each library were used for the construction of the assembly after quality trimming, with the exception of the diwm library. a total of, reads were assembled into isotigs of at least bp length. additional file contains the origin and the number of assembled reads. all isotigs were clustered in isogroups. the mean isotig length was bp, the mean coverage per base was reads, and the mean number of reads per isotig was. for isotigs longer than bp, the mean coverage increased to reads per base. furthermore, the mean number of reads per isotig reached. length of the isotig was positively correlated with the number of reads. d = douglas-fir, c/i = coastal/interior, n/w = needle/wood tissue, c/m/s = no/mild/severe drought stress. based on the results of the assembly, we constructed a set of puts as outlined in the methods section. of reads with a length >  bp initially labelled as singletons were mapped to isotigs and were considered as false positive singletons. therefore, the final put set consisted of sequences with an average sequence length of bp. as no reference sequence of douglas-fir was available, we used the put set as reference for the following analysis, including functional annotation and snp detection. functional annotation of the puts for functional annotation, we compared all puts to the ncbi nr database using blastx with an e-value cutoff of e− at least one blast hit was obtained with transcripts. if only the best hit of each transcript is considered, a total of different sequences were hit in the nr database. the largest number of hits was against picea sitchensis sequences, followed by vitis vinifera. in the subsequent analysis, blast2go assigned at least one go term to transcripts. for the three main go categories, transcripts were assigned a go term from the molecular function category, from the biological process, and from the cellular component category. to compare the distribution of go terms of the douglas-fir transcriptome with the distribution of go terms of transcriptomes from other species, we also applied blast2go to the arabidopsis thaliana and picea sitchensis sequences downloaded from tair and ncbi, respectively. we chose these two species for comparison because a. thaliana is a well studied model species with a well studied transcriptome and p. sitchensis is the species with most top blastx hits in our douglas-fir put set. figure and additional file show that the distributions of go terms at go level to for each of the three ontology classes are highly similar for all three species. identification of treatment-specific puts the isotigs were divided and clustered according to the origin of their reads. about one third of the isotigs contained reads from all three treatments and were therefore grouped in the cms group. the average length of isotigs in the cms group was bp. each of the cm, cs, and ms groups contained 6-8% of the isotigs with an average isotig length of to bp. the c, m, and s groups contained 14-15% of the isotigs in each case. the average lengths of those isotigs were between and bp. the search for specific keywords in the blastx results revealed that different isotigs coming from isogroups had a blastx hit containing one of the keywords related to stress. of those isotigs coming from isogroups were part of the m, s, or ms groups and will serve as top candidate genes in future studies. we expected that cms group sequences are more conserved than sequences assigned to the remaining groups because drought stress specific sequences may evolve faster or are of a more recent evolutionary origin than common or widely expressed genes. percentages of isotigs with a blastx hit against the nr, ara, and picea databases within the groups. c = control, m = mild stress, s = severe stress, cm = control and mild stress, cs = control and severe stress, ms = mild and severe stress, cms = control, mild and severe stress. number of isotigs with a blastx hit containing a keyword for each group. c = control, m = mild stress, s = severe stress, cm = control and mild stress, cs = control and severe stress, ms = mild and severe stress, cms = control, mild and severe stress. 1aba = abscisic acid. 2snf = sucrose non-fermenting. 3dreb = dehydration-responsive element-binding. 4erd = early responsive to dehydration. 5cipk = cbl-interacting protein kinase. 6cdpk = calcium-dependent protein kinase. 7cbl = calcineurin b-like protein. 8pks = phytochrome kinase substrate. to test this hypothesis, we determined the proportion of significant blastx hits within each group of isotigs against the nr, the ara, and the picea databases. most hits were observed in the cms group and the least number of hits in the m group. however, there is a highly significant correlation between the average length of isotigs and percent blast hits, and also between the total sequence length of each isotig group with the proportion of blast hits. hence, the differences in the proportion of blast hits among classes of isotigs are not a result of differential evolutionary conservation, but of the amount of sequence data in each class. results of similarity searches with blastx and functional annotation using blast2go subdivided by transcript length in bp. snp identification snp detection was performed with three different programs, gsmapper, ssahasnp, and bwa/samtools, to minimize the number of false positives. puts obtained from the assembly served as reference. the programs detected,, and snps, resulting in a total number of different snps. however, only snps were detected by all three tools. these snps were selected for further analysis because we consider them as most reliable true positive polymorphisms. these snps were distributed over different puts of different isogroups. most transcripts harbored only a single snp and transcripts contained more than three snps. a total of snps were detected in the most polymorphic put. in the transcripts with at least one snp and a significant match against the nr database, snps were classified as synonymous and as non-synonymous. in addition, we estimated the polymorphism level of the transcriptome by dividing the number of snps with the total number of different nucleotides in puts. if only the most reliable snps are considered, the estimated nucleotide diversity is very low. using all snps identified by newbler, bwa/samtools and ssahasnp separately, resulted in estimated polymorphism levels of %, %, and %, respectively. to investigate differences in the level of genetic diversity between coastal and interior douglas-firs, we divided the snps into groups depending on whether their alleles segregated in coastal or interior provenances, or in both. the majority of snps are polymorphic in both coastal and interior provenances, but coastal provenances have a higher number of provenance-specific alleles, as seen in the comparison of ci/c versus ci/i groups. partitioning of snps into groups depending on the origin of sequence reads. c: only reads of coastal libraries; i: only reads of interior libraries, ci: reads of coastal as well as interior libraries; ci/ci: both possible nucleotides were confirmed by reads of coastal and interior libraries; c/i: one of the possible nucleotides at the snp position was confirmed only by reads of coastal libraries, the other nucleotide only by reads of interior libraries; etc. number of snps with a specific composition of reads. ci/ci: variant and reference nucleotide appeared in reads from coastal and interior libraries; ci/c, ci/i: variant or reference nucleotide appeared only in reads of the coastal or interior libraries, the other one appeared in reads of both kind of libraries; c/i: variant or reference nucleotide appeared only in reads of the coastal libraries, the other one only in reads of the interior libraries; c/c, i/i: variant and reference nucleotides appeared only in reads of coastal or interior libraries. discussion sequencing and assembly next generation sequencing has now a major impact on the genome-wide analysis of transcriptomes in non-model species. to achieve a comprehensive characterization of the protein-coding genome of douglas-fir, we exposed young seedlings from different provenances to drought stress treatments and generated normalized cdna libraries to enrich for rare transcripts or genes not constitutively expressed. all libraries were assembled into a single assembly to maximize the number of reads per transcript and to improve the quality of assembly and snp annotation. there is a strong relationship between the number of reads and the length of a transcript, confirming the observation that longer transcripts consist of more reads than shorter transcripts. the number and average length of the reads of four libraries were below expectation, probably because of problems during the sequencing process. however, we did not exclude these libraries, because they contributed the same proportion of reads to the assembly as the other libraries. % of all reads were assembled into isotigs during assembly, which is a high proportion compared to similar transcriptome assemblies. for example, 88% of reads were assembled in melitaea cinxia and eucalyptus grandis, 78% in pandinus imperator, and 48% in pinus contorta. one cause for the high proportion in our data is the stringent preprocessing of reads, which excluded most uninformative reads prior to the assembly. the number of puts in the assembly exceeds the number of expected genes in conifer genomes, which ranges from to genes. nevertheless, the put set is smaller than the one obtained with pinus contorta, but larger than in other transcriptome sequencing projects. it should be noted that it is difficult to compare numbers of transcripts in different projects, because they are influenced by the genome and transcriptome sizes of the sequenced organism, the assembly method used, and the number of reads used for the assembly. overall, the high number of transcripts compared to the expected number of genes is likely a result of incompletely assembled genes. the average length of transcripts is bp, which is approximately half of the expected average gene length in eukaryotes. functional annotation of puts we annotated the puts by using stringent blastx searches against the nr database from ncbi with a cutoff e-value of e− assuming that each unique blast hit represents a different transcribed gene, we tagged genes of the douglas-fir genome. this number is similar to other projects in pinus contorta with tagged genes and is lower than the estimated total number of about transcribed genes in white spruce picea glauca. if we further assume that the number of unique blast hits equals the number of transcribed genes and that the transcriptome sizes of douglas-fir and white spruce are comparable, it seems that the put set generated in this study tags about two thirds of the protein-coding genes of douglas-fir. the missing third likely comprise cdnas that were excluded from assembly because of low quality; genes that are expressed at different developmental stages, growth conditions, or tissues and were thus not represented in the twelve cdna libraries despite the normalization process; and non-conserved genes, which are either lineage-specific and not yet present in the nr database or rapidly evolving genes with e-values > e− in the blastx comparisons. since about 75% of the transcripts do not result in significant blastx hits, one may speculate that a large proportion represents non-conserved genes. gene ontologies provide a standardized set of terms to describe genes and gene products consistently in different species and databases. go terms are widely used for annotation and for comparisons of gene products of different species. the similarity of the go annotation distributions of douglas-fir puts to the well-characterized transcriptome of a. thaliana and the protein sequences of p. sitchensis suggests that the put set broadly represents the douglas-fir transcriptome and can be viewed as being representative for further applications and investigations. identification of drought stress related genes dividing puts consisting of multiple reads by the origin of their reads is a simple, but useful method to identify potential treatment-specific sequences. about 50% of isotigs consist of reads from the single treatment groups c, m, or s. on average, they are shorter than isotigs containing reads from at least two different treatments. the isotigs of the m and s groups, but also the ms group, were considered to be top candidates for drought stress tolerance or resistance. however, since most treatment-specific isotigs consist of only two or three reads that originated from a total of twelve cdna libraries, we reasoned that the power of a statistical model to detect significant differences is low. therefore, we compared the extent of evolutionary conservation between groups as judged by the proportion of significant blast hits. this analysis is based on the notion that widely expressed genes are under stronger selective constraint than treatment-specific genes. under the assumption that constitutively expressed genes are expressed in all different treatments, we expected that cms group isotigs are more conserved than isotigs from the c, m, and s groups. since the libraries were normalized and cdna levels do not represent true expression levels, we restricted our analysis to presence-absence patterns. the differences in the fractions of blastx hits in single treatment groups and the cms group suggested that genes expressed in all three treatments are more conserved. however, if groups are corrected for the total amount of sequence data, cms group isotigs are not more conserved than treatment-specific isotigs, because the main determinant for a blast hit is isotig length which is longer in cms isotigs. this pattern was also observed in white spruce. in addition to testing the general hypothesis that treatment-specific genes are less conserved than widely-expressed genes, we also parsed blastx results for drought stress related keywords to find potential drought stress related puts. we expect that the transcripts with a blastx hit containing one of the keywords are probably involved in the douglas-firs response to drought. more than candidate puts are part of the cms group and only candidates are part of the m, s, and ms groups. this reflects that the response to drought seems to be mainly facilitated through changes in gene expression levels via up- or down-regulation. the small set of puts exclusively induced by drought stress appears to contribute to a specific drought response, but this needs to be further verified because their expression pattern may reflect a sampling artifact. even though the function of those puts may not be conserved across large evolutionary distances, the identified puts serve as top candidates for further analysis of sequence and expression variation in comparisons of differentially adapted douglas-fir provenances. analysis of genetic variation the construction of the cdna libraries representing genetically different provenances allowed the detection of snps for later analysis of patterns of genetic variation. the two most important results are the high proportion of shared polymorphisms and the strong influence of the snp calling algorithm on the number of detected snps. by taking a conservative approach and considering only those snps that were called by all three programs, only snps were obtained, which is only about half of the number detected with gsmapper, which identified the lowest number of snps. since the numbers of called snps differed highly between snp detection tools, our results indicate that snp calling from next generation sequencing data depend to a high degree on the software tools used. therefore, results should be interpreted with caution, if relying on a single snp detection approach only. to our knowledge there are no systematic studies yet that compared the accuracy of different snp callers with next-generation sequencing data. a comparison of the snp density of the most reliable snps with published data shows that the former is an underestimate of the true level of sequence variation in douglas-fir. the snp density is snp for every bp, whereas other studies estimated an average snp density from snp per bp to snp per bp using sanger sequencing protocols. the reasons for the large difference to the reported snp density are probably the stringency criteria used and the better quality of base-calling using sanger sequencing. if we take only the snps identified by bwa/samtools or ssahasnp in account, the snp density increases to snp per bp and snp per bp, respectively. nevertheless, our sequence data make a significant contribution to the number of douglas-fir snps available for further applications. until now, only around snps have been identified in douglas-fir. if only the most reliable snps are considered, a key result is the large number of snps whose alleles are segregating in both the coastal and interior provenances. in only 5% of snps the two alleles are specific to coastal and interior provenances, respectively. this high proportion of shared polymorphisms indicates either a high level of shared ancestral polymorphisms between the two main douglas-fir varieties, or recent, possibly pollen-mediated gene flow. the comparison of snps, in which only one of the two alleles is shared between coastal and interior provenances suggest a higher level of genetic diversity in coastal provenances because three times as many snps are polymorphic for both alleles in the coastal than in the interior accessions. this difference is also observed for snps which were called only in either the interior or coastal provenances because no reads were available from the other provenance, respectively. although these results are consistent with earlier studies on the genetic diversity of douglas-fir varieties, they are also certainly influenced by the different numbers of reads originating from coastal and interior cdna libraries. since there are 70% more reads from the coastal provenances, the probability of finding a polymorphism in these provenances is increased and needs to be accounted for in further conclusions. different numbers of reads can be accounted for by using methods for population genetic inference developed for next-generation sequencing that account for differences in read numbers from individuals or pools of individuals in estimating allele frequencies and population parameters. however, such an approach does not work in the present study because allele frequencies depend on the sampled individuals in a library, the gene expression level and the effect of normalization on read numbers. unbiased population genetic estimators like tajima’s Π can be calculated from data, but as the coverage at most snp positions is much smaller than the total number of individuals, the results are not reliable. the development of genotyping and resequencing arrays using the present set of puts to estimate snp allele frequencies and population genetic inference will allow accurate and unbiased estimates of nucleotide diversity.  in this study we established a catalogue of douglas-fir putative unique transcripts enriched for drought stress induced genes. although the real magnitude of the transcriptome is yet unknown, we estimate that the majority of the transcriptome has been tagged by the put set presented here. this is based on the results of the functional annotation and the comparison of the go term distributions with those of arabidopsis thaliana and picea sitchensis. by analyzing sequence variation in the puts we detected thousands of new snps. furthermore, we identified drought stress specific candidate sequences. taken together these data represent a useful resource for the next steps in the characterization of the douglas-fir genome and transcriptome and the association of genetic variation with phenotypic traits such as adaptation to different ecogeographic environments. 
the two alleles of about 60% of the snps are segregating simultaneously in both coastal and interior provenances, which indicates a high proportion of ancestral shared polymorphisms or a high level of gene flow between these two ecologically and phenotypically different varieties. both will serve as a useful resource for the further characterization of the genome and transcriptome of douglas-fir and for the analysis of genetic variation using genotyping or resequencing methods. however, genomic data and tools for studying genetic variation in natural populations to understand the genetic and physiological mechanisms of adaptation are currently missing for douglas-fir. functional annotation by blast searches and gene ontology mapping showed that the composition of functional classes is very similar to other plant transcriptomes and demonstrated that a large fraction of the douglas-fir transcriptome is tagged by the puts. the sequencing of twelve normalized cdna libraries resulted in million reads from which a set of putative unique transcripts was assembled. a total number of single nucleotide polymorphisms were detected by three snp detection tools. the libraries were constructed from needle and wood tissue of coastal and interior provenances subjected to drought stress experiments. this study represents a first step towards characterizing the douglas-fir transcriptome based on sequencing of twelve cdna libraries. however, only snps were identified by all three methods, indicating that snp detection depends on the particular method used. based on evolutionary conservation, we identified about candidate genes related to drought stress.
 douglas-fir extends over a wide range of contrasting environmental conditions, reflecting substantial local adaptation. for this reason, it is an interesting model species to study plant adaptation and the effects of global climate change such as increased temperatures and significant periods of drought on individual trees and the forest landscape in general. however, genomic data and tools for studying genetic variation in natural populations to understand the genetic and physiological mechanisms of adaptation are currently missing for douglas-fir. this study represents a first step towards characterizing the douglas-fir transcriptome based on sequencing of twelve cdna libraries. the libraries were constructed from needle and wood tissue of coastal and interior provenances subjected to drought stress experiments.  the sequencing of twelve normalized cdna libraries resulted in million reads from which a set of putative unique transcripts was assembled. functional annotation by blast searches and gene ontology mapping showed that the composition of functional classes is very similar to other plant transcriptomes and demonstrated that a large fraction of the douglas-fir transcriptome is tagged by the puts. based on evolutionary conservation, we identified about candidate genes related to drought stress. a total number of single nucleotide polymorphisms were detected by three snp detection tools. however, only snps were identified by all three methods, indicating that snp detection depends on the particular method used. the two alleles of about 60% of the snps are segregating simultaneously in both coastal and interior provenances, which indicates a high proportion of ancestral shared polymorphisms or a high level of gene flow between these two ecologically and phenotypically different varieties.  we established a catalogue of puts and large snp database for douglas-fir. both will serve as a useful resource for the further characterization of the genome and transcriptome of douglas-fir and for the analysis of genetic variation using genotyping or resequencing methods. 