BACKGROUND
transcriptome sequencing describes global trends in gene expression while also detailing alterations to biological pathways at the gene-specific level  <cit> .
high-throughput sequencing of rna  quantifies transcripts expressed from known genes with an enormous dynamic range and discovers transcriptional units that are biologically novel yet previously unannotated, or not fully characterized in available databases or gene expression arrays  <cit> .
the localization of transcript sequence to different areas of a gene  at base-pair resolution can detect instances of alternative splicing and transcript isoforms while sequencing of multiple rna species allows for detection of many genetic regulatory elements of biological significance  <cit> .
these advantages make rna-seq a suitable tool for examining the biology of preimplantation embryos.
while single cell transcriptome sequencing is becoming more common, sequencing of individual embryos is not typically performed due to a scarcity of biological material necessary for sequencing library preparation.
this limitation has previously been overcome by pooling of embryos, which precludes important aspects of genome biology from being examined; namely analysis of an individual sample’s genetic variation and identification of allele-specific expression .
the ability to ascertain expression of known and novel snp and to detect imbalanced allelic expression, in addition to discrete quantification of genome-wide transcript abundance, give rna-seq of individual embryos enormous utility in studying early developmental processes.
detection and categorization of snp within production animal systems has been performed extensively  <cit> .
use of these variants as markers and predictors of performance in a large variety of traits  is common  <cit> .
in the context of early development, classifying snp between samples of varying viability, sex, or breed allows for discovery of novel markers of fertility and characterization of critical regulatory mechanisms of embryonic development, such as epigenetic reprogramming or embryonic genome activation.
use of transcript sequence for variant detection has been performed with various assays  <cit> , but also recently in cattle using rna sequencing data  <cit> .
while identification of sequence variation identity at nucleotide resolution and position is valuable, determination of allelic imbalances , defined as a shift from a 50: <dig> expressed ratio, require transcript sequencing of single samples.
expression imbalances have been associated with variation in performance traits and disease processes  <cit> , and could be relevant to embryonic development  <cit> .
ai provides another means of genetic regulation and is a characteristic of imprinted genes, although detection of these with transcript sequencing remains controversial  <cit> .
to examine the feasibility of transcript sequencing of single embryos we performed rna-seq of bovine ivf-derived blastocysts.
particular importance was placed on the analysis options available for this type of data and consideration was given to sequence read processing, bovine genome assemblies available for alignment, and mapping strategy.
transcriptome profiling of replicates was found to be highly reproducible and genes and pathways associated with embryonic samples were highly expressed or enriched.
variant analysis was also performed with in silico validation of detected snp to define criteria for discerning true variation from sequencing artifacts.
similarly, statistical tests and skew thresholds were defined for classification of ai.
the ability to locate genetic variation on a global scale and also quantify the expression of allelic variants demonstrates the unique advantages afforded by sequencing of individual pre-implantation embryos.
RESULTS
sequencing library preparation
the initial obstacle to performing rna-seq from a single blastocyst is obtaining sufficient amounts of high quality rna for use in sequencing library preparation.
to accomplish this we used a standard column-based method of total rna isolation, including a dnase treatment, and eluted the rna into a small volume to achieve a concentration suitable for amplification and library preparation.
we obtained between  <dig>  and  <dig>  ng of total rna per embryo .
these amounts correspond closely to other reports using an identical methodology  <cit>  and the same rna extraction kit  <cit>  used presently.
examination of the 18s and 28s rrna fractions by micro-electrophoresis showed profiles similar to those previously observed in total rna derived from bovine blastocysts  <cit>  and possessing high rqi quality scores .
the amount of rna harvested from a single embryo is not sufficient for preparation of sequencing libraries using standard methods and necessitates amplification following cdna synthesis.
methods for rna amplification have been developed including pcr, in vitro transcription , and ribo-spia -based amplification methods.
we chose the latter approach given that it is recommended for the levels of rna input we obtained from blastocysts and requires only one round of amplification to produce sufficient amounts of material for downstream library preparation.the ribo-spia amplification method is based on an isothermic reaction where transcription initiation sites are primed with random and oligo-dt primers and synthesis of single-stranded cdna occurs via degradation of rna in a dna:rna primer hybrid to create a template permitting multiple instances of transcription  <cit> .
this method has been extensively validated in microarray  <cit> , large-scale rt-pcr  <cit> , and recently rna-seq  <cit>  studies.
after rna amplification starting from ~ <dig> ng of total rna, we obtained on average  <dig> μg of cdna.
amplification produced cdna fragments with an average size range of  <dig> to  <dig> bp, therefore physical or enzymatic fragmentation and size exclusion by gel electrophoresis was not necessary.
no difference in transcript coverage was observed when libraries were prepared from sheared or non-sheared cdna from the ovation system in other studies  <cit> .
the cdna output from the ovation amplification kit was then used for sequencing library preparation, at which point sequencing adapters are ligated to cdna fragments.
electropherograms obtained from amplified cdna and sequencing libraries indicated that high quality libraries were obtained .
previous attempts to generate rna-seq libraries from single bovine embryos were unsuccessful  <cit> , this probably resulted from the lower output of amplified material generated by the ivt-based approach, even when two rounds of amplification were performed.
conversely, the spia approach has been shown to produce higher outputs while maintaining a linear amplification pattern  <cit> .
normally, the high proportion of ribosomal rna  compared to mrna represents a problem for rna-seq studies, requiring steps aimed at eliminating the rrna fraction before library construction and sequencing.
however, human and mouse samples prepared using the ovation amplification system contain low percentages of rrna fragments .
alignment of unmapped reads to a non-coding rna database  resulted in a  <dig> % alignment rate, confirming that the bovine embryo libraries prepared using the ovation v <dig> system did not contain large proportions of rrna fractions.
while the exact mechanism by which rrna are not amplified is currently unknown, it is speculated that the secondary structure of rrna is likely responsible for inefficient primer binding resulting in low cdna conversion efficiency and minimization of rrna amplification.
it is also important to note that this amplification methodology excludes small rna populations below  <dig> bp in size, such as micrornas, meaning that a variety of small non-coding rnas with potential significance to the blastocyst are not included in the sequencing data.
as a confirmation of this, mirna transcripts were not detected in any of our samples.
overall, rna extraction from samples with typically restrictive amounts of rna using a column-based extraction and amplification with ribo-spia allowed for the preparation of high quality sequencing libraries from single bovine blastocysts.
sequence read processing
on average,  <dig> ±  <dig>  million single reads of  <dig> bp were generated per embryo with an illumina gaiix sequencer.
initial mapping of unprocessed reads to the bovine genome  allowing for two mismatches resulted in only  <dig> % of total reads mapping.
allowing for the possibility of low quality reads at the ends of sequences, we evaluated read quality using galaxy .
no evidence of low quality reads based on illumina’s scoring was observed spanning the  <dig> bp of sequence .
however, base sequence content across the reads indicated that the 5′ ends had a higher proportion of g and c nucleotides that deviated from the rest of the read.
this abnormal pattern spanned the first  <dig> bases sequenced.
a potential source of this problem could be that the ovation rna amplification system uses 9mers to perform the rt and are likely to have high gc content.
also, since amplified reads are not shared, and the sequencing primer is ligated to one end of a read, it would be expected by chance that half of the sequences generated begin sequencing at the 9mer annealing site.
moreover, it is possible that 9mer primer binding is not 100% efficient and thus could incorporate mistakes to the sequence.
thus, we hypothesized that removing the first  <dig> bp of the reads would improve mapping by removing bad sequence.
we tested if removing nucleotides at both or either end, and especially nine bp from the 5′ end, would enhance the proportion of reads mapping to the genome .
trimming  <dig> nucleotides at each end of the reads significantly improved the % of read alignment from  <dig> to 82%.
this was not the result of just removing  <dig> bp from the reads, since trimming  <dig> bp at the 3′ end results in a lower read alignment of 71%.
similarly, removing  <dig> bp from the 3′ end resulted in a low alignment percent , although this was higher than untrimmed read mapping probably as a result of the reduced mapping specificity of shortened reads.
trimming  <dig> bp and  <dig> bp from 5′ end resulted in a larger improvement in alignment percentages , which is likely the result of removing low quality information from the reads.
based on these results, we performed all other analyses with reads trimmed  <dig> bp at the 5′ end.
processing of reads prior to alignment by removal of this priming sequence not only significantly improved the fidelity of read alignment, but also the accuracy of snp analyses .
since the time of these experiments, nugen has incorporated a cdna fragmentation step after amplification, minimizing the chance that a fragment will be sequenced at the primer binding site and thus this trimming may not be necessary in future libraries.
however, our experience indicates that quality control of reads is important and removing affected nucleotide positions, although reducing the amount of total sequence, can greatly enhance the accuracy of the results.
mapping of reads to different genome annotations
multiple assemblies and annotations of the bovine genome exist from different sources .
since these annotations contain different numbers of genes and transcripts we compared general mapping statistics after mapping trimmed reads of each embryo to the annotated assemblies using the rna-seq algorithm of clc genomics allowing for  <dig> mismatches .
the differences in total reads aligned varied by as much as  <dig> million  when comparing ncbi-btau  <dig>  and ensembl-umd <dig>  .
the difference in the proportion of uniquely aligned reads was even greater between ncbi-btau  <dig>  and ensembl-umd <dig>  .
these large differences were observed even though both annotations contained a similar number of genes.
also, differences were observed when comparing different annotations of the same genome assembly .
the ensembl annotation of umd <dig>  resulted in higher proportion of reads mapping even though this annotation contained fewer genes than the ncbi annotation.
these mapping read statistics indicate ensembl-umd <dig>  is a more comprehensive genome annotation, not only in terms of identifying actual genes, but through its increased ability to uniquely localize reads.
it is possible that these annotations contain different proportions of embryonically expressed genes and thus alternative annotations may be better for different tissues.
increased alignment to annotated genes improves the robustness of rna-seq results, thus we selected the ensembl-umd <dig>  annotation for further analysis.
rna-seq mapping statistics
of the total reads, an average of  <dig> %  mapped to annotated portions of the ensembl-umd <dig>  genome.
among these,  <dig> % of the reads mapped to protein coding genes,  <dig> % to mitochondrial rrna and the other  <dig> % was distributed between miscellaneous rnas, pseudogenes, retrotransposed elements, and mitochondrial trnas .
the presence of high levels of mitochondrial rrna and trna was in agreement with the high level of expression of mitochondrial protein coding genes.
out of the  <dig> protein coding genes present in the mitochondrial dna,  <dig> were the most highly expressed genes among all protein coding genes and the other  <dig> ranked at positions  <dig>   <dig>   <dig> and  <dig>  among a total of  <dig>  protein coding genes.
this indicates that a large number of blastocysts transcripts are coded by the mitochondrial genome.
also, it is not clear why mitochondrial, but not nuclear, rrna is detected by the methodology used.
it could be that the three dimensional structure of mitochondrial rrna does not prevent random primer binding and is amplified.
on the other hand, the amount of mitochondrial rrna detected does not interfere with analysis of protein coding rnas, since as indicated above, most reads mapped to protein coding genes.
among reads mapping uniquely to protein coding genes , 60% were located to exons, including 6% to exon-exon boundaries.
39% were located to introns and  <dig> % to exon-intron boundaries .
the relatively high proportion of reads assigned to introns is not uncommon when the sequencing library preparation includes random priming of the mrna  <cit> .
this is not often seen when the rt reaction is performed using oligo-dt primers and only amplifying polyadenylated transcripts.
whether these intron sequences belong to un-processed transcripts or un-annotated exons is not clear and deserves further investigation.
towards the latter possibility, clc genomics identifies and quantifies the presence of putative exons and we identified an average of  <dig>  putative exons per sample.
analysis of the location and identity of these potentially novel coding sequences goes beyond the scope of this manuscript but warrants additional attention.
of the  <dig>  protein coding genes in umd <dig> ,  <dig>  ±  <dig>  were expressed per embryo with at least  <dig>  rpkm, a value slightly more conservative than the  <dig>  rpkm threshold suggested to represent above background expression levels  <cit> .
the maximum rpkm value was  <dig> , indicating a dynamic range of expression spanning  <dig> orders of magnitude.
the correlation of rpkm values considering all genes evaluated was high among individual blastocysts , indicating the consistency of this methodology .
milk somatic cell samples analyzed with a similar approach  had a high correlation among the replicates .
on the other hand, the correlation between embryos and somatic cells was extremely low .
these results indicate a high correlation between biological replicates but not between samples of different origin, suggesting that discrimination between samples is possible based on whole transcriptome analysis.
the number of genes detected in our study could be considered low when contrasted to those from other bovine blastocysts studies .
one important consideration is that our analysis results in number of genes expressed while other publications usually report number of transcripts.
given that multiple transcripts can results from the single gene, the comparison of number of expressed genes to number of expressed transcripts should be made cautiously.
also, it is important to note that the threshold used for considering a gene detected can have large implications in the number of detected genes.
we selected an rpkm >  <dig>  as a conservative approach to the recommended  <dig>   <cit> .
reducing the rpkm threshold results in detection of  <dig>  genes, which is in line with what is reported elsewhere  <cit> .
we compared our rna-seq results from single embryos to those generated by huang and khatib using pools of  <dig> normal bovine blastocysts .
raw data was downloaded and processed according to our bioinformatics pipeline.
at rpkm >  <dig>  a total of  <dig>  genes were detected in the pooled blastocysts dataset, a similar number of genes  was also detected in our dataset, when considering it expressed in at least one embryo.
among these,  <dig>  genes were commonly expressed in all datasets.
when we considered genes consistently expressed in all  <dig> unique embryo samples,  <dig>  transcripts were detected, out of which  <dig>   were also present in the huang and khatib dataset.
also, principal component analysis  and hierarchical clustering of all rpkm values indicated that pooled and single blastocysts clustered together and apart from milk somatic cell samples .
these results indicate a highly similar representation of blastocyst transcripts between the studies; in spite of different embryo production, rna extraction, amplification, and library preparation methodologies.
furthermore, the similarity in results highlights the robustness of rna-seq analysis.
the average median coverage of annotated exons with aligned reads was  <dig> x. this coverage corresponds closely to a depth necessary to detect 73% of moderately abundant transcripts  and identify gene isoforms and novel splice junctions  <cit> .
gene isoforms derived from alternative splicing are a significant contributor to extra-genomic variation  <cit>  and exon-exon junctions have been used to identify instances of alternative splicing  <cit> .
exon-exon junctions contributed on average  <dig> % of the total aligned exon reads and intron-exon reads, potentially representing pre-mrnas, constituted a very low proportion  of reads.
this subset of possible regulatory isoforms could lead to discovery of developmentally-related forms of gene regulation unrecognized in non-embryonic tissues and examination of these isoforms between embryonic developmental stages is needed to determine their significance.
functional annotation of blastocyst expressed genes
we performed functional annotation of genes based on level of expression.
genes expressed in all  <dig> blastocysts were sorted by average rpkm into 1) high; 2) medium high; 3) medium; 4) medium low; and 5) low expression level groups.
functional categories enriched in each group compared to the genome  are shown in table  <dig>  enrichment of constitutive cellular elements such as the cytoskeleton, ribosomes and mitochondria was prevalent in higher expression groups, indicating that the blastocysts are rapidly synthesizing proteins to sustain a high rate of cell division and growth.
also, all proteins required for oxidative phosphorylation were in the high expression level clusters, suggesting this mechanism of energy production is very active in blastocysts.
high and medium high brackets mostly related to cell structure and basic cellular biology functions, respectively.
in medium expressing groups, transcription and related processes  were prominent, while in the lower expression brackets, more nuanced aspects of transcription/translation regulation were found to be enriched .
processes of the lower expression brackets were more varied and functionally specific as evidenced by the presence of specific dna repair mechanisms and regulators of cell death, reflecting the greater variety of cellular functions represented outside of the constitutive structural functions found in the medium and high expression level groups.
to identify overexpressed genes in embryos the transcriptomes of blastocysts and milk somatic cells  were compared.
the choice of milk somatic cell as a reference for comparison was made based on the availability of the data and its expected dissimilarity in global gene expression.
as expected, we identified large differences between these samples  with a total of  <dig>  genes differentially expressed at an adjusted p value <  <dig> .
among these genes,  <dig> were overexpressed in the blastocyst samples by at least a  <dig> fold-change difference.
this group of blastocyst-overexpressed genes was subjected to gene ontology classification and functional analysis using david.
the top overrepresented biological processes corresponded to embryonic development, including blastocyst development .
the blastocyst stage embryo is composed of inner cell mass  and trophectoderm  cell lineages.
the icm is a pluripotent cell lineage that gives rise to all the tissues in the embryo proper, and from which embryonic stem cells are derived.
as such, stem cell maintenance was among the top overrepresented categories.
the te differentiates into placental tissues, and placental development was another one of the top overrepresented biological processes.
interplay between the icm and te ensures that these tissues are correctly specified.
specific transcription factors need to be expressed to differentiate embryonic cells into these specific tissues, and transcription factor activity was the most overrepresented molecular function.
among them, oct <dig>  nanog and sox <dig> are well known transcription factors associated with pluripotency and were all expressed in embryos, but not in milk somatic cells.
sall <dig>  a gene known to regulate oct <dig> expression in mice was detected.
cdx <dig> a transcription factor required for te development was highly expressed in embryos.
also, gata <dig>  a known regulator of trophoblast lineage transcription was present in all five samples, as well its upstream regulator tead <dig>  ifntau, the factor responsible for pregnancy recognition in cows, was also present.
grb <dig> is a known suppressor of nanog in the primitive endoderm  of the inner cell mass , and homozygous mutants of grb <dig> are believed to arrest shortly after implantation  <cit> .
grb <dig> and gata <dig> were expressed in all samples, which suggest that primitive endoderm differentiation may already be active in day  <dig> embryos.
gata <dig>  tead <dig>  and grb <dig> suppress expression of nanog and oct <dig> and are thought to diminish pluripotent gene expression to create and maintain extra embryonic lineages  <cit> .
establishment of new embryonic cell fates requires epigenetic changes, and de novo dna methylation genes  were highly expressed in embryos.
also, post-transcriptional regulatory mechanisms are important for regulating cell fate changes, and lin <dig>  a factor required for embryonic stem cell differentiation through regulation of mirna activity, was highly expressed in blastocysts.
to create the blastocoel cavity the trophectoderm cells form a polarized epithelium and transport ions across them, generating osmotic pressure that moves water into the blastocoel.
cell-cell junction related genes were highly expressed, including genes associated with gap junctions and desmosomes.
also, atpase activity mostly associated with ion pumps such as the na+/k +  and h +  atpases were among the top overrepresented molecular functions in blastocyst genes .
sex-specific gene expression profiles
differences in gene expression have been detected between blastocysts of different gender.
in an attempt to determine the sex of the embryos we analyzed expression of candidate sex-specific genes.
sry, the gene responsible for sex determination, was not detected in any embryo.
on the other hand, other genes associated to the y chromosome, such as eif1ay and ddx3y, were expressed in all embryos except blastocyst  <dig>  furthermore, xist, a gene involved in x-chromosome inactivation in female cells, was expressed at high levels only in blastocyst  <dig>  this suggests that blastocyst  <dig> is a female embryo while the others are male.
the high male to female ratio could be explained by the tendency of male embryos to grow faster under in vitro culture conditions and it is possible that selecting for more advanced embryos at the time of collection resulted in a male bias.
we also noticed that the clustering of samples based on all genes analyzed discriminated blastocyst  <dig> from the other embryos, suggesting that global levels of gene expression can discriminate embryo gender.
in support of this, one third of expressed genes were previously found to be differentially regulated when comparing male and female embryos using microarray analysis of bovine blastocysts produced with sexed semen  <cit> .
we thus performed a statistical analysis comparing blastocyst  <dig> to the rest.
a total of  <dig> genes were differentially expressed between the males and female blastocyst at an adjusted p-value of <  <dig>  and  <dig> fold change.
of these,  <dig> were overexpressed in the female embryo.
among them,  <dig> genes  were located to the x chromosome, a proportion much higher than expected by chance since the x chromosome only contains  <dig> % of annotated genes.
moreover, of the overexpressed genes in males,  <dig> of the  <dig> genes  were homologous to y-linked genes .
this is consistent with results from a microarray study comparing male and female bovine embryos where differentially expressed genes were enriched to the sex chromosomes  <cit> .
interestingly, the microarray study found most differences to be lower than  <dig> fold-changes, whereas our study indicated fold-changes ranging from  <dig>  to >  <dig>  these discrepancies are probably associated with the larger dynamic range of rna-seq versus microarrays and suggest that further analysis using rna-seq in female and male embryos would be informative.
we also compared our results to the list of differentially expressed genes reported by bermejo-alvarez et al.
<cit>  .
for genes reported overexpressed in female embryos, we found that among the  <dig> genes with a common identity to ours,  <dig>  had a fold change >  <dig>  in our data.
similarly, all matching genes reported upregulated in males by bermejo-alvarez et al.
were also higher in our male samples compared to the female one.
finally, among transcripts validated by qpcr in the aforementioned study, all the matching genes  differed in the same direction in our data, with  <dig> of them having a p-value <  <dig> .
the similarities in the results of these two studies reinforce the notion that global gene expression analysis can differentiate between embryos of different sexes.
no functional categories were enriched in the group of genes overexpressed in male embryos.
among female overexpressed genes, no go molecular process was significantly overexpressed while cell adhesion , glucose metabolism  and cell motility  were overrepresented categories among biological functions .
snp identification in single embryo transcriptomes
next generation transcriptome sequencing allows for identification and discovery of genetic variants located in transcribed regions of the genome.
genetic variation in gene coding sequences has a higher potential to affect phenotypic characteristics.
to investigate the usefulness of single embryo transcriptome data for detecting known and novel genetic variation we used clc genomic workbench snp identification tool.
reads were mapped to the umd <dig>  reference genome, and nucleotides within reads that differed from the reference were identified as snp.
an snp was considered homozygous in the sample if only a variant allele was present.
snp were considered heterozygous if both a variant and the reference nucleotide were detected at a given position.
snp validation was performed in-silico by matching putative snp positions to known bovine dbsnp entries .
snp were considered validated when a corresponding dbsnp entry was found and the variant nucleotide identities coincided exactly with our data.
the random chance of exact variant matching is only  <dig> % for heterozygous snp and 33% for homozygous snp.
initially, variants were called if the snp was covered by more than  <dig> reads and, in the case of heterozygous snp, if the lower expressed allele was present in at least  <dig> reads.
roughly, half of the snp that were identified matched a dbsnp variant position, and out of these more than 99% shared the same nucleotide alleles between our samples and the reference database.
when performing the validation it was noticed that a higher proportion of homozygous snps corresponded to already known snp compared to heterozygous snp.
however, when an snp was found in both our data and dbsnp, the allelic variants coincided in > 99% of the cases for both homozygous and heterozygous snp.
we interpreted these findings as a suggestion that the rate of false snp discovery is higher in heterozygous snp .
this could be attributed to a lower threshold required to call the lower variant in a heterozygous snp, since only  <dig> reads are required for the lower allele to be identified as a snp versus  <dig> reads for homozygous snp.
we therefore examined how snp coverage related to validation rate in order to estimate the minimum sequencing depth to accurately call new snp.
the proportion of snp matching dbsnp across a spectrum of coverage levels was calculated and found to increase appreciably from  <dig> to  <dig> reads and then plateau afterwards through  <dig> reads per variant .
therefore,  <dig> reads was selected as the minimum coverage level for novel snp .
increasing the threshold for snp discovery resulted in fewer sequencing errors being determined as snp.
it is still possible however that certain genes with high expression levels could accumulate a sufficient number of errors to reach the minimum threshold for snp identification.
to investigate this possibility, snp detected in mitochondrial genes were examined, because these genes were highly expressed and because the mitochondrial genome is assumed to be homozygous  <cit> .
we assumed that the presence of heterozygosity in the mitochondrial genome is the result of sequencing errors and therefore represents false snp discovery.
when heterozygous snp in mitochondrial genes were evaluated, a total of  <dig> unique snp with an average of  <dig> per sample were observed.
the frequency of the minor allele averaged  <dig>  ±  <dig> %, and was consistent for each of the embryos.
more than 95% of the snp in mitochondrial genes consisted of snp with the minor allele representing less than 15% of the reads.
based on this, 15% was chosen as the lower limit cutoff for minor alleles to limit the probability  that sequencing errors in highly expressed genes influence heterozygous snp discovery and allele specific expression analysis.
thus, coverage thresholds used should be carefully considered in order to minimize calling of sequencing artifacts in the interest of maximizing variation coverage.
it should also be noted that the potential for heteroplasmy in mtdna  <cit>  makes the lower threshold for snp detection very conservative.
the optimized parameters were only used for identification of de-novo snps, while snp detected with the less stringent set of parameters that matched dbsnps were considered valid.
using these criteria, a total of  <dig>  unique snp were identified,  <dig>  , of which were heterozygous in at least one sample .
a large proportion of detected snp  were present in the dbsnp database.
this suggests that snp identified from the transcriptome sequencing of single embryos are highly reliable and can be used to discover novel snp at genome coding sequences and for analysis of allele specific expression.
allele specific expression analysis
analysis of allele specific expression requires the presence of heterozygous snp in the transcribed portions of a gene.
among the  <dig>  heterozygous unique snp that were identified,  <dig>  known genes were represented.
allelic imbalance  was determined for these snp when the proportion of reads assigned to one allele was higher than 65%.
we found that forty percent of expressed heterozygous snp demonstrated an allelic imbalance .
a similar proportion was observed for novel snp, although in this category the highest allowed ratio was 85%.
the statistical significance of these imbalances was assessed in order to diminish the presence of false positives.
statistical significance of expressed ai bias was determined if the read distribution was significantly different than a theoretical 50: <dig> distribution, as determined by a chi <dig> test adjusted for fdr .
statistical thresholds reduced the number of snp with ai to  <dig>  within  <dig> unique genes .
the stringency of the criteria used to assess these significant biases, while minimizing false positives, restricts the ability to create a comprehensive list of snp given that only those with high coverage provided sufficient information for statistical analysis.
indeed, 82% of genes with significant ai were among the top 20th percentile of expression levels, demonstrating the necessity of sufficient read depth for ase detection.
increased sequencing depth would likely allow discrimination of ai in more genes and also provide more statistical power to determine the significance of low level ai.
in order to estimate the number of genes in blastocysts that could possibly express ai, we determined the proportion of expressed ai with statistical significance within a subset of high coverage genes.
of  <dig> genes containing snp with high coverage ,  <dig> were found to also have snp expressing significant ai .
we estimate that with sufficient sequencing depth, ~33% of expressed genes in bovine blastocyst genes could have ai with a major allele frequency > 65%.
differences in epigenetic reprogramming of maternal and paternal genomes during preimplantation development could result in allelic skewing of embryonic gene expression.
ai dependent on parental origin and cellular lineage may also be important for embryonic development  <cit> .
however, it is not known at this time if the preferentially expressed alleles have a common parental origin.
go enrichment analysis of genes containing snp with statistically significant ai identified major cellular constitutive elements as overrepresented categories.
these included go terms such as cytoskeleton, basolateral membrane, focal adhesion, non-membrane-bounded organelles, organelle lumen and macromolecular complex organization, among others.
the enrichment of genes to these specific cellular elements is not surprising given that they are highly expressed cellular components and, as stated above, detection of imbalance is greatly aided by high coverage or expression.
imprinting is one phenomenon that leads to allelic expression imbalance, with typically only one specific parental allele being expressed.
the approach used to detect imbalances is based on the expressed genes without knowledge of the individual embryo genotype, and therefore in cases of monoallelic expression it would be expected that no heterozygous alleles would be identified.
among  <dig> genes known to be imprinted in cattle, human and/or mouse,  <dig> were expressed in bovine blastocysts at rpkm >  <dig>   <cit> .
expressed heterozygous snp were found in  <dig> of these genes .
lack of heterozygous snp in other imprinted genes could indicate mono-allelic expression.
allelic biases  were observed in  <dig> out of the  <dig> detected heterozygous snp .
because of low imprinted gene expression and low coverage of the snp, only  <dig> candidate imprinted snp from  <dig> genes  had ai with statistical significance .
identification of imbalanced alleles in most of the known imprinted genes within our dataset further corroborates the utility of rna-seq in characterization of ai.
* chi-squared p <  <dig>  difference from 50: <dig> expression of both alleles.
CONCLUSIONS
this study reports the analysis of individual bovine embryo transcriptome sequencing, providing details for amplification procedures of low rna input samples and an analysis pipeline that examines differences in gene expression profiles, identifies novel snp and determines instances of allelic imbalance.
rna-seq analysis in single embryos allows for discrimination of embryo gender and provided the opportunity to characterize individual variability in gene expression.
snp analysis of individual samples demonstrates the use of rna-seq to identify embryo-specific variation for association studies and ase that could represent novel layers of developmental regulation subject to influence by ai.
our data suggests that ai is prevalent in bovine blastocysts.
rna-seq analysis of the individual embryonic transcriptome is feasible and presents valuable insights into gene expression, variation and regulation of the early developmental transcriptome.

 transcriptome sequencing describes global trends in gene expression while also detailing alterations to biological pathways at the gene-specific level. high-throughput sequencing of rna quantifies transcripts expressed from known genes with an enormous dynamic range and discovers transcriptional units that are biologically novel yet previously unannotated, or not fully characterized in available databases or gene expression arrays. the localization of transcript sequence to different areas of a gene at base-pair resolution can detect instances of alternative splicing and transcript isoforms while sequencing of multiple rna species allows for detection of many genetic regulatory elements of biological significance. these advantages make rna-seq a suitable tool for examining the biology of preimplantation embryos. while single cell transcriptome sequencing is becoming more common, sequencing of individual embryos is not typically performed due to a scarcity of biological material necessary for sequencing library preparation. this limitation has previously been overcome by pooling of embryos, which precludes important aspects of genome biology from being examined; namely analysis of an individual sample’s genetic variation and identification of allele-specific expression. the ability to ascertain expression of known and novel snp and to detect imbalanced allelic expression, in addition to discrete quantification of genome-wide transcript abundance, give rna-seq of individual embryos enormous utility in studying early developmental processes. detection and categorization of snp within production animal systems has been performed extensively. use of these variants as markers and predictors of performance in a large variety of traits is common. in the context of early development, classifying snp between samples of varying viability, sex, or breed allows for discovery of novel markers of fertility and characterization of critical regulatory mechanisms of embryonic development, such as epigenetic reprogramming or embryonic genome activation. use of transcript sequence for variant detection has been performed with various assays, but also recently in cattle using rna sequencing data. while identification of sequence variation identity at nucleotide resolution and position is valuable, determination of allelic imbalances, defined as a shift from a 50: expressed ratio, require transcript sequencing of single samples. expression imbalances have been associated with variation in performance traits and disease processes, and could be relevant to embryonic development. ai provides another means of genetic regulation and is a characteristic of imprinted genes, although detection of these with transcript sequencing remains controversial. to examine the feasibility of transcript sequencing of single embryos we performed rna-seq of bovine ivf-derived blastocysts. particular importance was placed on the analysis options available for this type of data and consideration was given to sequence read processing, bovine genome assemblies available for alignment, and mapping strategy. transcriptome profiling of replicates was found to be highly reproducible and genes and pathways associated with embryonic samples were highly expressed or enriched. variant analysis was also performed with in silico validation of detected snp to define criteria for discerning true variation from sequencing artifacts. similarly, statistical tests and skew thresholds were defined for classification of ai. the ability to locate genetic variation on a global scale and also quantify the expression of allelic variants demonstrates the unique advantages afforded by sequencing of individual pre-implantation embryos.  sequencing library preparation the initial obstacle to performing rna-seq from a single blastocyst is obtaining sufficient amounts of high quality rna for use in sequencing library preparation. to accomplish this we used a standard column-based method of total rna isolation, including a dnase treatment, and eluted the rna into a small volume to achieve a concentration suitable for amplification and library preparation. we obtained between and ng of total rna per embryo. these amounts correspond closely to other reports using an identical methodology and the same rna extraction kit used presently. examination of the 18s and 28s rrna fractions by micro-electrophoresis showed profiles similar to those previously observed in total rna derived from bovine blastocysts and possessing high rqi quality scores. the amount of rna harvested from a single embryo is not sufficient for preparation of sequencing libraries using standard methods and necessitates amplification following cdna synthesis. methods for rna amplification have been developed including pcr, in vitro transcription, and ribo-spia -based amplification methods. we chose the latter approach given that it is recommended for the levels of rna input we obtained from blastocysts and requires only one round of amplification to produce sufficient amounts of material for downstream library preparation.the ribo-spia amplification method is based on an isothermic reaction where transcription initiation sites are primed with random and oligo-dt primers and synthesis of single-stranded cdna occurs via degradation of rna in a dna:rna primer hybrid to create a template permitting multiple instances of transcription. this method has been extensively validated in microarray, large-scale rt-pcr, and recently rna-seq studies. after rna amplification starting from ~ ng of total rna, we obtained on average μg of cdna. amplification produced cdna fragments with an average size range of to bp, therefore physical or enzymatic fragmentation and size exclusion by gel electrophoresis was not necessary. no difference in transcript coverage was observed when libraries were prepared from sheared or non-sheared cdna from the ovation system in other studies. the cdna output from the ovation amplification kit was then used for sequencing library preparation, at which point sequencing adapters are ligated to cdna fragments. electropherograms obtained from amplified cdna and sequencing libraries indicated that high quality libraries were obtained. previous attempts to generate rna-seq libraries from single bovine embryos were unsuccessful, this probably resulted from the lower output of amplified material generated by the ivt-based approach, even when two rounds of amplification were performed. conversely, the spia approach has been shown to produce higher outputs while maintaining a linear amplification pattern. normally, the high proportion of ribosomal rna compared to mrna represents a problem for rna-seq studies, requiring steps aimed at eliminating the rrna fraction before library construction and sequencing. however, human and mouse samples prepared using the ovation amplification system contain low percentages of rrna fragments. alignment of unmapped reads to a non-coding rna database resulted in a % alignment rate, confirming that the bovine embryo libraries prepared using the ovation v system did not contain large proportions of rrna fractions. while the exact mechanism by which rrna are not amplified is currently unknown, it is speculated that the secondary structure of rrna is likely responsible for inefficient primer binding resulting in low cdna conversion efficiency and minimization of rrna amplification. it is also important to note that this amplification methodology excludes small rna populations below bp in size, such as micrornas, meaning that a variety of small non-coding rnas with potential significance to the blastocyst are not included in the sequencing data. as a confirmation of this, mirna transcripts were not detected in any of our samples. overall, rna extraction from samples with typically restrictive amounts of rna using a column-based extraction and amplification with ribo-spia allowed for the preparation of high quality sequencing libraries from single bovine blastocysts. sequence read processing on average, ± million single reads of bp were generated per embryo with an illumina gaiix sequencer. initial mapping of unprocessed reads to the bovine genome allowing for two mismatches resulted in only % of total reads mapping. allowing for the possibility of low quality reads at the ends of sequences, we evaluated read quality using galaxy. no evidence of low quality reads based on illumina’s scoring was observed spanning the bp of sequence. however, base sequence content across the reads indicated that the 5′ ends had a higher proportion of g and c nucleotides that deviated from the rest of the read. this abnormal pattern spanned the first bases sequenced. a potential source of this problem could be that the ovation rna amplification system uses 9mers to perform the rt and are likely to have high gc content. also, since amplified reads are not shared, and the sequencing primer is ligated to one end of a read, it would be expected by chance that half of the sequences generated begin sequencing at the 9mer annealing site. moreover, it is possible that 9mer primer binding is not 100% efficient and thus could incorporate mistakes to the sequence. thus, we hypothesized that removing the first bp of the reads would improve mapping by removing bad sequence. we tested if removing nucleotides at both or either end, and especially nine bp from the 5′ end, would enhance the proportion of reads mapping to the genome. trimming nucleotides at each end of the reads significantly improved the % of read alignment from to 82%. this was not the result of just removing bp from the reads, since trimming bp at the 3′ end results in a lower read alignment of 71%. similarly, removing bp from the 3′ end resulted in a low alignment percent, although this was higher than untrimmed read mapping probably as a result of the reduced mapping specificity of shortened reads. trimming bp and bp from 5′ end resulted in a larger improvement in alignment percentages, which is likely the result of removing low quality information from the reads. based on these results, we performed all other analyses with reads trimmed bp at the 5′ end. processing of reads prior to alignment by removal of this priming sequence not only significantly improved the fidelity of read alignment, but also the accuracy of snp analyses. since the time of these experiments, nugen has incorporated a cdna fragmentation step after amplification, minimizing the chance that a fragment will be sequenced at the primer binding site and thus this trimming may not be necessary in future libraries. however, our experience indicates that quality control of reads is important and removing affected nucleotide positions, although reducing the amount of total sequence, can greatly enhance the accuracy of the results. mapping of reads to different genome annotations multiple assemblies and annotations of the bovine genome exist from different sources. since these annotations contain different numbers of genes and transcripts we compared general mapping statistics after mapping trimmed reads of each embryo to the annotated assemblies using the rna-seq algorithm of clc genomics allowing for mismatches. the differences in total reads aligned varied by as much as million when comparing ncbi-btau and ensembl-umd. the difference in the proportion of uniquely aligned reads was even greater between ncbi-btau and ensembl-umd. these large differences were observed even though both annotations contained a similar number of genes. also, differences were observed when comparing different annotations of the same genome assembly. the ensembl annotation of umd resulted in higher proportion of reads mapping even though this annotation contained fewer genes than the ncbi annotation. these mapping read statistics indicate ensembl-umd is a more comprehensive genome annotation, not only in terms of identifying actual genes, but through its increased ability to uniquely localize reads. it is possible that these annotations contain different proportions of embryonically expressed genes and thus alternative annotations may be better for different tissues. increased alignment to annotated genes improves the robustness of rna-seq results, thus we selected the ensembl-umd annotation for further analysis. rna-seq mapping statistics of the total reads, an average of % mapped to annotated portions of the ensembl-umd genome. among these, % of the reads mapped to protein coding genes, % to mitochondrial rrna and the other % was distributed between miscellaneous rnas, pseudogenes, retrotransposed elements, and mitochondrial trnas. the presence of high levels of mitochondrial rrna and trna was in agreement with the high level of expression of mitochondrial protein coding genes. out of the protein coding genes present in the mitochondrial dna, were the most highly expressed genes among all protein coding genes and the other ranked at positions and among a total of protein coding genes. this indicates that a large number of blastocysts transcripts are coded by the mitochondrial genome. also, it is not clear why mitochondrial, but not nuclear, rrna is detected by the methodology used. it could be that the three dimensional structure of mitochondrial rrna does not prevent random primer binding and is amplified. on the other hand, the amount of mitochondrial rrna detected does not interfere with analysis of protein coding rnas, since as indicated above, most reads mapped to protein coding genes. among reads mapping uniquely to protein coding genes, 60% were located to exons, including 6% to exon-exon boundaries. 39% were located to introns and % to exon-intron boundaries. the relatively high proportion of reads assigned to introns is not uncommon when the sequencing library preparation includes random priming of the mrna. this is not often seen when the rt reaction is performed using oligo-dt primers and only amplifying polyadenylated transcripts. whether these intron sequences belong to un-processed transcripts or un-annotated exons is not clear and deserves further investigation. towards the latter possibility, clc genomics identifies and quantifies the presence of putative exons and we identified an average of putative exons per sample. analysis of the location and identity of these potentially novel coding sequences goes beyond the scope of this manuscript but warrants additional attention. of the protein coding genes in umd, ± were expressed per embryo with at least rpkm, a value slightly more conservative than the rpkm threshold suggested to represent above background expression levels. the maximum rpkm value was, indicating a dynamic range of expression spanning orders of magnitude. the correlation of rpkm values considering all genes evaluated was high among individual blastocysts, indicating the consistency of this methodology. milk somatic cell samples analyzed with a similar approach had a high correlation among the replicates. on the other hand, the correlation between embryos and somatic cells was extremely low. these results indicate a high correlation between biological replicates but not between samples of different origin, suggesting that discrimination between samples is possible based on whole transcriptome analysis. the number of genes detected in our study could be considered low when contrasted to those from other bovine blastocysts studies. one important consideration is that our analysis results in number of genes expressed while other publications usually report number of transcripts. given that multiple transcripts can results from the single gene, the comparison of number of expressed genes to number of expressed transcripts should be made cautiously. also, it is important to note that the threshold used for considering a gene detected can have large implications in the number of detected genes. we selected an rpkm > as a conservative approach to the recommended. reducing the rpkm threshold results in detection of genes, which is in line with what is reported elsewhere. we compared our rna-seq results from single embryos to those generated by huang and khatib using pools of normal bovine blastocysts. raw data was downloaded and processed according to our bioinformatics pipeline. at rpkm > a total of genes were detected in the pooled blastocysts dataset, a similar number of genes was also detected in our dataset, when considering it expressed in at least one embryo. among these, genes were commonly expressed in all datasets. when we considered genes consistently expressed in all unique embryo samples, transcripts were detected, out of which were also present in the huang and khatib dataset. also, principal component analysis and hierarchical clustering of all rpkm values indicated that pooled and single blastocysts clustered together and apart from milk somatic cell samples. these results indicate a highly similar representation of blastocyst transcripts between the studies; in spite of different embryo production, rna extraction, amplification, and library preparation methodologies. furthermore, the similarity in results highlights the robustness of rna-seq analysis. the average median coverage of annotated exons with aligned reads was x. this coverage corresponds closely to a depth necessary to detect 73% of moderately abundant transcripts and identify gene isoforms and novel splice junctions. gene isoforms derived from alternative splicing are a significant contributor to extra-genomic variation and exon-exon junctions have been used to identify instances of alternative splicing. exon-exon junctions contributed on average % of the total aligned exon reads and intron-exon reads, potentially representing pre-mrnas, constituted a very low proportion of reads. this subset of possible regulatory isoforms could lead to discovery of developmentally-related forms of gene regulation unrecognized in non-embryonic tissues and examination of these isoforms between embryonic developmental stages is needed to determine their significance. functional annotation of blastocyst expressed genes we performed functional annotation of genes based on level of expression. genes expressed in all blastocysts were sorted by average rpkm into 1) high; 2) medium high; 3) medium; 4) medium low; and 5) low expression level groups. functional categories enriched in each group compared to the genome are shown in table  enrichment of constitutive cellular elements such as the cytoskeleton, ribosomes and mitochondria was prevalent in higher expression groups, indicating that the blastocysts are rapidly synthesizing proteins to sustain a high rate of cell division and growth. also, all proteins required for oxidative phosphorylation were in the high expression level clusters, suggesting this mechanism of energy production is very active in blastocysts. high and medium high brackets mostly related to cell structure and basic cellular biology functions, respectively. in medium expressing groups, transcription and related processes were prominent, while in the lower expression brackets, more nuanced aspects of transcription/translation regulation were found to be enriched. processes of the lower expression brackets were more varied and functionally specific as evidenced by the presence of specific dna repair mechanisms and regulators of cell death, reflecting the greater variety of cellular functions represented outside of the constitutive structural functions found in the medium and high expression level groups. to identify overexpressed genes in embryos the transcriptomes of blastocysts and milk somatic cells were compared. the choice of milk somatic cell as a reference for comparison was made based on the availability of the data and its expected dissimilarity in global gene expression. as expected, we identified large differences between these samples with a total of genes differentially expressed at an adjusted p value <. among these genes, were overexpressed in the blastocyst samples by at least a fold-change difference. this group of blastocyst-overexpressed genes was subjected to gene ontology classification and functional analysis using david. the top overrepresented biological processes corresponded to embryonic development, including blastocyst development. the blastocyst stage embryo is composed of inner cell mass and trophectoderm cell lineages. the icm is a pluripotent cell lineage that gives rise to all the tissues in the embryo proper, and from which embryonic stem cells are derived. as such, stem cell maintenance was among the top overrepresented categories. the te differentiates into placental tissues, and placental development was another one of the top overrepresented biological processes. interplay between the icm and te ensures that these tissues are correctly specified. specific transcription factors need to be expressed to differentiate embryonic cells into these specific tissues, and transcription factor activity was the most overrepresented molecular function. among them, oct nanog and sox are well known transcription factors associated with pluripotency and were all expressed in embryos, but not in milk somatic cells. sall a gene known to regulate oct expression in mice was detected. cdx a transcription factor required for te development was highly expressed in embryos. also, gata a known regulator of trophoblast lineage transcription was present in all five samples, as well its upstream regulator tead ifntau, the factor responsible for pregnancy recognition in cows, was also present. grb is a known suppressor of nanog in the primitive endoderm of the inner cell mass, and homozygous mutants of grb are believed to arrest shortly after implantation. grb and gata were expressed in all samples, which suggest that primitive endoderm differentiation may already be active in day embryos. gata tead and grb suppress expression of nanog and oct and are thought to diminish pluripotent gene expression to create and maintain extra embryonic lineages. establishment of new embryonic cell fates requires epigenetic changes, and de novo dna methylation genes were highly expressed in embryos. also, post-transcriptional regulatory mechanisms are important for regulating cell fate changes, and lin a factor required for embryonic stem cell differentiation through regulation of mirna activity, was highly expressed in blastocysts. to create the blastocoel cavity the trophectoderm cells form a polarized epithelium and transport ions across them, generating osmotic pressure that moves water into the blastocoel. cell-cell junction related genes were highly expressed, including genes associated with gap junctions and desmosomes. also, atpase activity mostly associated with ion pumps such as the na+/k + and h + atpases were among the top overrepresented molecular functions in blastocyst genes. sex-specific gene expression profiles differences in gene expression have been detected between blastocysts of different gender. in an attempt to determine the sex of the embryos we analyzed expression of candidate sex-specific genes. sry, the gene responsible for sex determination, was not detected in any embryo. on the other hand, other genes associated to the y chromosome, such as eif1ay and ddx3y, were expressed in all embryos except blastocyst furthermore, xist, a gene involved in x-chromosome inactivation in female cells, was expressed at high levels only in blastocyst this suggests that blastocyst is a female embryo while the others are male. the high male to female ratio could be explained by the tendency of male embryos to grow faster under in vitro culture conditions and it is possible that selecting for more advanced embryos at the time of collection resulted in a male bias. we also noticed that the clustering of samples based on all genes analyzed discriminated blastocyst from the other embryos, suggesting that global levels of gene expression can discriminate embryo gender. in support of this, one third of expressed genes were previously found to be differentially regulated when comparing male and female embryos using microarray analysis of bovine blastocysts produced with sexed semen. we thus performed a statistical analysis comparing blastocyst to the rest. a total of genes were differentially expressed between the males and female blastocyst at an adjusted p-value of < and fold change. of these, were overexpressed in the female embryo. among them, genes were located to the x chromosome, a proportion much higher than expected by chance since the x chromosome only contains % of annotated genes. moreover, of the overexpressed genes in males, of the genes were homologous to y-linked genes. this is consistent with results from a microarray study comparing male and female bovine embryos where differentially expressed genes were enriched to the sex chromosomes. interestingly, the microarray study found most differences to be lower than fold-changes, whereas our study indicated fold-changes ranging from to > these discrepancies are probably associated with the larger dynamic range of rna-seq versus microarrays and suggest that further analysis using rna-seq in female and male embryos would be informative. we also compared our results to the list of differentially expressed genes reported by bermejo-alvarez et al.. for genes reported overexpressed in female embryos, we found that among the genes with a common identity to ours, had a fold change > in our data. similarly, all matching genes reported upregulated in males by bermejo-alvarez et al. were also higher in our male samples compared to the female one. finally, among transcripts validated by qpcr in the aforementioned study, all the matching genes differed in the same direction in our data, with of them having a p-value <. the similarities in the results of these two studies reinforce the notion that global gene expression analysis can differentiate between embryos of different sexes. no functional categories were enriched in the group of genes overexpressed in male embryos. among female overexpressed genes, no go molecular process was significantly overexpressed while cell adhesion, glucose metabolism and cell motility were overrepresented categories among biological functions. snp identification in single embryo transcriptomes next generation transcriptome sequencing allows for identification and discovery of genetic variants located in transcribed regions of the genome. genetic variation in gene coding sequences has a higher potential to affect phenotypic characteristics. to investigate the usefulness of single embryo transcriptome data for detecting known and novel genetic variation we used clc genomic workbench snp identification tool. reads were mapped to the umd reference genome, and nucleotides within reads that differed from the reference were identified as snp. an snp was considered homozygous in the sample if only a variant allele was present. snp were considered heterozygous if both a variant and the reference nucleotide were detected at a given position. snp validation was performed in-silico by matching putative snp positions to known bovine dbsnp entries. snp were considered validated when a corresponding dbsnp entry was found and the variant nucleotide identities coincided exactly with our data. the random chance of exact variant matching is only % for heterozygous snp and 33% for homozygous snp. initially, variants were called if the snp was covered by more than reads and, in the case of heterozygous snp, if the lower expressed allele was present in at least reads. roughly, half of the snp that were identified matched a dbsnp variant position, and out of these more than 99% shared the same nucleotide alleles between our samples and the reference database. when performing the validation it was noticed that a higher proportion of homozygous snps corresponded to already known snp compared to heterozygous snp. however, when an snp was found in both our data and dbsnp, the allelic variants coincided in > 99% of the cases for both homozygous and heterozygous snp. we interpreted these findings as a suggestion that the rate of false snp discovery is higher in heterozygous snp. this could be attributed to a lower threshold required to call the lower variant in a heterozygous snp, since only reads are required for the lower allele to be identified as a snp versus reads for homozygous snp. we therefore examined how snp coverage related to validation rate in order to estimate the minimum sequencing depth to accurately call new snp. the proportion of snp matching dbsnp across a spectrum of coverage levels was calculated and found to increase appreciably from to reads and then plateau afterwards through reads per variant. therefore, reads was selected as the minimum coverage level for novel snp. increasing the threshold for snp discovery resulted in fewer sequencing errors being determined as snp. it is still possible however that certain genes with high expression levels could accumulate a sufficient number of errors to reach the minimum threshold for snp identification. to investigate this possibility, snp detected in mitochondrial genes were examined, because these genes were highly expressed and because the mitochondrial genome is assumed to be homozygous. we assumed that the presence of heterozygosity in the mitochondrial genome is the result of sequencing errors and therefore represents false snp discovery. when heterozygous snp in mitochondrial genes were evaluated, a total of unique snp with an average of per sample were observed. the frequency of the minor allele averaged ± %, and was consistent for each of the embryos. more than 95% of the snp in mitochondrial genes consisted of snp with the minor allele representing less than 15% of the reads. based on this, 15% was chosen as the lower limit cutoff for minor alleles to limit the probability that sequencing errors in highly expressed genes influence heterozygous snp discovery and allele specific expression analysis. thus, coverage thresholds used should be carefully considered in order to minimize calling of sequencing artifacts in the interest of maximizing variation coverage. it should also be noted that the potential for heteroplasmy in mtdna makes the lower threshold for snp detection very conservative. the optimized parameters were only used for identification of de-novo snps, while snp detected with the less stringent set of parameters that matched dbsnps were considered valid. using these criteria, a total of unique snp were identified,, of which were heterozygous in at least one sample. a large proportion of detected snp were present in the dbsnp database. this suggests that snp identified from the transcriptome sequencing of single embryos are highly reliable and can be used to discover novel snp at genome coding sequences and for analysis of allele specific expression. allele specific expression analysis analysis of allele specific expression requires the presence of heterozygous snp in the transcribed portions of a gene. among the heterozygous unique snp that were identified, known genes were represented. allelic imbalance was determined for these snp when the proportion of reads assigned to one allele was higher than 65%. we found that forty percent of expressed heterozygous snp demonstrated an allelic imbalance. a similar proportion was observed for novel snp, although in this category the highest allowed ratio was 85%. the statistical significance of these imbalances was assessed in order to diminish the presence of false positives. statistical significance of expressed ai bias was determined if the read distribution was significantly different than a theoretical 50: distribution, as determined by a chi test adjusted for fdr. statistical thresholds reduced the number of snp with ai to within unique genes. the stringency of the criteria used to assess these significant biases, while minimizing false positives, restricts the ability to create a comprehensive list of snp given that only those with high coverage provided sufficient information for statistical analysis. indeed, 82% of genes with significant ai were among the top 20th percentile of expression levels, demonstrating the necessity of sufficient read depth for ase detection. increased sequencing depth would likely allow discrimination of ai in more genes and also provide more statistical power to determine the significance of low level ai. in order to estimate the number of genes in blastocysts that could possibly express ai, we determined the proportion of expressed ai with statistical significance within a subset of high coverage genes. of genes containing snp with high coverage, were found to also have snp expressing significant ai. we estimate that with sufficient sequencing depth, ~33% of expressed genes in bovine blastocyst genes could have ai with a major allele frequency > 65%. differences in epigenetic reprogramming of maternal and paternal genomes during preimplantation development could result in allelic skewing of embryonic gene expression. ai dependent on parental origin and cellular lineage may also be important for embryonic development. however, it is not known at this time if the preferentially expressed alleles have a common parental origin. go enrichment analysis of genes containing snp with statistically significant ai identified major cellular constitutive elements as overrepresented categories. these included go terms such as cytoskeleton, basolateral membrane, focal adhesion, non-membrane-bounded organelles, organelle lumen and macromolecular complex organization, among others. the enrichment of genes to these specific cellular elements is not surprising given that they are highly expressed cellular components and, as stated above, detection of imbalance is greatly aided by high coverage or expression. imprinting is one phenomenon that leads to allelic expression imbalance, with typically only one specific parental allele being expressed. the approach used to detect imbalances is based on the expressed genes without knowledge of the individual embryo genotype, and therefore in cases of monoallelic expression it would be expected that no heterozygous alleles would be identified. among genes known to be imprinted in cattle, human and/or mouse, were expressed in bovine blastocysts at rpkm >. expressed heterozygous snp were found in of these genes. lack of heterozygous snp in other imprinted genes could indicate mono-allelic expression. allelic biases were observed in out of the detected heterozygous snp. because of low imprinted gene expression and low coverage of the snp, only candidate imprinted snp from genes had ai with statistical significance. identification of imbalanced alleles in most of the known imprinted genes within our dataset further corroborates the utility of rna-seq in characterization of ai. * chi-squared p < difference from 50: expression of both alleles.  this study reports the analysis of individual bovine embryo transcriptome sequencing, providing details for amplification procedures of low rna input samples and an analysis pipeline that examines differences in gene expression profiles, identifies novel snp and determines instances of allelic imbalance. rna-seq analysis in single embryos allows for discrimination of embryo gender and provided the opportunity to characterize individual variability in gene expression. snp analysis of individual samples demonstrates the use of rna-seq to identify embryo-specific variation for association studies and ase that could represent novel layers of developmental regulation subject to influence by ai. our data suggests that ai is prevalent in bovine blastocysts. rna-seq analysis of the individual embryonic transcriptome is feasible and presents valuable insights into gene expression, variation and regulation of the early developmental transcriptome. 
