Systematic identification of long noncoding RNAs expressed during zebrafish embryogenesis

Assembly of a high-confidence embryonic transcriptome

To systematically discover noncoding transcripts with potential functions during early vertebrate development, we performed large-scale cDNA sequencing experiments across zebrafish embryogenesis. We chose eight time-points that mark important developmental stages (Fig. 1A): (1) shortly after fertilization (two- to four-cell stage); (2) at the time when zygotic transcription of the genome is initiated (1000-cell stage); (3–5) during blastula and gastrula stages (dome, shield, and bud stages), when cell fates are specified and large-scale cell movements occur; and (6–8) at late embryonic and early larval developmental stages, when organs are forming (28 h post fertilization [hpf], 48 hpf, and 120 hpf) (see overview in Fig. 1A). Polyadenylated RNA was purified from approximately 1000 embryos per time-point and converted into cDNA libraries for strand-specific, paired-end 76 bp sequencing on Illumina's HiSeq platform (see Methods; Parkhomchuk et al. 2009; Levin et al. 2010). On average, we obtained about 200–300 million reads per stage (more than two billion reads in total) (Supplemental Table 1). Eighty-eight percent of the reads passed initial quality thresholds, of which ∼80% could be aligned to the latest assembly (Zv9) of the zebrafish genome sequence (Methods; Supplemental Table 1).

We assembled transcripts using a step-wise protocol (Methods; Fig. 1A). Briefly, we used TopHat (Trapnell et al. 2009) to align all reads per time-point, including those that span splice junctions. We then reconstructed transcripts using two assemblers—Cufflinks (Trapnell et al. 2010) and Scripture (Guttman et al. 2010), resulting in the assembly of a total number of 316,373 nonredundant transcript isoforms from 143,626 loci across all embryonic stages.

We defined a “high-confidence” set of 56,535 embryonic transcripts, following a similar strategy as described by Cabili et al. (2011). Specifically, we developed a filtering pipeline aimed at reducing the number of transcripts that might be erroneously assembled or below significance thresholds (Supplemental Fig. 1A). We first required a transcript be assembled at least twice: either identified by both assemblers or in at least two developmental stages. Next, we removed transcripts there were likely to be assembly artifacts or run-on fragments or that did not pass our high-confidence thresholds (Methods; Supplemental Fig. 1A). This resulted in a final set of 56,535 embryonic transcripts from 28,912 loci (on average, 1.95 transcripts per locus) (Supplemental Fig. 1B,C), of which 50,904 were multi-exonic and 5631 were single exon transcripts. We will henceforth refer to this set as the “embryonic transcriptome,” and all subsequent analyses are based on it.

The embryonic transcriptome has high coverage, quality, and depth

To estimate the quality and coverage of our embryonic transcriptome, we compared it to the current RefSeq and Ensembl gene annotations. Compared to the 15,175 zebrafish RefSeq genes, our embryonic transcriptome provides more than a threefold increase in the number of identified transcripts (56,535) from nearly twice as many loci (28,912), suggesting that the increase in the number of individual transcripts is due to both novel isoforms of known genes and novel loci (Fig. 1B). Notably, of the 13,942 RefSeq genes that are expressed (FPKM [fragments per kilobase of exon per million fragments mapped] > 0) during the stages covered by our data set, 90% (12,527/13,942) have transcript evidence (exonic overlapping transcripts) in our embryonic transcriptome, and 70% of those (8751/12,527) are identical to RefSeq isoforms (Supplemental Fig. 2A). In addition, 3532 of our transcripts are variants of known RefSeq genes (novel isoforms and partial transcripts), many of which extend the existing exon–intron structures with additional 5′ or 3′ exons (Supplemental Fig. 2A).

Compared to the most recent Ensembl Zv9 gene models (52,873 transcripts in 31,711 loci), our embryonic transcriptome is of similar size (Supplemental Fig. 2B). The Ensembl gene set integrates transcript annotations from several sources and includes RNA-seq transcript models built from a total of 376 million reads derived from embryonic, larval, and adult zebrafish (Sanger Institute). Thus, the comparable transcriptome size and larger read numbers in our purely embryonic RNA-seq experiments again suggest that our embryonic transcriptome is of high depth. Moreover, we have transcript evidence for ∼74% of Ensembl gene loci of comparable (>160 nt) transcript sizes, corresponding to ∼68% of our embryonically identified loci (Fig. 1B). This high degree of overlap provides independent confirmation for a large fraction of our transcriptome.

