Skip to main content
  • Research article
  • Open access
  • Published:

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



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.


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.


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.


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 down-regulated. 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 up-regulated 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.

Figure 1
figure 1

Venn diagram showing the effects of BDE47 exposure on the detected genes. Genes were grouped on the basis of their expression changes under different BDE47 concentrations. Expressions were evaluated on the basis of the number of reads per gene model. The shaded areas represent the existence of no intersection between two subsets.

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 down-regulated 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 up-regulation 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].

Table 1 Most significant genes changed their expressions in zebrafish exposed to 500 μg/l BDE47 (FC > 10 or FC < 0.1)

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.

Table 2 Genes with significantly altered alternative splicing under 500 g/l BDE47 exposure

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.

Figure 2
figure 2

Functional enrichment analysis based on different gene expression patterns. (A) Expression profiles in color indicate significant ones (p < 0.05). Green, up-regulated; red, down-regulated. Profile number (up left), gene number (bottom left), and trend (line) in each profile were labelled. The significant terms were graphed using (B) GO annotation (p < 0.001, FDR < 0.05), and (C) KEGG pathway annotation (p < 0.05). Tawny block, pattern 4; teal block, pattern 3; red block, pattern 0.

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).

Figure 3
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.

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 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.

Table 3 Comparison of gene expression changes between deep sequencing and qRT-PCR data

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 down-regulated. 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.

Figure 4
figure 4

Morphological observation of 6 dpf zebrafish larvae under BDE47 treatment and control. (A) Normal histological patterns of larval retina (40X). (B) Larvae exposed to 500 μg/l BDE47 had retinal morphology distinct from the control (40X). (C) Normal body type of zebrafish larvae in the control group (4X). (D) Body curvature phenotype after 500 μg/l BDE47 exposure (4X). RPE, retinal pigment epithelium; PCL, photoreceptor cell layer; INL, inner nuclear layer; IPL, inner plexiform layer; GCL, ganglion cell layer.

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).


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).

Figure 5
figure 5

