BACKGROUND
understanding the genetic basis of phenotypic variation is a major challenge for modern biology.
phenotypes result from interactions between genes and environment during development and from interactions between genes and natural selection over evolutionary time.
two barriers to building a complete understanding of phenotypic variation are that  single genes rarely act alone in building a phenotype, so genome-wide information is needed, and  those organisms for which we have the best ecological knowledge are not generally those for which we have the best genomic knowledge.
sequencing entire genomes of non-model organisms is still out of reach for most researchers but sequencing smaller subsets of the genome, like transcriptomes, provides an attractive alternative.
transcriptomes correspond to the transcribed dna of an organism and therefore represent functional genomic data.
de novo assembly and annotation is easier for transcribed genes than for complete genomes because new sequences can be compared to conserved protein sequences and transcribed genes contain fewer repetitive elements.
recent advances in dna sequencing technology have reduced the cost and time associated with gathering large amounts of sequence data.
for example a single run of  <dig> gs flx technology can generate nearly  <dig> mbp  <cit> .
development of transcriptomes in "ecological model species" can provide access to functional and evolutionary analyses previously restricted to genetic model organisms.
a well characterized transcriptome can help identify genes underlying phenotypic variation in several ways.
candidate genes that have been identified and characterized in model organisms can be identified in transcriptome databases and tested for signatures of selection in wild populations.
genome-wide scans of selection can identify genes involved in adaptations to specific environments  <cit> .
rna sequencing  can produce short sequences that can be aligned to a reference transcriptome, and used as an assay for genome-wide rna expression  <cit> .
creating a reference transcriptome can therefore be an invaluable tool for deciphering the genetic architecture of adaptive traits in species for which complete genome sequence is not available.
several decades of research on the trinidadian guppy  have established the species as a model system in evolutionary biology, ecology, and animal behaviour  <cit> .
extensive documentation of parallel evolution along a repeated environmental gradient is one reason that guppies are an important evolutionary model.
waterfalls, characteristic of streams in northern trinidad, separate guppies into populations that differ in predation and other ecological factors .
comparative, common garden experiments and transplant studies have revealed parallel evolution of life history traits , male colour patterns  and behaviour  in response to the ecological differences above and below barrier waterfalls.
this research has provided textbook examples of the operation of natural selection .
despite their importance in ecology and evolutionary biology, few genetic or genomic tools have been available for guppies until very recently.
these resources currently include approximately  <dig>  expressed genes with annotation  <cit>  and a snp-based genetic map useful for qtl mapping  <cit>  and population genetic analyses  <cit> .
our goal here was to add to these resources by describing a de novo assembly of the guppy transcriptome and to test the resulting database for its completeness and its utility for several genome-scale analyses.
to assess the completeness of the transcriptome database, we compared it with those for other fish species and annotated it by searching protein and gene ontology databases.
we also identified putative single nucleotide polymorphisms  and sex-specific transcripts.
identifying sex-specific or sex-biased genes can shed light on sex determination, sexual dimorphism, and sex-specific selection  <cit> .
in addition, we identified candidate genes of particular interest to the guppy research community.
guppies are important models for studies of mate choice and other social behaviours, sexual selection and sexual dimorphism of colour patterns, and the effects of parasitism on fitness .
we therefore report transcripts that are candidates for genes implicated in behaviour, colour vision, skin colouration, and parasite resistance.
in addition to sequencing and annotating the transcriptome, we used illumina short-read sequencing to show that the assembled transcriptome can be a reliable reference for rna-seq analysis of differential gene expression, both in guppies and in a related species, the sailfin molly .
rna-seq does not depend on the hybridization chemistry of a microarray probe and therefore can detect expression variation over a large dynamic range  <cit> .
another advantage is that a priori knowledge about what genes will be expressed is not necessary, a particular asset in organisms for which a complete genome assembly is lacking.
despite clear benefits, efficient application of this technique to organisms without full genome assemblies has not been demonstrated.
to address this issue, we present data on gene expression response in juvenile guppies exposed to a natural predator, rivulus hartii.
anti-predator adaptation has been a major theme in guppy research for many years , and these data contribute to a larger research program intended to elucidate anti-predator behaviour and adaptation.
finally, we assessed the extent to which the guppy transcriptome could provide a reliable reference for gene expression studies in related species.
we mapped rna-seq data from a congeneric species  to the guppy reference transcriptome.
although microarray chips have been used for distantly related taxa with varying success  <cit>  the usefulness of a de novo reference transcriptome for studies of related taxa has not yet been determined.
RESULTS
transcriptome assembly
guppies used in the sequenced samples were from several natural populations encompassing the three major trinidadian river drainages; some were wild caught and the others reared in the lab under a range of conditions designed to maximize representation of expressed genes .
four separate normalized cdna libraries were created from pooled adult tissue: male brain, male body, female brain, and female body.
these libraries were sequenced on four separate plates, one library per plate, by  <dig> gs flx technology to produce  <dig> mbp of sequence.
a total of  <dig> , <dig> reads were assembled into contigs ;  <dig>  high quality reads were not assembled  and were excluded from further analysis.
the assembly produced  <dig>  contigs after excluding  <dig> contigs that were less than  <dig> bp long and less than  <dig> reads deep .
the average length of the remaining contigs was   <dig>  ±  <dig>  bp , and the n <dig> of the assembly was  <dig> bp.
the longest 10% of contigs were 892- <dig>  bp .
the average number of reads per contig was  <dig>  ±  <dig>  .
contigs in the top 10% of number of reads ranged from 72- <dig>  number of reads .
number of reads was significantly correlated with length of contig .
short-read files were deposited in the sequence read archive on genbank .
contigs and read files are also available from our website: http://www.bio.fsu.edu/kahughes/databases.html.
comparison to reference genomes
to our contig data, we added  <dig>  guppy expressed sequence tag  sequences available from genbank for a total of  <dig>  sequences.
we compared this full dataset to the medaka , three-spined stickleback , and zebrafish  unigene records  <cit> .
using blastn similarity searches , we found that  <dig>  sequences had matches in the medaka unigene database.
these sequences matched  <dig>  unique medaka records, covering  <dig> %  of the medaka transcriptome.
of these  <dig>  were matched by single guppy sequences and the remainder were matched by multiple guppy sequences;  <dig>  reciprocal best-hit matches were identified.
using tblastx similarity searches , we found that  <dig>  sequences had matches in the medaka unigene database, corresponding to  <dig>  unique medaka records, covering  <dig> %  of the medaka transcriptome.
of these  <dig>  matched single guppy reference sequences and the remainder matched multiple sequences.
similar results emerged when we searched the three-spined stickleback unigene database.
using blastn we found that  <dig>  sequences had matches in the stickleback database.
these sequences matched to  <dig>  unique stickleback records, covering  <dig> %  of the stickleback transcriptome.
of these  <dig>  matched to single guppy reference sequences and with the remainder matching multiple sequences;  <dig>  reciprocal best hit matches were found.
using tblastx we found that  <dig>  sequences had matches.
these sequences matched  <dig>  unique stickleback records, covering  <dig> %  of the transcriptome.
of these  <dig>  matched single guppy reference sequences, and with the remainder matched multiple sequences.
the zebrafish unigene records yielded more matches than did the other two databases, but they corresponded to a smaller percentage of the entire zebrafish transcriptome.
using blastn we found that  <dig>  sequences had matches in the zebrafish database.
these sequences matched to  <dig>  unique zebrafish records, covering  <dig> %  of the zebrafish transcriptome.
of these  <dig>  zebrafish sequences matched to single guppy reference sequences;  <dig>  reciprocal best hit matches were found.
similarly, using tblastx, we found that  <dig>  sequences had matches.
these sequences matched to  <dig>  unique zebrafish records, covering  <dig> %  of the transcriptome.
of these  <dig>  matched to single guppy reference sequences, and the remainder were matched by multiple reference sequences.
annotation
we annotated our database by first searching the swiss-prot  <cit>  and then the ncbi non-redundant  protein  <cit>  databases using blastx.
we found that  <dig>   sequences had matches, with  <dig>  unique records, in the swiss-prot database.
an additional  <dig>   sequences had matches in the nr database and matched to  <dig>  unique records.
in total,  <dig>   sequences were annotated and corresponded to  <dig>  unique swiss-prot or nr records.
if multiple guppy sequences matched the same record in either database, we grouped these sequences into "clusters" so that each cluster represented a unique match.
taxa with the most matches were human , mouse , and zebrafish , where percentages are based on the top hit for each annotated reference sequence.
in addition,  <dig> clusters matched records annotated as hypothetical proteins,  <dig> as uncharacterized proteins, and  <dig> as unknown proteins.
guppy sequences that had matches in either the swiss-prot or nr databases were annotated with gene ontology  annotations with the uniprot database  <cit> .
of these,  <dig>  of  <dig>   were annotated with go ids corresponding to  <dig>  unique matches in the uniprot database.
these unique matches were then grouped into generic go terms   <cit>  .
we found that  <dig>   records were annotated with a cellular component ,  <dig>   with a molecular function , and  <dig>   with a biological process .
representation of go categories in the guppy transcriptome set was similar to that in the zebrafish go gene association database with a few categories over- or underrepresented in each of the three main go categories .
after correcting for multiple tests, we found that  <dig> of the  <dig> comparisons were significantly over or underrepresented in comparison to the zebrafish records.
for example, in the biological-processes category, protein metabolic processes  and catabolic processes  were overexpressed in the guppy data, and multicellular organismal development  and embryonic development  were underrepresented.
we next investigated the utility of the guppy transcriptome data set for identifying candidate genes by searching for a subset of genes of particular interest to the guppy research community.
guppies have been used extensively in studies of mate choice, sexual selection, and parasite-mediated selection so we searched the transcriptome for candidate genes involved in visual communication, male ornaments, parasite resistance, and behavioural activation.
we found three clusters annotated as nonvisual opsins and eight as visual opsins, eight implicated in pigment synthesis ,  <dig> major histocompatibility complex genes , four clusters implicated in the regulation of behaviour or in behavioural activation, and two immediate early genes  .
presented are the gene description found in the protein database, the accession number, the number of contigs and est sequences in that annotation cluster, the database , the mean e-value, and mean percent coverage of the reference sequence by the contigs or est sequences.
male-specific expression
in the  <dig> assembly, we found  <dig> contigs that were assembled from male sequence reads only, which were therefore candidates for male-specific expression.
five of these contigs were annotated in the swiss-prot or nr database .
for one, the most significant match was a hypothetical protein and the second most significant match is therefore reported.
we chose these five annotated contigs plus one additional contig that had a high depth of coverage  to test for male specific expression using pcr amplification on whole body homogenate.
we confirmed that expression was male-specific for three of these contigs, whereas the remaining three were expressed in both males and females .
presented are the annotations with either the swiss-prot or nr databases, the number of reads used to assemble the contig, length of contig, the number of putative snps, and whether male-specific expression was confirmed with pcr.
snp discovery
we searched the contigs generated from  <dig> sequencing for single nucleotide polymorphisms using mosaik and gigabayes  <cit> .
we considered putative snps that had a bayesian probability of above  <dig>  and greater than 4× coverage.
we found  <dig>  putative snps meeting these criteria .
the mean coverage per snp was  <dig>  ±  <dig> .
a total of  <dig>  contigs had at least  <dig> snp.
the mean number of snps per contig was  <dig>  ±  <dig> .
four contigs had more than  <dig> snps, and  <dig> contigs had more than  <dig> snps.
we report the number of putative snps in male-specific and male-biased transcripts in table  <dig>
to test our prediction that gene classes thought to be under diversifying natural selection  would have higher numbers of snps, we compared the proportion of contigs without snps to the proportion of contigs with snps in generic go slim terms by means of a chi-squared test.
among the  <dig>  contigs possessing snps that were annotated with uniprot goids, we found no significant difference between the number of contigs with snps and the number without in any go slim term .
differential gene expression analysis with rna-seq
to investigate the utility of the guppy transcriptome for use in rna-seq based gene expression analysis, we conducted a small rna-seq experiment and mapped the resulting reads to a non-redundant version of our transcriptome data.
we compared two groups of juvenile fish, one reared in the presence of visual and chemical cues produced by rivulus hartii , and one reared identically but without predator cues.
each treatment group had two biological replicates, and each replicate consisted of mrna extracted from whole-head homogenate of a mean of two individuals per replicate.
we also sequenced one sample of whole-head homogenate taken from two sailfin mollies to determine the utility of the guppy transcriptome data for gene expression analysis in a related species.
the trinidadian guppy and sailfin molly are estimated to have diverged  <dig> mya  <cit> .
in total,  <dig> , <dig> high quality reads were obtained from for the guppy samples and  <dig> , <dig> high quality reads was obtained for the molly sample .
short read archive accession numbers is: study accession id: srp <dig>
to generate a non-redundant reference database for the rna-seq analysis, we performed a self blast search, using an e-value cut off  <dig> , of the entire sequence database .
if sequences were more than 90% identical and overlapped by more than 80% , they were grouped together.
the longest sequence of the group  was retained in the database, and the redundant sequences were removed.
the reduced reference database contained  <dig>  sequences .
we mapped rna-seq reads to the reduced database using bwa  <cit>  with default settings.
for each guppy sample, a mean of  <dig> , <dig> reads mapped to a unique sequence in the reference database with high confidence, corresponding to 52% of high quality reads.
over all four guppy samples,  <dig>  reference sequences  had reads align to them.
for the sailfin molly sample, > 40% of reads mapped to the guppy transcriptome, and reads mapped to nearly  <dig>  unique reference sequences .
to test for differences in counts between treatment groups in the predator-exposure experiment, we used generalized linear models  <cit>  and an empirical bayesian technique edger  <cit> .
both approaches model the distribution of count data as negative-binomial or poisson and account for differences between samples in the total number of reads .
we chose the negative binomial distribution in both analyses because our experiment included both biological and technical replication, and most genes were over-dispersed relative to the poisson expectation .
edger applies an empirical bayesian method to moderate dispersion estimates across reference sequences by borrowing information across all sequences in the analysis.
this moderation improves the reliability of inference in small- to moderate-sized experiments  <cit> .
to determine whether the data contained a signal of differential expression, we first evaluated the distribution of p-values obtained from generalized linear models applied to the counts for each transcript, as did bullard et al.
<cit> .
additional file  <dig> shows that this distribution is enriched for small p values, indicating that the data contain transcripts that are truly differently expressed in different treatment groups.
if expression were not truly different, we would expect the distribution of p values to be approximately uniform.
we used edger to identify reference sequences showing the strongest evidence for differential expression.
we found that the number of reference sequences classified as differently expressed depended strongly on the degree of moderation applied to the dispersion estimates.
with strong moderation ,  <dig> reference sequences were differently expressed at p <  <dig> ;  <dig> of these remained significant after correction for multiple testing by the method of benjamini and hochberg  <cit>  to control the false discovery rate .
with weak moderation , we found  <dig> reference sequences differently expressed at p <  <dig> ;  <dig> of these had fdr <  <dig>  after correction for multiple tests .
of the  <dig> reference sequences identified as differentially expressed with the more conservative model,  <dig> had higher counts in the predator-exposed treatment , and  <dig> had higher counts in the non-predator-exposed treatment .
annotations and counts for these  <dig> genes are given in table  <dig> and, figure  <dig>
shown is the log fold change  of the number of counts, p-value, corrected p-value , counts for samples not exposed to predators  and exposed to predators .
annotation is from either the swiss-prot or the nr databases with their accession number in brackets.
our ability to detect differential gene expression was not restricted to high-abundance transcripts.
in fact, figure  <dig> shows that differential expression was more likely to be detected at low to moderate abundance.
we also found little bias in favour of detecting differential expression for longer genes or reference sequences, as has been reported in some studies  <cit> .
the length of the reference sequence was very slightly but significantly correlated with the number of counts .
using the reduced reference database should have decreased the likelihood that a single read would map to multiple reference sequences.
nevertheless, we tested for bias that could result from eliminating any read with multiple hits by including all mapped reads in a reanalysis of the rna-seq data.
the only substantive change in the results was that one reference sequence  that did not meet the fdr <  <dig>  cut-off in the original analysis did meet that cut-off in the re-analysis.
discussion
using next generation sequencing technologies, we were able to sequence and annotate a reference transcriptome for the guppy.
approximately,  <dig> mbp were sequenced to create the reference transcriptome and then an additional  <dig> gbp were sequenced in the gene expression analysis.
our results represent the most extensive sequencing resource published for the guppy and indeed one of the largest sequencing projects for a non-model species .
by means of next generation sequence technologies, de novo transcriptome assemblies have been created for several such species in a variety of taxa: plants , insects , coral , a marine gastropod , fishes , and  <dig> species of birds  <cit> .
reference transcriptomes can be developed by assembly of the short sequence reads typical of next-generation sequencing technology into longer contiguous sequences more representative of a complete transcript.
the average length of contigs from our study  is comparable to those of other studies using similar technologies .
we obtained a greater depth of coverage per contig than did these other studies.
the average number of reads per contig was  <dig>  higher than that of other reference transcriptomes previously reported .
we were able to confirm the presence of  <dig>  contigs by pcr.
this validation and high degree of coverage suggest that our contigs are representative of real transcripts and not due to assembly error.
our reference transcriptome represents an extensive sampling of the guppy transcriptome.
by comparing it to those of species with well characterized genomes, we estimated that we have recovered roughly 40% of the entire transcriptome, corresponding to  <dig>  unique unigenes.
this number is probably an underestimate, however, because many guppy transcripts would probably not align to transcripts of species as divergent as the medaka, three-spined stickleback, or zebrafish .
we were able to match 37% of our database with functional annotations in the swiss-prot, nr protein, and go databases, a figure comparable to those from other studies using similar approaches .
while there are some differences between our reference database for the guppy and the available database for the zebrafish in go annotations, concordance in the overall distributions suggests that our library sampled widely across categories and provides a good representation of the transcriptome.
under-representation of a few go categories  is probably due to the choice of tissues in the two different sources of sequence data that we used.
the  <dig> data used to assemble the transcriptome were derived from adult tissue, whereas the ests downloaded from genbank represent many different developmental stages.
under-representation of genes involved in early development might be the result of the relatively larger amount of sequence derived from adults.
enriching the reference transcriptome with sequence derived from ovaries and developing embryos  could produce a more complete reference transcriptome
guppies have been the focus of evolutionary, ecological and behavioural research for over a century, and much of the motivation for developing a transcriptome resource was to provide tools to investigate genetic and genomic bases of adaptive variation.
we were able to identify a large number of genes in categories of particular interest to evolutionary and behavioural ecologists, supporting the utility of the transcriptome assembly in these fields.
we used a relatively liberal e-value cut-off, however, and, in some cases, the percent coverage of the subject sequence by the candidate gene sequence was relatively low, suggesting caution when these sequences are used in candidate gene studies.
coding-region snps are useful markers for mapping studies and for genome wide screens for signatures of selection, and we were able to identify putative snp markers in nearly  <dig>  unique contigs representing a broad sampling of functional categories.
use of sex-specific libraries for  <dig> sequencing was moderately successful in identifying genes with sex-specific expression.
we identified  <dig> candidate male-specific transcripts that were assembled from reads originating in male libraries; only 50% of transcripts tested were confirmed to be male-specific by a pcr assay.
a similar approach was taken by hale and colleagues  <cit> , who found that all  <dig> contigs unique to either sex's transcriptome assembly in sturgeon , were expressed in both sexes when tested by pcr.
in contrast, analysis of differential gene expression by mapping of short-read sequences to the reference transcriptome appears to be a robust strategy.
on average, 52% of reads mapped to the reference, and this coverage was sufficient to reveal differential expression for transcripts with as few as  <dig> mapped reads.
wolf and colleagues  <cit> , studying carrion crows , found that 69% of their short-sequence reads hit a zebra finch or chicken reference genome.
that study did not report using a mapping quality threshold, however, so the percentage of high-confidence mapped reads could be comparable to our own.
we studied gene expression differences in two treatment groups, one in which juvenile guppies were exposed to predators and one in which they were not.
at least two of the differently expressed reference sequences are good candidates for response to the predator exposure treatments: cerebellin- <dig> and cerebellin- <dig>  which were both expressed at higher levels in samples from predator exposed fish.
cerebellin- <dig> has been implicated in stress response: it stimulates the production of norepinephrine and increases andrenocortical secretion in rats  <cit> .
cerebellin- <dig> also has been implicated in regulation of neuronal processes  <cit> .
these results are encouraging in that we were able to obtain significant differences in transcript counts over a broad dynamic range even though the experiment was small.
results of this analysis were sensitive to parameters chosen for the empirical bayesian estimates of dispersion and somewhat sensitive to exclusion of reads that mapped to multiple reference sequences.
this sensitivity suggests that care should be taken in the choice of mapping and analysis parameters, especially in small experiments.
the rna-seq data also provided additional confirmation of the quality of the transcriptome assembly.
more than 70% of the transcriptome reference was matched by > <dig> reads in the rna-seq data.
because only a single age class of fish was used in this experiment, and only a small subset of tissues were sampled, this result suggests both that a large fraction of the assembled contigs in the data set are accurate, and that the rna-seq process recovered a high fraction of expressed genes.
in addition, the guppy transcriptome appears to be useful for rna-seq experiments in related species.
approximately 40% of reads from a sailfin molly sample mapped to the guppy data set.
CONCLUSIONS
like many organisms that are of great interest in ecological, evolutionary, and behavioural research, guppies and their close relatives lack complete genome sequences and most other genetic tools and resources.
here we show that next-generation sequencing provided a reliable, broad reference transcriptome that we have assembled and annotated.
this resource allowed us to identify candidate genes, putative snps, sex-specific gene expression, and differential gene expression at a genomic scale.
the reference transcriptome also proved useful for transcript mapping in a related species, the sailfin molly.

 understanding the genetic basis of phenotypic variation is a major challenge for modern biology. phenotypes result from interactions between genes and environment during development and from interactions between genes and natural selection over evolutionary time. two barriers to building a complete understanding of phenotypic variation are that single genes rarely act alone in building a phenotype, so genome-wide information is needed, and those organisms for which we have the best ecological knowledge are not generally those for which we have the best genomic knowledge. sequencing entire genomes of non-model organisms is still out of reach for most researchers but sequencing smaller subsets of the genome, like transcriptomes, provides an attractive alternative. transcriptomes correspond to the transcribed dna of an organism and therefore represent functional genomic data. de novo assembly and annotation is easier for transcribed genes than for complete genomes because new sequences can be compared to conserved protein sequences and transcribed genes contain fewer repetitive elements. recent advances in dna sequencing technology have reduced the cost and time associated with gathering large amounts of sequence data. for example a single run of gs flx technology can generate nearly mbp. development of transcriptomes in "ecological model species" can provide access to functional and evolutionary analyses previously restricted to genetic model organisms. a well characterized transcriptome can help identify genes underlying phenotypic variation in several ways. candidate genes that have been identified and characterized in model organisms can be identified in transcriptome databases and tested for signatures of selection in wild populations. genome-wide scans of selection can identify genes involved in adaptations to specific environments. rna sequencing can produce short sequences that can be aligned to a reference transcriptome, and used as an assay for genome-wide rna expression. creating a reference transcriptome can therefore be an invaluable tool for deciphering the genetic architecture of adaptive traits in species for which complete genome sequence is not available. several decades of research on the trinidadian guppy have established the species as a model system in evolutionary biology, ecology, and animal behaviour. extensive documentation of parallel evolution along a repeated environmental gradient is one reason that guppies are an important evolutionary model. waterfalls, characteristic of streams in northern trinidad, separate guppies into populations that differ in predation and other ecological factors. comparative, common garden experiments and transplant studies have revealed parallel evolution of life history traits, male colour patterns and behaviour in response to the ecological differences above and below barrier waterfalls. this research has provided textbook examples of the operation of natural selection. despite their importance in ecology and evolutionary biology, few genetic or genomic tools have been available for guppies until very recently. these resources currently include approximately expressed genes with annotation and a snp-based genetic map useful for qtl mapping and population genetic analyses. our goal here was to add to these resources by describing a de novo assembly of the guppy transcriptome and to test the resulting database for its completeness and its utility for several genome-scale analyses. to assess the completeness of the transcriptome database, we compared it with those for other fish species and annotated it by searching protein and gene ontology databases. we also identified putative single nucleotide polymorphisms and sex-specific transcripts. identifying sex-specific or sex-biased genes can shed light on sex determination, sexual dimorphism, and sex-specific selection. in addition, we identified candidate genes of particular interest to the guppy research community. guppies are important models for studies of mate choice and other social behaviours, sexual selection and sexual dimorphism of colour patterns, and the effects of parasitism on fitness. we therefore report transcripts that are candidates for genes implicated in behaviour, colour vision, skin colouration, and parasite resistance. in addition to sequencing and annotating the transcriptome, we used illumina short-read sequencing to show that the assembled transcriptome can be a reliable reference for rna-seq analysis of differential gene expression, both in guppies and in a related species, the sailfin molly. rna-seq does not depend on the hybridization chemistry of a microarray probe and therefore can detect expression variation over a large dynamic range. another advantage is that a priori knowledge about what genes will be expressed is not necessary, a particular asset in organisms for which a complete genome assembly is lacking. despite clear benefits, efficient application of this technique to organisms without full genome assemblies has not been demonstrated. to address this issue, we present data on gene expression response in juvenile guppies exposed to a natural predator, rivulus hartii. anti-predator adaptation has been a major theme in guppy research for many years, and these data contribute to a larger research program intended to elucidate anti-predator behaviour and adaptation. finally, we assessed the extent to which the guppy transcriptome could provide a reliable reference for gene expression studies in related species. we mapped rna-seq data from a congeneric species to the guppy reference transcriptome. although microarray chips have been used for distantly related taxa with varying success the usefulness of a de novo reference transcriptome for studies of related taxa has not yet been determined.  transcriptome assembly guppies used in the sequenced samples were from several natural populations encompassing the three major trinidadian river drainages; some were wild caught and the others reared in the lab under a range of conditions designed to maximize representation of expressed genes. four separate normalized cdna libraries were created from pooled adult tissue: male brain, male body, female brain, and female body. these libraries were sequenced on four separate plates, one library per plate, by gs flx technology to produce mbp of sequence. a total of, reads were assembled into contigs; high quality reads were not assembled and were excluded from further analysis. the assembly produced contigs after excluding contigs that were less than bp long and less than reads deep. the average length of the remaining contigs was ± bp, and the n of the assembly was bp. the longest 10% of contigs were 892- bp. the average number of reads per contig was ±. contigs in the top 10% of number of reads ranged from 72- number of reads. number of reads was significantly correlated with length of contig. short-read files were deposited in the sequence read archive on genbank. contigs and read files are also available from our website: http://www.bio.fsu.edu/kahughes/databases.html. comparison to reference genomes to our contig data, we added guppy expressed sequence tag sequences available from genbank for a total of sequences. we compared this full dataset to the medaka, three-spined stickleback, and zebrafish unigene records. using blastn similarity searches, we found that sequences had matches in the medaka unigene database. these sequences matched unique medaka records, covering % of the medaka transcriptome. of these were matched by single guppy sequences and the remainder were matched by multiple guppy sequences; reciprocal best-hit matches were identified. using tblastx similarity searches, we found that sequences had matches in the medaka unigene database, corresponding to unique medaka records, covering % of the medaka transcriptome. of these matched single guppy reference sequences and the remainder matched multiple sequences. similar results emerged when we searched the three-spined stickleback unigene database. using blastn we found that sequences had matches in the stickleback database. these sequences matched to unique stickleback records, covering % of the stickleback transcriptome. of these matched to single guppy reference sequences and with the remainder matching multiple sequences; reciprocal best hit matches were found. using tblastx we found that sequences had matches. these sequences matched unique stickleback records, covering % of the transcriptome. of these matched single guppy reference sequences, and with the remainder matched multiple sequences. the zebrafish unigene records yielded more matches than did the other two databases, but they corresponded to a smaller percentage of the entire zebrafish transcriptome. using blastn we found that sequences had matches in the zebrafish database. these sequences matched to unique zebrafish records, covering % of the zebrafish transcriptome. of these zebrafish sequences matched to single guppy reference sequences; reciprocal best hit matches were found. similarly, using tblastx, we found that sequences had matches. these sequences matched to unique zebrafish records, covering % of the transcriptome. of these matched to single guppy reference sequences, and the remainder were matched by multiple reference sequences. annotation we annotated our database by first searching the swiss-prot and then the ncbi non-redundant protein databases using blastx. we found that sequences had matches, with unique records, in the swiss-prot database. an additional sequences had matches in the nr database and matched to unique records. in total, sequences were annotated and corresponded to unique swiss-prot or nr records. if multiple guppy sequences matched the same record in either database, we grouped these sequences into "clusters" so that each cluster represented a unique match. taxa with the most matches were human, mouse, and zebrafish, where percentages are based on the top hit for each annotated reference sequence. in addition, clusters matched records annotated as hypothetical proteins, as uncharacterized proteins, and as unknown proteins. guppy sequences that had matches in either the swiss-prot or nr databases were annotated with gene ontology annotations with the uniprot database. of these, of were annotated with go ids corresponding to unique matches in the uniprot database. these unique matches were then grouped into generic go terms. we found that records were annotated with a cellular component, with a molecular function, and with a biological process. representation of go categories in the guppy transcriptome set was similar to that in the zebrafish go gene association database with a few categories over- or underrepresented in each of the three main go categories. after correcting for multiple tests, we found that of the comparisons were significantly over or underrepresented in comparison to the zebrafish records. for example, in the biological-processes category, protein metabolic processes and catabolic processes were overexpressed in the guppy data, and multicellular organismal development and embryonic development were underrepresented. we next investigated the utility of the guppy transcriptome data set for identifying candidate genes by searching for a subset of genes of particular interest to the guppy research community. guppies have been used extensively in studies of mate choice, sexual selection, and parasite-mediated selection so we searched the transcriptome for candidate genes involved in visual communication, male ornaments, parasite resistance, and behavioural activation. we found three clusters annotated as nonvisual opsins and eight as visual opsins, eight implicated in pigment synthesis, major histocompatibility complex genes, four clusters implicated in the regulation of behaviour or in behavioural activation, and two immediate early genes. presented are the gene description found in the protein database, the accession number, the number of contigs and est sequences in that annotation cluster, the database, the mean e-value, and mean percent coverage of the reference sequence by the contigs or est sequences. male-specific expression in the assembly, we found contigs that were assembled from male sequence reads only, which were therefore candidates for male-specific expression. five of these contigs were annotated in the swiss-prot or nr database. for one, the most significant match was a hypothetical protein and the second most significant match is therefore reported. we chose these five annotated contigs plus one additional contig that had a high depth of coverage to test for male specific expression using pcr amplification on whole body homogenate. we confirmed that expression was male-specific for three of these contigs, whereas the remaining three were expressed in both males and females. presented are the annotations with either the swiss-prot or nr databases, the number of reads used to assemble the contig, length of contig, the number of putative snps, and whether male-specific expression was confirmed with pcr. snp discovery we searched the contigs generated from sequencing for single nucleotide polymorphisms using mosaik and gigabayes. we considered putative snps that had a bayesian probability of above and greater than 4× coverage. we found putative snps meeting these criteria. the mean coverage per snp was ±. a total of contigs had at least snp. the mean number of snps per contig was ±. four contigs had more than snps, and contigs had more than snps. we report the number of putative snps in male-specific and male-biased transcripts in table to test our prediction that gene classes thought to be under diversifying natural selection would have higher numbers of snps, we compared the proportion of contigs without snps to the proportion of contigs with snps in generic go slim terms by means of a chi-squared test. among the contigs possessing snps that were annotated with uniprot goids, we found no significant difference between the number of contigs with snps and the number without in any go slim term. differential gene expression analysis with rna-seq to investigate the utility of the guppy transcriptome for use in rna-seq based gene expression analysis, we conducted a small rna-seq experiment and mapped the resulting reads to a non-redundant version of our transcriptome data. we compared two groups of juvenile fish, one reared in the presence of visual and chemical cues produced by rivulus hartii, and one reared identically but without predator cues. each treatment group had two biological replicates, and each replicate consisted of mrna extracted from whole-head homogenate of a mean of two individuals per replicate. we also sequenced one sample of whole-head homogenate taken from two sailfin mollies to determine the utility of the guppy transcriptome data for gene expression analysis in a related species. the trinidadian guppy and sailfin molly are estimated to have diverged mya. in total,, high quality reads were obtained from for the guppy samples and, high quality reads was obtained for the molly sample. short read archive accession numbers is: study accession id: srp to generate a non-redundant reference database for the rna-seq analysis, we performed a self blast search, using an e-value cut off, of the entire sequence database. if sequences were more than 90% identical and overlapped by more than 80%, they were grouped together. the longest sequence of the group was retained in the database, and the redundant sequences were removed. the reduced reference database contained sequences. we mapped rna-seq reads to the reduced database using bwa with default settings. for each guppy sample, a mean of, reads mapped to a unique sequence in the reference database with high confidence, corresponding to 52% of high quality reads. over all four guppy samples, reference sequences had reads align to them. for the sailfin molly sample, > 40% of reads mapped to the guppy transcriptome, and reads mapped to nearly unique reference sequences. to test for differences in counts between treatment groups in the predator-exposure experiment, we used generalized linear models and an empirical bayesian technique edger. both approaches model the distribution of count data as negative-binomial or poisson and account for differences between samples in the total number of reads. we chose the negative binomial distribution in both analyses because our experiment included both biological and technical replication, and most genes were over-dispersed relative to the poisson expectation. edger applies an empirical bayesian method to moderate dispersion estimates across reference sequences by borrowing information across all sequences in the analysis. this moderation improves the reliability of inference in small- to moderate-sized experiments. to determine whether the data contained a signal of differential expression, we first evaluated the distribution of p-values obtained from generalized linear models applied to the counts for each transcript, as did bullard et al.. additional file shows that this distribution is enriched for small p values, indicating that the data contain transcripts that are truly differently expressed in different treatment groups. if expression were not truly different, we would expect the distribution of p values to be approximately uniform. we used edger to identify reference sequences showing the strongest evidence for differential expression. we found that the number of reference sequences classified as differently expressed depended strongly on the degree of moderation applied to the dispersion estimates. with strong moderation, reference sequences were differently expressed at p <; of these remained significant after correction for multiple testing by the method of benjamini and hochberg to control the false discovery rate. with weak moderation, we found reference sequences differently expressed at p <; of these had fdr < after correction for multiple tests. of the reference sequences identified as differentially expressed with the more conservative model, had higher counts in the predator-exposed treatment, and had higher counts in the non-predator-exposed treatment. annotations and counts for these genes are given in table and, figure shown is the log fold change of the number of counts, p-value, corrected p-value, counts for samples not exposed to predators and exposed to predators. annotation is from either the swiss-prot or the nr databases with their accession number in brackets. our ability to detect differential gene expression was not restricted to high-abundance transcripts. in fact, figure shows that differential expression was more likely to be detected at low to moderate abundance. we also found little bias in favour of detecting differential expression for longer genes or reference sequences, as has been reported in some studies. the length of the reference sequence was very slightly but significantly correlated with the number of counts. using the reduced reference database should have decreased the likelihood that a single read would map to multiple reference sequences. nevertheless, we tested for bias that could result from eliminating any read with multiple hits by including all mapped reads in a reanalysis of the rna-seq data. the only substantive change in the results was that one reference sequence that did not meet the fdr < cut-off in the original analysis did meet that cut-off in the re-analysis. discussion using next generation sequencing technologies, we were able to sequence and annotate a reference transcriptome for the guppy. approximately, mbp were sequenced to create the reference transcriptome and then an additional gbp were sequenced in the gene expression analysis. our results represent the most extensive sequencing resource published for the guppy and indeed one of the largest sequencing projects for a non-model species. by means of next generation sequence technologies, de novo transcriptome assemblies have been created for several such species in a variety of taxa: plants, insects, coral, a marine gastropod, fishes, and species of birds. reference transcriptomes can be developed by assembly of the short sequence reads typical of next-generation sequencing technology into longer contiguous sequences more representative of a complete transcript. the average length of contigs from our study is comparable to those of other studies using similar technologies. we obtained a greater depth of coverage per contig than did these other studies. the average number of reads per contig was higher than that of other reference transcriptomes previously reported. we were able to confirm the presence of contigs by pcr. this validation and high degree of coverage suggest that our contigs are representative of real transcripts and not due to assembly error. our reference transcriptome represents an extensive sampling of the guppy transcriptome. by comparing it to those of species with well characterized genomes, we estimated that we have recovered roughly 40% of the entire transcriptome, corresponding to unique unigenes. this number is probably an underestimate, however, because many guppy transcripts would probably not align to transcripts of species as divergent as the medaka, three-spined stickleback, or zebrafish. we were able to match 37% of our database with functional annotations in the swiss-prot, nr protein, and go databases, a figure comparable to those from other studies using similar approaches. while there are some differences between our reference database for the guppy and the available database for the zebrafish in go annotations, concordance in the overall distributions suggests that our library sampled widely across categories and provides a good representation of the transcriptome. under-representation of a few go categories is probably due to the choice of tissues in the two different sources of sequence data that we used. the data used to assemble the transcriptome were derived from adult tissue, whereas the ests downloaded from genbank represent many different developmental stages. under-representation of genes involved in early development might be the result of the relatively larger amount of sequence derived from adults. enriching the reference transcriptome with sequence derived from ovaries and developing embryos could produce a more complete reference transcriptome guppies have been the focus of evolutionary, ecological and behavioural research for over a century, and much of the motivation for developing a transcriptome resource was to provide tools to investigate genetic and genomic bases of adaptive variation. we were able to identify a large number of genes in categories of particular interest to evolutionary and behavioural ecologists, supporting the utility of the transcriptome assembly in these fields. we used a relatively liberal e-value cut-off, however, and, in some cases, the percent coverage of the subject sequence by the candidate gene sequence was relatively low, suggesting caution when these sequences are used in candidate gene studies. coding-region snps are useful markers for mapping studies and for genome wide screens for signatures of selection, and we were able to identify putative snp markers in nearly unique contigs representing a broad sampling of functional categories. use of sex-specific libraries for sequencing was moderately successful in identifying genes with sex-specific expression. we identified candidate male-specific transcripts that were assembled from reads originating in male libraries; only 50% of transcripts tested were confirmed to be male-specific by a pcr assay. a similar approach was taken by hale and colleagues, who found that all contigs unique to either sex's transcriptome assembly in sturgeon, were expressed in both sexes when tested by pcr. in contrast, analysis of differential gene expression by mapping of short-read sequences to the reference transcriptome appears to be a robust strategy. on average, 52% of reads mapped to the reference, and this coverage was sufficient to reveal differential expression for transcripts with as few as mapped reads. wolf and colleagues, studying carrion crows, found that 69% of their short-sequence reads hit a zebra finch or chicken reference genome. that study did not report using a mapping quality threshold, however, so the percentage of high-confidence mapped reads could be comparable to our own. we studied gene expression differences in two treatment groups, one in which juvenile guppies were exposed to predators and one in which they were not. at least two of the differently expressed reference sequences are good candidates for response to the predator exposure treatments: cerebellin- and cerebellin- which were both expressed at higher levels in samples from predator exposed fish. cerebellin- has been implicated in stress response: it stimulates the production of norepinephrine and increases andrenocortical secretion in rats. cerebellin- also has been implicated in regulation of neuronal processes. these results are encouraging in that we were able to obtain significant differences in transcript counts over a broad dynamic range even though the experiment was small. results of this analysis were sensitive to parameters chosen for the empirical bayesian estimates of dispersion and somewhat sensitive to exclusion of reads that mapped to multiple reference sequences. this sensitivity suggests that care should be taken in the choice of mapping and analysis parameters, especially in small experiments. the rna-seq data also provided additional confirmation of the quality of the transcriptome assembly. more than 70% of the transcriptome reference was matched by > reads in the rna-seq data. because only a single age class of fish was used in this experiment, and only a small subset of tissues were sampled, this result suggests both that a large fraction of the assembled contigs in the data set are accurate, and that the rna-seq process recovered a high fraction of expressed genes. in addition, the guppy transcriptome appears to be useful for rna-seq experiments in related species. approximately 40% of reads from a sailfin molly sample mapped to the guppy data set.  like many organisms that are of great interest in ecological, evolutionary, and behavioural research, guppies and their close relatives lack complete genome sequences and most other genetic tools and resources. here we show that next-generation sequencing provided a reliable, broad reference transcriptome that we have assembled and annotated. this resource allowed us to identify candidate genes, putative snps, sex-specific gene expression, and differential gene expression at a genomic scale. the reference transcriptome also proved useful for transcript mapping in a related species, the sailfin molly. 