Compared to two recent RNA-seq–based transcriptome studies in zebrafish embryos (Aanes et al. 2011; Vesterlund et al. 2011), our data set is of significantly higher depth: a total of about 220 million (Vesterlund et al. 2011) and about 100 million (Aanes et al. 2011) mapped reads from four and six embryonic stages, respectively, versus about 1.5 billion mapped reads in our study (Supplemental Table 1). Moreover, we identify many more known and novel transcribed loci: about 4000 “novel transcribed regions” reported by Vesterlund et al. (2011) and Aanes et al. (2011) versus more than 9000 novel loci in our embryonic transcriptome with no previous annotations in RefSeq or Ensembl. This suggests that our embryonic transcriptome provides a highly comprehensive and more complete assembly than that previously available.

Our transcript assemblies are also consistent with chromatin marks known to be associated with promoters (Zhou et al. 2011) (see also below). The fraction of marked protein-coding loci of our embryonic transcriptome (44% for H3K4me3 only, 19% for H3K4me3 and H3K27me3) is nearly identical to the fraction of marked RefSeq loci (46% for H3K4me3 only, 16% for H3K4me3 and H3K27me3) (see Fig. 5). This suggests that (1) our embryonic transcriptome is of a quality comparable to RefSeq genes, and (2) many of our RNA-seq–based transcript structures contain complete 5′ ends.

Identification of a stringent set of embryonic lncRNAs

To identify mRNAs that exert their biological function as lncRNAs, we developed a highly stringent filtering pipeline aimed at removing transcripts with evidence for protein-coding potential (Methods; Fig. 2). We identified putative lncRNAs by considering their phylogenetic conservation across species, homology with known proteins and protein domains, and potential ORFs.

Four filters were used. First, we used PhyloCSF (phylogenetic coding substitution frequency; see Methods), to score the coding potential of transcripts using phylogenetic alignments (Lin et al. 2011). PhyloCSF exploits the fact that protein-coding sequences—but not lncRNAs and other sequences—tend to have a higher rate of synonymous versus nonsynonymous substitutions (Supplemental Figs. 3, 4A). We chose a PhyloCSF threshold of less than 20 because it retained the majority of RefSeq ncRNAs (Supplemental Figs. 3, 4A) but removed 96.2% of protein-coding RefSeq transcripts. This filter retained 4867 putative noncoding transcripts (Fig. 2).

Second, we removed transcripts that had similarity to known proteins or protein domains based on blastx, blastp, and HMMER (Pfam domains) (Eddy 2009). This filter retained 2531 putative noncoding transcripts (Fig. 2B). The excluded transcripts had not been captured by the PhyloCSF filter because they typically received low PhyloCSF scores due to poorly aligned sequences (complete branch lengths [CBLs] of zero) (Methods; Supplemental Fig. 4B).

Third, we removed any remaining transcript of uncertain coding potential by applying a maximal ORF filter. Consistent with the traditional cutoff for protein-coding transcripts (Okazaki et al. 2002), we excluded any transcript with a maximal ORF > 100 amino acids (aa). For transcripts that were not scored by PhyloCSF due to lacking sequence alignments (CBL = 0), we used a more stringent maximal ORF cutoff of 30 aa. The ORF filter retained 1301 transcripts.

Finally, to exclude potentially incomplete transcript structures, we removed any transcript that had sense exonic overlap with a protein-coding transcript. The resulting set contained 902 lncRNAs (mean PhyloCSF score of 5) (Fig. 2B; Supplemental Fig. 3).

Identification of antisense overlapping embryonic lncRNAs

Some putative noncoding transcripts had antisense exonic overlap with protein-coding genes. Examination of the range of PhyloCSF scores obtained for antisense strands of sense-coding transcripts revealed that transcripts with a high-scoring sense strand also tended to score relatively high on the antisense strand (Supplemental Fig. 4C). Thus, PhyloCSF scores of antisense exonic overlapping transcripts can be confounded by high coding potential on the opposite strand.

To address this issue and “rescue” noncoding antisense transcripts, we employed a modified filtering pipeline with four additional criteria (Fig. 2; for details, see Methods): (1) The putative noncoding transcript had a lower PhyloCSF score than the overlapping coding transcript; (2) its highest PhyloCSF score was obtained in the region of overlap (e.g., Supplemental Fig. 4D); (3) its PhyloCSF score was less than 300; and (4) the sense/antisense exonic overlap did not exceed 81% of the sense strand. This approach retained 231 multi-exonic antisense transcripts and resulted in a final stringent set of 1133 lncRNAs (Fig. 2B).

