High-throughput RNA sequencing reveals the effects of 2,2′,4,4′ -tetrabromodiphenyl ether on retina and bone development of zebrafish larvae

Background 2,2′,4,4′-Tetrabromodiphenyl ether (BDE47) is a prevalent environmental pollutant and has been demonstrated to be a serious toxicant in both humans and animals, but little is known about the molecular mechanism underlying its toxic effect on the early development of vertebrates. BDE47-treated zebrafish larvae were found to present the light-related locomotion reduction in our previous study, therefore, we aimed to use high throughput sequencing to investigate the possible reasons from a transcriptomic perspective. Results By exposing zebrafish embryos/larvae to 5 μg/l and 500 μg/l BDE47, we measured the influence of BDE47 on the mRNA expression profiles of zebrafish larvae until 6 days post-fertilization, using Illumina HiSeq 2000 sequencing. Differential expression analysis and gene enrichment analysis respectively revealed that a great number of genes, and gene sets based on two popular terminologies, were affected by the treatment of 500 μg/l BDE47. Among them, BDE47 caused changes in the retinal metabolism and related biological processes involving eye morphogenesis and visual perception, as confirmed by disordered photoreceptor arrangement and thickened bipolar cell layer of larval retina from histological observations. Other altered genes such as pth1a and collaborative cathepsin family exhibited disrupted bone development, which was consistent with the body curvature phenotype. The transcriptome of larvae was not significantly affected by the treatment of 5 μg/l BDE47, as well as the treatment of DMSO vehicle. Conclusions Our results suggest that high BDE47 concentrations disrupt the eye and bone development of zebrafish larvae based on both transcriptomic and morphological evidences. The abnormal visual perception may result in the alteration of dark adaption, which was probably responsible for the abnormal larval locomotion. Body curvature arose from enhanced bone resorption because of the intensive up-regulation of related genes. We also proposed the larval retina as a novel potential target tissue for BDE47 exposure. Electronic supplementary material The online version of this article (doi:10.1186/s12864-014-1194-5) contains supplementary material, which is available to authorized users.


Background
Polybrominated diphenyl ethers (PBDEs) are commonly used as flame retardants. Although the use of commercial mixtures such as pentabromodiphenyl ether and octabromodiphenyl ether in products has been banned or limited by the European Union, the U.S., and other countries, PBDE residues in the environment and in animal tissues still pose a serious threat [1,2]. Among PBDE congeners, 2,2′,4,4′-tetrabromodiphenyl ether (CAS: 5436-43-1) is the most frequently detected one [3], and previous studies have demonstrated it to be the most toxic [4,5]. Exposure to BDE47, an endocrine-disrupting chemical, can cause multiple adverse effects, and current research has focused on the examination of its endocrine-disrupting activity and neurobehavioral toxicity [3,6]. The best-known mechanism of BDE47 toxicity involves its ability to impair the homeostasis and function of thyroid hormone (TH) in animals because of its structural similarity to THs such as triiodothyronine (T3) and thyroxine (T4). Moreover, BDE47 can lead to various types of developmental neurotoxicity, including long-lasting behavioral alterations, by interfering with signal transduction pathways and by damaging neuron cell viability [7].
Recently, deep sequencing (next-generation RNA sequencing) has become a fast-growing high-throughput technology. Compared to microarrays, deep sequencing is independent of predefined probe sets, thereby permitting the detection of novel transcripts and alternative splicing [8]. However, thus far, while RNA sequencing has been increasingly used in toxicological studies of rodents and humans, it has only been employed for toxicological studies of aquatic species in a few studies [9,10]. In our previous study, 500 μg/L of BDE47 caused a significant reduction of larval locomotion specifically at the switch between light and dark periods; the effects seemed to be related to light stimulus [11]. The evidence for this was that a PBDE mixture, DE-71, containing BDE47 was found to cause biochemical changes in the larval eye as measured by optokinetic responses and phototaxis tests [12]. Therefore, to validate our hypothesis and understand the novel toxic mechanism of BDE47, we employed deep sequencing to obtain and analyze whole-transcriptome information for zebrafish larvae after exposure to BDE47.
In this study, the Illumina HiSeq sequencing platform was utilized to investigate the influence of PBDEs on the gene expression profiles of zebrafish larvae. We chose 6 day post-fertilization (dpf ) zebrafish larvae to keep aligned with our previous work [11], and zebrafish at this stage already have many well-developed organs, for example retina. First, this study discovered that in 6 dpf zebrafish, BDE47 exposure led to major adverse effects, including disrupting the visual perception and bone development of larvae. Furthermore, the reduction of zebrafish larval locomotion was probably caused by the disruption of retinol metabolism by light stimulus, and body curvature was probably caused by abnormal bone development induced by BDE47.

