Phenotypically Anchored mRNA and miRNA Expression Profiling in Zebrafish Reveals Flame Retardant Chemical Toxicity Networks

Abstract

The ubiquitous use of flame retardant chemicals (FRCs) in the manufacture of many consumer products leads to inevitable environmental releases and human exposures. Studying toxic effects of FRCs as a group is challenging since they widely differ in physicochemical properties. We previously used zebrafish as a model to screen 61 representative FRCs and showed that many induced behavioral and teratogenic effects, with aryl phosphates identified as the most active. In this study, we selected 10 FRCs belonging to diverse physicochemical classes and zebrafish toxicity profiles to identify the gene expression responses following exposures. For each FRC, we executed paired mRNA-micro-RNA (miR) sequencing, which enabled us to study mRNA expression patterns and investigate the role of miRs as posttranscriptional regulators of gene expression. We found widespread disruption of mRNA and miR expression across several FRCs. Neurodevelopment was a key disrupted biological process across multiple FRCs and was corroborated by behavioral deficits. Several mRNAs (e.g., osbpl2a) and miRs (e.g., mir-125b-5p), showed differential expression common to multiple FRCs (10 and 7 respectively). These common miRs were also predicted to regulate a network of differentially expressed genes with diverse functions, including apoptosis, neurodevelopment, lipid regulation and inflammation. Commonly disrupted transcription factors (TFs) such as retinoic acid receptor, retinoid X receptor, and vitamin D regulator were predicted to regulate a wide network of differentially expressed mRNAs across a majority of the FRCs. Many of the differential mRNA-TF and mRNA-miR pairs were predicted to play important roles in development as well as cancer signaling. Specific comparisons between TBBPA and its derivative TBBPA-DBPE showed contrasting gene expression patterns that corroborated with their phenotypic profiles. The newer generation FRCs such as IPP and TCEP produced distinct gene expression changes compared to the legacy FRC BDE-47. Our study is the first to establish a mRNA-miR-TF regulatory network across a large group of structurally diverse FRCs and diverse phenotypic responses. The purpose was to discover common and unique biological targets that will help us understand mechanisms of action for these important chemicals and establish this approach as an important tool for better understanding toxic effects of environmental contaminants.

Keywords: flame retardants (additives, reactives); gene expression; mRNA; micro-RNA (miRNA); network analysis; neurodevelopment; transcription factors (TFs); zebrafish.

FIGURE 1

Chemical name, structure, abbreviation, CAS number, physicochemical classes, EC₈₀ s and exposure concentrations (“Conc”) of flame retardant chemicals (FRCs) used in this study. “Conc” represents concentration used for exposures in the study. For FRCs, EC₈₀s [exposure concentrations demonstrating 80% morphological effects based on Truong et al. (2020)] were used as exposure concentrations, except BDE-47, TCEP and TCPP where limit concentration (85 μM) or TBBPA-DBPE where a TBBPA-matched concentration (4 μM) was used.

FIGURE 2

Phenotypic screening, mRNA sequencing and miR sequencing of 10 FRCs. (A) Phenotypic screening with lowest effect level (LEL) of 18 morphological endpoints and four behavioral endpoints; details about these endpoints are included in Supplementary Table 1. Data based on Truong et al. (2020). Exposure concentrations for this study are represented as “Conc.” Black lines with numbers indicate FRC classes based on Figure 1. Colored lines indicate phenotypic response groups: only behavior response; morphology + behavior response; only morphology response and no response. Embryonic photomotor response is included within morphology response and behavior response indicates only larval photomotor response. Numbers within cells represent LELs for each FRC/endpoint combination. (B) Log₂ fold changes of all increased and decreased differentially expressed mRNAs across FRCs. Cutoff was log₂ fold change ≥1.5 and p ≤ 0.05. Numbers represent number of genes with increased (↑) or decreased (↓) mRNA levels. (C) Log₂ fold changes of increased and decreased miR levels across different FRCs. Cutoff was p ≤ 0.05. Numbers represent number of miRs with increased (↑) or decreased (↓) levels.

