Transcriptome Analysis of the Zebrafish Model of Diamond-Blackfan Anemia from RPS19 Deficiency via p53-Dependent and -Independent Pathways

Diamond-Blackfan anemia (DBA) is a rare inherited bone marrow failure syndrome that is characterized by pure red-cell aplasia and associated physical deformities. It has been proven that defects of ribosomal proteins can lead to this disease and that RPS19 is the most frequently mutated gene in DBA patients. Previous studies suggest that p53-dependent genes and pathways play important roles in RPS19-deficient embryos. However, whether there are other vital factors linked to DBA has not been fully clarified. In this study, we compared the whole genome RNA-Seq data of zebrafish embryos injected with RPS19 morpholino (RPS19 MO), RPS19 and p53 morpholino simultaneously (RPS19+p53 MO) and control morpholino (control). We found that genes enriched in the functions of hematological systems, nervous system development and skeletal and muscular disorders had significant differential expression in RPS19 MO embryos compared with controls. Co-inhibition of p53 partially alleviates the abnormalities for RPS19-deficient embryos. However, the hematopoietic genes, which were down-regulated significantly in RPS19 MO embryos, were not completely recovered by the co-inhibition of p53. Furthermore, we identified the genome-wide p53-dependent and -independent genes and pathways. These results indicate that not only p53 family members but also other factors have important impacts on RPS19-deficient embryos. The detection of potential pathogenic genes and pathways provides us a new paradigm for future research on DBA, which is a systematic and complex hereditary disease.


Introduction
Diamond-Blackfan anemia (DBA) is a congenital anemia and broad developmental disease that develops soon after birth. The anemia results from a failure of erythropoiesis, with normal platelet and myeloid lineages, and it can be managed with steroids, blood transfusions, or stem cell transplantation [1]. DBA is the first known inherited disease that results from a defect in a structural ribosomal protein [2]. In humans, mutations in an increasing number of genes encoding RPs of the small (RPS19, RPS24, RPS17, RPS7, RPS10, RPS26) and large (RPL35A, RPL5, RPL11) ribosomal subunits have been shown to cause DBA. Moreover, RPS19 is the most frequently mutated ribosomal protein gene, which accounts for approximately 25% of DBA patients [3]. However, how the ribosomal protein mutations specifically affect hematopoiesis is still unclear.
The zebrafish (Danio rerio) is an excellent model organism for such studies because of its embryo transparency, high fecundity and fast development of organogenesis. It is notable that many of the key molecular players and events that drive organogenesis in zebrafish are evolutionarily and functionally conserved with other organisms, including mammalians [4]. Several groups have previously used zebrafish for modeling ribosomal protein knockdown in the embryo with morpholino technology or generating ribosomal protein mutants [5]. Recent studies have suggested that p53-mediated cell cycle arrest and/or apoptosis in erythroid cells could be the major factors in DBA development [3,5,6,7]. For example, regarding the RPL11 mutant in zebrafish, Danilova et al propose that the unique phenotype of DBA is the sum of several abnormally regulated molecular pathways mediated by the p53 protein family and p53-independent synergistic impacts on hematological and other cellular pathways affected in DBA [6]. In RPS7-deficient zebrafish embryos, p53 was activated, and its downstream target genes and biological events were induced, including apoptosis and cell cycle arrest. Furthermore, simultaneous knockdown of the p53 protein could partially reverse the abnormal phenotype of the morphants [3]. In the RPS29 mutant, genes that are up-regulated are enriched with genes that are upregulated by p53 after irradiation and the p53 knockdown almost completely rescues the RPS29 morphological and hematopoietic phenotypes, demonstrating that p53 mediates the effects of rps29 knockdown [5]. Knockdown of two other RP genes, RPS3A and RPL36A, result in severe morphological abnormalities with mild erythroid defects and elicited an activated p53 response. For the RPS19-deficient zebrafish, its phenotype is mediated by dysregulation of deltaNp63 and p53, and suppression of p53 and deltaNp63 alleviates the RPS19-deficient phenotypes [8]. At the same time, co-inhibition of p53 activity rescued the morphological abnormalities but did not alleviate erythroid aplasia in RPS19deficient zebrafish. Hence, both the p53-independent and p53dependent pathways could be responsible for the defective erythropoiesis in the zebrafish model of Diamond-Blackfan Anemia because of RPS19 deficiency [7].
In this study, we observed the hemoglobin synthesis defect at 48 hpf in RPS19 knockdown zebrafish embryos, and this defect can only be partially reversed by co-inhibition of p53 activity. To investigate the underlying regulatory mechanisms, we generated three types of zebrafish morphants by MO microinjection, including control morphants (control), RPS19 morphants (RPS19 MO), and RPS19 and p53 morphants (RPS19+p53 MO), and we performed transcriptome analysis of all of the pairs using the RNA-Seq technique. We found significant differentially expressed genes in RPS19 MO and RPS19+p53 MO compared with the controls. These genes are associated with the functions of cell cycle, hematological system and nervous system development and the skeletal and muscular disorders. Additionally, we determined the genome-wide p53-dependent and -independent genes and pathways. Our results demonstrated that members of the p53 network as well as other partners exert important impacts on RPS19-deficient embryos. The detection of potential pathogenic genes and pathways in this study will provide a new research paradigm for the study of DBA. Figure 1. Effectiveness of RPS19 morpholinos and hemoglobin staining of embryos coinjected with Rps19 mRNA and P53 MO using o-dianisidine. The Rps19:egfp construct was made by inserting a partial sequence of Rps19 cDNA (containing 60 bps from the 59 UTR ) with the Nterminus of egfp into modified pEGFP-N1 (the ATG codon of EGFP was removed). The sequence of RPS19 MO1 is a compliment of bp 1-24 of Rps19 cDNA. Embryos co-injected with 25 ng Rps19:egfp DNA and 5 ng control MO produced green fluorescent fusion protein (A), and expression of the fusion protein was inhibited by co-injection with 2 ng Rps19 Mo (B). O-staining results show a drastic reduction in the number of hemoglobin-stained blood cells when Rps19 is knockdown (C and D are the control, E and F are Rps19 knockdown) and partially reversed by co-injection of P53 morpholino (G and H). A, B, D, F and H are the lateral view; C, E and G are the ventral view. doi:10.1371/journal.pone.0071782.g001