Characterization of the transcriptome of zebrafish larvae
The original RNA-Seq data are available at the Gene Expression Omnibus database (accession number GSE59968). Sequencing produced approximately 48 M total reads from a blank control sample, 56 M total reads from a sample of 0.1% DMSO vehicle, 42 M total reads from a sample of 5 μg/L BDE47, and 38 M total reads from sample of 500 μg/L BDE47. Of the total number of reads in these four groups, the mapping/unique mapping rates were 89.05%/86.67%, 89.42%/87.11%, 90.21%/ 88.06%, and 89.15%/86.87%, respectively. Although the sequencing depth in this study was lower than that in some previous experiments [13], high unique mapping rates provided the basis for adequate analysis. 21098 transcripts were detected in at least one of four groups.
Excluding 753 non-coding RNA and 521 miscellaneous RNA, here we used 19824 mRNA transcripts in further steps.
Because TMM normalization was considered to be superior to reads per kilobase of the transcript per million mapped (RPKM) normalization in Illumina-based RNA sequencing [14], TMM was adopted as our normalization strategy. The most abundant transcripts in the four groups, including si, myhz1, ef1a, LOC100535672, try, pvalb2, myhz2, krt4, hsp90ab1, and pvalb1, varied slightly in content across all treatments. However, the 500 μg/l BDE47-treated group had more distinct impacts on the abundance of transcripts than the other 3 groups (Additional file 1: Table S1).

