BACKGROUND
studies of molecular evolution have until recently focused on nucleotide divergence, while studies of variation in gene expression profiles have mainly been restricted to a few model species such as drosophila and mice  <cit> .
this is because the technologies for studying gene expression have not been available  for non-model species  <cit> .
however, sequencing-based technologies for expression profiling can now be utilised to this end.
by counting the number of reads generated by sequencing of cdna from different genes in the transcriptome, one can get an estimate of the expression level of these genes in the particular tissues sampled  <cit> .
a complementary approach is to scan publicly available databases of expressed sequence tags  for the genes of interest.
in addition to microarrays, these strategies, called digital transcriptomics, are today the most commonly used methods for investigating expression patterns  <cit> .
digital transcriptomics has received a great deal of attention, but the use of these methods has been restricted in many species by the requirement of having a reference genome to evaluate and analyse the data.
the advent of massively parallel  sequencing is now starting to change this picture by providing a cost-effective way of generating large amount of sequence data in species where there is no prior knowledge of the genome sequence  <cit> .
next-generation sequencing technology generally generates millions of short sequence reads, each read being tens to hundreds of base pairs long, depending on the specific platform.
this enables detection of genes even with very low expression levels.
roche 454-sequencing  <cit> , in particular, generates reads that are long enough to be informative in the absence of a reference genome  <cit> .
here, we evaluate the use of 454-sequencing to investigate tissue specific gene expression profiles.
next-generation sequencing can be used to not only describe genome-wide patterns of gene expression, but also to characterise specific gene families or genetic pathways.
to illustrate this point, we use the ecologically important and widely studied genes of the major histocompatibility complex  for a more detailed analysis.
these genes are a very common focus of studies that take a candidate gene approach to investigate functionally important genetic variation in immune function  <cit> .
mhc genes are among the most variable of the vertebrate genomes  <cit> .
in particular, the classical mhc genes  exhibit an extraordinary level of polymorphism.
this polymorphism is strongly associated to the role of these genes in regulating and triggering the adaptive immune response.
studies have found links between nucleic acid variation in the mhc genes and resistance to parasites  <cit> , sexually selected ornaments  <cit> , mate choice  <cit> , maternal-foetal incompatibilities  <cit>  and local adaptation  <cit> .
typically studies of mhc variation have focused on sequence variation only in a few highly polymorphic regions of class i and class ii genes, while variation in other genes, regions and expression levels has largely been ignored.
the completion of the genome sequence of first the chicken   <cit>  and now of the zebra finch  have opened the door for in-depth studies of organisation and expression of mhc genes in birds.
there are striking differences in the way the adaptive immune defence operates in birds compared to mammals  <cit>  and it could be envisioned that such studies will reveal new insights in the evolution of vertebrate immunity.
the aim of the present study was to investigate tissue-specific gene expression patterns in the zebra finch.
with the sequencing of its genome, the zebra finch has taken a major step towards becoming an important model system for bird genomics  <cit> .
outside of some recent studies of gene expression in brain  <cit> , however, little is known about genome-scale, and organism-wide patterns of gene expression in song birds.
in this study we describe patterns of gene expression across six zebra finch tissues and explore the relationship between expression profiles of genes and characteristics of their molecular evolution.
to this end, we use a next-generation sequencing  digital transcriptomics approach known as rna-seq  <cit> .
this methodology was recently employed to study gene expression differentiation between two subspecies of crow   <cit> , but as far as we are aware, this is the first time that a bird transcriptome has been characterised in multiple tissues using an ngs rna-seq approach.
in addition to global patterns of gene expression, we highlight patterns of expression in the genes of the mhc.
because of the complex history of duplication among certain mhc genes, gene expression profiles have the potential to offer insight into the evolutionary fates of these duplicated genes.
importantly, characterizing the expression of mhc genes will also facilitate downstream studies of these genes in ecological contexts by identifying functionally important loci.
RESULTS
assembly of  <dig> sequencing reads
after trimming and removal of contaminant sequences a total of  <dig> , <dig> reads were available, with a mean read length of  <dig> nucleotides.
<dig>  of these reads were assembled  de-novo  into  <dig>  contigs with a mean contig length of  <dig> nucleotides  and a mean of  <dig> reads per contig.
the total length of all contigs was  <dig>  kb.
for read and contig statistics for each tissue separately see additional file 1: appendix s <dig>   <dig>  of the contigs showed signatures of multiple splice variants, as indicated by gaps in alignments between the contig and one or more of the reads that contribute to that contig.
expression levels for contigs
expression levels were highly variable between tissues and contigs.
a vast majority of contigs were made up by only a few reads  but some had indications of very high expression levels .
as would be expected  there was a positive relationship between contig length  and contig depth .
the tissue specificity of the expression  of contigs was negatively correlated with the overall  amount of expression .
this negative correlation could be the result of a sampling artefact during the calculation of τ. we found however, that the observed correlation was significantly  stronger than the mean simulated correlation  based on unbiased τ-values calculated from randomisations of re-sampled data.
this suggests that there is indeed a sampling bias in the calculation of τ but that it is not strong enough to alone explain our observed correlation.
interestingly, τ seems to have a somewhat bimodal distribution  with peaks around  <dig>  and  <dig> , indicating that the genes may group into two different classes of tissue specificity.
however, only a few contigs showed evidence of very low τ-values indicative of housekeeping genes.
outlier contigs with high expression levels
three contigs were found to have strikingly high overall expression levels; following blast searches against the chicken and zebra finch gene databases they were found to represent albumin , heat shock  <dig> kda protein  <dig> beta  and nadh dehydrogenase subunit  <dig> .
all of these are considered to be so-called housekeeping genes  and are also highly expressed in mammals.
some contigs were conspicuous in having very strong expression in one or few tissues.
genes represented by these include elongation factor 1-alpha , cytoplasmic beta-actin gene , haemoglobin alpha  and mhc class ii associated invariant chain ii .
one contig was found in high levels in testes  but was almost completely absent in other tissues, the blast search revealed that this originates from a contamination with dna from a freshwater planarian .
this is likely to have occurred in the laboratory that carried out the sequencing , since the genome of schmidtea mediterranea was being sequenced there at the same time as the zebra finch cdna preparation.
this contig, together with  <dig> other contigs resulting from contamination , was removed from the data before conducting downstream analyses.
coverage of the zebra finch transcriptome
<dig>  contigs from the de-novo assembly and  <dig>  of the non-assembled singletons gave significant blast hits against at least one predicted zebra finch gene.
since the contigs were generally much shorter than the total cdna length of the gene it was commonly found that several different contigs matched the same gene.
in total  <dig>  zebra finch transcripts present in the biomart database were found to correspond to the 454/est transcriptome contigs and singletons.
this represents 65% of the total characterised zebra finch transctiptome .
the transcripts are derived from  <dig>  different genes, suggesting that more than one splice variant was detected  for ~2% of the genes.
we also identified potential novel splice variants for  <dig> of the expressed zebra finch genes, as indicated by gaps in the alignments of the contig and the gene prediction.
on average 38% of the lengths of represented transcripts were covered by contig sequences and  <dig> transcripts were fully covered.
to further investigate the extent of transcriptomic coverage, we investigated the presence of known genes in various metabolic pathways and signalling cascades  <cit> .
for the metabolic pathways about 85% of the genes were represented and for signalling cascades we found around 60%.
<dig>   of all genes found were expressed in all investigated tissues and  <dig>   were expressed exclusively in one tissue .
out of the  <dig>  contigs that did not give any matches to known predicted zebra finch transcripts, most  still gave highly significant blast hits  against the zebra finch genome sequence, suggesting that these represent transcribed regions that have not yet been annotated.
the remaining  <dig>  contigs  may represent genes in regions of the genome that have not been sequenced and/or assembled in the current genome assembly, or additional contamination from other organisms that are not represented in genbank.
number of mhc related genes included in this study is also given.
the mean index of tissue specificity of expression  and its 95% ci  for each pathway is also given.
nmax represents the number of genes with maximal expression in each of the tissues and the number within brackets  is the number of genes which are expressed uniquely in that tissue.
analyses of gene expression profiles
we found a positive correlation between tissue specificity of expression  and the ratio of non-synonymous to synonymous substitution rate  when compared to the chicken orthologue of the gene in question .
there was also a negative correlation between total expression level of the gene and ω .
there was a weak positive correlation between the length of the gene and the total level of gene expression , and a negative correlation between gene length and τ .
there were differences in ω between the tissues in which genes were primarily expressed .
genes that were primarily expressed in the embryo had the lowest mean ω-value .
the expression specificity  of genes also varied significantly between tissues of maximal expression .
the lowest tissue specificity was found in genes with primary expression in embryo and muscle, while the highest τ was found in genes with maximal expression in skin and testes .
analysis of expression in relation to go-terms
there were  <dig> gene ontology  terms overrepresented  in genes with high levels of expression specificity .
these represent processes such as cellular and organelle movement and specific enzymatic processes .
some go terms overrepresented in genes with high tissue specificity are associated with reproduction  and immune defence .
genes identified as having low tissue specificity of gene expression were significantly overrepresented for  <dig> different go terms .
these terms generally represented functions such as protein synthesis and basal metabolic processes.
gene ontology terms overrepresented in genes primarily expressed in embryo were mostly associated with cell division and protein synthesis .
gene ontology associated with genes with highest expression in liver indicated functions of specific metabolic processes - reactions involving oxygen and energy related processes .
also genes primarily expressed in muscles were associated with go terms related to energy utilisation and especially the function of the mitochondria .
in genes with the highest expression levels in skin there was an overrepresentation of go terms related to cytoskeletal structures and cell proliferation .
of main interest in relation to mhc and immune function were genes with primary expression in spleen.
go terms associated with expression in this tissue include "leukocyte adhesion", "immune response", "cell surface receptor linked signal transduction" and "chemokine activity", but also several terms related to ribosomal activity .
lastly, there were a large number of go terms overrepresented in genes with maximal expression in testes, including for example "spermatogenesis" and "microtubule motor activity" .
validation of expression profiling: "housekeeping" genes
we specifically investigated expression patterns in four widely used housekeeping genes that have been shown to have similar levels of expression over a wide range of tissues and treatments in birds  <cit> .
two highly expressed genes, ubiquitin  and glyceraldehyde-3-phosphate dehydrogenase , were represented by  <dig>  and  <dig>  reads, respectively.
the tissue specificities of gene expression  for these were  <dig>  for ub and  <dig>  for gapdh  two genes with medium expression levels also had low gene expression variation between tissues.
ribosomal protein s <dig>  was found in  <dig> reads and had a τ-value of  <dig> , while 60s ribosomal protein l <dig>  was found in  <dig> reads with a τ-value of  <dig>  .
case study of expression profiling: mhc genes
we surveyed expression data for  <dig> mhc-related genes found in the zebra finch assembly and targeted bac sequencing .
we found evidence for expression of ten different mhc related genes in the zebra finch .
thus the coverage of these genes is comparable to the rest of the genome .
among these there was evidence for one expressed mhc class i loci but we did not find expression of any mhc class ii loci in the present dataset.
this is not to say that there are no expressed mhc class ii molecules in the zebra finch, but only that these genes are expressed at too low levels in the sampled tissues to be detected using our methodology.
the expression patterns of mhc genes were generally tissue specific , with the highest expression levels for most genes in spleen.
a detailed presentation of the expression for specific mhc genes can be found in a separate supplementary text .
one of the mhc genes, cd <dig>  presents a case of alternative splicing.
as is true in many other species, we found evidence for at least two differently spliced isoforms of this gene, represented by different contigs in our  <dig> sequence assembly .
the total number of reads in each tissue library after trimming  and the tissue specificity of gene expression  are also given .
discussion
we have used transcriptomic data from six different tissues, generated by 454-sequencing  <cit> , to investigate expression patterns of different zebra finch genes.
our results highlight, in a new evolutionary lineage, a number of trends in the evolution of gene expression profiles.
genes with a high degree of tissue specificity in expression levels also had high non-synonymous to synonymous rate of nucleotide substitutions , while genes with high overall expression levels had low dn/ds ratios.
thus genes with a more specialised function  seem to be evolving at a higher rate  than genes with a more general function .
these results recapitulate those of axelsson and co-workers  <cit>  who analysed chicken expression profiles in conjunction with sequence divergence data from chicken and zebra finch.
similar patterns of molecular evolution and expression specificity have also been found in mammals  <cit> .
one important consequence of this finding for future studies of gene expression is that genes under strong positive selection might be missed if rna from the appropriate tissues is not sequenced.
in other words, the genes that are likely to be relevant for explaining genetic variation in ecologically important processes such as host-parasite co-evolution or reproduction  <cit>  may be relatively less likely to be sequenced.
this is particularly relevant in the present study system because the vast majority of gene expression studies in passerine birds have focussed on a single tissue, the brain.
overall about 65% of the annotated zebra finch transcripts were covered by  <dig> sequencing in this study.
an analysis of genes in well characterised metabolic pathways and signalling cascades  <cit>  also corroborate this number.
there is also some indication that more than one splice variant  <cit>  was detected for some of the genes.
most of the contigs that did not match any of the annotated zebra finch transcripts still gave highly relevant hits against the zebra finch genome, suggesting that these represent novel genes that have yet to be annotated in the zebra finch genome.
a few contigs that did not match anywhere in the zebra finch genome could either be part of genetic regions that have not been sequenced in the present zebra finch genome assembly or may represent contamination from other organisms.
higher coverage transcriptome sequencing will be needed to complete the zebra finch transcriptome and to fully characterize splice-variants.
genes primarily expressed in embryo had low dn/ds ratios, while genes with the highest expression in testes showed high ratio.
low dn/ds ratios of embryonically expressed genes may represent stabilizing selection and high evolutionary constraint on core developmental and housekeeping genes  <cit> .
high dn/ds in testes-expressed and reproductive genes has also been observed in human versus chimpanzee comparisons  <cit> , in drosophila  <cit>  and in mice  <cit> .
such a pattern may be attributable to sexual selection acting on genes important for traits involved in reproduction.
high dn/ds values of genes expressed primarily in spleen is also concordant with previous studies showing high rate of evolution in genes involved in the immune system  <cit> .
several of the mhc genes investigated in this study had primary expression in spleen and high dn/ds ratios of these genes are often seen as an indication of balancing selection acting on them  <cit> .
in addition to performing genome wide analyses we also used the  <dig> transcriptome sequence data to investigate specific genes of interest.
in particular, special attention was given to genes of the major histocompatibillity complex .
we found evidence for expression of ten different genes associated with the major histocompatibility complex.
most of these were primarily expressed in spleen, although there were also high levels of expression in brain and liver for some.
expression in the spleen is hardly surprising given the function of spleen in the immune defence.
the expression of mhc genes in the brain however, was only relatively recently discovered in mammals  <cit>  and has not been previously described in birds.
it will be of interest to determine whether the role of the mhc in the brain is conserved across vertebrates.
furthermore, several gene ontology terms related to immune response were overrepresented in genes with primary expression in spleen.
some go terms related to immune response were also overrepresented in genes with high tissue specificity, indicating that many immune genes are expressed mainly in a few specialised tissues.
for a few known mhc genes we could not detect any expression.
this illustrates the fact that one may not necessarily find specific genes of interest in a next generation transcriptome sequencing dataset, especially if they are expressed at very low levels or only in specific tissues or life history stages.
on the other hand, ongoing development of next generation sequencing technologies means that deeper coverage will be obtained enabling gene finding of lowly expressed genes.
coverage of mhc genes was within the range of other well characterised groups of genes related to specific metabolic and signalling pathways.
these genes had medium levels of tissue specificity of expression, and there was a tendency for mhc genes to have higher levels of expression specificity .
in expression profiling it is preferable to use sequences from a non-normalized cdna library to avoid bias in the estimates of expression individual genes  <cit> .
in our case the only data available for gene expression in different tissues came from cdna libraries that were normalized to increase the abundance of rare transcripts  <cit> .
thus there is a risk that our expression estimates might be biased.
in particular the expression levels of rare transcripts are probably overestimated while the levels for very common transcripts should be underestimated.
this also means that estimates of tissue specificity of gene expression  may be underestimated for individual genes.
we argue, however, that the comparative analyses presented here can be performed using this dataset.
there are at least four lines of evidence that these analyses are valid.
1) there is still considerable variation in expression levels between the different genes and tissues in our study, with many genes only expressed in one or a few tissues.
2) the analysis concerning gene expression gave results in the predicted direction.
for example there was a positive relationship between specificity of gene expression and dn/ds ratio  <cit> .
3) the expression of most mhc genes was by far strongest in spleen which is what would be predicted for genes involved in immune defence.
further, go terms overrepresented for genes with maximal expression in a certain tissue seemed to correspond well to those expected given the biological functions of the different tissues.
4) the expression levels of several housekeeping genes seemed to be stable across the different tissues analysed here.
one potential explanation for our failure to find  <dig> of the  <dig> mhc genes surveyed is that the relatively short contigs generated here, in combination with oligo dt priming, produced a strong 3' bias in the  <dig> sequencing.
indeed, many of the  <dig> reads fell in the 3' untranslated region  of genes .
it is therefore possible that these mhc genes were expressed, but the sequence reads only included utr sequence.
to investigate this issue we collected information from the avian mhc genes where the 3' utr has been sequenced.
utr regions of avian mhc genes are not well-described at this point but we found 3' utr sequence data for mhc class iib from chicken, turkey, quail, new zealand robin , bengalese finch  and zebra finch .
for mhc class i we found data from chicken, turkey, quail, mallard duck  and great reed warbler , and we also included data from duck cd <dig> .
these sequences were blasted against all zebra finch  <dig> contigs and positive matches were verified by a reciprocal blast against the zebra finch genome and chicken transcriptome databases.
only two of our contigs matched the 3' utr mhc sequences, both representing the cd <dig>  transcript.
therefore it is unlikely that the failure to detect more mhc genes can be attributed solely to the short, and 3' utr biased, contigs we assembled.
new and improved methods for library preparation are now used to deal with this problem of 3' bias.
in general, the contigs produced using de-novo assembly of the 454-reads only partially covered the gene transcripts, with a mean contig length of only  <dig> nucleotides.
these data were produced using the first generation of the 454-sequencing system  for which maximal read lengths were only around  <dig> bp.
with application of the new generation of 454-sequencing , which generates more and longer reads, one would expect to get longer contigs and more contigs covering the whole of the gene coding sequence  <cit> .
on the other hand deeper coverage of the transcriptome, and expression data on more genes, would be obtained using illumina/solexa or abi sollid technology.
both of these approaches generates a much larger amount of reads compared to  <dig> sequencing but at a cost of much shorter reads.
they are thus particularly useful for species, like the zebra finch, that have a characterized genome sequence.
this study highlights the utility of next-generation sequencing data for expression pattern profiling.
the zebra finch genome sequence was recently released and this, together with the gene predictions available, has been very useful when analysing the data.
still, this methodology would also work well when addressing a non-model species without any prior genome information  <cit> .
in particular, the long read lengths of the new titanium 454-generation means that many expressed genes can be identified using comparative sequence analysis against genomes of distantly related species.
the combination of data on sequence and gene expression variation makes this strategy useful for future studies in novel species.
however, our study also shows that it may not always be possible to find and sequence specific genes of interest using whole-transcriptome sequencing.
for example, we did not find any mhc class ii, tap or tapasin sequences, even though there is no reason to believe that these are not present and expressed in the zebra finch genome.
it may be that gene capture methods  <cit>  or more efficient cdna normalization and random primed libraries are needed to be able to pick up specific and very rare transcripts.
another approach to improving the discovery of genes specifically involved in the immune system would be to boost an immune response prior to cdna sampling.
CONCLUSIONS
our analysis of the zebra finch transcriptome extends conserved patterns of gene expression profiles and molecular evolution to the avian lineage.
genes with low overall and tissue specific expression were shown to evolve at a higher rate than genes with high and unspecific expression levels.
such genes were also shown to be related to biological functions such as reproduction and immune response.
furthermore genes with primary expression in spleen were often related to the immune function .
our results highlight the usefulness of next-generation sequence data for investigating expression profiles in the genome as well as in specific candidate genes.
however, as illustrated by our survey of mhc genes, it is far from certain that all genes of interest will be present in a given transcriptome sequencing run.
therefore care must thus be taken to ensure sampling of the appropriate tissues and life stages if the aim of the sequencing run is to examine specific gene families or physiological pathways.

 studies of molecular evolution have until recently focused on nucleotide divergence, while studies of variation in gene expression profiles have mainly been restricted to a few model species such as drosophila and mice. this is because the technologies for studying gene expression have not been available for non-model species. however, sequencing-based technologies for expression profiling can now be utilised to this end. by counting the number of reads generated by sequencing of cdna from different genes in the transcriptome, one can get an estimate of the expression level of these genes in the particular tissues sampled. a complementary approach is to scan publicly available databases of expressed sequence tags for the genes of interest. in addition to microarrays, these strategies, called digital transcriptomics, are today the most commonly used methods for investigating expression patterns. digital transcriptomics has received a great deal of attention, but the use of these methods has been restricted in many species by the requirement of having a reference genome to evaluate and analyse the data. the advent of massively parallel sequencing is now starting to change this picture by providing a cost-effective way of generating large amount of sequence data in species where there is no prior knowledge of the genome sequence. next-generation sequencing technology generally generates millions of short sequence reads, each read being tens to hundreds of base pairs long, depending on the specific platform. this enables detection of genes even with very low expression levels. roche 454-sequencing, in particular, generates reads that are long enough to be informative in the absence of a reference genome. here, we evaluate the use of 454-sequencing to investigate tissue specific gene expression profiles. next-generation sequencing can be used to not only describe genome-wide patterns of gene expression, but also to characterise specific gene families or genetic pathways. to illustrate this point, we use the ecologically important and widely studied genes of the major histocompatibility complex for a more detailed analysis. these genes are a very common focus of studies that take a candidate gene approach to investigate functionally important genetic variation in immune function. mhc genes are among the most variable of the vertebrate genomes. in particular, the classical mhc genes exhibit an extraordinary level of polymorphism. this polymorphism is strongly associated to the role of these genes in regulating and triggering the adaptive immune response. studies have found links between nucleic acid variation in the mhc genes and resistance to parasites, sexually selected ornaments, mate choice, maternal-foetal incompatibilities and local adaptation. typically studies of mhc variation have focused on sequence variation only in a few highly polymorphic regions of class i and class ii genes, while variation in other genes, regions and expression levels has largely been ignored. the completion of the genome sequence of first the chicken and now of the zebra finch have opened the door for in-depth studies of organisation and expression of mhc genes in birds. there are striking differences in the way the adaptive immune defence operates in birds compared to mammals and it could be envisioned that such studies will reveal new insights in the evolution of vertebrate immunity. the aim of the present study was to investigate tissue-specific gene expression patterns in the zebra finch. with the sequencing of its genome, the zebra finch has taken a major step towards becoming an important model system for bird genomics. outside of some recent studies of gene expression in brain, however, little is known about genome-scale, and organism-wide patterns of gene expression in song birds. in this study we describe patterns of gene expression across six zebra finch tissues and explore the relationship between expression profiles of genes and characteristics of their molecular evolution. to this end, we use a next-generation sequencing digital transcriptomics approach known as rna-seq. this methodology was recently employed to study gene expression differentiation between two subspecies of crow, but as far as we are aware, this is the first time that a bird transcriptome has been characterised in multiple tissues using an ngs rna-seq approach. in addition to global patterns of gene expression, we highlight patterns of expression in the genes of the mhc. because of the complex history of duplication among certain mhc genes, gene expression profiles have the potential to offer insight into the evolutionary fates of these duplicated genes. importantly, characterizing the expression of mhc genes will also facilitate downstream studies of these genes in ecological contexts by identifying functionally important loci.  assembly of sequencing reads after trimming and removal of contaminant sequences a total of, reads were available, with a mean read length of nucleotides. of these reads were assembled de-novo into contigs with a mean contig length of nucleotides and a mean of reads per contig. the total length of all contigs was kb. for read and contig statistics for each tissue separately see additional file 1: appendix s of the contigs showed signatures of multiple splice variants, as indicated by gaps in alignments between the contig and one or more of the reads that contribute to that contig. expression levels for contigs expression levels were highly variable between tissues and contigs. a vast majority of contigs were made up by only a few reads but some had indications of very high expression levels. as would be expected there was a positive relationship between contig length and contig depth. the tissue specificity of the expression of contigs was negatively correlated with the overall amount of expression. this negative correlation could be the result of a sampling artefact during the calculation of τ. we found however, that the observed correlation was significantly stronger than the mean simulated correlation based on unbiased τ-values calculated from randomisations of re-sampled data. this suggests that there is indeed a sampling bias in the calculation of τ but that it is not strong enough to alone explain our observed correlation. interestingly, τ seems to have a somewhat bimodal distribution with peaks around and, indicating that the genes may group into two different classes of tissue specificity. however, only a few contigs showed evidence of very low τ-values indicative of housekeeping genes. outlier contigs with high expression levels three contigs were found to have strikingly high overall expression levels; following blast searches against the chicken and zebra finch gene databases they were found to represent albumin, heat shock kda protein beta and nadh dehydrogenase subunit. all of these are considered to be so-called housekeeping genes and are also highly expressed in mammals. some contigs were conspicuous in having very strong expression in one or few tissues. genes represented by these include elongation factor 1-alpha, cytoplasmic beta-actin gene, haemoglobin alpha and mhc class ii associated invariant chain ii. one contig was found in high levels in testes but was almost completely absent in other tissues, the blast search revealed that this originates from a contamination with dna from a freshwater planarian. this is likely to have occurred in the laboratory that carried out the sequencing, since the genome of schmidtea mediterranea was being sequenced there at the same time as the zebra finch cdna preparation. this contig, together with other contigs resulting from contamination, was removed from the data before conducting downstream analyses. coverage of the zebra finch transcriptome contigs from the de-novo assembly and of the non-assembled singletons gave significant blast hits against at least one predicted zebra finch gene. since the contigs were generally much shorter than the total cdna length of the gene it was commonly found that several different contigs matched the same gene. in total zebra finch transcripts present in the biomart database were found to correspond to the 454/est transcriptome contigs and singletons. this represents 65% of the total characterised zebra finch transctiptome. the transcripts are derived from different genes, suggesting that more than one splice variant was detected for ~2% of the genes. we also identified potential novel splice variants for of the expressed zebra finch genes, as indicated by gaps in the alignments of the contig and the gene prediction. on average 38% of the lengths of represented transcripts were covered by contig sequences and transcripts were fully covered. to further investigate the extent of transcriptomic coverage, we investigated the presence of known genes in various metabolic pathways and signalling cascades. for the metabolic pathways about 85% of the genes were represented and for signalling cascades we found around 60%. of all genes found were expressed in all investigated tissues and were expressed exclusively in one tissue. out of the contigs that did not give any matches to known predicted zebra finch transcripts, most still gave highly significant blast hits against the zebra finch genome sequence, suggesting that these represent transcribed regions that have not yet been annotated. the remaining contigs may represent genes in regions of the genome that have not been sequenced and/or assembled in the current genome assembly, or additional contamination from other organisms that are not represented in genbank. number of mhc related genes included in this study is also given. the mean index of tissue specificity of expression and its 95% ci for each pathway is also given. nmax represents the number of genes with maximal expression in each of the tissues and the number within brackets is the number of genes which are expressed uniquely in that tissue. analyses of gene expression profiles we found a positive correlation between tissue specificity of expression and the ratio of non-synonymous to synonymous substitution rate when compared to the chicken orthologue of the gene in question. there was also a negative correlation between total expression level of the gene and ω. there was a weak positive correlation between the length of the gene and the total level of gene expression, and a negative correlation between gene length and τ. there were differences in ω between the tissues in which genes were primarily expressed. genes that were primarily expressed in the embryo had the lowest mean ω-value. the expression specificity of genes also varied significantly between tissues of maximal expression. the lowest tissue specificity was found in genes with primary expression in embryo and muscle, while the highest τ was found in genes with maximal expression in skin and testes. analysis of expression in relation to go-terms there were gene ontology terms overrepresented in genes with high levels of expression specificity. these represent processes such as cellular and organelle movement and specific enzymatic processes. some go terms overrepresented in genes with high tissue specificity are associated with reproduction and immune defence. genes identified as having low tissue specificity of gene expression were significantly overrepresented for different go terms. these terms generally represented functions such as protein synthesis and basal metabolic processes. gene ontology terms overrepresented in genes primarily expressed in embryo were mostly associated with cell division and protein synthesis. gene ontology associated with genes with highest expression in liver indicated functions of specific metabolic processes - reactions involving oxygen and energy related processes. also genes primarily expressed in muscles were associated with go terms related to energy utilisation and especially the function of the mitochondria. in genes with the highest expression levels in skin there was an overrepresentation of go terms related to cytoskeletal structures and cell proliferation. of main interest in relation to mhc and immune function were genes with primary expression in spleen. go terms associated with expression in this tissue include "leukocyte adhesion", "immune response", "cell surface receptor linked signal transduction" and "chemokine activity", but also several terms related to ribosomal activity. lastly, there were a large number of go terms overrepresented in genes with maximal expression in testes, including for example "spermatogenesis" and "microtubule motor activity". validation of expression profiling: "housekeeping" genes we specifically investigated expression patterns in four widely used housekeeping genes that have been shown to have similar levels of expression over a wide range of tissues and treatments in birds. two highly expressed genes, ubiquitin and glyceraldehyde-3-phosphate dehydrogenase, were represented by and reads, respectively. the tissue specificities of gene expression for these were for ub and for gapdh two genes with medium expression levels also had low gene expression variation between tissues. ribosomal protein s was found in reads and had a τ-value of, while 60s ribosomal protein l was found in reads with a τ-value of. case study of expression profiling: mhc genes we surveyed expression data for mhc-related genes found in the zebra finch assembly and targeted bac sequencing. we found evidence for expression of ten different mhc related genes in the zebra finch. thus the coverage of these genes is comparable to the rest of the genome. among these there was evidence for one expressed mhc class i loci but we did not find expression of any mhc class ii loci in the present dataset. this is not to say that there are no expressed mhc class ii molecules in the zebra finch, but only that these genes are expressed at too low levels in the sampled tissues to be detected using our methodology. the expression patterns of mhc genes were generally tissue specific, with the highest expression levels for most genes in spleen. a detailed presentation of the expression for specific mhc genes can be found in a separate supplementary text. one of the mhc genes, cd presents a case of alternative splicing. as is true in many other species, we found evidence for at least two differently spliced isoforms of this gene, represented by different contigs in our sequence assembly. the total number of reads in each tissue library after trimming and the tissue specificity of gene expression are also given. discussion we have used transcriptomic data from six different tissues, generated by 454-sequencing, to investigate expression patterns of different zebra finch genes. our results highlight, in a new evolutionary lineage, a number of trends in the evolution of gene expression profiles. genes with a high degree of tissue specificity in expression levels also had high non-synonymous to synonymous rate of nucleotide substitutions, while genes with high overall expression levels had low dn/ds ratios. thus genes with a more specialised function seem to be evolving at a higher rate than genes with a more general function. these results recapitulate those of axelsson and co-workers who analysed chicken expression profiles in conjunction with sequence divergence data from chicken and zebra finch. similar patterns of molecular evolution and expression specificity have also been found in mammals. one important consequence of this finding for future studies of gene expression is that genes under strong positive selection might be missed if rna from the appropriate tissues is not sequenced. in other words, the genes that are likely to be relevant for explaining genetic variation in ecologically important processes such as host-parasite co-evolution or reproduction may be relatively less likely to be sequenced. this is particularly relevant in the present study system because the vast majority of gene expression studies in passerine birds have focussed on a single tissue, the brain. overall about 65% of the annotated zebra finch transcripts were covered by sequencing in this study. an analysis of genes in well characterised metabolic pathways and signalling cascades also corroborate this number. there is also some indication that more than one splice variant was detected for some of the genes. most of the contigs that did not match any of the annotated zebra finch transcripts still gave highly relevant hits against the zebra finch genome, suggesting that these represent novel genes that have yet to be annotated in the zebra finch genome. a few contigs that did not match anywhere in the zebra finch genome could either be part of genetic regions that have not been sequenced in the present zebra finch genome assembly or may represent contamination from other organisms. higher coverage transcriptome sequencing will be needed to complete the zebra finch transcriptome and to fully characterize splice-variants. genes primarily expressed in embryo had low dn/ds ratios, while genes with the highest expression in testes showed high ratio. low dn/ds ratios of embryonically expressed genes may represent stabilizing selection and high evolutionary constraint on core developmental and housekeeping genes. high dn/ds in testes-expressed and reproductive genes has also been observed in human versus chimpanzee comparisons, in drosophila and in mice. such a pattern may be attributable to sexual selection acting on genes important for traits involved in reproduction. high dn/ds values of genes expressed primarily in spleen is also concordant with previous studies showing high rate of evolution in genes involved in the immune system. several of the mhc genes investigated in this study had primary expression in spleen and high dn/ds ratios of these genes are often seen as an indication of balancing selection acting on them. in addition to performing genome wide analyses we also used the transcriptome sequence data to investigate specific genes of interest. in particular, special attention was given to genes of the major histocompatibillity complex. we found evidence for expression of ten different genes associated with the major histocompatibility complex. most of these were primarily expressed in spleen, although there were also high levels of expression in brain and liver for some. expression in the spleen is hardly surprising given the function of spleen in the immune defence. the expression of mhc genes in the brain however, was only relatively recently discovered in mammals and has not been previously described in birds. it will be of interest to determine whether the role of the mhc in the brain is conserved across vertebrates. furthermore, several gene ontology terms related to immune response were overrepresented in genes with primary expression in spleen. some go terms related to immune response were also overrepresented in genes with high tissue specificity, indicating that many immune genes are expressed mainly in a few specialised tissues. for a few known mhc genes we could not detect any expression. this illustrates the fact that one may not necessarily find specific genes of interest in a next generation transcriptome sequencing dataset, especially if they are expressed at very low levels or only in specific tissues or life history stages. on the other hand, ongoing development of next generation sequencing technologies means that deeper coverage will be obtained enabling gene finding of lowly expressed genes. coverage of mhc genes was within the range of other well characterised groups of genes related to specific metabolic and signalling pathways. these genes had medium levels of tissue specificity of expression, and there was a tendency for mhc genes to have higher levels of expression specificity. in expression profiling it is preferable to use sequences from a non-normalized cdna library to avoid bias in the estimates of expression individual genes. in our case the only data available for gene expression in different tissues came from cdna libraries that were normalized to increase the abundance of rare transcripts. thus there is a risk that our expression estimates might be biased. in particular the expression levels of rare transcripts are probably overestimated while the levels for very common transcripts should be underestimated. this also means that estimates of tissue specificity of gene expression may be underestimated for individual genes. we argue, however, that the comparative analyses presented here can be performed using this dataset. there are at least four lines of evidence that these analyses are valid. 1) there is still considerable variation in expression levels between the different genes and tissues in our study, with many genes only expressed in one or a few tissues. 2) the analysis concerning gene expression gave results in the predicted direction. for example there was a positive relationship between specificity of gene expression and dn/ds ratio. 3) the expression of most mhc genes was by far strongest in spleen which is what would be predicted for genes involved in immune defence. further, go terms overrepresented for genes with maximal expression in a certain tissue seemed to correspond well to those expected given the biological functions of the different tissues. 4) the expression levels of several housekeeping genes seemed to be stable across the different tissues analysed here. one potential explanation for our failure to find of the mhc genes surveyed is that the relatively short contigs generated here, in combination with oligo dt priming, produced a strong 3' bias in the sequencing. indeed, many of the reads fell in the 3' untranslated region of genes. it is therefore possible that these mhc genes were expressed, but the sequence reads only included utr sequence. to investigate this issue we collected information from the avian mhc genes where the 3' utr has been sequenced. utr regions of avian mhc genes are not well-described at this point but we found 3' utr sequence data for mhc class iib from chicken, turkey, quail, new zealand robin, bengalese finch and zebra finch. for mhc class i we found data from chicken, turkey, quail, mallard duck and great reed warbler, and we also included data from duck cd. these sequences were blasted against all zebra finch contigs and positive matches were verified by a reciprocal blast against the zebra finch genome and chicken transcriptome databases. only two of our contigs matched the 3' utr mhc sequences, both representing the cd transcript. therefore it is unlikely that the failure to detect more mhc genes can be attributed solely to the short, and 3' utr biased, contigs we assembled. new and improved methods for library preparation are now used to deal with this problem of 3' bias. in general, the contigs produced using de-novo assembly of the 454-reads only partially covered the gene transcripts, with a mean contig length of only nucleotides. these data were produced using the first generation of the 454-sequencing system for which maximal read lengths were only around bp. with application of the new generation of 454-sequencing, which generates more and longer reads, one would expect to get longer contigs and more contigs covering the whole of the gene coding sequence. on the other hand deeper coverage of the transcriptome, and expression data on more genes, would be obtained using illumina/solexa or abi sollid technology. both of these approaches generates a much larger amount of reads compared to sequencing but at a cost of much shorter reads. they are thus particularly useful for species, like the zebra finch, that have a characterized genome sequence. this study highlights the utility of next-generation sequencing data for expression pattern profiling. the zebra finch genome sequence was recently released and this, together with the gene predictions available, has been very useful when analysing the data. still, this methodology would also work well when addressing a non-model species without any prior genome information. in particular, the long read lengths of the new titanium 454-generation means that many expressed genes can be identified using comparative sequence analysis against genomes of distantly related species. the combination of data on sequence and gene expression variation makes this strategy useful for future studies in novel species. however, our study also shows that it may not always be possible to find and sequence specific genes of interest using whole-transcriptome sequencing. for example, we did not find any mhc class ii, tap or tapasin sequences, even though there is no reason to believe that these are not present and expressed in the zebra finch genome. it may be that gene capture methods or more efficient cdna normalization and random primed libraries are needed to be able to pick up specific and very rare transcripts. another approach to improving the discovery of genes specifically involved in the immune system would be to boost an immune response prior to cdna sampling.  our analysis of the zebra finch transcriptome extends conserved patterns of gene expression profiles and molecular evolution to the avian lineage. genes with low overall and tissue specific expression were shown to evolve at a higher rate than genes with high and unspecific expression levels. such genes were also shown to be related to biological functions such as reproduction and immune response. furthermore genes with primary expression in spleen were often related to the immune function. our results highlight the usefulness of next-generation sequence data for investigating expression profiles in the genome as well as in specific candidate genes. however, as illustrated by our survey of mhc genes, it is far from certain that all genes of interest will be present in a given transcriptome sequencing run. therefore care must thus be taken to ensure sampling of the appropriate tissues and life stages if the aim of the sequencing run is to examine specific gene families or physiological pathways. 