Results
Zebrafish Phenotypes of RPS19-deficient and Coinhibition of p53 Activity using MO Developmental defects are found in approximately 40% of DBA patients with mutations in RPS19. Zebrafish RPS19 is approximately 87.7% identical to the human homolog. As previously demonstrated, the RPS19 deficiency in zebrafish results in hematopoietic and developmental abnormalities that resemble DBA. The embryos injected with control MO did not display any morphological changes. RPS19 morphants showed an obvious ventrally bent tail and a reduction in the circulating blood cells compared with the control embryos. At 48 hpf, the hemoglobin staining results showed that hemoglobin-stained blood cells in the heart region were markedly decreased in RPS MO, which is partially rescued in RPS19+p53 MO (Fig. 1). Further in situ results showed that the gata1 expression level are comparable in control Mo, RPS MO and RPS19+p53 MO (Fig. S1), which is consistent with previous results [8]. In addition, we used HSCs and definitive hematopoiesis markers cmyb and runx1 to examine the expression after RPS19 morpholino (2 ng) injection or Rps19 (mo) and P53 (mo) co-injection. We did not find significant changes of the expression of cmyb and runx1 (Fig. S1), indicating the divergent functions of RPS19 and RPL22 in hematopoiesis [9]. The effectiveness of translational inhibition by RPS19 MO was confirmed by examining green fluorescent fusion protein under fluorescence microscopy ( Fig. 1).