Differential expression of genes with BDE47 treatments
The differentially expressed genes of zebrafish larvae treated with two BDE47 concentrations are shown as a Venn diagram ( Figure 1). By comparing BDE47 treatments with the DMSO solvent control, 2235 transcripts were affected after 500 μg/L exposure, of which 1338 were up-regulated and 897 were down-regulated, and 552 transcripts were affected after 5 μg/L exposure, of which 155 were up-regulated and 397 were downregulated. These values clearly indicate that exposure to a high BDE47 concentration resulted in a higher number of up-regulated genes. The intersections of the Venn diagram further revealed that less than 4% of upregulated genes in the 500 μg/L BDE47-treated group were also up-regulated in the 5 μg/L BDE47-treated group, and the expression trends of a total of 32 genes reversed with increasing BDE47 concentration. These features suggested that exposure to different BDE47 concentrations resulted in different types of biological effects in the zebrafish development process. The maximum FC of genes was 15.67 (zgc:194355) when exposed to 5 μg/L BDE47 and up to 115.04 (hsp70.2, LOC100535101) when exposed to 500 μg/L BDE47. Complete lists of differently expressed genes treated with both BDE47 concentrations are shown in Additional files 2 and 3: Tables S2 and S3.
The most significant genes (FC > 10 or FC < 0.1) of zebrafish larvae after exposure to 500 μg/L BDE47 are shown in Table 1. Excluding those genes that were totally inhibited in the exposure group and with a TMM value < 1 in the DMSO group, 49 up-regulated and 22 downregulated genes were present. Of the four most dramatically up-regulated genes, three transcripts encoded heat shock cognate (HSC) proteins: their mean FC was up to 100.59. Unlike canonical heat-shock proteins, the HSC proteins constitutively expressed and functioned during normal cellular processes; however, their upregulation still reflected immune stress such as inflammation, infection, and cancer [15]. In the range of 10 < FC < 50, the notable transcripts were those encoding cathepsin and PIM1 enzymes. The protein products of cathepsin-related transcripts were similar to those of human cathepsin V, a protease highly homologous to mouse cathepsin L [16]. The Pim1 oncogene is involved in multiple human cancers, and enhanced Pim1 expression would benefit vision formation of zebrafish [17].
Marginal effect of BDE47 on alternative splicing events One advantage of deep sequencing over microarrays is its ability to detect alternative splicing, which is a common event in higher vertebrates. Therefore, the effect of BDE47 on the alternative splicing of larval transcriptome was investigated. The results are shown in Table 2. Many diseases in animals are due to abnormal alternative splicing. However, BDE47 exhibited only a marginal effect on the alternative splicing events of zebrafish larvae despite seriously affecting its transcriptome. We noticed that a similar outcome was observed in the case of benzene exposure [18]. Two genes of tropomyosin, TPM3 and TPM4, were significantly changed their alternative splicing patterns in the 500 μg/L BDE47 exposure group. Tropomyosin, which has diverse isoforms, plays a critical role in regulating the function of actin filaments; in particular, TPM3 may have a unique and vital role in embryogenesis.
Functional enrichment based on different dynamic expression patterns All significant mRNAs were further classified into eight dynamic expression patterns under different concentrations, and the mRNA number and p-value of each pattern were calculated ( Figure 2). Three patterns were significant: pattern 4 (up-regulated; only exposed to 500 μg/l BDE47), pattern 3 (down-regulated; only exposed to 500 μg/l BDE47), and pattern 0 (constantly down-regulated; exposed to both 500 and 5 μg/l BDE47). The results reaffirmed that the exposure of zebrafish larvae to high BDE47 concentration exerted a significant influence on their transcriptome, and most genes remained stable expression under low BDE47 concentration.
According to the expression patterns of zebrafish larvae, functional enrichment analysis was performed, and biological process ontologies were adopted in GO analysis. Because of the excessive numbers of terms obtained, the false discovery rate (FDR) was used to calibrate the p-value. The complete GO analysis results (p < 0.001, FDR < 0.05) are shown in Figure 3b. Pattern 4 represented the biological functions and processes that were significantly induced by BDE47 treatment. The affected GO terms mainly involved stress responses and metabolic processes such as: response to stress (GO:0006950); small molecule metabolic process (GO:0044281); nerve development (GO:0021675); and steroid metabolic process (GO:0008202). Patterns 3 and 0 represented the biological functions and processes that were significantly inhibited by BDE47 treatment. In pattern 3, the GO analysis results reflected BDE47 damage to the development of nerve, vision, and skeletal system in zebrafish larvae, e.g., extracellular matrix organization (GO:0030198); visual perception (GO:0007601); collagen fibril organization (GO:0030199); and skeletal system development (GO:0001501). The most enriched terms in pattern 0 included: protein heterotrimerization (GO: 0070208); extracellular matrix organization (GO:0030198); complement activation, alternative pathway (GO:0006957); and blood coagulation extrinsic pathway (GO:0007598).
KEGG pathway analysis ( Figure 3c) screened out fewer terms than GO analysis, and reflected different information from that obtained by GO analysis. The significant terms in pattern 4 included: retinol metabolism (PATH:00830); chemical carcinogenesis (PATH:05204); steroid hormone biosynthesis (PATH:00140); metabolism of xenobiotics by cytochrome P450 (PATH:00980); and sulfur metabolism (PATH:00920). The affected pathway terms in patterns 3 and 0 included: protein digestion and absorption (PATH:04974); ECM-receptor interaction (PATH:04512); and glycosaminoglycan biosynthesis-chondroitin sulfate (PATH:00532). Patterns 0 and 3 had similar enrichment results. The gene regulatory network of BDE47 was constructed on the basis of the information provided by the KEGG pathway database and visualized ( Figure 3).

