Identification of Novel Transcribed Regions in Zebrafish (Danio rerio) Using RNA-Sequencing

Zebrafish (Danio rerio) has emerged as a model organism to investigate vertebrate development and human genetic diseases. However, the zebrafish genome annotation is still ongoing and incomplete, and there are still new gene transcripts to be found. With the introduction of massive parallel sequencing, whole transcriptome studies became possible. In the present study, we aimed to discover novel transcribed regions (NTRs) using developmental transcriptome data from RNA sequencing. In order to achieve this, we developed an in-house bioinformatics pipeline for NTR discovery. Using the pipeline, we detected 152 putative NTRs that at the time of discovery were not annotated in Ensembl and NCBI gene database. Four randomly selected NTRs were successfully validated using RT-PCR, and expression profiles of 10 randomly selected NTRs were evaluated using qRT-PCR. The identification of these 152 NTRs provide new information for zebrafish genome annotation as well as new candidates for studies of zebrafish gene function.


Introduction
Transcriptome analysis has become a key tool for understanding functional roles of genes involved in a variety of biological processes including early development [1].In 2005 Mathavan et al. [2] published a genome-wide microarray analysis of the embryonic zebrafish transcriptome using 12 different embryonic time points.The study revealed a highly dynamic and diverse transcriptional profile during embryogenesis and identified a previously unknown set of very early genes transcribed prior to the mid-blastula transition (MBT).Over the last decade, technologies for transcriptome analysis have improved dramatically, and high throughput RNA sequencing technologies (RNA-Seq) have revolutionized transcriptomics by detailed examination of cellular transcriptomes with high sensitivity, high dynamic range, and expression at a single-cell resolution [3,4].Unlike microarray expression analysis, transcriptome profiling using RNA-Seq may allow for discovery of novel transcribed regions (NTRs) since it is not limited by the availability of reference information.
Zebrafish is a valuable vertebrate model organism for studies in developmental biology [5] and for functional characterization of human disease genes, especially where the functional tissue such as brain is not readily available in human [6].However, the zebrafish genome annotation is not complete and still ongoing.Therefore, zebrafish sequencing data provide an opportunity for NTR discovery.Several whole transcriptome analyses using RNA-Seq have been performed in zebrafish [7][8][9][10][11].Besides unraveling expression profiles of transcriptome dynamics during early embryonic stages, a number of NTRs in annotated and unannotated regions of the zebrafish genome have also been described [7,8,12].The aim of this study is to discover novel transcribed regions (NTRs) using developmental transcriptome data from RNA sequencing.We have identified 152 putative NTRs that at the time of discovery were not annotated in Ensembl and NCBI gene database using an in-house bioinformatics pipeline to systematically in zebrafish early development by reanalyzing our previously obtained RNA-Seq data [8].

