BACKGROUND
among the most prevalent evolutionary transitions in plants is the shift from cross-fertilization to predominant self-fertilization among numerous angiosperm lineages  <cit> .
this change in mating system has important consequences for many aspects of the biology of selfing taxa including population genetic structure, colonizing ability, genome evolution and the morphology of flowers  <cit> .
the loss of floral mechanisms that reduce the incidence of self-fertilization results in high rates of autogamous selfing, leading to the evolution of the 'selfing syndrome'  <cit> .
although studies of the causes and consequences of cross- and self-fertilization in flowering plants have a long and venerable history, beginning with darwin's seminal work  <cit> , relatively little is known about the underlying molecular changes that accompany the transition from outcrossing to selfing  <cit>  and genomic analyses of related outcrossing and selfing plants are in their infancy.
recent technological advances in dna sequencing technology have removed several limitations associated with gathering large amounts of genomic data from non-model organisms  <cit> , providing opportunities for detailed investigation into the genomics of mating-system variation and evolution.
although assembling large eukaryotic genomes, de novo, may not yet be practical , sets of expressed genes or transcriptomes present a viable and attractive alternative to population genetic analyses of whole genome sequences.
transcriptomes represent a fraction of the total genome in size, contain fewer repetitive elements, and by selecting specific tissues they can be enriched for genes relevant to the particular aims of the research.
in addition, if the rna sample is not normalized the relative abundance of different reads has been shown to accurately reflect the expression level of transcripts in the tissue .
despite these potential advantages there remain a number of challenges for de novo transcriptome assembly, including gene duplication or paralogy, heterozygosity and alternative splicing, each of which require careful consideration.
there are relatively few studies to date involving de novo transcriptome assembly in non-model organisms.
so far the majority have used the roche  <dig> gs platform  which has the advantage of longer reads, but produces a fraction of the total amount of sequence produced per instrument per run compared with the illumina gaii platform  .
therefore, to maximize the coverage for rare transcripts, cdna samples are typically normalized.
as a result, studies using  <dig> are not able to estimate expression levels of different ests, as this requires deep sequencing of non-normalized cdna.
further, due to lower sequencing depth many transcripts are represented by a single read, and others by very few reads.
this can create problems in accurately distinguishing snps from errors, and in retrieving orthologous transcripts for sequence comparisons across experiments or species.
therefore, new methods are required to generate and assemble large datasets, many of which currently consist of substantially shorter reads.
here, we present de novo floral transcriptome assemblies using short read sequencing of four individual plants of neotropical eichhornia  species that vary in floral morphology and mating system.
the samples include three individuals of e. paniculata, including two from independently derived selfing populations and the third an outcrosser.
eichhornia paniculata is an annual diploid that has been the subject of detailed studies over the past two decades on the ecology and genetics of mating-system variation .
populations of e. paniculata are largely concentrated in northeastern brazil, with smaller foci in jamaica and cuba and isolated localities in nicaragua and mexico.
populations in brazil are largely outcrossing and possess the sexual polymorphism tristyly, which promotes cross-pollination among the three floral morphs .
morphological, genetic and biogeographical evidence indicates that tristyly in e. paniculata has broken down on multiple occasions resulting in independently derived selfing populations  <cit> .
populations in jamaica are largely composed of selfing variants of the mid-styled morph  in which short-level stamens are elongated to a position adjacent to mid-level stigmas resulting in autonomous self-fertilization.
in contrast, plants in mexico and nicaragua are selfing variants of the long-styled morph  with a different arrangement of their sexual organs .
although both variants possess the selfing syndrome, comparisons of molecular variation at  <dig> est-derived nuclear loci indicate a high level of differentiation consistent with their separate origins from different outcrossing ancestors .
our analysis included both of these selfing variants, an individual of an outcrossing l-morph from northeastern brazil, and a selfing individual of e. paradoxa, the sister species of e. paniculata  <cit> .
we included e. paradoxa to serve as an additional selfing phenotype and as an outgroup for future studies of molecular evolution.
the approaches we describe demonstrate the utility of short-read sequencing for quantifying variation in gene expression among related samples.
methods
sampling & rna preparation
we selected the four plants used in our study from glasshouse collections maintained at the university of toronto.
the plants were originally obtained by germinating open-pollinated seed collected from field populations at the following localities: outcrossing l-morph b <dig>  fortaleza, ceará, n.e.
brazil; selfing m-morph j <dig>  georges plain, westmoreland, jamaica; selfing l-morph n <dig>  rio las lajas, rivas, nicaragua; e. paradoxa, patos, paraíba, n.e.
brazil.
we collected fresh tissue from different stages of bud and flower maturation to sample as many of the genes expressed during development as possible.
to standardize sampling across all four genotypes flower buds were classified into six sizes , with multiple buds for each stage and each stage represented equally among all individuals.
following bud and flower removal, samples were immediately frozen in liquid nitrogen to avoid rna degradation.
we extracted rna from pooled bud samples from each individual using the invitrogen  trizol high salt precipitation extraction protocol.
we visualized rna extracts on a gel to provide an initial assessment of quality and quantified using a spectrophotometer.
sequencing
we provided  <dig> μg of total rna to the center for analysis of genome evolution & function  at the university of toronto  for sequencing.
the rna was sequenced using the illumina  mrna-seq, paired-end protocol on a genome analyzer, gaii, for  <dig> cycles.
this resulted in an average of ~ <dig>  ×  <dig> total reads per sample or ~ <dig>  ×  <dig> bp of sequence per sample .
it is common for the quality of bases from the 3' end of illumina reads to drop in quality, we therefore trimmed the 3' end of reads when the phred quality score dropped below q =  <dig>  for two consecutive bases.
in addition, we also trimmed all 5' and 3' stretches of ambiguous 'n' nucleotides.
trimming resulted in reads with a mean length of  <dig>  bp across all samples and a minimum length of  <dig> bp was applied during sequence trimming.
the number of reads and total amount of sequence generated for each sample of eichhornia.
values are shown for both raw sequence and the sequence after trimming from reads with low quality bases.
de novo assembly
we performed de novo assembly on each sample separately using software packages designed for short read sequence assembly including abyss  <cit> , edena  <cit> , velvet  <cit>  and oases .
to choose the optimal parameters for each method, we used a combination of blastx searches of the ncbi non-redundant protein database , and summary statistics of the assemblies .
consideration of the summary statistics led us to finally choose oases, which generated the longest assembled ests, with the best hits to nr in terms of low e-values.
oases is a program designed as an extension of velvet, specifically released for assembly of transcriptome sequences.
unlike the other software mentioned above, oases handles the uneven coverage of contigs due to variation in expression levels of the transcripts in the sample.
we assembled each sample using the same assembly parameters .
a consequence of the algorithm in the the version of oases we used was a tendency to generate identical or near-identical contigs, possibly due to allelic variants or sequencing errors.
to lower redundancy in the dataset we removed these by comparing each transcriptome assembly to itself using blast  <cit> .
any pair of contigs that were > 99% identical over 95% of the length of the shorter contig were collapsed into a single contig.
consensus transcriptome generation
to create a reference transcriptome we conducted a 'four-way' reciprocal blast  to identify all orthologous sequences.
the goal here was to identify sequences that may not show similarity to other known proteins or ests, but which are expressed in more than one sample.
this procedure allowed us to confirm a large proportion of our transcripts without having to rely on comparative searches to distantly related species.
in addition, we were able to generate longer consensus sequences when one of the reciprocal best blast sequences was longer than the others.
this was implemented using a custom biopython script  <cit>  and blast.
we compared each of the four individual redundancy-reduced transcriptome assemblies to each other using blastn .
reciprocal best blast hits found in more than two samples were then placed into groups and aligned using muscle  <cit>  to generate a consensus sequence.
we defined a number of criteria to identify orthologous sequences including minimum alignment length , minimum sequence identity , and minimum alignment proportion .
this last criterion was used to avoid alternatively spliced transcripts or incompletely aligned contigs being collapsed in an alignment.
after generating the consensus sequences with reciprocal blast we identified unaligned sequences that aligned well to the ortholog groups, but may not have been >  <dig> bp.
these sequences were incorporated into the consensus only when the contig extended the length of the consensus sequence, and had >  <dig> % identity over >  <dig> bp with no unalignable segments.
due to low coverage or repetitive elements within coding loci it is possible that separate contigs are fragments of a single protein.
to reduce fragmentation and recover longer coding sequences we aligned each contig to all unique oryza sativa  proteins using blastx.
we used o. sativa because it is the closest related plant for which an extensive set of protein sequences is available.
this allowed us to identify consensus sequences that probably belong to the same protein and assemble them into a single contig.
we aligned sequences that were potentially from the same protein enabling an elongated consensus to be generated.
only a small number of contigs were found to be potentially fragments of longer ests  and all of the alignments made in sequencher  <dig>  were verified manually to ensure that no gaps, or mismatches were introduced.
after we assembled the consensus of all potential orthologs we identified sequences that were not included in these groups, but had homologs in other species .
we compared each singleton against nr and those over the size threshold of  <dig> bp and with a strong blastx hit  were included in the reference sequence along with all potential orthologs identified with our reciprocal blast scheme.
to ensure that there was no remaining redundancy of transcripts in the consensus we used the same technique to reduce redundancy, as outlined above.
the assembled sequences have been deposited in the ncbi's transcriptome shotgun assembly  database.
genotype calling & snp detection
for each sample we mapped the original short reads, trimmed using the base qualities as outlined above, to the consensus transcriptome using the burrows-wheeler aligner  version  <dig> .
<dig>  <cit> .
by mapping the reads back to the consensus transcriptome we were able to more accurately estimate coverage and use counts of nucleotides at each position to call genotypes.
bwa allowed us to vary the number of mismatches between the reference and aligned reads.
we tested a number of parameter values for alignment and all analyses presented here use a value of n =  <dig> , where n is the fraction of missing alignments given 2% uniform base error rate.
to generate a final sequence for each individual we included all loci where coverage was on average greater than five fold across the locus.
furthermore, within loci only sites where coverage exceeded five were used to call genotypes, the other sites were marked as ambiguous.
to generate genotypes for each sample we used the genotype calling method implemented in the software maq  <cit> .
this method takes into account the counts of different bases at each site, as well as the quality of each base and the mapping quality of the sequence read.
maq uses a bayesian statistical model to compare the inferred genotype to the original reference.
to call genotypes we used a threshold 'consensus quality' cutoff of q >  <dig> , where q is the phred-scaled probability that the consensus genotype call is wrong  <cit> .
sites for which we could not determine the consensus, with at least this level of confidence, were marked as ambiguous.
to detect potential errors in read mapping we assumed that selfing genotypes were largely homozygous, and therefore the presence of heterozygous sites in multiple selfing genotypes may indicate errors in read mapping .
these loci were excluded from downstream analyses.
true heterozygosity was therefore estimated at loci where there was no evidence of read mapping errors.
we used a subset of the identified loci that were shared between all four samples to assess the number of snps between pairs of genotypes.
we also calculated nucleotide polymorphism values, θw  <cit>  for the three e. paniculata sequences
measurement of gene expression
in addition to generating a consensus sequence, we also used abundance of reads derived from each locus to estimate gene expression.
we calculated the number of fragments per kilobase per million fragments mapped  with the program cuffdiff from the package cufflinks v  <dig>   <cit> .
this program estimates confidence intervals around expression estimates of each transcript using a bayesian inference method and will identify significant differences in expression using a fdr control and benjamini-hochberg correction for multiple tests  <cit> .
to compare expression differences among selfing genotypes, we identified loci in which both e. paniculata selfing genotypes, or all three selfing genotypes , differed in expression relative to the outcrossing genotype of e. paniculata from brazil.
to estimate the overall similarity in expression we calculated the correlation of fpkm among pairs of plants after log transformation so that the data fitted a normal distribution.
to test for significant differences among correlation coefficients we bootstrapped the data  to estimate 95% confidence intervals.
functional annotation
to functionally annotate each gene and to assess the quality of our assembly we used the gene ontology  based annotation suite blast2go v <dig> .
<dig>  <cit> .
blast2go allowed us to identify similarity of the sequences in our reference transcriptome to known and predicted proteins, and to assign each of the sequences go terms that were associated with the proteins found by blastx.
we searched all  <dig>  sequences against the non-redundant protein database  with maximum e-value =  <dig> × 10- <dig>  blast2go assigned go terms using a pro-similarity-hit-filter of  <dig>  an annotation cut-off of  <dig> and a go weight of  <dig>  we conducted enrichment analyses to test for an excess or paucity of gene classes  in test sets relative to the whole reference transcriptome.
these included: 1) genes absent in the two selfing e. paniculata genotypes ; 2) genes absent in the outcrossing genotype , and 3) genes with low expression in the two selfing e. paniculata genotypes, and 4) genes with high expression in the two selfing e. paniculata genotypes.
additionally, we repeated these contrasts with e. paradoxa included in all sets with the jamaican and nicaraguan selfers.
however, the results of these analyses were not informative because there were no over-represented gene classes in each group, and the results are therefore not presented here.
all comparisons were implemented in blast2go, which uses a fisher's exact test to determine significance after controlling for multiple tests with a false discovery rate  =  <dig> .
assessment of the accuracy of est assembly
to assess the quality of our assembly, we compared the ests assembled in our consensus transcriptome with a set of  <dig> unique ests sequenced from e. paniculata with sanger sequencing.
these 'sanger-ests' were sampled from a cdna library generated from leaf and floral tissue .
we assembled and aligned forward and reverse strands of each est using sequencher  <dig>  and edited chromatographs and alignments manually.
we compared the two sets of ests to identify conflicts, in which pairs from each est set share a significant portion of their sequence but cannot be aligned over another overlapping portion.
this could be an indication of our illumina ests having been assembled incorrectly, creating chimeric sequences of distinct transcripts.
because the set of sanger ests was not exhaustive, we attempted to identify well-assembled ests from the subset with conflicting alignments by comparing them to the nr protein database to determine whether there is a known protein that covered the full length of the illumina est.
RESULTS
assembly & consensus transcriptome generation
de novo assemblies of each sample using oases resulted in an average of  <dig>   contigs per sample, totaling approximately  <dig>  mbp of sequence for each individual.
many of the sequences were small, resulting in an n <dig> size of  <dig> bp and mean contig size of  <dig> bp .
after removing very similar sequences and contigs that were shorter than  <dig> bp, the mean number of contigs per sample was  <dig> , totaling on average  <dig>  mbp/sample.
the four-way reciprocal blast scheme returned consensus sequences of all contigs present in at least two samples.
this resulted in a nearly two-fold reduction in the total number of contigs to  <dig> , along with an increase in the n <dig> to  <dig> bp and only a slight decrease in the total amount of sequence in the transcriptome, relative to each individual oases assembly .
we attempted to improve our consensus transcriptome by incorporating the remaining contigs into the consensus and by joining contigs that were fragments of the same protein.
these steps had only a minor effect on our consensus, increasing the n <dig> to  <dig> bp, the total length to  <dig>  mbp and decreasing the total number of contigs to  <dig> .
blast identified all the unincorporated contigs that had similarity to known proteins.
<dig> contigs over  <dig> bp had blastx hits with e <  <dig> × 10- <dig> to known proteins and were added to the consensus transcriptome.
after we removed all remaining redundancy, the final consensus transcriptome consisted of  <dig>  contigs, representing  <dig>  mbp of sequence with an n <dig> of  <dig> bp and a mean contig size of  <dig> bp.
the final consensus transcriptome resulted in a nearly two-fold increase in the n <dig>  was comprised of ~ <dig>  fewer contigs than any of the original assemblies, and represented a slightly larger total transcriptome length .
<dig> singletons over  <dig> bp with blastx hits where e <  <dig> × 10- <dig> are included in the final consensus
for the first stage all statistics are the average across the four independently sequenced genotypes.
genotype calling & snp detection
using our consensus transcriptome as a reference, we mapped the original short sequence reads for each sample with the software bwa- <dig> .
<dig>  on average there was  <dig>  reads covering each position in the reference with a standard deviation of  <dig>  .
we generated, on average,  <dig>  contigs for each sample.
<dig>  of the  <dig>  original contigs were found in all four samples and only  <dig> of the reference contigs were not recovered in any of the samples, likely due to low coverage.
few loci were unique to any single e. paniculata genotype,  <dig> sequences were shared by all e. paniculata samples and not found in e. paradoxa.
moreover,  <dig> loci were unique to e. paradoxa.
results of read mapping using the original trimmed short reads for each of the four eichhornia samples against the reference transcriptome
we identified heterozygous loci and potential read mapping errors as loci with one or more bases called as heterozygotes.
assuming that selfing genotypes are largely homozygous, the presence of heterozygous sites in multiple samples of selfers may indicate errors in read mapping .
we identified  <dig>  loci where there was no evidence of read mapping errors,  <dig> loci with some evidence for read mapping errors, and  <dig> loci were expressed in either one or zero selfers, precluding the application of this test.
for the loci in which there was no evidence for read mapping errors, the number of heterozygous loci was highest in the outcrossing e. paniculata genotype  compared to the two selfing genotypes from jamaica  and nicaragua .
eichhornia paradoxa had an intermediate number  of heterozygous loci.
to detect the number of snps between pairs of genotypes, we selected a conservative set of  <dig>  loci  that were expressed in all four individuals and were homozygous in all selfers .
the outcrossing brazilian and selfing nicaraguan genotypes had the fewest divergent sites .
intraspecific variation in e. paniculata was substantially lower than the divergence of each e. paniculata sample to e. paradoxa.
values reflect the polymorphism statistics for  <dig> loci that were present in all four samples and homozygous in the three selfing genotypes.
values below the diagonal show the total number of snps and values above the diagonal show the number snps per bp .
functional annotation
<dig>  of  <dig>  contigs  had similarity to known proteins in the ncbi nr database, with a cutoff of e <  <dig> × 10- <dig>   <dig> loci had alignments which covered more than 80% of the top protein hits and  <dig>  of the query sequences were at least 80% covered by their best protein hit .
blast2go assigned a functional annotation to  <dig>  of the loci .
within the broad go category 'cellular components', over a third of the sequences were localized to the plastid,  <dig> % to the mitochondrion and 17% to the nucleus .
a number of other cellular components made up the remaining  <dig>  % of the annotated loci.
within the broad go category 'biological process' the two most common type of genes were those involved in cellular  and metabolic  processes .
<dig> genes that were identified are involved in reproductive processes including flower development  and pollination .
gene expression
using fpkm to measure gene expression, we found significant correlation in expression among our samples .
as expected, the correlation of each of the three e. paniculata samples with e. paradoxa was lower  compared with the correlation of e. paniculata genotypes with one another.
the two independently derived selfing genotypes were slightly more correlated , but not significantly more so than either was to the outcrossing genotype from brazil .
there were  <dig> genes that were significantly up-regulated in all three selfing genotypes compared with the outcrosser,  <dig> of which were involved in reproduction.
a similar number of genes  were down-regulated in the selfers relative to the outcrosser,  <dig> of these genes were involved in pollination or flower development .
* p <  <dig> , ** p <  <dig> , ***p <  <dig> , ****p <  <dig>
<dig> measured as the mean ratio of expression  in each selfer relative to the outcrossing genotype
<dig> these genes were differentially expressed in the selfing genotypes of e. paniculata only
gene ontology enrichment tests
we investigated whether there was an excess or paucity of particular gene classes that were differentially expressed in the two selfing genotypes of e. paniculata compared to the outcrosser.
in the up-regulated genes of the e. paniculata selfers there were  <dig> go categories that were significantly enriched.
however, using a two-tailed fisher's exact test we found that only genes involved in photosynthesis  were significantly over-represented and the remaining  <dig> classes of genes were under-represented.
in the genes that were expressed at a lower level or were completely absent in the selfers,  <dig> of  <dig> and  <dig> of  <dig> go classes were significantly over-represented.
significantly, of the  <dig> classes of genes in the two e. paniculata selfers that were over-represented among the down-regulated genes, five are related to the regulation of cellular structure and development .
in the e. paniculata selfers, the genes involved in flower development were significantly under represented in both higher and lower expression genes.
assessment of est assembly accuracy
comparisons of our consensus transcriptome ests to a set of  <dig> sanger sequenced ests revealed  <dig> illumina ests  that conflicted with their best sanger est match.
these  <dig> ests each had regions that could not be aligned with the full length of the best blast hit from the sanger ests and also did not have full-length blast hits to known proteins in nr.
however, all  <dig> had highly significant hits to proteins in nr.
none of the  <dig> ests appeared to be chimeric assemblies of distinct proteins that may have resulted from errors during assembly.
moreover,  <dig> of the  <dig> ests had sequence flanking their protein-coding region with significant similarity to known 5' and 3' untranslated regions from other monocot genomes.
discussion
we assembled ~ <dig> mbp of transcriptome sequence in each of four individuals of two species of eichhornia.
the data represent an important genetic resource of nearly  <dig>  transcripts, many of which are common to all four samples.
further, using read mapping to a consensus transcriptome we have generated statistically informed genotypes for each individual in our study.
by choosing to extract rna from buds and mature flowers we were able to recover many genes involved in reproduction and floral development, some of which are likely to provide future insights on genetic changes to floral traits governing mating-system variation.
we now compare our assembly and analysis with previous de novo transcriptome sequencing projects and briefly review some of the challenges and interpretations specific to our assembly.
in addition, we also discuss the utility of short-read sequencing for characterizing genetic changes in transcriptomes and the expression level of different loci.
assembly & consensus transcriptome
our study is among the first published attempts at a de novo transcriptome assembly using short-read  sequencing.
although there are now many studies using next generation transcriptome sequencing, most have used the roche  <dig> platform .
as previously mentioned, this platform has the advantage of longer reads but at the expense of less sequence data per run.
longer reads may be critical for resolving assembly challenges associated with repetitive elements, such as gene duplications, allelic differences and alternative splicing .
however, despite using shorter reads, our assembly is comparable to other published transcriptomes, which at this time average  <dig>  contigs , a number similar to our own .
in addition, despite the longer read lengths , of previously published transcriptome studies using the roche  <dig> platform, our assembly generated contigs that were over twice as long averaging  <dig> bp.
one of the reasons that this is the case is because the illumina generates greater depth of sequencing thus ensuring more complete coverage of the transcriptome.
although summaries of the distribution of contig lengths are informative, the ultimate goal of transcriptome assembly is not long sequences, but accurate assembly of full-length sequences.
however, it is difficult to assess the success of an assembly without a priori knowledge of the transcriptome.
one metric that may be informative is the proportion of contigs that have significant similarity to known proteins.
it is difficult to compare this measure across studies because each reports slightly different results using different blast parameters and databases.
however, nearly 87% of our contigs had matches in nr and this value is as high or higher than all other comparable statistics reported in other de novo assemblies.
another useful metric is the proportion of the contig and its corresponding blast hit that align to one another .
<dig>  contigs cover greater than 75% of their best blast hit and  <dig>   contigs are fully covered by their best blast hit.
this means that although we assembled a large number of full-length proteins, many of the contigs appear to be fragments of larger proteins.
one explanation is that gene duplication or alternative splicing creates repetitive elements in the assembly and these cannot be resolved.
although we found a fraction of our ests  that had conflicts with sanger sequenced ests, this is likely an overestimate of assembly error because some of these conflicts could result from paralogy or alternative splicing.
there was no evidence that any of the  <dig> conflicting ests were chimeric assemblies of two or more proteins and most of the conflicts  appear to be the result of misaligned untranslated regions flanking coding regions.
it is possible that some of the discrepancies between the ests we assembled and sanger ests or known proteins are true differences, for example, paralogous transcripts or alternatively spliced isoforms.
genotype calling & snp detection
one of the major challenges in dealing with very short sequence reads is that they must be assembled into longer contigs based on overlap with other reads.
the algorithms used in many de novo assemblers, including velvet  <cit> , may misinterpret small differences between alleles  or gene duplicates as sequencing errors.
if this occurs they can be 'collapsed' or purged from the assembly.
although the program oases has begun to address this problem for transcriptome data, we chose to use read mapping to a consensus transcriptome because it allowed us to use allele frequencies at each site to statistically determine the genotype of each individual  <cit> .
from this approach we obtained, on average, more than  <dig> reads for each position to inform genotype calls.
this allowed us to generate sequences for ~ <dig>  loci/genotype,  <dig>   <dig> of which were found in all samples.
an additional benefit of using this approach is that we were able to identify heterozygous loci and potential read mapping errors.
of the  <dig>  consensus sequences,  <dig> heterozygous loci were identified.
as expected, the fraction of heterozygous loci was highest in the outcrossing genotype of e. paniculata from brazil.
the selfing genotypes from jamaica and nicaragua appear to retain some residual heterozygosity, despite their predominantly autogamous mating systems.
we also found evidence of read mapping errors in  <dig> loci where more than one selfer appeared heterozygous.
possible explanations for these read-mapping errors include, sequencing errors, alternative splicing and most likely paralogy.
with our current data and the available methods we have no way of determining their relative contribution to read mapping errors.
sequencing of multiple paralogous transcripts will generate short reads that are similar or identical to many other reads derived from different loci.
as a result, when there has not been sufficient divergence between gene copies, reads may be erroneously mapped to the reference sequence.
one possible explanation for false heterozygosity is that repetitive elements, for example, conserved motifs in gene families, may be difficult or impossible to assemble into long contigs using current technology.
as a result, a fraction of the original short reads may not have been assembled by oases, or included in our consensus transcriptome.
however, if they share enough similarity with a paralog in the transcriptome they may be incorrectly mapped.
this could explain, in part, why the original oases assemblies contain so many short contigs .
if there is divergence between paralogous loci, incorrectly mapped reads may create a signature similar to heterozygosity.
for future analyses it is critical that potentially paralogous sequences are identified because evolutionary inferences from non-orthologous genes are misleading.
although our approach does not allow us to unambiguously characterize all paralogous sequences it has provided a useful method for detecting single copy transcripts.
using a conservative set of  <dig> transcripts for which there was no evidence of paralogy, based on homozygosity in all three selfers and presence in all four genotypes, we determined the number of snps between each pair of genotypes.
as expected, e. paniculata samples were more differentiated from e. paradoxa than with one another .
the patterns of divergence among e. paniculata genotypes reflect relationships previously reported .
specifically, the nicaraguan selfing genotype is more similar to the outcrossing brazilian genotype than it is to the selfing jamaican genotype, despite more similar biogeographical origins and mating systems.
this suggests that the nicaraguan population is more closely related, or more recently derived from the brazilian population.
further, when we calculate nucleotide polymorphism across  <dig> sequences for the three e. paniculata genotypes the value we obtained  is comparable to our previously published species-wide estimate of total diversity , based on  <dig> nuclear-derived est loci assayed in samples of  <dig> individuals from  <dig> populations.
this evidence supports the validity of our snp detection method.
expression & enrichment
there was a weak trend indicating that the selfing genotypes of e. paniculata were more correlated in gene expression; however, this difference was not significant using bootstrapping to generate 95% confidence intervals.
all three genotypes of e. paniculata retain highly correlated gene expression despite phenotypic divergence and geographic isolation.
the slight elevation in correlated expression between the two selfing genotypes of e. paniculata may be caused by a small number of genes that are differentially expressed in both selfing genotypes , although overall patterns of expression during flower development appear to remain largely conserved.
this may be because we combined all stages of flower development in our assays and a more careful dissection of expression in each stage individually could potentially reveal different patterns of gene expression in selfers and outcrossers.
enrichment tests of our annotated transcriptomes demonstrated that genes that were differentially expressed in selfers exhibit a paucity of particular gene classes.
this can be interpreted as the conservation of expression of these gene functions, which are rarely differentially expressed.
of the  <dig> go categories found to be significantly enriched among all differentially expressed genes only  <dig> were found to be over-abundant, and  <dig> of these  <dig> categories were over-represented in genes absent or expressed at lower levels in the two selfers.
therefore, it appears that many of the differences common to the selfing lineages of e. paniculata are associated with reductions of gene expression in floral tissue.
this may be related to the convergence of floral traits in these two lineages, both of which have much smaller, less pigmented flowers, with reduced pollen production compared to the outcrossing genotype.
<dig> of the  <dig> gene classes that are expressed at a low level in the selfers are under-represented but  <dig> of the  <dig> over-represented go categories were associated with the regulation of cell development and structure.
selfing flowers display modified stamen positions and floral instability including twisted, fused or missing perianth parts, whereas outcrossing plants rarely display these floral modifications  <cit> .
it is possible that the changes in gene expression we have documented influence the regulation of cell growth and division and are responsible for changes to floral morphology that characterize selfing populations.
if so, these regulatory loci could be used as a set of candidate genes to investigate aspects of the evolution of the selfing syndrome.
differentially expressed floral genes in selfers
by sampling different stages of floral development up to and including anthesis we were able to sequence and annotate a large number of florally expressed genes.
in total  <dig> genes with the go annotation 'reproduction' were identified, which is a large fraction of the number reported for arabidopsis thaliana in which  <dig> genes for reproduction have been documented  <cit> .
the lower number of genes in our annotation is not unexpected because we did not include tissue from reproductive stages after flowering, such as fruit and seed development.
within the go category 'reproduction' we found  <dig> genes involved in floral development, similar to the number that has been annotated in a. thaliana .
of particular significance are the floral development genes that are differentially expressed in the three selfing genotypes , several of which affect structures that are modified in the selfers.
anther development and filament elongation are influenced by afb <dig>  <cit> , arf <dig>  arf <dig>  <cit> , bam <dig>  <cit> , gid1c,  <cit> myb <dig>  <cit> , pgp <dig>  <cit>  and ptr <dig>  <cit>  genes in a. thaliana, and pollen maturation and pollen tube growth are altered by afb <dig>  csdl <dig>  <cit> , cer <dig>  <cit> , pop <dig>  <cit> , myb <dig>  <cit> and tir <dig>  <cit> .
erl <dig> plays an important roles in normal anther lobe formation and anther cell differentiation  <cit>  and mutants in the erecta gene family have reduced lateral organ size and abnormal flower development, including defects in petal polar expansion and carpel elongation  <cit> .
we also identified a number of genes involved in flowering time including dcl <dig>  <cit> , pgi <dig>  <cit> , shk <dig>  <cit> , suf <dig>  <cit> .
lastly, all of the differentially expressed genes that influence ovule development were significantly up-regulated, including arf <dig>  arf <dig>  <cit> , ova <dig>  <cit> , pep  <cit> .
significantly, most of the candidate genes discussed above cause alterations to attractive structures  and male function  consistent with the relaxation of selection for showy flowers, reduced allocation to pollen production and the loss of herkogamy  through filament elongation of stamens.
in contrast, the requirement for functional ovules to maintain seed fertility in selfers may explain the apparent absence of changes to gene expression in female traits.
CONCLUSIONS
we have shown that short-read sequencing can be used to characterize the transcriptomes of multiple individuals for use in comparative studies.
we were able to assemble as many contigs as other sequencing methods by de novo assembly, but our contigs were on average substantially longer.
by comparing sequences from related individuals we generated a consensus transcriptome.
this allowed us to make snp genotype calls and provided a method for detecting paralogous sequences.
discerning among copies of paralogous sequences presents a major challenge to the current technology and requires either technological or analytical solutions to discern among different members of gene families or duplicates.
however, despite these complications we annotated > 80% of contigs and identified thousands of putative orthologs, many of which are differentially expressed.
we identified  <dig> genes that were differentially expressed in selfers and which have developmental functions that suggest a role in the evolution of the selfing syndrome.
this sequencing effort has generated a valuable resource of coding dna for a non-model species.
the transcriptome sequences will help in future studies of changes in the genetic architecture involved in the transition from outcrossing to selfing and also for identifying the genes controlling heterostyly.
authors' contributions
the study was conceived by rwn and schb.
the molecular lab work was carried out by rwn.
sequence assembly and annotation and expression analyses were carried out by rwn.
ms contributed programming expertise and to methods for implementing genotyping.
all authors read and approved the final manuscript.

 among the most prevalent evolutionary transitions in plants is the shift from cross-fertilization to predominant self-fertilization among numerous angiosperm lineages. this change in mating system has important consequences for many aspects of the biology of selfing taxa including population genetic structure, colonizing ability, genome evolution and the morphology of flowers. the loss of floral mechanisms that reduce the incidence of self-fertilization results in high rates of autogamous selfing, leading to the evolution of the'selfing syndrome'. although studies of the causes and consequences of cross- and self-fertilization in flowering plants have a long and venerable history, beginning with darwin's seminal work, relatively little is known about the underlying molecular changes that accompany the transition from outcrossing to selfing and genomic analyses of related outcrossing and selfing plants are in their infancy. recent technological advances in dna sequencing technology have removed several limitations associated with gathering large amounts of genomic data from non-model organisms, providing opportunities for detailed investigation into the genomics of mating-system variation and evolution. although assembling large eukaryotic genomes, de novo, may not yet be practical, sets of expressed genes or transcriptomes present a viable and attractive alternative to population genetic analyses of whole genome sequences. transcriptomes represent a fraction of the total genome in size, contain fewer repetitive elements, and by selecting specific tissues they can be enriched for genes relevant to the particular aims of the research. in addition, if the rna sample is not normalized the relative abundance of different reads has been shown to accurately reflect the expression level of transcripts in the tissue. despite these potential advantages there remain a number of challenges for de novo transcriptome assembly, including gene duplication or paralogy, heterozygosity and alternative splicing, each of which require careful consideration. there are relatively few studies to date involving de novo transcriptome assembly in non-model organisms. so far the majority have used the roche gs platform which has the advantage of longer reads, but produces a fraction of the total amount of sequence produced per instrument per run compared with the illumina gaii platform. therefore, to maximize the coverage for rare transcripts, cdna samples are typically normalized. as a result, studies using are not able to estimate expression levels of different ests, as this requires deep sequencing of non-normalized cdna. further, due to lower sequencing depth many transcripts are represented by a single read, and others by very few reads. this can create problems in accurately distinguishing snps from errors, and in retrieving orthologous transcripts for sequence comparisons across experiments or species. therefore, new methods are required to generate and assemble large datasets, many of which currently consist of substantially shorter reads. here, we present de novo floral transcriptome assemblies using short read sequencing of four individual plants of neotropical eichhornia species that vary in floral morphology and mating system. the samples include three individuals of e. paniculata, including two from independently derived selfing populations and the third an outcrosser. eichhornia paniculata is an annual diploid that has been the subject of detailed studies over the past two decades on the ecology and genetics of mating-system variation. populations of e. paniculata are largely concentrated in northeastern brazil, with smaller foci in jamaica and cuba and isolated localities in nicaragua and mexico. populations in brazil are largely outcrossing and possess the sexual polymorphism tristyly, which promotes cross-pollination among the three floral morphs. morphological, genetic and biogeographical evidence indicates that tristyly in e. paniculata has broken down on multiple occasions resulting in independently derived selfing populations. populations in jamaica are largely composed of selfing variants of the mid-styled morph in which short-level stamens are elongated to a position adjacent to mid-level stigmas resulting in autonomous self-fertilization. in contrast, plants in mexico and nicaragua are selfing variants of the long-styled morph with a different arrangement of their sexual organs. although both variants possess the selfing syndrome, comparisons of molecular variation at est-derived nuclear loci indicate a high level of differentiation consistent with their separate origins from different outcrossing ancestors. our analysis included both of these selfing variants, an individual of an outcrossing l-morph from northeastern brazil, and a selfing individual of e. paradoxa, the sister species of e. paniculata. we included e. paradoxa to serve as an additional selfing phenotype and as an outgroup for future studies of molecular evolution. the approaches we describe demonstrate the utility of short-read sequencing for quantifying variation in gene expression among related samples. methods sampling & rna preparation we selected the four plants used in our study from glasshouse collections maintained at the university of toronto. the plants were originally obtained by germinating open-pollinated seed collected from field populations at the following localities: outcrossing l-morph b fortaleza, ceará, n.e. brazil; selfing m-morph j georges plain, westmoreland, jamaica; selfing l-morph n rio las lajas, rivas, nicaragua; e. paradoxa, patos, paraíba, n.e. brazil. we collected fresh tissue from different stages of bud and flower maturation to sample as many of the genes expressed during development as possible. to standardize sampling across all four genotypes flower buds were classified into six sizes, with multiple buds for each stage and each stage represented equally among all individuals. following bud and flower removal, samples were immediately frozen in liquid nitrogen to avoid rna degradation. we extracted rna from pooled bud samples from each individual using the invitrogen trizol high salt precipitation extraction protocol. we visualized rna extracts on a gel to provide an initial assessment of quality and quantified using a spectrophotometer. sequencing we provided μg of total rna to the center for analysis of genome evolution & function at the university of toronto for sequencing. the rna was sequenced using the illumina mrna-seq, paired-end protocol on a genome analyzer, gaii, for cycles. this resulted in an average of ~ × total reads per sample or ~ × bp of sequence per sample. it is common for the quality of bases from the 3' end of illumina reads to drop in quality, we therefore trimmed the 3' end of reads when the phred quality score dropped below q = for two consecutive bases. in addition, we also trimmed all 5' and 3' stretches of ambiguous'n' nucleotides. trimming resulted in reads with a mean length of bp across all samples and a minimum length of bp was applied during sequence trimming. the number of reads and total amount of sequence generated for each sample of eichhornia. values are shown for both raw sequence and the sequence after trimming from reads with low quality bases. de novo assembly we performed de novo assembly on each sample separately using software packages designed for short read sequence assembly including abyss, edena, velvet and oases. to choose the optimal parameters for each method, we used a combination of blastx searches of the ncbi non-redundant protein database, and summary statistics of the assemblies. consideration of the summary statistics led us to finally choose oases, which generated the longest assembled ests, with the best hits to nr in terms of low e-values. oases is a program designed as an extension of velvet, specifically released for assembly of transcriptome sequences. unlike the other software mentioned above, oases handles the uneven coverage of contigs due to variation in expression levels of the transcripts in the sample. we assembled each sample using the same assembly parameters. a consequence of the algorithm in the the version of oases we used was a tendency to generate identical or near-identical contigs, possibly due to allelic variants or sequencing errors. to lower redundancy in the dataset we removed these by comparing each transcriptome assembly to itself using blast. any pair of contigs that were > 99% identical over 95% of the length of the shorter contig were collapsed into a single contig. consensus transcriptome generation to create a reference transcriptome we conducted a'four-way' reciprocal blast to identify all orthologous sequences. the goal here was to identify sequences that may not show similarity to other known proteins or ests, but which are expressed in more than one sample. this procedure allowed us to confirm a large proportion of our transcripts without having to rely on comparative searches to distantly related species. in addition, we were able to generate longer consensus sequences when one of the reciprocal best blast sequences was longer than the others. this was implemented using a custom biopython script and blast. we compared each of the four individual redundancy-reduced transcriptome assemblies to each other using blastn. reciprocal best blast hits found in more than two samples were then placed into groups and aligned using muscle to generate a consensus sequence. we defined a number of criteria to identify orthologous sequences including minimum alignment length, minimum sequence identity, and minimum alignment proportion. this last criterion was used to avoid alternatively spliced transcripts or incompletely aligned contigs being collapsed in an alignment. after generating the consensus sequences with reciprocal blast we identified unaligned sequences that aligned well to the ortholog groups, but may not have been > bp. these sequences were incorporated into the consensus only when the contig extended the length of the consensus sequence, and had > % identity over > bp with no unalignable segments. due to low coverage or repetitive elements within coding loci it is possible that separate contigs are fragments of a single protein. to reduce fragmentation and recover longer coding sequences we aligned each contig to all unique oryza sativa proteins using blastx. we used o. sativa because it is the closest related plant for which an extensive set of protein sequences is available. this allowed us to identify consensus sequences that probably belong to the same protein and assemble them into a single contig. we aligned sequences that were potentially from the same protein enabling an elongated consensus to be generated. only a small number of contigs were found to be potentially fragments of longer ests and all of the alignments made in sequencher were verified manually to ensure that no gaps, or mismatches were introduced. after we assembled the consensus of all potential orthologs we identified sequences that were not included in these groups, but had homologs in other species. we compared each singleton against nr and those over the size threshold of bp and with a strong blastx hit were included in the reference sequence along with all potential orthologs identified with our reciprocal blast scheme. to ensure that there was no remaining redundancy of transcripts in the consensus we used the same technique to reduce redundancy, as outlined above. the assembled sequences have been deposited in the ncbi's transcriptome shotgun assembly database. genotype calling & snp detection for each sample we mapped the original short reads, trimmed using the base qualities as outlined above, to the consensus transcriptome using the burrows-wheeler aligner version.. by mapping the reads back to the consensus transcriptome we were able to more accurately estimate coverage and use counts of nucleotides at each position to call genotypes. bwa allowed us to vary the number of mismatches between the reference and aligned reads. we tested a number of parameter values for alignment and all analyses presented here use a value of n =, where n is the fraction of missing alignments given 2% uniform base error rate. to generate a final sequence for each individual we included all loci where coverage was on average greater than five fold across the locus. furthermore, within loci only sites where coverage exceeded five were used to call genotypes, the other sites were marked as ambiguous. to generate genotypes for each sample we used the genotype calling method implemented in the software maq. this method takes into account the counts of different bases at each site, as well as the quality of each base and the mapping quality of the sequence read. maq uses a bayesian statistical model to compare the inferred genotype to the original reference. to call genotypes we used a threshold'consensus quality' cutoff of q >, where q is the phred-scaled probability that the consensus genotype call is wrong. sites for which we could not determine the consensus, with at least this level of confidence, were marked as ambiguous. to detect potential errors in read mapping we assumed that selfing genotypes were largely homozygous, and therefore the presence of heterozygous sites in multiple selfing genotypes may indicate errors in read mapping. these loci were excluded from downstream analyses. true heterozygosity was therefore estimated at loci where there was no evidence of read mapping errors. we used a subset of the identified loci that were shared between all four samples to assess the number of snps between pairs of genotypes. we also calculated nucleotide polymorphism values, θw for the three e. paniculata sequences measurement of gene expression in addition to generating a consensus sequence, we also used abundance of reads derived from each locus to estimate gene expression. we calculated the number of fragments per kilobase per million fragments mapped with the program cuffdiff from the package cufflinks v. this program estimates confidence intervals around expression estimates of each transcript using a bayesian inference method and will identify significant differences in expression using a fdr control and benjamini-hochberg correction for multiple tests. to compare expression differences among selfing genotypes, we identified loci in which both e. paniculata selfing genotypes, or all three selfing genotypes, differed in expression relative to the outcrossing genotype of e. paniculata from brazil. to estimate the overall similarity in expression we calculated the correlation of fpkm among pairs of plants after log transformation so that the data fitted a normal distribution. to test for significant differences among correlation coefficients we bootstrapped the data to estimate 95% confidence intervals. functional annotation to functionally annotate each gene and to assess the quality of our assembly we used the gene ontology based annotation suite blast2go v.. blast2go allowed us to identify similarity of the sequences in our reference transcriptome to known and predicted proteins, and to assign each of the sequences go terms that were associated with the proteins found by blastx. we searched all sequences against the non-redundant protein database with maximum e-value = × 10- blast2go assigned go terms using a pro-similarity-hit-filter of an annotation cut-off of and a go weight of we conducted enrichment analyses to test for an excess or paucity of gene classes in test sets relative to the whole reference transcriptome. these included: 1) genes absent in the two selfing e. paniculata genotypes; 2) genes absent in the outcrossing genotype, and 3) genes with low expression in the two selfing e. paniculata genotypes, and 4) genes with high expression in the two selfing e. paniculata genotypes. additionally, we repeated these contrasts with e. paradoxa included in all sets with the jamaican and nicaraguan selfers. however, the results of these analyses were not informative because there were no over-represented gene classes in each group, and the results are therefore not presented here. all comparisons were implemented in blast2go, which uses a fisher's exact test to determine significance after controlling for multiple tests with a false discovery rate =. assessment of the accuracy of est assembly to assess the quality of our assembly, we compared the ests assembled in our consensus transcriptome with a set of unique ests sequenced from e. paniculata with sanger sequencing. these'sanger-ests' were sampled from a cdna library generated from leaf and floral tissue. we assembled and aligned forward and reverse strands of each est using sequencher and edited chromatographs and alignments manually. we compared the two sets of ests to identify conflicts, in which pairs from each est set share a significant portion of their sequence but cannot be aligned over another overlapping portion. this could be an indication of our illumina ests having been assembled incorrectly, creating chimeric sequences of distinct transcripts. because the set of sanger ests was not exhaustive, we attempted to identify well-assembled ests from the subset with conflicting alignments by comparing them to the nr protein database to determine whether there is a known protein that covered the full length of the illumina est.  assembly & consensus transcriptome generation de novo assemblies of each sample using oases resulted in an average of contigs per sample, totaling approximately mbp of sequence for each individual. many of the sequences were small, resulting in an n size of bp and mean contig size of bp. after removing very similar sequences and contigs that were shorter than bp, the mean number of contigs per sample was, totaling on average mbp/sample. the four-way reciprocal blast scheme returned consensus sequences of all contigs present in at least two samples. this resulted in a nearly two-fold reduction in the total number of contigs to, along with an increase in the n to bp and only a slight decrease in the total amount of sequence in the transcriptome, relative to each individual oases assembly. we attempted to improve our consensus transcriptome by incorporating the remaining contigs into the consensus and by joining contigs that were fragments of the same protein. these steps had only a minor effect on our consensus, increasing the n to bp, the total length to mbp and decreasing the total number of contigs to. blast identified all the unincorporated contigs that had similarity to known proteins. contigs over bp had blastx hits with e < × 10- to known proteins and were added to the consensus transcriptome. after we removed all remaining redundancy, the final consensus transcriptome consisted of contigs, representing mbp of sequence with an n of bp and a mean contig size of bp. the final consensus transcriptome resulted in a nearly two-fold increase in the n was comprised of ~ fewer contigs than any of the original assemblies, and represented a slightly larger total transcriptome length. singletons over bp with blastx hits where e < × 10- are included in the final consensus for the first stage all statistics are the average across the four independently sequenced genotypes. genotype calling & snp detection using our consensus transcriptome as a reference, we mapped the original short sequence reads for each sample with the software bwa-. on average there was reads covering each position in the reference with a standard deviation of. we generated, on average, contigs for each sample. of the original contigs were found in all four samples and only of the reference contigs were not recovered in any of the samples, likely due to low coverage. few loci were unique to any single e. paniculata genotype, sequences were shared by all e. paniculata samples and not found in e. paradoxa. moreover, loci were unique to e. paradoxa. results of read mapping using the original trimmed short reads for each of the four eichhornia samples against the reference transcriptome we identified heterozygous loci and potential read mapping errors as loci with one or more bases called as heterozygotes. assuming that selfing genotypes are largely homozygous, the presence of heterozygous sites in multiple samples of selfers may indicate errors in read mapping. we identified loci where there was no evidence of read mapping errors, loci with some evidence for read mapping errors, and loci were expressed in either one or zero selfers, precluding the application of this test. for the loci in which there was no evidence for read mapping errors, the number of heterozygous loci was highest in the outcrossing e. paniculata genotype compared to the two selfing genotypes from jamaica and nicaragua. eichhornia paradoxa had an intermediate number of heterozygous loci. to detect the number of snps between pairs of genotypes, we selected a conservative set of loci that were expressed in all four individuals and were homozygous in all selfers. the outcrossing brazilian and selfing nicaraguan genotypes had the fewest divergent sites. intraspecific variation in e. paniculata was substantially lower than the divergence of each e. paniculata sample to e. paradoxa. values reflect the polymorphism statistics for loci that were present in all four samples and homozygous in the three selfing genotypes. values below the diagonal show the total number of snps and values above the diagonal show the number snps per bp. functional annotation of contigs had similarity to known proteins in the ncbi nr database, with a cutoff of e < × 10- loci had alignments which covered more than 80% of the top protein hits and of the query sequences were at least 80% covered by their best protein hit. blast2go assigned a functional annotation to of the loci. within the broad go category'cellular components', over a third of the sequences were localized to the plastid, % to the mitochondrion and 17% to the nucleus. a number of other cellular components made up the remaining % of the annotated loci. within the broad go category'biological process' the two most common type of genes were those involved in cellular and metabolic processes. genes that were identified are involved in reproductive processes including flower development and pollination. gene expression using fpkm to measure gene expression, we found significant correlation in expression among our samples. as expected, the correlation of each of the three e. paniculata samples with e. paradoxa was lower compared with the correlation of e. paniculata genotypes with one another. the two independently derived selfing genotypes were slightly more correlated, but not significantly more so than either was to the outcrossing genotype from brazil. there were genes that were significantly up-regulated in all three selfing genotypes compared with the outcrosser, of which were involved in reproduction. a similar number of genes were down-regulated in the selfers relative to the outcrosser, of these genes were involved in pollination or flower development. * p <, ** p <, ***p <, ****p < measured as the mean ratio of expression in each selfer relative to the outcrossing genotype these genes were differentially expressed in the selfing genotypes of e. paniculata only gene ontology enrichment tests we investigated whether there was an excess or paucity of particular gene classes that were differentially expressed in the two selfing genotypes of e. paniculata compared to the outcrosser. in the up-regulated genes of the e. paniculata selfers there were go categories that were significantly enriched. however, using a two-tailed fisher's exact test we found that only genes involved in photosynthesis were significantly over-represented and the remaining classes of genes were under-represented. in the genes that were expressed at a lower level or were completely absent in the selfers, of and of go classes were significantly over-represented. significantly, of the classes of genes in the two e. paniculata selfers that were over-represented among the down-regulated genes, five are related to the regulation of cellular structure and development. in the e. paniculata selfers, the genes involved in flower development were significantly under represented in both higher and lower expression genes. assessment of est assembly accuracy comparisons of our consensus transcriptome ests to a set of sanger sequenced ests revealed illumina ests that conflicted with their best sanger est match. these ests each had regions that could not be aligned with the full length of the best blast hit from the sanger ests and also did not have full-length blast hits to known proteins in nr. however, all had highly significant hits to proteins in nr. none of the ests appeared to be chimeric assemblies of distinct proteins that may have resulted from errors during assembly. moreover, of the ests had sequence flanking their protein-coding region with significant similarity to known 5' and 3' untranslated regions from other monocot genomes. discussion we assembled ~ mbp of transcriptome sequence in each of four individuals of two species of eichhornia. the data represent an important genetic resource of nearly transcripts, many of which are common to all four samples. further, using read mapping to a consensus transcriptome we have generated statistically informed genotypes for each individual in our study. by choosing to extract rna from buds and mature flowers we were able to recover many genes involved in reproduction and floral development, some of which are likely to provide future insights on genetic changes to floral traits governing mating-system variation. we now compare our assembly and analysis with previous de novo transcriptome sequencing projects and briefly review some of the challenges and interpretations specific to our assembly. in addition, we also discuss the utility of short-read sequencing for characterizing genetic changes in transcriptomes and the expression level of different loci. assembly & consensus transcriptome our study is among the first published attempts at a de novo transcriptome assembly using short-read sequencing. although there are now many studies using next generation transcriptome sequencing, most have used the roche platform. as previously mentioned, this platform has the advantage of longer reads but at the expense of less sequence data per run. longer reads may be critical for resolving assembly challenges associated with repetitive elements, such as gene duplications, allelic differences and alternative splicing. however, despite using shorter reads, our assembly is comparable to other published transcriptomes, which at this time average contigs, a number similar to our own. in addition, despite the longer read lengths, of previously published transcriptome studies using the roche platform, our assembly generated contigs that were over twice as long averaging bp. one of the reasons that this is the case is because the illumina generates greater depth of sequencing thus ensuring more complete coverage of the transcriptome. although summaries of the distribution of contig lengths are informative, the ultimate goal of transcriptome assembly is not long sequences, but accurate assembly of full-length sequences. however, it is difficult to assess the success of an assembly without a priori knowledge of the transcriptome. one metric that may be informative is the proportion of contigs that have significant similarity to known proteins. it is difficult to compare this measure across studies because each reports slightly different results using different blast parameters and databases. however, nearly 87% of our contigs had matches in nr and this value is as high or higher than all other comparable statistics reported in other de novo assemblies. another useful metric is the proportion of the contig and its corresponding blast hit that align to one another. contigs cover greater than 75% of their best blast hit and contigs are fully covered by their best blast hit. this means that although we assembled a large number of full-length proteins, many of the contigs appear to be fragments of larger proteins. one explanation is that gene duplication or alternative splicing creates repetitive elements in the assembly and these cannot be resolved. although we found a fraction of our ests that had conflicts with sanger sequenced ests, this is likely an overestimate of assembly error because some of these conflicts could result from paralogy or alternative splicing. there was no evidence that any of the conflicting ests were chimeric assemblies of two or more proteins and most of the conflicts appear to be the result of misaligned untranslated regions flanking coding regions. it is possible that some of the discrepancies between the ests we assembled and sanger ests or known proteins are true differences, for example, paralogous transcripts or alternatively spliced isoforms. genotype calling & snp detection one of the major challenges in dealing with very short sequence reads is that they must be assembled into longer contigs based on overlap with other reads. the algorithms used in many de novo assemblers, including velvet, may misinterpret small differences between alleles or gene duplicates as sequencing errors. if this occurs they can be'collapsed' or purged from the assembly. although the program oases has begun to address this problem for transcriptome data, we chose to use read mapping to a consensus transcriptome because it allowed us to use allele frequencies at each site to statistically determine the genotype of each individual. from this approach we obtained, on average, more than reads for each position to inform genotype calls. this allowed us to generate sequences for ~ loci/genotype, of which were found in all samples. an additional benefit of using this approach is that we were able to identify heterozygous loci and potential read mapping errors. of the consensus sequences, heterozygous loci were identified. as expected, the fraction of heterozygous loci was highest in the outcrossing genotype of e. paniculata from brazil. the selfing genotypes from jamaica and nicaragua appear to retain some residual heterozygosity, despite their predominantly autogamous mating systems. we also found evidence of read mapping errors in loci where more than one selfer appeared heterozygous. possible explanations for these read-mapping errors include, sequencing errors, alternative splicing and most likely paralogy. with our current data and the available methods we have no way of determining their relative contribution to read mapping errors. sequencing of multiple paralogous transcripts will generate short reads that are similar or identical to many other reads derived from different loci. as a result, when there has not been sufficient divergence between gene copies, reads may be erroneously mapped to the reference sequence. one possible explanation for false heterozygosity is that repetitive elements, for example, conserved motifs in gene families, may be difficult or impossible to assemble into long contigs using current technology. as a result, a fraction of the original short reads may not have been assembled by oases, or included in our consensus transcriptome. however, if they share enough similarity with a paralog in the transcriptome they may be incorrectly mapped. this could explain, in part, why the original oases assemblies contain so many short contigs. if there is divergence between paralogous loci, incorrectly mapped reads may create a signature similar to heterozygosity. for future analyses it is critical that potentially paralogous sequences are identified because evolutionary inferences from non-orthologous genes are misleading. although our approach does not allow us to unambiguously characterize all paralogous sequences it has provided a useful method for detecting single copy transcripts. using a conservative set of transcripts for which there was no evidence of paralogy, based on homozygosity in all three selfers and presence in all four genotypes, we determined the number of snps between each pair of genotypes. as expected, e. paniculata samples were more differentiated from e. paradoxa than with one another. the patterns of divergence among e. paniculata genotypes reflect relationships previously reported. specifically, the nicaraguan selfing genotype is more similar to the outcrossing brazilian genotype than it is to the selfing jamaican genotype, despite more similar biogeographical origins and mating systems. this suggests that the nicaraguan population is more closely related, or more recently derived from the brazilian population. further, when we calculate nucleotide polymorphism across sequences for the three e. paniculata genotypes the value we obtained is comparable to our previously published species-wide estimate of total diversity, based on nuclear-derived est loci assayed in samples of individuals from populations. this evidence supports the validity of our snp detection method. expression & enrichment there was a weak trend indicating that the selfing genotypes of e. paniculata were more correlated in gene expression; however, this difference was not significant using bootstrapping to generate 95% confidence intervals. all three genotypes of e. paniculata retain highly correlated gene expression despite phenotypic divergence and geographic isolation. the slight elevation in correlated expression between the two selfing genotypes of e. paniculata may be caused by a small number of genes that are differentially expressed in both selfing genotypes, although overall patterns of expression during flower development appear to remain largely conserved. this may be because we combined all stages of flower development in our assays and a more careful dissection of expression in each stage individually could potentially reveal different patterns of gene expression in selfers and outcrossers. enrichment tests of our annotated transcriptomes demonstrated that genes that were differentially expressed in selfers exhibit a paucity of particular gene classes. this can be interpreted as the conservation of expression of these gene functions, which are rarely differentially expressed. of the go categories found to be significantly enriched among all differentially expressed genes only were found to be over-abundant, and of these categories were over-represented in genes absent or expressed at lower levels in the two selfers. therefore, it appears that many of the differences common to the selfing lineages of e. paniculata are associated with reductions of gene expression in floral tissue. this may be related to the convergence of floral traits in these two lineages, both of which have much smaller, less pigmented flowers, with reduced pollen production compared to the outcrossing genotype. of the gene classes that are expressed at a low level in the selfers are under-represented but of the over-represented go categories were associated with the regulation of cell development and structure. selfing flowers display modified stamen positions and floral instability including twisted, fused or missing perianth parts, whereas outcrossing plants rarely display these floral modifications. it is possible that the changes in gene expression we have documented influence the regulation of cell growth and division and are responsible for changes to floral morphology that characterize selfing populations. if so, these regulatory loci could be used as a set of candidate genes to investigate aspects of the evolution of the selfing syndrome. differentially expressed floral genes in selfers by sampling different stages of floral development up to and including anthesis we were able to sequence and annotate a large number of florally expressed genes. in total genes with the go annotation'reproduction' were identified, which is a large fraction of the number reported for arabidopsis thaliana in which genes for reproduction have been documented. the lower number of genes in our annotation is not unexpected because we did not include tissue from reproductive stages after flowering, such as fruit and seed development. within the go category'reproduction' we found genes involved in floral development, similar to the number that has been annotated in a. thaliana. of particular significance are the floral development genes that are differentially expressed in the three selfing genotypes, several of which affect structures that are modified in the selfers. anther development and filament elongation are influenced by afb, arf arf, bam, gid1c, myb, pgp and ptr genes in a. thaliana, and pollen maturation and pollen tube growth are altered by afb csdl, cer, pop, myb and tir. erl plays an important roles in normal anther lobe formation and anther cell differentiation and mutants in the erecta gene family have reduced lateral organ size and abnormal flower development, including defects in petal polar expansion and carpel elongation. we also identified a number of genes involved in flowering time including dcl, pgi, shk, suf. lastly, all of the differentially expressed genes that influence ovule development were significantly up-regulated, including arf arf, ova, pep. significantly, most of the candidate genes discussed above cause alterations to attractive structures and male function consistent with the relaxation of selection for showy flowers, reduced allocation to pollen production and the loss of herkogamy through filament elongation of stamens. in contrast, the requirement for functional ovules to maintain seed fertility in selfers may explain the apparent absence of changes to gene expression in female traits.  we have shown that short-read sequencing can be used to characterize the transcriptomes of multiple individuals for use in comparative studies. we were able to assemble as many contigs as other sequencing methods by de novo assembly, but our contigs were on average substantially longer. by comparing sequences from related individuals we generated a consensus transcriptome. this allowed us to make snp genotype calls and provided a method for detecting paralogous sequences. discerning among copies of paralogous sequences presents a major challenge to the current technology and requires either technological or analytical solutions to discern among different members of gene families or duplicates. however, despite these complications we annotated > 80% of contigs and identified thousands of putative orthologs, many of which are differentially expressed. we identified genes that were differentially expressed in selfers and which have developmental functions that suggest a role in the evolution of the selfing syndrome. this sequencing effort has generated a valuable resource of coding dna for a non-model species. the transcriptome sequences will help in future studies of changes in the genetic architecture involved in the transition from outcrossing to selfing and also for identifying the genes controlling heterostyly. authors' contributions the study was conceived by rwn and schb. the molecular lab work was carried out by rwn. sequence assembly and annotation and expression analyses were carried out by rwn. ms contributed programming expertise and to methods for implementing genotyping. all authors read and approved the final manuscript. 