Validation of sequencing data
To validate the sequencing data of BDE47 treatment, the FCs of relative mRNA expression levels of several selected genes were determined by real-time quantitative Table 1 Most significant genes changed their expressions in zebrafish exposed to 500 μg/l BDE47 (FC > 10 or FC < 0.1) PCR (RT-qPCR). Because, in the sequencing data, most of the genes significantly changed their expressions under exposure to high BDE47 concentration but showed no significant changes of their expressions under low concentration, only a RT-qPCR assay using 500 μg/L BDE47 was performed. The genes are listed in Table 3, and they include ch25h, LOC799791, rx2, fetub, sesn2, opn1sw1, hsp70, cyp2k6, pth1a, and lepb. The selected genes had significant expression changes with enough dimensions and they participated in biological functions we were concerned with, including xenobiotic metabolism (cyp2k6), immune response (sesn2, hsp70, ch25h, and LOC799791), endocrine regulation (fetub, pth1a, and lepb), eye development (rx2), and photosensitivity (opn1sw1). After RT-qPCR amplification, they were confirmed to have expression patterns similar to those observed by deep sequencing. Even though the FC values of pth1a and lepb were lower than their sequencing data, the expression changes were still dramatically significant.

Transcriptomic effects of DMSO vehicle
To study the possible effects of DMSO on zebrafish larvae, we compared the transcriptomes of the DMSO vehicle and blank control. Four hundred and fifty-four transcripts were significantly expressed (p < 0.05), of which 214 transcripts were up-regulated and 240 transcripts were downregulated. Significantly changed transcripts accounted for about 2.2% of the total detected transcripts, and no abnormal phenotype was observed during larval development (data not shown). The effects of DMSO vehicle on transcriptome in our study were far weaker than those in a previous study of Turner et al. [19]. The most up-regulated gene was crfb9 (FC = 17.19), which encoded the cytokine receptor and its human homolog displayed inhibitory properties on IL-22 effects. The expressions of 31 genes were totally inhibited, such as cyp11a1, gtf3ab, and msh4. Functional enrichment analysis revealed that DMSO had a major effect on inducing the immune responses of larvae. Detailed results are shown in Additional file 4: Table S4. Morphological and histological abnormality caused by BDE47 exposure To confirm the relationship between locomotor changes and disrupted larvae vision, histological sections of larval retina were prepared to investigate the effects of BDE47 on eye histology ( Figure 4A-4B). In the sections exposed to 500 μg/L BDE47, a clear morphological alteration of photoreceptor cells was observed, and in the remaining sections, no significant change (data not shown) was observed compared to the control. At 6 dpf, disorderly arranged rods and cones (photoreceptors) were observed. A bipolar cell layer (BCL) was also distributed more sparsely than in the control, and the thickness of BCL increased noticeably. In general, this structural alteration of the retina closely resembled the effects caused by polychlorinated biphenyls [20], indicating that some common mechanisms probably exist among polyhalogenated aromatic hydrocarbons. Another apparent effect of 500 μg/l BDE47 exposure on zebrafish larvae was the body curvature phenotype. Almost all larvae exhibited different degrees of curvature (88.9 ± 2.9% of malformation rate at 6 dpf); however, it was particularly evident in some individuals ( Figure 4C-4D). A close examination revealed that the curvature was related to the failed inflation of the swim bladder, and we hypothesized that the curvature may result in the malformation of the swim bladder. The abnormality of the notochord and swim bladder could have implications for larval locomotion. Besides, there was no significant difference in hatching rates in both of the BDE 47 exposure groups (data not shown).

