BACKGROUND
until recently, genome wide sequencing has been cost limited to relatively few species.
newer sequencing technologies, offering rapid high throughput sequencing, reduction of sequencing cost, and complementary software have changed the situation dramatically, making these methods more readily considered as standard molecular tools for biology and medicine  <cit> .
transcriptome sequencing in particular can now be undertaken by individual research groups to generate genomic information for novel or "orphan" species of scientific interest  <cit> .
one such species is the red fox , including its farm bred coat color variant, the silver fox.
the red fox is one of the most widely distributed mammalian species across the globe  <cit>  - inhabiting the northern hemisphere naturally and australia by introduction.
foxes live in a wide range of temperature zones, geographical regions, wild and urban environments  <cit>  and are close relatives of the domestic dog .
the grey wolf  and the fox diverged from a common ancestor about  <dig> million years ago  <cit> , and the dog was domesticated from the grey wolf at least  <dig>  years ago  <cit> .
in contrast, the red fox -- or, more precisely, the silver fox -- has only recently been domesticated as part of a deliberate scientific experiment at the institute of cytology and genetics  of the russian academy of sciences in the second part of the twentieth century .
starting from a farm-bred population of silver foxes, the icg selectively bred foxes for tameness for over  <dig> generations  <cit> .
foxes from this tame strain exhibit human-directed behavior resembling that of domestic dogs http://cbsu.tc.cornell.edu/ccgr/behaviour/index.htm <cit> .
in parallel, the icg also undertook a selective breeding program to develop another strain of foxes that are aggressive and difficult to handle.
both strains of foxes are maintained at the icg as outbred populations.
the genetic basis of these tame and aggressive behavioral phenotypes was confirmed in experiments that included cross-fostering of pups, transplantation of embryos and cross-breeding experiments  <cit> .
these fox strains present a robust model, intermediate between rodent and primate in biological complexity, for the study of genetics of affiliative versus aggressive behaviors.
by exploiting genomic resources developed for the domestic dog - the only canid species with a sequenced genome, to date - we recently began molecular genetic studies of these selectively bred tame and aggressive fox strains.
to this end, we developed crossbred f <dig>  back-cross and intercross pedigrees segregating behavioral differences derived from the founder strains; a meiotic linkage map of the fox genome based on  <dig> canine-derived microsatellite markers  <cit> ; and, using these resources, mapped several quantitative trait loci  in the fox genome associated with behavior  <cit> .
although genetic mapping with canine derived microsatellite markers allowed identification of several behavioral loci in fox genome, this type of genetic analysis has inherent limitations in its mapping power.
single nucleotide polymorphism  markers would be much more powerful for fine mapping these behavioral loci, and identification of haplotypes cosegregating with qtls that were under selection in the founder strains.
unfortunately, the majority of dog snps are not informative in foxes and fox-specific snps need to be identified to move this study forward.
to this end we undertook deep transcriptome sequencing of pre-frontal brain samples from one tame and one aggressive fox.
because pre-frontal cortex is a brain area relatively easy to identify, samples from different individuals can be readily collected in a relatively consistent manner.
furthermore, the role of this brain region in aspects of behavior is well established  <cit> , and thus particularly relevant to the study.
comparison of gene expression profiles between a tame and an aggressive individual offers insights into functional differences in the pre-frontal brain cortex.
comparison of transcriptome sequences of the same genes between the tame and aggressive fox samples, has identified a large set of informative snp markers and begun a catalogue of gene-specific sequence variants between the two strains.
RESULTS
<dig> gs flx titanium sequencing the transcriptome of the silver fox prefrontal cortex
sequencing of non-normalized cdna libraries from prefrontal cortex of one tame and one aggressive male fox after image analysis, signal processing and data filtering yielded  <dig> , <dig> reads  from the tame fox sample and  <dig> , <dig> reads  from the aggressive fox sample .
distribution of filtered read lengths in each sample is presented in figure  <dig>  from the tame sample, the average read length was  <dig> bp , and the median  <dig> bp.
from the aggressive, the average was  <dig> bp , and the median  <dig> bp.
the primary sequencing data have been deposited in the genebank sequence read archive  under accession number .
the number of reads is tabulated for both samples , both individually and all together, with the total number of bases for each set that passed filtering criteria and were used for mapping against the canine refseq database, together with analytic results of the corresponding fox transcriptome assemblies.
isogroups: an isogroup is a collection of contigs containing reads with implied connections between them, analogous to individual genes; isotigs: an isotig is meant to be analogous to an individual transcript or splice variant; contigs: a contig represents the consensus sequence of a self-consistent set of overlapping reads, roughly analogous to exons; singletons: reads that did not overlap with any other reads in the input; repeats: reads that were either inferred to be repetitive early in the assembly process or determined to partially overlap a contig; outliers: reads that were identified by the gs de novo assembler as problematic, and were excluded from the final contigs; too short: the trimmed read was too short  to used in the computation.
.
assembly of fox transcriptome sequencing reads
fox sequencing reads that passed quality filtering were assembled de novo in three different manners: i) all reads together; ii) reads from the tame sample alone; iii) reads from the aggressive sample alone .
over 70% of reads were fully assembled, and close to a further 10% of reads were partly assembled in both the assembly of all reads together and in each individual assembly run .
blat of fox contigs to the dog genome
among all contigs identified in the assembly of all reads together we selected the  <dig>  contigs that were at least  <dig> bp long .
blat of these contigs against the canine genome sequence assembly  successfully located  <dig>  contigs.
the output file of blat results was then visualized on the ucsc genome browser .
alignment of fox contigs against the dog genome allowed quick visual examination of the genes represented in the fox transcriptome, revealed gene structure and possible splice variants.
blastx of fox isotigs against the swissprot data base
to estimate the number of genes well represented in this fox transcriptome, all isotigs identified in the assembly of all fox reads together  were aligned against the swissprot database using blastx.
of these,  <dig>  isotigs  yielded a hit to one or more proteins in the database.
only the best hit for each isotig was selected.
successfully aligned isotigs belonged to  <dig>  isogroups .
in total,  <dig>  genes were identified, of which  <dig>   were in more than one isogroup.
the total number of genes identified in the fox transcriptome was significantly lower than the number of isogroups identified in the assembly of fox reads.
putative homologs of canine genes in the fox transcriptome
because a fox reference genome sequence assembly does not exist, and the dog is the fox's closest relative with a sequenced genome, the canine refseq database was used to map fox reads to specific genes by nucleotide sequence alignment.
all  fox reads that passed filtering for quality control , from the six sequencing runs, were mapped against the canine refseq database .
in total,  <dig> % of these reads mapped either fully  or partially  to a known gene in the canine refseq database , and a total of  <dig>  canine genes  were mapped as homologs of at least one fox read .
after excluding from the remaining reads those that were short, chimeric or mapped repetitively,  <dig> % of the total number of reads did not map to the canine refseq database .
completeness of the fox transcriptome
because non-normalized cdna libraries were sequenced, it was expected that the most abundant transcripts in the fox transcriptome would be represented by a higher number of sequencing reads than would rarer transcripts.
to estimate whether fox cdna libraries were sequenced to a sufficient depth we mapped i) reads from individual  <dig> flx sequencing runs  and ii) reads from several sequencing runs combined together, to make datasets of several different sizes, against the dog refseq database in separate experiments .
mapping all available reads from the tame  and aggressive  samples in separate experiments identified  <dig>  and  <dig>  of dog gene orthologs, respectively.
the number of genes identified in the canine refseq database by mapping  <dig>   <dig>  and  <dig>  sequencing runs  was very similar in each case , suggesting that the data plot  was approaching its saturation point.
the plot likewise showed that increasing the depth of sequencing data from  <dig> to  <dig> half full runs identified only about 5% more genes.
the average breadth of coverage of dog refseq transcripts that were hit by at least one fox read was  <dig> % .
the average depth of coverage for dog refseq sequences that had at least one hit was  <dig>  for tame and  <dig>  for aggressive.
in determining breadth of coverage, we used a conservative approach and counted only reads which were uniquely mapped to a single transcript by the  <dig> gs reference mapper.
because the reads corresponding to multiple splice variants were counted only for a single splice variant , it is likely that the average depth of coverage of dog transcripts was underestimated.
the fact that the breadth of coverage of dog transcripts mapped by reads from two individual samples or both samples was almost identical suggests that an average breadth of coverage at ~66% is the practical maximum we can expect to obtain using our mapping approach.
identification of fox snps
sequence differences between dog and fox and among fox chromosomes were identified by mapping fox reads against the canine genome sequence assembly, and contigs from the assembly of all fox reads together.
•snps identified by mapping fox sequencing reads against the dog genome
to identify nucleotide differences between the dog and fox genomes and between the two foxes sequenced, we mapped all fox reads against the repeat-masked dog genomic sequence  .
in total,  <dig> % of reads mapped fully  or partially  to the dog genome;  <dig> % of fox reads mapped chimerically, with non-overlapping parts of the read aligned to distant locations in the dog genome.
reads that mapped to multiple locations  and reads that were too short to map reliably  were excluded from further analysis .
the mapping run created a total of  <dig>  dog-fox contigs containing  <dig> , <dig> bases.
comparison of sequences from the four fox chromosomes represented in the present study with the one canine reference chromosome  identified a total of  <dig>  sequence differences, including single nucleotide differences, indels and 2- <dig> nucleotide differences.
for  <dig>  of these differences only 1-5% of fox reads were in disagreement with the dog sequence.
these latter differences were deemed less likely to be true nucleotide differences and were excluded from further analysis.
this left  <dig>  differences representing one difference  for every  <dig> bases of dog-fox contigs.
of these  <dig>  differences,  <dig>  were fox-versus-dog differences, where at least 95% of fox reads were in agreement with each other, and differed from the dog sequence.
this represented one difference per  <dig> mapped bases.
an example of one such dog-fox nucleotide difference is presented in supplementary figure 2a .
in total,  <dig>  differences  were identified among fox reads that aligned to the dog genome, an example of which is presented in supplementary figure 2b .
the set of differences where the fox reads differed from each other were filtered to obtain  <dig>  high confidence snps.
• snps identified by mapping fox reads against the assembled fox transcriptome sequence
a total of  <dig>  fox contigs that were at least  <dig> bp long, from the assembly of all reads together , were used as a reference for mapping fox reads.
in total,  <dig>  differences between fox reads were identified.
these polymorphisms included snps, indels, and 2- <dig> nucleotide polymorphisms.
after applying filtering criteria   <dig>  high-confidence fox snps were identified.
contigs containing high-confidence snps were selected and localized in the dog genome by blat, and the output file was visualized on the ucsc genomic browser.
• non-redundant set of high-confidence fox-specific snps
high confidence fox snps were selected using parameters described in the methods section.
snps identified by both methods  were checked for redundancy and duplicates were deleted.
the non-redundant set contained  <dig>  snps.
snps were classified into four categories according to the number of alleles observed in each sample .
snps were classified by the observed genotypes in the tame and aggressive individuals, as either homozygous in both , heterozygous in both, or homozygous in one individual and heterozygous in the other.
in total,  <dig>  of these high confidence fox snps were localized in the dog genome and the inferred location of  <dig>  of these snps on fox chromosomes was assigned using the comparative fox/dog map .
the informativeness of each snp in each sample is also represented on fox chromosome graphs .
figure  <dig> shows the distribution and zygosity of snps in a segment of vvu <dig> that has conserved synteny to dog chromosome  <dig>  this region corresponds to part of a qtl interval identified by genetic mapping of behavioural phenotypes in fox experimental pedigrees  <cit> .
note that the aggressive and tame individuals are homozygous for opposite snp alleles in this interval.
in the tame animal, the region of homozygosity extends even further.
snp validation
to establish parameters for filtering high-confidence fox snps  we selected a subset of snps for validation by sanger sequencing.
three sets of snps were selected.
the first set included  <dig> snps in which the minor allele was present in at least 25% of the sequencing reads per sample and each allele was present in more than three reads in one individual ; the second set included  <dig> snps in which the minor allele was present in at least 25% of sequencing reads per sample and by exactly three reads in one individual ; the third set included  <dig> snps in which the minor allele was present in at least three reads between the two individuals .
snps that mapped to the canine x chromosome and showed heterozygosity in at least one sample were filtered out  and not included in the snp validation study.
of  <dig> snps, all but one  were validated by sanger sequencing .
allele heterozygosity, identified by sanger and  <dig> sequencing for each snp in each individual, differed for two snps in the first dataset , five snps in the second dataset , and two snps in the third data set .
for all these snps except snp #3- <dig> one of the two individuals was homozygous by  <dig> sequencing, but heterozygous by sanger sequencing .
the allele identified by sanger sequencing corresponded to the second allele observed in another individual.
the opposite situation was observed for snp #3- <dig>  with two snp alleles identified in the tame individual by  <dig> sequencing but only one by sanger sequencing .
all snps in the first two data sets were validated by sanger sequencing, indicating that the selection parameters applied identified high-confidence snps reliably.
over 90% of snps in the third data set  were also validated.
sanger resequencing of snps did find, however, that a second snp allele was not always identified in the  <dig> sequencing datasets at the coverage cutoff used for selection of high-confidence snps .
indeed, all the minor snp alleles not identified by  <dig> sequencing were among those with the lowest read coverage depth .
distribution of fox snps among different types of genomic regions
blat of snp-bearing fox amplicons to the dog genome localized a subset of snps to non-coding regions of the dog genome.
to evaluate the distribution of fox snps among known coding and non-coding regions we selected  <dig> random snps that passed defined stringent filtering parameters  and localized these snps in the dog genome  using the ucsc genome browser.
in total,  <dig> snps were identified in exons or utrs,  <dig> in introns or non-coding regions, and  <dig> snps were not localized unambiguously in the dog genome .
these results demonstrate that although the majority of high-confidence snps identified in the fox transcriptome correspond to translated sequences, our reads also contain sequences from non-coding genomic regions or regions not currently known to be transcribed in the dog genome.
<dig> randomly selected high-confidence fox snps were mapped onto the dog genome sequence assembly  to identify their distribution among different types of genome sequence.
comparison of gene expression between the two fox samples
a total of  <dig> genes had a tenfold difference of expression between the tame and aggressive individuals with p <  <dig> .
<dig> genes had at least a twofold difference of expression and p <  <dig>  .
genes with differences in expression between the two individuals were classified into biofunctional groups, using ingenuity ipa http://www.ingenuity.com <cit> .
the  <dig> most statistically significant biofunction groups for both the tame and the aggressive upregulated genes are shown in figure  <dig>  networks were also analyzed, and a network of  <dig> genes was found in the set of upregulated genes from the aggressive individual that contains genes involved in the behavior, nervous system development and function and cardiovascular development and function biofunction groups .
validation of gene expression by rt-qpcr
eleven genes that showed expression differences between the two fox samples, in analysis of transcriptome data, were selected for validation by rt-qpcr.
genes were chosen that had more than one exon and did not belong to a highly conserved gene family.
differences in expression of selected genes varied between two samples in terms of number of reads, fold differences, and statistical significance of observed expression differences.
six genes in this set had higher expression levels  in the transcriptome of the aggressive fox, versus five in the tame individual.
primers were designed to successfully amplify unique products for nine of these genes, but not successfully for pkhd1l <dig> and btnl <dig> - the latter two genes were thus excluded from further analysis.
the relative expression level was estimated for these nine genes  by calculating the difference in expression between "tame" and "aggressive" samples.
three genes  showed less than a two-fold difference between rt-qpcr results and transcriptome analysis.
five genes  showed more than a two-fold lower expression difference in rt-qpcr experiments than in the transcriptome analysis, although the difference between the "tame" and "aggressive" samples was in the same direction as in the transcriptome analysis.
for four of these genes the difference in expression between the "tame" and "aggressive" samples detected by rt-qpcr was less than two-fold.
one gene, slitrk <dig>  showed a difference of expression by rt-qpcr in the opposite direction from that by transcriptome analysis .
expression difference of slitrk <dig> in the transcriptome analysis had less significant support  among all genes analysed in rt-qpcr experiments.
nine genes, identified as differentially expressed between the tame and aggressive fox by transcriptome sequence analysis, were selected for expression evaluation by relative rt-qpcr.
results are sorted by z test p-value.
except for one gene , which had the least significant z test p-value, rt-qpcr consistently confirmed a difference of expression between the tame and aggressive samples in the same direction as had transcriptome analysis .
discussion
this study was designed to expand genomic resources for the red  fox, an emerging animal model for genetic research whose genome has not yet been sequenced.
using the  <dig> gs flx titanium platform we sequenced non-normalized cdna libraries from brain samples of two farm-bred foxes with markedly different genetically determined behavioral phenotypes.
samples for sequencing were selected from prefrontal cortex because this region is easily recognisable in the brain; can thus be collected in a consistent manner; and is related to behavior.
the role of different subdivisions of prefrontal cortex in regulating various aspects of cognition, memory, emotions, and social behavior is reviewed in clark et al.,  <dig>  <cit> .
in humans, a link between brain damage of prefrontal cortex and aggressive behavior has been reported  <cit> .
a total of  <dig>  gb of transcriptome sequencing data for the tame fox and  <dig>  gb for the aggressive fox were obtained.
the depth of sequencing chosen for this study was based on previous experience sequencing non-normalized cdna libraries  <cit> .
the fox reads were assembled in three different combinations: i) only reads from the tame fox, ii) only reads from the aggressive fox; iii) all reads together.
the percent of reads assigned into contigs in each of the three assemblies was very similar .
the assembly software used, the  <dig> gs de novo assembler , takes into account the existence of gene splice variants.
newbler first assembles sets of self-consistent reads into contigs, which are the rough equivalent of exons.
it then builds isotigs  from contigs that are consistently connected by subpopulations of reads.
when there are splice variants, contigs can be connected via different paths that are supported by read data, leading to multiple isotigs being generated.
isotigs are grouped into isogroups, which are the rough equivalent of genes.
in the context of a transcriptome assembly by newbler, contigs are the sections of isogroups where corresponding isotig sequences are the same.
consequently, the assembly produces two sets of sequences: a theoretically non-redundant set of contigs, and a set of isotigs which are candidate long- or full-length transcripts from multiple splice variants of genes, and contains some partially redundant sequence segments.
newbler is currently the only transcripome assembler that attempts this task, and is one of the mostly widely used assemblers for de novo assembly of organisms with few genomic recourses  <cit> .
comparison of transcriptome assemblies from different studies is inherently problematic because different assembly algorithms handle the complexities of the transciptome in a variety of ways.
cap <dig>  which was used in the assembly of transcriptomes of two mammalian species -- the bank vole  and the domestic ferret   <cit>  -- does not split reads into multiple contigs as newbler does  <cit> .
this makes the number of contigs non-comparable.
one measure of comparison is the span of contigs, the total number of bases in the all contigs in an assembly.
this number is available for the bank vole, as the number of contigs reported,  <dig> , multiplied by the average length of the contigs,  <dig>  for an assembly span of  <dig>  mb.
assembly of the fox transcriptome spans  <dig>  mb.
relative to the bank vole, the fox transcriptome was sequenced to approximately  <dig>  times greater depth and the assembly span was higher by  <dig>  times.
the number of isogroups, isotigs, and contigs identified in the assembly of all fox reads together was significantly higher than the number of isogroups, isotigs, and contigs identified in the separate assemblies of tame and aggressive reads .
singleton reads in the individual assemblies, when assembled together, created new contigs, isotigs and isogroups.
furthermore, sequence variations between individuals increased the number of contigs and isotigs in isogroups, when these variations were sufficiently large for the assembly program to assume them to be different exons, or real differences in splice variants of some genes between the two individuals.
the average isotig length and n <dig> isotig length were also longer in the assembly of all reads together, as would be expected from an increase in the amount of sequencing data .
because the fox genome has not yet been sequenced and a fox refseq database is thus not available, we used swissprot and the dog refseq database to identify gene orthologs in the fox transcriptome.
blastx of fox isotigs from the assembly of all reads together against swissprot database and mapping of all fox reads against dog refseq database identified very similar number of genes .
the number of identified genes in the fox transcriptome was significantly lower than the number of isogroups  found in the assembly of fox reads.
the total number of isogroups was also larger than the total number of genes in mammalian genome.
one hypothesis is that the amount of sequencing data used in the assembly was not sufficient to cover the full lengths of all fox gene transcripts, therefore some poorly covered genes might have been partitioned into two or more isogroups.
consistent with this hypothesis, blastx of fox isotigs against swissprot identified  <dig>  genes, of which  <dig>   were found in more than one isogroup .
in addition, a subset of fox reads  was mapped to non-gene regions in the dog genome.
these regions are expected not to be present in the dog refseq database.
we expect that in the fox transcriptome these regions will each be represented by an isogroup with only one isotig assigned to it.
also, not all genes that are transcribed are described in the databases we used.
the genes not in public databases are expected to be genes which would have the more rare, and less covered sequences, and subsequently less complex isogroups.
this low level of expression would imply that they would have been less likely to be found and described previously.
this may explain why only  <dig> %  of all isogroups with a single isotig were successfully aligned using blastx.
at the same time,  <dig> %  of isogroups containing more than one isotig  yielded a hit to the swissprot database.
the completeness of the fox transcriptome was tested by mapping fox reads combined variously into datasets of different sizes against the dog refseq database .
one striking result of this experiment was the identification of a very similar number of dog gene orthologs in transcriptomes from tame and aggressive samples .
the total number of genes identified in both transcriptomes was only approximately 5% larger than a number of genes identified in individual samples .
these data demonstrate that, without increasing the depth of sequencing by at least an order of magnitude, it is unlikely that a significantly higher number of rare transcripts would be identified in these transcriptome libraries.
the number of genes identified, in this study, from fox pre-frontal brain transcriptome is similar to that identified from mouse cortex   <cit> .
for each accession number in the dog refseq database to which at least one fox read mapped, we estimated the breadth of coverage of dog transcripts by fox reads.
figure  <dig> demonstrates the bivariate coverage of identified transcripts: roughly 75% of dog transcripts had less than 90% coverage and ~25% of dog transcripts had coverage over 90%.
we examined the depth of coverage of dog transcripts with relatively low coverage.
the slight rise in the middle of the graph, from around 40-60% of coverage contains, as expected, mainly transcripts with a low depth of coverage, but it also contains a significant portion of the transcripts, responsible for the rise above a simple fit curve, that have a large number of fox reads mapped to them.
it suggests that in many cases it is not the depth of coverage that is limiting the breadth of coverage.
the transcripts in the canine refseq database can originate from different tissues and it is likely that for some genes, brain transcripts only partially overlap with transcripts from other tissues.
possible sequence differences between canine and fox transcripts would also reduce the coverage.
the breadth of coverage of dog transcripts mapped by reads from two individual samples or from both samples together was almost identical.
these data again suggest that a moderate further increase in the amount of sequencing data is unlikely to increase the breadth of coverage of dog transcripts.
the dog is the closest relative of the fox with a sequenced genome  <cit> .
the dog and fox, however, have very different karyotypes, the dog having  <dig> mostly acrocentric chromosomes, whereas the red fox has  <dig> metacentric chromosomes  <cit> .
the cytological relationship between the dog and fox genomes is well understood: cytogenetic methods  <cit>  and alignment of the fox meiotic map against of the dog genome  <cit>  have demonstrated that most fox chromosomes correspond to two or three canine chromosomes.
the canine genome sequence assembly canfam <dig> was used to localize fox sequencing reads and contigs in the dog genome.
the approximate location of corresponding regions in the fox genome could be then inferred using the comparative dog-fox genetic and cytogenetic maps.
mapping fox reads to the canine genome allowed the sequence divergence between the dog and fox genomes to be estimated.
one sequence difference  between fox and dog was identified per  <dig> mapped bases.
the number of differences between fox chromosomes was one difference per  <dig>  mapped bases.
the level of divergence between fox and dog coding regions is comparable with the level of divergence between sheep and cattle  <cit> .
mapping fox reads against both the canine genome and the assembled fox transcriptome sequence identified  <dig>  high-confidence fox-specific snps.
the filtering criteria were established by sanger resequencing of snps with different numbers of reads per allele.
the number of snps identified in this study significantly exceeded those identified in the bank vole transcriptome , although filtering parameters in our study  significantly exceed the filtering parameters used in the vole transcriptome study   <cit> .
these comparisons demonstrate that a greater depth of transcriptome sequencing enables the discovery of a large number of snp's with higher confidence.
our filtering parameters yielded 100% confirmation of the polymorphisms in the validation study  but differences in individual genotype  identified by  <dig> sequencing were different from sanger sequencing in some cases.
to validate snps by sanger sequencing, we used cdna prepared from exactly the same rna samples which were used for  <dig> sequencing .
therefore, we expect allele specific bias to be equal in the samples used in both  <dig> and sanger sequencing experiments.
all snps for which the second allele in the same sample was not identified by  <dig> sequencing were among snps with the lowest read coverage depth.
over 90% of snps from those with the lowest sequence coverage included in the validation study  were also validated by sanger sequencing.
this filter parameter is still stricter than that in the bank vole study  <cit> .
therefore, for a small risk of false positives, we can also use some of these lower confidence snp's to increase the total set of snps by about 50% to around  <dig>  snps.
taking into account the significant cost difference between  <dig> and sanger technologies, our data suggest that, at the small cost of obtaining a relatively small number of false polymorphisms, snps from this category would be worth including in a snp array.
a small number of false snps can be easily detected during post-processing analysis of genotyping calls and filtered out prior to further analysis.
for regions in the fox genome associated with traits of interest it would also be useful to invest in validation of snps represented by an even lower number of reads.
although the false positive rate for these snps is expected to be higher, validation experiments can still provide a significant number of snps for these critical regions.
the  <dig> technology has been successfully used for developing snp resources for multiple organisms with limited genomic resources  <cit> .
initially, a large number of snp's that mapped to the canine x chromosome were heterozygous in at least one male individual.
to test whether this represented errors in our data or in the assembly of the canine x chromosome, the contigs bearing these snps were mapped onto the human x chromosome assembly.
all fox contigs containing snps that mapped to both the human and the dog x chromosome outside of the recognized pseudoautosomal region were homozygous for each individual fox.
to gain a better understanding of the high heterozygosity of snps assigned to cfax, we mapped approximately  <dig> bp region surrounding  <dig> randomly selected snps that were reported as heterozygous on the x back to the dog genome.
we found that, in 34% of cases the area surrounding the snp could not be located unambiguously to the x, or to a single area on the x; 50% of snps had high homology to the x but also had significant homology with another area of the genome.
this casts doubt on the x chromosome snps showing heterozygosity in our male individuals.
we subsequently filtered out all snps on the cfax that were not in the pseudoautosomal region and were heterozygous in either individual.
of all the fox high-confidence snps identified  <dig> % were localised in the dog genome, and the inferred location in the fox genome was identified for  <dig> % .
heterozygosity of snps in the two fox samples was compared and several continuous regions along fox chromosomes where snps were fixed for opposite alleles between tame and aggressive sample were revealed.
these regions warrant investigation in a larger number of samples from tame and aggressive fox strains for evaluation of putative selection sweeps.
one extended region of reduced heterozygosity in the tame sample was observed in a region of vvu <dig> corresponding to part of a previously identified qtl interval  <cit> .
overall, the number of high confidence snps and their locations on fox chromosomes indicated that now we have a set of snps sufficient for genome wide mapping in these populations.
detailed evaluation of the location in the dog genome of a randomly selected set of  <dig> fox snps revealed that approximately 77% of snps are located in transcribed regions and distributed approximately equally between coding regions and utrs.
approximately 9% of snps are located in introns and 14% in non-coding regions.
similar results were also identified in the bank vole transcriptome study  <cit> .
as has become apparent recently, the mammalian transcriptome is more complex than previously recognized  <cit> .
comparison of gene expression between the two fox samples revealed only a relatively small number of genes whose expression level differed by two-fold or more .
the biofunction groups over-represented in the set of differentially expressed genes differed between the two samples.
genes associated with neurological diseases were over represented in the tame sample.
interestingly, the second biofunction group over-represented in the aggressive sample represents genes associated with cardiovascular diseases.
the links between the cardiovascular system and behavior have been increasingly recognized recently  <cit> .
to validate differences in gene expression that were identified by transcriptome sequencing, we selected a subset of nine genes for expression analysis by rt-qpcr.
exactly the same two rna samples  which had been used for transcriptome sequencing were used for rt-qpcr validation experiments.
eight of nine differentially expressed genes selected for the validation study showed differences in expression in the same direction by both methods.
all genes that were validated by rt-qpcr showed differential expression between the two samples in transcriptome analysis with p <  <dig> .
although the consistent validation by rt-qpcr of results obtained by transcriptome analysis increases confidence in the observed differences in expression between the two samples, it will clearly be necessary to test these observed differences in a larger set of samples before coming to conclusions regarding biological relevance.
among the genes that were validated by rt-qpcr, htr2c is particularly noteworthy.
htr2c plays an important role in serotonergic and dopaminergic signalling  <cit> ; and is differentially expressed in specific brain regions of two strains of rats - one of which was selected for and the other against aggressive behavior  <cit> .
increased level of htr2c mrna in the frontal cortex and hippocampus of tame rats was observed  <cit> .
notably, the direction of difference in htr2c expression between the tame and aggressive fox individuals in the present study was consistent with that between rat strains in the popova et al.
study.
none of the genes reported in other studies  <cit>  investigating comparative gene expression in canid brain samples were significantly different between the two fox samples in the present study.
further resequencing and deeper analysis of brain transcriptome, including identification of non-synonymous mutations, comparison of splice variants between tame and aggressive samples, and evaluation of differentially expressed genes in larger independent sample sets  will be necessary to characterize critical differences in these fox strains associated with markedly different behavioral phenotypes.
in any experimental system there will be a certain amount of noise that might stem from a variety of sources.
these biological replicates would let us explore the differences in the underlying biological system we are interested in and validate the differences in gene expression.
CONCLUSIONS
using new sequencing technologies it is now possible to develop genomic resources for species without a sequenced genome.
that dramatically facilitates the study of interesting but previously inaccessible biological phenomena and helps to establish new animal models for genetic research.
overall, transcriptome sequencing from the prefrontal cortex of two fox individuals identified over  <dig>  high-confidence fox-specific snps, fox orthologs of over  <dig>  dog genes, and yielded new insights into potentially important differences in expression of genes in the pre-frontal cortex between tame and aggressive foxes.

 until recently, genome wide sequencing has been cost limited to relatively few species. newer sequencing technologies, offering rapid high throughput sequencing, reduction of sequencing cost, and complementary software have changed the situation dramatically, making these methods more readily considered as standard molecular tools for biology and medicine. transcriptome sequencing in particular can now be undertaken by individual research groups to generate genomic information for novel or "orphan" species of scientific interest. one such species is the red fox, including its farm bred coat color variant, the silver fox. the red fox is one of the most widely distributed mammalian species across the globe - inhabiting the northern hemisphere naturally and australia by introduction. foxes live in a wide range of temperature zones, geographical regions, wild and urban environments and are close relatives of the domestic dog. the grey wolf and the fox diverged from a common ancestor about million years ago, and the dog was domesticated from the grey wolf at least years ago. in contrast, the red fox -- or, more precisely, the silver fox -- has only recently been domesticated as part of a deliberate scientific experiment at the institute of cytology and genetics of the russian academy of sciences in the second part of the twentieth century. starting from a farm-bred population of silver foxes, the icg selectively bred foxes for tameness for over generations. foxes from this tame strain exhibit human-directed behavior resembling that of domestic dogs http://cbsu.tc.cornell.edu/ccgr/behaviour/index.htm. in parallel, the icg also undertook a selective breeding program to develop another strain of foxes that are aggressive and difficult to handle. both strains of foxes are maintained at the icg as outbred populations. the genetic basis of these tame and aggressive behavioral phenotypes was confirmed in experiments that included cross-fostering of pups, transplantation of embryos and cross-breeding experiments. these fox strains present a robust model, intermediate between rodent and primate in biological complexity, for the study of genetics of affiliative versus aggressive behaviors. by exploiting genomic resources developed for the domestic dog - the only canid species with a sequenced genome, to date - we recently began molecular genetic studies of these selectively bred tame and aggressive fox strains. to this end, we developed crossbred f back-cross and intercross pedigrees segregating behavioral differences derived from the founder strains; a meiotic linkage map of the fox genome based on canine-derived microsatellite markers; and, using these resources, mapped several quantitative trait loci in the fox genome associated with behavior. although genetic mapping with canine derived microsatellite markers allowed identification of several behavioral loci in fox genome, this type of genetic analysis has inherent limitations in its mapping power. single nucleotide polymorphism markers would be much more powerful for fine mapping these behavioral loci, and identification of haplotypes cosegregating with qtls that were under selection in the founder strains. unfortunately, the majority of dog snps are not informative in foxes and fox-specific snps need to be identified to move this study forward. to this end we undertook deep transcriptome sequencing of pre-frontal brain samples from one tame and one aggressive fox. because pre-frontal cortex is a brain area relatively easy to identify, samples from different individuals can be readily collected in a relatively consistent manner. furthermore, the role of this brain region in aspects of behavior is well established, and thus particularly relevant to the study. comparison of gene expression profiles between a tame and an aggressive individual offers insights into functional differences in the pre-frontal brain cortex. comparison of transcriptome sequences of the same genes between the tame and aggressive fox samples, has identified a large set of informative snp markers and begun a catalogue of gene-specific sequence variants between the two strains.  gs flx titanium sequencing the transcriptome of the silver fox prefrontal cortex sequencing of non-normalized cdna libraries from prefrontal cortex of one tame and one aggressive male fox after image analysis, signal processing and data filtering yielded, reads from the tame fox sample and, reads from the aggressive fox sample. distribution of filtered read lengths in each sample is presented in figure from the tame sample, the average read length was bp, and the median bp. from the aggressive, the average was bp, and the median bp. the primary sequencing data have been deposited in the genebank sequence read archive under accession number. the number of reads is tabulated for both samples, both individually and all together, with the total number of bases for each set that passed filtering criteria and were used for mapping against the canine refseq database, together with analytic results of the corresponding fox transcriptome assemblies. isogroups: an isogroup is a collection of contigs containing reads with implied connections between them, analogous to individual genes; isotigs: an isotig is meant to be analogous to an individual transcript or splice variant; contigs: a contig represents the consensus sequence of a self-consistent set of overlapping reads, roughly analogous to exons; singletons: reads that did not overlap with any other reads in the input; repeats: reads that were either inferred to be repetitive early in the assembly process or determined to partially overlap a contig; outliers: reads that were identified by the gs de novo assembler as problematic, and were excluded from the final contigs; too short: the trimmed read was too short to used in the computation.. assembly of fox transcriptome sequencing reads fox sequencing reads that passed quality filtering were assembled de novo in three different manners: i) all reads together; ii) reads from the tame sample alone; iii) reads from the aggressive sample alone. over 70% of reads were fully assembled, and close to a further 10% of reads were partly assembled in both the assembly of all reads together and in each individual assembly run. blat of fox contigs to the dog genome among all contigs identified in the assembly of all reads together we selected the contigs that were at least bp long. blat of these contigs against the canine genome sequence assembly successfully located contigs. the output file of blat results was then visualized on the ucsc genome browser. alignment of fox contigs against the dog genome allowed quick visual examination of the genes represented in the fox transcriptome, revealed gene structure and possible splice variants. blastx of fox isotigs against the swissprot data base to estimate the number of genes well represented in this fox transcriptome, all isotigs identified in the assembly of all fox reads together were aligned against the swissprot database using blastx. of these, isotigs yielded a hit to one or more proteins in the database. only the best hit for each isotig was selected. successfully aligned isotigs belonged to isogroups. in total, genes were identified, of which were in more than one isogroup. the total number of genes identified in the fox transcriptome was significantly lower than the number of isogroups identified in the assembly of fox reads. putative homologs of canine genes in the fox transcriptome because a fox reference genome sequence assembly does not exist, and the dog is the fox's closest relative with a sequenced genome, the canine refseq database was used to map fox reads to specific genes by nucleotide sequence alignment. all fox reads that passed filtering for quality control, from the six sequencing runs, were mapped against the canine refseq database. in total, % of these reads mapped either fully or partially to a known gene in the canine refseq database, and a total of canine genes were mapped as homologs of at least one fox read. after excluding from the remaining reads those that were short, chimeric or mapped repetitively, % of the total number of reads did not map to the canine refseq database. completeness of the fox transcriptome because non-normalized cdna libraries were sequenced, it was expected that the most abundant transcripts in the fox transcriptome would be represented by a higher number of sequencing reads than would rarer transcripts. to estimate whether fox cdna libraries were sequenced to a sufficient depth we mapped i) reads from individual flx sequencing runs and ii) reads from several sequencing runs combined together, to make datasets of several different sizes, against the dog refseq database in separate experiments. mapping all available reads from the tame and aggressive samples in separate experiments identified and of dog gene orthologs, respectively. the number of genes identified in the canine refseq database by mapping and sequencing runs was very similar in each case, suggesting that the data plot was approaching its saturation point. the plot likewise showed that increasing the depth of sequencing data from to half full runs identified only about 5% more genes. the average breadth of coverage of dog refseq transcripts that were hit by at least one fox read was %. the average depth of coverage for dog refseq sequences that had at least one hit was for tame and for aggressive. in determining breadth of coverage, we used a conservative approach and counted only reads which were uniquely mapped to a single transcript by the gs reference mapper. because the reads corresponding to multiple splice variants were counted only for a single splice variant, it is likely that the average depth of coverage of dog transcripts was underestimated. the fact that the breadth of coverage of dog transcripts mapped by reads from two individual samples or both samples was almost identical suggests that an average breadth of coverage at ~66% is the practical maximum we can expect to obtain using our mapping approach. identification of fox snps sequence differences between dog and fox and among fox chromosomes were identified by mapping fox reads against the canine genome sequence assembly, and contigs from the assembly of all fox reads together. •snps identified by mapping fox sequencing reads against the dog genome to identify nucleotide differences between the dog and fox genomes and between the two foxes sequenced, we mapped all fox reads against the repeat-masked dog genomic sequence. in total, % of reads mapped fully or partially to the dog genome; % of fox reads mapped chimerically, with non-overlapping parts of the read aligned to distant locations in the dog genome. reads that mapped to multiple locations and reads that were too short to map reliably were excluded from further analysis. the mapping run created a total of dog-fox contigs containing, bases. comparison of sequences from the four fox chromosomes represented in the present study with the one canine reference chromosome identified a total of sequence differences, including single nucleotide differences, indels and 2- nucleotide differences. for of these differences only 1-5% of fox reads were in disagreement with the dog sequence. these latter differences were deemed less likely to be true nucleotide differences and were excluded from further analysis. this left differences representing one difference for every bases of dog-fox contigs. of these differences, were fox-versus-dog differences, where at least 95% of fox reads were in agreement with each other, and differed from the dog sequence. this represented one difference per mapped bases. an example of one such dog-fox nucleotide difference is presented in supplementary figure 2a. in total, differences were identified among fox reads that aligned to the dog genome, an example of which is presented in supplementary figure 2b. the set of differences where the fox reads differed from each other were filtered to obtain high confidence snps. • snps identified by mapping fox reads against the assembled fox transcriptome sequence a total of fox contigs that were at least bp long, from the assembly of all reads together, were used as a reference for mapping fox reads. in total, differences between fox reads were identified. these polymorphisms included snps, indels, and 2- nucleotide polymorphisms. after applying filtering criteria high-confidence fox snps were identified. contigs containing high-confidence snps were selected and localized in the dog genome by blat, and the output file was visualized on the ucsc genomic browser. • non-redundant set of high-confidence fox-specific snps high confidence fox snps were selected using parameters described in the methods section. snps identified by both methods were checked for redundancy and duplicates were deleted. the non-redundant set contained snps. snps were classified into four categories according to the number of alleles observed in each sample. snps were classified by the observed genotypes in the tame and aggressive individuals, as either homozygous in both, heterozygous in both, or homozygous in one individual and heterozygous in the other. in total, of these high confidence fox snps were localized in the dog genome and the inferred location of of these snps on fox chromosomes was assigned using the comparative fox/dog map. the informativeness of each snp in each sample is also represented on fox chromosome graphs. figure shows the distribution and zygosity of snps in a segment of vvu that has conserved synteny to dog chromosome this region corresponds to part of a qtl interval identified by genetic mapping of behavioural phenotypes in fox experimental pedigrees. note that the aggressive and tame individuals are homozygous for opposite snp alleles in this interval. in the tame animal, the region of homozygosity extends even further. snp validation to establish parameters for filtering high-confidence fox snps we selected a subset of snps for validation by sanger sequencing. three sets of snps were selected. the first set included snps in which the minor allele was present in at least 25% of the sequencing reads per sample and each allele was present in more than three reads in one individual; the second set included snps in which the minor allele was present in at least 25% of sequencing reads per sample and by exactly three reads in one individual; the third set included snps in which the minor allele was present in at least three reads between the two individuals. snps that mapped to the canine x chromosome and showed heterozygosity in at least one sample were filtered out and not included in the snp validation study. of snps, all but one were validated by sanger sequencing. allele heterozygosity, identified by sanger and sequencing for each snp in each individual, differed for two snps in the first dataset, five snps in the second dataset, and two snps in the third data set. for all these snps except snp #3- one of the two individuals was homozygous by sequencing, but heterozygous by sanger sequencing. the allele identified by sanger sequencing corresponded to the second allele observed in another individual. the opposite situation was observed for snp #3- with two snp alleles identified in the tame individual by sequencing but only one by sanger sequencing. all snps in the first two data sets were validated by sanger sequencing, indicating that the selection parameters applied identified high-confidence snps reliably. over 90% of snps in the third data set were also validated. sanger resequencing of snps did find, however, that a second snp allele was not always identified in the sequencing datasets at the coverage cutoff used for selection of high-confidence snps. indeed, all the minor snp alleles not identified by sequencing were among those with the lowest read coverage depth. distribution of fox snps among different types of genomic regions blat of snp-bearing fox amplicons to the dog genome localized a subset of snps to non-coding regions of the dog genome. to evaluate the distribution of fox snps among known coding and non-coding regions we selected random snps that passed defined stringent filtering parameters and localized these snps in the dog genome using the ucsc genome browser. in total, snps were identified in exons or utrs, in introns or non-coding regions, and snps were not localized unambiguously in the dog genome. these results demonstrate that although the majority of high-confidence snps identified in the fox transcriptome correspond to translated sequences, our reads also contain sequences from non-coding genomic regions or regions not currently known to be transcribed in the dog genome. randomly selected high-confidence fox snps were mapped onto the dog genome sequence assembly to identify their distribution among different types of genome sequence. comparison of gene expression between the two fox samples a total of genes had a tenfold difference of expression between the tame and aggressive individuals with p <. genes had at least a twofold difference of expression and p <. genes with differences in expression between the two individuals were classified into biofunctional groups, using ingenuity ipa http://www.ingenuity.com. the most statistically significant biofunction groups for both the tame and the aggressive upregulated genes are shown in figure networks were also analyzed, and a network of genes was found in the set of upregulated genes from the aggressive individual that contains genes involved in the behavior, nervous system development and function and cardiovascular development and function biofunction groups. validation of gene expression by rt-qpcr eleven genes that showed expression differences between the two fox samples, in analysis of transcriptome data, were selected for validation by rt-qpcr. genes were chosen that had more than one exon and did not belong to a highly conserved gene family. differences in expression of selected genes varied between two samples in terms of number of reads, fold differences, and statistical significance of observed expression differences. six genes in this set had higher expression levels in the transcriptome of the aggressive fox, versus five in the tame individual. primers were designed to successfully amplify unique products for nine of these genes, but not successfully for pkhd1l and btnl - the latter two genes were thus excluded from further analysis. the relative expression level was estimated for these nine genes by calculating the difference in expression between "tame" and "aggressive" samples. three genes showed less than a two-fold difference between rt-qpcr results and transcriptome analysis. five genes showed more than a two-fold lower expression difference in rt-qpcr experiments than in the transcriptome analysis, although the difference between the "tame" and "aggressive" samples was in the same direction as in the transcriptome analysis. for four of these genes the difference in expression between the "tame" and "aggressive" samples detected by rt-qpcr was less than two-fold. one gene, slitrk showed a difference of expression by rt-qpcr in the opposite direction from that by transcriptome analysis. expression difference of slitrk in the transcriptome analysis had less significant support among all genes analysed in rt-qpcr experiments. nine genes, identified as differentially expressed between the tame and aggressive fox by transcriptome sequence analysis, were selected for expression evaluation by relative rt-qpcr. results are sorted by z test p-value. except for one gene, which had the least significant z test p-value, rt-qpcr consistently confirmed a difference of expression between the tame and aggressive samples in the same direction as had transcriptome analysis. discussion this study was designed to expand genomic resources for the red fox, an emerging animal model for genetic research whose genome has not yet been sequenced. using the gs flx titanium platform we sequenced non-normalized cdna libraries from brain samples of two farm-bred foxes with markedly different genetically determined behavioral phenotypes. samples for sequencing were selected from prefrontal cortex because this region is easily recognisable in the brain; can thus be collected in a consistent manner; and is related to behavior. the role of different subdivisions of prefrontal cortex in regulating various aspects of cognition, memory, emotions, and social behavior is reviewed in clark et al.,. in humans, a link between brain damage of prefrontal cortex and aggressive behavior has been reported. a total of gb of transcriptome sequencing data for the tame fox and gb for the aggressive fox were obtained. the depth of sequencing chosen for this study was based on previous experience sequencing non-normalized cdna libraries. the fox reads were assembled in three different combinations: i) only reads from the tame fox, ii) only reads from the aggressive fox; iii) all reads together. the percent of reads assigned into contigs in each of the three assemblies was very similar. the assembly software used, the gs de novo assembler, takes into account the existence of gene splice variants. newbler first assembles sets of self-consistent reads into contigs, which are the rough equivalent of exons. it then builds isotigs from contigs that are consistently connected by subpopulations of reads. when there are splice variants, contigs can be connected via different paths that are supported by read data, leading to multiple isotigs being generated. isotigs are grouped into isogroups, which are the rough equivalent of genes. in the context of a transcriptome assembly by newbler, contigs are the sections of isogroups where corresponding isotig sequences are the same. consequently, the assembly produces two sets of sequences: a theoretically non-redundant set of contigs, and a set of isotigs which are candidate long- or full-length transcripts from multiple splice variants of genes, and contains some partially redundant sequence segments. newbler is currently the only transcripome assembler that attempts this task, and is one of the mostly widely used assemblers for de novo assembly of organisms with few genomic recourses. comparison of transcriptome assemblies from different studies is inherently problematic because different assembly algorithms handle the complexities of the transciptome in a variety of ways. cap which was used in the assembly of transcriptomes of two mammalian species -- the bank vole and the domestic ferret -- does not split reads into multiple contigs as newbler does. this makes the number of contigs non-comparable. one measure of comparison is the span of contigs, the total number of bases in the all contigs in an assembly. this number is available for the bank vole, as the number of contigs reported,, multiplied by the average length of the contigs, for an assembly span of mb. assembly of the fox transcriptome spans mb. relative to the bank vole, the fox transcriptome was sequenced to approximately times greater depth and the assembly span was higher by times. the number of isogroups, isotigs, and contigs identified in the assembly of all fox reads together was significantly higher than the number of isogroups, isotigs, and contigs identified in the separate assemblies of tame and aggressive reads. singleton reads in the individual assemblies, when assembled together, created new contigs, isotigs and isogroups. furthermore, sequence variations between individuals increased the number of contigs and isotigs in isogroups, when these variations were sufficiently large for the assembly program to assume them to be different exons, or real differences in splice variants of some genes between the two individuals. the average isotig length and n isotig length were also longer in the assembly of all reads together, as would be expected from an increase in the amount of sequencing data. because the fox genome has not yet been sequenced and a fox refseq database is thus not available, we used swissprot and the dog refseq database to identify gene orthologs in the fox transcriptome. blastx of fox isotigs from the assembly of all reads together against swissprot database and mapping of all fox reads against dog refseq database identified very similar number of genes. the number of identified genes in the fox transcriptome was significantly lower than the number of isogroups found in the assembly of fox reads. the total number of isogroups was also larger than the total number of genes in mammalian genome. one hypothesis is that the amount of sequencing data used in the assembly was not sufficient to cover the full lengths of all fox gene transcripts, therefore some poorly covered genes might have been partitioned into two or more isogroups. consistent with this hypothesis, blastx of fox isotigs against swissprot identified genes, of which were found in more than one isogroup. in addition, a subset of fox reads was mapped to non-gene regions in the dog genome. these regions are expected not to be present in the dog refseq database. we expect that in the fox transcriptome these regions will each be represented by an isogroup with only one isotig assigned to it. also, not all genes that are transcribed are described in the databases we used. the genes not in public databases are expected to be genes which would have the more rare, and less covered sequences, and subsequently less complex isogroups. this low level of expression would imply that they would have been less likely to be found and described previously. this may explain why only % of all isogroups with a single isotig were successfully aligned using blastx. at the same time, % of isogroups containing more than one isotig yielded a hit to the swissprot database. the completeness of the fox transcriptome was tested by mapping fox reads combined variously into datasets of different sizes against the dog refseq database. one striking result of this experiment was the identification of a very similar number of dog gene orthologs in transcriptomes from tame and aggressive samples. the total number of genes identified in both transcriptomes was only approximately 5% larger than a number of genes identified in individual samples. these data demonstrate that, without increasing the depth of sequencing by at least an order of magnitude, it is unlikely that a significantly higher number of rare transcripts would be identified in these transcriptome libraries. the number of genes identified, in this study, from fox pre-frontal brain transcriptome is similar to that identified from mouse cortex. for each accession number in the dog refseq database to which at least one fox read mapped, we estimated the breadth of coverage of dog transcripts by fox reads. figure demonstrates the bivariate coverage of identified transcripts: roughly 75% of dog transcripts had less than 90% coverage and ~25% of dog transcripts had coverage over 90%. we examined the depth of coverage of dog transcripts with relatively low coverage. the slight rise in the middle of the graph, from around 40-60% of coverage contains, as expected, mainly transcripts with a low depth of coverage, but it also contains a significant portion of the transcripts, responsible for the rise above a simple fit curve, that have a large number of fox reads mapped to them. it suggests that in many cases it is not the depth of coverage that is limiting the breadth of coverage. the transcripts in the canine refseq database can originate from different tissues and it is likely that for some genes, brain transcripts only partially overlap with transcripts from other tissues. possible sequence differences between canine and fox transcripts would also reduce the coverage. the breadth of coverage of dog transcripts mapped by reads from two individual samples or from both samples together was almost identical. these data again suggest that a moderate further increase in the amount of sequencing data is unlikely to increase the breadth of coverage of dog transcripts. the dog is the closest relative of the fox with a sequenced genome. the dog and fox, however, have very different karyotypes, the dog having mostly acrocentric chromosomes, whereas the red fox has metacentric chromosomes. the cytological relationship between the dog and fox genomes is well understood: cytogenetic methods and alignment of the fox meiotic map against of the dog genome have demonstrated that most fox chromosomes correspond to two or three canine chromosomes. the canine genome sequence assembly canfam was used to localize fox sequencing reads and contigs in the dog genome. the approximate location of corresponding regions in the fox genome could be then inferred using the comparative dog-fox genetic and cytogenetic maps. mapping fox reads to the canine genome allowed the sequence divergence between the dog and fox genomes to be estimated. one sequence difference between fox and dog was identified per mapped bases. the number of differences between fox chromosomes was one difference per mapped bases. the level of divergence between fox and dog coding regions is comparable with the level of divergence between sheep and cattle. mapping fox reads against both the canine genome and the assembled fox transcriptome sequence identified high-confidence fox-specific snps. the filtering criteria were established by sanger resequencing of snps with different numbers of reads per allele. the number of snps identified in this study significantly exceeded those identified in the bank vole transcriptome, although filtering parameters in our study significantly exceed the filtering parameters used in the vole transcriptome study. these comparisons demonstrate that a greater depth of transcriptome sequencing enables the discovery of a large number of snp's with higher confidence. our filtering parameters yielded 100% confirmation of the polymorphisms in the validation study but differences in individual genotype identified by sequencing were different from sanger sequencing in some cases. to validate snps by sanger sequencing, we used cdna prepared from exactly the same rna samples which were used for sequencing. therefore, we expect allele specific bias to be equal in the samples used in both and sanger sequencing experiments. all snps for which the second allele in the same sample was not identified by sequencing were among snps with the lowest read coverage depth. over 90% of snps from those with the lowest sequence coverage included in the validation study were also validated by sanger sequencing. this filter parameter is still stricter than that in the bank vole study. therefore, for a small risk of false positives, we can also use some of these lower confidence snp's to increase the total set of snps by about 50% to around snps. taking into account the significant cost difference between and sanger technologies, our data suggest that, at the small cost of obtaining a relatively small number of false polymorphisms, snps from this category would be worth including in a snp array. a small number of false snps can be easily detected during post-processing analysis of genotyping calls and filtered out prior to further analysis. for regions in the fox genome associated with traits of interest it would also be useful to invest in validation of snps represented by an even lower number of reads. although the false positive rate for these snps is expected to be higher, validation experiments can still provide a significant number of snps for these critical regions. the technology has been successfully used for developing snp resources for multiple organisms with limited genomic resources. initially, a large number of snp's that mapped to the canine x chromosome were heterozygous in at least one male individual. to test whether this represented errors in our data or in the assembly of the canine x chromosome, the contigs bearing these snps were mapped onto the human x chromosome assembly. all fox contigs containing snps that mapped to both the human and the dog x chromosome outside of the recognized pseudoautosomal region were homozygous for each individual fox. to gain a better understanding of the high heterozygosity of snps assigned to cfax, we mapped approximately bp region surrounding randomly selected snps that were reported as heterozygous on the x back to the dog genome. we found that, in 34% of cases the area surrounding the snp could not be located unambiguously to the x, or to a single area on the x; 50% of snps had high homology to the x but also had significant homology with another area of the genome. this casts doubt on the x chromosome snps showing heterozygosity in our male individuals. we subsequently filtered out all snps on the cfax that were not in the pseudoautosomal region and were heterozygous in either individual. of all the fox high-confidence snps identified % were localised in the dog genome, and the inferred location in the fox genome was identified for %. heterozygosity of snps in the two fox samples was compared and several continuous regions along fox chromosomes where snps were fixed for opposite alleles between tame and aggressive sample were revealed. these regions warrant investigation in a larger number of samples from tame and aggressive fox strains for evaluation of putative selection sweeps. one extended region of reduced heterozygosity in the tame sample was observed in a region of vvu corresponding to part of a previously identified qtl interval. overall, the number of high confidence snps and their locations on fox chromosomes indicated that now we have a set of snps sufficient for genome wide mapping in these populations. detailed evaluation of the location in the dog genome of a randomly selected set of fox snps revealed that approximately 77% of snps are located in transcribed regions and distributed approximately equally between coding regions and utrs. approximately 9% of snps are located in introns and 14% in non-coding regions. similar results were also identified in the bank vole transcriptome study. as has become apparent recently, the mammalian transcriptome is more complex than previously recognized. comparison of gene expression between the two fox samples revealed only a relatively small number of genes whose expression level differed by two-fold or more. the biofunction groups over-represented in the set of differentially expressed genes differed between the two samples. genes associated with neurological diseases were over represented in the tame sample. interestingly, the second biofunction group over-represented in the aggressive sample represents genes associated with cardiovascular diseases. the links between the cardiovascular system and behavior have been increasingly recognized recently. to validate differences in gene expression that were identified by transcriptome sequencing, we selected a subset of nine genes for expression analysis by rt-qpcr. exactly the same two rna samples which had been used for transcriptome sequencing were used for rt-qpcr validation experiments. eight of nine differentially expressed genes selected for the validation study showed differences in expression in the same direction by both methods. all genes that were validated by rt-qpcr showed differential expression between the two samples in transcriptome analysis with p <. although the consistent validation by rt-qpcr of results obtained by transcriptome analysis increases confidence in the observed differences in expression between the two samples, it will clearly be necessary to test these observed differences in a larger set of samples before coming to conclusions regarding biological relevance. among the genes that were validated by rt-qpcr, htr2c is particularly noteworthy. htr2c plays an important role in serotonergic and dopaminergic signalling; and is differentially expressed in specific brain regions of two strains of rats - one of which was selected for and the other against aggressive behavior. increased level of htr2c mrna in the frontal cortex and hippocampus of tame rats was observed. notably, the direction of difference in htr2c expression between the tame and aggressive fox individuals in the present study was consistent with that between rat strains in the popova et al. study. none of the genes reported in other studies investigating comparative gene expression in canid brain samples were significantly different between the two fox samples in the present study. further resequencing and deeper analysis of brain transcriptome, including identification of non-synonymous mutations, comparison of splice variants between tame and aggressive samples, and evaluation of differentially expressed genes in larger independent sample sets will be necessary to characterize critical differences in these fox strains associated with markedly different behavioral phenotypes. in any experimental system there will be a certain amount of noise that might stem from a variety of sources. these biological replicates would let us explore the differences in the underlying biological system we are interested in and validate the differences in gene expression.  using new sequencing technologies it is now possible to develop genomic resources for species without a sequenced genome. that dramatically facilitates the study of interesting but previously inaccessible biological phenomena and helps to establish new animal models for genetic research. overall, transcriptome sequencing from the prefrontal cortex of two fox individuals identified over high-confidence fox-specific snps, fox orthologs of over dog genes, and yielded new insights into potentially important differences in expression of genes in the pre-frontal cortex between tame and aggressive foxes. 