Genes differentially transcribed related to vision formation and eye development in zebrafish larvae. Hierarchical-clustering-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.

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 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.


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.



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 (NH3-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 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 first-strand 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 ( and Kyoto encyclopedia of genes and genomes (KEGG) pathway ( 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.

Histological observation

First, fresh zebrafish larvae samples were fixed overnight in formalin/acetic acid/70%alcohol (18:1:1) at 4°C. Next, they were dehydrated using an alcohol dilution series (70%, 80%, 90%, 95%, followed by 100%). Finally, the samples were embedded into 100% paraffin. Serial sections were cut to a thickness of 3 μm and stained with hematoxylin-eosin staining. The slides were examined using a Motic BA310 microscope (Motic, Germany) with a Moticam Pro digital camera.

Availability of supporting data

The RNA sequencing data have been deposited in Gene Expression Omnibus with accession number GSE59968, which are available at


  1. Chernyak SM, Rice CP, Quintal RT, Begnoche LJ, Hickey JP, Vinyard BT. Time trends (1983–1999) for organochlorines and polybrominated diphenyl ethers in rainbow smelt (Osmerus mordax) from Lakes Michigan, Huron, and Superior, USA. Environ Toxicol Chem. 2005;24:1632–41.

    Article  PubMed  CAS  Google Scholar 

  2. Carlsson P, Herzke D, Wedborg M, Gabrielsen GW. Environmental pollutants in the Swedish marine ecosystem, with special emphasis on polybrominated diphenyl ethers (PBDE). Chemosphere. 2011;82:1286–92.

    Article  PubMed  CAS  Google Scholar 

  3. Birnbaum LS, Staskal DF. Brominated flame retardants: cause for concern? Environ Health Persp. 2004;112:9–17.

    Article  CAS  Google Scholar 

  4. Staskal DF, Diliberto JJ, DeVito MJ, Birnbaum LS. Toxicokinetics of BDE 47 in female mice: Effect of dose, route of exposure, and time. Toxicol Sci. 2005;83:215–23.

    Article  PubMed  CAS  Google Scholar 

  5. Darnerud PO. Toxic effects of brominated flame retardants in man and in wildlife. Environ Int. 2003;29:841–53.

    Article  PubMed  CAS  Google Scholar 

  6. Richardson VM, Staskal DF, Ross DG, Diliberto JJ, DeVito MJ, Birnbaum LS. Possible mechanisms of thyroid hormone disruption in mice by BDE 47, a major polybrominated diphenyl ether congener. Toxicol Appl Pharmacol. 2008;226:244–50.

    Article  PubMed  CAS  Google Scholar 

  7. Costa LG, Giordano G. Developmental neurotoxicity of polybrominated diphenyl ether (PBDE) flame retardants. Neurotoxicology. 2007;28:1047–67.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

  8. Scholz S. Zebrafish embryos as an alternative model for screening of drug-induced organ toxicity. Arch Toxicol. 2013;87:767–9.

    Article  PubMed  CAS  Google Scholar 

  9. Huang QS, Dong SJ, Fang C, Wu XL, Ye T, Lin Y. Deep sequencing-based transcriptome profiling analysis of Oryzias melastigma exposed to PFOS. Aquat Toxicol. 2012;120:54–8.

    Article  PubMed  CAS  Google Scholar 

  10. Olsvik PA, Berg V, Lyche JL. Transcriptional profiling in burbot (Lota lota) from Lake Mjosa-A Norwegian Lake contaminated by several organic pollutants. Ecotox Environ Safe. 2013;92:94–103.

    Article  CAS  Google Scholar 

  11. Zhao J, Xu T, Yin DQ. Locomotor activity changes on zebrafish larvae with different 2,2′,4,4′-tetrabromodiphenyl ether (PBDE-47) embryonic exposure modes. Chemosphere. 2014;94:53–61.

    Article  PubMed  CAS  Google Scholar 

  12. Chen LG, Huang YB, Huang CJ, Hu B, Hu CY, Zhou BS. Acute exposure to DE-71 causes alterations in visual behavior in zebrafish larvae. Environ Toxicol Chem. 2013;32:1370–5.

    Article  PubMed  CAS  Google Scholar 

  13. Vesterlund L, Jiao H, Unneberg P, Hovatta O, Kere J. The zebrafish transcriptome during early development. BMC Dev Biol. 2011;11: 10.1186/1471-213X-11-30

  14. Dillies MA, Rau A, Aubert J, Hennequet-Antier C, Jeanmougin M, Servant N, et al. A comprehensive evaluation of normalization methods for Illumina high-throughput RNA sequencing data analysis. Brief Bioinform. 2013;14:671–83.

    Article  PubMed  CAS  Google Scholar 

  15. Stricher F, Macri C, Ruff M, Muller S. HSPA8/HSC70 chaperone protein Structure, function, and chemical targeting. Autophagy. 2013;9:1937–54.

    Article  PubMed  CAS  Google Scholar 

  16. Bromme D, Li ZQ, Barnes M, Mehler E. Human cathepsin V functional expression, tissue distribution, electrostatic surface potential, enzymatic characterization, and chromosomal localization. Biochemistry. 1999;38:2377–85.

    Article  PubMed  CAS  Google Scholar 

  17. Yin J, Shine L, Raycroft F, Deeti S, Reynolds A, Ackerman KM, et al. Inhibition of the Pim1 Oncogene Results in Diminished Visual Function. PLoS One 2013;7. 10.1371/journal.pone.0052177

  18. Thomas R, McHale CM, Lan Q, Hubbard AE, Zhang LP, Vermeulen R, et al. Global gene expression response of a population exposed to benzene: A pilot study exploring the use of RNA-sequencing technology. Environ Mol Mutagen. 2013;54:566–73.

    Article  PubMed  CAS  Google Scholar 

  19. Turner C, Sawle A, Fenske M, Cossinsy A. Implications of the solvent vehicles dimethylformamide and dimethylsulfoxide for establishing transcriptomic endpoints in the zebrafish embryo toxicity test. Environ Toxicol Chem. 2012;31:593–604.

    Article  PubMed  CAS  Google Scholar 

  20. Wang YP, Hong Q, Qin DN, Kou CZ, Zhang CM, Guo M, et al. Effects of embryonic exposure to polychlorinated biphenyls on zebrafish (Danio rerio) retinal development. J Appl Toxicol. 2012;32:186–93.

    Article  PubMed  CAS  Google Scholar 

  21. Chen LG, Huang CJ, Hu CY, Yu K, Yang LH, Zhou BS. Acute exposure to DE-71: Effects on locomotor behavior and developmental neurotoxicity in zebrafish larvae. Environ Toxicol Chem. 2012;31:2338–44.

    Article  PubMed  CAS  Google Scholar 

  22. Chen XJ, Huang CJ, Wang XC, Chen JF, Bai CL, Chen YH, et al. BDE-47 disrupts axonal growth and motor behavior in developing zebrafish. Aquat Toxicol. 2012;120–121:35–44.

    Article  PubMed  CAS  Google Scholar 

  23. Padilla S, Hunter DL, Padnos B, Frady S, MacPhail RC. Assessing locomotor activity in larval zebrafish: Influence of extrinsic and intrinsic variables. Neurotoxicol Teratol. 2011;33:624–30.

    Article  PubMed  CAS  Google Scholar 

  24. Lema SC, Schultz IR, Scholz NL, Incardona JP, Swanson P. Neural defects and cardiac arrhythmia in fish larvae following embryonic exposure to 2,2′,4,4′-tetrabromodiphenyl ether (PBDE 47). Aquat Toxicol. 2007;82:296–307.

    Article  PubMed  CAS  Google Scholar 

  25. Ala-Laurila P, Cornwall MC, Crouch RK, Kono M. The action of 11-cis-retinol on cone opsins and intact cone photoreceptors. J Biol Chem. 2009;284:10422–32.

    Article  Google Scholar 

  26. Saari JC. Vitamin A Metabolism in Rod and Cone Visual Cycles. Annu Rev Nutr. 2012;32:125–46.

    Article  PubMed  CAS  Google Scholar 

  27. Furukawa T, Morrow EM, Cepko CL. Crx, a novel otx-like homeobox gene, shows photoreceptor-specific expression and regulates photoreceptor differentiation. Cell. 1997;91:531–41.

    Article  PubMed  CAS  Google Scholar 

  28. Babb SG, Kotradi SM, Shah B, Chiappini-Williamson C, Bell LN, Schmeiser G, et al. Zebrafish R-cadherin (Cdh4) controls visual system development and differentiation. Dev Dynam. 2005;233:930–45.

    Article  CAS  Google Scholar 

  29. Murray TM, Rao LG, Divieti P, Bringhurst FR. Parathyroid hormone secretion and action: Evidence for discrete receptors for the carboxyl-terminal region and related biological actions of carboxyl-terminal ligands. Endocr Rev. 2005;26:78–113.

    Article  PubMed  CAS  Google Scholar 

  30. Delaisse JM, Eeckhout Y, Vaes G. In vivo and in vitro evidence for the involvement of cysteine proteinases in bone resorption. Biochem Bioph Res Co. 1984;125:441–7.

    Article  CAS  Google Scholar 

  31. Alt B, Reibe S, Feitosa NM, Elsalini OA, Wendl T, Rohr KB. Analysis of origin and growth of the thyroid gland in zebrafish. Dev Dynam. 2006;235:1872–83.

    Article  Google Scholar 

  32. Holtrop ME, Raisz LG, Simmons HA. The effects of parathyroid hormone, colchicine, and calcitonin on the ultrastructure and the activity of osteoclasts in organ culture. J Cell Biol. 1974;60:346–55.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

  33. Wang LK, Feng ZX, Wang X, Wang XW, Zhang XG. DEGseq: an R package for identifying differentially expressed genes from RNA-seq data. Bioformatics. 2010;26:136–8.

    Article  CAS  Google Scholar 

  34. Zhou XX, Wu WW, Li H, Cheng YM, Wei N, Zong J, et al. Transcriptome analysis of alternative splicing events regulated by SRSF10 reveals position-dependent splicing modulation. Nucleic Acids Res. 2014;42:4019–30.

    Article  PubMed Central  PubMed  CAS  Google Scholar 

Download references


The authors thank Novel Bioinformatics Corporation, Shanghai, for excellent technical assistance. This study was supported by the National Natural Science Foundation of China [21477086] and the Collaborative Innovation Center for Regional Environmental Quality. It was also in part funded through support from the Swedish Research Council to the project [D0691301].

Author information

Authors and Affiliations


Corresponding author

Correspondence to Daqiang Yin.

Additional information

Competing interests

The authors declare that they have no competing interests.

Authors’ contributions

TX, JZ, and DY conceived the study and performed chemical exposure, RNA sequencing, and RT-qPCR experiments. TX and ZJ analyzed the data. TX 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.

Additional files

Additional file 1: Table S1.

The most abundant transcripts in four groups.

Additional file 2: Table S2.

The differentially expressed genes of zebrafish larvae treated by 500 μg/L BDE47 (p < 0.05).

Additional file 3: Table S3.

The differentially expressed genes of zebrafish larvae treated by 5 μg/L BDE47 (p < 0.05).

Additional file 4: Table S4.

The differentially expressed genes of zebrafish larvae between DMSO solvent control and negative control (p < 0.05).

Additional file 5: Table S5.

Primers used for quantitative PCR amplification.

Rights and permissions

Open Access  This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made.

The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder.

To view a copy of this licence, visit

The Creative Commons Public Domain Dedication waiver ( applies to the data made available in this article, unless otherwise stated in a credit line to the data.

Reprints and permissions

About this article

Check for updates. Verify currency and authenticity via CrossMark

Cite this article

Xu, T., Zhao, J., Yin, D. et al. High-throughput RNA sequencing reveals the effects of 2,2′,4,4′ -tetrabromodiphenyl ether on retina and bone development of zebrafish larvae. BMC Genomics 16, 23 (2015).

Download citation

  • Received:

  • Accepted:

  • Published:

  • DOI: