A quantitative reference transcriptome for Nematostella vectensis early embryonic development: a pipeline for de novo assembly in emerging model systems
© Tulin et al.; licensee BioMed Central Ltd. 2013
Received: 4 January 2013
Accepted: 14 March 2013
Published: 3 June 2013
The de novo assembly of transcriptomes from short shotgun sequences raises challenges due to random and non-random sequencing biases and inherent transcript complexity. We sought to define a pipeline for de novo transcriptome assembly to aid researchers working with emerging model systems where well annotated genome assemblies are not available as a reference. To detail this experimental and computational method, we used early embryos of the sea anemone, Nematostella vectensis, an emerging model system for studies of animal body plan evolution. We performed RNA-seq on embryos up to 24 h of development using Illumina HiSeq technology and evaluated independent de novo assembly methods. The resulting reads were assembled using either the Trinity assembler on all quality controlled reads or both the Velvet and Oases assemblers on reads passing a stringent digital normalization filter. A control set of mRNA standards from the National Institute of Standards and Technology (NIST) was included in our experimental pipeline to invest our transcriptome with quantitative information on absolute transcript levels and to provide additional quality control.
We generated >200 million paired-end reads from directional cDNA libraries representing well over 20 Gb of sequence. The Trinity assembler pipeline, including preliminary quality control steps, resulted in more than 86% of reads aligning with the reference transcriptome thus generated. Nevertheless, digital normalization combined with assembly by Velvet and Oases required far less computing power and decreased processing time while still mapping 82% of reads. We have made the raw sequencing reads and assembled transcriptome publically available.
Nematostella vectensis was chosen for its strategic position in the tree of life for studies into the origins of the animal body plan, however, the challenge of reference-free transcriptome assembly is relevant to all systems for which well annotated gene models and independently verified genome assembly may not be available. To navigate this new territory, we have constructed a pipeline for library preparation and computational analysis for de novo transcriptome assembly. The gene models defined by this reference transcriptome define the set of genes transcribed in early Nematostella development and will provide a valuable dataset for further gene regulatory network investigations.
KeywordsTranscriptome Gene regulatory network Nematostella embryonic development Body plan evolution Next-generation sequencing Illumina HiSeq Trinity Oases RNA-seq
External RNA Controls Consortium
fragments per kilobase of exon per million fragments mapped
filtered sea water
gene regulatory network
- N50 GO:
next generation sequencing
National Institute of Standards and Technology
non-redundant BLAST database
(quantitative) polymerase chain reaction.
Nematostella vectensis, the starlet sea anemone, offers many advantages as a model system for the evolution of animal developmental programs. As an anthozoan cnidarian, it is strategically positioned as an outgroup to Bilateria [1–3] and is well situated to reveal the early steps in the evolution of the bilaterian body plan. Two of these evolutionary steps are likely to include the formation of a secondary body axis and a mesodermal germ layer which are both essential, defining characteristics of a bilaterian animal. Embryonic dorsal-ventral patterning and mesodermal development have been studied in many bilaterian models yet the origins of these significant body plan innovations are not well understood. Initial studies of gene expression in Nematostella and non-anthozoan cnidarians have revealed that genes important to bilaterian mesoderm specification are expressed in the endoderm of the sea anemone, and suggests that the bilaterian mesoderm may have originated from the endoderm of diploblastic ancestors [4–6]. Genes encoding factors involved in dorsal-ventral axis specification in Bilaterians are likewise asymmetrically expressed in Nematostella, indicating the possibility that a secondary axis was present in the Cnidarian-Bilaterian ancestor [7, 8]. Defining the mechanisms controlling Nematostella development will help address these questions about the early evolutionary steps that led to bilaterian body plans with three germ layers and bilateral symmetry.
Gene regulatory networks (GRN) provide predictive models of gene regulation, as in the several examples that now exist for normal animal development (for example, Drosophila, sea urchin [10, 11], ascidians , chick , and zebrafish ). To gain a comprehensive view of the control system, it is necessary to identify all genes whose products make up the regulatory network. This applies to our current research efforts but is also generally applicable to studies of virtually any regulatory system. Advanced sequencing platforms now allow us to do this through RNA-seq techniques. Yet, deep RNA-seq brings challenges in analysis reflecting the scale and complexity of transcriptomes, the primary problem being adequate assembly of RNA-seq reads in order to define a reference set of gene models [15–17]. Transcriptome assembly can be achieved using a reference-based strategy, a de novo strategy or a combination of the two. The main drawback to using a genome reference for assembly is that it relies on the quality of the reference genome being used . This is a particular problem for emerging model systems with recently completed genomes because misassemblies, poor annotation and large gaps in coverage plague the genome assemblies of all but a few of the major model systems . There is also a challenge in assigning reads that align equally well to multiple places in the genome. The aligner must decide to either exclude these reads which can result in gaps or to choose which alignments to retain which could lead to wrong assignments or predictions of a transcript in a region that has no transcription.
A comprehensive GRN for early embryonic development in Nematostella will enable researchers to investigate the extent to which the bilaterian regulatory toolkit is present in this representative cnidarian, down to the level of precise signaling systems and transcription factor cis-regulatory interactions. By harnessing the power of high-throughput sequencing and perturbation techniques, we aim to build the sea anemone GRN in an unbiased and efficient manner that will serve as a GRN construction pipeline for other model systems to follow.
The current Nematostella genome assemblies [20, 21] fall into the category of young genome models that are still incomplete and contain gaps thus making the reference-based method alone insufficient for our needs. Taking these and all of the above complications into account and considering our goal to define an experimental and computational pipeline for emerging model systems, we elected to use the de novo assembly approach. This approach will be especially useful for evo-devo researchers aiming to harness the power of next-generation sequencing to bring their research into the genomics era; a trend already underway, for example Parhyale, Oncopeltus, sponge , and sea urchin .
The scale of reads, random and non-random sequencing errors, and inherent transcript complexity due to alternate transcription start sites or splice junctions all pose challenges for de novo transcriptome assembly. Indeed, the scale of the problem is only set to increase with the expanding capacity for transcriptome sequencing from advances in next-generation sequencing (NGS) platforms. In the last few years several assembly algorithms have been released to meet these challenges: Trans-ABySS , SOAPdenovo , Velvet/Oases [27, 28], and Trinity . The millions of short reads produced from NGS platforms result in millions of overlapping sequences. Short-read de novo assemblers exploit these overlaps to reconstruct the original transcripts by using the de Bruijn graph data structure, which encodes overlapping k-mers as adjacent vertices. Assembly algorithms then compute paths through the de Bruijn graph that correspond to valid assemblies of the sequence reads.
The Trinity assembler breaks the sequence data into many de Bruijn graphs in order to capture transcript complexity resulting from alternative splicing, gene duplications or allelic variation . Trinity consists of three modules called Inchworm, Chrysalis, and Butterfly. Inchworm assembles the RNA-seq reads into transcripts and reports only the unique portions of alternate transcripts. Chrysalis clusters the Inchworm contigs so that each cluster represents all known transcripts variants for each gene or related genes and then constructs De Bruijn graphs for each cluster. All reads are segregated to one of these separate graphs. Butterfly then processes these separate graphs in parallel by tracing a path through each one and reports full length transcripts separately for alternate splice forms and paralogs. The Oases assembler uploads a preliminary assembly created by Velvet, which was originally designed for genome assembly. Oases corrects this assembly using a range of k-mers to create separate assemblies, which are then combined into one. The longer k-mers perform better on high expression transcripts and the shorter k-mers have an advantage on low expression transcripts . While the multiple k-mer approach has been found to result in an increase of longer transcripts, it can also lead to an accumulation of incorrect assemblies or artificially fused transcripts .
In this study we designed a next-generation sequencing and analysis pipeline to produce a minimally biased and quantitative reference transcriptome. The resulting transcriptome represents the first 24 h of Nematostella development and will be the basis for further gene regulatory network studies. The experimental and computational pipeline will be used by us and others to produce transcriptomes for other model systems, particularly those evo-devo models that do not yet have an annotated genome but would benefit from an in depth molecular analysis.
Nematostella vectensis adults following normal culture at 18°C were spawned with a 9-h cycle of light at 25°C in an incubator. Male and female spawning adults were in separate bowls and egg sacs were removed to a fresh bowl and fertilized with sperm from male bowls for 10 minutes. The egg sacs were then dejellied with a 4% cysteine solution (pH 7.4) in 50% filtered sea water (FSW) for 8 minutes and rinsed five times with 50% FSW. All embryo processing was performed in an 18°C room and the embryos were cultured from the time of fertilization for 0, 6, 12, 18 or 24 h (five timepoints). An additional 24-h sample was prepared in the same way from a separate spawning event. Cultured embryos were transferred to an eppendorf tube, allowed to settle, gently spun to a pellet and the supernatant removed, approximately 600 embryos per sample. The embryo pellet was immediately immersed in 100 μl of lysis buffer from the Invitrogen mRNA DIRECT kit (Invitrogen, Life Technologies, Grand Island, NY, USA) and homogenized with a Kontes Pellet Pestle (distributed by Thermo Fisher Scientific, Pittsburgh, PA, USA) attached to a 12 V/700 rpm drill. Another 100 μl of lysis buffer was used to rinse the Kontes Pellet Pestle tip and collected in the same tube. Samples were then stored at −80°C until all timepoints had been collected.
To thawed lysates, a third aliquot of 100 μl of lysis buffer was added and then the normal protocol for the Invitrogen mRNA DIRECT kit was followed using 50 μl Dynabeads per sample and low adhesion microcentrifuge tubes, following the manufacturer’s recommendations. The mRNA yields were between 108 ng and 344 ng total per sample. The mRNA was used as starting material for the ScriptSeq V.1 kit from Epicentre (Epicentre Biotechnologies, Madison, WI, USA). A total of 9.0 μl of mRNA corresponding to between 74 ng to 233 ng per sample was combined with 1.0 μl of a 1:10 dilution of External RNA Controls Consortium (ERCC) spike-in control RNA for the first reaction (available from Invitrogen/Life Technologies). The protocol was followed exactly, using 12 cycles total of PCR in the amplification step with Phusion High Fidelity polymerase (available from Therm Scientific) and barcoded Illumina-compatible primers 1 to 6 from Epicentre.
The libraries were size selected with a 2% Pippin prep gel (from Sage Science, Beverly, MA, USA) for 450 bp and checked on a Agilent 2100 Bioanalyzer with a high sensitivity DNA chip (from Agilent Technologies, Santa Clara, CA, USA) and then by qPCR. The samples were combined and run on a single lane of the Illumina High Seq 1000 with version III chemistry with 200 cycles of paired end sequencing plus indexing reads. All raw read files are available on the Woods Hole Data Archive at http://hdl.handle.net/1912/5613, DOI [DOI:http://10.1575/1912/5613].
Quality control was implemented using a combination of Bowtie2 (version 2.0.0-beta6), Basic Local Alignment Search Tool (BLAST), btrim (build date 9 September 2011), and the FASTX-Toolkit. First, we computed overrepresented k-mers in the raw sequence data and ran BLAST against the set of non-redundant (nr) sequences from the National Center for Biotechnology Information (NCBI). The top several BLAST hits were analyzed for homology with ribosomal or mitochondrial RNA. Using the BLAST data, we constructed a catalog of rRNA and mtRNA sequences which could serve as a reference set to filter the overly abundant non-protein coding RNA prior to de novo transcriptome assembly. A Bowtie2 index was built from the rRNA and mtRNA sequences of Montastraea franksi, Savalia savaglia, Actiniaria, Nematostella vectensis, and Clathrina. Sequence reads successfully aligning by Bowtie2 to this set were removed. The overrepresented k-mers also contained adapter sequences that remained in the sequence. We retrieved the exact Illumina adapter sequences and used the software tool btrim to clip adapters.
Sequence reads demonstrating low complexity (containing only one or two unique bases) are likely due to technical artifacts and were removed. Next, the GC content distribution was computed for the set of all reads. GC content biases in the first 13 bases of Illumina RNA-seq data are known to exist due to random hexamer priming . This bias may cause an imbalance in read coverage and persist through the assembly process, which can affect the quality of assembly and quantification levels. Because we had extremely high sequence coverage, we removed the bias by simply trimming the start of the reads. Using the FASTX-Toolkit, we removed the initial 13 bases from the reads at each timepoint. Finally, btrim was also used to adaptively trim low quality bases from the end of the read. Adaptive trimming is performed by sliding a window of 5 bp from the end of the read to the start, removing bases and shifting the sliding window by 1 base if the average quality score is less than 30 until the average quality score is at least 30.
Digital normalization, Velvet and Oases
Digital normalization is a method to reduce the total number of reads to be assembled, thereby also reducing the computing power and time required for assembly. It preferentially removes high abundance reads but retains read complexity in order to remove errors but preserve low abundance transcripts prior to assembly. All links to digital normalization software are available electronically through http://ged.msu.edu/papers/2012-diginorm/. Raw paired-end read files were first interleaved into pairs using a python script, available at http://github.com/ged-lab/khmer/tree/2012-paper-diginorm/sandbox. Then, three rounds of digital normalization were applied to remove overabundant and erroneous reads. These depend on the khmer software package, available at http://github.com/ged-lab/khmer/. The khmer software also relies on the screed package for loading sequences, available at http://github.com/ged-lab/screed/ (khmer and screed are ©2010 Michigan State University, and are free software available for distribution, modification, redistribution under the BSD license). The digital normalized files were assembled with Velvet (version 1.2.03) and Oases (version 0.2.06). The details of the execution commands are available in Additional file 1. The most current recommendations for use of digital normalization for de novo transcriptome assembly recommend using only one round of normalization instead of three. Fewer low abundance transcripts may be lost by foregoing further rounds of digital normalization at the expense of increased computing time and power to assemble the greater number of remaining reads.
Trinity assembly and quantification
Spike-in standard curve R values for the ordinary least squares regression
To compute overexpressed GO terms in our transcriptome, we used BLASTx 2.2.26+, BLOSUM62 similarity matrix, Blast2GO database version August 2011, and pipeline B2G4Pipe version 2.3.5. The definition of each GO term is determined by the GO Consortium: http://www.geneontology.org/, and can be found using the EMBL-European Bioinformatics Institute QuickGO: http://www.ebi.ac.uk/QuickGO/, or the Gene Ontology Normal Usage Tracking System, GONUTS: http://gowiki.tamu.edu/wiki/index.php/Main_Page. Definitions for all GO terms presented in this paper can be found in Additional file 4.
Library preparation and quality control of reads
The polyA-RNA enriched sample was then processed with the Script Seq kit, version 1, from Epicentre. The main advantages of this kit are the resulting directional reads, short preparation time (4 h), and low input requirements (as low as 50 ng mRNA). The adapter-ligated libraries were then size selected for uniformity at 450 bp using a Pippin Prep gel electrophoresis apparatus from Sage Science and combined in one lane on an Illumina HiSeq 1000 to produce 2 × 100 bp paired-end reads.
Quality control read attrition
Raw paired end reads
QC step 1
QC step 2
1.74E + 07
1.29E + 07
1.24E + 07
1.98E + 07
1.46E + 07
1.40E + 07
1.44E + 07
9.17E + 06
8.76E + 06
2.88E + 07
1.96E + 07
1.89E + 07
3.61E + 07
2.57E + 07
2.46E + 07
2.75E + 06
1.87E + 06
1.77E + 06
1.19E + 08
8.39E + 07
8.05E + 07
De novo assembly with Trinity assembler
No. of distinct transcripts passing filters
Total number of assembled transcripts
≥1 BLASTx hit (Nematostella vectensis), e <5e-5
≥1 BLASTx hit (non-redundant), e <5e-5
≥1 BLASTx hit (non-redundant) and 80% similarity
≥1 BLASTx hit (Nematostella vectensis) and 80% similarity
Molecules per embryo (MPE) >100
MPE >100, ≥1 BLASTx hit (non-redundant)
Transcript families from transcripts with MPE >100, ≥1 BLASTx hit (non-redundant)
Trinity assembly compared to digital normalization/Velvet/Oases
Digital normalization, Oases-Velvet
Aligned concordantly 0 times
Aligned concordantly 1 time
Of reads aligning concordantly or discordantly 0 times:
Aligned discordantly 1 time
Aligned 0 times
Aligned exactly 1 time
Aligned >1 times
Overall alignment rate
Determining adequate coverage of transcriptomes is more challenging than determining coverage of genomes because every transcript species (including splice variants or those using alternate transcription start sites) is present at a different level across a large range. Several studies have used an independent assembly of randomly selected subsets of their reads to compare the rate of new transcript discovery, determine the lower abundance limit of detection and compare the average length of isotigs. While analyzing the sea urchin embryonic transcriptome, Tu et al. compared a 20 M read subset with a 2 M read subset and a 0.2 M read subset, and found that 20 M reads were sufficient to reliably detect levels of transcripts at 400 molecules/embryo, which they estimate as the lower limit for proteins of developmental significance (such as transcription factors, which may be functionally relevant at levels as low as 10 copies of transcript per cell) . In their analysis of the milkweed transcriptome, Ewen-Campen et al. created eight subsets of reads, assembled them separately and used BLASTx to compare gene discovery rates . They found that the rate of new transcript discovery plateaued at 1.5 M reads, although the N50 isotig length continued to increase when using 2 M to 17 M read subsets. After 2 rounds of quality control filtering of our reads, we were left with 80,537,812 paired and 5,362,854 unpaired reads, a depth which has been shown to produce good sensitivity in these other systems for identifying all protein-coding transcripts expressed in the early embryo.
To restate, these previous studies indicate that with the volume of reads coming off the latest Illumina platforms (250 million to 400 million reads/lane), and only multiplexing 6 samples in a lane, we should be beyond the necessary coverage to represent all relevant regulatory transcripts. The best indication of the quality of our assembly is that we have been able to use it as a reference to map reads from more recent experiments in our lab at a median 93% rate (with 90% of samples mapping 90% of their reads to the Trinity-assembled reference).
Digital normalization followed by Oases assembly
To evaluate more closely the quality of our de novo transcriptome assembly, we compared Trinity with an alternate strategy that combines a digital normalization step  with the assemblers Velvet  and Oases . Digital normalization is a computational normalization method that preferentially removes high abundance reads but retains read complexity in order to remove errors and preserve low abundance transcripts prior to assembly. The quality controlled reads were assembled using Velvet and Oases and then mapped back to the resulting assembly (commands in Additional file 1). The main advantage of this method is it greatly decreases the computing power and time required to process millions of reads. We also tested the Amazon Elastic Cloud Computing Service (Amazon EC2), http://aws.amazon.com/ec2/, to perform the digital normalization and Velvet-Oases assembly. This method proves to be a great alternative to using a home institution’s core computers if the institution does not have sufficient computing power for running an assembler or if the computers are expensive to rent, slow or unreliable. Whereas the computation of Trinity required over 50 h and 100 GB of RAM in addition to the quality control steps, the pipeline using digital normalization, Velvet and Oases can all be run in a single day using an XL computer rented from Amazon. We found that this alternative assembly approach gave competitive results when mapping our Nematostella vectensis RNA-seq reads. The overall mapping success rate was 82.17% for digital normalization-Oases as compared to 85.90% for Trinity (Table 4).
Quantification of Trinity-assembled transcriptome using known RNA ‘spike-ins’
Blast2GO analysis reveals genes involved in gene regulation
Transcript family analysis
The transcripts inferred by the Trinity assembler were tested for sequence similarity with known genes. Specifically, the sequence of each transcript was translated and BLASTx was run using the NCBI non-redundant RefSeq database of protein sequences (nr). We combined sequences sharing at least 1 BLAST result with an e-value <5e-5 and a percent identity >80% into a transcript family. The e-value threshold filtered out low confidence BLAST results, while the percent identity filter requires a large percentage of the transcript to match. We systematically tested multiple identity thresholds and found that 80% sequence match yielded a realistic transcript family distribution. Given the natural tradeoff between identification of transcript variants versus paralogous genes, overly stringent requirements of similarity result in an underrepresentation of true homologous relationships; conversely, when similarity thresholds are set too low, distinct transcripts are erroneously grouped together. Because the assembled Nematostella transcripts are also included in nr, the sequence similarity threshold is required to be set high. The transcript families therefore represent a mixture of fully assembled gene transcripts, pieces of transcripts, paralogs, or multiple splice forms, entirely compatible with our overarching goals. For BLAST analysis with sequence databases of very different sizes (as in our case, nr vs Nematostella) e-values are not an appropriate measure of similarity; in this case a similarity threshold is more informative. Thus, at a similarity threshold of 80% we were able to annotate 48,235 transcripts from the nr database, compared to only 47,193 transcripts when using the Nematostella genome alone (Table 3). To do transcript family analysis we used the BLAST hits from the nr database. We inferred 13,293 total transcript families from the 61,835 assembled transcripts with at least 1 BLAST hit. When restricting the transcript family analysis to transcripts expressed over 100 molecules per embryo (MPE) and at least 1 BLAST hit (MPE >100, ≥1 BLAST hit) across the 5 timepoints, we observed a total of 4,055 transcript families from the 6,169 transcripts passing the filter (Table 3). These computations likely represents an underestimate of the true number of genes expressed due to an inability to assemble very lowly expressed transcripts and, in a few cases, grouping paralogous genes together.
An estimation of the percentage of the genome transcribed during these time periods was computed by taking the length of the longest transcript in each transcript family and dividing that by the length of the genome taken from the estimate in the Nematostella genome paper (450 Mbp) . An important caveat is that this is likely an underestimate because some transcripts are not full length. The percent of the genome transcribed above 100 molecules per embryo according to this calculation is 0.368%. The average transcript length for all assembled transcripts is 622.53 bp. When only taking transcripts expressed after spike-in control correction above 0 molecules per embryo, the average transcript length is 456.18 bp. By timepoint, the average of assembled and expressed transcripts is: (0 h) 479.95 bp, (6 h) 487.82 bp, (12 h) 490.54 bp, (18 h) 494.03 bp, and (24 h) 441.29 bp.
The goals of the project discussed in this paper were (1) to identify all the protein-coding genes expressed during the first 24 h of Nematostella vectensis development and in so doing (2) to detail a modern, cost-effective, efficient and quantitative series of experimental and computational methods that together make up a transcriptome pipeline for non-model organisms. A key component of our pipeline is the inclusion of NIST RNA spike-in standards for quantification. This allows us to get around the problem of normalizing data to estimate gene expression levels , and provides an absolute measure of transcript abundance per embryo.
Many evolutionary developmental biology ‘evo-devo’ research projects have revealed candidate genes in non-model organisms leading to intriguing hypotheses regarding the conservation, or conversely, invention of pathways controlling development [5, 44, 45]. However, to answer the questions these hypotheses have generated, it is not only the gene homology, presence, absence or spatial localization that needs to be known. To say that a developmental program or subcircuit has been conserved or evolved in a specific way, the cis-regulatory network connections between all the regulatory genes involved must be at a minimum known and validated. Candidate genes identified from BLAST analysis will typically only make up a small fraction of the regulatory genes in any pathway. With the advent of next-generation sequencing platforms, identifying all the protein coding genes expressed at any given time during embryonic development is now within the reach of any model system where embryos can be acquired. The lack of a sequenced, annotated genome is no longer a major setback to GRN analysis.
The use of polyadenylated spike-in RNAs provides quantitative information on the absolute abundance of transcripts per embryo. It is important to note the difference between this method of standardization and normalization approaches. The ERCC spike-ins allow us to build a standard curve, in our case a 92-point standard curve. As the quantities of the spike-ins are known, this allows us to infer from the standard curve absolute quantities. Note that since spike-ins are added at the beginning of the library preparation procedures, any variation in preparation efficiencies (that is, technical noise) is in theory accounted for by the spike-ins. Thus, even without absolute quantitation, the use of spike-ins allows direct comparison between samples without the distorting effects of normalization to minimize the effects of technical variation. Further, quantitation by spike-ins also allows us to know the limits of our ability to detect and quantify lowly expressed transcripts. Since low expressed transcripts account for many of the problems in bioinfomatics analysis, our 100 molecules per embryo cut-off allows us to focus our analysis on those transcripts expressed at biologically relevant levels which are also within the linear range of our standard curve. Increasing the sequencing depth and being less conservative with mapping stringency could improve our ability to quantify these lowly expressed transcripts.
The embryo of the sea anemone Nematostella vectensis provides an important evo-devo model for understanding early animal development, particularly in relation to the question of how initial patterns of differential gene expression emerge along orthogonal body axes. Given Nematostella’s position among cnidarians and the molecular evidence thus far, it is possible that a bilaterally-symmetric pattern formation network stretches back to before the Cambrian to a time preceding the Cnidaria-Bilateria bifurcation. However, to make this argument we need a mechanistic understanding of early development in both cnidarian and canonical bilaterian models systems. Moreover, in light of the compare-and-contrast nature of these studies, we need to move away from a candidate gene approach as such methods clearly bias towards the ‘discovery’ of similarities as opposed to differences between regulatory networks. With the advent of genomics, we can now attempt exhaustive de novo approaches to define regulatory networks, though challenges in handling RNA-seq data sets still exist. In this report, we have undertaken preliminary steps in defining the Nematostella gene regulatory network for early pattern formation by building a comprehensive model of gene expression through 24 h of development. This quantitative reference transcriptome will help us identify, in a minimally biased manner, the most relevant genes to the pattern formation control system. The regulatory network for pattern formation in Nematostella will in turn provide a powerful basis for comparison with early development networks from canonical bilaterians.
In summary, we have presented our quantitative reference transcriptome for Nematostella vectensis early embryogenesis, which is available on the Woods Hole Data Archive at http://hdl.handle.net/1912/5613 [DOI:http://10.1575/1912/5613]. Additionally, our de novo transcriptome pipeline, based on the Trinity assembler, has been designed to meet the needs of the evo-devo community.
We give special thanks to Liliana Florea for her valuable contributions to the design of the Trinity de novo assembly pipeline used in this paper. We thank Kasia Hammar for animal care and technical assistance; Likit Preeyan and C Titus Brown for help with digital normalization and for credit on the Amazon EC2 system; Antje Fischer, Rebecca Helm, Freya Goetz, and Casey Dunn for discussions on experimental protocols; and C Titus Brown also for comments on the manuscript.
- Dunn CW, Hejnol A, Matus DQ, Pang K, Browne WE, Smith SA, Seaver E, Rouse GW, Obst M, Edgecombe GD, Sørensen MV, Haddock SHD, Schmidt-Rhaesa A, Okusu A, Kristensen RM, Wheeler WC, Martindale MQ, Giribet G: Broad phylogenomic sampling improves resolution of the animal tree of life. Nature 2008, 452:745–749.PubMedView Article
- Hejnol A, Obst M, Stamatakis A, Ott M, Rouse GW, Edgecombe GD, Martinez P, Baguna J, Bailly X, Jondelius U, Wiens M, Muller WEG, Seaver E, Wheeler WC, Martindale MQ, Giribet G, Dunn CW: Assessing the root of bilaterian animals with scalable phylogenomic methods. Proc R Soc B 2009, 276:4261–4270.PubMedView Article
- Mallatt J, Craig CW, Yoder MJ: Nearly complete rRNA genes assembled from across the metazoan animals: effects of more taxa, a structure-based alignment, and paired-sites evolutionary models on phylogeny reconstruction. Mol Phylogenet Evol 2010, 55:1–17.PubMedView Article
- Martindale MQ, Pang K, Finnerty JR: Investigating the origins of triploblasty: “mesodermal” gene expression in a diploblastic animal, the sea anemone Nematostella vectensis (phylum, Cnidaria; class, Anthozoa). Development 2004, 131:2463–2474.PubMedView Article
- Wikramanayake AH, Hong M, Lee PN, Pang K, Byrum CA, Bince JM, Xu R, Martindale MQ: An ancient role for nuclear beta-catenin in the evolution of axial polarity and germ layer segregation. Nature 2003, 426:446–450.PubMedView Article
- Fritzenwanker JH, Saina M, Technau U: Analysis of forkhead and snail expression reveals epithelial-mesenchymal transitions during embryonic and larval development of Nematostella vectensis . Dev Biol 2004, 275:389–402.PubMedView Article
- Matus DQ, Thomsen GH, Martindale MQ: Dorso/ventral genes are asymmetrically expressed and involved in germ-layer demarcation during cnidarian gastrulation. Curr Biol 2006, 16:499–505.PubMedView Article
- Rentzsch F, Anton R, Saina M, Hammerschmidt M, Holstein TW, Technau U: Asymmetric expression of the BMP antagonists chordin and gremlin in the sea anemone Nematostella vectensis : implications for the evolution of axial patterning. Dev Biol 2006, 296:375–387.PubMedView Article
- Stathopoulos A, Levine M: Genomic regulatory networks and animal development. Dev Cell 2005, 9:449–462.PubMedView Article
- Davidson EH, Levine MS: Properties of developmental gene regulatory networks. Proc Natl Acad Sci USA 2008, 105:20063–20066.PubMedView Article
- Davidson EH: The Regulatory Genome. London, UK: Academic Press; 2006.
- Imai KS, Levine M, Satoh N, Satou Y: Regulatory blueprint for a chordate embryo. Science 2006, 312:1183–1187.PubMedView Article
- Betancur P, Bronner-Fraser M, Sauka-Spengler T: Assembling neural crest regulatory circuits into a gene regulatory network. Annu Rev Cell Dev Biol 2010, 26:581–603.PubMedView Article
- Swiers G, Patient R, Loose M: Genetic regulatory networks programming hematopoietic stem cells and erythroid lineage specification. Dev Biol 2006, 294:525–540.PubMedView Article
- Tu Q, Cameron RA, Worley KC, Gibbs RA, Davidson EH: Gene structure in the sea urchin Strongylocentrotus purpuratus based on transcriptome analysis. Genome Res 2012, 22:2079–2087.PubMedView Article
- Martin J, Wang Z: Next-generation transcriptome assembly. Nat Rev Genet 2011, 12:671–682.PubMedView Article
- Brown CT, Howe A, Zhang Q, Pyrkosz A, Brom T: A reference-free algorithm for computational normalization of shotgun sequencing data. [http://arxiv.org/abs/1203.4802]
- Miller JR, Koren S, Sutton G: Assembly algorithms for next-generation sequencing data. Genomics 2010, 95:315–327.PubMedView Article
- Salzberg SL, Yorke JA: Beware of mis-assembled genomes. Bioinformatics 2005, 21:4320–4321.PubMedView Article
- Ryan JF, Burton PM, Mazza ME, Kwong GK, Mullikin JC, Finnerty JR: The cnidarian-bilaterian ancestor possessed at least 56 homeoboxes: evidence from the starlet sea anemone, Nematostella vectensis . Genome Biol 2006, 7:R64.PubMedView Article
- Putnam NH, Srivastava M, Hellsten U, Dirks B, Chapman J, Salamov A, Terry A, Shapiro H, Lindquist E, Kapitonov VV, Jurka J, Genikhovich G, Grigoriev IV, Lucas SM, Steele RE, Finnerty JR, Technau U, Martindale MQ, Rokhsar DS: Sea anemone genome reveals ancestral eumetazoan gene repertoire and genomic organization. Science 2007, 317:86–94.PubMedView Article
- Zeng V, Villanueva KE, Ewen-Campen BS, Alwes F, Browne WE, Extavour CG: De novo assembly and characterization of a maternal and developmental transcriptome for the emerging model crustacean Parhyale hawaiensis . BMC Genomics 2011, 12:581.PubMedView Article
- Ewen-Campen B, Shaner N, Panfilio KA, Suzuki Y, Roth S, Extavour CG: The maternal and early embryonic transcriptome of the milkweed bug Oncopeltus fasciatus . BMC Genomics 2011, 12:61.PubMedView Article
- Conaco C, Neveu P, Zhou H, Arcila ML, Degnan SM, Degnan BM, Kosik KS: Transcriptome profiling of the demosponge Amphimedon queenslandica reveals genome-wide events that accompany major life cycle transitions. BMC Genomics 2012, 13:209.PubMedView Article
- Robertson G, Schein J, Chiu R, Corbett R, Field M, Jackman SD, Mungall K, Lee S, Okada HM, Qian JQ, Griffith M, Raymond A, Thiessen N, Cezard T, Butterfield YS, Newsome R, Chan SK, She R, Varhol R, Kamoh B, Prabhu A-L, Tam A, Zhao Y, Moore RA, Hirst M, Marra MA, Jones SJM, Hoodless PA, Birol I: De novo assembly and analysis of RNA-seq data. Nat Methods 2010, 7:909–912.PubMedView Article
- Li R, Zhu H, Ruan J, Qian W, Fang X, Shi Z, Li Y, Li S, Shan G, Kristiansen K, Li S, Yang H, Wang J, Wang J: De novo assembly of human genomes with massively parallel short read sequencing. Genome Res 2010, 20:265–272.PubMedView Article
- Zerbino DR, Birney E: Velvet: algorithms for de novo short read assembly using de Bruijn graphs. Genome Res 2008, 18:821–829.PubMedView Article
- Schulz MH, Zerbino DR, Vingron M, Birney E: Oases: robust de novo RNA-seq assembly across the dynamic range of expression levels. Bioinformatics 2012, 28:1086–1092.PubMedView Article
- Grabherr MG, Haas BJ, Yassour M, Levin JZ, Thompson DA, Amit I, Adiconis X, Fan L, Raychowdhury R, Zeng Q, Chen Z, Mauceli E, Hacohen N, Gnirke A, Rhind N, di Palma F, Birren BW, Nusbaum C, Lindblad-Toh K, Friedman N, Regev A: Full-length transcriptome assembly from RNA-seq data without a reference genome. Nat Biotechnol 2011, 29:644–652.PubMedView Article
- Zhao Q-Y, Wang Y, Kong Y-M, Luo D, Li X, Hao P: Optimizing de novo transcriptome assembly from short-read RNA-seq data: a comparative study. BMC Bioinformatics 2011,12(Suppl 14):S2.View Article
- Hansen KD, Brenner SE, Dudoit S: Biases in Illumina transcriptome sequencing caused by random hexamer priming. Nucleic Acids Res 2010, 38:e131-e131.PubMedView Article
- Venables WN, Ripley BD: Modern Applied Statistics with S. Berlin, Germany: Springer Verlag; 2002.View Article
- Ren X, Liu T, Dong J, Sun L, Yang J, Zhu Y, Jin Q: Evaluating de bruijn graph assemblers on 454 transcriptomic data. PLoS ONE 2012, 7:e51188.PubMedView Article
- Vijay N, Poelstra JW, Künstner A, Wolf JBW: Challenges and strategies in transcriptome assembly and differential gene expression quantification. A comprehensive in silico assessment of RNA-seq experiments. Mol Ecol 2013, 22:620–634.PubMedView Article
- Baker SC, Bauer SR, Beyer RP, Brenton JD, Bromley B, Burrill J, Causton H, Conley MP, Elespuru R, Fero M, Foy C, Fuscoe J, Gao X, Gerhold DL, Gilles P, Goodsaid F, Guo X, Hackett J, Hockett RD, Ikonomi P, Irizarry RA, Kawasaki ES, Kaysser-Kranich T, Kerr K, Kiser G, Koch WH, Lee KY, Liu C, Liu ZL, Lucas A, et al.: The External RNA Controls Consortium: a progress report. Nat Methods 2005, 2:731–734.PubMedView Article
- Devonshire AS, Elaswarapu R, Foy CA: Evaluation of external RNA controls for the standardisation of gene expression biomarker measurements. BMC Genomics 2010, 11:662.PubMedView Article
- External RNA Controls Consortium: Proposed methods for testing and selecting the ERCC external RNA controls. BMC Genomics 2005, 6:150.View Article
- Li B, Dewey CN: RSEM: accurate transcript quantification from RNA-seq data with or without a reference genome. BMC Bioinformatics 2011, 12:323.PubMedView Article
- Conesa A, Götz S, García-Gómez JM, Terol J, Talón M, Robles M: Blast2GO: a universal tool for annotation, visualization and analysis in functional genomics research. Bioinformatics 2005, 21:3674–3676.PubMedView Article
- Alexa A, Rahnenfuher J: topGO: enrichment analysis for gene ontology. [http://www.bioconductor.org/packages/2.11/bioc/html/topGO.html]
- Marlow H, Roettinger E, Boekhout M, Martindale MQ: Functional roles of Notch signaling in the cnidarian Nematostella vectensis . Dev Biol 2012, 362:295–308.PubMedView Article
- Röttinger E, Dahlin P, Martindale MQ: A framework for the establishment of a Cnidarian gene regulatory network for “endomesoderm” specification: the inputs of β-catenin/TCF signaling. PLoS Genet 2012, 8:e1003164.PubMedView Article
- Tarazona S, García-Alcalde F, Dopazo J, Ferrer A, Conesa A: Differential expression in RNA-seq: a matter of depth. Genome Res 2011, 21:2213–2223.PubMedView Article
- Yamada A, Pang K, Martindale MQ, Tochinai S: Surprisingly complex T-box gene complement in diploblastic metazoans. Evol Dev 2007, 9:220–230.PubMedView Article
- Adamska M, Larroux C, Adamski M, Green K, Lovas E, Koop D, Richards GS, Zwafink C, Degnan BM: Structure and expression of conserved Wnt pathway components in the demosponge Amphimedon queenslandica . Evol Dev 2010, 12:494–518.PubMedView Article
- Giresi PG: Chromatin profiles of human cells in health and disease using FAIRE. Ann Arbor, MI: ProQuest; 2012.
This article is published under license to BioMed Central Ltd. This is an Open Access article distributed under the terms of the Creative Commons Attribution License (http://creativecommons.org/licenses/by/2.0), which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.