BACKGROUND
white clover  is an allotetraploid forage legume with a genetically determined gametophytic self incompatibility system  <cit> .
this reproductive feature renders it an obligate outbreeder, resulting in high levels of inter- and intragenotypic heterozygosity.
t. occidentale d.e.
coombe, and t. pallescens schreber are the currently extant species that seem to best represent the diploid progenitors of the species  <cit> , although more recent analyses  <cit>  seem to suggest that, while the former is a good candidate, the latter exhibits a greater divergence from the sub-genome to which it has been ascribed.
regardless of the identity of the progenitors, the earliest genetic mapping experiments revealed that the two constituent sub-genomes of white clover exhibit highly conserved macrosynteny, with eight homeolog pairs of linkage groups  <cit> .
this highly conserved type of allopolyploidy, combined with high intragenotypic heterozygosity, is problematic for the development of sequence-based resources for white clover.
this is because the presence of two meiotically independent, but very similar genomes results in a much reduced ability to distinguish between homolog and homeolog copies of genes.
the earliest snp discovery and validation studies in white clover represent a good illustration of the problems associated with distinguishing these two types of variation.
sawbridge et al.
<cit>  generated a collection of over  <dig>  expressed sequence tags  that were clustered to generate a non-redundant set of over  <dig>  consensus sequences.
these were subsequently utilised for snp discovery, exploiting the fact that the underlying cdna libraries were constructed in multiple, heterogeneous genotypes of the white clover variety grasslands huia  <cit> .
white clover varieties are generated via polycrosses of multiple parents and thus are relatively genetically heterogeneous  leading to an expectation of multiple allelic variants for individual loci in the consensus sequences generated by clustering.
from the est data, a total of  <dig>  candidate snps were identified across  <dig> loci  <cit> .
on validating snps from a small subset of these loci in mapping populations, almost half of the snp assays generated monomorphic patterns in f <dig> progeny despite prior validation as polymorphic markers  <cit> .
the authors attributed this to the frequent clustering of sequences from highly similar homeolog copies of genes in the consensus sequences generated in the earlier part of the process, and a subsequent inability to efficiently distinguish between allelic and homeolog specific variation amongst the clustered sequences when designing assays.
basing snp assays on homeolog specific, rather than allelic variants, results in monomorphic assays because alternative homeolog variants are ubiquitously present in segregating populations.
one potential approach for distinguishing homologous, homeologous  sequence variation in est-derived sequences from outbreeders is the generation, cloning and sequencing of pcr amplicons based on the ests.
this approach, which allows unambiguous reconstruction of the component haplotypes at a locus  was examined by hand et al.
<cit> .
for four mapping parents, a sufficient number of pcr amplicons were sequenced to represent all forms of variation associated with over 7 kb of sequence derived from ten stress response genes represented by nine ests.
the same regions were also amplified and sequenced in the proposed diploid progenitors of white clover, t. occidentale and t. pallescens.
comparison of the amplicon sequences generated for each gene in this clone set allowed the identification of all three types of variation in the amplicons .
most interestingly, it was possible to use the proposed ancestor sequences to assign putative homeolog specific variants  into ancestral sub-genomes .
this represents a robust way of both differentiating between allelic and non-allelic variants and also grouping them into consistent sub-genome classes , but the process would be difficult to apply to many thousands of sequences.
homeolocus organisation over tens/hundreds of kilobases has been examined in greater detail for a restricted number of loci in white clover.
hand et al.
<cit>  sequenced eight bacterial artificial chromosome  clones representing four homeologous regions in genic areas of white clover .
a total of 173 kb of overlapping sequence between the o and p’ sub-genomes was generated.
eighteen homeologous genes were identified in the overlapping sub-genome specific sequences and these exhibited conservation of order and orientation.
a further eight genes apparently specific to one or the other sub-genomes were also observed, although this observation is complicated by the fact that homeolog copies of these genes may have been present in the alternative sub-genome just outside the windows of comparison.
levels of nucleotide identity in exonic regions of homeolog pairs ranged from 86% to 100%, with an average of 97%.
in introns, nucleotide identity ranged between 66% and 100% with an average of 89%.
the constantly increasing throughput of next generation sequencing  technology is revolutionising our ability to describe both genome structure and variation.
recently, numerous draft whole genome sequences have been developed in higher plants, based either in part or entirely on ngs-based sequencing strategies.
plant species, many of which have large complex genomes with a history of ancient and recent duplications are still challenging to sequence.
ngs is also rendering the transcriptome more accessible, with technologies such as illumina-based rna-seq allowing sampling to depths capable of detecting relatively low-abundance transcripts in specific tissues and environmental conditions.
this is a useful tool in both functional genomics and the annotation of genome sequences, but also can also provide a potentially stand-alone resource in the face of difficult genomes where whole genome sequencing is problematic.
another potential use of large scale transcript sequencing datasets is in high throughput, genome-wide snp discovery.
this was previously possible using sanger-derived est sequences, but the cost and comparatively low throughput nature of this strategy meant that, in comparison to ngs-based strategies, it was really only effective in relatively well-resourced model and crop plant species, where est data could be aggregated from numerous sources.
the real potential of ngs for snp discovery relies heavily on an ability to both discover and validate snps in an automated, high throughput manner.
while the process of single nucleotide variant detection in allopolyploids is technically no more complex that the same process in homozygous diploids, the presence of multiple independently-segregating but potentially very similar sub-genomes requires extra layers of analysis in order to distinguish allelic and non-allelic variation.
as illustrated by the above studies, this is certainly the case for allotetraploid white clover.
while the approaches described in these studies  <cit>  were both elegant and informative, the processes required to partition genetic variation into allelic and non-allelic types are unsuited to the validation of the thousands of snps required for downstream applications such as association mapping and genomic selection.
recently, trick et al.
<cit>  described a strategy for ngs-based snp detection in allotetraploid brassica napus which allowed the discrimination of allelic single nucleotide variants  and sub-genome-specific or homeolog-specific variants .
the progenitor a and c sub-genomes of b. napus are very similar, with many genes expected to exist as homeolog pairs, the transcript sequences of which are expected to differ by only approximately  <dig> % on average.
as a result, such homeolog pairs will tend to co-assemble or cluster, even when relatively high similarity thresholds  are used.
the strategy is based on the use of doubled haploid  lines, in which the individual sub-genomes have been rendered homozygous.
allelic variants in a sub-genome will be fixed for one or the other alternative nucleotide.
nucleotide variants that distinguish homeolog pairs  originating from different subgenomes will still be present.
when reads from a homeolocus in a single dh line are aligned relative to a reference sequence, nucleotide variants apparent between constituent reads of the alignment are either inter-homeolog variants or real allelic variants that differ between the subgenomes.
in most alignment software, the consensus sequence generated from the alignments will represent such positions with ambiguity codes indicating the presence of more than one nucleotide.
the key point in the strategy is the comparison of the consensus sequences for alignments from different dh genotypes.
thus, any set of genotypes compared in this manner will always have the same consensus nucleotide interhomeolog variant positions.
for positions where real allelic variation is the source of the ambiguity code in one genotype, comparison to a second genotype may reveal that the same position has a normal nucleotide code rather than an ambiguity code.
this is likely to be the footprint of allelic variation at this position.
since, in essence, all that is required is the comparison of the consensus sequences resulting from at least two genotypes for each locus  this approach is amenable to being scaled up for numerous loci.
platforms such as illumina can produce tens or hundred of millions of 100nt transcript reads per genotype, providing a source of genotype-specific sequences to align to a reference sequence, at a depth of coverage that allows the robust detection of nucleotide variants in the alignments.
the levels of nucleic acid identity for the subgenomes of white clover are similar to that of b. napus, so a similar strategy could be implemented with suitable homozygous genotypes.
to our knowledge, no doubled haploid lines of white clover have been developed.
however, inbred lines based on starting material in which the single-locus gametophytic self-incompatibility system has been overcome, presumably by the presence of one or more self-fertility  alleles, have been developed  <cit> .
originally, fourth and fifth generation inbred lines of four distinct origins  were developed by single seed descent.
when these lines were intercrossed they yielded f <dig> hybrid progeny that consistently exhibited heterosis for a range of characteristics, supporting the generally inbred nature of the lines.
two of the lines were subsequently used in the first molecular marker based genetic linkage map of white clover  <cit> , which was based on an f <dig> population of a cross between a fourth generation derivative of one line and a fifth generation derivative of another ).
marker segregation data revealed a high retention of heterozygosity, with up to 60% of simple sequence repeat  markers exhibiting bands suggesting retention of heterozygosity in the parents.
in the study described here, we sought to develop publicly accessible est and snp resources for white clover.
in the first instance we exploited the longer read capabilities  of 454-pyrosequencing to generate an est-like reference transcript set of approximately  <dig>  sequences in a white clover inbred line.
this was then used as a reference for two rounds of snp discovery using shorter  illumina reads, the first of which was based on generating approximately  <dig> million sequences from mixed cdna samples of  <dig> white clover genotypes from a wide geographical range.
in the second round of snp discovery, we used the eighth generation derivatives of two of the sets of white clover inbred lines described above to mimic the genomic constitution that would be expected for doubled haploid lines.
these lines have undergone three to four further generations of inbreeding than the lines described by michaelson-yeates et al.
<cit> , and consequently, we started out from an assumption that the residual heterozygosity reported earlier  <cit>  would be diminished.
trick et al.
<cit>  applied the maq aligner  <cit>  and custom perl scripts to identify and classify snps between dh lines at various read depths.
in this paper we present a similar pipeline based on the combined use of bwa  <cit> , genome analysis toolkit  <cit>  and samtools  <cit>  both to identify single nucleotide variants, and tentatively classify their origin as allelic or homeolog specific.
RESULTS
de novo sequence assembly of  <dig> reads
to provide a reference transcript set for white clover, cdna from mixed above-ground material of a single pot-grown plant of the inbred line s <dig> was subjected to sequencing on the roche gs_flx platform using titanium sequencing chemistry, generating a total of  <dig>  cleaned reads .
assembly of these reads using mira  <cit>  resulted in  <dig>  contigs.
these  <dig>  contigs were built from  <dig>  reads .
of these, a total of  <dig> contigs incorporated more than  <dig> reads per contig, and  <dig> of these incorporated more than  <dig>  reads per contig.
a total of  <dig>  contigs were built from two reads only.
on average  <dig>  reads were built into each kb contig length.
the  <dig>  contigs were clustered using the cdhit-est program of the cd-hit package  <cit>  by setting a similarity threshold of  <dig>  and a word size of  <dig>  the clustering resulted in a non-redundant set of  <dig>  contigs.
from these,  <dig> contigs shorter than  <dig> nt were removed.
the remaining contigs were subjected to blast comparison against trifolium subterraneum chloroplast sequence .
two hundred and sixty-one contigs gave hits with more than 90% similarity at an e-value threshold of >10e- <dig>  these contigs were removed and the remaining  <dig>   were used for further analysis.
the length distribution of the assembly following these processing steps is summarised in figure  <dig>
sequences of  <dig> reads that were not incorporated into contigs were recovered by using the information given in the info_debrislist.txt file provided by the mira assembler.
the sequences were cleaned to remove vector and adapter contamination, and low quality and low complexity terminal sequencess  using the seqclean script  <cit> , resulting in a set of  <dig>  cleaned  <dig> reads.
of these sequences,  <dig>  were shorter than  <dig> nt and were removed.
the remaining  <dig>  sequences were clustered by cd-hit at a similarity threshold level  <dig>  resulting in  <dig>  cluster-representative sequences.
these sequences were subjected to a second clustering step along with the  <dig>  contigs from the previous sequence assembly at a  <dig>  similarity threshold level.
the second clustering step retained  <dig>  singlets.
eighty-three chloroplast specific sequences were identified by blast comparison with the t. subterraneum chloroplast sequence and removed.
the remaining set of  <dig>  non-redundant singlet sequences was used for further analysis.
the mean length of the non-redundant singlet sequence set was  <dig>  nt .
hereafter we refer to the combined set of contigs and non-redundant singlets  as the reference transcript set.
characterisation of the reference transcript sequences
in order to assess the overall quality of the  <dig>  contigs and  <dig>  non-redundant singlets, they were compared to the medicago truncatula gene index  transcript database using the tblastx program of the blast package  <cit> .
the top hits of the tblastx results were filtered into  <dig> categories using a perl script:
“class i” hits: e-value < 1e- <dig>  score >200;
“class ii” hits: e-value < 1e- <dig>  score >100;
“class iii” hits: e-value >1e- <dig>
seventy-seven percent of the contig sequences fell into the top two categories  while the remaining 23%  were in class iii.
for the singlet sequences the situation was reversed, with 66% of sequences  falling onto class iii,  <dig> % into class ii and only  <dig> sequences  in class i .
while this may be partially due to a difference in the gene content of the closely related model and forage legumes, it probably also indicates a generally lower level of sequence quality and reliability for the singlet sequences, and the appropriate caution was exercised in their use for the remainder of the study.
the sequences from the reference transcript set  were subjected to functional annotation by comparison to a custom database of all available viridiplantae protein sequences using the blast2go pipeline  <cit> .
using this approach  <dig> %  of the contigs, but only  <dig> %  of the singlets, could be annotated.
only  <dig> %  of the contigs did not yield significant blastx results compared to  <dig> % of the singlets.
the total number of go annotations was higher in the contigs than in the singlets, but the distribution of gene ontology categories  were largely similar in contigs and singlets .
a summary of the go annotation is presented in additional file  <dig>  and a fasta file of all  <dig>  reference transcript sequences is presented in additional file  <dig>
snp discovery in a mixed sample composed of highly divergent genotypes
in order to exploit the reference transcript set to identify snp variation across the white clover genome we used two approaches.
in the first instance, we generated approximately  <dig> million illumina gaii reads from pooled cdna of  <dig> very divergent heterozygous allotetraploid genotypes representing a pan-global geographic distribution of white clover, and thus likely to capture a significant proportion of the overall variation present in the cultivated genepool of the species.
we aligned these sequences to the reference transcript set to identify snp positions.
after cleaning,  <dig> , <dig> illumina reads  were available for alignment to the reference transcript set.
alignment, mapping and variant calling were performed as described in the methods section.
single nucleotide variants were identified by applying a read depth threshold of  <dig> and a base quality threshold of  <dig>  the process resulted in the identification of  <dig>  snps, .
of the  <dig>  snp variants  <dig>  were due to ambiguity codes in the read consensus of aligned illumina reads, the remainder  were departures from the reference.
in all, snps were detected in  <dig>  reference sequences .
the average read depth at the snp positions was  <dig>   in the case of the contigs and  <dig>   in the case of the singlets.
hereafter we refer the mapping assembly representing this dataset as the dl assembly.
snp discovery and classification using white clover inbred lines
in the second instance, we generated approximately  <dig> million reads for each of two further white clover inbred lines on the illumina hiseq  <dig> platform and aligned these to the reference transcript set.
for identification of snps between the two white clover inbred lines, cdna from individual plants was used.
one of the lines, sc was nearly identical to the s <dig> line used for the reference transcript assembly, as it was sourced from the same batch of seeds produced by a single self-fertilisation event.
line j <dig> was derived from a different, genetically distinct background, with the expectation of reasonable levels of nucleotide polymorphism relative to sc confirmed by some preliminary analysis using two aflp assays .
illumina reads  were aligned to the reference transcript sequences using the pipeline described more fully in the methods section.
from a merged alignment that contained illumina reads from both genotypes  we separately called snps for each read group in positions where there was read depth of at least  <dig>  applying a base quality threshold of  <dig>  this produced two lists of variants , one for each genotype.
conceptually, these variant lists contain positions where variants  occur in both genotypes and positions where variants are observed in one genotype but not the other.
we compared the two variant lists, and for the latter class we examined the merged alignment in regard to the opposing, non-variant genotype, and accepted the non-variant status where there was a read depth of at least  <dig> in the non-variant read group.
this produced a “high confidence” set of variant calls over the two genotypes that were supported by a minimum read depth of  <dig> in each genotype.
this variant set comprised instances in which variants were observed in either both genotypes, or exclusively in one genotype or the other.
in addition, for completeness, we also recorded all instances in which there was a variant in one genotype and not the other, where the second genotype had a read depth of lower than  <dig>  including those for which there was no read coverage .
finally, we recorded instances in there was no apparent variation in either sc or j <dig>  but where there was an apparent departure from the reference sequence.
this analysis identified  <dig>  independent snps in alignments to  <dig>  sequences from the reference transcript set sequences .
no variants satisfying our filtering parameters were found in the remaining reference transcript set of  <dig>  contig sequences and  <dig>  singlet sequences.
the positions and nucleotide composition of the variants identified in both the dl and il assemblies are listed in additional file  <dig>  variant positions are also embedded into the headers of the appropriate sequences in the fasta file presented in additional file  <dig>
in the case of the reference contigs, the mean read depth at the snp positions was  <dig>   for the j <dig> read type and  <dig>   for the sc read type .
perhaps surprisingly, given our assumption that they might have an inherently lower quality, similar short read coverage could be detected for the reference singlets, where the average read depth was  <dig>  for the j <dig> read type  and  <dig>  for the sc read type .
subsequent to variant position identification, we employed the set of principles outlined by trick et al.
<cit>  to tentatively partition snps into classes attributable in origin to homeolog-specific or allelic variation.
using custom perl scripts, the nucleotide consensus of both read type at each variant position from the variant lists for the two inbred lines was compared regarding ambiguity type and relation to the reference base, and variants were assigned to different snp categories  as follows:
inter-homeolog snps
representative of homeolog specific variation.
these comprise variants with the same ambiguity consensus sequence in both the j <dig> and sc read types .
hemi snps
representative of allelic variation.
in these cases, no variant can be detected in one of the read types , but the opposing read type exhibits variation  at a given position .
a subcategory of this class occurs for variants where the two read types have different ambiguous consensus sequences .
we categorized these as asymmetric hemi-snps .
simple snps
this class of variants represents inter-genotype polymorphism where the two read types display different unambiguous consensus nucleotides .
this class could arise from several different sources, eg., effective diploid status at that locus or co-positional allelic variation in both homeolog genomes.
irregular snps
in addition to the three main snp classes, further minor variant classes were identified that could not easily be assigned to the above classes, but may represent real polymorphism, and thus we recorded these in a single category.
this included:  instances in which an ambiguity code was observed in one read type, but low  read depth occurred in the opposing read type;  a class of variant position where both read types have the same unambiguous consensus sequence which differs from the reference sequence - in most cases the reference base was ambiguous at such position e.g.
m/c/c;  positions where one or both read types exhibit consensus of more than two constituent bases .
from the total of  <dig>  variants identified,  <dig>   were categorised as hemi-snps, representative of allelic variation.
a further  <dig>   were categorised as inter-homeolog snps and  <dig>   were classified as simple snps.
the remaining  <dig>  variants  fell into the irregular snp category.
nearly three quarters of the latter  were of the variety in which there was an ambiguity code representing a variant in one genotype, but low or no coverage in the other, making it impossible to assign to a definitive category using our threshold criteria.
the majority of the remainder involved positions where j <dig> and sc exhibited the same nucleotide but differed from the reference base.
these results are summarised in table  <dig>  detailed information regarding position on the reference sequence, snp type and read type specific consensus nucleotide for each snp is given in additional file  <dig>
one interesting feature of this study is that the reference sequence was generated from the same inbred lineage as one of the test genotypes .
as such, there are additional expectations in terms of the types of snp configurations that should be observable in the partitioning analysis relative to a situation where the reference data were aggregated from numerous different sources.
the most obvious expectation relates to the simple snps, in that the test genotype sc should share the same nucleotide variant as the reference sequence in this situation .
in fact, this expectation was met in the vast majority  of the  <dig>  simple snp calls .
a second, slightly more complex set of expectations relates to the hemi-snp calls.
part of the complexity in this case comes from the fact that we have allowed ambiguity codes to be used in the reference transcript assembly.
because of this, variant positions will sometimes be represented as ambiguity codes in the reference contigs, but this will not always be the case, due to the lower sequence depth of the  <dig> sequence dataset, and in these cases, one variant or the other will be represented.
this is not an issue for the singleton sequences, where such ambiguity codes cannot occur.
theoretically, ambiguity codes appearing in the reference contigs should be matched ambiguity codes in sc at that position, since the two genotypes are near-identical.
a more common occurrence is where an ambiguity code is observed in sc but not in the reference sequence.
in this case the nucleotide from the reference sequence should represent one of the two constituent bases of the ambiguity code.
when these factors are taken into account, of the  <dig>  hemi-snps identified,  <dig>   fell into four non-problematic subcategories  that accorded to expectations of identity between the reference sequence and sc, while the remaining  <dig>   fell into categories that were in some way problematic, due to either apparent discord between the reference nucleotide and that in sc, the occurrence of three nucleotide variants rather than two, or some combination thereof.
the latter categories may represent a mixture of sequencing, alignment and variant calling error, deviations from homozygosity, triallelism and paralog incorporation.
*the column contains one example from all of the possible similar configurations.
bold setting indicates non-problematic subcategories as described in the main text.
code h-ax indicates no variant nucleotides in sc, variant in j5;
code h-bx indicates variant nucleotides in sc, no variant in j5;
code ah-x indicates asymmetric variation in both sc and j <dig>
similarly, for the ih snps, there is an expectation of congruence between all three genotypes .
deviations from this expectation occurred in only  <dig>  of  <dig>  ih snp calls.
while, given our model assumptions, the origin of hemi-and ih snps is relatively clear in terms of allelic vs. homeolog specific variation, the origin of simple snps is perhaps less clear.
instances in which simple snps are arising from loci for which no clear homeolog copy exists , or from the transcriptional suppression of one homeolog copy in favour of another are likely to be represented in our dataset by reference sequences in which only simple snps between the inbred lines are present.
from  <dig> reference transcript sequences containing simple snps,  <dig> contained only simple snps , meaning that up to one quarter of contigs containing simple snps could come from the aforementioned sources.
in the majority of occasions, simple snps occur in contigs also containing hemi- and ih-snps, indicating that both homeologs are present and expressed in the genotypes studied.
hand et al.
<cit>  found that approx 22% of snp positions that they surveyed for ten independent genes between a number of mapping population parents contained polymorphisms in both sub-genomes, generally exhibiting the same nucleotide variants .
for many of the positions appearing as simple snps in our data, it’s likely that, in the process of homogenisation associated with inbreeding, opposing allelic variants were fixed in the inbred lines  mapping assembly using  <dig> randomly chosen sequences from the reference transcript set .
our hypothesis was that homozygous loci would consistently exhibit no more than two detectable haplotypes per genotype in the alignments, corresponding to the ancestral sub-genomes, while loci exhibiting retention of heterozygosity would tend to exhibit three to four haplotypes due to combined allelic and inter-genomic variation.
we developed a rigorous set of rules for this visual examination that took into account the fact that the alignments consisted of sequences of only  <dig> nt, limiting the ability to identify “in-phase” variants to instances in which these occur on the same constituent read of an alignment.
in the majority of cases, a nearly equal distribution of two haplotypes could be detected, of which the sc-derived haplotype was similar to that of the reference sequence .
after inspecting the mapping assemblies of the  <dig> contigs, more than two haplotypes were detected in only  <dig>  of the examined contigs, with three or more haplotypes occurring in both genotypes on  <dig> occasions, in j <dig> alone on  <dig> occasions, and in sc alone on  <dig> occasions.
although we hesitate to use this approach to assign a specific estimate for the level of retained heterozygosity on a genome wide level, it supports the hypothesis of an allotetraploid state, consisting largely of two extensively homozygous homeologues for both  test genotypes.
in addition to demonstrating that residual heterozygosity is limited in the inbred lines, this analysis suggests that, in general, the expectation that homeolog copies will tend co-assemble in the reference transcript assembly and subsequent mapping assemblies is valid.
available data suggest that the average nucleotide identity levels in white clover are as high as 97%  <cit> .
with iterative clustering steps at a 95% similarity level during the de novo assembly phase, our process should generally achieve co-assembly of homeolog copies, although nucleotide similarity levels as low as 86% were previously reported  <cit> , and such homeolog pairs may not co-assemble at this similarity level.
less stringent clustering parameters would remediate this effect, but increase the frequency of the incorporation of paralog copies.
while not an in depth examination of paralog co-assembly, the examination of the  <dig> contigs suggests co-assembly of multiple paralogous copies was also not widespread, since this would also yield more than two apparent haplotypes.
comparison of variants from divergent and inbred lines
in order to identify the total number of independent snp positions identified in this study, we made a direct comparison between the variant lists arising from the divergent  and inbred  line mapping assemblies.
comparing the  <dig>  snps from the inbred lines to the  <dig>  snps from the divergent lines revealed that  <dig>  snps were positionally common between the two datasets, with  <dig>  variants occurring only in the inbred line dataset and  <dig>  variants specific to the divergent lines .
thus, in the course of this study a total of  <dig>  independent variant positions were identified in  <dig>  reference transcript sequences  comprising a total sequence length of  <dig> .57 kb .
this corresponds to an average density of  <dig>  snps per kb, or one snp per  <dig>  bp of sequence.
the density distribution of snps on the reference sequences is shown in figure  <dig>
examining the dl assembly at the variant positions that were exclusively detected in the il assembly, we concluded that the main factors responsible for a lack of corresponding variants were either total absence of short read coverage or an insufficient read depth in the dl assembly.
however, in  <dig>  cases, positions that were found to be variant in the il assembly exhibited no variation in the dl assembly at read depth levels of greater than  <dig>  snps were exclusively detected in the dl dataset at  <dig>  positions.
of these,  <dig>  were located in non-variant positions in the il dataset supported by a read depth of greater than  <dig>  and the remaining  <dig>  positions exhibited low read depth and/or low base quality values that prevented the calling of variants in the il assembly.
the co-incidence of a large number of variants between the dl and il assemblies offered the opportunity to compare the frequency of occurrence of the major snp categories from the il analysis in the dl dataset.
variant positions categorised in the il assembly as ih snps are more representative of ancestral variation and thus should have representation approaching 100% in comparable positions in the dl assembly.
conversely, hemi-snps  are representative of allelic variation, and thus one might expect that there are allelic variants represented in the il assembly that are not present in the variation captured in the dl assembly.
in order to test this, we identified all variant positions in the il assembly supported by read depths of  <dig>   <dig> and  <dig> in each individual inbred line, for which there was matching read depth support in the dl assembly.
we then identified whether variants were also present at these positions within the dl assembly.
at the highest stringency level, which required read depth of at least ten in each inbred line and the dl assembly, there were  <dig>  ih snps from the il assembly for which there was sequence coverage in the dl mapping assembly.
of these,  <dig>  positions  were also variant in the dl assembly .
by comparison, of the  <dig>  hemi-snp positions from the il assembly at this read depth level,  <dig>   were also variant in the dl assembly.
thus, at these very well supported positions, ~90% of ih snps were represented in the dl assembly as opposed to ~70% of hemi-snps.
at lower read depths , representation of individual classes was slightly lower, but the difference of ~20% between the categories was maintained.
simple snps were consistently represented at lower levels still .
these significant differences in the representation of categories derived from allelic and homeolog specific variation are in accordance with the expectations described above, and also offer support for the ability of the analysis in the inbred lines to successfully partition these two types of variation.
ih: inter-homeolog snps, he: hemi-snps, s: simple snps.
*number of reference sequence positions where short read coverages were at or above the applied read depth threshold  for both read types in the il assembly, as well as in the dl assembly.
**number and percentage  of dl variants with ambiguity consensus code in matching positions.
note that no base quality filtering was applied for this comparison when searching for matching dl variants.
one potential problem in using transcript data for snp discovery in allopolyploids is the possibility of transcriptional suppression of genes from different sub-genomes in different lines.
on an individual variant basis this could result in the mis-identification of interhomeolog variation as allelic variation.
interestingly, the pipeline described in this study can address this problem at the discovery stage, since the identification of ih snps in alignments to any reference sequence indicates the expression of both homeolog copies of that sequence in both inbred lines.
in addition, the identification of he snps from both inbred lines for a particular reference sequence also indicates the expression of both homeologs.
in this study, of the  <dig>  sequences originally identified as containing variants using the il assembly, 74%  showed clear expression of both  homeolog copies in both of the inbred lines using these criteria.
amongst the remaining 26% of reference sequences, many variants falling below the read depth criteria  were identified, indicating probable expression of both homeologs.
excluding sequences containing irregular category variants for this reason, only 11% of reference sequences  clearly exhibit the footprint of possible homeolog-specific transcriptional suppression in one inbred line.
in these cases, no ih snps were identified and all he and s variants were derived from one line, indicating possible suppression in the other.
this suggests that homeolog-specific transcript suppression is not a prevalent feature of this dataset, and that the analysis pipeline employed can identify potential cases, reducing the need for subsequent validation to account for this phenomenon.
potential for haplotype reconstruction and progenitor comparison
as previously outlined, hand et al.
<cit>  identified and sequenced four pairs of partially overlapping bacs representing the o and p’ subgenomes of white clover.
in a total overlap of 173 kb, eighteen homeologous pairs of genes were identified.
at the time of writing, these represented the best set of matched homeolog gene sequences available from white clover, and we decided to compare our results to them, with the specific objective of assessing the potential for our pipeline to successfully partition sequences into their respective sub-genomes.
we performed blast-based comparisons of the published sequences of these eight bac clones to the contigs from the reference transcript set, and identified nine sequences from the latter which gave good full-length hits to predicted coding sequences on the bacs, allowing them to be identified as unequivocal homologs/homeologs.
we then examined the reference transcript sequences for snps identified using the il mapping assembly.
for three of these sequences, no snps had been identified.
for one sequence we found only two snps, in the irregular category, indicating that they were supported by insufficient coverage.
the remaining five transcript sequences contained snps that had been identified at sufficient coverage to be categorised as hemi-, ih- or simple-snps.
one of these contained only  <dig> simple snps.
the remaining four transcript sequences contained ih- and hemi-snps as well as the other categories .
these alignments present the first opportunity to examine the sub-genomic origin of ih snps identified by our analysis.
of the  <dig> ih snps in the  <dig> transcript sequences,  <dig> had equivalent sequence coverage in the bac-derived gene sequences.
nucleotide variation between the bac-derived sequences from alternate sub-genomes matched ih variation in the inbreds on  <dig> occasions.
at the remaining  <dig> comparable ih snp locations, no variation was observed between bac-derived sequences.
this figure is not unreasonable given that hand et al.
<cit>  found that there was a high degree of coincidence between allelic snp and homeolog specific variation, and a general concordance of nucleotide variants at sites where this was occurring.
they quantified the proportions of co-locating allelic snps and hsvs, and found that, in their dataset, on average, 28% of hsv locations experienced simultaneous allelic snp variation.
this rose to 75% in one individual case.
the effect of this is that, rather than there being a clearly identifiable set of homeolog specific variants that reflect the progenitor state of each sub-genome, homeolog specific variation in white clover is much more variable.
this conclusion has specific ramifications for our attempts to partition the observed snp variation in our inbred lines into allelic and homeolog specific variants, and the general utility of this information.
the most salient feature from a utility point of view is that our analysis should not be viewed as a comprehensive catalogue of homeolog specific variation for the sequences surveyed.
the process we implement will specifically capture instances of interhomeolog variation between these two lines.
given the limitations inherent in exploring the sequence variability at isolated nucleotide positions, we decided to extend the analysis to encompass multiple snp positions simultaneously by attempting to reconstruct the underlying haplotypes for the inbred lines across two of the gene sequences.
we chose the transcript sequences that, given the snp density and distribution , and the read coverage in j <dig> and sc, it was most likely that we would be able to manually reconstruct the haplotypes for the full length of the underlying sequences.
reference transcript sequences wcd- <dig> and wcd- <dig> both possess  <dig> snps in the il assembly .
for both sequences this allowed the reconstruction of two uninterrupted, contiguous haplotypes from each inbred line, presumably corresponding to the two sub-genomes, spanning just over 900 bp in wcd- <dig> and  <dig>  kb in wcd- <dig>  these were subsequently aligned to the bac-derived o and p’ genome specific coding sequences of genes encoding dehydration responsive element binding protein  and acyl-coa oxidase respectively .
interestingly, for the wcd-2781/acyl coa oxidase comparison, the initial alignment shows that, in addition to snp variation, a large indel event spanning ~180 bp exists between the o and p’ sub-genome copies, and this is present in the haplotypes reconstructed from the inbred lines, allowing immediate identification of the sub-genomic origin of the four reconstructed haplotype sequences.
for both this and the wcd-12699/dreb <dig> comparison, sequences were edited to remove non-overlapping regions and realigned to generate simple unrooted phenetic trees to test whether the reconstructed haplotypes would cluster with the different sub-genome specific sequences from the bacs.
in both cases the resulting trees exhibited two distinct clusters , each with one j <dig> and one sc haplotype clustering with one of the bac-derived sequences representing the o and p’ subgenomes.
although limited to two sequences, this demonstrates the potential for both haplotype reconstruction and progenitor-based comparison to further elucidate homeolog relationships based on the type of data generated during this study.
the read lengths of ngs technologies are generally increasing, and even at the time of writing, increased single end read lengths, and the ability for paired end sequencing could make haplotype reconstruction more feasible, especially given the reduction in complexity offered by inbreeding-derived homogeneity.
hand et al.
<cit>  pointed to the need to identify a core set of snps in white clover that are specific to individual sub-genomes and do not coincide with homeolog specific variants.
they suggest a strategy of “extensive sequence-based sampling of germplasm collections using highly parallel technology” combined with the elegant progenitor-based comparison method they describe, largely driven by t. occidentale, since t. pallescens seems to be less similar to the p’ sub-genome the results above suggest a complementary strategy.
were it possible to assemble a reasonable-sized collection of divergent inbred lines , they could be subjected to the pipeline described in this study.
inter-homeolog snp calls would in this case be based on only those ambiguity codes occurring in all of the genotypes under examination, eliminating homeolog specific variant sites at which allelic variation was also occurring, thus producing a core set of ubiquitously present homeolog specific variants.
the success of the strategy would largely be determined by the number of lines examined, but, even using advanced generation inbred lines derived only from the four genotypes described originally  <cit> , it is probable that a significant number of sites exhibiting allelic variation could be discounted.
obviously, such a process cannot partition the variants into specific sub-genomes, but combining it with the progenitor comparison approach and more automated methods of haplotype reconstruction from mapping assemblies would significantly extend the utility of the approach in this respect.
CONCLUSIONS
the main goal of this study was to develop a significant publicly available snp resource for white clover.
we have developed, and partially annotated a reference transcript set of over  <dig>  sequences for white clover, and used these to identify over  <dig>  independent snps in approximately 45% of the reference transcript sequences.
in recognition that the utility of these snps is partially based on whether they are truly allelic or sub-genomic in origin, we have attempted to partition a significant portion of them into categories related to these origins, and examined the most significant parameters likely to effect the accuracy of these categorisations.
the snp data processing and analysis pipeline presented in this study can easily be adapted to snp discovery projects in other allopolyploid species.
based on the principles outlined by trick et al.
<cit>  it differs in that the perl-script enhanced comparison of genotype specific variant lists is obtained and further processed by different utilities of the samtools package rather than the somewhat older maq, overcoming several problems related to implementing this type of analysis with the latter.
samtools takes alignments in the generic sam  format, which is natively supported by various short read aligners.
unlike maq, samtools is under active development, is more likely to be able to handle newer file formats from different sequencing platforms, and produces output formats that are adopted by a wide range of downstream application software.
an additional advantage of using samtools is that short reads from different genotypes can easily be mapped and aligned onto the same set of reference sequences, while variant calling can be performed both globally as well as specifically for any particular genotype, or for a set of genotypes from combined alignments.
this can be effectively utilised in the separation of allelic and sub-genomic variants in two or more genotypes, and we have demonstrated how this type of analysis could be extended to identify a core set of homeolog specific snp variants that would be useful for distinguishing the sub-genomic origin of white clover transcripts, while simultaneously providing an even more comprehensive collection of potential allelic variants for genome-based analysis in white clover.

 white clover is an allotetraploid forage legume with a genetically determined gametophytic self incompatibility system. this reproductive feature renders it an obligate outbreeder, resulting in high levels of inter- and intragenotypic heterozygosity. t. occidentale d.e. coombe, and t. pallescens schreber are the currently extant species that seem to best represent the diploid progenitors of the species, although more recent analyses seem to suggest that, while the former is a good candidate, the latter exhibits a greater divergence from the sub-genome to which it has been ascribed. regardless of the identity of the progenitors, the earliest genetic mapping experiments revealed that the two constituent sub-genomes of white clover exhibit highly conserved macrosynteny, with eight homeolog pairs of linkage groups. this highly conserved type of allopolyploidy, combined with high intragenotypic heterozygosity, is problematic for the development of sequence-based resources for white clover. this is because the presence of two meiotically independent, but very similar genomes results in a much reduced ability to distinguish between homolog and homeolog copies of genes. the earliest snp discovery and validation studies in white clover represent a good illustration of the problems associated with distinguishing these two types of variation. sawbridge et al. generated a collection of over expressed sequence tags that were clustered to generate a non-redundant set of over consensus sequences. these were subsequently utilised for snp discovery, exploiting the fact that the underlying cdna libraries were constructed in multiple, heterogeneous genotypes of the white clover variety grasslands huia. white clover varieties are generated via polycrosses of multiple parents and thus are relatively genetically heterogeneous leading to an expectation of multiple allelic variants for individual loci in the consensus sequences generated by clustering. from the est data, a total of candidate snps were identified across loci. on validating snps from a small subset of these loci in mapping populations, almost half of the snp assays generated monomorphic patterns in f progeny despite prior validation as polymorphic markers. the authors attributed this to the frequent clustering of sequences from highly similar homeolog copies of genes in the consensus sequences generated in the earlier part of the process, and a subsequent inability to efficiently distinguish between allelic and homeolog specific variation amongst the clustered sequences when designing assays. basing snp assays on homeolog specific, rather than allelic variants, results in monomorphic assays because alternative homeolog variants are ubiquitously present in segregating populations. one potential approach for distinguishing homologous, homeologous sequence variation in est-derived sequences from outbreeders is the generation, cloning and sequencing of pcr amplicons based on the ests. this approach, which allows unambiguous reconstruction of the component haplotypes at a locus was examined by hand et al.. for four mapping parents, a sufficient number of pcr amplicons were sequenced to represent all forms of variation associated with over 7 kb of sequence derived from ten stress response genes represented by nine ests. the same regions were also amplified and sequenced in the proposed diploid progenitors of white clover, t. occidentale and t. pallescens. comparison of the amplicon sequences generated for each gene in this clone set allowed the identification of all three types of variation in the amplicons. most interestingly, it was possible to use the proposed ancestor sequences to assign putative homeolog specific variants into ancestral sub-genomes. this represents a robust way of both differentiating between allelic and non-allelic variants and also grouping them into consistent sub-genome classes, but the process would be difficult to apply to many thousands of sequences. homeolocus organisation over tens/hundreds of kilobases has been examined in greater detail for a restricted number of loci in white clover. hand et al. sequenced eight bacterial artificial chromosome clones representing four homeologous regions in genic areas of white clover. a total of 173 kb of overlapping sequence between the o and p’ sub-genomes was generated. eighteen homeologous genes were identified in the overlapping sub-genome specific sequences and these exhibited conservation of order and orientation. a further eight genes apparently specific to one or the other sub-genomes were also observed, although this observation is complicated by the fact that homeolog copies of these genes may have been present in the alternative sub-genome just outside the windows of comparison. levels of nucleotide identity in exonic regions of homeolog pairs ranged from 86% to 100%, with an average of 97%. in introns, nucleotide identity ranged between 66% and 100% with an average of 89%. the constantly increasing throughput of next generation sequencing technology is revolutionising our ability to describe both genome structure and variation. recently, numerous draft whole genome sequences have been developed in higher plants, based either in part or entirely on ngs-based sequencing strategies. plant species, many of which have large complex genomes with a history of ancient and recent duplications are still challenging to sequence. ngs is also rendering the transcriptome more accessible, with technologies such as illumina-based rna-seq allowing sampling to depths capable of detecting relatively low-abundance transcripts in specific tissues and environmental conditions. this is a useful tool in both functional genomics and the annotation of genome sequences, but also can also provide a potentially stand-alone resource in the face of difficult genomes where whole genome sequencing is problematic. another potential use of large scale transcript sequencing datasets is in high throughput, genome-wide snp discovery. this was previously possible using sanger-derived est sequences, but the cost and comparatively low throughput nature of this strategy meant that, in comparison to ngs-based strategies, it was really only effective in relatively well-resourced model and crop plant species, where est data could be aggregated from numerous sources. the real potential of ngs for snp discovery relies heavily on an ability to both discover and validate snps in an automated, high throughput manner. while the process of single nucleotide variant detection in allopolyploids is technically no more complex that the same process in homozygous diploids, the presence of multiple independently-segregating but potentially very similar sub-genomes requires extra layers of analysis in order to distinguish allelic and non-allelic variation. as illustrated by the above studies, this is certainly the case for allotetraploid white clover. while the approaches described in these studies were both elegant and informative, the processes required to partition genetic variation into allelic and non-allelic types are unsuited to the validation of the thousands of snps required for downstream applications such as association mapping and genomic selection. recently, trick et al. described a strategy for ngs-based snp detection in allotetraploid brassica napus which allowed the discrimination of allelic single nucleotide variants and sub-genome-specific or homeolog-specific variants. the progenitor a and c sub-genomes of b. napus are very similar, with many genes expected to exist as homeolog pairs, the transcript sequences of which are expected to differ by only approximately % on average. as a result, such homeolog pairs will tend to co-assemble or cluster, even when relatively high similarity thresholds are used. the strategy is based on the use of doubled haploid lines, in which the individual sub-genomes have been rendered homozygous. allelic variants in a sub-genome will be fixed for one or the other alternative nucleotide. nucleotide variants that distinguish homeolog pairs originating from different subgenomes will still be present. when reads from a homeolocus in a single dh line are aligned relative to a reference sequence, nucleotide variants apparent between constituent reads of the alignment are either inter-homeolog variants or real allelic variants that differ between the subgenomes. in most alignment software, the consensus sequence generated from the alignments will represent such positions with ambiguity codes indicating the presence of more than one nucleotide. the key point in the strategy is the comparison of the consensus sequences for alignments from different dh genotypes. thus, any set of genotypes compared in this manner will always have the same consensus nucleotide interhomeolog variant positions. for positions where real allelic variation is the source of the ambiguity code in one genotype, comparison to a second genotype may reveal that the same position has a normal nucleotide code rather than an ambiguity code. this is likely to be the footprint of allelic variation at this position. since, in essence, all that is required is the comparison of the consensus sequences resulting from at least two genotypes for each locus this approach is amenable to being scaled up for numerous loci. platforms such as illumina can produce tens or hundred of millions of 100nt transcript reads per genotype, providing a source of genotype-specific sequences to align to a reference sequence, at a depth of coverage that allows the robust detection of nucleotide variants in the alignments. the levels of nucleic acid identity for the subgenomes of white clover are similar to that of b. napus, so a similar strategy could be implemented with suitable homozygous genotypes. to our knowledge, no doubled haploid lines of white clover have been developed. however, inbred lines based on starting material in which the single-locus gametophytic self-incompatibility system has been overcome, presumably by the presence of one or more self-fertility alleles, have been developed. originally, fourth and fifth generation inbred lines of four distinct origins were developed by single seed descent. when these lines were intercrossed they yielded f hybrid progeny that consistently exhibited heterosis for a range of characteristics, supporting the generally inbred nature of the lines. two of the lines were subsequently used in the first molecular marker based genetic linkage map of white clover, which was based on an f population of a cross between a fourth generation derivative of one line and a fifth generation derivative of another ). marker segregation data revealed a high retention of heterozygosity, with up to 60% of simple sequence repeat markers exhibiting bands suggesting retention of heterozygosity in the parents. in the study described here, we sought to develop publicly accessible est and snp resources for white clover. in the first instance we exploited the longer read capabilities of 454-pyrosequencing to generate an est-like reference transcript set of approximately sequences in a white clover inbred line. this was then used as a reference for two rounds of snp discovery using shorter illumina reads, the first of which was based on generating approximately million sequences from mixed cdna samples of white clover genotypes from a wide geographical range. in the second round of snp discovery, we used the eighth generation derivatives of two of the sets of white clover inbred lines described above to mimic the genomic constitution that would be expected for doubled haploid lines. these lines have undergone three to four further generations of inbreeding than the lines described by michaelson-yeates et al., and consequently, we started out from an assumption that the residual heterozygosity reported earlier would be diminished. trick et al. applied the maq aligner and custom perl scripts to identify and classify snps between dh lines at various read depths. in this paper we present a similar pipeline based on the combined use of bwa, genome analysis toolkit and samtools both to identify single nucleotide variants, and tentatively classify their origin as allelic or homeolog specific.  de novo sequence assembly of reads to provide a reference transcript set for white clover, cdna from mixed above-ground material of a single pot-grown plant of the inbred line s was subjected to sequencing on the roche gs_flx platform using titanium sequencing chemistry, generating a total of cleaned reads. assembly of these reads using mira resulted in contigs. these contigs were built from reads. of these, a total of contigs incorporated more than reads per contig, and of these incorporated more than reads per contig. a total of contigs were built from two reads only. on average reads were built into each kb contig length. the contigs were clustered using the cdhit-est program of the cd-hit package by setting a similarity threshold of and a word size of the clustering resulted in a non-redundant set of contigs. from these, contigs shorter than nt were removed. the remaining contigs were subjected to blast comparison against trifolium subterraneum chloroplast sequence. two hundred and sixty-one contigs gave hits with more than 90% similarity at an e-value threshold of >10e- these contigs were removed and the remaining were used for further analysis. the length distribution of the assembly following these processing steps is summarised in figure  sequences of reads that were not incorporated into contigs were recovered by using the information given in the info_debrislist.txt file provided by the mira assembler. the sequences were cleaned to remove vector and adapter contamination, and low quality and low complexity terminal sequencess using the seqclean script, resulting in a set of cleaned reads. of these sequences, were shorter than nt and were removed. the remaining sequences were clustered by cd-hit at a similarity threshold level resulting in cluster-representative sequences. these sequences were subjected to a second clustering step along with the contigs from the previous sequence assembly at a similarity threshold level. the second clustering step retained singlets. eighty-three chloroplast specific sequences were identified by blast comparison with the t. subterraneum chloroplast sequence and removed. the remaining set of non-redundant singlet sequences was used for further analysis. the mean length of the non-redundant singlet sequence set was nt. hereafter we refer to the combined set of contigs and non-redundant singlets as the reference transcript set. characterisation of the reference transcript sequences in order to assess the overall quality of the contigs and non-redundant singlets, they were compared to the medicago truncatula gene index transcript database using the tblastx program of the blast package. the top hits of the tblastx results were filtered into categories using a perl script: “class i” hits: e-value < 1e- score >200; “class ii” hits: e-value < 1e- score >100; “class iii” hits: e-value >1e- seventy-seven percent of the contig sequences fell into the top two categories while the remaining 23% were in class iii. for the singlet sequences the situation was reversed, with 66% of sequences falling onto class iii, % into class ii and only sequences in class i. while this may be partially due to a difference in the gene content of the closely related model and forage legumes, it probably also indicates a generally lower level of sequence quality and reliability for the singlet sequences, and the appropriate caution was exercised in their use for the remainder of the study. the sequences from the reference transcript set were subjected to functional annotation by comparison to a custom database of all available viridiplantae protein sequences using the blast2go pipeline. using this approach % of the contigs, but only % of the singlets, could be annotated. only % of the contigs did not yield significant blastx results compared to % of the singlets. the total number of go annotations was higher in the contigs than in the singlets, but the distribution of gene ontology categories were largely similar in contigs and singlets. a summary of the go annotation is presented in additional file and a fasta file of all reference transcript sequences is presented in additional file snp discovery in a mixed sample composed of highly divergent genotypes in order to exploit the reference transcript set to identify snp variation across the white clover genome we used two approaches. in the first instance, we generated approximately million illumina gaii reads from pooled cdna of very divergent heterozygous allotetraploid genotypes representing a pan-global geographic distribution of white clover, and thus likely to capture a significant proportion of the overall variation present in the cultivated genepool of the species. we aligned these sequences to the reference transcript set to identify snp positions. after cleaning,, illumina reads were available for alignment to the reference transcript set. alignment, mapping and variant calling were performed as described in the methods section. single nucleotide variants were identified by applying a read depth threshold of and a base quality threshold of the process resulted in the identification of snps,. of the snp variants were due to ambiguity codes in the read consensus of aligned illumina reads, the remainder were departures from the reference. in all, snps were detected in reference sequences. the average read depth at the snp positions was in the case of the contigs and in the case of the singlets. hereafter we refer the mapping assembly representing this dataset as the dl assembly. snp discovery and classification using white clover inbred lines in the second instance, we generated approximately million reads for each of two further white clover inbred lines on the illumina hiseq platform and aligned these to the reference transcript set. for identification of snps between the two white clover inbred lines, cdna from individual plants was used. one of the lines, sc was nearly identical to the s line used for the reference transcript assembly, as it was sourced from the same batch of seeds produced by a single self-fertilisation event. line j was derived from a different, genetically distinct background, with the expectation of reasonable levels of nucleotide polymorphism relative to sc confirmed by some preliminary analysis using two aflp assays. illumina reads were aligned to the reference transcript sequences using the pipeline described more fully in the methods section. from a merged alignment that contained illumina reads from both genotypes we separately called snps for each read group in positions where there was read depth of at least applying a base quality threshold of this produced two lists of variants, one for each genotype. conceptually, these variant lists contain positions where variants occur in both genotypes and positions where variants are observed in one genotype but not the other. we compared the two variant lists, and for the latter class we examined the merged alignment in regard to the opposing, non-variant genotype, and accepted the non-variant status where there was a read depth of at least in the non-variant read group. this produced a “high confidence” set of variant calls over the two genotypes that were supported by a minimum read depth of in each genotype. this variant set comprised instances in which variants were observed in either both genotypes, or exclusively in one genotype or the other. in addition, for completeness, we also recorded all instances in which there was a variant in one genotype and not the other, where the second genotype had a read depth of lower than including those for which there was no read coverage. finally, we recorded instances in there was no apparent variation in either sc or j but where there was an apparent departure from the reference sequence. this analysis identified independent snps in alignments to sequences from the reference transcript set sequences. no variants satisfying our filtering parameters were found in the remaining reference transcript set of contig sequences and singlet sequences. the positions and nucleotide composition of the variants identified in both the dl and il assemblies are listed in additional file variant positions are also embedded into the headers of the appropriate sequences in the fasta file presented in additional file in the case of the reference contigs, the mean read depth at the snp positions was for the j read type and for the sc read type. perhaps surprisingly, given our assumption that they might have an inherently lower quality, similar short read coverage could be detected for the reference singlets, where the average read depth was for the j read type and for the sc read type. subsequent to variant position identification, we employed the set of principles outlined by trick et al. to tentatively partition snps into classes attributable in origin to homeolog-specific or allelic variation. using custom perl scripts, the nucleotide consensus of both read type at each variant position from the variant lists for the two inbred lines was compared regarding ambiguity type and relation to the reference base, and variants were assigned to different snp categories as follows: inter-homeolog snps representative of homeolog specific variation. these comprise variants with the same ambiguity consensus sequence in both the j and sc read types. hemi snps representative of allelic variation. in these cases, no variant can be detected in one of the read types, but the opposing read type exhibits variation at a given position. a subcategory of this class occurs for variants where the two read types have different ambiguous consensus sequences. we categorized these as asymmetric hemi-snps. simple snps this class of variants represents inter-genotype polymorphism where the two read types display different unambiguous consensus nucleotides. this class could arise from several different sources, eg., effective diploid status at that locus or co-positional allelic variation in both homeolog genomes. irregular snps in addition to the three main snp classes, further minor variant classes were identified that could not easily be assigned to the above classes, but may represent real polymorphism, and thus we recorded these in a single category. this included: instances in which an ambiguity code was observed in one read type, but low read depth occurred in the opposing read type; a class of variant position where both read types have the same unambiguous consensus sequence which differs from the reference sequence - in most cases the reference base was ambiguous at such position e.g. m/c/c; positions where one or both read types exhibit consensus of more than two constituent bases. from the total of variants identified, were categorised as hemi-snps, representative of allelic variation. a further were categorised as inter-homeolog snps and were classified as simple snps. the remaining variants fell into the irregular snp category. nearly three quarters of the latter were of the variety in which there was an ambiguity code representing a variant in one genotype, but low or no coverage in the other, making it impossible to assign to a definitive category using our threshold criteria. the majority of the remainder involved positions where j and sc exhibited the same nucleotide but differed from the reference base. these results are summarised in table  detailed information regarding position on the reference sequence, snp type and read type specific consensus nucleotide for each snp is given in additional file one interesting feature of this study is that the reference sequence was generated from the same inbred lineage as one of the test genotypes. as such, there are additional expectations in terms of the types of snp configurations that should be observable in the partitioning analysis relative to a situation where the reference data were aggregated from numerous different sources. the most obvious expectation relates to the simple snps, in that the test genotype sc should share the same nucleotide variant as the reference sequence in this situation. in fact, this expectation was met in the vast majority of the simple snp calls. a second, slightly more complex set of expectations relates to the hemi-snp calls. part of the complexity in this case comes from the fact that we have allowed ambiguity codes to be used in the reference transcript assembly. because of this, variant positions will sometimes be represented as ambiguity codes in the reference contigs, but this will not always be the case, due to the lower sequence depth of the sequence dataset, and in these cases, one variant or the other will be represented. this is not an issue for the singleton sequences, where such ambiguity codes cannot occur. theoretically, ambiguity codes appearing in the reference contigs should be matched ambiguity codes in sc at that position, since the two genotypes are near-identical. a more common occurrence is where an ambiguity code is observed in sc but not in the reference sequence. in this case the nucleotide from the reference sequence should represent one of the two constituent bases of the ambiguity code. when these factors are taken into account, of the hemi-snps identified, fell into four non-problematic subcategories that accorded to expectations of identity between the reference sequence and sc, while the remaining fell into categories that were in some way problematic, due to either apparent discord between the reference nucleotide and that in sc, the occurrence of three nucleotide variants rather than two, or some combination thereof. the latter categories may represent a mixture of sequencing, alignment and variant calling error, deviations from homozygosity, triallelism and paralog incorporation. *the column contains one example from all of the possible similar configurations. bold setting indicates non-problematic subcategories as described in the main text. code h-ax indicates no variant nucleotides in sc, variant in j5; code h-bx indicates variant nucleotides in sc, no variant in j5; code ah-x indicates asymmetric variation in both sc and j similarly, for the ih snps, there is an expectation of congruence between all three genotypes. deviations from this expectation occurred in only of ih snp calls. while, given our model assumptions, the origin of hemi-and ih snps is relatively clear in terms of allelic vs. homeolog specific variation, the origin of simple snps is perhaps less clear. instances in which simple snps are arising from loci for which no clear homeolog copy exists, or from the transcriptional suppression of one homeolog copy in favour of another are likely to be represented in our dataset by reference sequences in which only simple snps between the inbred lines are present. from reference transcript sequences containing simple snps, contained only simple snps, meaning that up to one quarter of contigs containing simple snps could come from the aforementioned sources. in the majority of occasions, simple snps occur in contigs also containing hemi- and ih-snps, indicating that both homeologs are present and expressed in the genotypes studied. hand et al. found that approx 22% of snp positions that they surveyed for ten independent genes between a number of mapping population parents contained polymorphisms in both sub-genomes, generally exhibiting the same nucleotide variants. for many of the positions appearing as simple snps in our data, it’s likely that, in the process of homogenisation associated with inbreeding, opposing allelic variants were fixed in the inbred lines mapping assembly using randomly chosen sequences from the reference transcript set. our hypothesis was that homozygous loci would consistently exhibit no more than two detectable haplotypes per genotype in the alignments, corresponding to the ancestral sub-genomes, while loci exhibiting retention of heterozygosity would tend to exhibit three to four haplotypes due to combined allelic and inter-genomic variation. we developed a rigorous set of rules for this visual examination that took into account the fact that the alignments consisted of sequences of only nt, limiting the ability to identify “in-phase” variants to instances in which these occur on the same constituent read of an alignment. in the majority of cases, a nearly equal distribution of two haplotypes could be detected, of which the sc-derived haplotype was similar to that of the reference sequence. after inspecting the mapping assemblies of the contigs, more than two haplotypes were detected in only of the examined contigs, with three or more haplotypes occurring in both genotypes on occasions, in j alone on occasions, and in sc alone on occasions. although we hesitate to use this approach to assign a specific estimate for the level of retained heterozygosity on a genome wide level, it supports the hypothesis of an allotetraploid state, consisting largely of two extensively homozygous homeologues for both test genotypes. in addition to demonstrating that residual heterozygosity is limited in the inbred lines, this analysis suggests that, in general, the expectation that homeolog copies will tend co-assemble in the reference transcript assembly and subsequent mapping assemblies is valid. available data suggest that the average nucleotide identity levels in white clover are as high as 97%. with iterative clustering steps at a 95% similarity level during the de novo assembly phase, our process should generally achieve co-assembly of homeolog copies, although nucleotide similarity levels as low as 86% were previously reported, and such homeolog pairs may not co-assemble at this similarity level. less stringent clustering parameters would remediate this effect, but increase the frequency of the incorporation of paralog copies. while not an in depth examination of paralog co-assembly, the examination of the contigs suggests co-assembly of multiple paralogous copies was also not widespread, since this would also yield more than two apparent haplotypes. comparison of variants from divergent and inbred lines in order to identify the total number of independent snp positions identified in this study, we made a direct comparison between the variant lists arising from the divergent and inbred line mapping assemblies. comparing the snps from the inbred lines to the snps from the divergent lines revealed that snps were positionally common between the two datasets, with variants occurring only in the inbred line dataset and variants specific to the divergent lines. thus, in the course of this study a total of independent variant positions were identified in reference transcript sequences comprising a total sequence length of.57 kb. this corresponds to an average density of snps per kb, or one snp per  bp of sequence. the density distribution of snps on the reference sequences is shown in figure  examining the dl assembly at the variant positions that were exclusively detected in the il assembly, we concluded that the main factors responsible for a lack of corresponding variants were either total absence of short read coverage or an insufficient read depth in the dl assembly. however, in cases, positions that were found to be variant in the il assembly exhibited no variation in the dl assembly at read depth levels of greater than snps were exclusively detected in the dl dataset at positions. of these, were located in non-variant positions in the il dataset supported by a read depth of greater than and the remaining positions exhibited low read depth and/or low base quality values that prevented the calling of variants in the il assembly. the co-incidence of a large number of variants between the dl and il assemblies offered the opportunity to compare the frequency of occurrence of the major snp categories from the il analysis in the dl dataset. variant positions categorised in the il assembly as ih snps are more representative of ancestral variation and thus should have representation approaching 100% in comparable positions in the dl assembly. conversely, hemi-snps are representative of allelic variation, and thus one might expect that there are allelic variants represented in the il assembly that are not present in the variation captured in the dl assembly. in order to test this, we identified all variant positions in the il assembly supported by read depths of and in each individual inbred line, for which there was matching read depth support in the dl assembly. we then identified whether variants were also present at these positions within the dl assembly. at the highest stringency level, which required read depth of at least ten in each inbred line and the dl assembly, there were ih snps from the il assembly for which there was sequence coverage in the dl mapping assembly. of these, positions were also variant in the dl assembly. by comparison, of the hemi-snp positions from the il assembly at this read depth level, were also variant in the dl assembly. thus, at these very well supported positions, ~90% of ih snps were represented in the dl assembly as opposed to ~70% of hemi-snps. at lower read depths, representation of individual classes was slightly lower, but the difference of ~20% between the categories was maintained. simple snps were consistently represented at lower levels still. these significant differences in the representation of categories derived from allelic and homeolog specific variation are in accordance with the expectations described above, and also offer support for the ability of the analysis in the inbred lines to successfully partition these two types of variation. ih: inter-homeolog snps, he: hemi-snps, s: simple snps. *number of reference sequence positions where short read coverages were at or above the applied read depth threshold for both read types in the il assembly, as well as in the dl assembly. **number and percentage of dl variants with ambiguity consensus code in matching positions. note that no base quality filtering was applied for this comparison when searching for matching dl variants. one potential problem in using transcript data for snp discovery in allopolyploids is the possibility of transcriptional suppression of genes from different sub-genomes in different lines. on an individual variant basis this could result in the mis-identification of interhomeolog variation as allelic variation. interestingly, the pipeline described in this study can address this problem at the discovery stage, since the identification of ih snps in alignments to any reference sequence indicates the expression of both homeolog copies of that sequence in both inbred lines. in addition, the identification of he snps from both inbred lines for a particular reference sequence also indicates the expression of both homeologs. in this study, of the sequences originally identified as containing variants using the il assembly, 74% showed clear expression of both homeolog copies in both of the inbred lines using these criteria. amongst the remaining 26% of reference sequences, many variants falling below the read depth criteria were identified, indicating probable expression of both homeologs. excluding sequences containing irregular category variants for this reason, only 11% of reference sequences clearly exhibit the footprint of possible homeolog-specific transcriptional suppression in one inbred line. in these cases, no ih snps were identified and all he and s variants were derived from one line, indicating possible suppression in the other. this suggests that homeolog-specific transcript suppression is not a prevalent feature of this dataset, and that the analysis pipeline employed can identify potential cases, reducing the need for subsequent validation to account for this phenomenon. potential for haplotype reconstruction and progenitor comparison as previously outlined, hand et al. identified and sequenced four pairs of partially overlapping bacs representing the o and p’ subgenomes of white clover. in a total overlap of 173 kb, eighteen homeologous pairs of genes were identified. at the time of writing, these represented the best set of matched homeolog gene sequences available from white clover, and we decided to compare our results to them, with the specific objective of assessing the potential for our pipeline to successfully partition sequences into their respective sub-genomes. we performed blast-based comparisons of the published sequences of these eight bac clones to the contigs from the reference transcript set, and identified nine sequences from the latter which gave good full-length hits to predicted coding sequences on the bacs, allowing them to be identified as unequivocal homologs/homeologs. we then examined the reference transcript sequences for snps identified using the il mapping assembly. for three of these sequences, no snps had been identified. for one sequence we found only two snps, in the irregular category, indicating that they were supported by insufficient coverage. the remaining five transcript sequences contained snps that had been identified at sufficient coverage to be categorised as hemi-, ih- or simple-snps. one of these contained only simple snps. the remaining four transcript sequences contained ih- and hemi-snps as well as the other categories. these alignments present the first opportunity to examine the sub-genomic origin of ih snps identified by our analysis. of the ih snps in the transcript sequences, had equivalent sequence coverage in the bac-derived gene sequences. nucleotide variation between the bac-derived sequences from alternate sub-genomes matched ih variation in the inbreds on occasions. at the remaining comparable ih snp locations, no variation was observed between bac-derived sequences. this figure is not unreasonable given that hand et al. found that there was a high degree of coincidence between allelic snp and homeolog specific variation, and a general concordance of nucleotide variants at sites where this was occurring. they quantified the proportions of co-locating allelic snps and hsvs, and found that, in their dataset, on average, 28% of hsv locations experienced simultaneous allelic snp variation. this rose to 75% in one individual case. the effect of this is that, rather than there being a clearly identifiable set of homeolog specific variants that reflect the progenitor state of each sub-genome, homeolog specific variation in white clover is much more variable. this conclusion has specific ramifications for our attempts to partition the observed snp variation in our inbred lines into allelic and homeolog specific variants, and the general utility of this information. the most salient feature from a utility point of view is that our analysis should not be viewed as a comprehensive catalogue of homeolog specific variation for the sequences surveyed. the process we implement will specifically capture instances of interhomeolog variation between these two lines. given the limitations inherent in exploring the sequence variability at isolated nucleotide positions, we decided to extend the analysis to encompass multiple snp positions simultaneously by attempting to reconstruct the underlying haplotypes for the inbred lines across two of the gene sequences. we chose the transcript sequences that, given the snp density and distribution, and the read coverage in j and sc, it was most likely that we would be able to manually reconstruct the haplotypes for the full length of the underlying sequences. reference transcript sequences wcd- and wcd- both possess snps in the il assembly. for both sequences this allowed the reconstruction of two uninterrupted, contiguous haplotypes from each inbred line, presumably corresponding to the two sub-genomes, spanning just over 900 bp in wcd- and  kb in wcd- these were subsequently aligned to the bac-derived o and p’ genome specific coding sequences of genes encoding dehydration responsive element binding protein and acyl-coa oxidase respectively. interestingly, for the wcd-2781/acyl coa oxidase comparison, the initial alignment shows that, in addition to snp variation, a large indel event spanning ~180 bp exists between the o and p’ sub-genome copies, and this is present in the haplotypes reconstructed from the inbred lines, allowing immediate identification of the sub-genomic origin of the four reconstructed haplotype sequences. for both this and the wcd-12699/dreb comparison, sequences were edited to remove non-overlapping regions and realigned to generate simple unrooted phenetic trees to test whether the reconstructed haplotypes would cluster with the different sub-genome specific sequences from the bacs. in both cases the resulting trees exhibited two distinct clusters, each with one j and one sc haplotype clustering with one of the bac-derived sequences representing the o and p’ subgenomes. although limited to two sequences, this demonstrates the potential for both haplotype reconstruction and progenitor-based comparison to further elucidate homeolog relationships based on the type of data generated during this study. the read lengths of ngs technologies are generally increasing, and even at the time of writing, increased single end read lengths, and the ability for paired end sequencing could make haplotype reconstruction more feasible, especially given the reduction in complexity offered by inbreeding-derived homogeneity. hand et al. pointed to the need to identify a core set of snps in white clover that are specific to individual sub-genomes and do not coincide with homeolog specific variants. they suggest a strategy of “extensive sequence-based sampling of germplasm collections using highly parallel technology” combined with the elegant progenitor-based comparison method they describe, largely driven by t. occidentale, since t. pallescens seems to be less similar to the p’ sub-genome the results above suggest a complementary strategy. were it possible to assemble a reasonable-sized collection of divergent inbred lines, they could be subjected to the pipeline described in this study. inter-homeolog snp calls would in this case be based on only those ambiguity codes occurring in all of the genotypes under examination, eliminating homeolog specific variant sites at which allelic variation was also occurring, thus producing a core set of ubiquitously present homeolog specific variants. the success of the strategy would largely be determined by the number of lines examined, but, even using advanced generation inbred lines derived only from the four genotypes described originally, it is probable that a significant number of sites exhibiting allelic variation could be discounted. obviously, such a process cannot partition the variants into specific sub-genomes, but combining it with the progenitor comparison approach and more automated methods of haplotype reconstruction from mapping assemblies would significantly extend the utility of the approach in this respect.  the main goal of this study was to develop a significant publicly available snp resource for white clover. we have developed, and partially annotated a reference transcript set of over sequences for white clover, and used these to identify over independent snps in approximately 45% of the reference transcript sequences. in recognition that the utility of these snps is partially based on whether they are truly allelic or sub-genomic in origin, we have attempted to partition a significant portion of them into categories related to these origins, and examined the most significant parameters likely to effect the accuracy of these categorisations. the snp data processing and analysis pipeline presented in this study can easily be adapted to snp discovery projects in other allopolyploid species. based on the principles outlined by trick et al. it differs in that the perl-script enhanced comparison of genotype specific variant lists is obtained and further processed by different utilities of the samtools package rather than the somewhat older maq, overcoming several problems related to implementing this type of analysis with the latter. samtools takes alignments in the generic sam format, which is natively supported by various short read aligners. unlike maq, samtools is under active development, is more likely to be able to handle newer file formats from different sequencing platforms, and produces output formats that are adopted by a wide range of downstream application software. an additional advantage of using samtools is that short reads from different genotypes can easily be mapped and aligned onto the same set of reference sequences, while variant calling can be performed both globally as well as specifically for any particular genotype, or for a set of genotypes from combined alignments. this can be effectively utilised in the separation of allelic and sub-genomic variants in two or more genotypes, and we have demonstrated how this type of analysis could be extended to identify a core set of homeolog specific snp variants that would be useful for distinguishing the sub-genomic origin of white clover transcripts, while simultaneously providing an even more comprehensive collection of potential allelic variants for genome-based analysis in white clover. 