RNA preparation and sequencing
Embryo collection, RNA extraction and RNA-Seq have been described previously [8].All sequence data (accession numbers ERP000635 and PRJEB9889) are available at the European Nucleotide Archive (http://www.ebi.ac.uk/ena/) website.In summary, total RNA was extracted from four developmental stages (1-cell (0.75 hpf), 16-cell (1.5 hpf), 512-cells (2.75 hpf) and 50% epiboly (5.25 hpf)) and subsequently rRNA depleted using RiboMinus ™ Eukaryote Kit for RNA-Seq (Life Technologies).We performed in total three runs of RNA-Seq on the SOLiD System Sequencing platform (ABI, Applied Biosystems), with two biological replicate runs and one technical replicate run.The libraries were sequenced using 50 base pair (bp) reads.The research protocol was approved by the local Ethical Board, Stockholms Norra Djurförsöksetiska Nämnd (application number N413/11, N230/10 and S170/08).

Reads mapping and filtering
RNA-Seq data were analyzed using Tophat v2.0.4 [13] that applies Bowtie v.0.12.8 to handle color space reads generated by the SOLiD (ABI, Applied Biosystems).Briefly, after quality check [8] the sequencing reads from three runs at four stages were individually aligned to the zebrafish reference genome danRer7/Zv9 assembly from Ensembl.Only uniquely mapped reads were used for NTR discovery.Read filtering was performed using SAMtools [14].

Detection of putative NTRs
NTRs were detected using in-house pipeline composed of BEDTools [15] modules (mergeBed, slopBed, intersectBed) and customized scripts.The pipeline conducts two types of tasks, fragment construction based on the uniquely mapped RNA-Seq reads and formation of clusters (Fig 1 ).A fragment is defined as a number of adjacent reads overlapped by one or more bp at ends on the same strand or a read with detected splicing junction sites.A cluster is a group of linked fragments on the same strand.Output of the pipeline could be controlled by two parameters, D1 and D2.D1 is the distance between a cluster and any annotated transcript, while D2 is the distance between two adjacent clusters (Fig 1A ).
In order to find independent NTRs, fragments with a distance D1 away from any known transcripts on the same strand at either end were excluded from clustering since they could be additional parts of those known transcripts.Ensembl genome annotation release 79 was used as a primary reference to detect putative NTRs.In addition, we also checked all detected NTRs against NCBI zebrafish (Danio rerio) annotation releases 103 and 104 for further known transcript filtering.
For clustering, any two fragments with a distance less than D2 on the same strand were considered as linked.A series of linked fragments formed a cluster.A cluster could be defined as a putative NTR if 1) it included at least two splicing junction sites detected by TopHat; 2) it was not annotated in the databases mentioned above.