Genomic characterization of embryonic lncRNAs

According to their genomic location, our 1133 embryonic lncRNAs are partitioned into 397 lincRNAs without overlap with any genes, 184 intronic overlapping lncRNAs, and 566 antisense exonic overlapping lncRNAs (Fig. 3). Intronic overlapping lncRNAs are defined as loci that have no exon–exon overlap with another locus, i.e., there is no overlap between the mature lncRNA with exons of the overlapping locus. Intronic overlapping lncRNAs are in either sense or antisense orientation with respect to the overlapping gene and can be further partitioned into 105 intronic contained lncRNAs (incs; the lncRNA is contained within the transcribed region of another locus), 60 completely overlapping lncRNAs (concs; the other locus is contained within the transcribed region of the lncRNA locus), and 19 partially overlapping lncRNAs (poncs; neither incs nor poncs but with at least one exon of the lncRNA contained within an intron of another locus).

Some lncRNAs may function as precursors for the generation of sRNAs (ENCODE Project Consortium et al. 2007; Wilusz et al. 2008). To identify sRNA-precursor lncRNAs, we compared our lncRNA transcripts to a set of sRNAs present in 2-d-old zebrafish (Methods; Cifuentes et al. 2010). We identified 41 lncRNAs that appear to function as precursors for the production of miRNAs (16), snoRNAs (nine), or sRNAs of unknown categories (20) (Supplemental Table 2). Four lncRNAs of the latter category contained a vast number of sRNAs throughout the entire transcript. For example, the zebrafish ortholog of the abundant nuclear lncRNA MALAT1 (also called NEAT2) was cleaved throughout its transcript and gave rise to multiple sRNAs (Supplemental Fig. 5). Consistent with this observation, MALAT1 has previously been shown to be associated with Ago2 (also known as EIF2C2), a known component of the sRNA processing machinery (Weinmann et al. 2009). This analysis indicates that the large majority of our lncRNAs are not processed into sRNAs.

Zebrafish lncRNAs are shorter, less conserved and expressed at lower levels than are protein-coding genes

Previous studies in mammals have shown that lncRNAs are shorter, less conserved, and expressed at significantly lower levels than are protein-coding genes (Guttman et al. 2010; Cabili et al. 2011). To determine whether embryonic lncRNAs have similar features, we analyzed the structure, expression level, and conservation of our lncRNAs (Fig. 4). We found that zebrafish lncRNAs were on average about one-third of the length of protein-coding transcripts (mean length of 1113 nt for lncRNAs versus 3352 nt for coding transcripts) (Fig. 4Aa). Moreover, lncRNAs had fewer exons per transcript (about 2.8) than the average protein-coding gene (about 11) (Fig. 4Ab). These properties are comparable to the estimated transcript length and exon number of human lincRNAs (on average, ∼1 kb and 2.9 exons, respectively) (Cabili et al. 2011). Notably, zebrafish embryonic lncRNAs were expressed on average at about 10-fold lower levels than protein-coding genes (Fig. 4B), consistent with the low expression levels of their mammalian counterparts (Guttman et al. 2010; Cabili et al. 2011).

To assess the level of conservation of lncRNAs, we used the CBL score, a measure of the fraction of phylogenetic teleost alignments present over the region of interest (Methods). In agreement with signatures of conservation in mammalian lncRNAs (Ponjavic et al. 2007, 2009; Guttman et al. 2009, 2010; Ørom et al. 2010), a few lncRNAs were clearly conserved across fish species (Fig. 4C; for two conserved examples, see Supplemental Fig. 6A). However, the majority of zebrafish lncRNA loci had low CBL scores, indicating a lack of sequence alignments over many noncoding regions (Fig. 4C). The conservation of zebrafish lncRNAs as reflected by CBL scores was substantially lower than the conservation of protein-coding genes and was comparable to the conservation of intronic sequences (Supplemental Fig. 6B).

lncRNA genes carry chromatin marks associated with developmental regulators

To assess to which extent lncRNA genes carry chromatin marks that are known to be associated with protein-coding genes (Vastenhouw et al. 2010; Zhou et al. 2011), we performed chromatin immunoprecipitation assays in shield stage embryos followed by deep sequencing (ChIP-seq). We tested for the presence of trimethylated lysine 4 on histone 3 (H3K4me3), a known marker of promoters, and trimethylated lysine 27 on histone H3 (H3K27me3), a repressive histone modification. We restricted our analysis to lincRNAs and intronic overlapping lncRNAs since unambiguous assignment of marks to antisense exonic overlapping transcripts can be confounded by the overlapping genes.