mapping fox reads against the fox transcriptome assembly and the dog genome identified over high confidence fox-specific snps. an at least 2-fold expression difference between the two samples was observed for genes, fewer than 3% of the total number of genes identified in the fox transcriptome. two strains of the silver fox, with markedly different behavioral phenotypes, have been developed by long-term selection for behavior. to understand the genetic differences underlying these behavioral phenotypes fox-specific genomic resources are needed. transcriptome sequencing significantly expanded genomic resources available for the fox, a species without a sequenced genome. this result demonstrates the utility of this approach for developing genomic resources in species with limited genomic information. in a very cost efficient manner this yielded a large number of fox-specific snp markers for genetic studies and provided significant insights into the gene expression profile of the fox pre-frontal cortex; expression differences between the two fox samples; and a catalogue of potentially important gene-specific sequence variants. cdna from mrna from pre-frontal cortex of a tame and an aggressive fox was sequenced using the roche flx titanium platform. foxes from the tame strain exhibit friendly behavior towards humans, paralleling the sociability of canine puppies, whereas foxes from the aggressive strain are defensive and exhibit aggression to humans. fox transcripts for approximately genes were identified using swissprot and the dog refseq databases.
 two strains of the silver fox, with markedly different behavioral phenotypes, have been developed by long-term selection for behavior. foxes from the tame strain exhibit friendly behavior towards humans, paralleling the sociability of canine puppies, whereas foxes from the aggressive strain are defensive and exhibit aggression to humans. to understand the genetic differences underlying these behavioral phenotypes fox-specific genomic resources are needed.  cdna from mrna from pre-frontal cortex of a tame and an aggressive fox was sequenced using the roche flx titanium platform. over 80% of the fox reads were assembled into contigs. mapping fox reads against the fox transcriptome assembly and the dog genome identified over high confidence fox-specific snps. fox transcripts for approximately genes were identified using swissprot and the dog refseq databases. an at least 2-fold expression difference between the two samples was observed for genes, fewer than 3% of the total number of genes identified in the fox transcriptome.  transcriptome sequencing significantly expanded genomic resources available for the fox, a species without a sequenced genome. in a very cost efficient manner this yielded a large number of fox-specific snp markers for genetic studies and provided significant insights into the gene expression profile of the fox pre-frontal cortex; expression differences between the two fox samples; and a catalogue of potentially important gene-specific sequence variants. this result demonstrates the utility of this approach for developing genomic resources in species with limited genomic information. 