Transcriptome Profile Analysis
As previously demonstrated, the RPS19 deficiency in zebrafish results in hematopoietic and developmental abnormalities that resemble DBA. To determine the effect of RPS19 MO on the whole transcriptome and to delineate the function of p53 in RPS19-deficient embryos, we analyzed the transcript profiles using RNA-Seq. Three mRNA-Seq libraries were generated, including zebrafish embryos with control morpholino (control), RPS19 morpholino knockdown (RPS19 MO), and RPS19 and p53 morpholino knockdown simultaneously (RPS19+p53 MO). These libraries were sequenced using the Illumina Hi-seq 2000 Genome Analyzer platform with paired-end 100 base-pair tags to a depth of 35-60 million reads. We mapped the sequencing data to the latest zebrafish genome assembly version 2010 (Zv9). A total of 17-25 million reads could be mapped to the genome, representing 40%-46% of all of the generated reads (Table S1 in Tables S1). FPKM was used to estimate the transcript abundance. Overlapped and distinct gene expression among three samples was shown in Fig. 2. Unsupervised hierarchical clustering of genome-wide expression profiles (FPKM.1) was performed to examine the relationship among these three samples (Fig. 2). We observed that, compared with the control, the pattern of transcriptome profile of RPS19 MO is more similar to that of RPS19+p53 MO. The transcriptome difference between the profiles of RPS19 MO and the control is more significant than that of RPS19+p53 MO and the control. These results implied that simultaneous knockdown of p53 by co-injecting a p53 MO could partially reverse the abnormal phenotype in the RPS19 morphants.