In silico confirmation of putative NTRs
A variety of in silico methods were used to confirm the presence of putative NTRs.Besides using junction information to clarify hypothetical exon boundaries within clusters, we also used gene models provided by GENSCAN [16] and open reading frame (ORF) predicted using FGENESH software (http://www.softberry.com) to confirm the structure of potential protein- coding genes.Zebrafish CAGE data [11] were used as additional supporting evidence.Moreover, we also looked for conservation of NTRs using human proteins mapped by tBLASTn and RefSeq gene of other species, which were download from the UCSC genome browser (http://hgdownload.soe.ucsc.edu/goldenPath/danRer7/database/xenoRefGene.txt.gz,http:// hgdownload.soe.ucsc.edu/goldenPath/danRer7/database/blastHg18KG.txt.gz).

Experimental validation of putative NTRs
Total RNA was extracted from zebrafish embryos at different developmental stages using Trireagent (Sigma Aldrich).cDNA was synthesized using SuperScript III First-Strand Synthesis SuperMix for reverse transcription polymerase chain reaction (RT-PCR) according to the manufacturer's protocol (Invitrogen, Life Technologies).Half a microgram of the total RNA was taken as input material.HotStarTaq plus DNA polymerase (Qiagen) was used in combination with specific PCR primers (S1 Table ) to amplify regions of the putative transcripts.The fragments were subsequently cloned into pCRII-TOPO vector (Invitrogen, Life Technologies) and validated using Sanger sequencing (Eurofin Genomics, Ebersberg, Germany).cDNA from 50% epiboly was used as template.The primer for validation was designed by using Primer3 [17,18], and they were located at the first and the last exons of each predicted isoform.

Expression evaluation of putative NTRs
The expression profiles of the discovered 152 NTRs were first evaluated based on the RNA-Seq from the biological replicate samples at all four studied developmental stages.The Cuffdiff module from Cufflinks 2.0.0 [19] was employed to estimate expression levels in terms of fragments per kilobase of transcript per million mapped reads (FPKM).Gplots package in R (http://CRAN.R-project.org/package=gplots) was employed to explore the expression profiles of the 152 NTRs.The expression levels (FPKM) were scaled in each NTR using heatmap2 function.Subsequently, 10 NTRs were randomly selected for validation using quantitative reverse transcription polymerase chain reaction (qRT-PCR).
For the qRT-PCR validation, RNA extraction and cDNA synthesis were performed as described above.The developmental stages used in cDNA synthesis were 1-cell, 16-cell, 512-cell and 50% epiboly.The primers (S1 Table ) were designed to span the exon splicing junctions, 150-200 bp in length, by using Primer3 [17,18].The gene expression analysis was performed in duplicates on three biological replicate samples of these four developmental stages.Expression of bactin2 in each developmental stage was used as the control.Fast SYBR Green Master Mix (Thermo Fisher Scientific) was used for qRT-PCR according to manufacturers protocol and the experiments were run on the Applied Biosystems 7500 Real-Time PCR system.All amplicons were validated using Sanger sequencing.

RNA-Seq data mapping
Each of the three runs of previously obtained RNA-Seq data from four zebrafish embryonic developmental stages (1-cell, 16-cell, 512-cell and 50% epiboly) [8] was realigned individually using TopHat v2.0.4 [13] against the zebrafish reference genome Zv9.The realignment was done using TopHat default parameters corresponding to SOLiD sequencing data and no genome annotation was used.In total, we obtained more than 200 million raw sequencing reads for each developmental stage from all three runs.Table 1 shows the raw read amounts of the four developmental stages in each run and the mapping rates for each stage based on the total read numbers.About 62-76 million reads could be mapped to the zebrafish reference genome Zv9, accounting for 26%-34% of the total reads from different stages.At each developmental stage, 37-55 million reads were uniquely mapped and subsequently used for NTR discovery.

Discovery of NTRs
The workflow of NTR discovery is described in Fig 1B .In order to increase the read coverage and depth, we pooled uniquely mapping reads from all four stages and all three runs.With the pooled reads (about 166 million in total) as the input, we obtained 487937 fragments with either overlapped reads and/or reads containing splicing junction sites.Among them, 254508 fragments were located in annotated transcripts regions (Ensembl genome annotation release 79).To avoid potential extensions of annotated transcripts, we further filtered out 42120 fragments with D1 1 kb, i.e. 1 kb or less away from any known transcripts.The remaining 191309 fragments formed 60194 clusters with separating distances less or equal to 5 kb between fragments (D2 to 5kb).We calculated numbers of splicing junction within each cluster and excluded those clusters without or having only one splicing junction site to avoid potential random error (Fig 1B).After this filtering step, we obtained 648 NTRs (S2 Table ).Among the 648 NTRs, 449 were predicted by GENSCAN and 105 had putative homologues in mammals.
During the course of the study, the zebrafish genome annotations have been continuously updated.Therefore, some of the NTRs predicted based on Ensembl genome annotation release 79 were subsequently annotated in NCBI zebrafish annotation release 103, validating our approach for those genes.However, even with the exclusion of these newly annotated regions, there were 180 NTRs remaining unannotated (Fig 1B).Among those NTRs, 28 were annotated in NCBI zebrafish annotation release 104 (Zv10 zebrafish reference genome) currently.Therefore, as a final result, we have identified 152 NTRs that had not been previously annotated in the databases (Table 2).These NTRs are distributed on all 25 chromosomes in the zebrafish genome (Fig 2 ).

Validation of putative NTRs
Before experimental validation, structures of NTRs built based on RNA-Seq data were checked in silico to investigate overlap with any predicted open reading frames (ORF) or potential gene models.The expression of NTRs was validated using RT-PCR with pooled RNA samples from 50% epiboly.We used Sanger sequencing to confirm the sequences of amplified PCR products.
Four NTRs, NTR50, NTR88, NTR103, and NTR145 were randomly selected for experimental validation.In general, the validation results showed high similarity with the predicted

Expression of NTRs
The expression of the discovered 152 NTRs was evaluated using RNA-Seq data from the biological replicates run using the Cufflinks.

RNA-Seq analysis and pipeline development
In our previous study [8] we used Bioscope that is developed specifically for color space reads from the Applied Biosystems Inc. (ABI).Bioscope could provide 65-73% mapping rates.However, it provides junction information only for known transcripts, not for unannotated regions.Thus, for NTR discovery, we chose TopHat because it is capable of detecting novel splice sites, even though having a lower mapping rate compared with Bioscope.TopHat uses an efficient read-mapping algorithm designed to align reads from a RNA-Seq experiment to a reference genome without relying on known splice sites [13].
For NTR discovery, we did not use Cufflinks package that focuses on finding new isoforms of annotated transcripts as a previous study did [12].In our study, we aimed to find NTRs that were not annotated.Although Cufflinks could provide information of transcribed fragments in intragenic regions, the program would not consider the separated fragments as linked transcript even if they are located close by to each other without overlapping.Our pipeline is more flexible handling these fragments.It assesses those reads with adjustable parameters D1 and D2 (Fig 1A).The potential random errors will be restricted by numbers of splicing junction in a putative NTR and other predication methods.
There are two input parameters, D1 and D2, in our pipeline for NTR discovery (Fig 1A).The values we used for D1 and D2 were arbitrary.D1 was used to avoid potential extension to an annotated gene by additional exon/s.We set D1 to 1 kb, so that any fragments located at most 1 kb away from any known transcripts were excluded.D1 could be increased to reduce the probability for fragments being detected as extensions of known transcripts, but with a larger value of D1 we could risk a failure to observe short transcripts near known transcripts.We used D2 to separate any two clusters on the same strand.Based on Ensembl genome annotation version 79, the mean of intronic lengths in the zebrafish genome is approximated 3 kb (standard deviation, SD 8.2 kb) and the mean of intragenic lengths is 64 kb (SD 112 kb).Since sizes of both intragenic and intronic regions vary to a high extent, it is difficult to choose a generic value suitable to all transcripts.The results we present were based on a setting of D2 to 5 kb (Fig 1B).However, we also tested a setting of D2 to 10 kb, which resulted in a slightly reduction of the total cluster number, from 60194 to 46376 (S5 Fig) .The discovery of a larger number of putative NTRs (672 vs. 648) with the setting of D2 to 10 kb was due to less clusters with single fragments being removed (S5 Fig) .Increasing D2 could cause fusion of two putative NTRs, while decreasing it might risk getting more incomplete NTRs.However, the selected settings of D2 to 5 kb worked well based on the outcome of both in silico and experimental validations.

NTRs discovery and validation
We pooled RNA-Seq data from all three runs and four developmental stages to enhance detection of NTR at discovery step since they were all from early developmental stages, but we validated the presence of NTRs at individual developmental stage by independent experiments separately.
We excluded multiple mapped reads to avoid complexity they could introduce in both fragment formation and clustering steps, although by doing so we risked losing information.However, by using in silico and/or experimental validations the influence of multiple mapped reads on the structures of predicted NTRs could be compensated for.We excluded 2448 NTRs with only one fragment over 50 bp without junction observed (Fig 1B).These single fragment NTRs could be candidates for novel noncoding gene transcripts, however this is beyond the scope of this study.Moreover, we also excluded 1461 clusters containing only one splicing junction (Fig 1B).As consequence, the number of NTRs we report here is smaller compared to those from previous studies [7,8,12], however the NTRs discovered in the present study could be more reliable because they were subject to several filtering steps with stricter criteria.On the other hand, most previously reported NTRs also included a large number of isoforms for individual NTRs.Therefore, the numbers could be larger as well.
In this study, we randomly selected and successfully validated 4 NTRs by using RT-PCR.Interestingly, among these validated NTRs, 3 NTRs (Fig 3B , 3C and 3D) contained a very large intron-like structure in each, which was detected based on junction information and confirmed by validation.The large intron-like structure could split those NTRs into two parts by D2 without junction information.Therefore, the junction information was very valuable.
Furthermore, validation also rendered evidence of non-predicted transcripts with alternative splicing sites (Fig 3C and 3D).Further validations are needed to clarify boundaries at both 5' and 3' ends of the discovered NTRs.

Regulation of discovered NTRs during development
In this study we have generated a list of NTRs that are putative candidate genes that may provide novel information on early embryo development, as well as being novel zebrafish genes with homology to human or other mammalian genes.
The majority of the discovered NTRs were expressed at low levels in early developmental stages (S1 Transcripts originating from the zygote genome have been shown to encode for factors involved in biological processes such as differentiation, pattern formation and cell morphology among others and thus some of the 152 NTRs may fall into these categories [7,8,12].For example, the expression of NTR57 measured by two platforms (RNA-Seq and qRT-PCR) showed a similar pattern, down regulated at 16-cell, upregulated at 512-cell stage, and further downregulated at 50% epiboly, suggesting that this NTR may be involved in MBT and deactivated after MBT.However, the standard error of NTR57 gene expression in 512-cell is very large among the biological replicates (Fig 4B).This standard error may reflect a biological variance between different individual female zebrafish, or it may reflect embryo developmental competence differences.Since the embryos were sampled at MBT (512-cell stage) it is not possible to get information on whether a high or low relative expression of NTR57 could be associated with any change in developmental competence or embryo survival past the 512-cell stage.It would be of particular interest to investigate the NTRs with relative high expression values, such as NTR157, NTR8, NTR61 and NTR94 in future studies.The zebrafish model system benefits from the continuous discovery of novel genes and more gene orthologues with candidate human disease genes.Although the putative biological relevance of the NTRs discovered is difficult to determine based only on the information obtained within this study, the validated transcripts give some indication on the biological relevance of the NTR candidates.Future studies will be needed in order to determine the function and characteristics of each of these.

Conclusion
With the increasing number of RNA-Seq data from a large variety of species, there is a vast amount of novel information that can be found using relatively simple in-house bioinformatics pipelines.We here described the development of an in-house bioinformatics pipeline for NTR discovery based on RNA-Seq data sets.Using this pipeline we discovered 152 putative NTRs that had not been previously annotated.Four randomly selected NTRs were successfully validated experimentally using RT-PCR.Using qRT-PCR we showed that the expression levels of 10 NTRs varied during the four early developmental stages investigated, suggesting that these NTRs may have a specific role in zebrafish early development.However, the characterization and function of these NTRs was not within the scope of present study and will require further investigation.
The 152 discovered NTRs provide new information for zebrafish genome characterization and zebrafish developmental studies as well as new gene candidates for developmental studies, thus further increasing the value of the zebrafish model system for the scientific community.

Fig 1 .
Fig 1. Systematic workflow for the identification of NTRs using RNA-Seq datasets.A. RNA-Seq reads form fragments (shown in green, blue, red and purple), which further generate clusters for NTR identification.D1: the distance between a putative NTR and any annotated transcribed regions; D2: the distance between two putative NTRs; B. Systematic workflow of NTR identification.The numbers for fragments are circled in blue, while numbers for clusters are in red circles.Singleton fragments are one-fragment clusters with length over 50 bp.doi:10.1371/journal.pone.0160197.g001 Cuffdiff module.S1 Fig shows the expression profiles of these 152 NTRs.Except for NTR20 and NTR10 that were highly expressed in 50% epiboly, all the other NTRs showed relatively low expression levels (FPKM 10) in the four studied early developmental stages.The average expression values were below 1 FPKM for more than 90% of the NTRs.Five NTRs (NTR8, NTR41, NTR61, NTR81 and NTR115) were upregulated at 50% epiboly by more than 1 FPKM compared with the 512-cell stage (S2A Fig).Several other NTRs also demonstrated a clear upregulated pattern after MBT albeit at lower expression levels (S2B Fig).Many NTRs showed a peak in expression levels at the 512-cell stage (S2D, S2E and S2F Fig).Out of the six selected expression patterns, two displayed relatively high expression values (S2A and S2D Fig).The NTRs displaying an increase after MBT (NTR61, NTR81 and NTR115) may be associated with processes important for organismal and anatomical structure development (S2A Fig) [7, 8].The NTRs showing relatively high maternal expression were expressed throughout MBT to subsequently diminish in 50% epiboly stage (S2D Fig).We randomly selected 10 NTRs for gene expression profile validation.According to our RNA-Seq data analysis, 8 out of the 10 NTRs were downregulated as the embryo development progressed from 1-cell to 50% epiboly (Fig 4A).Two exceptions were NTR36, with expression level showing 30-fold change from 1-cell to 50% epiboly, and NTR133, which was only detected as expressed at 50% epiboly.All 10 NTRs were validated using qRT-PCR (S3 Fig).The expression levels of NTR36, NTR57, NTR76, NTR88 and NTR151 were low at all studied stages (S4 Fig).The expression measured by qRT-PCR replicated a similar expression profile

Fig 2 .
Fig 2. Position of discovered NTRs in the zebrafish genome.The USCS Zebrafish Genome Graphs tool was used to generate the figure.The 152 NTRs with approximate genomic positions (vertical lines in black) were detected on each chromosome (represented by grey bars) of the zebrafish genome.The NTRs marked with crosses were validated using RT-PCR, and those marked with asterisks were evaluated using qRT-PCR.doi:10.1371/journal.pone.0160197.g002

Fig 4 .
Fig 4. Gene expression profiles of 10 randomly selected NTRs. A. The upper panel shows the gene expression fold change of the randomly selected 10 NTRs calculated from RNA-Seq data, except NTR 36 whose fold change was over 29 at 50% epiboly, and NTR133 which was not expressed at 1-cell stage; B. The lower panel shows the gene expression fold changes of the 10 randomly selected NTRs calculated from qRT-PCR.The 1-cell stage was used as basal line for fold change calculations.Developmental stages are given on the X-axis.MBT-Midblastula transition.doi:10.1371/journal.pone.0160197.g004 Fig), which might be a reason for them remaining previously unannotated.However, when validating using qRT-PCR, all 10 randomly selected NTRs were detected, albeit at different expression levels (Fig 4B, S4 Fig).In addition, the expression patterns shown for the different NTRs were similar when comparing the gene expression profiles between RNA-Seq and qRT-PCR (Fig 4A and 4B).Furthermore, the expression profiles show that the majority of the investigated NTRs peak in expression levels at MBT or at 50% epiboly (post-MBT).(S1 and S2 Figs) This indicates that processes during MBT induce or activate the expression of the NTRs, although the majority of the 152 NTRs are present already at 1 cell stage (S1 and S2 Figs).It has been shown in previous studies that there is an overall increase in expression post-MBT [7, 8].

Table 1 .
Read counts (in million) in four developmental stages.A sum of reads in million from all three runs at corresponding developmental stages.**Proportion against to the total reads for each developmental stage. *

Table 2 .
Putative NTRs without annotation in NCBI Zebrafish Annotation.

Table 2 .
(Continued) structures of NTRs (Fig3).One isoform of NTR50, NTR50_1, was validated as amplicon1 in Fig 3A.We were not able to detect the predicted isoform NTR50_2 in 50% epiboly.Sanger sequencing confirmed both isoforms of NTR88 (Fig 3B), in agreement with the predicted structures.The two isoforms were supported by ORF prediction as well.For NTR103, besides two predicted isoforms being confirmed, we also detected two additional isoforms, amplicon2 and amplicon3 (Fig3C).For NTR145, we found three additional alternative splicing patterns at the 5'-end of the NTR (Fig3D).

Table 2 .
(Continued)The start and end position were defined according to the coverage of RNA-Seq reads.
#ORF prediction was performed with GENSCAN program.§ The NTR was marked as "Yes" if there were reads covered in the promoter region according to zebrafish CAGE data.doi:10.1371/journal.pone.0160197.t002