approximately million sequencing reads were generated per embryo and known bovine genes were found to be expressed, with a high correlation of expression levels between samples. expressed biallelic snp variants with allelic imbalances were observed in snp, where one allele represented between 65-95% of a variant’s transcripts. this technology provides the opportunity to assay the developing embryo, but the paucity of biological material available from individual embryos has made this a challenging prospect. we report here the first application of rna-seq for the analysis of individual blastocyst gene expression, snp detection, and characterization of allele specific expression. transcriptomic data was analyzed to identify snp in expressed genes, and individual snp were examined to characterize allele specific expression. rna was extracted from single bovine blastocysts, amplified, and analyzed using high-throughput sequencing. this study represents the first application of rna-seq technology in single bovine embryos allowing a representation of the embryonic transcriptome and the analysis of transcript sequence variation to describe specific allele expression.
 use of rna-seq presents unique benefits in terms of gene expression analysis because of its wide dynamic range and ability to identify functional sequence variants. this technology provides the opportunity to assay the developing embryo, but the paucity of biological material available from individual embryos has made this a challenging prospect.  we report here the first application of rna-seq for the analysis of individual blastocyst gene expression, snp detection, and characterization of allele specific expression. rna was extracted from single bovine blastocysts, amplified, and analyzed using high-throughput sequencing. approximately million sequencing reads were generated per embryo and known bovine genes were found to be expressed, with a high correlation of expression levels between samples. transcriptomic data was analyzed to identify snp in expressed genes, and individual snp were examined to characterize allele specific expression. expressed biallelic snp variants with allelic imbalances were observed in snp, where one allele represented between 65-95% of a variant’s transcripts.  this study represents the first application of rna-seq technology in single bovine embryos allowing a representation of the embryonic transcriptome and the analysis of transcript sequence variation to describe specific allele expression. rna-seqbovineembryoblastocystsnpase 