Transcriptome Analysis of RPS19-deficient Embryos
To observe the transcriptome changes of RPS19-deficient embryos, we compared the differences of the whole genome RNA-Seq data between RPS19 MO and the control using the software Cufflinks. There were 47 up-regulated genes and 312 down-regulated genes in RPS19 MO compared with the control sample (fold-change.2.0, p-Value,0.05). These genes were subjected to Ingenuity Pathway Analysis (IPA) to identify the enrichment of genes in specific functional groups and pathways (IPA, http://www.ingenuity.com). The IPA accepts human UniGene IDs as one of the identifiers for data uploading and analysis. Hence, we mapped these differentially expressed genes to their human homologs using the HomoloGene database. Then, human homologs of up-and down-regulated genes of the RPS19deficient embryos were analyzed by using IPA tools in which the gene sets are enriched for a specific function/pathway and the enrichment is represented as a ratio. The main biological functions that were enriched by differentially expressed genes are associated with genetic disorder, neurological disease, cellular growth and proliferation, cancer and cell death (p-Value,0.05). The number of enriched genes for the above GO terms is more than 50, which accounts for at least 14% of the whole set of differentially expressed genes (Fig. 3). The significantly changed genes are associated with signaling pathways, including the coagulation system, acute phase response signaling, LXR/RXR activation, intrinsic prothrombin activation pathway, and the taurine and hypotaurine metabolism (Table S2 in Tables S1). The networks constructed by differentially expressed genes are associated with cellular development, cellular function and maintenance, development disorder, hematological system development and function, and neurological disease, each of which received an enrichment score of more than 20 (Fig. 4). Since the networks with functions of 'cellular development', 'cellular function and maintenance', 'Nervous system development and function', and 'cell death and survival' are inter-connected to form a big and complex regulatory network that is controlled by a few central factors and we observed expression changes in many genes and they could not be attributed to gene mutations (the probability of multi-gene mutations appear at the same time is very low), we assume this change to be a result of abnormal expression of a few central factors that affect the central nodes in the networks. In the network analysis, we found central nodes for which the human homologs are SST, OXA1L, CLDND1, MAL, Hsp70, GCG, PLC, SOX1, VAX1, POMC, and TSH.

Transcriptome Analysis of RPS19+p53 MO Embryos
To further study the role of p53 in RPS19-deficient embryos, we analyzed the RNA-Seq data of RPS19+p53 MO. Differential expression analysis is conducted by Cufflinks on the transcriptome data of RPS19+p53 MO and control/RPS19 MO. Compared with the control, RPS19+p53 MO is detected with 34 upregulated genes and 113 down-regulated genes (fold-change.2.0, p-Value,0.05). These genes were subjected to Ingenuity Pathways Analysis (IPA) to identify the enrichment of genes in specific functional categories and canonical pathways. We found that the main biological functions of these genes are associated with cellular growth and proliferation, hereditary disorder, cellular development, tissue development and small molecule biochemistry (p-Value,0.05) (Fig. 5). The top significant canonical pathways are LXR/RXR activation, acute phase response signaling and coagulation system (Table S3 in Tables S1). The main networks constructed by these differentially expressed genes are involved in small molecule biochemistry, cell cycle, cellular development, and hematological disease. The above networks all obtained enrichment scores of more than 20 (Fig. 6). In addition, the networks with functions of 'hematological disease', 'small molecule biochemistry', 'molecular transport', and 'cellular development' are found to be connected. We also found several central nodes in these networks, the human homologs of which are TGM2, A2M, LDL, HDL, Collagen, AGT, SERPINF1, NROB2, CETP, HBZ, hemoglobin, CLDND1, calpain, CD81, and G-protein beta.  However, by comparing the transcriptome data of RPS19+p53 MO and RPS19 MO, we found only 17 up-regulated genes and 11 down-regulated genes that were significantly expressed in RPS19+p53 MO embryos (fold-change.2.0, p-Value,0.05), which is much less than the number of significant differentially expressed genes of RPS19 versus the control. These results demonstrate that co-inhibition of p53 cannot rescue the abnormalities completely in RPS19-deficient embryos. The main enriched biological functions of these genes, by identified IPA, were associated with cancer, cell death, cellular growth and proliferation, cell signaling and small molecule biochemistry. Moreover, the canonical pathways are enriched on methane metabolism, DNA double-strand break repair by homologous recombination, stilbene, coumarine and lignin biosynthesis (data not shown). These functions might be corrected by p53 MO for RPS19-deficient embryos, which are significantly influenced by the expression of p53.

Changed Hematopoietic Genes in RPS19 MO and RPS19+p53 MO Embryos
To observe the defects of the hematopoietic functions in RPS19deficiency embryos, we detected the differentially expressed hematopoietic genes by searching the keyword 'hematopoietic' while using the software AmiGO in the Gene Ontology database (http://www.geneontology.org/). We observed 12 GO terms that are associated with 'hematopoietic' and 177 genes in the latest annotated zebrafish genome. We matched our data to these genes, which have the GO term keyword 'hematopoietic'. In the RPS19- deficient embryos, we found that no hematopoietic gene was upregulated significantly, but 5 hematopoietic genes were downregulated highly and significantly, such as bmper, fzd5, nrp1a, sema3d, and tbx1. The specific bio functions of these 5 genes are 'hematopoiesis'(bmper), 'T cell differentiation in thymus'(fzd5), and 'thymus development' (nrp1a, sema3d, tbx1). To determine whether co-inhibition of p53 would alleviate hematopoietic defect, we also matched the differential expressed genes in RPS19+p53 MO versus RPS19 MO to the 177 genes associated with hematopoietic in The Gene Ontology database. Interestingly, there was only one hematopoietic gene, fzd5, which was downregulated in RPS19 MO, but was significantly up-regulated in RPS19+p53 MO (Table S4 in Tables S1). Previous studies showed that in zebrafish, fzd5 was associated with not only eye and retina development, but also canonical Wnt signaling, T cell differentiation in thymus and early liver formation [10,11] and thus participated in hematopoiesis [12,13]. As the expression of fzd5 is reversed by co-inhibition of p53, hematopoiesis and other functions of fzd5 are reversed by p53 MO. In combination with our hemoglobin staining results of RPS19+p53 MO, we concluded that the hematopoietic defects could not be remedied completely through co-inhibition of p53 in RPS19-deficient embryos.

p53-dependent and -independent Genes and Pathways
Next, we observed the p53-dependent/independent pathways in RPS19-deficient embryos by comparing the transcriptome profiles of control, RPS19 MO and RPS19+p53 MO, using the software Cufflinks. To determine the p53-dependent genes, we firstly searched the up/down-regulated genes in RPS19 MO (compared with control). We then identified genes that their expression levels were reversed (down/up-regulated) in embryos co-injected with RPS19 MO+p53 MO, and defined them as p53 dependent genes. Here we presume that genes whose expression levels did not return to the normal level (control), but did have partial reverse effect compared with pathological state (RPS19 MO), are p53-dependent genes. To determine the p53-independent genes, we presume genes that are up/down-regulated in RPS19 MO (compared with control), and up/down-regulated in RPS19+p53 MO (compared with control), are p53-independent genes. No matter the degree of changes in these genes in RPS19+p53 MO reach the level of them in RPS19 MO or not, as long as the change trends are the same, we accept them as p53independent genes. (Fig. 7).
Based on the above strategy, we observed 12 potential p53dependent genes, in which 1 gene was up-regulated and 11 genes were down-regulated ( Table 1). The biological functions enriched by these genes are protein folding, integral to membrane and ATP binding. We also found 110 potential p53-independent genes, of which 21 genes were up-regulated and 89 genes were downregulated. The up-regulated genes are mainly associated with extracellular regions, structural molecule activity, response to stress and intracellular. The down-regulated genes are mainly associated with negative regulation of endopeptidase activity, negative regulation of peptidase activity, and negative regulation of hydrolase activity (p-Value,0.01, number of gene products.2) ( Table 2). There were 20 of the regulated pathways that were determined to be p53-dependent (Table S5 in Tables S1) and 76 that were determined to be p53-independent (Table S6 in Tables S1). Real-time PCR and in situ hybridization results are mostly consistent with the above sequencing results except for mt2. The real-time PCR results of mt2, fzd5, hsp70l, fn1b, optc, klf11a and srf1 are shown in Fig S2. Interestingly, the mt2 expression in sequencing results is different from the real-time PCR results. We postulated that the discrepancy may arise from: a) the design and optimization of gene-specific real-time primer pairs for mt2; b) the PCR amplification steps bear inherent biases; c) the associated tools for computational analysis are in their infancy, thus may contain errors. It is crucial to identify all existing isoforms in the long run for accurate quantification of a transcriptome. Notably, most of the genes we tested are consistent with the initial RNA-seq experiment data and thus we are confident about the quality of our data. The in situ hybridization data of fzd5, alas2, nrp1a and tagln are shown in Fig S3.

