BACKGROUND
innovations in sequencing technology have allowed deep sequencing of complementary dna , known as ribonucleic acid sequencing , enabling transcriptome assembly and identification of coding and non-coding transcripts across many cell types  <cit> .
transcriptome profiling studies have been undertaken in zebrafish, using polyadenylated  selected messenger rna .
these studies have reported identification of thousands of maternal genes and identified the earliest set of embryonic transcripts.
they also identified a large number of novel transcribed regions in annotated and unannotated regions of the zebrafish genome  <cit> .
in xenopus, several deep-sequencing studies have created different libraries of small rnas from oocytes, eggs, gastrula, liver and skin  <cit> .
a gastrula stage polyadenylated  selected rna-seq profile was used to identify transcribed loci, to enhance gene annotation and to analyze spatial regulation of gene expression  <cit> .
recently, similar polya+ libraries of multiple stages of development were published  <cit> .
for the analysis of transcriptome dynamics it is important to appreciate that, like many other vertebrates, the xenopus maternal-to-zygotic transition involves two important processes: first, destabilization of a subset of maternal mrnas; second, onset of transcription at the mid-blastula transition   <cit> .
studies in xenopus laevis identified distinct phases of maternal, late embryonic and larval gene expression during the course of embryogenesis, whereas microarray analysis in xenopus tropicalis identified several developmentally important maternal mrnas that are regulated by changes in their adenylation during oogenesis and early development  <cit> .
cytoplasmic polyadenylation is essential for the meiotic maturation of the oocyte as it mediates translational activation of mrnas encoding mos kinase and mitotic cyclins involved in early rapid synchronous cell divisions  <cit> .
several maternally polyadenylated mrnas lose their polyadenylated tails after fertilization.
in most cases, this is mediated by an embryonic deadenylation element  in the 3′untranslated region  of the mrna, which binds embryonic deadenylation element - binding protein   <cit> .
processes that regulate mrna deadenylation and degradation are temporally uncoupled.
deadenylated rnas are as stable as their polyadenylated counterparts until the blastula stage, several hours after fertilization  <cit> .
for developmental analysis it is important to establish the dynamics and scale of maternal transcript destabilization on a genome-wide level and to identify the full complement of embryonic transcripts, including as of yet unannotated and long non-coding rnas , the analysis of which will be facilitated by transcriptome profiling using polya+ or total rna-sequencing.
here, we present results from polya+ and ribosomal rna depleted total-rna  sequencing.
our study distinguishes changes in polyadenylation and abundance, which is critical for the the analysis of early transcript dynamics and proper identification of maternal and embryonically induced transcripts.
the embryonic genome shows a gradual cascade of activation, which involves only a third of the number of genes expressed in the oocyte.
by expanding and updating our previously published xenopus tropicalis experimentally validated  annotation pipeline, we also identified  <dig>  new transcripts resulting in a total collection of  <dig>  gene models.
these new transcripts do not overlap with gene models in the xenopus model organism database   <cit> .
using stringent filtering criteria and manual curation,  <dig> transcripts were identified as “stand-alone” lncrna transcripts.
we characterize these transcripts in terms of exon number, transcript length, conservation and expression pattern during embryogenesis and thus anticipate that our catalogue of coding and long non-coding transcripts will enable more developmental and genomic studies directed towards dissecting their functional roles.
RESULTS
deep sequencing of polya+ and total rna libraries
to systematically analyze the transcriptome during early development, we performed polya+ and total rna  sequencing experiments across xenopus tropicalis embryogenesis based on three biological replicates :  oocytes ;  early blastula ;  late blastula ;  gastrula ;  neurula  and  an early larval stage .
we verified abundance of transcripts in biological replicates with real time rt- pcr  using random hexamers.
we not only find minimum variance in transcript abundance among replicates, but also the stage-dependent expression dynamics are similar for the replicates .
total rna was independently extracted for each biological replicate, subjected to polya+ or ribosomal rna depleted total rna  enrichment protocols, quality-controlled, pooled and converted into complementary dna sequencing libraries for the illumina genome analyzer platform .
gene expression was calculated as reads per kilobase of exon model per million mapped reads  and shows a comparable median distribution across sequencing libraries .
heatmap representation of the pearson correlation coefficients reveal the similarity within the early  and the late  transcriptomes respectively .
there are major changes in the pa transcriptome marking meiotic maturation and fertilization  and the maternal-to-zygotic transition .
the total rna  profiles of the early developmental stages correlate relatively well with each other, especially between oocyte and stage  <dig>  most likely due to the presence of stable maternal rnas and the early embryo being transcriptionally quiescent.
correlation between the stages is higher for total rna than the polya+ data, most likely due to changes in polyadenylation of maternal mrna.
to rule out any bias in correlation arising from low expression values, we filtered the data for a threshold of  <dig> rpkm in oocyte .
the pearson correlation heatmap of the filtered data shows a similar profile .
the correlation between same stages in different data sets , while moderate, is highly significant , reflecting the representation of most transcripts in both types of libraries .
a spearman’s rank order correlation analysis strongly underscores the similarities in the total rna and polya+ data .
the data are also in good agreement with previously published ribonucleic acid sequencing  and microarray data .
abundance and polyadenylation state of maternal and maternal-embryonic transcripts
to investigate the maternal contribution to the transcriptome of the developing embryo, we used polya+ and total rna-seq data to classify maternal and embryonic transcripts.
a set of  <dig>  transcripts were called as maternally expressed as they met a filtering criteria of rpkm greater than or equal to  <dig> in oocyte  .
this set of transcripts includes maternal-embryonic genes, which are transcribed in both oocytes and embryos.
our maternal subset also includes well known mouse maternal genes like c-mos, zp <dig>  zp <dig> and slbp   <cit> .
the polya+ and total rna-seq pool of transcripts are likely to have different complexity, and while their rpkm expression levels cannot be compared directly, the ratio of these two measures  may reflect the relative polyadenylation state.
this would allow us to examine the polyadenylated state of all transcripts during early development.
using this strategy we initially compared the pa/rz rpkm ratios for genes like cdk <dig>  kif <dig> and hes <dig> , which are known to be deadenylated after fertilization  <cit> , and indeed, their pa/rz ratios faithfully reflect their fertilization-induced deadenylated state .
genome-wide average pa/rz ratios show an overall abundance of polyadenylated maternal gene products at stage  <dig> relative to earlier stages .
these abundant polya+ transcripts arise either from polyadenylation of the maternally derived message, or new transcription .
we validated the adenylation states of transcripts using rt-qpcr in combination with oligo <dig> primers and random hexamers .
with one exception , polya+ ribonucleic acid sequencing  data and rt-qpcr with oligo <dig> primers correlated very well.
to analyze the deadenylated transcripts more systematically, we filtered the log ratio of two measures  to be less than or equal to - <dig>  for oocyte, stage  <dig> and stage  <dig> respectively.
this filtering gave us sets of genes that are highly enriched in rz data compared to the polya+ data .
interestingly we find stage  <dig> to be highly enriched in deadenylated transcripts.
this may reflect the fertilization-induced deadenylated state, which has been reported in the literature  <cit> .
in order to gain an insight in to the functional categories of deadenylated transcripts at this stage, we compared enrichment of gene ontology  terms of biological processes  using david  <cit> .
we extracted the stage-specific functional annotation charts from david and compared them using clusterprofiler, an r package for comparing gene clusters, with a p-value cut off of  <dig>   <cit> .
interestingly most of the stage  <dig> transcripts show an enrichment of biological processes related to cell-cycle control and regulation .
for example important cell-cycle regulators like aurka, cdc25c and birc <dig>  belong to stage  <dig> rz-enriched fraction.
several chromatin reorganizers and modifiers like ezh <dig>  cbx <dig> and hdac <dig> are also enriched in the stage  <dig> rz- enriched fraction.
to assess patterns of genome-wide polyadenylation during early development in more detail, we used k-means clustering .
clusters  <dig>   <dig>  and  <dig> are groups of maternally abundant polyadenylated transcripts and include well known genes like ccnb <dig>  aurkb, mos, emi <dig> and maskin.
these genes show peak expression in the oocyte and are deadenylated or degraded in a stage-specific manner .
cluster  <dig> represents 12% of all the maternally loaded transcripts which are relatively deadenylated.
this cluster includes histone variants like hist1h2ad, hist1h2al, which are known to exist as deadenylated transcripts  <cit> , emi <dig> which is well studied for its role in unfertilized eggs where it, along with its partner mos causes arrest at metaphase of meiosis ii   <cit> .
cluster  <dig> transcripts are polyadenylated during early development.
a notable gene in this cluster is celf <dig>  which codes for embryonic deadenylation element - binding protein , known to mediate sequence-specific mrna deadenylation  <cit> .
cluster  <dig> represents a group of genes that seem to be loaded as relatively deadenylated messages in the oocyte and are then polyadenylated post-fertilization or post-mbt .
overall 59%  of transcripts are deadenylated during oocyte maturation and early post-fertilization development, whereas 57%  show a higher relative polyadenylation state in late blastulae compared to early blastulae.
motif analysis of 3′ ends of transcripts clustered in figure 2b show a significant enrichment of deadenylation and polyadenylation elements  in several clusters .
to gain insight into the fate and temporal expression patterns of maternally-abundant polyadenylated transcripts after the blastula stage, we compared the polya+ data from six stages  using k-means clustering .
cluster  <dig> includes aurkb, a mitotic serinethreonine kinase, which declines in abundance post-mbt.
we find that genes like tcf <dig> and oct <dig> from cluster  <dig> have different profiles of abundance during development.
oct <dig>  a homologue of the mammalian pluripotency factor oct3/ <dig>  peaks in abundance at late gastrula  and declines drastically thereafter   <cit> .
on the other hand tcf <dig>  a gene encoding a helix-loop-helix transcription factor responsible for mesoderm and axis formation as well as anterior forebrain development via repression of wnt/beta-catenin targets, dramatically peaks at blastula and then exists as a stable polyadenylated transcript up to organogenesis   <cit> .
this analysis shows that the abundance of many maternally loaded polyadenylated transcripts declines after late blastula.
progressive activation of the embryonic genome
an embryonic set of  <dig>  transcripts was obtained by filtering the data for a 10-fold increase in any of the stages compared to oocyte expression.
this set includes transcripts which are transcribed around or immediately after the mid-blastula transition .
concomitant polyadenylation could stabilize the newly synthesized transcripts and mark them for translation.
to examine whether this is the case we compared the relative adenylation status of embryonic transcripts expressed at stage  <dig> to that of maternally abundant stage  <dig> transcripts.
this comparison revealed that a large number of stage  <dig> embryonic transcripts are highly polyadenylated compared to their maternal counterparts, although the polya+ state is also more variable in the embryonic subset of the late blastula transcriptome than in the maternal-embryonic subset .
the polyadenylation distribution of this variable embryonic transcriptome at stage  <dig> is robust and does not change with filtering criteria of embryonic transcripts  .
to explore this variable transcriptomic landscape during embryogenesis, we used k-means clustering to group genes according to their expression  or according to their expression relative to maximum gene expression  .
this revealed that many genes are activated at the mbt, but most are more strongly induced at later stages .
comparing the numbers of expressed transcripts relative to the total embryonic pool of transcripts reveals that about 30% of the embryonic genes are induced by fourfold or more during blastula stages and another 36% by late gastrulation .
clustering gene expression relative to maximum gene expression reveals exquisitely well-defined clusters of dynamic gene expression .
genes which peak in the late blastula  include nodal  <dig> /5/ <dig> and sia1/ <dig>  respectively signalling and transcription factors important for the spemannmangold organizer   <cit> .
gastrula expression represented by cluster  <dig> contains 8% of genes peaking in expression .
this cluster is dominated by high expression of genes involved in mesendoderm specification and patterning such as eomes, chrd, gsc, and t .
cluster  <dig> and  <dig> shows genes peaking during neurulation and include genes like neurog, pax <dig>  nkx2- <dig>  wnt <dig> and myod.
cluster  <dig> and cluster  <dig> together represent clusters of genes peaking in expression at stage  <dig> of development and include transcription factors involved in hematopoietic development like gata <dig> and tal <dig>  celf <dig> is an embryonic gene which belongs to the celf family of genes that code for rna binding proteins involved in deadenylation of mrnas.
it is known to be exclusively expressed in the nervous system including domains in the brain, spinal cord, optic and otic vesicles  <cit> .
irx <dig>  a homeobox transcription factor known for its role in neural patterning also peaks during organogenesis  <cit> .
as the embryonic genome is progressively taking control of development, we found 50% of the developmentally induced genes reach maximum expression late in development .
to correlate these temporal profiles with gene function, enrichment for go terms of bp were examined using david  <cit> .
using gene names as unique identifiers, we found a significant enrichment of stage-specific processes .
we then extracted the stage-specific functional annotation charts from david and compared them using clusterprofiler, an r package for comparing gene clusters, with a p-value cut off of  <dig>    <cit> .
cluster  <dig>  is a small cluster of  <dig> genes and shows enrichment of terms like nucleosome assembly, a term with  <dig> genes  .
the stage 12-specific genes  show enrichment for the term gastrulation and include well-known genes  .
in conclusion, apart from confirming genes enriched in terms with known functions, our go analysis also provides a framework of hypothesis for several genes with unknown function.
experimental validation of gene models and analysis of novel transcripts
to improve gene annotation and identify potentially novel transcripts, we updated our previously published xenopus tropicalis experimentally validated  annotation pipeline  <cit> .
using more sequencing data and the latest genome build, we performed guided transcript assembly with cufflinks using all our polya+ and total rna-seq data with jgi  <dig>  annotation as reference  <cit>  and combined the cufflinks transcripts with expressed sequence tags  clusters .
both histone h <dig> lysine  <dig> tri-methylation  and rna polymerase ii  chromatin immuno-precipitation sequencing  data were used to validate or update the 5′ ends of the gene models as described previously  <cit>  .
the annotation pipeline resulted in a collection of  <dig>  xenopus tropicalis spliced gene models out of which  <dig>  were validated or updated by the xenopus tropicalis experimentally validated  pipeline.
from these validated models,  <dig>   can be detected by rna-seq and 65% have h3k4me <dig> enrichment at the annotated start site.
several thousand gene models were updated and/or reannotated leading to addition of 5′, 3′ or internal exons  known gene model update see additional file 1: figures s <dig> and s7a).
in addition  <dig>  spliced transcripts were newly identified on basis of rna-seq and/or est evidence.
as a by-product of our gene annotation pipeline, we find evidence for a total of  <dig>  single exon unspliced gene models.
these unspliced single exon gene models are filtered out early on in the pipeline and have not been analyzed further .
these single exon models include malat <dig> , a known single exon lncrna conserved in mammals, zebrafish and xenopus <cit> .
from the expression data it appears to be most abundant at the at neurula stage, suggestive of a specialized stage-dependent regulatory role .
compared to our published annotation pipeline, where we only validated and updated known gene models, which are mostly protein-coding, the new implementation is more inclusive in annotating both coding and long non-coding rnas.
to identify new gene models we looked for cufflinks transcripts lacking any overlap with gene models from the xenopus model organism database   <cit> .
we found  <dig>  gene models to be new and non-overlapping , figure 5a, for a complete list with genomic co-ordinates see additional file 3: page –ngm).
out of this set,  <dig> gene models are supported by both expression data  and a 5′ h3k4me <dig> modification peak , see additional file 1: figure s <dig>  for a complete list with genomic co-ordinates co-ordinates see additional file 3: page – ngm-vv), meaning that they are likely stand-alone transcripts.
to validate some of these transcripts, we looked at their relative abundance by performing rt-qpcr of rna from several developmental stages .
we find stage-specific expression .
to assess the coding potential of the various gene model subsets , ngm, ngm-vv) we used open reading frame  length and the log likelihood ratio  of codon bias .
we find that 65% of new transcripts that are not overlapping with known genes , show a codon bias score comparable to non-coding genomic sequences .
also, 70% of these ngm transcripts have a maximum orf length comparable to non-coding genomic sequences .
by contrast, transcripts from annotated protein-coding genes  and all gm have large orf lengths and higher codon bias score, the cumulative distribution of which is well-separated from that of non-coding genomic dna.
however it should be noted that the subset of new non-overlapping transcripts that has strong experimental validation support  may contain both coding and non-coding transcripts.
60% have an orf length of less than  <dig> amino acids and also the codon bias distribution suggests this subset represents a mixed population of coding and non-coding rnas .
to enrich for long non-coding rnas , we curated the ngm-vv subset both bioinformatically and manually.
since high-confidence lncrnas have been shown to lack an orf length greater than  <dig> amino acids  <cit> , we used this as a first step to enrich for putative lncrnas.
inspection of this subset of gene models revealed that it was enriched for the 5′ untranslated region  exons of downstream gene models without h3k4me <dig> peak and a number of other artifacts.
we therefore excluded new gene models that were upstream of the known genes without a h3k4me <dig> peak, or did not meet more stringent expression and splicing criteria .
the resulting  <dig> gene models were screened using blastn and blastp against homology to known protein coding sequences and an array of other problems on the genome browser.
for example models in regions with many gaps were excluded.
also, several intronic transcripts were selected against as they may not represent independent transcription units.
this screening resulted in a set of  <dig> transcripts.
we compared our different subsets of transcripts with the new transcripts identified by tan et al.,  <cit> .
the tan study reports the identification of  <dig>  novel transcripts, a set of transcripts which collapses to a set of  <dig>  unique models that map to the jgi <dig>  genome assembly and are not included in the jgi <dig>  annotation.
only  <dig> of these models overlap with the ngm set of  <dig>  transcripts.
all of the  <dig> transcripts we selected were identified in the tan study, however six of these were linked to protein-coding genes in their multi-stage ribonucleic acid sequencing  data sets.
therefore, our stringent filtering and curation approach has generally enriched for transcripts identified as new gene models in both studies and we conclude that the remaining  <dig> transcript models are likely stand-alone transcription units.
this subset is referred to as ngm-vvo: manually curated new gene models with h3k4me <dig> peaks and rna-seq evidence, and a longest orf length of  <dig> amino acids.
in terms of assessing coding metrics, this set has a codon bias score comparable to random genomic sequences and their cumulative distribution is well separated from known protein-coding mrna .
an example is shown in figure 6b .
compared to the protein-coding transcripts, the manually curated new gene models with orf less than  <dig> amino acids  transcripts have low exon number and relatively shorter transcripts , similar to what has been reported for lncrnas  <cit> .
based on our manual curation, coding potential metric and gene-structure analysis, we conclude that ngm-vvo subset represents a set of high-confidence lncrnas.
mammalian and zebrafish lncrna show a mean expression varying from a 3-fold to a 10-fold difference from their protein-coding counterparts  <cit> .
to investigate the expression of ngm-vvo lncrnas during embryogenesis, we looked at their polya+ mrna profiles during embryogenesis.
we find that the median expression level of lncrnas during embryogenesis is one-third of their protein-coding counterparts .
to investigate expression patterns of ngm-vvo transcripts, we performed unsupervised hierarchical clustering of polya+ and total-rna expression profiles during embryogenesis.
just like the protein-coding genes, some ngm-vvo transcripts have maternal expression.
many are developmentally regulated and show a stage-specific peak in expression .
the gencode v <dig> catalogue of human lncrnas has looked into conservation of human lncrnas using phastcons analysis  <cit> .
our lncrna conservation results are in line with these analyses, since we find ngmvvo exons to be less conserved than annotated proteincoding mrna, but more conserved than the random genomic sequences .
the evolutionary constraints on their sequence and their developmentally regulated expression may be an indication of their stage-specific functionality.
discussion
our results present temporal profiles of maternal and embryonic transcripts during early development.
we report a total number of  <dig>  non-redundant xenbase transcripts expressed in any of the six assayed stages from oocyte to tailbud embryos and mapped to xenopus tropicalis genome assembly   <dig> ).
in our data set the maternal transcriptome consists of over  <dig>  transcripts that show differential adenylation and of these 46% are abundant in the oocyte polya+data .
this is interesting in view of the fact that the oocyte serves as a reservoir of stable maternal transcripts which drive early development in the absence of embryonic transcription.
to better understand the dynamics of polyadenylated vs. deadenylated mrna, we compared the ratio of polya+ and total rna.
this comparison gave us a tool to examine the dynamics of transcript abundance and polyadenylation during early development.
we observed fertilization-induced deadenylation of several cell cycle regulators like cdk <dig>  kif <dig> as already reported  <cit> .
also, it is interesting to note that there is an exclusive pool of relatively deadenylated transcripts, which in our analysis accounts for 12% of the maternal genes and includes well known non-adenylated transcripts like the histone mrnas .
transcription has been reported to start at the mid-blastula stage  <cit> , although a number of genes are transcribed before this stage.
we find little evidence of pre-mbt transcription at stage  <dig> in our data.
many maternal transcripts are gaining polyadenylation during post-fertilization development and may appear as false-positives in an analysis of early transcription if only polya+ messages are considered.
between stage  <dig> and stage  <dig>  some genes are activated early as described for several nodal genes  <cit> .
the embryonic transcript abundance is stage-specific.
about 30% of the embryonic genes are induced four-fold or more by late blastula and another 35% by late gastrulation .
50% of the genes peak late in development .
our go analysis of embryonic genes provides confirmation of genes with known functions as well as provides a framework for hypothesis for several genes with unknown functions .
we have generated an updated annotation pipeline for xenopus tropicalis experimentally validated  gene models, featuring a total of  <dig>  transcripts.
of these, we find  <dig>  gene models to be new, however many of these may be linked to known transcripts and are not independently generated.
the ngm-vvo subset of  <dig> transcripts is a high-confidence set of lncrnas, which shares many of the characteristic features of lncrnas such as low exon number, relatively short length and overall low expression during embryogenesis  <cit> .
they are decorated with h3k4me <dig> histone modification at the 5′ end, evidence that these high-confidence lncrnas are transcribed from their own promoter.
like their protein-coding counterparts, the expression profile of high confidence lncrnas  is stage-specific and temporally restricted.
it proved surprisingly difficult to identify this high-confidence set of lncrnas.
this is because any selection for novel or unannotated transcripts enriches resulting subsets for annotation problems  and assembly problems .
our high-confidence approach may however underestimate the true number of lncrnas that are expressed during embryogenesis.
first of all, lncrna may be transcribed from complex loci and not all may meet our criteria of stand-alone transcripts.
second, many lncrnas are expressed at very low levels.
inclusion of more rna-seq data is therefore likely to identify more lncrnas.
on the other hand, true stand-alone lncrna transcription units, produced from their own promoter, may be far less common than frequently assumed, and the majority of “new transcripts” may arise as by-products of known genes or be transcribed from highly complex loci.
also, rna-seq alignment tools produce artifacts and multiple, sometimes incorrect, models for the same locus.
therefore, approaches that integrate expression and histone modification data are essential to curate transcription units.
functional analyses of the novel lncrnas are required to elucidate their potential roles in pre-mbt transcriptional repression, gastrulation, neurulation and organogenesis.
our catalogue of high confidence stand-alone lncrnas with sequence conservation and stage-specific expression provides a prioritized resource for studies in lncrna function during vertebrate development.
CONCLUSIONS
we provide a comprehensive survey of the xenopus tropicalis transcriptome using polya+ and ribosomal-rna depleted total rna expression data.
these results provide insights into the maternal and embryonic components of expression and polyadenylation dynamics through-out early embryogenesis.
in addition, our improved annotation has led to the discovery of new transcripts which constitute subset of high-confidence stand-alone lncrnas.
together, these data provide a rich developmentally relevant resource, integration of which will enable new genomic and genetic studies in the near future.
materials and methods
animal procedures
x.tropicalis embryos were obtained with in vitro fertilization from three separate crosses .
briefly, both females and males were primed with  <dig> units of human chorionic gonadotropin .
four to six hours before embryo collection, female frogs were injected with  <dig> units of hcg.
forty-five minutes after the onset of egg laying, embryos were collected and dejellied in 3% cysteine hydrochloride  in 10% mmr.
embryos were then cultured in 10% mmr at room temperature and were staged according to nieuwkoop and faber   <cit> .
embryos from three separate clutches were harvested and frozen at -80°c until rna isolation.
stage vi oocytes were harvested by treating ovarian follicles with collagenase .
rna preparation and sequencing
oocyte and embryos from nieuwkoop-faber stages  <dig> -  <dig> were collected and total rna was isolated using trizol and the qiagen rneasy kit.
subsequently, polyadenylated rna was selected by enriching with the oligotex mrna kit .
to ensure complete removal of ribosomal rna , polyadenylated mrna was subjected to an additional round of oligotex treatment.
total rna was subjected to depletion of ribosomal rna  using ribozero epicenter low input kit.
two important quality control measures were taken to confirm removal of ribosomal rna .
first, the ribosomal rna  -depleted sample was tested on a rna-chip  in comparison with non-depleted total rna.
absence of 28s and 18s peaks in the ribosomal rna  depleted sample confirmed good depletion.
second, rt-qpcr with primers against 28s, 5s and gapdh was performed.
28s rna levels were less than 5%, typically around 1% after depletion, whereas gapdh was typically at more than 80% of the levels before depletion.
for sequencing, cdna was prepared for both polya+ and rz samples with random hexamer primers using superscript iii  and the second strand was made with dna polymerase i, dna ligase and t <dig> dna polymerase.
the purified double-stranded cdna was used for illumina sample preparation.
all quality control qpcr reactions were performed on a mylq single-color reader real-time pcr detection system  using iq sybr green supermix .
the three biological replicates were checked by rt-qpcr and pooled for sample preparation and sequencing.
these samples were then processed according to the manufacturer’s protocol .
shortly, adapter sequences were linked to the complementary dna  samples, the library was size selected , and amplified by polymerase chain reaction .
the subsequent sequencing was carried out on genome analyzer .
rna-seq expression analysis
on average, we obtained about 16- <dig> million reads per stage .
out of the total reads about 50-60% could be aligned to the genome assembly  of the xenopus tropicalis genome sequence.
to allow a quantitative comparison all reads were normalized before analysis.
the transcript list contains all the genes that are expressed  in at least one stage.
the rpkm per gene is the mean of all rpkm of all the non-redundant exons of all isoforms per gene.
the total list contains around  <dig>  genes of which only  <dig> are not detected as expressed in any stage.
all unknown/unnamed gene names have been changed to include the genomic position for reasons of identification.
alignment was performed using burrows-wheeler aligner , reads mapping to multiple positions  were not included in the rpkm calculation  <cit> .
xtev  gene annotation pipeline
gene models  were downloaded from xenbase  and est clusters mapped to the jgi  <dig>  x.tropicalis genome were supplied by mike gilchrist .
spliced transcription units were generated from the rna-seq data.
all reads were mapped to jgi  <dig>  using tophat v <dig> .
<dig>  <cit> .
the tophat output was filtered to keep only new splice sites with evidence of at least  <dig> spliced reads.
the filtered tophat output was used with cufflinks v <dig> .
<dig> to perform transcript assembly  <cit> .
the experimental annotation pipeline consists of several steps: 1) collect gene models; 2) update with experimental data; 3) validate and/or update gene models with rna-seq data; 4) validate and/or update transcription start-site  with h3k4me <dig> and/or rnapii chip-seq data and 5) choose the most likely model .
all xenbase gene models sharing at least  <dig> exon were considered as multiple models of a single gene.
the est clusters and transcripts determined by cufflinks were used to update the gene models with extra putative exons, mainly at the 5′ and 3′ end of genes.
the number of rna-seq reads was determined for all exons of all models.
if 1/ <dig> of the exons of a model contained at least  <dig> rna-seq reads the model was considered as expressed.
if the first exon of a gene model overlapped with a h3k4me <dig> peak, the tss was considered as validated.
if no single model of a gene had a validated tss, we looked for evidence of a tss upstream.
in this case there had to be a h3k4me <dig> peak upstream of a gene model, with no different gene models in between, and the mean rnapii level of the region between the upstream h3k4me <dig> peak and the 5′ exon of the gene had to be at least  <dig>  of the mean rnapii level of the gene body.
furthermore, all gene models were checked for evidence of a downstream h3k4me <dig> peak, which can indicate a putative alternative tss.
for each single gene the most likely model was chosen according to the following criteria, in order of decreasing importance: a validated tss, number of expressed exons, number of exons.
of the new models, which are not present in the xenbase jgi  <dig>  annotation, only spliced transcripts were included.
analysis of coding potential
coding potential of rna sequences was determined using maximal orf length and codon bias metrics, as described  <cit> .
the codon bias metric is based on unequal usage of synonymous codons.
briefly, triplet frequencies were determined in non-coding genomic x.tropicalis dna , whereas x.tropicalis codon frequencies were downloaded from http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species= <dig>  log likelihood ratios  for each codon were calculated based on the codon frequency conditional of the encoded amino acid, such that for each codon i coding for amino acid ai, llri = log <dig>  in which ci and ni correspond to the likelihood of codon i conditional on amino acid ai in coding and non-coding sequences respectively .
the total llr score is determined by summing llri values in all  <dig> bp windows in six potential reading frames.
after computing a score for windows, the max llr score was defined as the maximum score observed in all windows of the transcript.
quantitative rt-qpcr for known and new gene model  validation
validation of known and novel transcripts was performed on total-rna which was subjected to depletion of ribosomal rna  using ribozero epicenter low input kit.
total rna was then dnase treated and column purified to remove any contaminating genomic dna.
cdna was prepared with oligo <dig> primers or random hexamers using superscript iii .
qpcr reactions were performed on a mylq single-color reader real-time pcr detection system  using iq sybr green supermix .
with stage  <dig> as reference, fold change was calculated by normalizing ct values in stages  <dig>   <dig>  and  <dig> against odc <dig> gene using the 2−ΔΔct method   <cit> .
bioinformatic and manual curation of ngm-vv subset
ngm-vv subset is collection of  <dig> new gene models.
as a first step, we filtered these models for orf length less than or equal to  <dig> amino acids.
this resulted in a set of  <dig> gene models, which were then screened using following criteria: 1) absence of a downstream gene  with a x or u annotation for h3k4me <dig> ; 2) rpkm of all exons should be greater than or equal  <dig>  and 3) evidence of splicing in our data.
this resulted in a set of  <dig> gene models.
these models were then manually curated using blastn and blastp to filter against homology to known protein coding sequences .
conservation  analysis
for conservation analysis all gene models were mapped to jgi v <dig>   using blat.
the average phastcons score per gene model was calculated using the conservation track  of the ucsc genome browser.
data availability
the data have been deposited in ncbi′s gene expression omnibus  <cit>  and are accessible through geo series accession number gse <dig>  xtev gene models are available at: http://veenstra.ncmls.nl/genomedata.asp.
statement on animal use
animal care and use was in accordance with national and european guidelines and standard operating procedures approved by the institutional animal care and use committee .
abbreviations
polya+: polyadenylated; rt-qpcr: real time rt-pcr; cdna: complementary dna; rna-seq: ribonucleic acid sequencing; mrna: messenger rna; mbt: mid-blastula transition; eden: embryonic deadenylation element; eden-bp: embryonic deadenylation element - binding protein; utr: untranslated] region; xtev: xenopus tropicalis experimentally validated; lncrnas: long non-coding rnas; orf: open reading frame; pa: polya+; ribozero, rz: ribosomal rna depleted total-rna; rpkm, see materials and methods: reads per kilobase of exon model per million mapped reads; go: gene ontology; bp: biological processes; est: expressed sequence tags; gm: gene models; ngm: new gene models; ngm-vv: new gene models with validation support; ngm-vvo: manually curated new gene models with orf less than  <dig> amino acids; llr: log likelihood ratio; jgi: joint genome institute; rrna: ribosomal rna; bwa: burrows-wheeler aligner; pcr: polymerase chain reaction; h3k4me3: histone h <dig> lysine  <dig> tri-methylation; rz: ribosomal rna depleted total rna; rnapii: rna polymerase ii; chip-seq: chromatin immuno-precipitation sequencing; tss: transcription start-site.
competing interests
the authors declare that they have no competing interests.
authors’ contributions
ssp, ugj and gjcv conceived and designed the study, ugj performed the rna-seq experiments, ssp did the analysis and performed the validation experiments.
sjvh performed transcriptome assembly, designed and assembled the xenopus tropicalis experimentally validated  pipeline and performed the conservation analysis.
gjcv performed orf length and codon bias analysis.
ssp and gjcv wrote the manuscript.
all authors read and approved the final manuscript.
supplementary material
additional file 1
supplementary tables  and figures .
click here for file
additional file 2
multi-page list of genes expressed in
xenopus tropicalis
.
click here for file
additional file 3
multi-page xtev annotation pipeline gene list with genomic co-ordinates.
click here for file
additional file 4
single-page list of codon usage frequencies used for calculating codon bias score.
click here for file
additional file 5
single-page list of primers used for real-time qpcr validation of highly conserved lncrnas from the ngm subset.
click here for file
acknowledgements
we thank dr.mike gilchrist for providing the est data .
this work was supported by grants of the national institutes of health  and nwo, the netherlands organization of scientific research .

 innovations in sequencing technology have allowed deep sequencing of complementary dna, known as ribonucleic acid sequencing, enabling transcriptome assembly and identification of coding and non-coding transcripts across many cell types. transcriptome profiling studies have been undertaken in zebrafish, using polyadenylated selected messenger rna. these studies have reported identification of thousands of maternal genes and identified the earliest set of embryonic transcripts. they also identified a large number of novel transcribed regions in annotated and unannotated regions of the zebrafish genome. in xenopus, several deep-sequencing studies have created different libraries of small rnas from oocytes, eggs, gastrula, liver and skin. a gastrula stage polyadenylated selected rna-seq profile was used to identify transcribed loci, to enhance gene annotation and to analyze spatial regulation of gene expression. recently, similar polya+ libraries of multiple stages of development were published. for the analysis of transcriptome dynamics it is important to appreciate that, like many other vertebrates, the xenopus maternal-to-zygotic transition involves two important processes: first, destabilization of a subset of maternal mrnas; second, onset of transcription at the mid-blastula transition. studies in xenopus laevis identified distinct phases of maternal, late embryonic and larval gene expression during the course of embryogenesis, whereas microarray analysis in xenopus tropicalis identified several developmentally important maternal mrnas that are regulated by changes in their adenylation during oogenesis and early development. cytoplasmic polyadenylation is essential for the meiotic maturation of the oocyte as it mediates translational activation of mrnas encoding mos kinase and mitotic cyclins involved in early rapid synchronous cell divisions. several maternally polyadenylated mrnas lose their polyadenylated tails after fertilization. in most cases, this is mediated by an embryonic deadenylation element in the 3′untranslated region of the mrna, which binds embryonic deadenylation element - binding protein. processes that regulate mrna deadenylation and degradation are temporally uncoupled. deadenylated rnas are as stable as their polyadenylated counterparts until the blastula stage, several hours after fertilization. for developmental analysis it is important to establish the dynamics and scale of maternal transcript destabilization on a genome-wide level and to identify the full complement of embryonic transcripts, including as of yet unannotated and long non-coding rnas, the analysis of which will be facilitated by transcriptome profiling using polya+ or total rna-sequencing. here, we present results from polya+ and ribosomal rna depleted total-rna sequencing. our study distinguishes changes in polyadenylation and abundance, which is critical for the the analysis of early transcript dynamics and proper identification of maternal and embryonically induced transcripts. the embryonic genome shows a gradual cascade of activation, which involves only a third of the number of genes expressed in the oocyte. by expanding and updating our previously published xenopus tropicalis experimentally validated annotation pipeline, we also identified new transcripts resulting in a total collection of gene models. these new transcripts do not overlap with gene models in the xenopus model organism database. using stringent filtering criteria and manual curation, transcripts were identified as “stand-alone” lncrna transcripts. we characterize these transcripts in terms of exon number, transcript length, conservation and expression pattern during embryogenesis and thus anticipate that our catalogue of coding and long non-coding transcripts will enable more developmental and genomic studies directed towards dissecting their functional roles.  deep sequencing of polya+ and total rna libraries to systematically analyze the transcriptome during early development, we performed polya+ and total rna sequencing experiments across xenopus tropicalis embryogenesis based on three biological replicates: oocytes; early blastula; late blastula; gastrula; neurula and an early larval stage. we verified abundance of transcripts in biological replicates with real time rt- pcr using random hexamers. we not only find minimum variance in transcript abundance among replicates, but also the stage-dependent expression dynamics are similar for the replicates. total rna was independently extracted for each biological replicate, subjected to polya+ or ribosomal rna depleted total rna enrichment protocols, quality-controlled, pooled and converted into complementary dna sequencing libraries for the illumina genome analyzer platform. gene expression was calculated as reads per kilobase of exon model per million mapped reads and shows a comparable median distribution across sequencing libraries. heatmap representation of the pearson correlation coefficients reveal the similarity within the early and the late transcriptomes respectively. there are major changes in the pa transcriptome marking meiotic maturation and fertilization and the maternal-to-zygotic transition. the total rna profiles of the early developmental stages correlate relatively well with each other, especially between oocyte and stage most likely due to the presence of stable maternal rnas and the early embryo being transcriptionally quiescent. correlation between the stages is higher for total rna than the polya+ data, most likely due to changes in polyadenylation of maternal mrna. to rule out any bias in correlation arising from low expression values, we filtered the data for a threshold of rpkm in oocyte. the pearson correlation heatmap of the filtered data shows a similar profile. the correlation between same stages in different data sets, while moderate, is highly significant, reflecting the representation of most transcripts in both types of libraries. a spearman’s rank order correlation analysis strongly underscores the similarities in the total rna and polya+ data. the data are also in good agreement with previously published ribonucleic acid sequencing and microarray data. abundance and polyadenylation state of maternal and maternal-embryonic transcripts to investigate the maternal contribution to the transcriptome of the developing embryo, we used polya+ and total rna-seq data to classify maternal and embryonic transcripts. a set of transcripts were called as maternally expressed as they met a filtering criteria of rpkm greater than or equal to in oocyte. this set of transcripts includes maternal-embryonic genes, which are transcribed in both oocytes and embryos. our maternal subset also includes well known mouse maternal genes like c-mos, zp zp and slbp. the polya+ and total rna-seq pool of transcripts are likely to have different complexity, and while their rpkm expression levels cannot be compared directly, the ratio of these two measures may reflect the relative polyadenylation state. this would allow us to examine the polyadenylated state of all transcripts during early development. using this strategy we initially compared the pa/rz rpkm ratios for genes like cdk kif and hes, which are known to be deadenylated after fertilization, and indeed, their pa/rz ratios faithfully reflect their fertilization-induced deadenylated state. genome-wide average pa/rz ratios show an overall abundance of polyadenylated maternal gene products at stage relative to earlier stages. these abundant polya+ transcripts arise either from polyadenylation of the maternally derived message, or new transcription. we validated the adenylation states of transcripts using rt-qpcr in combination with oligo primers and random hexamers. with one exception, polya+ ribonucleic acid sequencing data and rt-qpcr with oligo primers correlated very well. to analyze the deadenylated transcripts more systematically, we filtered the log ratio of two measures to be less than or equal to - for oocyte, stage and stage respectively. this filtering gave us sets of genes that are highly enriched in rz data compared to the polya+ data. interestingly we find stage to be highly enriched in deadenylated transcripts. this may reflect the fertilization-induced deadenylated state, which has been reported in the literature. in order to gain an insight in to the functional categories of deadenylated transcripts at this stage, we compared enrichment of gene ontology terms of biological processes using david. we extracted the stage-specific functional annotation charts from david and compared them using clusterprofiler, an r package for comparing gene clusters, with a p-value cut off of. interestingly most of the stage transcripts show an enrichment of biological processes related to cell-cycle control and regulation. for example important cell-cycle regulators like aurka, cdc25c and birc belong to stage rz-enriched fraction. several chromatin reorganizers and modifiers like ezh cbx and hdac are also enriched in the stage rz- enriched fraction. to assess patterns of genome-wide polyadenylation during early development in more detail, we used k-means clustering. clusters and are groups of maternally abundant polyadenylated transcripts and include well known genes like ccnb aurkb, mos, emi and maskin. these genes show peak expression in the oocyte and are deadenylated or degraded in a stage-specific manner. cluster represents 12% of all the maternally loaded transcripts which are relatively deadenylated. this cluster includes histone variants like hist1h2ad, hist1h2al, which are known to exist as deadenylated transcripts, emi which is well studied for its role in unfertilized eggs where it, along with its partner mos causes arrest at metaphase of meiosis ii. cluster transcripts are polyadenylated during early development. a notable gene in this cluster is celf which codes for embryonic deadenylation element - binding protein, known to mediate sequence-specific mrna deadenylation. cluster represents a group of genes that seem to be loaded as relatively deadenylated messages in the oocyte and are then polyadenylated post-fertilization or post-mbt. overall 59% of transcripts are deadenylated during oocyte maturation and early post-fertilization development, whereas 57% show a higher relative polyadenylation state in late blastulae compared to early blastulae. motif analysis of 3′ ends of transcripts clustered in figure 2b show a significant enrichment of deadenylation and polyadenylation elements in several clusters. to gain insight into the fate and temporal expression patterns of maternally-abundant polyadenylated transcripts after the blastula stage, we compared the polya+ data from six stages using k-means clustering. cluster includes aurkb, a mitotic serinethreonine kinase, which declines in abundance post-mbt. we find that genes like tcf and oct from cluster have different profiles of abundance during development. oct a homologue of the mammalian pluripotency factor oct3/ peaks in abundance at late gastrula and declines drastically thereafter. on the other hand tcf a gene encoding a helix-loop-helix transcription factor responsible for mesoderm and axis formation as well as anterior forebrain development via repression of wnt/beta-catenin targets, dramatically peaks at blastula and then exists as a stable polyadenylated transcript up to organogenesis. this analysis shows that the abundance of many maternally loaded polyadenylated transcripts declines after late blastula. progressive activation of the embryonic genome an embryonic set of transcripts was obtained by filtering the data for a 10-fold increase in any of the stages compared to oocyte expression. this set includes transcripts which are transcribed around or immediately after the mid-blastula transition. concomitant polyadenylation could stabilize the newly synthesized transcripts and mark them for translation. to examine whether this is the case we compared the relative adenylation status of embryonic transcripts expressed at stage to that of maternally abundant stage transcripts. this comparison revealed that a large number of stage embryonic transcripts are highly polyadenylated compared to their maternal counterparts, although the polya+ state is also more variable in the embryonic subset of the late blastula transcriptome than in the maternal-embryonic subset. the polyadenylation distribution of this variable embryonic transcriptome at stage is robust and does not change with filtering criteria of embryonic transcripts. to explore this variable transcriptomic landscape during embryogenesis, we used k-means clustering to group genes according to their expression or according to their expression relative to maximum gene expression. this revealed that many genes are activated at the mbt, but most are more strongly induced at later stages. comparing the numbers of expressed transcripts relative to the total embryonic pool of transcripts reveals that about 30% of the embryonic genes are induced by fourfold or more during blastula stages and another 36% by late gastrulation. clustering gene expression relative to maximum gene expression reveals exquisitely well-defined clusters of dynamic gene expression. genes which peak in the late blastula include nodal /5/ and sia1/ respectively signalling and transcription factors important for the spemannmangold organizer. gastrula expression represented by cluster contains 8% of genes peaking in expression. this cluster is dominated by high expression of genes involved in mesendoderm specification and patterning such as eomes, chrd, gsc, and t. cluster and shows genes peaking during neurulation and include genes like neurog, pax nkx2- wnt and myod. cluster and cluster together represent clusters of genes peaking in expression at stage of development and include transcription factors involved in hematopoietic development like gata and tal celf is an embryonic gene which belongs to the celf family of genes that code for rna binding proteins involved in deadenylation of mrnas. it is known to be exclusively expressed in the nervous system including domains in the brain, spinal cord, optic and otic vesicles. irx a homeobox transcription factor known for its role in neural patterning also peaks during organogenesis. as the embryonic genome is progressively taking control of development, we found 50% of the developmentally induced genes reach maximum expression late in development. to correlate these temporal profiles with gene function, enrichment for go terms of bp were examined using david. using gene names as unique identifiers, we found a significant enrichment of stage-specific processes. we then extracted the stage-specific functional annotation charts from david and compared them using clusterprofiler, an r package for comparing gene clusters, with a p-value cut off of. cluster is a small cluster of genes and shows enrichment of terms like nucleosome assembly, a term with genes. the stage 12-specific genes show enrichment for the term gastrulation and include well-known genes. in conclusion, apart from confirming genes enriched in terms with known functions, our go analysis also provides a framework of hypothesis for several genes with unknown function. experimental validation of gene models and analysis of novel transcripts to improve gene annotation and identify potentially novel transcripts, we updated our previously published xenopus tropicalis experimentally validated annotation pipeline. using more sequencing data and the latest genome build, we performed guided transcript assembly with cufflinks using all our polya+ and total rna-seq data with jgi annotation as reference and combined the cufflinks transcripts with expressed sequence tags clusters. both histone h lysine tri-methylation and rna polymerase ii chromatin immuno-precipitation sequencing data were used to validate or update the 5′ ends of the gene models as described previously. the annotation pipeline resulted in a collection of xenopus tropicalis spliced gene models out of which were validated or updated by the xenopus tropicalis experimentally validated pipeline. from these validated models, can be detected by rna-seq and 65% have h3k4me enrichment at the annotated start site. several thousand gene models were updated and/or reannotated leading to addition of 5′, 3′ or internal exons known gene model update see additional file 1: figures s and s7a). in addition spliced transcripts were newly identified on basis of rna-seq and/or est evidence. as a by-product of our gene annotation pipeline, we find evidence for a total of single exon unspliced gene models. these unspliced single exon gene models are filtered out early on in the pipeline and have not been analyzed further. these single exon models include malat, a known single exon lncrna conserved in mammals, zebrafish and xenopus. from the expression data it appears to be most abundant at the at neurula stage, suggestive of a specialized stage-dependent regulatory role. compared to our published annotation pipeline, where we only validated and updated known gene models, which are mostly protein-coding, the new implementation is more inclusive in annotating both coding and long non-coding rnas. to identify new gene models we looked for cufflinks transcripts lacking any overlap with gene models from the xenopus model organism database. we found gene models to be new and non-overlapping, figure 5a, for a complete list with genomic co-ordinates see additional file 3: page –ngm). out of this set, gene models are supported by both expression data and a 5′ h3k4me modification peak, see additional file 1: figure s for a complete list with genomic co-ordinates co-ordinates see additional file 3: page – ngm-vv), meaning that they are likely stand-alone transcripts. to validate some of these transcripts, we looked at their relative abundance by performing rt-qpcr of rna from several developmental stages. we find stage-specific expression. to assess the coding potential of the various gene model subsets, ngm, ngm-vv) we used open reading frame length and the log likelihood ratio of codon bias. we find that 65% of new transcripts that are not overlapping with known genes, show a codon bias score comparable to non-coding genomic sequences. also, 70% of these ngm transcripts have a maximum orf length comparable to non-coding genomic sequences. by contrast, transcripts from annotated protein-coding genes and all gm have large orf lengths and higher codon bias score, the cumulative distribution of which is well-separated from that of non-coding genomic dna. however it should be noted that the subset of new non-overlapping transcripts that has strong experimental validation support may contain both coding and non-coding transcripts. 60% have an orf length of less than amino acids and also the codon bias distribution suggests this subset represents a mixed population of coding and non-coding rnas. to enrich for long non-coding rnas, we curated the ngm-vv subset both bioinformatically and manually. since high-confidence lncrnas have been shown to lack an orf length greater than amino acids, we used this as a first step to enrich for putative lncrnas. inspection of this subset of gene models revealed that it was enriched for the 5′ untranslated region exons of downstream gene models without h3k4me peak and a number of other artifacts. we therefore excluded new gene models that were upstream of the known genes without a h3k4me peak, or did not meet more stringent expression and splicing criteria. the resulting gene models were screened using blastn and blastp against homology to known protein coding sequences and an array of other problems on the genome browser. for example models in regions with many gaps were excluded. also, several intronic transcripts were selected against as they may not represent independent transcription units. this screening resulted in a set of transcripts. we compared our different subsets of transcripts with the new transcripts identified by tan et al.,. the tan study reports the identification of novel transcripts, a set of transcripts which collapses to a set of unique models that map to the jgi genome assembly and are not included in the jgi annotation. only of these models overlap with the ngm set of transcripts. all of the transcripts we selected were identified in the tan study, however six of these were linked to protein-coding genes in their multi-stage ribonucleic acid sequencing data sets. therefore, our stringent filtering and curation approach has generally enriched for transcripts identified as new gene models in both studies and we conclude that the remaining transcript models are likely stand-alone transcription units. this subset is referred to as ngm-vvo: manually curated new gene models with h3k4me peaks and rna-seq evidence, and a longest orf length of amino acids. in terms of assessing coding metrics, this set has a codon bias score comparable to random genomic sequences and their cumulative distribution is well separated from known protein-coding mrna. an example is shown in figure 6b. compared to the protein-coding transcripts, the manually curated new gene models with orf less than amino acids transcripts have low exon number and relatively shorter transcripts, similar to what has been reported for lncrnas. based on our manual curation, coding potential metric and gene-structure analysis, we conclude that ngm-vvo subset represents a set of high-confidence lncrnas. mammalian and zebrafish lncrna show a mean expression varying from a 3-fold to a 10-fold difference from their protein-coding counterparts. to investigate the expression of ngm-vvo lncrnas during embryogenesis, we looked at their polya+ mrna profiles during embryogenesis. we find that the median expression level of lncrnas during embryogenesis is one-third of their protein-coding counterparts. to investigate expression patterns of ngm-vvo transcripts, we performed unsupervised hierarchical clustering of polya+ and total-rna expression profiles during embryogenesis. just like the protein-coding genes, some ngm-vvo transcripts have maternal expression. many are developmentally regulated and show a stage-specific peak in expression. the gencode v catalogue of human lncrnas has looked into conservation of human lncrnas using phastcons analysis. our lncrna conservation results are in line with these analyses, since we find ngmvvo exons to be less conserved than annotated proteincoding mrna, but more conserved than the random genomic sequences. the evolutionary constraints on their sequence and their developmentally regulated expression may be an indication of their stage-specific functionality. discussion our results present temporal profiles of maternal and embryonic transcripts during early development. we report a total number of non-redundant xenbase transcripts expressed in any of the six assayed stages from oocyte to tailbud embryos and mapped to xenopus tropicalis genome assembly ). in our data set the maternal transcriptome consists of over transcripts that show differential adenylation and of these 46% are abundant in the oocyte polya+data. this is interesting in view of the fact that the oocyte serves as a reservoir of stable maternal transcripts which drive early development in the absence of embryonic transcription. to better understand the dynamics of polyadenylated vs. deadenylated mrna, we compared the ratio of polya+ and total rna. this comparison gave us a tool to examine the dynamics of transcript abundance and polyadenylation during early development. we observed fertilization-induced deadenylation of several cell cycle regulators like cdk kif as already reported. also, it is interesting to note that there is an exclusive pool of relatively deadenylated transcripts, which in our analysis accounts for 12% of the maternal genes and includes well known non-adenylated transcripts like the histone mrnas. transcription has been reported to start at the mid-blastula stage, although a number of genes are transcribed before this stage. we find little evidence of pre-mbt transcription at stage in our data. many maternal transcripts are gaining polyadenylation during post-fertilization development and may appear as false-positives in an analysis of early transcription if only polya+ messages are considered. between stage and stage some genes are activated early as described for several nodal genes. the embryonic transcript abundance is stage-specific. about 30% of the embryonic genes are induced four-fold or more by late blastula and another 35% by late gastrulation. 50% of the genes peak late in development. our go analysis of embryonic genes provides confirmation of genes with known functions as well as provides a framework for hypothesis for several genes with unknown functions. we have generated an updated annotation pipeline for xenopus tropicalis experimentally validated gene models, featuring a total of transcripts. of these, we find gene models to be new, however many of these may be linked to known transcripts and are not independently generated. the ngm-vvo subset of transcripts is a high-confidence set of lncrnas, which shares many of the characteristic features of lncrnas such as low exon number, relatively short length and overall low expression during embryogenesis. they are decorated with h3k4me histone modification at the 5′ end, evidence that these high-confidence lncrnas are transcribed from their own promoter. like their protein-coding counterparts, the expression profile of high confidence lncrnas is stage-specific and temporally restricted. it proved surprisingly difficult to identify this high-confidence set of lncrnas. this is because any selection for novel or unannotated transcripts enriches resulting subsets for annotation problems and assembly problems. our high-confidence approach may however underestimate the true number of lncrnas that are expressed during embryogenesis. first of all, lncrna may be transcribed from complex loci and not all may meet our criteria of stand-alone transcripts. second, many lncrnas are expressed at very low levels. inclusion of more rna-seq data is therefore likely to identify more lncrnas. on the other hand, true stand-alone lncrna transcription units, produced from their own promoter, may be far less common than frequently assumed, and the majority of “new transcripts” may arise as by-products of known genes or be transcribed from highly complex loci. also, rna-seq alignment tools produce artifacts and multiple, sometimes incorrect, models for the same locus. therefore, approaches that integrate expression and histone modification data are essential to curate transcription units. functional analyses of the novel lncrnas are required to elucidate their potential roles in pre-mbt transcriptional repression, gastrulation, neurulation and organogenesis. our catalogue of high confidence stand-alone lncrnas with sequence conservation and stage-specific expression provides a prioritized resource for studies in lncrna function during vertebrate development.  we provide a comprehensive survey of the xenopus tropicalis transcriptome using polya+ and ribosomal-rna depleted total rna expression data. these results provide insights into the maternal and embryonic components of expression and polyadenylation dynamics through-out early embryogenesis. in addition, our improved annotation has led to the discovery of new transcripts which constitute subset of high-confidence stand-alone lncrnas. together, these data provide a rich developmentally relevant resource, integration of which will enable new genomic and genetic studies in the near future. materials and methods animal procedures x.tropicalis embryos were obtained with in vitro fertilization from three separate crosses. briefly, both females and males were primed with units of human chorionic gonadotropin. four to six hours before embryo collection, female frogs were injected with units of hcg. forty-five minutes after the onset of egg laying, embryos were collected and dejellied in 3% cysteine hydrochloride in 10% mmr. embryos were then cultured in 10% mmr at room temperature and were staged according to nieuwkoop and faber. embryos from three separate clutches were harvested and frozen at -80°c until rna isolation. stage vi oocytes were harvested by treating ovarian follicles with collagenase. rna preparation and sequencing oocyte and embryos from nieuwkoop-faber stages - were collected and total rna was isolated using trizol and the qiagen rneasy kit. subsequently, polyadenylated rna was selected by enriching with the oligotex mrna kit. to ensure complete removal of ribosomal rna, polyadenylated mrna was subjected to an additional round of oligotex treatment. total rna was subjected to depletion of ribosomal rna using ribozero epicenter low input kit. two important quality control measures were taken to confirm removal of ribosomal rna. first, the ribosomal rna -depleted sample was tested on a rna-chip in comparison with non-depleted total rna. absence of 28s and 18s peaks in the ribosomal rna depleted sample confirmed good depletion. second, rt-qpcr with primers against 28s, 5s and gapdh was performed. 28s rna levels were less than 5%, typically around 1% after depletion, whereas gapdh was typically at more than 80% of the levels before depletion. for sequencing, cdna was prepared for both polya+ and rz samples with random hexamer primers using superscript iii and the second strand was made with dna polymerase i, dna ligase and t dna polymerase. the purified double-stranded cdna was used for illumina sample preparation. all quality control qpcr reactions were performed on a mylq single-color reader real-time pcr detection system using iq sybr green supermix. the three biological replicates were checked by rt-qpcr and pooled for sample preparation and sequencing. these samples were then processed according to the manufacturer’s protocol. shortly, adapter sequences were linked to the complementary dna samples, the library was size selected, and amplified by polymerase chain reaction. the subsequent sequencing was carried out on genome analyzer. rna-seq expression analysis on average, we obtained about 16- million reads per stage. out of the total reads about 50-60% could be aligned to the genome assembly of the xenopus tropicalis genome sequence. to allow a quantitative comparison all reads were normalized before analysis. the transcript list contains all the genes that are expressed in at least one stage. the rpkm per gene is the mean of all rpkm of all the non-redundant exons of all isoforms per gene. the total list contains around genes of which only are not detected as expressed in any stage. all unknown/unnamed gene names have been changed to include the genomic position for reasons of identification. alignment was performed using burrows-wheeler aligner, reads mapping to multiple positions were not included in the rpkm calculation. xtev gene annotation pipeline gene models were downloaded from xenbase and est clusters mapped to the jgi x.tropicalis genome were supplied by mike gilchrist. spliced transcription units were generated from the rna-seq data. all reads were mapped to jgi using tophat v.. the tophat output was filtered to keep only new splice sites with evidence of at least spliced reads. the filtered tophat output was used with cufflinks v. to perform transcript assembly. the experimental annotation pipeline consists of several steps: 1) collect gene models; 2) update with experimental data; 3) validate and/or update gene models with rna-seq data; 4) validate and/or update transcription start-site with h3k4me and/or rnapii chip-seq data and 5) choose the most likely model. all xenbase gene models sharing at least exon were considered as multiple models of a single gene. the est clusters and transcripts determined by cufflinks were used to update the gene models with extra putative exons, mainly at the 5′ and 3′ end of genes. the number of rna-seq reads was determined for all exons of all models. if 1/ of the exons of a model contained at least rna-seq reads the model was considered as expressed. if the first exon of a gene model overlapped with a h3k4me peak, the tss was considered as validated. if no single model of a gene had a validated tss, we looked for evidence of a tss upstream. in this case there had to be a h3k4me peak upstream of a gene model, with no different gene models in between, and the mean rnapii level of the region between the upstream h3k4me peak and the 5′ exon of the gene had to be at least of the mean rnapii level of the gene body. furthermore, all gene models were checked for evidence of a downstream h3k4me peak, which can indicate a putative alternative tss. for each single gene the most likely model was chosen according to the following criteria, in order of decreasing importance: a validated tss, number of expressed exons, number of exons. of the new models, which are not present in the xenbase jgi annotation, only spliced transcripts were included. analysis of coding potential coding potential of rna sequences was determined using maximal orf length and codon bias metrics, as described. the codon bias metric is based on unequal usage of synonymous codons. briefly, triplet frequencies were determined in non-coding genomic x.tropicalis dna, whereas x.tropicalis codon frequencies were downloaded from http://www.kazusa.or.jp/codon/cgi-bin/showcodon.cgi?species= log likelihood ratios for each codon were calculated based on the codon frequency conditional of the encoded amino acid, such that for each codon i coding for amino acid ai, llri = log in which ci and ni correspond to the likelihood of codon i conditional on amino acid ai in coding and non-coding sequences respectively. the total llr score is determined by summing llri values in all bp windows in six potential reading frames. after computing a score for windows, the max llr score was defined as the maximum score observed in all windows of the transcript. quantitative rt-qpcr for known and new gene model validation validation of known and novel transcripts was performed on total-rna which was subjected to depletion of ribosomal rna using ribozero epicenter low input kit. total rna was then dnase treated and column purified to remove any contaminating genomic dna. cdna was prepared with oligo primers or random hexamers using superscript iii. qpcr reactions were performed on a mylq single-color reader real-time pcr detection system using iq sybr green supermix. with stage as reference, fold change was calculated by normalizing ct values in stages and against odc gene using the 2−ΔΔct method. bioinformatic and manual curation of ngm-vv subset ngm-vv subset is collection of new gene models. as a first step, we filtered these models for orf length less than or equal to amino acids. this resulted in a set of gene models, which were then screened using following criteria: 1) absence of a downstream gene with a x or u annotation for h3k4me; 2) rpkm of all exons should be greater than or equal and 3) evidence of splicing in our data. this resulted in a set of gene models. these models were then manually curated using blastn and blastp to filter against homology to known protein coding sequences. conservation analysis for conservation analysis all gene models were mapped to jgi v using blat. the average phastcons score per gene model was calculated using the conservation track of the ucsc genome browser. data availability the data have been deposited in ncbi′s gene expression omnibus and are accessible through geo series accession number gse xtev gene models are available at: http://veenstra.ncmls.nl/genomedata.asp. statement on animal use animal care and use was in accordance with national and european guidelines and standard operating procedures approved by the institutional animal care and use committee. abbreviations polya+: polyadenylated; rt-qpcr: real time rt-pcr; cdna: complementary dna; rna-seq: ribonucleic acid sequencing; mrna: messenger rna; mbt: mid-blastula transition; eden: embryonic deadenylation element; eden-bp: embryonic deadenylation element - binding protein; utr: untranslated] region; xtev: xenopus tropicalis experimentally validated; lncrnas: long non-coding rnas; orf: open reading frame; pa: polya+; ribozero, rz: ribosomal rna depleted total-rna; rpkm, see materials and methods: reads per kilobase of exon model per million mapped reads; go: gene ontology; bp: biological processes; est: expressed sequence tags; gm: gene models; ngm: new gene models; ngm-vv: new gene models with validation support; ngm-vvo: manually curated new gene models with orf less than amino acids; llr: log likelihood ratio; jgi: joint genome institute; rrna: ribosomal rna; bwa: burrows-wheeler aligner; pcr: polymerase chain reaction; h3k4me3: histone h lysine tri-methylation; rz: ribosomal rna depleted total rna; rnapii: rna polymerase ii; chip-seq: chromatin immuno-precipitation sequencing; tss: transcription start-site. competing interests the authors declare that they have no competing interests. authors’ contributions ssp, ugj and gjcv conceived and designed the study, ugj performed the rna-seq experiments, ssp did the analysis and performed the validation experiments. sjvh performed transcriptome assembly, designed and assembled the xenopus tropicalis experimentally validated pipeline and performed the conservation analysis. gjcv performed orf length and codon bias analysis. ssp and gjcv wrote the manuscript. all authors read and approved the final manuscript. supplementary material additional file 1 supplementary tables and figures. click here for file additional file 2 multi-page list of genes expressed in xenopus tropicalis. click here for file additional file 3 multi-page xtev annotation pipeline gene list with genomic co-ordinates. click here for file additional file 4 single-page list of codon usage frequencies used for calculating codon bias score. click here for file additional file 5 single-page list of primers used for real-time qpcr validation of highly conserved lncrnas from the ngm subset. click here for file acknowledgements we thank dr.mike gilchrist for providing the est data. this work was supported by grants of the national institutes of health and nwo, the netherlands organization of scientific research.