the potential for using haplotype reconstruction and progenitor genome comparison to assign haplotypes to specific ancestral sub-genomes of white clover is demonstrated for sequences corresponding to genes encoding dehydration responsive element binding protein and acyl-coa oxidase. allelic variantsallotetraploidhaplotype reconstructionnext generation sequencingsingle nucleotide polymorphismssub-genometranscript assemblytrifoliumwhite clover the apparent existence of highly similar homeolog copies for the majority of genes in white clover is problematic for the development of genome-based resources in the species. robust methods for categorising single nucleotide variants as allelic or homeolog-specific in large transcript datasets are required. this is especially true for the development of genetic markers based on single nucleotide polymorphisms, since it is difficult to distinguish between homeolog-specific and allelic variants. we discuss the potential to extend the analysis to identify a “core set” of ancestrally derived homeolog specific variants in white clover. we illustrate one potential approach in this study. in total, independent snps in reference sequences were discovered, approximately three quarters of which were categorised as representing allelic or homeolog-specific variation using two inbred lines. white clover is an allotetraploid species possessing two highly collinear ancestral sub-genomes. we used 454-pyrosequencing sequencing to generate ~ transcript sequences from an 8th generation white clover inbred line. we subsequently performed illumina sequencing on three further white clover samples, generating million transcript reads from a mixed sample comprising divergent white clover genotypes, and million reads on two further eighth generation white clover inbred lines. mapping these reads to the reference transcript set allowed us to develop a significant snp resource for white clover, and to partition the snps from the inbred lines into categories reflecting allelic or homeolog-specific variation. this represents a significant resource for white clover genomics and genetics studies. these were assembled and partially annotated to yield a reference transcript set comprising sequences.
 white clover is an allotetraploid species possessing two highly collinear ancestral sub-genomes. the apparent existence of highly similar homeolog copies for the majority of genes in white clover is problematic for the development of genome-based resources in the species. this is especially true for the development of genetic markers based on single nucleotide polymorphisms, since it is difficult to distinguish between homeolog-specific and allelic variants. robust methods for categorising single nucleotide variants as allelic or homeolog-specific in large transcript datasets are required. we illustrate one potential approach in this study.  we used 454-pyrosequencing sequencing to generate ~ transcript sequences from an 8th generation white clover inbred line. these were assembled and partially annotated to yield a reference transcript set comprising sequences. we subsequently performed illumina sequencing on three further white clover samples, generating million transcript reads from a mixed sample comprising divergent white clover genotypes, and million reads on two further eighth generation white clover inbred lines. mapping these reads to the reference transcript set allowed us to develop a significant snp resource for white clover, and to partition the snps from the inbred lines into categories reflecting allelic or homeolog-specific variation. the potential for using haplotype reconstruction and progenitor genome comparison to assign haplotypes to specific ancestral sub-genomes of white clover is demonstrated for sequences corresponding to genes encoding dehydration responsive element binding protein and acyl-coa oxidase.  in total, independent snps in reference sequences were discovered, approximately three quarters of which were categorised as representing allelic or homeolog-specific variation using two inbred lines. this represents a significant resource for white clover genomics and genetics studies. we discuss the potential to extend the analysis to identify a “core set” of ancestrally derived homeolog specific variants in white clover. allelic variantsallotetraploidhaplotype reconstructionnext generation sequencingsingle nucleotide polymorphismssub-genometranscript assemblytrifoliumwhite clover 