An Overview of Our Work
In this study, we provided a transcriptome profile of the zebrafish model of DBA, in which the DBA gene RPS19 is knocked-down by MO, for genome-wide analysis using the RNA-Seq technique. We determined differentially expressed genes with statistical significance and enriched pathways and networks by these genes among RPS19 MO, RPS19+p53 MO and control embryos, and we identified p53-dependent and p53-independent genes and pathways by comparing each pair of the three samples. Our data illustrate that not only p53 is the key factor contributing to the significant abnormalities of RPS19-deficient embryos but also some other important factors and pathways participate in regulating the abnormal phenotypes of RPS19-deficient embryos.

The Reliability of RNA-Seq Data
These samples were sequenced using the Illumina Hi-Seq 2000 Genome Analyzer platform with paired-end 100 base-pair tags to a depth of 35-60 million reads, which is sufficient sequence coverage for transcriptome profiling [14]. We mapped these reads to the zebrafish genome assembly version 2010 (Zv9) using TopHat. Approximately 17-25 million reads could be mapped to the genome, which represents 40%-46% of all of the generated reads. Thus, we believe that the whole genome RNA-Seq data of RPS19-deficient zebrafish transcriptome is of high quality and the genome-wide effects of p53 in RPS19-deficient embryos is reasonable and of scientific significance. It should be noted that the sorted erythrocytes could be a better source for RNA-seq sinc DBA is a red cell disease. However, we could not have enough