Discussion
According to the results of our study and previous ones, the locomotion of zebrafish was found to be less active after exposure to PBDE. Several plausible hypotheses were suggested, including disruption of the transport and absorption of the neurotransmitter acetylcholine [21], inhibition of the axonal growth of primary and secondary motor neurons [22], and deformation of the larvae [23,24]. Microscopic examination confirmed the spinal curvature of larvae under exposure to 500 μg/L BDE47. However, these above hypotheses do not explain the specific locomotor reduction during switching between light and dark conditions [11]. We believe that the most relevant cause is larval responses to light. Therefore, based on deep sequencing data, our interpretation focuses on vision. The expression levels of dozens of retina-related genes were changed by exposure to BDE47 ( Figure 5). In addition, GO or KEGG terms such as visual perception and retinol metabolism had a relatively high significance in our enrichment analysis results (Figure 2).
Considering the essential roles of the conversion between retinol and retinal in the visual cycle and eye development, the up-regulation of retinol metabolism suggested that BDE47 probably affected these processes. The configurations of retinal (retinaldehyde) determine their functions during the visual cycle: in the traditional rod visual cycle, 11-cis retinal is regarded as the only chromophore. However, some evidences have revealed that in the newly discovered cone visual cycle, the retina could also use all-trans retinal directly [25,26]. The expression changes of several retinol dehydrogenase (RDH)-encoding genes were examined ( Figure 5). The results revealed that expressions of 11-cis RDH genes such as rdh5, rdh10b, and LOC555864 were up-regulated, as well as an all-trans RDH rdh12l. Therefore, BDE47 probably increased the photosensitive ability of zebrafish larvae and shortened their dark adaptation procedure, resulting in lack of responsiveness to changed light conditions.
The expression changes of genes in terms of visual perception indicated the possibility of abnormal eye development in larvae, which was confirmed by histological sections ( Figure 5). For example, the down-regulated expressions of the cone-rod homeobox-encoding gene Crx (FC = 0.48) and the R-cadherin-encoding gene CDH4 (FC = 0.49) were found in BDE47 exposed zebarfish larvae. The deficit of the Crx function led to the failure of forming the photoreceptor outer segment (OS) and the inhibition of expressing OS-specific genes [27]. Cdh4 controls visual system development and functions of cell differentiation and axon migration of retina in cooperation with Crx [28]. Their simultaneous down-regulated expressions suggested the impaired development of the retina, especially of photoreceptor cells.
BDE47 was reported to cause abnormal dorsal curvature of the trunk and tail from 72 hpf [22], and this observation was reproduced in our experiments. We found some evidences from our sequencing results. With the exception of hsp70.2, pth1a had the most intensive expression change among all transcripts (FC = 112.61) in 500 μg/L BDE47 treated larvae. The gene encodes zebrafish parathyroid hormone (PTH), which regulates calcium and phosphorus concentrations, vitamin D metabolism, and bone cell activity [29]. Increased PTH could elevate the distribution of calcium in blood and restrain absorption of calcium into bones, as well as induce cathepsin proteins [30]. Besides those transcripts in Table 1, several transcripts encoding other cathepsin members were up-regulated (FC < 10) their expressions in the complete differential expression analysis results. Because of the key function of cathepsins in bone resorption, these changes indicated that the formation of larval spinal curvature was probably caused by the restricted absorption of calcium in the bone via the cooperation of PTH and cathepsins. Our functional enrichment analysis (Figure 3b) also reflected the adverse impacts of BDE47 on the larval skeletogeny process.
PBDEs are known for their ability to interrupt the production, transport, and metabolism of T3 and T4. However, our studies showed no evidence that expressions of TH genes were affected by BDE47, including genes encoding TH receptors, transthyretin, and T3 deiodinases. Different effects of BDE47 were seen in adults and larvae. During development, the zebrafish thyroid forms the first follicle from around 60 hpf [31], and all the observed TH changes occurred in zebrafish Figure 3 Gene regulatory network of 6 dpf zebrafish larvae under BDE47 treatments. The abbreviations on arrows between two nodes reflect their regulatory relationship. A, activation; b, binding/association; c, compound; dep, dephosphorylation; e, expression; ind, indirect effect; inh, inhibition; p, phosphorylation. Node colors indicate gene expression pattern. Red, pattern 0; violet, pattern 1; green, pattern 2; teal, pattern 3; tawny, pattern 4; yellow, pattern 5; lilac, pattern 6; lavender, pattern 7.
larvae older than that stage. These observations pose some interesting questions: when and why did BDE47 begin to disrupt zebrafish THs? What role did BDE47 play on the total thyroid system of the zebrafish larvae? The only affected hormone in our experiment was calcitonin, which is also secreted from the thyroid; the expression of its gene calca was moderately up-regulated (FC = 2.37). Calcitonin can both respond to PTH, and antagonize the biological effects of PTH [32]. Therefore, in 6 dpf larvae, the primary influence of BDE47 on the thyroid was also relevant to the disturbance of bone formation.