here we present a de novo assembly of the guppy transcriptome using sequence reads, and we evaluate potential uses of this transcriptome, including detection of sex-specific transcripts and deployment as a reference for gene expression analysis in guppies and a related species. we identified transcripts that, in juveniles, are regulated differently in the presence and absence of an important predator, rivulus hartii, including two genes implicated in stress response. in addition, we evaluated the use of the guppy reference transcriptome for gene expression analyses in a congeneric species, the sailfin molly. for each sample in the rna-seq study, >50% of high-quality reads mapped to unique sequences in the reference database with high confidence. guppies have been model organisms in ecology, evolutionary biology, and animal behaviour for over years. we generated approximately mbp of mrna sequence data from male brain, male body, female brain, and female body. we annotated nearly 40% of this reference transcriptome by searching protein and gene ontology databases. we show that next-generation sequencing provided a reliable and broad reference transcriptome. using this annotated transcriptome database, we identified candidate genes of interest to the guppy research community, putative single nucleotide polymorphisms, and male-specific expressed genes. this resource allowed us to identify candidate gene variants, snps in coding regions, and sex-specific gene expression, and permitted quantitative analysis of differential gene expression. the resulting, reads assembled into contigs, creating a reference transcriptome for the guppy with an average read depth of 28×. we also showed that our reference transcriptome can be used for rna-sequencing-based analysis of differential gene expression. an annotated transcriptome and other genomic tools will facilitate understanding the genetic and molecular bases of adaptation and variation in a vertebrate species with a uniquely well known natural history. over 40% of reads from the sailfin molly sample aligned to the guppy transcriptome.
 next-generation sequencing is providing researchers with a relatively fast and affordable option for developing genomic resources for organisms that are not among the traditional genetic models. here we present a de novo assembly of the guppy transcriptome using sequence reads, and we evaluate potential uses of this transcriptome, including detection of sex-specific transcripts and deployment as a reference for gene expression analysis in guppies and a related species. guppies have been model organisms in ecology, evolutionary biology, and animal behaviour for over years. an annotated transcriptome and other genomic tools will facilitate understanding the genetic and molecular bases of adaptation and variation in a vertebrate species with a uniquely well known natural history.  we generated approximately mbp of mrna sequence data from male brain, male body, female brain, and female body. the resulting, reads assembled into contigs, creating a reference transcriptome for the guppy with an average read depth of 28×. we annotated nearly 40% of this reference transcriptome by searching protein and gene ontology databases. using this annotated transcriptome database, we identified candidate genes of interest to the guppy research community, putative single nucleotide polymorphisms, and male-specific expressed genes. we also showed that our reference transcriptome can be used for rna-sequencing-based analysis of differential gene expression. we identified transcripts that, in juveniles, are regulated differently in the presence and absence of an important predator, rivulus hartii, including two genes implicated in stress response. for each sample in the rna-seq study, >50% of high-quality reads mapped to unique sequences in the reference database with high confidence. in addition, we evaluated the use of the guppy reference transcriptome for gene expression analyses in a congeneric species, the sailfin molly. over 40% of reads from the sailfin molly sample aligned to the guppy transcriptome.  we show that next-generation sequencing provided a reliable and broad reference transcriptome. this resource allowed us to identify candidate gene variants, snps in coding regions, and sex-specific gene expression, and permitted quantitative analysis of differential gene expression. 