The Effects of RPS19 MO and RPS19+p53 MO on Zebrafish Embryos
Zebrafish response to RPS19 deficiency has many features resembling DBA, including impaired ribosome biogenesis, increased apoptosis, developmental abnormalities, and defective hematopoiesis [8]. With the results of the mapping, we finished the gene annotation and differential expression analysis using Cufflinks. For a comparison of the transcriptome of RPS19 MO and control embryos, we observed that the number of down-regulated genes is much more than that of the up-regulated genes in RPS19 MO. This result suggests that in RPS19-deficient embryos, more genes have been inhibited and cannot function properly. IPA provides us more information about gene ontology clusters, enriched pathways and networks of these differentially expressed genes. We found that these genes are mainly associated with cell death, developmental disorder, skeletal and muscular disorders, hematological system development and function, nervous system development and function, and tissue morphology. These results are in agreement with previous studies on ribosomal protein defects and the phenotype of RPS19-deficient zebrafish in this study [8,15]. Overall, our results indicate that, as a ubiquitous protein, ribosomal protein RPS19 deficiency has very strong effects on zebrafish, and many core factors and pathways are involved in the abnormalities.
The p53 pathway is a key pathway in ribosomal protein defect embryos, with additional roles in regulating other signaling molecules. p53 is always transcriptionally up-regulated when ribosomal protein suffers haplo-insufficiency [5,6,16,17]. However, it is still unclear why the mutations in the ubiquitously expressed RPS19 gene specifically affect erythropoiesis, and whether the transcriptome could return to its original state after p53 knockdown in RPS19-deficient embryos. To clarify the influence of the co-inhibition of p53 in RPS19-deficient embryos, we used RPS19 control Mo, Rps19 Mo, and Rps19 Mo+P53 Mo three groups for the comparative studies. We deduced the effects of p53 MO alone by comparing Rps19 Mo+P53 Mo with Rps19 Mo since we used the same dosage of Rps19 Mo. Our O-staining results showed that the phenotype of Rps19 Mo we obtained in the study can be rescued by the P53 and RPS19 double knockdown, confirming the effectiveness of the P53 Mo knockdown. This indicates that our current study and analysis method do not affect the robustness of the conclusions we made. We observed from the heatmap, which is the clustering of the transcriptome of the three samples, that there is a large difference among the profiles of these samples. Compared with RPS19 MO, the distance between RPS19+p53 MO and the control is smaller, suggesting that, after p53 knockdown, the abnormalities of RPS19-deficient embryos have not been retrieved completely but instead have been alleviated partially. Additionally, the number of differentially expressed genes for RPS19+p53 versus the control is much smaller, and the down-regulated genes are still the majority. IPA analysis demonstrated that the differentially expressed genes are mainly associated with functions of hematological system development and function, cell cycle, cellular development, cellular growth and proliferation, lipid metabolism, molecular transport, and small molecule biochemistry. These results suggest that these functions, regulated by some other important factors, could not rescued by the inhibition of p53 in RPS19-deficienct embryos.

Less Alleviation Effect of p53 on Hematological Defects
To further study the alleviation effects of p53 on hematological defects caused by RPS19 deficiency, we screened all the genes that are associated with hematopoiesis in the database (The Gene Ontology), and performed matching analysis. We found that, in the differentially expressed genes between RPS19 MO and RPS19+p53 MO, only one gene associated with hematopoiesis was up-regulated in RPS19+p53 MO (fold-change.2.0, pvalue,0.05). We supposed that this gene was alleviated by p53 MO.While in RPS19-deficient embryos there are more abnormally down-regulated genes, that are not rescued by p53 MO (Table S4 in Tables S1).This observation indicates that p53 knockdown could not completely rescue the developmental and functional defects of the hematological system in RPS19-deficiency zebrafish. This is expected since hematopoiesis is governed by multiple genes and some of them are not fully controlled by p53. Both our hemoglobin staining observation and RNA-Seq data analysis support the notion that the hematopoiesis is not fully functional in RPS19+p53 MO. p53-dependent and -independent Genes and Pathways In this study we developed a strategy to determine the p53dependent and -independent genes and pathways from the information on the transcriptome using RNA-Seq data of RPS19 MO, RPS19+p53 MO and the control. Our data shows that the number of p53-independent genes and pathways is significantly more than the number of p53-dependent genes and pathways. This result is in agreement with our hypothesis that many genes and pathways are outside the control of p53. The upregulated genes link to the functions of structural molecule activity, the extracellular region, response to stress, the intracellular region and more. At the same time, down-regulated genes are associated with the regulation of metabolic processes, the regulation of transcription, sequence-specific DNA binding transcription factor activity and more.