Of all lncRNA promoter regions that were assessed, 29% were marked with H3K4me3 (both H3K4me3-only and H3K4me3/H3K27me3) (Fig. 5A). Notably, the fraction of H3K4me3-marked zebrafish lncRNA genes was similar to the 24% of human lincRNA genes that have a K4-K36 domain (Cabili et al. 2011), but was smaller than the fraction (63%) of marked zebrafish protein-coding genes (Fig. 5A).

To consider the possibility that the discrepancy between the fraction of H3K4me3-marked lncRNA and protein-coding loci could be due to the lower expression levels of lncRNA loci, we restricted our analysis to protein-coding genes expressed at shield stage and at expression levels similar to lncRNAs. Even under these conditions, the discrepancy between H3K4me3-positive noncoding (34%) and coding (74%) loci remained (Fig. 5B). This suggests that (1) the different expression levels of noncoding and protein-coding loci are not the primary cause of the different fractions of H3K4me3-marked loci, and (2) similarly to protein-coding genes (Vastenhouw et al. 2010), noncoding loci are marked with H3K4me3 largely independently of their expression status.

Interestingly, 7% of lincRNA and intronic overlapping lncRNA loci were marked by both H3K4me3 and H3K27me3 at shield stage (Fig. 5). Since Gene Ontology (GO)–term analysis of protein-coding genes marked with both H3K4me3 and H3K27me3 at shield stage revealed enrichment for developmental and regulatory functions (Supplemental Table 3), lncRNA loci may be important developmental regulators.

Nearest neighbor analysis of lncRNA genes

Previous studies have shown that mammalian lncRNAs are preferentially located next to genes with developmental functions (Dinger et al. 2008; Mercer et al. 2008; Guttman et al. 2009; Ponjavic et al. 2009; Cabili et al. 2011). We therefore analyzed the GO terms of genes that overlap with or are neighbors of zebrafish lncRNAs. We found significant enrichments (P < 0.05) of transcription factor activity, fate specification, and embryonic development and morphogenesis for genes that overlap with antisense exonic lncRNAs (Supplemental Fig. 7A; Supplemental Table 4) but not for neighbors of lincRNAs and intronic overlapping lncRNAs (Supplemental Table 4).

The mere physical proximity of lncRNAs and genes with developmental functions does not necessarily imply a functional link between the protein-coding gene and the lncRNA. For example, recent studies in the mouse did not detect a strong correlation between the expression levels of most lncRNAs and their neighbors (Guttman et al. 2011). Consistent with this study and with data from human lincRNAs (Cabili et al. 2011), we did not detect a higher degree of expression correlation for the majority of lncRNAs and their neighbors (or overlapping genes) than for protein–protein gene pairs or randomly assigned gene pairs (Supplemental Fig. 7B). The only exceptions were sense intronic overlapping lncRNAs, which tended to positively correlate in expression with the overlapping genes (Supplemental Fig. 7B). Such overlapping lncRNAs might resemble enhancer-associated lncRNAs (De Santa et al. 2010; Kim et al. 2010; Ørom et al. 2010; Wang et al. 2011).

To test whether lncRNA genes are preferentially located near protein-coding genes of certain evolutionary ages, we analyzed the phylostratographic classes (Domazet-Lošo and Tautz 2010) of genes that neighbor or overlap lncRNAs. Comparison with enrichments observed in a control set of protein–protein gene neighbors revealed no significant enrichment of particular evolutionary age groups for our lncRNA neighbors (Supplemental Fig. 7C).

Collectively, our analysis suggests that the neighbors of zebrafish lncRNAs belong to various classes of protein-coding genes of both ancient and more recent evolutionary origin and generally do not correlate in their expression with the neighboring lncRNAs.

Temporal expression profiles of lncRNAs

Our high-resolution time-series of RNA-seq experiments allowed us to follow the expression dynamics of lncRNAs and protein-coding genes as development proceeds. Comparison of independently clustered expression profiles of noncoding and protein-coding loci (Methods) revealed that both types of loci could be grouped into three broad classes (Fig. 6A): (1) loci whose transcripts were parentally supplied—these transcripts were present in the two- to four-cell-stage embryo (cleavage stages) and rapidly decayed after the first few hours of embryogenesis; (2) loci whose expression peaked during blastula and gastrula stages (dome, shield, and bud stages)—these transcripts were absent or only present at low levels during the early cleavage stages and were zygotically transcribed; and (3) loci that were only induced 1 d after fertilization during the process of organogenesis.

