Transcriptomic Changes in Zebrafish Embryos and Larvae Following Benzo[a]pyrene Exposure

Abstract

Benzo[a]pyrene (BaP) is an environmentally relevant carcinogenic and endocrine disrupting compound that causes immediate, long-term, and multigenerational health deficits in mammals and fish. Previously, we found that BaP alters DNA methylation patterns in developing zebrafish, which may affect gene expression. Herein, we performed a genome-wide transcriptional analysis and discovered differential gene expression and splicing in developing zebrafish. Adult zebrafish were exposed to control or 42.0 ± 1.9 µg/l BaP for 7 days. Eggs were collected and raised in control conditions or continuously exposed to BaP until 3.3 and 96 h post-fertilization (hpf). RNA sequencing (RNA-Seq) was conducted on zebrafish embryos and larvae. Data were analyzed to identify differentially expressed (DE) genes (changed at the gene or transcript variant level) and genes with differential exon usage (DEU; changed at the exon level). At 3.3 hpf, BaP exposure resulted in 8 DE genes and 51 DEU genes. At 96 hpf, BaP exposure altered expression in 1153 DE genes and 159 DEU genes. Functional ontology analysis by Ingenuity Pathway Analysis revealed that many disease pathways, including organismal death, growth failure, abnormal morphology of embryonic tissue, congenital heart disease, and adverse neuritogenesis, were significantly enriched for the DE and DEU genes, providing novel insights on the mechanisms of action of BaP-induced developmental toxicities. Collectively, we discovered substantial transcriptomic changes at the gene, transcript variant, and exon levels in developing zebrafish after early life BaP waterborne exposure, and these changes may lead to long-term adverse physiological consequences.

Keywords: RNA-Seq; alternative splicing; benzo[a]pyrene; development; gene expression; zebrafish.

FIG 1.

BaP exposure strategy. Parental zebrafish were acclimated in new culture conditions for 7 days before BaP (42.0 ± 1.9 µg/l) or control exposure for 11 days. Fertilized eggs were collected during exposure days 7–11. Eggs were continuously exposed to BaP until 3.3 or 96 hpf before being collected for gene expression analysis. Arrows indicate collection of offspring for exposure or RNA sequencing (RNA-Seq).

FIG. 2.

BaP induced DE of genes and transcripts in zebrafish larvae. Smear plots show relative expression of genes (A) or transcripts (B) in BaP-treated zebrafish larvae compared with control at 96 hpf. X-axis is the log2 value of read counts per million (CPM). Y-axis is log2 fold change (FC). The blue horizontal lines indicate 2- or −2-fold. Black dots represent non-significant genes/transcripts, whereas red dots indicate significant DE genes/transcripts (P < .05). (C) PCA was performed on the RNA-Seq data from 96 hpf zebrafish samples. The percent variability attributed to the first 2 principal components is displayed on the X and Y-axes. Full color version available online

FIG. 3.

Predicted mechanisms by which BaP leads to activation of organismal death. IPA on gene expression results from 96 hpf zebrafish predicted that early life BaP exposure activated organismal death via transforming growth factor beta (tgf beta; A), bone morphogenetic protein 2 (bmp2; B), and growth differentiation factor 2 (gdf2; C). Red or green indicates genes were up- or down-regulated by BaP, respectively. Orange or blue indicates pathways or interactions that were predicted to be activated or deactivated, respectively, by BaP. Yellow lines indicate where our data were inconsistent with IPA predictions. Full color version available online

FIG. 4.

Predicted mechanisms by which BaP leads to inhibited neuritogenesis, reduced body size, and cardiovascular defects. IPA on gene expression results from 96 hpf zebrafish larvae predicted that early life BaP exposure inhibited neuritogenesis via transforming growth factor beta-1 (tgfb1) (A) and reduced the size of body via apolipoprotein E (apoe) (B). BaP was predicted to activate the pathways of atrial septal defect (C) and bleeding (D). Full color version available online

FIG. 5.

BaP effects on the expression of BaP marker genes and genes involved in the AHR signaling pathway. IPA revealed that BaP exposure in zebrafish during early development altered the expression of known marker genes of BaP exposure, including cyp1a1, cyp1a2, cyp1b1, and nqo1. AHR was predicted to be activated by BaP and influence the expression of AHR-responsive genes, including cyp1a1, cyp1a2, cyp1b1, cyp3a4, gsta1, gsta5, and nqo1. Full color version available online

FIG. 6.

Upstream regulators that were predicted to be activated by BaP. IPA predicted that BaP activated hepatocyte nuclear factor 1-alpha (hnf1a) (A), nuclear factor, erythroid 2-like 2 (nfe2l2; B), peroxisome proliferator-activated receptor gamma, coactivator 1 alpha (ppargc1a) (C), and PXR ligand-PXR-Retinoic acid RXRα complex (D) in 96 hpf zebrafish larvae, which may be causing the observed gene expression changes. Full color version available online.

FIG. 7.

DEU in developmental genes. DEXseq analysis on RNA-Seq data revealed that early life BaP exposure significantly altered exon usage of llgl2 (A), myhb (B), Tnc (C), and skiv2l2 (D) in 96 hpf zebrafish larvae. The y-axis represents normalized read counts of exons (exon usage), and the x-axis shows individual exons within a gene. The bars below the x-axis represent exons, and the lines between the bars represent introns. The numbers at the bottom are genomic locations of the gene. Exons highlighted in purple had a FDR < 0.05 in the DEU analysis of BaP and control samples. Red arrows indicate the significant DEU exons that had a FDR < 0.05 and an absolute fold change > 1.5. The fold change is the individual exon expression in BaP samples verses the control expression. Full color version available online.

FIG. 8.

IPA identified pathways affected by abnormal alternative splicing. IPA analysis on genes with DEU revealed that the AHR signaling (A) and development of body axis (B) pathways were affected by BaP exposure in 96 hpf zebrafish larvae.

FIG. 9.

Overlap of DEU and differential expressed (DE) genes. Venn-diagrams show the overlap of genes affected by expression and/or splicing between BaP and control zebrafish larvae (96 hpf; A). IPA was used to identify the canonical pathways (B), top functions (C), and top disease pathways (D) that are affected by BaP at splicing and expression levels. In panel B, bars represent the percentage of overlaping genes (in black) within a pathway. Dots indicate the negative log10 of the P-values. Larger −log(P-value) means that the pathway is more significant. The threshold for significance is marked in the graph as a dotted-line at 1.3 (−log(0.05)). In panels C and D, bars represent −log(P-value), and dotted-lines indicate the significance threshold at 1.3 (−log(0.05)). Eight AHR-responsive genes were affected by splicing and expression with BaP treatment (E). Distribution of DE and DEU genes across zebrafish chromosomes is shown in (F). Full color version available online.