Conclusions
In conclusion, our study provides an outline of the transcript changes in the zebrafish model of DBA. We can conclude from this study that DBA is a systemic and complex disease due to ribosomal protein defects. DBA causes certain downstream phenotypes, including dysfunction of basic biological and physiological process such as transcription, translation, cellular metabolism, and pathways and networks related to many other disorders. The p53 network is one of the primary pathways in DBA disease. Meanwhile, there are many other p53-independent factors and pathways playing important roles in ribosomal protein defect. Our work has laid foundation in zebrafish DBA model and it will be interesting to see the clinical relevance of zebrafish DBA model in comparison with human DBA data.

Zebrafish Care
Breeding wild-type zebrafish (Danio rerio) (AB type) were maintained, and embryos were raised under standard library conditions [18]. Zebrafish embryos were kept in a 28.5uC incubator, and the stages (hours post-fertilization) in this study were as described [19]. All of the studies using zebrafish were approved by the Animal Care and Use Committee of Huazhong University of science and technology.

Microinjection of Morpholino RNA
The rps19 Morpholino (5-CACTGTTACACCACCTGG-CATCTTG-3), P53 Morpholino (5-GCGCCATTGCTTTG-CAAGAATTG-3), and RPS19 control MO (5-CACTcTTA-gACgCACCTGcCATgTTG-3) were obtained from Gene-Tools, LLC. Zebrafish embryos at the one-cell stage were injected with the MOs using a Microinjector (WPI SYS-PV830). Based on the literature and our initial injection trials, 2 ng MO and control Mo was chosen as the optimal concentration. Injected embryos were grown at 28.5uC and were observed under a microscope. The effectiveness of translation inhibition by RPS19 MO (2 ng/ embryo) was used by examining in vivo the rps19:egfp green fluorescent fusion protein under fluorescence microscopy. We found that at this dosage of 2 ng/embryo, most of the GFP expression was knocked-down.

Hemoglobin Staining
Hemoglobin in zebrafish embryos was analyzed using odianisidine (Sigma) as described [20]. Low-power images were collected using an Olympus microscope with a digital camera (OLYMPUS IX71). Images were imported into Adobe Photoshop CS2 9.0.2 for orientation and figure preparation. O-dianisidine staining results confirmed the effectiveness of the morpholino knockdown while the embryo morphology looks normal.

RNA Preparation
Immediately after harvesting, 40,50 pooled embryos at 48 hpf from different experiment replicates were snap-frozen in liquid nitrogen and stored at 280uC. Total RNA was extracted from the pooled embryos using TRIzol (Invitrogen) according to the manufacturer's instructions. RNA concentrations were determined using NanoDrop 2000 (Thermo Scientific). The integrity of RNA samples was determined using 1.2% Agarose gel electrophoresis, followed by removal of the residual genomic DNA with RNasefree DNaseI (Ambion).