our study is among the first to demonstrate the use of illumina short read sequencing for de novo transcriptome assembly in non-model species, and the first to implement this technology for comparing floral transcriptomes in outcrossing and selfing plants. by sequencing mrna from tissues sampled at various stages of flower development, our goal was to sequence and assemble the floral transcriptome and identify differential patterns of gene expression. increased inbreeding has important consequences for floral morphology, population genetic structure and genome evolution. many of the differentially expressed genes affect floral traits commonly altered in selfing plants and these represent a set of potential candidate genes for investigating the evolution of the selfing syndrome. all four genotypes had highly correlated gene expression, but the three e. paniculata genotypes were more correlated with one another than each was to e. paradoxa. the shift from cross-fertilization to predominant self-fertilization is among the most common evolutionary transitions in the reproductive biology of flowering plants. here, we use short-read sequencing to assemble, de novo, the floral transcriptomes of three genotypes of eichhornia paniculata, including an outcrosser and two genotypes from independently derived selfers, and a single genotype of the sister species e. paradoxa. our analysis identified genes associated with floral development, of which were differentially expressed in selfing lineages relative to the outcrosser.
 the shift from cross-fertilization to predominant self-fertilization is among the most common evolutionary transitions in the reproductive biology of flowering plants. increased inbreeding has important consequences for floral morphology, population genetic structure and genome evolution. the transition to selfing is usually characterized by a marked reduction in flower size and the loss of traits involved in pollinator attraction and the avoidance of self-fertilization. here, we use short-read sequencing to assemble, de novo, the floral transcriptomes of three genotypes of eichhornia paniculata, including an outcrosser and two genotypes from independently derived selfers, and a single genotype of the sister species e. paradoxa. by sequencing mrna from tissues sampled at various stages of flower development, our goal was to sequence and assemble the floral transcriptome and identify differential patterns of gene expression.  our mbp assembly resulted in ~ contigs that averaged ~ bp in length. all four genotypes had highly correlated gene expression, but the three e. paniculata genotypes were more correlated with one another than each was to e. paradoxa. our analysis identified genes associated with floral development, of which were differentially expressed in selfing lineages relative to the outcrosser. many of the differentially expressed genes affect floral traits commonly altered in selfing plants and these represent a set of potential candidate genes for investigating the evolution of the selfing syndrome.  our study is among the first to demonstrate the use of illumina short read sequencing for de novo transcriptome assembly in non-model species, and the first to implement this technology for comparing floral transcriptomes in outcrossing and selfing plants. 