A dual platform approach to transcript discovery for the planarian Schmidtea mediterranea to establish RNAseq for stem cell and regeneration biology

Abstract

The use of planarians as a model system is expanding and the mechanisms that control planarian regeneration are being elucidated. The planarian Schmidtea mediterranea in particular has become a species of choice. Currently the planarian research community has access to this whole genome sequencing project and over 70,000 expressed sequence tags. However, the establishment of massively parallel sequencing technologies has provided the opportunity to define genetic content, and in particular transcriptomes, in unprecedented detail. Here we apply this approach to the planarian model system. We have sequenced, mapped and assembled 581,365 long and 507,719,814 short reads from RNA of intact and mixed stages of the first 7 days of planarian regeneration. We used an iterative mapping approach to identify and define de novo splice sites with short reads and increase confidence in our transcript predictions. We more than double the number of transcripts currently defined by publicly available ESTs, resulting in a collection of 25,053 transcripts described by combining platforms. We also demonstrate the utility of this collection for an RNAseq approach to identify potential transcripts that are enriched in neoblast stem cells and their progeny by comparing transcriptome wide expression levels between irradiated and intact planarians. Our experiments have defined an extensive planarian transcriptome that can be used as a template for RNAseq and can also help to annotate the S. mediterranea genome. We anticipate that suites of other 'omic approaches will also be facilitated by building on this comprehensive data set including RNAseq across many planarian regenerative stages, scenarios, tissues and phenotypes generated by RNAi.

Figure 1. Overview of pipeline obtaining planarian transcripts.

A schematic flow chart showing the work-flow for defining the planarian transcript set using a dual platform approach. 20 size-matched worms for each time point were cut into head, middle and tail pieces. Regenerating head, tail and middle fragments for RNA preparation were collected and frozen at 6, 12, 24, 36, 48, 72, 96, 120, 144 and 168 hours of regeneration. Total RNA was prepared from these time points, pooled and libraries made for both platforms.

Figure 2. 454 read length distribution.

A plot showing the distribution of the raw planarian 454 reads generated from one pico-titre plate run with Titanium chemistry.

Figure 3. Increasing isogroup/isotig numbers and contig length with increasing read number.

A. The number of Isotigs and Isogroups continue to increase with increasing number of reads added to the assembly, however Isogroup (gene) discovery appears to be a approaching a maximum whereas Isotig discovery (transcript isoforms) continues to increase. This suggests further sequencing would be beneficial to define more alternately spliced transcripts. B. The lengths of contigs also increase as more reads are added to the assembly. While many contigs are not any getting longer (average contig length) further sequencing is increasing the N50 perhaps by joining some transcripts together to form single contigs. Both analyses suggest that further sequencing of the mixed stage library may generate more transcripts and continue to increase the average length of contigs.

Figure 4. Contig length distributions.

Contig length distribution for the following datasets (A) Sm454EST (B) SmABI before iterative mapping (C) SmABI after iterative mapping (D) Sm454ESTABI merged dataset. The subplots for each figure focus on contigs with length up to 2,500 bp. The short read approach using SOLiD clearly produces many short putative transcripts, many of these may be transcribed regions that do not encode proteins or which splice reads were not identified. Iterative mapping resulted in fewer short transcribed regions after filtering and longer contigs (B vs C) Combining the datasets led to an increase in contig lengths compared to either platform alone.

Figure 5. Defining transcript variants in more detail.

A schematic flow chart showing the iterative mapping approach. New split reads found in the second round of mapping are in red.

Figure 6. Open reading frame length distribution.

Length distribution of the longest putative protein that could be coded for by each transcript for the Sm454ESTABI dataset (black). 11,599 open reading frames had both a putative start and stop codon (red), 9,076 missed the start (green) codon and 1,757 missed the stop codon (orange). The remaining 2620 transcripts either had neither a start or stop.

Figure 7. Cross-species BLASTX results.

The BLASTX cross-species results for the Sm454ESTABI dataset. All species with a significant hit (probability < 1e-5 red bar, 1e-10 orange bar) are shown. Also shown are the number of times the top hit was to a particular species (black bar and number in brackets).

Figure 8. Gene ontology hits.

The number of gene ontology (GO) hits for the second tier categories: (A) Biological process (B) Molecular function (C) Cellular component. More detailed GO term assignment is available as an additional file.

Figure 9. Distribution of ORF and transcript lengths with and without homology based annotation.

Frequency plots of open reading frame length (A) and total transcript length separated by whether a transcript had homology (black) to another species by BLAST, PFAM and/or GO assignment or did not have any homology (orange) while shorter coding regions and transcript lengths are less likely to have homology this analysis suggests that many short transcripts and of the short ORFs (in the 0–50 aa range) are likely to be non-coding transcripts of the untranslated regions of non-coding transcripts.

Figure 10. RNA seq expression values of Whole vs Irradiated expression.

Plots show Log 10 transformed RPKM values for transcripts detected above an RPKM of >1 (Log 10 RPKM>0) in at least one sample. We compared the two replicates of RNAseq of irradiated worms (A) and intact (mock irradiated) to each other demonstrating a high level of congruency with differences in expression values only apparent at low expression levels. We then compared the replicates of irradiated worms against the replicates of intact worms (C) and are able to detect many differentially expressed transcripts, many of which decrease after irradiation. These represent a group of genes potentially expressed in planarian stem cell and their recent progeny.