RNA-Seq Data Analysis
Reads were processed and aligned to the UCSC zebrafish reference genome (build Zv9/danRer7, Jul. 2010) using TopHat (version 1.3.3) [21]. TopHat incorporates the Bowtie v0.12.7 algorithm to perform the alignment. Briefly, TopHat initially removes a portion of reads based on quality information accompanying each read, and maps the qualified reads to the reference genome. The reference index was built using Bowtie with a fasta file for the whole genome of zebrafish downloaded from UCSC (http://genome.ucsc.edu/). The parameters were set by default, but the number of threads to align the reads was set to 6. The aligned read files processed by TopHat were used by Cufflinks (version 1.2.1) software for further analysis, including assembling transcripts, estimating their abundances, and testing for differential expression and regulation in RNA-Seq samples [22]. To calculate gene expression intensity, the read counts were normalized to fragments per kilobase of transcript per million mapped reads (FPKM) according to the gene length and the total mapped reads [23]. Confidence intervals for FPKM estimates were calculated using a Bayesian inference method [24]. Once all of the short read sequences were assembled with Cufflinks, the output.GTF files were sent to Cuffcompare along with a reference.GTF annotation file downloaded from UCSC. This process classified each transcript as known or novel. Cuffcompare produces a combined.GTF file that is passed to Cuffdiff along with the original alignment (.BAM) files produced by TopHat. Cuffdiff then re-estimates the abundance of transcripts listed in the.GTF file using alignments from the.BAM file, and concurrently tests for differential expression. The expression testing is performed at the level of transcripts, primary transcripts and genes [25].

Cluster Analysis
Average linkage hierarchical clustering of gene expression intensity was performed in this study. The distance (d XY ) between genes and samples was measured by the Pearson correlation coefficient (r XY ), which is d XY~1 {r XY . The Pearson correlation coefficient is calculated through r~1 n{1 and n is the count of samples, X i and Y i are observed values of the two variables, and X X and Y Y are mean values of the two variables. Computation and visualization were achieved using the heatmap plus package in R. Genes with FPKM greater than 1 in at least one of the 6 samples were chosen to perform cluster analysis.

IPA Analysis
Gene interaction networks and signal pathways were generated using Ingenuity Pathway Analysis software (http://www.ingenuity. com/). IPA is biological data analysis software from IngenuityH Systems; it analyzes data from a variety of experimental platforms and provides accurate biological insight into the interactions between genes, proteins, chemicals, pathways, cellular phenotypes, and disease processes.
The differentially expressed genes processed by Cuffdiff and screened with our criteria (fold change.2, p-Value,0.05, we set the p-value to 0.05, which means that the error of discovery of our differential expression genes is 0.05 for statistical significance. It is a meaningful measurement similar to FDR and routinely used in such studies.) Genes identified to be significant were then submitted to IPA for biological function, canonical pathway, and interaction network analysis. Because the current edition of IPA cannot accept and identify the gene ID of zebrafish yet, we must match these differentially expressed genes of zebrafish to their homologs of human genes. The homolog genes database that we used was HomoloGene, which was downloaded from NCBI. After converting to human homolog genes, we submitted these homologs to IPA. IPA identified a set of genes that are enriched for a specific function or pathway, and the enrichment is represented as a ratio. To find the significance of the enrichment in a specific function, IPA calculates the significance value based on the measure of involvement of the gene in the input data set to their respective molecular function/signaling pathways [26]. The significance of the networks is calculated using Fisher's exact test, and the p-value is the executed negative logarithmic transformation.
Quantitative PCR (qPCR) 2 mg total RNA was used for reverse transcription with TransScript First-strand cDNA Synthesis SuperMix (TransGen). cDNA was then diluted 1:5 with RNase free H 2 O and 2 mL of this cDNA was used in the qPCR reaction. Quantitative PCR (qPCR) was performed using GoTaqH qPCR Master Mix (Promega) on an Applied Biosystems 7900 Real-Time PCR System (ABI 7900). The data were normalized against zebrafish b-actin. The sequences of the primers used for qPCR are given in Table S7 in Tables S1.

Whole Mount in situ Hybridization
For in situ hybridization, the following genes were used as probes: nrp1a, alas2, tagln and fzd5. cDNAs were amplified by RT-PCR and the products were cloned into pTA2 vector (TOYOBO). The sequences of the primers used for probe are given in Table S8 in Tables S1. Antisense probes were synthesized using T7 RNA polymerases (Promega) and labeled with DIG RNA labeling Mix (Roche). Embryos at 48 hpf were fixed in 4% PFA/PBS overnight at 4uC, dehydrated in 100% methanol and stored at 220uC before use. Rehydrated embryos were prehybridized for 3 h at 65uC followed by overnight hybridization with DIG-labeled probes. After washing, embryos were incubated with 2% blocking solution, followed by incubation with anti-Digoxlgenin-Ap (Roche), 1:4000 in blocking solution. After extensive washing, samples are visualized by incubating with NBT/BCIP AP stock solution (Roche) in AP buffer. Photographs were taken by using a stereo fluorescence microscope (OLYMPUS SZX2-ILLB).

Data Access
The raw sequence data can be accessed from the Gene Expression Omnibus, accession number is GSE45699 (available upon request or online when paper get accepted), and all of the genomic resources generated in this study are provided as Supplementary Files. Tables S1 File includes Tables S1-S8. (DOC)