Conclusions
In general, exposure to BDE47 resulted in the body curvature phenotype and abnormal locomotion. Distinct from the reported explanation focused on neurological factors, our transcriptomic analysis revealed that the phenotypes resulted from the exposure of BDE47 were due to disrupt the vision and bone development in zebrafish larvae. Light was a critical factor regulating fish behavior. Hence, altered photosensitivity and locomotion were anticipated to have significant consequences on the prey, avoidance, and survival of zebrafish larvae. This study proposes the larval retina as a potential target organ of BDE47 exposure. Although the toxicological effects identified in this study resulted from high levels of exposure, it is likely that our findings will be helpful for understanding the underlying mechanisms and developmental consequences of the toxicity of BDE47 and other PBDE congeners in environmental situations featuring lower concentrations, but longer terms and more complex pollutant mixtures, and therefore merit further investigation.

Zebrafish
Wild-type Tuebingen zebrafish (Danio rerio) were placed in a recirculating filtration system using water treated by reverse osmosis (ESEN, China). Zebrafish were fed with live brine shrimp twice a day, and aquarium flake food (OSI, Hayward, CA) once a day during the spawning period. The system was maintained at 28.5°C with a 14-h light/10-h dark cycle. The main water quality parameters, such as pH, temperature, conductivity, and ammonia-N (NH 3 -N) were monitored weekly. Adult fish aged from 4 to 6 months were chosen for spawning. The embryos were collected within the first hour of lighting and rinsed with the water from filtration system. All related protocols were approved by the Animal Ethics Committee of Tongji University.
Toxicological exposure and morphological examination BDE47 (purity > 99%, Accustandard, New Haven, CT) was dissolved in dimethyl sulfoxide (DMSO) vehicle (purity > 98%, Amresco, Solon, OH). The final concentrations of BDE47 were 0, 5, and 500 μg/L of medium with 0.1% vehicle, respectively. A sample without a vehicle was also set as a negative control for zebrafish development. To maintain consistency with our previous behavioral assay [11], exposure in an incubator started from 3 to 4 hours post-fertilization (hpf) at 28.5°C. Exposure solutions were replaced 50% daily to minimize the changes of BDE47 concentrations. In each treatment, approximately 60 larvae exposed until 6 dpf were used for deep sequencing and RT-qPCR. Deep sequencing was performed once with the total RNAs isolated from the 6 dpf larvae exposed with BDE47 and DMSO started from 3 to 4 hpf whereas RT-qPCR assays were performed three times with the total RNA isolated from three different batches of the 6 dpf larvae exposed with BDE47 and DMSO, respectively. Approximately 30 embryos/larvae were used in triplicate for continuous morphological monitoring using an Olympus IX71 reversed microscope (Olympus, Japan) with a DP72 digital camera.