to analyze the dynamics and diversity of coding and non-coding transcripts during development, both polyadenylated mrna and ribosomal rna-depleted total rna were harvested across six developmental stages and subjected to high throughput sequencing. decay of a subset of maternal mrnas and new transcription define the maternal-to-zygotic transition, but the full complement of polyadenylated and deadenylated coding and non-coding transcripts has not yet been assessed in xenopus embryos. many maternal genes show peak expression in the oocyte and include genes which are known to be the key regulators of events like oocyte maturation and fertilization. of all the transcripts that increase in abundance between early blastula and larval stages, about 30% of the embryonic genes are induced by fourfold or more by the late blastula stage and another 35% by late gastrulation. using a gene model validation and discovery pipeline, we identified novel transcripts and putative long non-coding rnas.
 dynamics of polyadenylation vs. deadenylation determine the fate of several developmentally regulated genes. decay of a subset of maternal mrnas and new transcription define the maternal-to-zygotic transition, but the full complement of polyadenylated and deadenylated coding and non-coding transcripts has not yet been assessed in xenopus embryos.  to analyze the dynamics and diversity of coding and non-coding transcripts during development, both polyadenylated mrna and ribosomal rna-depleted total rna were harvested across six developmental stages and subjected to high throughput sequencing. the maternally loaded transcriptome is highly diverse and consists of both polyadenylated and deadenylated transcripts. many maternal genes show peak expression in the oocyte and include genes which are known to be the key regulators of events like oocyte maturation and fertilization. of all the transcripts that increase in abundance between early blastula and larval stages, about 30% of the embryonic genes are induced by fourfold or more by the late blastula stage and another 35% by late gastrulation. using a gene model validation and discovery pipeline, we identified novel transcripts and putative long non-coding rnas. these lncrna transcripts were stringently selected as spliced transcripts generated from independent promoters, with limited coding potential and a codon bias characteristic of noncoding sequences. many lncrnas are conserved and expressed in a developmental stage-specific fashion.  these data reveal dynamics of transcriptome polyadenylation and abundance and provides a high-confidence catalogue of novel and long non-coding rnas. xenopus tropicalisrna-seqmaternal and embryonic transcriptomepolyadenylationdeadenylationmbtcodon biaslong-noncoding rnas 