here we investigate tissue specific gene expression patterns in the zebra finch with special emphasis on the genes of the major histocompatibility complex. expression of mhc genes in particular, corresponds well with expression patterns in other vertebrates. several mhc genes, including mhc class i also showed expression in brain. overall, the data show predicted patterns of gene expression profiles and molecular evolution in the zebra finch genome. furthermore, for all genes with highest levels of expression in spleen there was an overrepresentation of several gene ontology terms related to immune function. a total of zebra finch transcripts were represented in this est data, indicating a transcriptome coverage of about 65%. we found evidence for expression of different genes related to the mhc. in line with this, there was also a negative correlation between overall expression levels and expression specificity of contigs. mhc genes showed relatively tissue specific expression levels and were in general primarily expressed in spleen. almost million 454-sequencing reads from cdna of six different tissues were assembled and analysed. our study highlights the usefulness of next-generation sequence data for quantifying gene expression in the genome as a whole as well as in specific candidate genes. in order to understand patterns of adaptation and molecular evolution it is important to quantify both variation in gene expression and nucleotide sequence divergence. gene expression profiling in non-model organisms has recently been facilitated by the advent of massively parallel sequencing technology. there was a positive correlation between the tissue specificity of gene expression and non-synonymous to synonymous nucleotide substitution ratio of genes, suggesting that genes with a specialised function are evolving at a higher rate than genes with a more general function.
 in order to understand patterns of adaptation and molecular evolution it is important to quantify both variation in gene expression and nucleotide sequence divergence. gene expression profiling in non-model organisms has recently been facilitated by the advent of massively parallel sequencing technology. here we investigate tissue specific gene expression patterns in the zebra finch with special emphasis on the genes of the major histocompatibility complex.  almost million 454-sequencing reads from cdna of six different tissues were assembled and analysed. a total of zebra finch transcripts were represented in this est data, indicating a transcriptome coverage of about 65%. there was a positive correlation between the tissue specificity of gene expression and non-synonymous to synonymous nucleotide substitution ratio of genes, suggesting that genes with a specialised function are evolving at a higher rate than genes with a more general function. in line with this, there was also a negative correlation between overall expression levels and expression specificity of contigs. we found evidence for expression of different genes related to the mhc. mhc genes showed relatively tissue specific expression levels and were in general primarily expressed in spleen. several mhc genes, including mhc class i also showed expression in brain. furthermore, for all genes with highest levels of expression in spleen there was an overrepresentation of several gene ontology terms related to immune function.  our study highlights the usefulness of next-generation sequence data for quantifying gene expression in the genome as a whole as well as in specific candidate genes. overall, the data show predicted patterns of gene expression profiles and molecular evolution in the zebra finch genome. expression of mhc genes in particular, corresponds well with expression patterns in other vertebrates. 