Total RNA isolation and sequencing experiments
The total RNA was extracted using TRIzol reagent (Invitrogen, Carlsbad, CA) and purified by the RNeasy Figure 5 Genes differentially transcribed related to vision formation and eye development in zebrafish larvae. Hierarchicalclustering-analysis-based transcription levels were performed on 25 related genes showing significant differential expression (p < 0.05) in larvae. Gene tree (left) and condition tree (top) were obtained using Pearson's uncentered distance metric calculated from all log10 transcription ratios (exposed/controls). Color scale from green to red indicate log10 ratios from −1.00 (10-fold down-regulation) to +1.00 (10-fold up-regulation). C: control; S: vehicle (0.1% DMSO); 5: 5 μg/l BDE47 treatment; 500: 500 μg/l BDE47 treatment.
Mini Kit (Qiagen, Germany) according to the manufacturer's instructions. The quality of the total RNA was examined using an Agilent 2200 TapeStation (Agilent Technologies, Santa Clara, CA) and then RNA library and sequencing was performed by the Illumina HiSeq 2000 platform (Illumina, San Diego, CA). In detail, the total RNA of zebrafish larvae was purified using Sera-Mag Oligo (dT) beads (Thermo Scientific, Indianapolis, IN) after DNase I treatment. The mRNA was fragmented by heating and then treating with sodium acetate. The cleaved RNA fragments were transcribed into firststrand cDNA using reverse transcriptase, followed by second-strand cDNA synthesis employing DNA polymerase and RNase H and cDNA purification using the QIAquick Gel Extraction kit (Qiagen, Valencia, CA). The double-stranded cDNA was further subjected to end repair employing T4 DNA polymerase, the Klenow fragment, and T4 polynucleotide kinase followed by a poly(A)-tailing procedure using Klenow exo-polymerase and ligation with an adapter using T4 DNA ligase. Adapter-ligated fragments were separated by agarose gel electrophoresis, and the desired range of cDNA fragments (200 ± 25 bp) was excised from the gel and enriched by PCR to construct the final cDNA library. After validation with Agilent 2200, the cDNA library was pair-end sequenced on a flow cell using Illumina HiSeq 2000.

Differential expression analysis and gene enrichment analysis
After obtaining sequences from the raw sequencing data, clean reads were mapped to the reference zebrafish genome (Zv9 assembly) and H. sapiens (hg19 assembly) using Basic Local Alignment Search Tool (BLAST) to increase the subsequent annotation efficiency. Mapped data were then normalized using the trimmed mean of M-values (TMM) normalization method. The statistical analysis of differential expression was performed by the DEGseq package in the R language [33], and "significant" was defined as significance level p < 0.05 and a fold change (FC) of 2. Gene functions were defined and annotated using gene ontology (GO) terminology (http://www.geneontology.org) and Kyoto encyclopedia of genes and genomes (KEGG) pathway (http://www.genome.jp/kegg/ pathway.html). To statistically identify differential gene sets between the control and treatments, the Fisher exact test was used for GO analysis (p < 0.001) and KEGG pathway analysis (p < 0.05). Differential alternative splicing events were calculated by employing the methodology of Zhou et al. [34], and the significance level was set as p < 0.05. The gene regulatory network was constructed on the basis of the information obtained from the KEGG pathway database and visualized by Cytoscape, which is a common open source software package.

Real-time quantitative PCR
After RNA quality control by gel electrophoresis and biophotometer (Eppendorf, Germany), the total RNA was reverse-transcribed into cDNA with Superscript III reverse transcriptase (Invitrogen), random hexamers, and oligo(dT) primers. SYBR Green I qPCR was performed using a 7500 Fast Real-Time PCR system (Applied Biosystems, Foster City, CA). The PCR amplification mixture contained 0.2 μL platinum Taq polymerase, 2 μL PCR buffer (10×), 0.5 μL forward primers, 0.5 μL reverse primers, 1 μL SYBR Green dye (20×), and 1 μL reverse transcription products to make a final volume of 20 μL. PCR conditions were 40 cycles of 95°C for 10 s, 60°C for 30 s, and 70°C for 30 s, and melting curves were generated after the last cycle. The threshold cycle values for selected genes and the housekeeping RPL13a gene were used to acquire the amount of RNA equivalents after the comparison with other alternatives like β-actin and 16S RNA. Primers for RT-qPCR were designed using the Primer Express software (Applied Biosystems, Foster City, CA); they are shown in Additional file 5: Table S5.
performed the histological sections of zebrafish larvae. JZ observed and recorded the morphologies of zebrafish larvae. TX, QZ, and BD drafted the initial manuscript and all members contributed to the preparation of the final manuscript. All authors read and approved the final version of this manuscript.