We discovered two differences between the expression patterns of protein-coding and noncoding loci. First, lncRNAs were more likely to be parentally supplied than were protein-coding mRNAs (see Fig. 6A). Any locus was classified as “parentally provided” for which at least 10% of its total expression across all eight embryonic stages was derived from the two- to four-cell stage. Of all transcripts present in our catalog, ∼42% of lncRNAs classified as parentally provided, compared with only ∼34% of protein-coding transcripts (Fisher's exact test, P < 10⁻⁰⁵). These observations suggest that parentally provided transcripts are specifically enriched in lncRNAs.

Second, the changes in a transcript's expression level between two consecutive stages were more pronounced for lncRNAs than for protein-coding genes. This observation suggests that lncRNAs have a more restricted temporal expression than do coding RNAs. To further test this hypothesis, we calculated a Shannon entropy-based specificity score per locus as a measure of expression level divergence during embryogenesis (Methods). All three classes of lncRNAs (lincRNAs, intronic overlapping, and antisense exonic overlapping lncRNAs) showed an increased temporal specificity compared with protein-coding genes (Fig. 6B). To rule out that this effect was caused by an increase in noise due to the lower expression levels of lncRNAs, we also sampled protein-coding loci from the same expression quantiles as lncRNAs (Methods). Although protein-coding loci that were expressed at low levels tended to be more restricted in time than were highly expressed protein-coding loci, they were significantly less restricted than were lncRNAs (P < 10⁻⁴) (Fig. 6B). Together, these analyses reveal high temporal specificity during development as a novel property of lncRNAs.

Assigning function through expression correlation

The lack of annotated features makes the assignment of functions to lncRNAs a more challenging task than for proteins. Therefore, functional predictions for mammalian lncRNAs have often been based on “guilt-by-association” analyses (Dinger et al. 2008; Guttman et al. 2009; Cabili et al. 2011). We therefore analyzed the correlation between the expression dynamics of each protein-coding gene with the expression dynamics of each lncRNA locus. We performed gene set enrichment analysis (GSEA) (Methods; Mootha et al. 2003; Subramanian et al. 2005; Guttman et al. 2009) and associated GO terms and lncRNAs. This analysis identified several groups of lncRNAs associated with protein-coding gene sets of distinct functional categories such as signaling (cluster 2), development (cluster 6), and cell cycle (cluster 8) (Fig. 6C; for a complete list of enriched GO terms in each cluster, see also Supplemental Table 5). Interestingly, about one-third of our lncRNAs were associated with clusters enriched in developmental functions (clusters 4–7). These results indicate that many embryonic lncRNAs are putative developmental regulators.

LncRNAs show tissue-specific and subcellularly restricted expression patterns

To determine whether embryonic lncRNAs were expressed in specific tissues, we performed RNA in situ hybridization for a selected set of lncRNAs. Thirty-two lncRNAs were amplified from cDNA. While the majority of lncRNAs did not reveal strong or tissue-specific expression (Supplemental Fig. 8), several lncRNAs showed both expression in specific cell types and distinct subcellular RNA localization patterns (Fig. 7 and data not shown). Examples for cell type–specific expression patterns included lncRNAs that were loaded into the fertilized embryo by cytoplasmic streaming (e.g., hoxAa_lncRNA) (Fig. 7Ai), a lncRNA expressed in developing somites (myo18a-lncRNA) (Fig. 7Aii), and lncRNAs with distinct expression patterns in the developing nervous system (Fig. 7Aiii,iv).

Several zebrafish lncRNAs showed distinct subcellular localization patterns, supporting and extending previous localization studies in the mouse brain (Mercer et al. 2008; Ponjavic et al. 2009). For example, we observed nuclear enrichment of some lncRNAs in early cleavage stage embryos, including chromatin association in mitotically dividing cells (e.g., hoxAa_lncRNA) (Fig. 7Bi). Other lncRNAs such as mprip_lncRNA were found to accumulate at the cytoplasmic side of yolk syncytial layer nuclei at the bud stage (Fig. 7Bii). A particular striking example for a subcellularly localized lncRNA was myo18a-lncRNA, which was enriched specifically at the myoseptum, the boundary where myotubes of adjacent somites meet (Fig. 7Aii,Biii). The myo18a-lncRNA localization pattern was distinct from many muscle-specific mRNAs that are ubiquitously expressed in somites (e.g., myzh1.1 [myosin heavy chain]), but resembled the mRNA localization of the protein-coding gene dystrophin, a known marker of the myoseptum (Fig. 7Biii, left panel; Bassett 2003). Dystrophin is a key component of the protein complex that connects the cytoskeleton of the muscle fiber to the extracellular matrix, and its deficiency causes severe myopathies (Koenig et al. 1988). Intriguingly, a potential function of myo18a_lncRNA in cell-–cell contact formation was supported by our expression-based GSEA approach, which associated myo18a_lncRNA with functions in “cell adhesion” and “structural molecule activity” (Fig. 6C, cluster 3). Collectively, these results reveal that several embryonic lncRNAs are expressed not only in specific tissues but also in specific subcellular domains.

RNA-seq of embryonic time course

Wild-type zebrafish embryos (TLAB) were staged according to standard procedures. About 1000 embryos were collected per stage (two to four cell, 1000 cell, dome, shield, bud, 28 hpf, 48 hpf, and 120 hpf) within a tight time window of ∼10 min; it was ensured that all embryos were at the same developmental stage. Total RNA was isolated using the standard TRIzol (Invitrogen) protocol. Genomic DNA was removed by DNase treatment and confirmed by qPCR assay. Two rounds of PolyA⁺-RNA purification were performed for each sample, using the PolyA(Purist)-MAG kit (Ambion). The quality of the RNA and lack of contaminating ribosomal RNA were confirmed using the Agilent 2100 Bioanalyzer. Strand-specific libraries for 76-bp paired-end sequencing were prepared according to a modified UTP-method (Parkhomchuk et al. 2009), as detailed by Levin et al. (2010). Libraries were sequenced on the GA-analyzer (shield stage library) and on the Illumina HiSeq 2000 (all stages), at a depth of 200–300 million reads per library (for statistics on read counts, see Supplemental Table 1).

Transcriptome assembly

RNA-seq–derived reads were aligned independently for each developmental stage with TopHat (version 1.2.1) (Trapnell et al. 2009). To aid these alignments, all known transcript annotations (Ensembl, RefSeq, and mRNAs from UCSC danRer7 [Zv9]) were pooled and used as an additional junction set (AJS) for each TopHat run. The junction outputs from individual stage-specific TopHat runs were pooled and added to the AJS (augmented AJS) to allow TopHat to use junction information from all stages. TopHat was rerun on each of the stages using the augmented AJS. The output of this second run comprised the final alignment and junction set for transcript assembly.

Transcriptomes were assembled with two different assemblers: Cufflinks (version 1.0.3) (Trapnell et al. 2010) and Scripture (version R4) (Guttman et al. 2010). The resulting transcripts were pooled, and a transcript was only considered reliable if it had support either from both assemblers or from at least two stages. Transcripts <160 bp were excluded, as these were most likely sequencing or assembly artifacts.

Multi-exonic transcripts

Multi-exonic transcripts were merged with Cuffcompare, and all transcripts classified as repeat were discarded. Scripture's strategy is to call all possible isoforms, including some that are most likely wrong and have no Cufflinks support. Therefore, all Scripture-only isoforms lacking Cufflinks support were excluded whenever Cufflinks had an assembled transcript for this locus. Furthermore, any di-exonic antisense transcript only supported by Scripture was removed since these transcripts are likely artifacts due to Scripture's lack of strand-aware library support.

Single exon transcripts

Single exon transcripts were subjected to additional scrutiny: They had to be significantly enriched in read coverage by Scripture (multiple testing corrected P < 0.01) (Guttman et al. 2010) and had to have at least one supporting transfrag from Cufflinks. Cufflinks uses library strand information and can therefore correctly assign the strand for single exon transcripts, while Scripture relies only on splice junctions and therefore cannot determine the strand-orientation of single exons. Transcripts classified by Cuffcompare (Trapnell et al. 2010) as contained (c), exon–intron fragment (e), exonic overlap (o), RNA pol II run-on (p), and repeat (r) were removed. Moreover, any single exon that was within a range of 500 bp in the sense direction relative to a multi-exonic transcript was removed.

Finally, single exon and multi-exonic transcripts were merged with Cuffcompare, discarding all contained and redundant isoforms.

ncRNA classification

Classification of each transcript as either coding or noncoding was determined using a step-wise filtering pipeline. First, all candidates were scored with PhyloCSF (Lin et al. 2011) to determine their coding potential (see PhyloCSF section). All transcripts that scored less than 20 were retained as potential noncoding candidates, and transcripts with PhyloCSF scores greater than 50 were considered proteins. The remaining transcripts (20 < PhyloCSF < 50) were initially classified as an ambiguous “gray” set. Second, the putatively noncoding transcripts and the “gray” set transcripts were repeat-masked and subjected to blastx, blastp, and HMMER (versus Pfam-A and Pfam-B) (Eddy 2009). For blastp and HMMER, the transcripts were translated (stop to stop codon, due to possible incomplete assemblies that could result in incomplete ORFs lacking the ATG start codon) in all three sense frames. Any transcript with an E-value less than 10⁻⁴ in any of the three search algorithms was considered as protein-coding.

Not all candidates were alignable to regions in the other four fish species. PhyloCSF-based coding potential predictions are less reliable for transcripts with no alignments over their entire region (see Supplemental Fig. 4B). Therefore, a maximal ORF cutoff was imposed. For candidates without alignments (CBL = 0, see Comparative Genomics Analysis of Conservation and Coding Potential) this cutoff was set to 30 aa, for all remaining transcripts (CBL > 0) it was set to 100 aa.

Finally, transcripts were removed that have exonic sense overlap with either Ensembl or RefSeq protein-coding genes or with the protein-coding gene set of the embryonic transcriptome.

Antisense rescue pipeline

PhyloCSF scores of antisense transcripts can be confounded by high-scoring protein-coding genes on the opposite strand, necessitating an alternative strategy for this set. The antisense-rescue pipeline was similar to the general pipeline (see above), with the exception that the PhyloCSF threshold was set to 300. In addition, the highest scoring region (HSR) for the putative antisense lncRNA had to have overlap with a protein-coding transcript on the opposite strand, and the PhyloCSF score for the protein-coding gene had to be higher than the PhyloCSF score for the lncRNA (see Supplemental Fig. 4D). Finally, after manual inspection, a threshold of maximal 81% was set for the sense/antisense exonic overlap with the protein-coding gene to remove a small number of likely artifactual transcripts that had substantial overlap and likely stemmed from errors in strand-calling during assembly.

Classification of lncRNAs

The resulting set of lncRNAs was subdivided into three categories: (1) lncRNAs without any overlap with other loci classify as lincRNAs (intergenic lncRNAs); (2) lncRNAs with intronic overlap are expressed from loci that have overlap (exon–intron or intron–intron but not exon–exon) with another transcribed locus (i.e., there is no overlap between the mature lncRNA with exons of the overlapping locus). They can be in either sense or antisense orientation with respect to the overlapping gene and can be further partitioned into intronic contained lncRNAs (incs; the lncRNA is contained within the transcribed region of another locus), completely overlapping lncRNAs (concs; the other locus is contained within the transcribed region of the lncRNA locus), and partially overlapping lncRNAs (poncs; neither incs nor poncs, but with at least one exon of the lncRNA contained within an intron of another locus). (3) Exonic antisense overlapping lncRNAs have exonic overlap with an exon of a protein-coding transcript on the opposite strand.

Comparative genomics analysis of conservation and coding potential

The RNA-seq transcripts were analyzed in MULTIZ whole-genome alignments of zebrafish with four other fish species (Tetraodon, Fugu, Stickleback, Medaka), generated by UCSC (http://genome.ucsc.edu/cgi-bin/hgTrackUi?db=danRer7&g=cons8way).

sRNA analysis

sRNAs expressed in 2-d-old wild-type zebrafish larvae (two biological replicates) were obtained from Cifuentes et al. (2010) and mapped to Zv9 using Bowtie (Langmead et al. 2009). The number of sRNAs overlapping lncRNA loci were counted. Transcripts with at least five uniquely mapped overlapped sRNAs were annotated according to known sRNA classes (miRNA precursor, snoRNA precursor, MALAT1-like transcripts, transcripts of unknown sRNA types).

Expression analysis

The expression level of each locus was assessed using Cuffdiff (Trapnell et al. 2010) in its time-series mode with upper quantile normalization. To visualize developmental expression profiles via heatmaps, expression levels were normalized to get relative expression levels over the developmental time-course (sum to an expression of one over all stages). LncRNAs and protein-coding loci were clustered separately using k-means (10 clusters) with a distance matrix constructed from the Pearson correlation.

The temporal specificity score over N time-points (N = 8 embryonic stages) was defined as 1 − H(g)/log2(N), where H(g) is the Shannon entropy expressed in bits of the expression vector of gene g. To compare lncRNA loci to protein-coding loci of similar expression levels, the expression levels of each locus were summed over all time-points and sorted into 100 quantiles. The Shannon entropy was then calculated for samples of proteins equal in number to the lncRNAs and from the same expression bins (10,000 repetitions to estimate the dispersion and to calculate the P-value).

In situ expression analysis

To analyze the expression pattern and localization of RNAs, 300- to 800-nt-long partial or full sequences of mRNAs were amplified from cDNA by PCR and cloned into the pSC-vector (Strataclone) according to standard procedures (primer sequences are available upon request). Thirty-two lncRNAs were cloned and tested by in situ hybridization experiments for their expression patterns from shortly after fertilization to 2-d-old larvae. Digoxigenin (DIG)-labeled antisense RNA probes were generated by in vitro transcription with T3 or T7 RNA polymerases, using plasmid-encoded T3 or T7 polymerase binding sites. In situ hybridization of zebrafish embryos of different embryonic stages was performed according to standard procedures (Thisse and Thisse 2008), using immunohistochemical detection of the DIG-labeled RNA–RNA hybrids by an anti-DIG Alkaline-Phosphatase coupled antibody, followed by nonfluorescent detection with BCIP/NBT. DNA was visualized by incubation of stained embryos with the DNA-dye OliGreen (Invitrogen, used at 1/400 in PBST). Images were processed with Photoshop and ImageJ.

Chromatin mark analysis: ChIP-seq for H3K4me3 and H3K27me3 at shield stage

ChIP was performed as previously described (Vastenhouw et al. 2010). Antibodies used were H3K4me3 (Millipore no. 07-473) and H3K27me3 (Millipore no. 07-449). For analysis on the Illumina Hi-Seq platform, sequencing libraries were prepared according to Illumina protocols.

Peak calling for chromatin marks was done using Scripture's ChIP-seq module (Guttman et al. 2010). This module scans fixed-size windows across the genome and computes read coverage and a multiple hypothesis corrected P-value for the observed coverage. For both H3K4me3 and H3K27me3, 500- and 1000-bp windows were scanned to account for both short regions with high read coverage and for larger regions with lower read coverage. All windows that were covered at a significant level (P < 0.01) were merged into “peaks.” The ends of the peaks were finally trimmed until coverage at the ends is at least the average peak coverage. To account for systematic biases—e.g., due to open chromatin—peaks were filtered using input genomic DNA sequence by requiring that every peak called contained a 500-bp window with a library score at least threefold higher than the input genomic sequence score.

Peaks were intersected with promoter regions (±500 bp relative to the transcriptional start site [TSS]) of our transcripts. To obtain protein-coding loci of similar expression levels as lncRNA loci at shield stage, the same strategy was used as for expression analysis (see above), except that ranking of protein-coding loci was based exclusively on their expression levels at shield stage.

Gene set enrichment analysis

The expression level of each lncRNA locus was correlated with all protein-coding loci, similar to (Guttman et al. 2009). For each lncRNA locus, a list of correlation-based ranked protein-coding loci was constructed and subjected to GSEA (Mootha et al. 2003; Subramanian et al. 2005). An association matrix between lncRNA loci and GO terms was constructed, using a false-discovery rate threshold of 0.01. Rows (lncRNA loci) and columns (GO terms) were clustered (k-means, 10 clusters), resulting in distinct subsets of lncRNAs associated with functional GO terms. To determine the enrichment level of positively associated GO terms for each cluster with respect to other clusters, positively correlated GO terms were ranked according to a binominal test.

Nearest neighbor analysis

For each lincRNA locus the nearest protein-coding neighbor within <10 kb was identified. For antisense overlapping and intronic overlapping lncRNAs, overlapping gene(s) were identified. This resulted in a list of lncRNA loci/protein-coding loci pairs. Similar to the method described by Cabili et al. (2011), Pearson correlation was used to explore the expression-based relationship between these pairs. The results were subdivided based on (1) the class of the lncRNA and (2) the orientation of the lncRNA locus relative to the neighbor/overlapping protein-coding locus. The list of pairs (lncRNA loci/protein-coding loci) also formed the basis for GO term enrichment analysis using GOstat (Beissbarth and Speed 2004) and phylostratographic analysis of nearest neighbors (see below).

Phylostratographic analysis of nearest neighbors

Phylostratographic classes for zebrafish genes were obtained from Domazet-Lošo and Tautz (2010). LncRNA neighboring/overlapping protein-coding loci (see above) were tested for enrichment in certain phylostratographic classes by a sampling procedure that used the protein population as a null model.