FIGURE 3

Major gene ontology (GO) processes across FRCs based on (A) differentially expressed mRNAs and (B) target mRNAs of differentially expressed miRs. GO was analyzed using human orthologs within Metacore. For panel (B), differentially expressed miRs were imported into TFmiR and both experimentally validated, and computationally predicted gene targets were imported into Metacore for GO analysis. Data is represented as -log (FDR p value) for each term; a value of ∼1.3 represents FDR p threshold of 0.05. Numbers within cells represent the significant -log (FDR p value). (C) Expression of representative mRNAs known to regulate nervous system development and neurotransmitter activity. Horizontal lines with numbers (1–7) represent FRC class based on Figure 1.

FIGURE 4

Heatmap representing log₂ fold changes of (A) top mRNAs and (B) top miRs across FRCs. Up to 5 genes with highest increase and decrease in mRNA or miR expression levels were selected and their fold changes were plotted for all FRCs. Each column dendrogram color represents a cluster.

FIGURE 5

mRNA-miR interactions across FRCs based on experimentally validated predictions. (A) FRC-mRNA-miR co-regulatory network based on TFmiR gene-miR interactions. Only mRNA-miR pairs with reverse expression directions were considered for the network. FRCs, miRs, miRs that were differentially expressed across 4 or more FRCs, with larger nodes denoting higher number of FRCs. mRNAs are represented within the connected lines. (B) Selected miRs that were decreased across multiple FRCs, with fold changes and major functions (based on GO analysis) of their anti-correlated mRNA targets in specific FRCs. Numbers represent FRC classes based on Figure 1. TBPH and TCPP not represented since there were no miR disruptions.

FIGURE 6

Comparison of TBBPA-DBPE and TCEP for neurotoxic effects. TBBPA-DBPE exposures showed only LPR phenotypes while TCEP showed no phenotype. (A) Heatmap representing mRNA expression for the two FRCs; colored bars on the rows represent gene clusters. (B) GO processes for unique differential mRNAs with increased (Cluster 1) and decreased (Cluster 2) in the TBBPA-DBPE exposures. (C) mRNA-miR network for the two FRCs. FRC, genes. miRs are represented within connected lines.

FIGURE 7

mRNA-transcription factor (TF) interactions across FRCs according to experimentally validated predictions. (A) FRC-mRNA-TF regulatory network based on TFmiR gene-TF interactions. FRCs, TFs. (black nodes with white text)-TFs that were differentially expressed across 5 or more FRCs, with larger nodes denoting higher number of FRCs. (B) Heatmap representing fold changes of selected TFs across all FRCs. Both zebrafish and human orthologs are provided. (C) Heatmap representing mRNA-TF combinations that were co-altered across multiple FRCs. Numbers within cells represent number of FRCs that a specific pair was altered in. Red arrows represent the FRCs for specific pairs that are discussed in the manuscript. All data based on TFmiR experimentally validated predictions.

FIGURE 8

Representative interactions between mRNA, miR and TFs that show a feed forward loop (FFL, dotted circle), a TF regulating a miR (dotted box) and a miR regulating a TF (solid circle) within IPP exposures. TFs, miRs, mRNA. These interactions were selected from the mRNA-miR-TF co-regulatory network for IPP. All interaction data based on TFmiR experimentally validated predictions.

FIGURE 9

mRNA expression for TBBPA vs its derivative, TBBPA-DBPE. (A) Heatmap representing log₂ fold changes of mRNA expression; colored bars on the rows represent gene clusters. Numbers within bars represent cluster numbers for specific mRNA sets. (B) GO processes for unique mRNAs that were increased or decreased uniquely in either FRC.

FIGURE 10

mRNA expression for BDE-47, aryl phosphates and chlorinated phosphates. (A,C) Heatmaps representing log₂ fold changes of mRNA expression; colored bars on the rows represent mRNA clusters. Numbers within bars represent cluster numbers for specific mRNA sets. (B,D) GO processes for unique mRNAs that were increased or decreased uniquely in various clusters.