Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish
- Institute of Marine Research, Norway
- University of Oslo, Norway
- Northwest Fisheries Science Center, National Marine Fisheries Service, United States
- University of Agder, Norway
Abstract
Crude oil spills are a worldwide ocean conservation threat. Fish are particularly vulnerable to the oiling of spawning habitats, and crude oil causes severe abnormalities in embryos and larvae. However, the underlying mechanisms for these developmental defects are not well understood. Here, we explore the transcriptional basis for four discrete crude oil injury phenotypes in the early life stages of the commercially important Atlantic haddock (Melanogrammus aeglefinus). These include defects in (1) cardiac form and function, (2) craniofacial development, (3) ionoregulation and fluid balance, and (4) cholesterol synthesis and homeostasis. Our findings suggest a key role for intracellular calcium cycling and excitation-transcription coupling in the dysregulation of heart and jaw morphogenesis. Moreover, the disruption of ionoregulatory pathways sheds new light on buoyancy control in marine fish embryos. Overall, our chemical-genetic approach identifies initiating events for distinct adverse outcome pathways and novel roles for individual genes in fundamental developmental processes.
https://doi.org/10.7554/eLife.20707.001eLife digest
Accidental oil spills are a worldwide threat to ocean life. Fish eggs and larvae are especially vulnerable; therefore oil spills in areas where fish spawn are of great concern. Fish embryos exposed to crude oil grow slower than normal as larvae and juveniles and often show defects in the heart, face and jaw. However, the underlying mechanisms behind these defects are largely unknown.
Working with the Atlantic haddock (Melanogrammus aeglefinus), Sørhus et al. have now examined fish embryos and larvae that had been exposed to crude oil, and identified those genes that were more active or less active than normal. The findings add further support to the idea that exposure to crude oil causes heart and face defects because it interferes with how the cells that develop into these structures use calcium ions. Signals sent via calcium ions are not only important for the contraction of muscle cells, but they are also essential for regulation of some genes. So, by interfering with the circulation of calcium ions, crude oil can have consequences for both how muscles work and how genes are regulated.
Sørhus et al. also report two previously uncharacterized defects. Firstly, genes that help to regulate the ion and water content of the tissues were highly affected in young fish exposed to crude oil. Some of the genes were more active than normal, while others were less active. This finding in particular would explain why oil-exposed embryos often accumulate fluids, and suggests that the larvae may have altered buoyancy too. Secondly, the oil-exposed embryos showed signs of a shortage of cholesterol and other fatty molecules. This is most likely because they absorbed less material from their yolk, which could also explain why larvae exposed to crude oil grow more slowly than normal.
Finally, in the future, these newly identified genes connected to crude oil toxicity could be used as diagnostic markers to confirm oil-induced injury in fish, and monitor the health of fish populations in the ocean.
https://doi.org/10.7554/eLife.20707.002Introduction
Catastrophic oil spills, rising water temperatures, ocean acidification, and other large-scale anthropogenic forcing pressures impact the health and survival of myriad marine species in ways that are often unknown. Mechanistic relationships between environmental stress and adverse health outcomes are most readily studied in model laboratory organisms. However, domesticated experimental models may be relatively insensitive to real-world environmental change. Emerging genomic technologies, including high-throughput RNA sequencing (RNA-seq), are providing new opportunities to profile physiological stress in wild, non-model species at a transcriptional level. This approach is premised on the full anchoring of gene expression to physiological or morphological injury phenotypes.
The early life stages of marine fish are particularly vulnerable to pollution and other stressors. However, developmental analyses in wild species can be challenging due to limited access to embryos and larvae. Also limited are molecular and cellular tools for imaging-specific structures (e.g. fluorescent protein-expressing transgenes). This includes, for example, a lack of species-specific probes for visualizing gene expression via in situ hybridization. On the other hand, one of the major vertebrate models for studying developmental genetics is a fish. Zebrafish have been a focus for high-throughput experimental techniques for more than three decades. This has yielded a wealth of information on embryonic gene expression patterns, with comparable data often available in chick and mouse embryos. The zebrafish platform therefore provides a powerful mechanistic context for anticipating environmental health impacts in marine fish spawning habitats.
Here, we use an RNA-seq approach to assess the effects of crude oil on the early life stages of a cold-water marine species, Atlantic haddock (Melanogrammus aeglefinus). Crude oil spills such as the 1989 Exxon Valdez (Prince William Sound), 2002 Prestige (Spain), and 2010 Deepwater Horizon (Gulf of Mexico) continue to threaten fisheries worldwide. Haddock are commercially valuable in Norway and other North Atlantic countries, and they spawn in areas proposed for future crude oil production (e.g. the Lofoten archipelago in Nordland). Similar to many other fish species, haddock embryos are vulnerable to developmental defects from crude oil exposure (Norwegian Sea crude; [Sørhus et al., 2016b]). Moreover, we recently showed that RNA-seq applied to normally developing haddock clearly anchored organogenesis phenotypes to the expression of genes involved in determination and differentiation (Sørhus et al., 2016a).
Crude oils are complex chemical mixtures, and fish early life stages are especially vulnerable to component polycyclic aromatic hydrocarbons (PAHs) and their alkylated homologues (Carls and Meador, 2009; Adams et al., 2014). Crude oil-derived PAHs containing three rings disrupt the normal form and function of the embryonic heart, and circulatory failure causes a range of secondary defects (Incardona et al., 2004, 2005). For individual heart muscle cells, the cardiotoxic mechanism involves a blockade of the repolarizing potassium efflux and a reduction in intracellular calcium cycling (Brette et al., 2014). The consequent disruption of excitation-contraction (E-C) coupling leads to rhythm and contractility defects at the scale of the developing heart (Incardona et al., 2009, 2014; Sørhus et al., 2016b). Mechanisms underpinning morphological defects in other embryonic tissues are still poorly understood.
Based on conventional measures of cardiac function and embryolarval anatomy, zebrafish and haddock respond to Norwegian Sea crude oil in ways that are similar and dissimilar. Both species show characteristic abnormalities including bradycardia, reduced chamber contractility, and fluid accumulation in the vicinity of the heart (edema). This suggests an understanding of zebrafish developmental genetics will inform the interpretation of changing messenger RNA (mRNA) levels in crude oil-exposed haddock as determined by RNA-seq. This is particularly true for tissue-specific patterns of gene expression that are highly conserved across vertebrates—for example, genes involved in cardiac organogenesis.
Relative to zebrafish, however, haddock are sensitive to much lower concentrations of crude oil and also display a distinct suite of craniofacial defects that cannot be attributed to circulatory failure (Sørhus et al., 2016b). There are several reasons to expect divergent effects of crude oil on marine fish embryos and larvae. These are attributable to differences in physiology and life history. For example, accumulation of cardiac edema is a canonical form of crude oil toxicity in both freshwater and marine species. Yet marine embryos are hyposmotic to the surrounding water and hence expected to lose water with circulatory failure. This suggests that PAHs may have distinct impacts on ionoregulation and related processes. Also, unlike zebrafish, many pelagic marine embryos are buoyant and have a characteristic morphology not found in species with demersal (sinking) eggs and larvae. Shelbourne first described a relationship between this unique morphology of pelagic fish larvae, osmoregulation and buoyancy control in the mid-twentieth century (Shelbourne, 1955, 1956, 1957), but there has been little progress in the decades since, particularly at a molecular scale.
Understanding cause-effect relationships between exposure to environmental contaminants like crude oil and adverse impacts on organismal health are critical for the construction of adverse outcome pathways (AOPs). The development and application of AOPs is an ongoing movement to improve risk assessments. AOPs are derived from detailed toxicological cause-and-effect relationships that span multiple levels of biological organization, ideally from molecular initiating events to species, community or ecosystem scale responses of regulatory concern (e.g. reduction in a fisheries abundance target). AOPs are widely used in risk assessments for both human and environmental (ecological) health (Ankley et al., 2010; Kramer et al., 2011; Villeneuve et al., 2014; Garcia-Reyero, 2015) Our long-term aim is to develop AOPs specific to oil spills and fish populations, premised on well-studied early life stage toxicity. AOPs based on crude oil cardiotoxicity in developing fish are already fairly well constructed (Incardona and Scholz, 2016) but currently lack details at the molecular level at several steps, particularly in relation to cardiac dysmorphogenesis. We anticipate that identification of changes in gene expression associated with oil-induced developmental defects will further complete these AOPs and expand the molecular toolkit for quantifying oil spill impacts.
In the present study, we used visible developmental abnormalities as phenotypic anchors for evaluating changes in haddock gene expression. We sequenced the full haddock transcriptome at several time points during and after embryonic and larval crude oil exposures. This approach allowed us to explore gene regulation in association with three distinct phenotypes: (1) heart form and function defects, (2) craniofacial deformities, and (3) fluid balance abnormalities and the characteristic pelagic larval form. We also identify perturbations in cholesterol homeostasis linked to poor absorption of yolk as a novel form of crude oil toxicity in marine fish embryos. Our findings are interpreted in the context of highly conserved gene regulation in zebrafish and other vertebrates.
Results
Structure of pelagic larvae and visible phenotypes associated with crude oil exposure
At a rearing temperature of 7°C, haddock embryos began hatching at 12 days post-fertilization (dpf). Unlike zebrafish that complete gastrulation (epiboly) before segmentation (somitogenesis) begins, haddock embryos begin forming anterior somites at about 50% epiboly (3 dpf). They subsequently reach the tailbud stage by 6 dpf (~30 somites), have a regular heartbeat by 8 dpf and completion of organogenesis at 10 dpf (Fridgeirsson, 1978; Hall et al., 2004; Sørhus et al., 2016a). Haddock yolk sac larvae have the characteristic morphology associated with ichthyoplankton from pelagic marine habitats (Figure 1A), marked by a large marginal finfold that surrounds the larva nearly completely on both the dorsal and ventral sides. The outer epidermis is thus separated from the brain, main body axis muscles, and internal organs by a voluminous subdermal space. This space is filled with extracellular matrix (Morrison, 1993) and is continuous with an avascular yolk sac sinus, with connections between the dorsal space and the ventral yolk sac in the vicinity of the pectoral fin (Shelbourne, 1955). The subdermal space acts as a reservoir for water balance in order to maintain larval buoyancy (Shelbourne, 1955, 1956, 1957), with specialized cells regulating ion and water balance (ionocytes or mitochondria rich cells, MRCs) distributed throughout the adjacent epidermis (Shelbourne, 1957; Hirose et al., 2003; Hiroi et al., 2005).
Figure 1 with 1 supplement see all
Terminal phenotypes after high dose exposure.
Control (A) and exposed (B) three days post hatch (dph) larvae (6 days post embryonic exposure). Open arrowheads in (A) indicate the marginal finfold surrounding the larvae and the white asterisk indicate the location of the connection between the dorsal space and the ventral yolk sac in the vicinity of the pectoral fin. In (B) the black arrowhead indicates severely reduced craniofacial outgrowth, while the black arrow indicates yolk sac edema. The ventricle and atrium in control (C) and embryonically exposed (D) animals are indicated by black and white arrows, respectively. (E) Normal craniofacial structure in control, and (F) moderate and (G) severe craniofacial defects in exposed animals. (H) Normal marginal finfold in control, (I) exposed animals with severe reduction of anterior marginal finfold (open arrowheads). Yolk mass (*) in control (J) and embryonically exposed larvae (K). (L) Control and (M) exposed 18 dph larvae. Open arrowheads indicate increased anterior marginal finfold, black arrowhead indicates reduced upper jaw outgrowth, and black arrow indicates edema formation in the peritoneal cavity in oil-exposed larvae (M). Scale bar: 0.2 mm (C,D; E–G; H–K) and 1 mm (A,B and L,M).
Haddock embryos were continuously exposed to crude oil at environmentally relevant total PAH concentrations of 6.7 ± 0.2 µg/L (high dose) and 0.58 ± 0.05 µg/L (low dose), and intermittently at 6.1 µg/L (pulse dose). This yielded internal total PAH doses of 3.0 ± 1.3 µg/g wet weight, 0.19 ± 0.02 µg/g, and 0.22 ± 0.06 µg/g, respectively (see [Sørhus et al., 2016b] for experimental details). Although the low and pulsed exposures led to similar total PAH accumulation in embryos, phenotypes were slightly more severe in the latter due to relatively higher exposure concentrations during critical windows of early development (Sørhus et al., 2016b). The embryonic exposure began at 2 dpf and ended shortly after the end of organogenesis at 10 dpf, just prior to hatch. Larvae were exposed to the same regimen from day of hatch to 18 days post hatch (dph). As expected, tissue PAH accumulation was lower than for embryos, with 0.81 ± 0.18 µg/g, 0.086 ± 0.015 µg/g, and 0.081 ± 0.024 µg/g, for high, low, and pulse doses, respectively. Except where indicated below, phenotypes were quantified as described in detail previously (Sørhus et al., 2016b), with 96% of high-dose animals showing abnormal phenotypes, ranging to ~60% for pulse dose and ~35% for the low dose. Representative terminal phenotypes in high dose are shown in Figure 1 for the embryonic exposure at 3 dph and larval exposure at 18 dph. Grossly, as in other species, crude oil exposure led to defects in cardiac function and morphology and accumulation of edema around the heart and in the yolk sac (Figure 1B). Defects in cardiac morphology included a failure to properly loop the atrial and ventricular chambers from a linear to an adjacent orientation, and reduced size of the ventricle (Figure 1C,D). In addition, oil-exposed haddock embryos displayed craniofacial defects in their larval stages that ranged in severity (Figure 1E–G), from marked reductions in upper jaw/skull base structures (Figure 1F) to near complete lack of upper and lower jaws (Figure 1G). Moreover, the anterior portion of the dorsal marginal finfold was collapsed or missing and the hindbrain ventricle typically failed to fill with cerebrospinal fluid in embryonically exposed larvae with severe edema (Figure 1H,I). Finally, in more severely affected embryos, a failure of yolk absorption was obvious at 3 dph (Figure 1J,K). Even in mildly affected embryos, yolk absorption was reduced after hatch as assessed by measuring the two-dimensional area of the yolk mass in lateral images (control yolk area control 0.63 ± 0.06 mm2, low-dose group 0.90 ± 0.11, high-dose group 1.2 ± 0.3; mean ± s.d., ANOVA p<0.0001). In contrast, there was no measurable difference in yolk area at day of hatch. After larval exposure, the primary morphological defects were reduced jaw growth and edema accumulation (Figure 1L,M), the latter in the peritoneal cavity. In contrast to embryos, the dorsal anterior subdermal space accumulated fluid in larvae and did not collapse.
Abnormal phenotypes relating to the formation of edema, heart function and morphogenesis, craniofacial structure, and yolk absorption manifested at different developmental time points during embryonic exposure and afterwards when embryos were transferred to clean water for hatching (Figure 2). Samples were collected for transcriptome profiling at four embryonic stages (E1-4, Figure 3A) and at two stages post-hatch (E5-6, Figure 3A). At 6 dpf (E2 sampling point; ~30 somites), embryos exposed to the high dose were indistinguishable from controls (Figure 2A). By 8 dpf, small accumulations of edema could be observed in the yolk sac of oil-exposed embryos, but their hindbrain ventricles were ‘inflated’ with fluid (Figure 2B). By 10 dpf (E3 sampling point), edema was evident in most embryos, typically filling the space above the yolk between the anterior of the head and the tail, and hindbrain ventricles lacked fluid (Figure 2C). Similarly, at 6 dpf/E2, heart development appeared unaffected in oil-exposed embryos and was at the un-rotated midline cone stage (Figure 2D). By 9 dpf (one day before sample E3), hearts in both control and high-dose-exposed embryos had rotated and were beginning to loop, but ventricular walls already appeared slightly thinner (Figure 2E) and heart rate was significantly slower (20 + 6 beats/min compared to 26 ± 3 beats/min for controls; [Sørhus et al., 2016b]). By day of hatch (E5 sampling), a high percentage (54%) of exposed embryos showed un-looped hearts with small, silent ventricles (Figure 2F).
Figure 2
Appearance of phenotypes over time.
In each panel control and high-dose-exposed embryos are shown on the left and right, respectively. (A–C) Lateral overview of whole embryos showing accumulation of edema (anterior to the left). (A) 6 dpf/E2 sampling point. (B) 8 dpf (between E2 and E3 sampling points). Heart (h) and liver bud (l) are indicated. White arrowheads indicate outer margins of the yolk sac membranes; asterisk indicates small pocket of edema. Black arrowheads indicate the hindbrain ventricle. (C) 10 dpf/E3 sampling point. Arrowheads same as (B); asterisks indicate expanded yolk sac edema. (D–E) High-magnification ventral views of the heart (anterior at top). (D) 6 dpf/E2. Dashed turquoise lines indicate outer border and lumen of midline cardiac cone. (E) 9 dpf (between E2 and E3). Arrows indicate the atrium (a), ventricle (v) and bulbus arteriosus (ba). (F) 0 dph (E5 sampling point). Chambers indicated as in (E). (G–I) Lateral views of the developing head (anterior to the left). (G) 8 dpf (between E2 and E3). (H) 0 dph (E5). Arrow indicates abnormal lower jaw cartilages in oil-exposed larva. (I) 3 dph (E6 sampling point). Red bars indicate difference in eye diameter between control and exposed larvae.
Figure 3 with 2 supplements see all
Exposure regimes and differentially expressed genes (DEGs) during embryonic development.
(A) Embryos were exposed to a continuous high dose (black line; 6.7 ± 0.2 μg/L TPAH), a pulsed dose (red line; 0.09 ± 0.02–6.8 ± 1.0 μg/L TPAH) and a continuous low dose (blue line; 0.58 ± 0.05 μg/L TPAH) of crude oil. Photos indicate normal developmental progress at each of six sampling time points (E1–E6). Venn diagrams show shared and exclusive DEGs for each of the three oil exposures at E1–E6. (B) Venn diagrams illustrating the number of shared and exclusive DEGs at each stage in development up to hatching for the three exposure regimes.
Onset of craniofacial abnormalities took a longer course. At 8 dpf (2 days before sample E3), head structures appeared identical in control and exposed embryos (Figure 2G). At hatch (E5 sample), jaw structures appeared somewhat abnormal (Figure 2H) but became much more strikingly severe by 3 dph (E6 sample; Figure 2I). Notably, the eyes appeared smaller in exposed embryos by 3 dph (Figure 2I). We did not quantify this effect, because it was demonstrated earlier in zebrafish that small eyes result from loss of cardiac function by either genetic or chemical means (Incardona et al., 2004), and hence, this phenotype is not specific to crude oil toxicity.
Other organs and structures were apparently unaffected by oil exposure. For example, the development of the liver and lateral line neuromasts progressed normally even in the most severely impacted larvae that were exposed as embryos (Figure 1—figure supplement 1A–D).
Oil-induced changes in gene expression during embryonic development
Relative to unexposed controls, differently expressed genes (DEGs) in oil-exposed fish were defined as having significantly (p<0.05) higher or lower levels of expression. The number of exclusive and shared DEGs varied across exposure regime and haddock developmental age (Figure 3B). After 24 hr of oil exposure (sampling stage E1; 3 dpf), relatively few genes were differentially expressed, and most were significantly downregulated (Supporting dataset 1, Sørhus et al., 2017). From sampling point E2 (6 dpf) through E6 (three days post hatch, dph), however, the number of DEGs was substantial, particularly in response to the high dose exposure. With the exception of the initial sampling point (E1), a total of 28 DEGs were shared across all embryonic stages (Table 1). As expected, the largest category (nine DEGs) included genes associated with stress response and xenobiotic metabolism. The remaining genes play a role in tyrosine catabolism, myofibrillar establishment and cardiac tissue repair, central nervous system (CNS) function and degeneration, ATP metabolism, and cholesterol synthesis.
Table 1
Genes expressed at all stages during and after embryonic exposure (E2–E6) in high dose group. SP, swissprot; GB, genebank; IE; increased expression; DE; decreased expression.
Cod ID | Swissprot annotation | SP ID | GB ID | Category | Regulation |
|---|---|---|---|---|---|
ENSGMOG00000018302 | Fumarylacetoacetase | faaa | fah | Tyrosine metabolism | IE |
ENSGMOG00000000318 | Cytochrome P450 1A1 | cp1a1 | cyp1a1 | xenobiotic metabolism and stress | IE |
ENSGMOG00000012518 | Glutathione S-transferase P | gstp1 | gstp1 | xenobiotic metabolism and stress | IE |
ENSGMOG00000016016 | Glutathione S-transferase omega-1 | gsto1 | gsto1 | xenobiotic metabolism and stress | IE |
ENSGMOG00000018752 | 3-hydroxyanthranilate 3,4-dioxygenase | 3hao | haao | xenobiotic metabolism and stress | IE |
ENSGMOG00000006796 | 3-beta-hydroxysteroid-Delta(8),Delta(7)-isomerase | ebp | ebp | xenobiotic metabolism and stress | IE |
ENSGMOG00000007636 | Glutamine synthetase | glna | glul | xenobiotic metabolism and stress | IE |
ENSGMOG00000015234 | Heat shock protein HSP 90-alpha | h90a1 | hsp90a.1 | xenobiotic metabolism and stress | IE |
ENSGMOG00000012029 | Peptidyl-prolyl cis-trans isomerase | ppia | - | xenobiotic metabolism and stress | IE |
ENSGMOG00000000218 | Ammonium transporter Rh type A OS=Mus | rhag | rhag | xenobiotic metabolism and stress | Mainly IE |
ENSGMOG00000003353 | Ferritin, middle subunit | frim | - | xenobiotic metabolism and stress | IE |
ENSGMOG00000018206 | Filamin-C | flnc | Flnc | myofibrillar establishment and repair | IE |
ENSGMOG00000001317 | Iron-sulfur cluster assembly enzyme ISCU, mitochondrial | iscu | Iscu | cardiac defects | IE |
ENSGMOG00000010446 | Fatty acid-binding protein, heart | fabph | fabp3 | cardiac defects and repair | IE |
ENSGMOG00000007115 | Lanosterol 14-alpha demethylase | cp51a | cyp51a1 | Cholesterol syntheis | IE |
ENSGMOG00000005565 | Squalene monooxygenase | erg1 | Sqle | Cholesterol syntheis | IE |
ENSGMOG00000018991 | Farnesyl pyrophosphate synthase | fpps | fdps | Cholesterol syntheis | IE |
ENSGMOG00000005774 | 3-hydroxy-3-methylglutaryl-coenzyme A reductase | hmdh | hmgcr | Cholesterol syntheis | IE |
ENSGMOG00000015657 | Epididymal secretory protein | npc2 | npc2 | Cholesterol syntheis | IE |
ENSGMOG00000001249 | Putative adenosylhomocysteinase | sahh3 | ahcyl2 | cardiac defects | DE |
ENSGMOG00000013374 | Peptide Y OS=Dicentrarchus | py | - | CNS function and development | IE |
ENSGMOG00000014820 | Complement C1q-like protein | c1ql2 | c1ql2 | CNS function and development | IE |
ENSGMOG00000017148 | Augurin-A OS=Danio rerio | augna | zgc:112443 | CNS function and development | IE |
ENSGMOG00000001072 | C-4 methylsterol oxidase | erg25 | sc4mol | CNS function and development | IE |
ENSGMOG00000013980 | Fatty acid-binding protein, brain | fabp7 | fabp7 | CNS function and development | IE |
ENSGMOG00000014938 | Maltase-glucoamylase, intestinal | mga | mgam | ATP metabolism | IE |
ENSGMOG00000003530 | ADP/ATP translocase | adt3 | slc25a6 | ATP metabolism | DE |
ENSGMOG00000006172 | IEF0762 protein C6orf58 homolog | cf058 | - | not known | IE |
Oil-induced changes in gene expression during larval development
Haddock larvae under the same three exposure regimes were transcriptionally profiled at five distinct developmental stages (L1-5; Figure 4A). Relative to embryos, transcriptional responses to crude oil-exposed haddock larvae were more subtle at 1, 2, and 9 dph (L1-3). This was followed by marked changes in gene expression at 14 dph (L4) for the high dose and for all treatments at 18 dph (L5) (Supporting dataset 2, Sørhus et al., 2017). In the high-dose group, nearly 1000 of the >3000 DEGs at L4 and L5 were shared (Figure 4B). However, for the high-dose treatment, only eight genes were shared across all larval stages. As expected, five of these genes are involved in xenobiotic metabolism (Table 2).
Figure 4 with 2 supplements see all
Exposure regimes and differentially expressed genes (DEGs) during larval development.
(A) Larvae were exposed to a continuous high dose (black line; 7.6 ± 0.7 μg/L TPAH), a pulsed dose (red line; 0.3 ± 0.3–6.1 ± 0.5 μg/L TPAH), and a continuous low dose (blue line; 0.65 ± 0.08 μg/L TPAH) of crude oil. Photos indicate normal developmental progress at each of five sampling time points (L1–L5). Venn diagrams show shared and exclusive DEGs for each of the three oil exposures at L1–5. (B) Venn diagrams illustrating the number of shared and exclusive DEGs at each larval stage for the three exposure regimes.
Table 2
Genes expressed at all stages during larval exposure (L1–L5) in high-dose group. SP, swissprot; GB, genebank; IE; increased expression; DE; decreased expression.
Cod ID | Swissprot annotation | SP ID | GB ID | Category | Regulation |
|---|---|---|---|---|---|
ENSGMOG00000009114 | Aryl hydrocarbon receptor repressor | ahrr | ahrr | Xenobiotic metabolism | IE |
ENSGMOG00000020141 | Cytochrome P450 1B1 | cp1b1 | cyp1b1 | Xenobiotic metabolism | IE |
ENSGMOG00000006842 | Cytochrome P450 1B1 | cp1b1 | cyp1b1 | Xenobiotic metabolism | IE |
ENSGMOG00000019790 | Cytochrome P450 1B1 | cp1b1 | cyp1b1 | Xenobiotic metabolism | IE |
ENSGMOG00000000318 | Cytochrome P450 1A1 | cp1a1 | cyp1a1 | Xenobiotic metabolism | IE |
ENSGMOG00000014967 | Keratinocyte growth factor | fgf7 | fgf7 | Myocardial development and tissue repair | IE |
ENSGMOG00000020500 | Forkhead box protein Q1 | foxq1 | foxq1 | Transcription factor | IE |
ENSGMOG00000000218 | Ammonium transporter Rh type A | rhag | rhag | Gas transport | IE |
General patterns of gene expression in response to crude oil
Read count data from the RNA-Seq closely matched expected abundances based on tissue-specific expression patterns for orthologous genes in zebrafish, available in the expression database at www.zfin.org (Supplementary file 1A), and generally correlated with the overall mass of the contributing tissues. Genes known to have tightly restricted cardiac expression generally had read count values below 100, with bmp10 and kcnh2 just above detection limits (10 reads). At the onset of expression in zebrafish, bmp10 transcripts are detected by in situ hybridization in perhaps fewer than 100 cells (Laux et al., 2013). In contrast, bmp4 is more widely expressed in zebrafish at the segmentation stage, including the eye, tailbud, and epidermis in addition to the heart. As expected, this gene had a correspondingly higher read count (265) in haddock. Genes expressed more strongly throughout the entire heart tube had read counts above 60 (e.g. nkx2.5 at 85 and fhl2 at 65). The cardiac-specific Serca2 isoform (atp2a2) had a read count of 312, while the isoform expressed in skeletal muscle (atp2a1) had a read count of 9114. Similarly, read counts for the atrial-specific isoform of myosin heavy chain (myh6) and the major skeletal muscle isoform myh1 were 176 and 2543, respectively. Genes expressed in the neural tube, a tissue mass much larger than the heart but less than the myotomes, had intermediate read counts that also fit with their relative expression patterns. For example, pax6 and nkx2.2 had read counts of 1773 and 333, respectively, with pax6 expressed in a fairly wide dorsal domain of the neural tube, and nkx2.2 expressed in a narrower ventral region.
We characterized pathways affected by oil exposure using three methods: extensive manual curation, KEGG Pathway Mapping, and Ingenuity Pathway Analysis (IPA; see Materials and methods). As detailed in the following sections, our manual curation identified specific patterns of gene expression in the context of cardiotoxicity, craniofacial deformities, disrupted ion and water balance, and disrupted cholesterol homeostasis. These same pathways were identified with statistical rigor using IPA and KEGG. Moreover, IPA demonstrated enrichment for these pathways at developmental time points that preceded the onset of visible phenotypes, and a lack of enrichment for pathways associated with structures that were phenotypically normal (Table 3).
Table 3
Time course of pathway enrichment relating to affected and unaffected developmental and functional phenotypes.
Phenotype† | Development stage* | |||||
|---|---|---|---|---|---|---|
3 dpf/E1 | 6 dpf/E2 | 10 dpf/E3 | 11 dpf/E4 | 0 Dph/E5 | three Dph/E6 | |
Cardiovascular | 0 (0/8) | 22.4 (11/49) | 5.7 (4/70) | 7.0 (4/57) | 4.7 (2/43) | 2.1 (1/48) |
Craniofacial | 0 (0/8) | 12.2 (6/49) | 10 (7/70) | 5.3 (3/57) | 7.0 (3/43) | 2.1 (1/48) |
Liver | 12.5 (1/8) | 0 (0/49) | 5.7 (4/70) | 8.8 (5/57) | 0 (0/43) | 0 (0/48) |
Eye | 0 (0/8) | 4.1 (2/49) | 20 (14/70) | 48.6 (17/35) | 51.2 (22/43) | 50.0 (24/48) |
Osmoregulation | -- | 43.3 (13/30) | 29.3 (12/41) | 15.0 (3/20) | 16 (4/25) | -- |
Cholesterol | 0/30 | 0 (0/27) | 27.1 (13/48) | 31.3 (10/32) | 25.5 (12/47) | -- |
Lipid | 0/30 | 40.7 (11/27) | 35.4 (17/48) | 50.0 (16/32) | 48.9 (23/47) | -- |
-
*Percentage of total enriched pathways (absolute values).
-
†Numbers of affected pathways representing Cardiovascular, Craniofacial, Liver and Eye were extracted from the combined Development category in IPA results; numbers of pathways representing osmoregulation/ion transport were extracted from the Molecular Transport category; numbers of pathways affecting Cholesterol/sterol metabolism and other non-cholesterol lipids (Lipid) were extracted from the Lipid Metabolism category.
At all stages, the IPA subcategory of Organismal Development or Embryonic Development (henceforth combined as Development) was in the top 5 Diseases and Bio Functions category under Physiological System Development and Function with p values ranging from 10−3 to 10−19 (Tables 2 and 3). Counts of the number of pathways specifically involving the cardiovascular system showed that no pathways were affected at E1 (3 dpf, 50% epiboly/cardiac progenitor stage), while the heart represented 22% of the affected Development pathways at 6 dpf/E2, the cardiac cone stage at which there was no visible phenotype (Table 3). The number of cardiovascular pathways fell to 5.7% at 10 dpf/E3, by which point hearts were visibly abnormal, falling to only one or two affected pathways at hatching stages (E5 and E6).
Pathways specifically related to head, face, or skull development were similarly enriched at all stages except 3dpf/E3, representing 12% and 10% at 6 dpf/E2 and 10 dpf/E3, prior to visible differences in head structures. Importantly, at 6 dpf/E2, prior to the onset of both visible cardiac and craniofacial defects, the top 10 enriched pathways under Development included three involving head development and two involving heart development (p values 10−6 to 10−12; Supplementary file 1B). In contrast, pathways relating to liver development were enriched at 5.7% and 8.8% at only two time points, 6dpf/E3 and 10 dpf/E4, and these did not appear in the top 10. Moreover, inspection of individual DEGs associated with those pathways showed genes involved in lipid transport rather than bona fide regulators of liver development (see below). The single pathway relating to the liver at 3 dpf/E1 was represented by a single gene, cyp1a. At these stages, these lipid transport genes are most strongly expressed in the yolk syncytial layer. Notably, IPA also detected larger scale enrichment of eye genes, almost all down-regulated, accounting for roughly 50% of developmental pathways during pigmentation of the retina, but prior to obvious differences in eye sizes after hatch.
Genes associated with osmoregulation were identified by IPA under the Molecular Transport category (Diseases and Bio Functions, Molecular and Cellular Functions). We quantified pathways relating to specific ions (e.g., Na+, K+), inorganic ions, and metals (Table 3). At E1/3 dpf Molecular Transport was not in the top five affected pathways, but became enriched at 43% at E2/6 dpf, with primarily down-regulation of genes prior to the onset of visible edema. These Molecular Transport pathways remained significantly enriched (29%, 15%, 16%) until onset of hatch (E5). By 3 dph/E6, Molecular Transport pathways dropped below the top 5.
Pathways related to cholesterol and other lipids (e.g. phospholipids, fatty acids) followed a similar pattern as osmoregulation. As for Molecular Transport, Lipid Metabolism was consistently in the top five Molecular and Cellular Functions category. Pathway enrichment was overall at the highest levels for Lipid Metabolism. We separately quantified individual pathways relating to sterols (e.g. cholesterol synthesis, cholesterol transport, steroid biogenesis) and other fundamental (non-signaling) lipids (e.g. fatty acid synthesis and transport, glycerolipids, phospholipids) (Table 3). General lipid metabolism pathways were highly enriched at 6 dpf/E2 (41%) and remained high (35–50%) until 3 dph/E6 when Lipid Metabolism pathways fell below the top 5. Cholesterol metabolism pathways were first enriched at 10 dpf/E3 at 27%, remaining at 31% and 26% until 3 dph/E6, when they also fell below the top 5. Notably, all Lipid Metabolism pathways were enriched prior to measureable detection of reduced yolk absorption at 3 dph.
Pathway enrichment was dose-dependent and clearly associated with the frequencies of abnormal phenotypes (Supplementary file 1C). For example, the combined general Development categories at 6 dpf included 203, 121, and 0 pathways (among the top five general categories) for the high, pulse, and low doses, respectively. At this stage, Molecular Transport pathways were enriched at levels of 117, 80, and 0 for the high, pulse, and low doses, respectively. At 10 dpf/E3, the Lipid Metabolism category included 115, 28, and 5 pathways for the high, pulse, and low doses, respectively. Finally, Cardiotoxicity pathways were prominent for nearly all times points for each dose (Supplementary file 1C). For high, pulse, and low doses, numbers of enriched pathways at 6 dpf/E2 were, respectively, 63, 33, and 17; at 10 dpf/E3, 79, 10, and 5; at 11 dpf/E4, 92, 40, and 9; at 0 dph/E5, 69, 46, and 26; and at 3 dph/E6, 34, 19, and 12.
We identified 10 individual genes (Supplementary file 1D) and KEGG pathways (Figure 3—figure supplement 1 and Figure 4—figure supplement 1) that were the most highly responsive (highest positive or negative fold change) to the high oil treatment regime relative to unexposed control fish. Briefly, in both embryos and larvae, DEGs involved in xenobiotic metabolism and stress response were highly represented. Genes involved in the development and function of neural networks and cholesterol/steroid biosynthesis were upregulated, while genes involved in intracellular calcium signaling were primarily downregulated.
In order to investigate both pathways with numerous and few genes, we chose two different approaches for KEGG pathways analysis (1) Total: Pathways with the highest number of DEGs ≥ 2 FC (Figure 3—figure supplement 1A and Figure 4—figure supplement 1A) and (2) Normalised: Pathways with the largest fraction of DEGs ≥2 FC/ Total number of genes in pathway) (Figure 3—figure supplement 1B and Figure 4—figure supplement 1B).
During embryonic exposure pathways associated with PAH metabolism were represented among the most affected. Indicative of disrupted osmoregulation and ion channel blockade, secretion pathways and calcium signaling showed decreased expression at the earliest stages. Further, we observed increased expression in steroid metabolism and biosynthesis pathways suggesting an effect on cholesterol metabolism (Figure 3—figure supplement 1B). Post exposure, we observed increased expression of several pathways suggestive of an inflammatory response (protein digestion and degradation, influenza A, antigen processing and presentation), while expression of genes in phototransduction pathway was decreased (Figure 3—figure supplement 1A).
During larval exposure at the first three sampling stages a small number of genes, and thus, pathways were regulated and most were participating in PAH metabolism. Consistent with total number of DEGs (Figure 4B), stage L4 and L5 included pathways with higher number of DEGs ≥ 2 FC. Most noticeable was the decreased expression in calcium signaling pathway and hypertrophic cardiomyopathy (HCM) pathway at 14 dph and decreased expression in pancreatic secretion and protein digestion and absorption (Figure 4—figure supplement 1A) and steroid biosynthesis pathways (Figure 4—figure supplement 1B) at 18 dph.
Finally, four genes stood out as unique for (1) being highly upregulated in both embryos and larvae, (2) their non-affiliation with a larger network or pathway, and (3) their potential connections to visible phenotypes. These included collagen and calcium-binding EGF-like domain 1 (ccbe1), the ammonia transporter rhag, forkhead box transcription factor foxq1, and fibroblast growth factor fgf7 (Table 1, Supplementary file 1D). The following sections identify specific patterns of gene expression in the context of cardiotoxicity, craniofacial deformities, disrupted ion and water balance, and disrupted cholesterol homeostasis.
Genes associated with defects in cardiac function and morphogenesis
As noted above, early formation of the heart was not affected by crude oil exposure, but morphological defects followed after functional defects were first observed at 9 dpf (bradycardia). In addition, morphology became more severely impacted over time, with later defects including failure of looping and poor ventricular growth becoming prominent by 0 dph/E5. There are several possible etiologies for ventricular size reduction. For example, although the precise mechanism by which intracellular calcium regulates embryonic cardiomyocyte proliferation is still unknown, a disruption of calcium cycling could reduce proliferation, thereby yielding fish with smaller hearts (Rottbauer et al., 2001; Ebert et al., 2005). We therefore focused on genes involved in cardiac morphogenesis. The earliest alteration was a fourfold increase in the signaling molecule, bmp10 at 6 dpf/E2 while the heart was at the midline cone stage, and appeared unaffected in oil-exposed embryos (Figure 5, Supplementary file 1E). IPA also identified Bmp signaling as a significantly enriched pathway at this time point, under the Organismal Development category. Elevation of bmp10 was followed by the upregulation of the cardiac transcription factor nkx25, to a threefold increase at 10 dpf/E4 and a nearly sixfold increase at 11dpf/E5, when the heart was beating regularly. At hatch (0 dph), the expression of the transcription factor tbx3 was elevated eightfold. Lastly, atrial natriuretic factor (nppa), a key homeostatic regulator of contractility, was downregulated by 2.3-fold in larvae at 3 dph (Figure 5, Supplementary file 1E). Notably, overexpression of bmp10, nkx25, or tbx3 is associated with serious heart defects in other vertebrates (Chen et al., 2006; Ribeiro et al., 2007; Tu et al., 2009).
Figure 5
DEGs involved in cardiogenesis.
Regulation of genes involved in cardiogenesis during and after embryonic exposure. Purple: increased expression, red: decreased expression in exposed group.
Crude oil exposures caused functional defects in the developing haddock heart, in the form of bradycardia, ventricular asystole and decreased contractility in embryos and partial atrio-ventricular conduction blockade in larvae (Sørhus et al., 2016b). This is consistent with disruption of the rhythmic fluxes of Ca2+ and K+ ions that regulate E-C coupling in heart muscle cells (Brette et al., 2014). We therefore focused on genes associated with cardiomyocyte membrane potential and intracellular Ca2+ cycling—for example, sarcoplasmic reticulum calcium ATPases (SERCAs) and the ryanodine receptor (RyR) ([Sørhus et al., 2016a], Supplementary file 1F and 1G ). We found three paralogs for at2a2 (Serca2) that were present at very different read count values (~300, 900, and 4000 at 6 dpf/E2). The two more abundant paralogs were transiently down-regulated in oil-exposed embryos at 6 dpf/E2, prior to the onset of functional and morphological defects, while the third paralog was down-regulated at 0 dph/E5 (Figure 6, Supplementary file 1H). Similarly, there were four nac1 paralogs that were all low abundance, and one was transiently down-regulated with the at2a2 genes at 6 dpf. Finally, the kcnh2 gene contributing to the repolarizing K+ current was down-regulated ninefold at hatch/E5. There were effects on a few other E-C coupling genes, but these had very high read counts, and are therefore likely to be associated with skeletal muscle. These included two at2a1 (serca1) paralogs that had opposite responses, and atpa.
A different picture emerged from the larval exposure. Changes in expression of E-C coupling genes occurred after the onset of functional defects (AV block arrhythmia). No changes in cardiac E-C coupling genes were observed at the L3/9 dph time point when larvae showed AV block. Six days later (L4), there was fourfold down-regulation of a nac1 paralog with the lowest read count value and an at2a2 paralog with the highest value. At this stage there was up-regulation of two kcnj2 paralogs (encoding potassium channels), two high abundance at2a1 paralogs, and a low abundance scn2a paralog. At 18 dph/L5, one paralog each of at2a1 and kcnj2 remained elevated, while atpa was elevated, and the kcnj12 potassium channel gene were down-regulated, the latter strongly (~6 fold).
Genes associated with craniofacial abnormalities
Craniofacial structures that shape the head include cartilage derived from neural crest cells and muscles that develop from paraxial mesoderm. Neural crest cells migrate from the anterior neural tube to form the pharyngeal arches with both dorsal (upper jaw) and ventral (lower jaw) patterning. They then differentiate into chondrocytes (Knight and Schilling, 2006; Simões-Costa and Bronner, 2015) and grow by processes such as convergence-extension (Shwartz et al., 2012; Kamel et al., 2013). Concurrently, mesodermal cells differentiate into patterned muscle in appropriate association with partner cartilage. Several lines of evidence suggest multidirectional signaling between all associated tissues, including endoderm (i.e. pharyngeal pouches), mesoderm, and overlying ectoderm (Minoux and Rijli, 2010; Medeiros and Crump, 2012; Kamel et al., 2013; Kong et al., 2014). Studies on zebrafish craniofacial mutants have primarily focused on the neural crest cell lineage, with less attention to muscle development or interactions between developing muscle and cartilage (Lin et al., 2013). Defects in oil-exposed haddock were marked by a dose-dependent reduction in more anterior cartilages (Figure 1E–G). This affected the basicranium most severely, with progressive loss of more posterior arch derivatives. Where present, craniofacial cartilage appeared small and distorted, without transformation to dorsal or ventral fates. This morphometry superficially aligns to several zebrafish mutants affecting either neural crest cell (Kimmel et al., 2001; Nissen et al., 2006; Lu and Carson, 2009; Kamel et al., 2013) or muscle development (Hinits et al., 2011; Shwartz et al., 2012). We therefore interpreted developmental changes in haddock gene expression in the context of these well-characterised zebrafish mutants.
The expression patterns of 12 genes with known roles in neural crest cell-dependent craniofacial development were significantly altered in the highest exposure regime (Figure 6). Read counts for these genes were all relatively low, consistent with highly restricted tissue-specific expression patterns (Supporting dataset 1, Sørhus et al., 2017). At 6 dpf, prior to visible craniofacial malformation, we observed lower expression levels of foxi1 (pharyngeal pouches) wnt9b (ectoderm), fgfr2 (chondrocytes), and fgfr3 (chondrocytes) compared to control (Figure 7A). The downregulation of wnt9b and fgfr2 persisted to 10 dpf, together with a downregulation of tgfb3 and two sox9b paralogues, the latter also expressed in neural crest cell-derived chondrocytes. Conversely, edn1, dlx3b, and dlx5a were upregulated at 10 dpf. In zebrafish embryos the two dlx genes are normally expressed in endoderm and arch neural crest cell-derived mesenchyme (Talbot et al., 2010). At 11 dpf, expression levels were down for fgfr2, fzd7a (a receptor for Wnt9b; chondrocytes), and tgfb3 and up for edn1 and dlx3b. None of these genes were differentially expressed after hatching (Figure 7A, Supplementary file 1I).
Figure 6
DEGs involved in E–C coupling.
Embryonic developmental samples (E1–6) were collected during (black lettering) and after (blue lettering) crude oil exposure. Oil exposure was continuous across the larval sampling points (L1–5).
Figure 7
DEGs involved in craniofacial development.
(a) Regulation of genes involved in craniofacial development during and after embryonic exposure. (b) Regulation of myosin heavy chain genes. Purple: increased expression, red: decreased expression in exposed group.
Genes controlling craniofacial muscle patterning are poorly characterised; however, muscle determination factors (e.g., myod, myf6 [Li et al., 2014]) appeared unaffected in oil-exposed haddock. Nevertheless, expression levels were significantly reduced for genes involved in the terminal differentiation of skeletal muscle cells, including several myosin heavy chain genes (myh) (Figure 7B, Supplementary file 1J). These included myh4, myh9, and myh10 paralogues at 6 dpf, myh9 and another paralogue of myh4 and myh10 at 10 dpf, and myh3, myh4 and myh9 at 11 dpf. Only myh3 remained significantly downregulated after hatching (3 dph) relative to controls. Notably, expression of myh1, encoding the major fast myosin heavy chain gene expressed in the body musculature (Thisse et al., 2001), was largely unaffected except for a small reduction at 11 dpf. Other myosin genes specific to muscle groups in the head and trunk (Peng et al., 2002; Elworthy et al., 2008), on the other hand, showed more complex differential expression patterns.
Genes associated with ion and water regulatory imbalance
Fluid accumulation in the form of edema is a hallmark indication of crude oil toxicity in fish embryos. Although patterns of edema formation vary across freshwater and marine fish species (Incardona and Scholz, 2016), it nearly always involves anatomical compartments adjacent to the heart and the yolk sac. However, marine embryos are hyposmolar to surrounding seawater, and they should therefore lose water along a diffusion gradient if osmoregulation is disrupted as a consequence of heart and circulatory failure. In fish embryos and yolk sac larvae, osmoregulation is controlled by MRCs in the epidermis and yolk sac membrane that actively secrete NaCl (specifically Cl-) to maintain an appropriate water and ion balance (Hiroi et al., 2005). Genetic and pharmacologic studies have shown that circulation is required to maintain embryonic MRC cell function. For example, total body osmolality increased in seawater-adapted tilapia embryos with a ~50% reduction in total cardiac output (Miyanishi et al., 2013). Therefore, edema formation in marine species with oil-induced circulatory defects is not a consequence of water moving into the embryo (as it is for freshwater species) but rather water moving along an internal osmotic gradient from peripheral tissues to the vicinity of the heart and yolk sac. Accordingly, dorsal anterior finfold defects in edematous embryos and larvae (Figure 1H,I) represent a visible phenotypic anchor for ionoregulatory disruption, a third distinct oil-induced adverse outcome pathway.
Our analysis focused on key ionoregulatory proteins in MRCs and their associated genes, including Na+/K+ ATPase subunits (at1 genes, e.g. at1a1-a4, at1b1-b4), a urea transporter (ut1), a Na+/K+/2Cl- co-transporter (s12a2), the sodium-hydrogen exchanger Nhe3 (slc9a3), and a chloride channel, the latter an ortholog to the human cystic fibrosis transmembrane conductance regulator (cftr) (Hirose et al., 2003). The disruption of MRC function in oil-exposed haddock embryos corresponded to significantly lower levels of at1a1-3, at1b2-3, ut1, s12a2, sl9a3, and cftr (Figure 8, Supplementary file 1K). This downregulation primarily spanned a developmental window between 6 dpf and hatching (Figure 8, Supplementary file 1K). We also found significant transcriptional modifications of genes encoding other pumps, channels, and transporters specific to the nervous system and other tissues, including the aquaporins (aqp genes) that rapidly transport water across cell membranes (Supplementary file 1L). For example, there was a pronounced decrease in the expression of the primary neuronal water channel, aqp4, from 6 dpf to hatching and a strong upregulation of aqp12 at hatch (Figure 8). Crude oil exposures therefore appear to cause osmotic stress in the developing embryonic nervous system.
Figure 8
DEGs involved in osmoregulation.
E1–E6: Embryonic exposure, L1–L5: Larval exposure. Black letters: during exposure, blue letters: after exposure.
During larval exposure, edema accumulated in different compartments from embryos, and there were corresponding differences in expression of genes related to ion and water balance. At late stages of larval exposure, edema accumulated in the peritoneal cavity, and the dorsal marginal finfold appeared increased rather than decreased as in embryos. Fewer ion transport genes were affected, with increased expression observed for only at1a3 at 14 dph, and one at1b2 paralogue and at1b3 at 14 and 18 dph. Similarly, aquaporin genes were affected differently. Whereas aqp4 was unaltered in larvae, expression of aqp7 and aqp9 was increased while aqp3 and aqp12 were decreased (Figure 8, Supplementary file 1K).
A novel adverse outcome pathway: disruption of cholesterol homeostasis
Cholesterol is an essential structural component required for maintaining both the integrity and the fluidity of all metazoan cell membranes. Cholesterol is sourced from de novo cellular synthesis and from the uptake of external lipoprotein cholesterol from the circulation (Bjorkhem and Meaney, 2004). During fish development, cholesterol is mobilized from the yolk and distributed to cells during embryonic and larval yolk sac stages. Later, after the yolk is absorbed and larvae begin exogenous feeding, cholesterol is transported from the intestines. Crude-oil-exposed haddock embryos and larvae with the most severe morphological abnormalities were visibly unable to effectively mobilize yolk (Figure 1J,K). Moreover, larvae from the highest exposure concentration had less visible food in their stomachs relative to controls. These observations together suggest that chemical components of crude oil may deprive developing tissues of externally available cholesterol.
Of the 28 genes differentially regulated at all developmental time points, 5 are involved in cholesterol synthesis and feedback control (Table 1). These include 3-hydroxy-3-methylglutaryl-coenzyme A reductase, an enzyme encoded by hmdh that plays a primary feedback regulation role in the cholesterol biosynthetic pathway (Brown and Goldstein, 2009). Although reduced yolk absorption was not physically measureable in exposed embryos until after 3 dph, genes controlling cholesterol synthesis were upregulated much earlier, prior to visible cardiac circulation (6 dpf/E2). We also detected complex regulation of apolipoproteins during and after exposure, with mainly down-regulation of apob paralogs before first heartbeat (6 dpf/E2) and up-regulation of apoa4, apod apoeb and apoc2 after initiation of cardiac circulation (Supporting dataset 1, Sørhus et al., 2017). Scavenger receptor class B 1(encoded by scarb1), a transcytotic receptor for cholesterol-containing high-density lipoprotein (Acton et al., 1996), was also down-regulated in the exposed groups in haddock at 6 dpf (Supporting dataset 1, Sørhus et al., 2017). Pathway analysis was also consistent with a significant effect on cholesterol homeostasis (Figure 3—figure supplement 1).
In the larval exposure, we also detected increased expression of hdmh, erg1, cp51a, and npc2 (encoding the proteins squalene epoxidase, and cytochrome P450 51A, Niemann-Pick disease, type C2, respectively) at the latest stages examined (14 and 18 dph). Conversely, pathways involved in digestion – that is, pancreatic secretion, protein digestion, and protein absorption – were suppressed. This includes the downregulation of genes encoding digestive enzymes such as trypsin and chymotrypsin (Figure 4—figure supplement 1). The stomachs of oil-exposed larvae at the final time point were relatively empty, and the associated loss of food-derived nutrients likely triggered the observed increase in endogenous cellular cholesterol synthesis.
Unaltered gene expression in relation to visibly normal organs: lateral line and liver
Whereas abnormal phenotypic traits corresponded to differential gene expression, genes associated with normal traits were unchanged. For example, the lateral line and liver appeared normal in the most severely affected embryos (Figure 1—figure supplement 1). Consistent with this, markers for the lateral line (protein atonal homolog 1, atoh1) (Cai and Groves, 2015), liver growth (hepatocyte growth factor, met) and differentiation (genes encoding wnt2 and 2b protein (wnt2, wnt2b) (Wilkins and Pack, 2013), hematopoietically-expressed homeobox protein (hhex), and protein heg (heg)) (Supplementary file 1M) were not significantly modified. While some markers for liver differentiation, including genes encoding transferrin (tfr1) and fatty acid binding protein (fa10a) were differentially expressed, the changes were subtle and not consistent throughout development. Similarly, the related KEGG pathways that are inclusive of these genes were relatively unaffected by oil exposure at all time points. As noted above, IPA failed to identify significant enrichment of pathways related to phenotypically normal structures.
Discussion
Overall, we observed tight anchoring of temporal gene expression patterns to measurable phenotypes in crude oil-exposed haddock. First, the global changes in gene expression observed in the embryonic and larval exposures matched the general nature and severity of phenotypes. Embryonic exposure to crude oil or component cardiotoxic PAHs produces a coarse chemical phenocopy of the loss-of-function zebrafish mutants affecting heart function or development (Incardona et al., 2004). Many aspects of the oil toxicity phenotype are secondary to a loss of circulation—that is, defects in non-cardiac tissues, such as the eye, that require circulation for normal organogenesis (Incardona et al., 2004). In contrast, larval stages are primarily a period of rapid growth after major organogenesis is complete, and the injury phenotype in larvae is less severe. Consistent with this, embryos showed a larger number of DEGs than larvae, with a preponderance of down-regulation. Second, we identified specific changes in the expression of key genes involved in the function or morphogenesis of individual tissues and organs with visible abnormalities. Given unaltered gene expression and lack of statistically enriched pathways associated with apparently unaffected structures such as the liver, the DEGs in oil-exposed haddock indicate a disruption of specific developmental processes, as opposed to non-specific effects (e.g. general developmental delay).
This study demonstrates the ability to resolve changes in tissue-specific genes in a pool of total RNA from embryos and larvae, even for organs such as the heart that contribute a very small fraction to total tissue mass. A key finding is the general correlation of read count data with tissue specific patterns previously characterized in model species. Our findings have important implications for the utility of RNA-Seq and other quantitative measures of mRNA abundance in whole embryo or larval samples. For example, this demonstrates the feasibility of developing real-time monitoring tools based on quantification of gene expression in environmental samples collected following an oil spill. In addition, our extensive manual curation of the transcriptome groundtruths the utility of applications like IPA for use with non-model, non-mammalian organisms. Moreover, the changes in gene expression identified here represent significant information to be added to existing cardiotoxicity AOPs and novel AOPs associated with disruption of osmoregulation and lipid metabolism.
Two major initiating events for crude-oil-associated cardiac defects during fish development are chemical blockade of IKr repolarizing potassium currents, (encoded by kcnh2) and disruption of intracellular calcium handling, the latter culminating in sarcoplasmic reticulum (SR) calcium depletion through effects on either RyR or SERCA2 (encoded by ryr2 and at2a2, respectively) (Brette et al., 2014). In the fully formed heart, these pharmacologic effects impair cardiac function by inducing arrhythmia and reducing contractility (Incardona et al., 2009, 2014). However, rhythm and contractility defects during heart development lead to morphological defects (Andrés-Delgado and Mercader, 2016). In haddock embryos, these include poor chamber looping and outgrowth of the ventricle (Sørhus et al., 2016b). Our data demonstrate a transcriptional cascade that is tightly linked to these defects in cardiac function (cardiomyocyte intracellular calcium cycling) and form (heart chamber growth) through bmp10.
While chemical blockade of calcium cycling alone would be sufficient to induce the ventricular arrhythmias observed in oil-exposed embryolarval haddock, other elements of the E-C coupling physiological cascade were also selectively modified at the mRNA level. As shown previously using qPCR (Sørhus et al., 2016b), RNA-seq revealed a downregulation of genes encoding the Na/Ca exchanger (nac1) and IKr (kcnh2) in haddock embryos. Notably however, kcnh2 downregulation was only observed at later time points, in response to the highest oil exposures that caused the most severe phenotypes. There was no consistent decrease in the mRNA for a larger suite of proteins involved in cardiac E-C coupling. Assuming a broader transcriptional response in the heart was not masked by more abundant, normal expression of these genes in larger non-cardiac tissues, other non-specific mechanisms were unlikely to contribute to the formation of misshapen hearts. Moreover, the changing expression of key E-C coupling genes is a close match to the cardiac arrhythmia phenotype in both embryos and larvae. This includes a marked downregulation (>5 fold) of kcnj12, which encodes a subunit of the repolarizing IK1 current and causes the same types of ventricular arrhythmias as a reduction of kcnh2 (Domenighetti et al., 2007). At the same time, the up-regulation of E-C coupling genes following the chemical induction of arrhythmia in the larval exposures suggest that the more mature larval heart mounts a compensatory response.
Intracellular calcium has multiple direct roles in regulating gene expression, including the process of excitation-transcription (E-T) coupling (Wamhoff et al., 2006). Our findings suggest that E-T coupling may link reduced cardiomyocyte calcium cycling to structural defects in the haddock heart. Among vertebrates, bmp10 is expressed exclusively in the early tubular hearts of zebrafish, mouse, and chick embryos. The normal function of Bmp10 in the developing heart is primarily to drive ventricular cardiomyocyte proliferation during trabeculation (Grego-Bessa et al., 2007), a relatively late process during embryogenesis (around hatching stages in fish). Both loss of and excess Bmp10 leads to severe abnormalities in ventricular development in mouse (Chen et al., 2004). In mice lacking the RyR-associated Fkbp12 protein, disruption of SR calcium handling leads to ventricular defects through elevated bmp10 transcription (Shou et al., 1998; Chen et al., 2004), probably through calcium-dependent activation of myocardin (Wamhoff et al., 2004, 2006), the transcriptional activator of bmp10 (Huang et al., 2012). Moreover, Bmp10 is the most potent Bmp family member, showing greater resistance to Bmp antagonists (e.g. Noggin) than Bmp4 (Lichtner et al., 2013), the primary cardiac Bmp family member at early stages. While bmp10 normally functions at late stages of cardiogenesis, bmp4 is normally expressed at the cardiac cone stage in zebrafish. At this stage, bmp4 levels shift from radially symmetric to elevated on the left side of the cone, to drive proper looping (Chen et al., 2004). Loss of this asymmetry with ectopic bmp4 leads to un-looped hearts. Therefore, the premature up-regulation of a more potent family member, bmp10, at the cone stage is very likely to underlie the looping defects observed in oil-exposed embryos. Further evidence for bmp10 overexpression initiating abnormal cardiac morphogenesis is the secondary up-regulation of the Bmp10 target gene nkx25 (Chen et al., 2004). In zebrafish, nkx25 overexpression or loss of function (Tu et al., 2009) yields a reduced ventricle, and nkx25 must be down-regulated or antagonised in specific regions of the ventricle in order to form specialised conduction cells through the repressor action of Tbx3 (Hoogaars et al., 2004). The higher levels of tbx3 that follow upregulation of nkx25 and subsequent downregulation of anf thus likely reflect an imbalance between myocardial and non-myocardial cell fates within the ventricle. Thus, normal heart development in zebrafish requires tight control over bmp10, nkx25, and tbx3 expression, and all three genes were dysregulated in oil-exposed haddock. The observed ventricular and looping defects may represent chemical phenocopies of the fkbp12 mutant, wherein reduced intracellular calcium transients are linked to altered bmp10 expression by E-T coupling, thereby changing cell fate and chamber formation in the developing heart. Calcium-mediated E-T coupling may also be a feedback mechanism for altering the expression of genes that encode repolarizing potassium channels.
Although the haddock with craniofacial deformities superficially resemble certain zebrafish mutants, associated changes in gene expression suggest a more complex developmental perturbation than previously described. As is the case with Bmp10, Edn1 is a strong morphogen that must be tightly regulated, as both too much and too little lead to craniofacial defects (Sato et al., 2008; Clouthier et al., 2010). Higher levels of edn1 observed here are thus highly likely to be related to the craniofacial defects, which is supported by the subsequent up-regulation of its target dlx genes. However, the phenotype does not appear to reflect changes in dorsal-ventral patterning, as expected for perturbation of edn1-dependent NCC identity. Most studies of craniofacial development in zebrafish and other vertebrates have focused on NCC-derived cartilaginous precursors. However, craniofacial skeletal elements develop in synchrony with their associated muscles (Schilling and Kimmel, 1997), and defects in muscle formation or function produce phenotypes that are in many ways indistinguishable from primary cartilage defects and appear more similar to the phenotypes observed here (Hinits et al., 2011; Shwartz et al., 2012). Mesodermal precursors of craniofacial muscle cells express edn1, the downregulation of which is associated with terminal muscle differentiation (Choudhry et al., 2011). A failure to downregulate edn1 is consistent with failure to up-regulate or maintain myh gene expression. It is unknown whether a failure of skeletal muscle to terminally differentiate would lead to reduced expression of genes associated with NCC development and cartilage growth and survival—that is foxi1, Wnt pathway genes (wnt9b, fzd7a, mycn), Fgf receptors and tgfb3. On the other hand, the early expression of foxi1 has no clear linkage to craniofacial muscle development in the literature, but it is required for NCC survival indirectly by driving Fgf signaling (Edlund et al., 2014). Sorting out the precise pathways of craniofacial malformation will thus require a concerted spatial localization of these DEGs.
The regulation of genes controlling ion and water balance, combined with the collapse of the dorsal marginal finfold, provides a novel insight into the origins of edema in marine fish. Conversely, the genetic elements of this phenotype provide a starting point to study the molecular basis of buoyancy control in pelagic fish larvae. Our findings suggest that crude oil may disrupt MRC function, leading to decreased expression of MRC channel and transporter genes. These effects could be direct, indirect, or both. First, ion and water balance in embryos require cardiac circulation (Miyanishi et al., 2013). Flow is essential for osmoregulatory counter-current exchange in the gills and kidney (Somero, 1998; Grosell, 2011), and similar principals should apply to the yolk sac epidermis. Second, epidermal MRCs are likely directly exposed to the highest PAH concentrations, because the highest Cyp1a upregulation in response to oil occurs in the skin of embryos and yolk sac larvae (Sørhus et al., 2016b). PAHs or their metabolites could block solute carriers in a similar manner as shown for cardiac calcium and potassium channels (Brette et al., 2014; Incardona et al., 2014). Lastly, MRC function could be impaired if the metabolic cost of PAH degradation competes with ion transport. Although the effects of oil exposure on salt and water balance in fish embryos have not been examined previously, exposing water-soluble oil fraction to adult Pacific herring (Clupea pallasi) showed an effect on ion homeostasis (Kennedy and Farrell, 2005).
Genes involved in ion and water balance were differentially regulated in embryos and larvae, in close correspondence to the loss and expansion of the dorsal subdermal space, respectively. Edema flows along the path of least resistance and accumulates in expandable spaces. In haddock embryos, fluid moves from the dorsal subdermal space to the yolk sac. At the larval stage, the presence of fluid in the dorsal finfold suggests that a developed peritoneal cavity and body wall provide greater resistance than the thin, permeable yolk sac membrane. As a consequence, the central nervous system is likely deprived of water during embryonic development and turgidly stressed at the larval stages. Changes in cell volume can modify intra- and extracellular ionic concentrations, and thus the electrophysiological properties of neurons (Pasantes-Morales et al., 2000). This may underlie the observed down-regulation of the main water channel in the brain, aquaporin 4 (aqp4), during and after embryonic exposure. Although not a focus of the current study, there are likely important mechanistic connections between a loss of embryonic MRC function, disrupted osmoregulation in the brain, and pathophysiological impacts on neuronal development.
Our findings also provide new insights into fundamental transcriptional mechanisms of lipid mobilization from yolk. The delivery of yolk-derived cholesterol not been widely studied in early fish embryos, which are distinct from other vertebrates in that the yolk mass is contained separately from the vasculature by the yolk syncytial layer (YSL, or periblast in older literature). The YSL is a multicellular membrane that forms during gastrulation and has been studied mostly for its role in early pattern formation (Carvalho and Heisenberg, 2010). At later embryonic and larval stages, the YSL transports yolk-derived nutrients into the circulation (e.g, [Poupard et al., 2000]). Although cholesterol is exported from the yolk prior to established circulation (Fraher et al., 2016), the cellular basis for this is not known—for example, by trancytosis or direct membrane transport. Cellular sterol levels are tightly controlled by membrane-bound transcription factors, which are cleaved and activated when membrane cholesterol levels fall, leading to transcription of hmdh, the rate-limiting enzyme for cholesterol synthesis (Brown and Goldstein, 2009). Therefore, the brisk up-regulation of intrinsic cholesterol biosynthetic genes, especially hmdh, prior to a visible mobilization of yolk, indicates the importance of yolk-derived cholesterol for embryonic tissues. In zebrafish and other fish, apolipoprotein genes (e.g. ApoB and ApoE) as well as Scarb1 are first expressed in YSL (Poupard et al., 2000; Thisse et al., 2001; Otis et al., 2015). The downregulation of multiple apob paralogs and scarb1 during early development suggest a specific defect in packaging and transporting of lipoprotein-cholesterol in the YSL, possibly involving Scarb1-dependent transcytosis. At stages subsequent to heartbeat initiation, the upregulation of intrinsic cholesterol biosynthetic genes likely reflects cholesterol deprivation as a consequence of the heart failing to deliver lipoproteins from the yolk (embryos) and intestine (larvae). The identification of broader lipid metabolism pathways by IPA also suggests that oil exposure leads to more global derangements relating to either poor yolk absorption or dysfunction of the YSL. It is well known both that embryonic oil exposure leads to later growth impairment (e.g. [Heintz, 2000; Incardona and Scholz, 2016]) and that lipids provide critical fuel for marine fish larvae, particularly at the first-feeding stage (Tocher et al., 2003). The consequences of disordered lipid metabolism for larval physiology and survival should be a focus for future studies. Importantly, the induction of cholesterol synthetic genes is a promising and novel indicator of crude oil toxicity to fish embryos.
Finally, we consider the four individual genes that were consistently upregulated across all developmental stages: ccbe1, rhag, foxq1, and fgf7. Elevated tissue pressure is a signal for lymphangiogenesis (Schulte-Merker et al., 2011) and the marked increase in ccbe1 expression is consistent with a compensatory formation of lymphatics to alleviate the physical effects of edema (Planas-Paz et al., 2012). In zebrafish, the secreted Ccbe1 protein appears to function exclusively in lymphangiogenesis (Hogan et al., 2009; Le Guen et al., 2014) by enhancing the activation of VEGF-C (Le Guen et al., 2014). For future spills, quantitative measures of ccbe1 upregulation should serve as very sensitive indicators of edema formation in crude oil-exposed fish embryos.
The Rh proteins are primarily structural components of erythrocyte membranes but were also recently identified as ammonia transporters (Weiner and Verlander, 2014). The rhag, rhbg, and rhcg genes are important for excretion of nitrogenous waste in fish (Braun et al., 2009). The strong upregulation of rhag observed in haddock embryos and larvae might be a consequence of MRC dysfunction. The increase in rhag expression corresponded to a downregulation of urea transporters ut1 and ut2 in embryos but not larvae, suggesting metabolic defects relating to amino acid or protein catabolism and an increased demand for nitrogen excretion. Alternatively, rhag may play a novel role in embryolarval physiology.
The last two highly responsive genes, foxq1 and fgf7, may be linked to the craniofacial injury phenotype based on prior work in other species. In zebrafish and frog embryos, foxq1 is expressed in the pharyngeal region (Choi et al., 2006; Planchart and Mattingly, 2010). In chick embryos, fgf7 is first expressed in the pharyngeal endoderm and head mesoderm (Kumar and Chapman, 2012). However, foxq1 has been shown to be a downstream target of AhR in other tissues (Faust et al., 2013), and crude oil exposures strongly induced fgf7 expression in the livers of juvenile polar cod (Andersen et al., 2015). Both genes are promising new biomarkers for future studies, but tissue localization during craniofacial development in embryonic fish is needed to confirm a role in this adverse outcome pathway. In zebrafish exposed to the dioxin TCDD, foxq1 was up-regulated in the branchial arches (Planchart and Mattingly, 2010), but TCDD-induced craniofacial defects have been shown to be entirely secondary to cardiotoxicity (Lanham et al., 2014).
Lastly, our findings show how transcriptomics can inform chemical genetics and environmental forensics in non-model organisms. First, we used known spatial mRNA distributions in model species (primarily zebrafish) to more accurately phenotypically anchor the transcriptome data for crude oil-exposed haddock. This accelerates the pace of discovery, particularly given difficulties in applying zebrafish methods for in situ hybridization to wild species. Second, our results reveal transcriptional aspects of chemical injury in fish, in response to a global pollution threat. Although endocrine disrupting compounds act on steroid hormone receptors or epigenetically modify DNA or histones (Walker and Gore, 2011), it has been much less clear whether crude oil, which acts on protein targets (e.g. [Brette et al., 2014]), also influences mRNA levels as part of a widely studied developmental injury phenotype.
In conclusion, our findings greatly expand our understanding of crude oil impacts to fish early life stages at a molecular level. The scientific focus on highly vulnerable fish embryos and larvae began with the 1989 Exxon Valdez disaster in Prince William Sound, Alaska. Major crude oil spills have continued worldwide—for example, the 2002 Prestige spill in Spain, the 2007 Hebei Spirit spill in Korea, and the 2010 Deepwater Horizon blowout in the northern Gulf of Mexico. An enduring challenge spans these and other spills: namely, the accurate determination of fisheries losses as a consequence of oiled spawning habitats. Our results more clearly define the different categories of developmental injury. We have also identified responsive genes that hold promise as sensitive molecular indicators of physiological stress. These can be developed into novel assessment tools for diagnostic use in future spill zones. Lastly, differential patterns of gene expression in oil-exposed haddock provide insight into fundamental but still poorly understood developmental processes in marine fish, including calcium-mediated excitation-transcription coupling, fluid balance, lipid mobilization, and buoyancy.
Materials and methods
For detailed procedures on animal collection, maintenance, exposure regime and analytical chemistry see (Sørhus et al., 2015, 2016b). All animal experiments within the study were approved by NARA, the governmental Norwegian Animal Research Authority (http://www.fdu.no/fdu/, reference number 2012/275334-2). All methods were performed in accordance with approved guidelines.
Animal collection, maintenance, and exposure set-up and regime
Request a detailed protocolThe samples transcriptome profiled here are the exact same samples described in Sørhus et al. (2016b), with the methodology for animal collection, maintenance, and crude oil exposure provided therein. Briefly, fertilized eggs were collected from tanks with wild captured Atlantic haddock (collected spring 2013) and transferred to indoor egg incubators (7°C). At experiment start, eggs from the egg incubators were transferred to 50 L exposure tanks (7.8°C). Embryonic and larval stages of Atlantic haddock were exposed to two doses, low and high dose, in addition to a pulsed high dose (see Figures 3 and 4 for details). Heidrun oil blend was artificially weathered by distillation and then pumped into the dispersion system using a HPLC pump (Shimadzu, LC-20AD Liquid Chromatograph Pump) with a flow of 5 μL/min together with a flow of seawater of 180 mL/min. The system described in Nordtug et. al 2011 (Nordtug et al., 2011) generates an oil dispersion with oil droplets in the low μm range with a nominal oil load of 26 mg/L (stock solution). Water samples were collected from each exposure tank at the beginning of each experiment and subjected for detailed PAH analysis. At end of exposure, pooled samples of eggs and larvae were extracted by solid liquid extraction and purified by solid phase extraction prior to analysis by GC-MS/MS (Sørensen et al., 2015) to reveal the PAH content in animals. After end of exposure, the animals were transferred to new tanks with clean water. Images of the animals were collected from 12 and 8 stages during and after embryonic and larval exposure, respectively. Videos from 60 (embryonic) and 48 (larval) individual embryos/larvae per treatment per stage were collected from 9 dpf, 0 dph, and 3 dph (embryonic exposure) and 2 dph and 9 dph (larval exposure). The experiments included four replicate tanks for each dose, and pooled samples from three replicate tanks for each dose were subjected for sequencing (see details below).
Total RNA and cDNA preparation
Request a detailed protocolTotal RNA was isolated from frozen pools of embryos and larvae using Trizol reagent (Invitrogen) per manufacturer instructions. This included a DNase treatment step using a TURBO DNA-free kit (Life Technologies Corporation). RNA was quantified using a Nanodrop spectrophotometer (NanoDrop Technologies), and confirmed using a 2100 Bioanalyzer (Agilent Technologies). cDNA was subsequently generated using SuperScript VILO cDNA Synthesis Kit (Life Technologies Corporation), according to the manufacturer’s instructions. The cDNA was normalized to obtain a concentration of 50 ng/µL.
Real-time qPCR
Request a detailed protocolSix responsive genes from the transcriptome were validated by real-time quantitative PCR (qPCR) (Figure 3—figure supplement 2 [Embryonic exposure], Figure 4—figure supplement 2 [Larval exposure]). Specific primers and probes (Supplementary file 1N) for a reference gene (ef1α) and the six DEGs (cp1a, wnt11, kcnh2, nac1, cac1c, and at2a2) were designed with either Primer Express Software (Applied Biosystems) or Eurofins qPCR probe and primer design software (Eurofins Scientific), according to the manufacturer’s guidelines. The two methods generally yielded the same quantitative trends. Primer and probe sequences are shown in Supplementary file 1J. TaqMan PCR assays were performed in duplicate, using 384-well optical plates on an ABI Prism Fast 7900HT Sequence Detection System (Applied Biosystems) with settings as follows: 50°C for 2 min, 95°C for 20 s, followed by a 40 cycles of 95°C for 1 s and 60°C for 20 s. For each 10 μL PCR reaction, a 2 μL cDNA 1:40 dilution was mixed with 200 nM fluorogenic probe, 900 nM sense primer, and 900 nM antisense primer in 1xTaqMan Fast Advanced Master Mix (Applied Biosystems). Gene expression was calculated relative to the exposure time zero sample (2 dpf and 0 dph in embryonic and larval exposures, respectively) using the ΔΔΔCt method, generating reference residuals (Edmunds et al., 2014) from ef1a and at2a2.
Extraction of mRNA, RNA sequencing, and bioinformatics
Request a detailed protocolcDNA library preparation and sequencing was performed by the Norwegian Sequencing Centre (NSC, Oslo, Norway) using the Illumina TruSeq RNA Sample Preparation Kit. A total of 132 samples were sequenced and 126 were subjected for analysis: Three (control, low and high dose) or two (pulse) biological replicates for each dose from six stages during and after embryonic exposure and three biological replicates for each dose from five stages during larval exposure (see Figures 3 and 4). Paired-end libraries were sequenced on the Illumina HiSeq2500 platform. The raw data are available from the Sequence Read Archive (SRA) at NCBI (Accession ID: PRJNA328092).
The high sequence similarity between the two species justified use of the cod template, and we chose a verified model over a reference-free de novo transcriptome approach to avoid fragmentation noise and false positives from un-collapsed genes. The average sequence similarity between mapped haddock reads and the previously verified cod genome was 98.4%. Moreover, of the 20,954 annotated cod genes, there were 18,990 (90.6%) corresponding haddock genes with at least 10 reads in one sample. Thus, the RNA sequencing data were mapped to the coding sequences of the cod genes (Star et al., 2011) using the Bowtie aligner (RRID:SCR_005476) (Langmead et al., 2009) and annotated as described in Sørhus et al. (2016a). Samtools idxstat (RRID:SCR_002105) (Li et al., 2009) was used to extract the number of mapped reads which were then normalised to the total number of mapped sequences. NOISeqBIO (RRID:SCR_003002) (Tarazona et al., 2011, 2015) was used to identify differentially expressed genes (DEGs, threshold of 0.95). The total number of reads averaged 41.39 million per sample and the mapping efficiency averaged 32.69%, giving an overall average of 13.51 million mapped paired end reads (125 bp) for each sample. Given the absence of untranslated regions (UTRs) and mitochondrial genes from the reference cod genome, reads with a UTR sequence and all reads for mitochondrial genes were excluded from the haddock analyses.
Only genes with 10 reads or more in at least one of the samples were included for further analysis. Kyoto Encyclopedia of Genes and Genomes (KEGG) pathways analysis (RRID:SCR_012773) (Kanehisa et al., 2012) was performed by mapping the KEGG annotated DEGs from NOISeqBIO to KEGG pathways as described in the KEGG Mapper tool. Heat maps were generated from fold change data in R (R Core Team , 2013) and Venn diagrams were created using the web-tool, Venn (http://bioinformatics.psb.ugent.be/webtools/Venn). Individual genes involved in cardiac and craniofacial development, osmoregulation and lipid metabolism, or not directly linked to KEGG pathways were curated manually in an intensive process that took a full year. In a previous effort, we characterized the transcriptome of normally developing haddock embryos and larvae (Sørhus et al., 2016a). Through extensive literature searches (PubMed, Web of Science, and Google Scholar) and reading, lists of roughly 150 key genes involved in cardiac and craniofacial development and cardiac function were assembled. After obtaining the oil exposure RNA-Seq dataset, these lists were expanded by further literature searches. Genes of interest were identified as ones that showed strong dose-dependency and had potential phenotypic association based on the literature. We relied heavily on the expression database at the Zebrafish Information Network (www.zfin.org) to determine whether individual genes of interest were expressed in the relevant tissues at the appropriate developmental stage to be associated with phenotypes, helping to narrow in on key linkages. This was performed both by searching individual gene names and by generating lists of all genes expressed in a specific tissue at a time relevant to the phenotypes (e.g. all genes expressed in the ethmoid plate between segmentation and hatch). These lists were cross-referenced to the list of DEGs to identify candidates. If zfin.org lacked expression data, we searched for any papers describing tissue localization by in situ hybridization in other model or non-model species.
Qiagen’s Ingenuity Pathway Analysis (IPA version 01–04) (RRID:SCR_008653) was used subsequent to our manual curation. The fold change and p values were extracted from the original NOISeqBIO output. Fold change values were multiplied by −1 to flip the direction (opposite the convention used by IPA), and any fold change values between −1 and 1 was set to 1, as non-significant values and those below 10 reads were originally set to 0.5 and 0.75. The cod genes were then BLASTed (version 2.3.0+) against the ENSEMBL zebrafish (GRCz10) and human (GCRh38) databases. The top match for each gene, filtered based on e-value (cutoff 10−5), was used to build a mapping table of the genes to the IPA database. For genes that lacked mapping at that point, SwissProt information was used to manually map. The mapping information was then combined with the fold change and p-value data and uploaded to IPA. This resulted in 17608 of 20954 genes (84%) mapping to the IPA database. Data were then analyzed using a fold change cutoff of 2.0 for both up- and down-regulation, considering relationships only when confidence was equal to the experimentally observed level. For this analysis, IPA uses a built-in right-tailed Fisher Exact Test to calculate p-values for significant functions and pathways. After each time point was analyzed individually in this manner, data outputs were used in comparison analyses in which the control dataset was compared to each dose at each time point. The output of the comparison analyses included rankings of top pathways in a variety of categories and subcategories. These rankings were expanded in the IPA software interface and inspected manually to generate Table 3 and Supplementary file 1B and 1C.
Supporting dataset 1: Expression of all genes in control, low dose, pulse dose, high dose during and after embryonic exposure. Highlighted FC: Prob ≥ 0.95; Not highlighted FC: Prob 0.8–0.95; FC = 0.5: Not significant; FC = 0.75: Less than 10 reads in both treatment and control. SP, swissprot; GB, genebank; E1, 3 dpf; E2, 6 dpf; E3, 10 dpf; E4, 11 dpf; E5, 0 dph; E6, 3 dph; C, control; L, low dose; P, pulse dose; H, high dose; FC; fold change; prob; probability BioNoiseq. (Sørhus et al., 2017).
Supporting dataset 2: Expression of all genes in control, low dose, pulse dose, high dose during larval exposure. Highlighted FC: Prob ≥ 0.95; Not highlighted FC: Prob 0.8–0.95; FC = 0.5: Not significant; FC = 0.75: Less than 10 reads in both treatment and control. SP, swissprot; GB, genebank; L1, 1 dph; L2, 3 dph; L3, 9 dph; L4, 14 dph; L5, 18 dph; C, control; L, low dose; P, pulse dose; H, high dose; FC; fold change; prob; probability BioNoiseq. (Sørhus et al., 2017).
Data availability
-
Data from: Novel adverse outcome pathways revealed by chemical genetics in a developing marine fishAvailable at Dryad Digital Repository under a CC0 Public Domain Dedication.
-
Sequence dataPublicly available at the NCBI Sequence Read Archive (accession no: PRJNA328092).
References
-
Identification of compounds in heavy fuel oil that are chronically toxic to rainbow trout embryos by effects-driven chemical fractionationEnvironmental Toxicology and Chemistry 33:825–835.https://doi.org/10.1002/etc.2497
-
Interplay between cardiac function and heart developmentBiochimica et Biophysica Acta (BBA) - Molecular Cell Research 1863:1707–1716.https://doi.org/10.1016/j.bbamcr.2016.03.004
-
Adverse outcome pathways: a conceptual framework to support ecotoxicology research and risk assessmentEnvironmental Toxicology and Chemistry 29:730–741.https://doi.org/10.1002/etc.34
-
Brain cholesterol: Long secret life behind a barrierArteriosclerosis, Thrombosis, and Vascular Biology 24:806–815.https://doi.org/10.1161/01.ATV.0000120374.59826.1b
-
Cholesterol feedback: from Schoenheimer's bottle to Scap's MELADLJournal of Lipid Research 50 Suppl:S15–S27.https://doi.org/10.1194/jlr.R800054-JLR200
-
The role of atonal factors in mechanosensory cell specification and functionMolecular Neurobiology 52:1315–1329.https://doi.org/10.1007/s12035-014-8925-0
-
A perspective on the toxicity of petrogenic PAHs to developing fish embryos related to environmental chemistryHuman and Ecological Risk Assessment: An International Journal 15:1084–1098.https://doi.org/10.1080/10807030903304708
-
The yolk syncytial layer in early zebrafish developmentTrends in Cell Biology 20:586–592.https://doi.org/10.1016/j.tcb.2010.06.009
-
Overexpression of bone morphogenetic protein 10 in myocardium disrupts cardiac postnatal hypertrophic growthJournal of Biological Chemistry 281:27481–27491.https://doi.org/10.1074/jbc.M604818200
-
Developmental expression of FoxJ1.2, FoxJ2, and FoxQ1 in Xenopus tropicalisGene Expression Patterns 6:443–447.https://doi.org/10.1016/j.modgep.2005.11.007
-
Alcama mediates Edn1 signaling during zebrafish cartilage morphogenesisDevelopmental Biology 349:483–493.https://doi.org/10.1016/j.ydbio.2010.11.006
-
Regulation of facial morphogenesis by endothelin signaling: insights from mice and fishAmerican Journal of Medical Genetics Part A 152A:2962–2973.https://doi.org/10.1002/ajmg.a.33568
-
Chronic angiotensin II stimulation in the heart produces an acquired long QT syndrome associated with IK1 potassium current downregulationJournal of Molecular and Cellular Cardiology 42:63–70.https://doi.org/10.1016/j.yjmcc.2006.09.019
-
Toward enhanced MIQE compliance: reference residual normalization of qPCR gene expression dataJournal of Biomolecular Techniques 25:54–60.https://doi.org/10.7171/jbt.14-2502-003
-
AhR-mediated changes in global gene expression in rat liver progenitor cellsArchives of Toxicology 87:681–698.https://doi.org/10.1007/s00204-012-0979-z
-
Rit Fiskeideildar1–68, Embryonic development of five species of gadoid fishes in icelandic waters, Rit Fiskeideildar, Marine Research Institute.
-
Are adverse outcome pathways here to stay?Environmental Science & Technology 49:3–9.https://doi.org/10.1021/es504976d
-
Notch signaling is essential for ventricular chamber developmentDevelopmental Cell 12:415–429.https://doi.org/10.1016/j.devcel.2006.12.011
-
Stages of embryonic development in the Atlantic cod Gadus morhuaJournal of Morphology 259:255–270.https://doi.org/10.1002/jmor.10222
-
Molecular biology of major components of chloride cellsComparative Biochemistry and Physiology Part B: Biochemistry and Molecular Biology 136:593–620.https://doi.org/10.1016/S1096-4959(03)00287-2
-
Ccbe1 is required for embryonic lymphangiogenesis and venous sproutingNature Genetics 41:396–398.https://doi.org/10.1038/ng.321
-
Myocardin regulates BMP10 expression and is required for heart developmentJournal of Clinical Investigation 122:3678–3691.https://doi.org/10.1172/JCI63635
-
Cardiac arrhythmia is the primary response of embryonic pacific herring (Clupea pallasi) exposed to crude oil during weatheringEnvironmental Science & Technology 43:201–207.https://doi.org/10.1021/es802270t
-
Aryl hydrocarbon receptor-independent toxicity of weathered crude oil during fish developmentEnvironmental Health Perspectives 113:1755–1762.https://doi.org/10.1289/ehp.8230
-
Defects in cardiac function precede morphological abnormalities in fish embryos exposed to polycyclic aromatic hydrocarbonsToxicology and Applied Pharmacology 196:191–205.https://doi.org/10.1016/j.taap.2003.11.026
-
KEGG for integration and interpretation of large-scale molecular data setsNucleic Acids Research 40:D109–D114.https://doi.org/10.1093/nar/gkr988
-
Ion homeostasis and interrenal stress responses in juvenile Pacific Herring, Clupea pallasi, exposed to the water-soluble fraction of crude oilJournal of Experimental Marine Biology and Ecology 323:43–56.https://doi.org/10.1016/j.jembe.2005.02.021
-
Specification and morphogenesis of the zebrafish larval head skeletonDevelopmental Biology 233:239–257.https://doi.org/10.1006/dbio.2001.0201
-
Cranial neural crest and development of the head skeletonAdvances in Experimental Medicine and Biology 589:120–133.https://doi.org/10.1007/978-0-387-46954-6_7
-
Adverse outcome pathways and ecological risk assessment: bridging to population-level effectsEnvironmental Toxicology and Chemistry 30:64–76.https://doi.org/10.1002/etc.375
-
Cloning and expression analysis of Fgf5, 6 and 7 during early chick developmentGene Expression Patterns 12:245–253.https://doi.org/10.1016/j.gep.2012.05.002
-
Cardiac myocyte-specific AHR activation phenocopies TCDD-induced toxicity in zebrafishToxicological Sciences 141:141–154.https://doi.org/10.1093/toxsci/kfu111
-
The Sequence Alignment/Map format and SAMtoolsBioinformatics 25:2078–2079.https://doi.org/10.1093/bioinformatics/btp352
-
MyoD and Myf6 gene expression patterns in skeletal muscle during embryonic and posthatch development in the domestic duck (Anas platyrhynchos domestica)Journal of Animal Breeding and Genetics 131:194–201.https://doi.org/10.1111/jbg.12057
-
Role of cardiac natriuretic peptides in seawater adaptation of medaka embryos as revealed by loss-of-function analysisAJP: Regulatory, Integrative and Comparative Physiology 304:R423–R434.https://doi.org/10.1152/ajpregu.00384.2012
-
Histology of the Atlantic cod, Gadus morhua: an atlasPart 4. Eleuthereoembryo and larva, Histology of the Atlantic cod, Gadus morhua: an atlas, In, Canadian Special Publication of Fisheries and Aquatic Sciences.
-
Method for generating parameterized ecotoxicity data of dispersed oil for use in environmental modellingMarine Pollution Bulletin 62:2106–2113.https://doi.org/10.1016/j.marpolbul.2011.07.015
-
Signaling events during swelling and regulatory volume decreaseNeurochemical Research 25:1301–1314.https://doi.org/10.1023/A:1007652330703
-
Myosin heavy chain expression in cranial, pectoral fin, and tail muscle regions of zebrafish embryosMolecular Reproduction and Development 63:422–429.https://doi.org/10.1002/mrd.10201
-
Mechanoinduction of lymph vessel expansionThe EMBO Journal 31:788–804.https://doi.org/10.1038/emboj.2011.456
-
2,3,7,8-Tetrachlorodibenzo-p-dioxin upregulates FoxQ1b in zebrafish jaw primordiumChemical Research in Toxicology 23:480–487.https://doi.org/10.1021/tx9003165
-
BookR: A language and environment for statistical computingVienna, Austria: R Foundation for Statistical Computing.
-
Musculoskeletal patterning in the pharyngeal segments of the zebrafish embryoDevelopment 124:2945–2960.
-
Lymphatic vascular morphogenesis in development, physiology, and diseaseThe Journal of Cell Biology 193:607–618.https://doi.org/10.1083/jcb.201012094
-
The effect of water conservation on the structure of marine fish embryos and larvaeJournal of the Marine Biological Association of the United Kingdom 35:275.https://doi.org/10.1017/S0025315400009103
-
FKBP12-deficient mice exhibit congenital dilated cardiomyopathy and altered ryanodine receptor functionBiophysical Journal 74:A353.
-
Cold Ocean Physiology33–57, Adaptation to cold and depth: contrasts between polar and deep-sea animals, Cold Ocean Physiology, Cambridge University Press.
-
Optimization and comparison of miniaturized extraction techniques for PAHs from crude oil exposed Atlantic cod and haddock eggsAnalytical and Bioanalytical Chemistry 408:1023–1032.https://doi.org/10.1007/s00216-015-9225-x
-
Differential expression in RNA-seq: a matter of depthGenome Research 21:2213–2223.https://doi.org/10.1101/gr.124321.111
-
Metabolism and functions of lipids and fatty acids in teleost fishReviews in Fisheries Science 11:107–184.https://doi.org/10.1080/713610925
-
Adverse outcome pathway development II: best practicesToxicological Sciences 142:321–330.https://doi.org/10.1093/toxsci/kfu200
-
Transgenerational neuroendocrine disruption of reproductionNature Reviews Endocrinology 7:197–207.https://doi.org/10.1038/nrendo.2010.215
-
Excitation-transcription coupling in arterial smooth muscleCirculation Research 98:868–878.https://doi.org/10.1161/01.RES.0000216596.73005.3c
-
Ammonia transport in the kidney by Rhesus glycoproteinsAJP: Renal Physiology 306:F1107–F1120.https://doi.org/10.1152/ajprenal.00013.2014
-
Zebrafish models of human liver development and diseaseComprehensive Physiology 3:1213–1230.https://doi.org/10.1002/cphy.c120021
Decision letter
-
Marianne BronnerReviewing Editor; California Institute of Technology, United States
In the interests of transparency, eLife includes the editorial decision letter and accompanying author responses. A lightly edited version of the letter sent to the authors after peer review is shown, indicating the most substantive concerns; minor comments are not usually included.
Thank you for submitting your article "Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Marianne Bronner as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers have extensively discussed the reviews with one another and the Reviewing Editor has drafted this decision to help you prepare a revised submission.
As you will see from the individual reviews below, all of the reviewers had extensive issues with the analysis and lack of detail in the description of method. Even more importantly, however, they also felt that the there was a great deal of over-interpretation and "overselling of the work", rendering many statements unsupported by the data. We feel that this is fixable but will require considerable rewriting and re-review. If you feel you are able to revise the paper accordingly, we would be willing to consider a significantly revised version for re-review.
Reviewer #1:
In this article the authors expand on previous experiments, where developing haddock are exposed to low levels of crude oil, to test for impacts on sensitive developmental endpoints. Developmental impacts are described in detail in previous reports, and this study reports transcriptomic changes during development following three oil exposure regimes.
I have two main difficulties with this paper. First, there is insufficient description of methods for gene expression analyses. Second, there is much speculation and unsupported statements such that much of the manuscript reads as over-statement.
Gene expression analysis:
Though the experimental design includes three different exposure regimes, specific gene transcriptional responses are not correlated in any way to dosing regime (high, low, pulse). The authors ought to show gene expression dose responses, at least for the genes for which they are trying to build cause-effect relationships with adverse outcomes. More convincing arguments for cause-effect would be supported by data showing that dysregulation of the gene precedes the appearance of the phenotype, and that dysregulation is more pronounced in treatments with more pronounced phenotypic effects (e.g., higher doses). The lack of detail reported here could make a reader concerned that perhaps interpretation of gene-specific responses are not as clear-cut as represented in the text. Are Figures 4, 5, 6, 7 showing data from any dose, or just high dose? Furthermore, there is no description of how emergence of phenotypes varied with dosing regime. Text (starting paragraph two of the Results section) that describes developmental phenotypes does not distinguish between effects that differ between doses, and all of the images from exposed animals in Figure 1 appear to be only from the high dose.
Other aspects of gene expression analysis that are problematic:
There is no description of the statistical model. What was compared to what? How many tests were there? Was there multiple test correction given the size of the gene set (e.g., false discovery rate correction)? Was this experiment analyzed as a fully specified statistical model, or were individual treatments compared to individual control treatments separately from others, and if so was there multiple test correction for that? There is insufficient information provided here to be able to evaluate the quality and robustness of the gene expression analysis. Given that gene expression analysis is a focus of the paper this is a problem.
41.9 million reads (subsection “Extraction of mRNA, RNA sequencing and bioinformatics”) – is this per sample?
32.69% mapping efficiency (subsection “Extraction of mRNA, RNA sequencing and bioinformatics”) – this is extremely low! Why is this so low if sequence similarity between haddock and reference (cod) genome is so high? How was successful mapping assessed – how were mapping parameters set? This should be a red flag that something went wrong in the mapping, or that there was significant contamination issues.
After mapping, this left 31.51 million reads "for each group". I'm not sure what this means "for each group". It is generally accepted that a minimum of ~10M reads PER SAMPLE (that is, per experimental unit) is adequate for RNA-seq.
"three biological replicates per stage” – section "Total RNA and cDNA preparation" states that RNA was from pools of embryos. So are these more accurately 3 replicate pools?
KEGG analysis: what were the criteria used to include a KEGG pathway in figures S4/S5? At least one of the criteria should be statistically significant enrichment. From my look at these figures, there are no obvious pathway-level connections to the phenotypes/physiological pathways on which the authors hang their hat. E.g., these figures do not have any terms that obviously relate to cholesterol homeostasis. Nor does Figure 3—figure supplement 1 indicate any enrichment of pathways having to do with ionoregulation, or craniofacial development. So what is the analysis here?
Selected gene sets: The authors make assertions about physiological mechanisms given responses of specific genes with inferred functions that make a nice story. However it is unclear what is the null expectation from this type of analysis.
E.g. Subsection “Genes associated with defects in cardiac function and morphogenesis”: "We therefore focused on genes associated with cardiomyocyte membrane potential and intracellular Ca2+ cycling" – how were these genes selected?
In the same section: "The expression of a few E-C coupling genes" – was there an unbiased collection of E-C coupling genes, then test for significant enrichment of these genes?
In the same section: "genes involved in cardiac morphogenesis" – again, how were these collected/curated as a gene set? Please provide the query gene set.
Subsection “Genes associated with craniofacial abnormalities”: "We therefore interpreted developmental changes in haddock gene expression in the context of these well-characterised zebrafish mutants." – so the "tester" set here was all genes shown in zebrafish mutants to be involved in neural crest cell or muscle development? Please include this set of genes. The 12 genes indicated below this section – is this a significant enrichment?
Subsection “Genes associated with ion and water regulatory imbalance”: "Our analysis focused on key ionoregulatory proteins in MRCs and their associated genes, including Na+/K+ ATPase subunits (at1 genes, e.g. at1a1-a4, at1b1-b4), a urea transporter (ut1), a Na+/K+/2Cl- co-transporter (s12a2), the sodium hydrogen exchanger Nhe3 (slc9a3), and a chloride channel" – this is clearly not a complete or objective curation of ionoregulatory genes.
Is there an unbiased query gene set for any of these analyses? This whole "genes associated with phenotypes" section – is there objective analysis? What makes this anything more than just story telling given some cherry-picked genes? The authors should at least use language that proposes some hypotheses, rather than writing text so as to imply strong conclusions of cause-effect adverse outcome pathways. This is a problem that permeates the reporting of results and discussion, and contributes in a large way to overstatement of the results that are expanded on below.
Unsupported claims and overstatement issues:
There are pervasive issues throughout the manuscript, often associated with claims or conclusions that are (or should be) in fact stated as hypotheses. It is often unclear whether certain statements are supported by their data, supported by the literature, or are speculation or informed hypotheses. Some examples follow:
Subsection “Genes associated with defects in cardiac function and morphogenesis”: "Overall, the transcriptional response to disrupted E-C coupling was not a simple pattern of compensation" – how do they know that this is a transcriptional response to E-C coupling? Is expression not being measured in whole animals? What proportion of overall tissue mass is accounted for by the heart? I think the authors are drawing conclusions where they ought to be proposing hypotheses. This is an issue in MANY places throughout the manuscript
Subsection “Genes associated with craniofacial abnormalities”: "Overall, the genes identified above for both neural crest and muscle lineages represent a more complex pattern of dysregulation than previously been reported for any of the individual zebrafish craniofacial mutants. This suggests that crude oil is acting on novel developmental processes." – It is unclear what is supporting this assertion (1st sentence), and how the assertion leads to the conclusion (2nd sentence)
Subsection “Genes associated with ion and water regulatory imbalance”: "and they should therefore lose water along a diffusion gradient if osmoregulation is disrupted as a consequence of heart and circulatory failure." – Does oil exposure cause ionoregulatory dysfunction and/or water loss in developing marine fish? Are these statements speculation/hypotheses, or are they supported in the literature? There are several prominent studies on the effects of oil exposure on osmoregulation in fish, but none of these are referred to.
In the same section: "a third distinct oil-induced adverse outcome pathway." – how do they know that this is a pathway distinct from cardiac dysfunction? E.g., rather than just a secondary manifestation of circulatory failure?
"Crude oil exposures therefore appear to cause osmotic stress in the developing embryonic nervous system" – the authors are implying an overall impact on ionoregulatory abilities? Is there any support for this here or in the literature? E.g., any studies showing dysfunction of net sodium or chloride flux upon oil exposure? Do the authors mean to propose this as a hypothesis?
Subsection: "A novel adverse outcome pathway: disruption of cholesterol homeostasis" – this subheading, and the following paragraphs, claim novelty, and appear to claim precedent for this discovery. Has oil exposure impacts on yolk utilization not previously been reported in the literature? I can find a publication from the 1990s after just a couple minutes searching. Furthermore, the studies presented here do not provide direct evidence for impacts on yolk mobilization or cholesterol homeostasis. E.g., how do the authors know that the apparently increased yolk size in Figure 1K isn't just a result of fluid accumulation stemming from observed edema?
In the same section: "Pathway analysis was also consistent with a significant effect on cholesterol homeostasis” – The associated figure is is missing. Also, staring at Figure 3—figure supplement 1, I see no pathways that include any terms obviously related to cholesterol homeostasis. Nor for that matter does Figure 3—figure supplement 1 indicate any enrichment of pathways having to do with ionoregulation, or craniofacial development,
Subsection “Unaltered gene expression in relation to visibly normal organs: lateral line and liver.”: "Similarly, the related KEGG pathways that are inclusive of these genes were relatively unaffected by oil exposure at all time points" – but what does this really mean? KEGG pathways (Figure 3—figure supplements 1 and 2), as far as I can tell, don't implicate many (any?) of the mechanisms that the authors are building their story on – e.g., ionoregulation, cholesterol homeostasis, craniofacial development. If this paragraph was intended to serve as a test or confirmation that their conclusions/assertions in preceding paragraphs have merit, then this is unconvincing – and problematic for supporting the assertions that the other functions (ionoregulation, cholesterol homeostasis, craniofacial development) ARE supported by pathway-level analysis.
Discussion section:
"We identified specific changes in the expression of key genes involved in the function or morphogenesis of individual tissues and organs with visible abnormalities. Given unaltered gene expression associated with apparently unaffected structures such as the liver, the DEGs in oil-exposed haddock indicate a disruption of specific developmental processes, as opposed to non-specific effects (e.g., general developmental delay)." – this appears to me as a gross overstatement. There was no organ-specific analysis of gene expression. The authors measured gene expression in whole animals. They therefore are not in a position to make assertions about expression in the heart or liver, though of course they are free to propose hypotheses.
"Our data demonstrate a transcriptional cascade that is tightly linked to these defects in cardiac function (cardiomyocyte intracellular calcium cycling) and form (heart chamber growth)." – there are 4 genes that are slightly down-regulated at one timepoint preceding the emergence of altered cardiac phenotypes (6 dpf), and none deferentially expressed immediately post-organogenesis (10 dpf) when the proposed E-C uncoupling should be apparent. Does this represent the discovery/demonstration of a "transcriptional cascade"?
"Crude oil likely disrupts normal MRC function" – what is the evidence for this?
"MRC function could be impaired by a high metabolic cost of PAH degradation." – why? Evidence or rationale to support this?
"In haddock embryos, fluid moves from the dorsal subdermal space to the yolk sac. At the larval stage, the permeable yolk sac membrane is replaced by the more resistant peritoneal cavity and body wall, causing fluid to move into the dorsal finfold and adjacent tissues." – I assume that the authors mean to write this as a proposal or hypothesis that is consistent with their observations?
Discussion section paragraph nine: much of this is highly speculative. There are plenty of transcriptomics studies of oil exposure during fish development, that include edema as an endpoint. If the authors want to make these assertions perhaps they should check those other studies for altered regulation of VEGF-C.
Discussion paragraph ten: This entire paragraph reads as excessive speculation
Discussion paragraph eleven: "but tissue localization during craniofacial development in embryonic fish is needed to confirm a role in this adverse outcome pathway" – look up Planchart & Mattingly 2010 TCDD upregulates FOXQ1 in zebrafish jaw primordium
Discussion section: "First, we used known spatial mRNA distributions in model species (primarily zebrafish) to more accurately phenotypically anchor the transcriptome data for crude oil-exposed haddock" – this sounds fancy but there is no description of this in the methods
Reviewer #2:
In their manuscript "Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish", Sorhus et al. characterize the impact on the transcriptome of PAHs during embryonic/larval development of Atlantic haddock. The authors examined gene expression profiles from fish exposed using three different exposure paradigms, two chronic (low =.58 µg/L & high = 6.7 µg/L) and one intermittent (6.1 µg/L per pulse). The authors examine genes that underlie one of four phenotypes observed in exposed embryos/larvae and show that changes in these genes may lead to the phenotypes observed. The authors' results will be of interest to the readership of eLife. However, several concerns exclude publication of the manuscript in its current form.
Major concerns:
1) Given that haddock is not a model system, better explanation of the developmental time windows will be critical for most readers to understanding the context of the embryonic and larval developmental stages. In zebrafish 50% epiboly occurs at 5.5 hpf while in haddock it occurs at 3 dpf.
2) Authors need to add Alcian images of the new facial phenotypes not listed in their previous work, the upper jaw is not shown in previously nor is the basicranium.
3) The authors need to be more explicit in their logic in moving from a broad overview of the transcriptome and to genes that directly impact the phenotypes listed when those genes are not the most highly responsive. Why discuss the broader transcriptional changes (subsection “General patterns of gene regulation in response to crude oil”)?
4) Without the ability to create genetic chimeras or cell-type specific transgenics to directly test the tissue target of PAHs effect on facial development, the authors need to soften their stance that PAH disrupts muscle development thereby affecting skeletal development. Crump and Schilling have shown that in zebrafish edn1 is expressed in and required in NCCs. The authors cannot rule out a direct NCC-PAH impact based only on their data (Discussion section).
5) What do the authors mean by the "transcriptional response to disrupted E-C coupling was not a simple pattern of compensation" (subsection “Genes associated with defects in cardiac function and morphogenesis”)?
6) Several the figures need to be reviewed. Two figure supplements are missing from the manuscript, though they are described in the text (subsection “A novel adverse outcome pathway: disruption of cholesterol homeostasis”). Figure 6B shows that myh1 is down regulated at 0 dph not 11 dpf as described in the text (subsection “Genes associated with craniofacial abnormalities”). The supporting datasets need to be organized with the supplemental table 1 for readability. It is difficult to relate these datasets with the table in the current form.
Reviewer #3:
Overall assessment: This is an interesting study that is well done overall and advances our understanding of oil impacts on developing marine fish. Strengths include the use of an environmental relevant non-model organism, the experimental design involving exposure at both embryonic and larval stages, use of multiple modes of exposure (two concentrations and a pulsed exposure), sampling at multiple time points, and a rich set of data on gene expression that is anchored to phenotypes.
Substantive concerns:
1) In the title, Abstract, and manuscript, the authors overreach when they invoke "chemical genetics" and "adverse outcome pathways" (AOPs). Chemical genetics involves high-throughput screening of libraries of individual compounds, a much more precise approach than the exposure to a complex chemical mixture performed here. Although nowhere defined by the authors, an AOP describes the entire sequence of events from a molecular initiating event, across multiple levels of biological organization, to an adverse outcome; it is synonymous with "mechanism of action" (see e.g. Ankley 2010 Envir. Tox. Chem. and Villeneuve 2011 Envir. Tox. Chem.). When the authors refer to AOPs they are actually referring only to adverse outcomes, not the pathways. Their gene expression data may help to inform our understanding of some AOPs associated with oil exposure, but they certainly have not "revealed novel AOPs." And contrary to the claim in the Abstract, it is not clear that they have "identified initiating events"-the specific chemical-protein interactions that lead to the gene expression and phenotypic changes that they report.
2) The authors' approach, which is stated explicitly (Results section), is to interpret changes in haddock gene expression in relation to known zebrafish mutants, i.e. they focused only on specific genes known to be involved in development of the tissues affected by oil. While this is valuable, is there a more objective approach that might be used to identify unexpected associations between changes in gene expression and specific phenotypes? It seems as though they have not taken full advantage of the unbiased RNA-seq dataset in this regard.
3) The authors over-interpret the connections between gene expression patterns and phenotypes, claiming cause-effect relationships in haddock from what are only associations. For example, Discussion section paragraph two claims a "tight linkage" between gene expression and cardiac defects. Whether the gene expression changes are causal or are secondary to the phenotypic changes is not clear. The authors could strengthen their arguments by being more explicit about the temporal and concentration-dependent associations between gene expression patterns and phenotypes, e.g. by adding a measure of phenotypic progression to Figures 4–7.
4) The experimental design is not described sufficiently. The methods (paragraph one) refer to a previous paper, but it is not clear whether these samples are from the same experiment described in that paper, or just used similar exposures. Indeed, the oil concentrations in that paper are expressed differently than in the current manuscript. Even if the methods are the same, this paper should be a description of the experimental design, including exposure conditions, numbers of replicates, numbers of pooled embryos in each replicate, etc.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for submitting your article "Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish" for consideration by eLife. Your article has been reviewed by three peer reviewers, and the evaluation has been overseen by a Reviewing Editor and Marianne Bronner as the Senior Editor. The reviewers have opted to remain anonymous.
The reviewers are generally happy with your revisions but request some clarifications in methodology and also raise some minor comments. I ask you to correct these points prior to passing on your manuscript for production. Individual reviews are below.
Reviewer #1:
Details of the experimental design and analysis are still not adequate. The authors refer the reader to the Sorhus 2016b publication for all aspects of experimental design. The authors should provide in the current manuscript at least the basics of experimental design in the methods section, including the number of doses, the specific developmental stages examined within each dose, and the number of biological replicates (replicate pools) within each dose*stage treatment. Furthermore, I have read Sorhus 2016b and it is not obvious WHAT were the exact contrasts without much effort on the part of the reader. The authors should make it easy to understand the basic of the experimental design.
Also, it is still not clear to me what was the unit of replication. I asked for clarification on the nature of replication previously. The revised manuscript reads "cDNA library preparation and sequencing was performed by the Norwegian Sequencing Centre (NSC, Oslo, Norway) on one pool from three replicate tanks per stage for each dose (low, pulse and high) plus control using the Illumina TruSeq RNA Sample Preparation Kit".
A treatment is a dose*stage. Were there replicate pools of embryos assayed per treatment? As far as I can tell, six developmental stages were profiled at each of 3 doses and a control, resulting in 24 treatments. I find 24 "samples" sequenced in the SRA (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP060012). This would suggest one pool per treatment, unless there are multiple samples embedded within each SRA entry. If there is just one sequenced pool per dose*stage treatment, then I recommend rejection of the manuscript for lack of replication of experimental units. If I misunderstand this, and there is in fact replication of pools of embryos within each dose*stage treatment, then I recommend the following:
Subsection “Structure of pelagic larvae and visible phenotypes associated with crude oil exposure”: "with 96% of high dose animals showing abnormal phenotypes, ranging to ~ 60% for pulse dose and ~ 35% for the low dose."
The authors should consider preparing a visual that summarizes this distribution of phenotypes across doses. E.g., stacked bar chart (or pie chart), including proportion of animals showing each phenotype, with dose as series. And perhaps separate plots for different developmental stages.
Subsection “Oil-induced changes in gene expression during embryonic development.”: "p > 0.05" Is it not more standard to notate as p<0.05?
Subsection “General patterns of gene expression in response to crude oil”: "At all stages, the subcategory of Organismal Development or Embryonic Development (henceforth combined as Development) was in the top 5 Diseases and Bio Functions category under Physiological System Development and Function with p values ranging from 10-3 to 10-19 " Are categories/subcategories of functions from IPA? Authors should state this.
In the same subsection: "Pathway enrichment was dose-dependent and clearly associated with the frequencies of abnormal phenotypes (Supplementary file 1C)." It would appear that reviewers do not have access to these supplementary files. This is frustrating.
Subsection “Genes associated with defects in cardiac function and morphogenesis”: "overexpression of bmp10, nkx25, or tbx3 is associated with serious heart defects in other vertebrates." This needs a citation
Discussion section: "Two major initiating events for crude oil-associated cardiac defects during fish development are chemical blockade of IKr repolarizing potassium currents, (encoded by kcnh2) and disruption of intracellular calcium handling, the latter culminating in sarcoplasmic reticulum (SR) calcium depletion through effects on either RyR or SERCA2 (encoded by ryr2 and at2a2, respectively). In the fully formed heart, these pharmacologic effects impair cardiac function by inducing arrhythmia and reducing contractility." Citations are needed in this section.
Subsection “Extraction of mRNA, RNA sequencing and bioinformatics”: "150 key genes involved in cardiac and craniofacial development and cardiac function were assembled." Please provide a table of these genes, and relevant citations describing their relationships to key phenotypes. The authors state in their rebuttal "The lists are provided in a new table, Supplementary file 1E." It would appear that I unfortunately do not have access to review these files (or I'm somehow just looking in the wrong place?). If not already done, authors please make sure these files are detailed with the phenotype relationship and citations to all relevant literature (e.g., defend criteria for including a gene in your curated set).
In the same subsection: "The contigs" is this the reference cod sequence, or the haddock RNA-seq read sequences?
The authors state in their rebuttal "Although the effects of oil exposure on salt and water balance in fish embryos have not been examined," This is perhaps true, but it has been examined in adults, and this may be worth noting in the manuscript. E.g.: Kennedy CJ, Farrell AP (2005) Ion homeostasis and interrenal stress responses in juvenile Pacific herring, Clupea pallasi, exposed to the water-soluble fraction of crude oil. Journal of Experimental Marine Biology and Ecology 323, 43-56.
Regarding the author's rebuttal "As for direct evidence of disrupted cholesterol homeostasis, we are not sure what is more direct than up-regulation of HMG-CoA-reductase" Up regulation of a gene is not a direct measure of altered cholesterol homeostasis.
Reviewer #2:
In their revised manuscript "Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish", Sorhus et al. characterize the impact on the transcriptome of PAHs during embryonic/larval development of Atlantic haddock. The authors examined gene expression profiles from fish exposed using three different exposure paradigms, two chronic (low =.58 µg/L & high = 6.7 µg/L) and one intermittent (6.1 µg/L per pulse). The authors examine genes that underlie one of four phenotypes observed in exposed embryos/larvae and show that changes in these genes may lead to the phenotypes observed. The authors have addressed concerned all concerns raised, therefore the results in the revised manuscript will be of interest to the readership of eLife and is ready for publishing in its current form.
Reviewer #3:
I would quibble with a few of the responses to my original comments, but I don't view these issues as serious enough to derail publication of the paper. Nevertheless, I point them out for the authors' consideration.
1) Regarding what is an "initiating event" in an adverse outcome pathway (AOP), I disagree with the authors' claim (response to reviewer #3, point 1) that bmp10 upregulation is such an initiating event. It is an early event, certainly. But the true initiating event is the chemical-protein interaction that causes bmp10 to be upregulated. That event has not been identified in this paper.
2) With regard to my questioning of their statement claiming that it has not been clear that oil influences mRNA levels as part of a developmental phenotype (Discussion section): The authors pointed out that the several previous papers I noted as showing oil affecting mRNA expression in fish did not involve measurements made in embryos. However, they did not provide this explanation in the revised manuscript, despite the fact that two of the three reviewers raised the same question. In addition, in their response the authors mention a more recent paper, published (27 June) prior to the original submission of this manuscript (30 Aug), that does include transcriptomic analysis of fish embryos exposed to oil (Xu et al. ES&T). It is unfortunate that the authors did not take advantage of the opportunity to better explain their original statement or compare their results with the prior work in embryos.
3) The manuscript would be improved by a clearer summary of the experimental design so that the readers don't have to dig out the other paper. And the point about replication is important – this needs clarification.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Thank you for resubmitting your work entitled "Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish" for further consideration at eLife. Your revised article has been favorably evaluated by Marianne Bronner as the Senior editor, a Reviewing editor, and one reviewer.
The manuscript is in principle ready for publication pending a minor but critical revision which is to include relevant information rather than referring the readers to other papers. The reviewers have asked for this repeatedly and I am only prepared to accept your paper with the inclusion of this additional information. This is an easy change and I hope you will make it quickly. Below is the comment from the reviewer.
Reviewer #1:
The revisions, in general are fine. However, I should note that I have been generally frustrated by the struggle with the authors to include all relevant information, accessibly, within this manuscript. Far too many papers are published these days that have incomplete description of methods. Also, many papers refer the reader to other papers to find out details of the methods. This results in un-necessary additional work for the reader – especially un-necessary since all journals these days have supplemental sections that allow authors to include all relevant information, without cluttering up the main manuscript. The authors are still insisting on sending the review on a hunt through their previous papers to find relevant information ("A thorough functional description of the genes including citations are provided in Sørhus et al. 2016a."). This is not too much to ask. We should all strive to make our research more transparent, and more reproducible.
https://doi.org/10.7554/eLife.20707.024Author response
As you will see from the individual reviews below, all of the reviewers had extensive issues with the analysis and lack of detail in the description of method. Even more importantly, however, they also felt that the there was a great deal of over-interpretation and "overselling of the work", rendering many statements unsupported by the data. We feel that this is fixable but will require considerable rewriting and re-review. If you feel you are able to revise the paper accordingly, we would be willing to consider a significantly revised version for re-review.
We appreciate the careful reviews and constructive comments from each of the three referees. We would like to initially provide a response that addresses many of the larger concerns of reviewers 1 & 3. Both helpfully pointed out that we failed to provide sufficient detail on our bioinformatics approach. Our original submission contained only a very brief description of our approach – i.e., “Individual genes involved in cardiac and craniofacial development, osmoregulation and lipid metabolism, or not directly linked to KEGG pathways were curated through manual searches of the literature (via PubMed and Web of Science) and expression databases at the Zebrafish Information Network (www.zfin.org).” We have explained our methods in much more detail in this revised submission.
Unlike conventional toxicogenomics, which relies primarily on bioinformatics to identify response pathways, our approach was more based on developmental genetics. We felt this was necessary because of the uncertainty inherent in cross-species extrapolations from a commercial ‘omics platform (e.g., designed for mammals) to a non-model fish (Atlantic haddock). In a genetic screen, investigators typically select the most robust (e.g., fully penetrant) phenotypes for in-depth analyses. In the original manuscript we focused on results from the high dose treatment group, because virtually 100% of the animals showed the abnormal phenotypes. This allowed for robust linkages between specific phenotypes and gene expression, and also kept the paper to a manageable length.
Prior to the current study, we characterized the normal developmental transcriptome of haddock (Sørhus et al., 2016; Dev. Biol. 411:301). This necessitated an extensive evaluation of the literature for vertebrate cardiac and craniofacial development to identify the most important genes for heart development/function as well as craniofacial/bone development (~ 150 genes for each category). The basic read count data were then used to identify all DEGs, and initial software-based statistical bioinformatics was used to compare changing patterns of gene expression at different developmental time points. Profiling the normal transcriptome was very labor-intensive; two members of our team worked on the project for the better part of a year.
This earlier gene identification and prioritization effort served as a cornerstone of our current study. This approach yielded very clear linkages between critical genes for which a loss or gain of function by mutation led to phenotypes resembling those observed in oil-exposed haddock. Moreover, we made novel observations of transcriptional regulation for clusters of functionally related genes, such as those involved in cholesterol synthesis and transport. This is a major new insight into developmental oil toxicity in fish, and will lead to new metrics for measuring injury at both the gross morphological (e.g., yolk absorption) and molecular scales.
The above explanation aside, we appreciate the reviewers’ comments related to a more conventional toxicogenomic analysis. Accordingly, during the revision phase, we have reanalyzed our data using Qiagen’s Ingenuity Pathway Analysis (IPA). We are pleased to report that each of the relationships identified in our original submission has been confirmed with statistical rigor. Therefore, our manual curation of the haddock transcriptome reinforces our confidence in an informatics application such as IPA to detect statistical enrichments of causal injury response pathways. This is uncommon, and we believe the paper is therefore strengthened. We also reiterate the importance of manual curation, as IPA also identified pathways that were misleading and apparently based on very tenuous associations in the literature.
Reviewer #1:
In this article the authors expand on previous experiments, where developing haddock are exposed to low levels of crude oil, to test for impacts on sensitive developmental endpoints. Developmental impacts are described in detail in previous reports, and this study reports transcriptomic changes during development following three oil exposure regimes.
I have two main difficulties with this paper. First, there is insufficient description of methods for gene expression analyses. Second, there is much speculation and unsupported statements such that much of the manuscript reads as over-statement.
Gene expression analysis:
Though the experimental design includes three different exposure regimes, specific gene transcriptional responses are not correlated in any way to dosing regime (high, low, pulse). The authors ought to show gene expression dose responses, at least for the genes for which they are trying to build cause-effect relationships with adverse outcomes. More convincing arguments for cause-effect would be supported by data showing that dysregulation of the gene precedes the appearance of the phenotype, and that dysregulation is more pronounced in treatments with more pronounced phenotypic effects (e.g., higher doses). The lack of detail reported here could make a reader concerned that perhaps interpretation of gene-specific responses are not as clear-cut as represented in the text. Are Figures 4, 5, 6, 7 showing data from any dose, or just high dose? Furthermore, there is no description of how emergence of phenotypes varied with dosing regime. Text (starting paragraph two of the Results section) that describes developmental phenotypes does not distinguish between effects that differ between doses, and all of the images from exposed animals in Figure 1 appear to be only from the high dose.
As noted above, we focused on the high dose treatment group because virtually 100% of the animals had the full range of abnormal phenotypes. This was clarified in the figures. However, as the reviewer suggests, we may not have provided a sufficiently detailed description of each phenotype – this was published previously (Sørhus et al., 2016; Sci. Rep. 6:31058). To address this, we have revised Figure 2 to more clearly relate changes in gene expression to the developmental onset of each phenotype. We also provided more information on the dose-dependency of gene expression in a new table (Supplementary file 1C), including the top IPA-identified pathways (e.g., Molecular and Cellular Function, etc.) for each dose and developmental time point.
Other aspects of gene expression analysis that are problematic:
There is no description of the statistical model. What was compared to what? How many tests were there? Was there multiple test correction given the size of the gene set (e.g., false discovery rate correction)? Was this experiment analyzed as a fully specified statistical model, or were individual treatments compared to individual control treatments separately from others, and if so was there multiple test correction for that? There is insufficient information provided here to be able to evaluate the quality and robustness of the gene expression analysis. Given that gene expression analysis is a focus of the paper this is a problem.
Each treatment (oil-exposed or unexposed controls) consisted of three pools.
We have used NOISeqBIO for the RNA-seq analysis where the 3 pools per treatment per lifestage analysed were treated as replicates. NOISeqBIO returns a DE (differential expression) probability that is equivalent to FDR (False discovery rate) adjusted P-values. The following reference has been added in subsection “Real time qPCR” to make this more clear:
Tarazona S, Furio-Tari P, Turra D, Pietro AD, Nueda MJ, Ferrer A, Conesa A (2015) Data quality aware analysis of differential expression in RNA-seq with NOISeq R/Bioc package. Nucleic Acids Research 43:e140.
41.9 million reads (subsection “Extraction of mRNA, RNA sequencing and bioinformatics”) – is this per sample?
Yes. To clarify, we changed the following sentence:
“The total number of reads averaged 41.39 million per sample and the mapping efficiency averaged 32.69%, giving an overall average of 13.51 million mapped paired end reads (125 bp) for each sample.”
32.69% mapping efficiency (subsection “Extraction of mRNA, RNA sequencing and bioinformatics”) – this is extremely low! Why is this so low if sequence similarity between haddock and reference (cod) genome is so high? How was successful mapping assessed – how were mapping parameters set? This should be a red flag that something went wrong in the mapping, or that there was significant contamination issues.
We mapped the haddock data to Atlantic cod because the genome of the latter is well annotated and nearly complete. The cod genes are listed with ensemble-IDs, allowing other investigators to extract the sequences to which we’ve mapped our haddock sequences. Also, the cod genome has served as a reference in several previous projects. Of 20613 sequences, 19000 cod genes map over 10 reads to haddock in at least one sample. We agree that 32% mapping efficiency appears low, but this is unlikely to reflect missing genes or a contamination issue. Rather, we strictly counted the reads that mapped exclusively to exons. Untranslated regions and alternatively spliced transcripts were not included. Also, the cod genome model is relatively fragmented, to the extent that RNA-seq in cod yields approximately 40% mapping to the exons of the model (Star et. al, 2011; Nature 477:207-210). Lastly, we expect some sequence variation across the two species.
After mapping, this left 31.51 million reads "for each group". I'm not sure what this means "for each group". It is generally accepted that a minimum of ~10M reads PER SAMPLE (that is, per experimental unit) is adequate for RNA-seq.
The number of reads per sample was high – i.e., the lowest number of reads for any sample after mapping was more than 10 M reads.
"three biological replicates per stage” – section "Total RNA and cDNA preparation" states that RNA was from pools of embryos. So are these more accurately 3 replicate pools?
This has been clarified in the Methods as follows: “cDNA library preparation and sequencing was performed by the Norwegian Sequencing Centre (NSC, Oslo, Norway) on one pool from three replicate tanks per stage for each dose (low, pulse and high) plus control using the Illumina TruSeq RNA Sample Preparation Kit”.
KEGG analysis: what were the criteria used to include a KEGG pathway in figures S4/S5? At least one of the criteria should be statistically significant enrichment. From my look at these figures, there are no obvious pathway-level connections to the phenotypes/physiological pathways on which the authors hang their hat. E.g., these figures do not have any terms that obviously relate to cholesterol homeostasis. Nor does Figure 3—figure supplement 1 indicate any enrichment of pathways having to do with ionoregulation, or craniofacial development. So what is the analysis here?
Please see our overarching response to this issue above (Editor’s comments). The additional analyses we’ve performed using IPA validates our manual curation (and vice versa), and also addresses the issue of statistical enrichment.
Selected gene sets: The authors make assertions about physiological mechanisms given responses of specific genes with inferred functions that make a nice story. However it is unclear what is the null expectation from this type of analysis.
E.g. Subsection “Genes associated with defects in cardiac function and morphogenesis”: "We therefore focused on genes associated with cardiomyocyte membrane potential and intracellular Ca2+ cycling" – how were these genes selected?
In the same section: "The expression of a few E-C coupling genes" – was there an unbiased collection of E-C coupling genes, then test for significant enrichment of these genes?
In the same section: "genes involved in cardiac morphogenesis" – again, how were these collected/curated as a gene set? Please provide the query gene set.
It seems more appropriate to consider null hypotheses in relation to specific structures or specific signaling pathways. For example, the hypothesis that genes governing phenotypically altered and unaltered structures would be equally responsive to crude oil. As stated in the subsection entitled “Unaltered gene expression in relation to visibly normal organs: lateral line and liver”, we essentially rejected this null hypothesis. For signaling pathways, a null hypothesis might be that Hedgehog signaling and Hox genes are affected in ways that are similar to BMP signaling. Our data do not support this hypothesis because developing haddock do not have phenotypes associated with disrupted Hedgehog or Hox signaling, and we observed minimal to no change in the transcription of these genes. As noted above, our priority gene lists were derived from an extensive evaluation of the developmental biology literature, and have since been confirmed through rigorous statistical analysis. The lists are provided in a new table, Supplementary file 1E.
Subsection “Genes associated with craniofacial abnormalities”: "We therefore interpreted developmental changes in haddock gene expression in the context of these well-characterised zebrafish mutants." – so the "tester" set here was all genes shown in zebrafish mutants to be involved in neural crest cell or muscle development? Please include this set of genes. The 12 genes indicated below this section – is this a significant enrichment?
The gene set is now provided in Supplementary file 1E.
Subsection “Genes associated with ion and water regulatory imbalance”: "Our analysis focused on key ionoregulatory proteins in MRCs and their associated genes, including Na+/K+ ATPase subunits (at1 genes, e.g. at1a1-a4, at1b1-b4), a urea transporter (ut1), a Na+/K+/2Cl- co-transporter (s12a2), the sodium hydrogen exchanger Nhe3 (slc9a3), and a chloride channel" – this is clearly not a complete or objective curation of ionoregulatory genes.
Is there an unbiased query gene set for any of these analyses? This whole "genes associated with phenotypes" section – is there objective analysis? What makes this anything more than just story telling given some cherry-picked genes? The authors should at least use language that proposes some hypotheses, rather than writing text so as to imply strong conclusions of cause-effect adverse outcome pathways. This is a problem that permeates the reporting of results and discussion, and contributes in a large way to overstatement of the results that are expanded on below.
For simplicity, we provided data on a small list of the most important actors in ionocyte function, as determined from a current review of the literature on iono- and osmoregulation in fish. The revised manuscript now provides a secondary analysis in IPA, statistically confirming alterations in a broader set of genes involved in ion transport.
Unsupported claims and overstatement issues:
There are pervasive issues throughout the manuscript, often associated with claims or conclusions that are (or should be) in fact stated as hypotheses. It is often unclear whether certain statements are supported by their data, supported by the literature, or are speculation or informed hypotheses. Some examples follow:
Subsection “Genes associated with defects in cardiac function and morphogenesis”: "Overall, the transcriptional response to disrupted E-C coupling was not a simple pattern of compensation" – how do they know that this is a transcriptional response to E-C coupling? Is expression not being measured in whole animals? What proportion of overall tissue mass is accounted for by the heart? I think the authors are drawing conclusions where they ought to be proposing hypotheses. This is an issue in MANY places throughout the manuscript
The role of transcription in the control of E-C coupling was discussed in detail in a previous paper that fully describes all aspects of the cardiac phenotypes (Sørhus et al., 2016; Sci. Rep. 6:31058). For brevity, we did not repeat that in-depth description in our original submission. However, it is clear from the reviewer’s comment that more clarification is needed, and we have revised the text in the “Genes associated with defects in cardiac function and morphogenesis” section accordingly. The current study specifically tested a hypothesis concerning transcriptional responses to disrupted E-C coupling. Very few studies have measured the responses of E-C coupling genes (i.e., ion channels) to pharmacologic blockade – these are cited in the original and revised version of our manuscript. Based on this limited information, we hypothesized a compensatory up-regulation of either the cognate gene for the blocked channel or other channel genes that function in the same pathway. This predicts, for example, an upregulation of the kcnh2 gene in response to potassium channel blockade. Because we observed the opposite, and did not see a compensatory upregulation of genes controlling repolarization, we concluded that the overall response was not a simple pattern of compensation. This suggests that crude oil is having a dual effect on E-C coupling at both the protein and gene levels.
In terms tissue specificity, as discussed in responses above, we assessed many genes that are expressed solely in the heart at a given developmental stage. Our RNA-Seq data demonstrate that these genes can be quantified despite the very contribution of the heart to total embryo mass. To address issues of tissue specificity, contribution of different tissues to the mass of embryos and larvae, and read count data, we added a new section and a table (Supplementary file 1D) to the Results.
Subsection “Genes associated with craniofacial abnormalities”: "Overall, the genes identified above for both neural crest and muscle lineages represent a more complex pattern of dysregulation than previously been reported for any of the individual zebrafish craniofacial mutants. This suggests that crude oil is acting on novel developmental processes." – It is unclear what is supporting this assertion (1st sentence), and how the assertion leads to the conclusion (2nd sentence)
Again, we wanted to keep this simple for brevity. There are some predictions about gene relationships based on the analysis of zebrafish mutants in terms of up- or down-regulation of signaling molecules, receptors, and downstream effectors. As is the case for E-C coupling genes, the pattern does not match the predictions based on single gene mutations in zebrafish. The oil-induced phenotype and the pattern of gene dysregulation are not a clear phenocopy of any one zebrafish craniofacial mutant. We have added text at this point in the manuscript to clarify.
Subsection “Genes associated with ion and water regulatory imbalance”: "and they should therefore lose water along a diffusion gradient if osmoregulation is disrupted as a consequence of heart and circulatory failure." – Does oil exposure cause ionoregulatory dysfunction and/or water loss in developing marine fish? Are these statements speculation/hypotheses, or are they supported in the literature? There are several prominent studies on the effects of oil exposure on osmoregulation in fish, but none of these are referred to.
We provided the justification for this statement in the original Discussion based on studies in saltwater medaka. Although the effects of oil exposure on salt and water balance in fish embryos have not been examined, there is a known role for cardiac circulation. Based on knockdown studies that disrupt circulation by causing heart defects, marine embryos can be expected to lose water and gain salt, with total osmolality increasing. We therefore infer that fluid accumulation in the vicinity of the yolk sac (edema) must represent free water moving from other compartments within the embryo. The visible phenotype is consistent with this – i.e., a reduction of the marginal subdermal space and a loss of cerebrospinal fluid (hindbrain ventricle). We also suggest several additional hypotheses regarding other mechanisms by which crude oil could disrupt ionoregulation more directly.
In the same section: "a third distinct oil-induced adverse outcome pathway." – how do they know that this is a pathway distinct from cardiac dysfunction? E.g., rather than just a secondary manifestation of circulatory failure?
We are drawing a distinction from the canonical cardiotoxicity pathways relating to sublethal heart malformation. Disturbances of salt and water balance, whether they are indirectly due to circulatory disruption or a consequence of direct PAH toxicity to ionocytes, are likely to represent adverse outcomes that have been previously unrecognized in this field, over three decades of research. For example, transcriptional changes observed in CNS-specific osmoregulatory genes may reflect a dysregulation of CNS development. Importantly, we also see changes in ion transport genes relatively early in development, prior to visible cardiac dysfunction. We have clarified this in the revised text.
"Crude oil exposures therefore appear to cause osmotic stress in the developing embryonic nervous system" – the authors are implying an overall impact on ionoregulatory abilities? Is there any support for this here or in the literature? E.g., any studies showing dysfunction of net sodium or chloride flux upon oil exposure? Do the authors mean to propose this as a hypothesis?
The visible phenotype, in tandem with specific patterns of gene expression, suggests a whole-embryo disruption of ionoregulation, and therefore a novel AOP. As in mammals, fluid accumulation in freshwater fish is a consequence of water retention. Given that marine embryos are hyposmolar to their environment, the only plausible way for edema to accumulate is a gain of electrolytes from the surrounding seawater. This is consistent with the observed reduction in the expression of genes that code for the transporters the keep salt out of the embryo. The gain of ions most likely is largest over the yolk sac due to surface area. This then leads to higher ion levels in the yolk sac sinus, and a shift of free water from other distal tissues, such as the dorsal subdermal space. We have rephrased parts of the text in this section. E.g.: “Our findings suggest that crude old may disrupt MRC function…”
Subsection: "A novel adverse outcome pathway: disruption of cholesterol homeostasis" – this subheading, and the following paragraphs, claim novelty, and appear to claim precedent for this discovery. Has oil exposure impacts on yolk utilization not previously been reported in the literature? I can find a publication from the 1990s after just a couple minutes searching. Furthermore, the studies presented here do not provide direct evidence for impacts on yolk mobilization or cholesterol homeostasis. E.g., how do the authors know that the apparently increased yolk size in Figure 1K isn't just a result of fluid accumulation stemming from observed edema?
As we clarify in this revised submission, edema accumulates well before there are detectable changes in the haddock yolk mass. Further, it seems unlikely that free water would move into the yolk platelets. Rather, edema accumulates in the yolk sac sinus, external to the yolk syncytial layer. We are aware of only one examination of yolk absorption, by our collaborators Mark Carls and Jeep Rice (published in 1999 on oil-exposed herring embryos). However, they did not make a connection to circulatory defects or lipid metabolism as we have in the present study. As for direct evidence of disrupted cholesterol homeostasis, we are not sure what is more direct than up-regulation of HMG-CoA-reductase. This is the rate-limiting enzyme for cholesterol synthesis, and its expression is extremely sensitive to sterol levels.
In the same section: "Pathway analysis was also consistent with a significant effect on cholesterol homeostasi." – The associated figure is missing. Also, staring at Figure 3—figure supplement 1 I see no pathways that include any terms obviously related to cholesterol homeostasis. Nor for that matter does Figure 1—figure supplement 1 indicate any enrichment of pathways having to do with ionoregulation, or craniofacial development,
We have addressed this through the new IPA analysis and associated tables.
Subsection “Unaltered gene expression in relation to visibly normal organs: lateral line and liver.”: "Similarly, the related KEGG pathways that are inclusive of these genes were relatively unaffected by oil exposure at all time points" – but what does this really mean? KEGG pathways (Figure 3—figure supplement 1 and 2), as far as I can tell, don't implicate many (any?) of the mechanisms that the authors are building their story on – e.g., ionoregulation, cholesterol homeostasis, craniofacial development. If this paragraph was intended to serve as a test or confirmation that their conclusions/assertions in preceding paragraphs have merit, then this is unconvincing – and problematic for supporting the assertions that the other functions (ionoregulation, cholesterol homeostasis, craniofacial development) ARE supported by pathway-level analysis.
We have addressed this through the new IPA analysis and associated tables.
Discussion section:
"We identified specific changes in the expression of key genes involved in the function or morphogenesis of individual tissues and organs with visible abnormalities. Given unaltered gene expression associated with apparently unaffected structures such as the liver, the DEGs in oil-exposed haddock indicate a disruption of specific developmental processes, as opposed to non-specific effects (e.g., general developmental delay)." – this appears to me as a gross overstatement. There was no organ-specific analysis of gene expression. The authors measured gene expression in whole animals. They therefore are not in a position to make assertions about expression in the heart or liver, though of course they are free to propose hypotheses.
We respectfully disagree, as captured (in part) in the discussion about null hypotheses above. While sophisticated tools are available for model species (e.g., organ-specific GFP expression) to isolate specific cell populations (e.g., by flow cytometry), this is practically impossible for a non-model species. A key aim of our study was to determine whether organ-specific changes in gene expression are detectable and consistent with the development of abnormal structures (phenotypic anchors). At issue is whether changes in gene expression, as visualized by whole-mount in situ hybridization, parallel quantitative changes in transcript abundance. Many previous studies of individual genes have confirmed this. As the reviewer pointed out earlier, the challenge is one of signal-to-noise, or whether gene regulation specific to a very small organ can be detected against a background of more abundant transcripts from a much larger tissue such as skeletal muscle. We have already confirmed this in our assessment of gene expression in normally developing haddock (Sørhus et al., 2016; Dev. Biol. 411:301). As structures became visible in the embryo over time (somites, heart, eye, etc.), we measured corresponding changes in the abundance of RNAs that are known to be associated with those structures via in situ hybridization. In this paper we show that this is also the case for these same structures that develop abnormally in response to crude oil. Consider, for example, bmp10. In situ hybridization in the three major vertebrate models (mouse, chick and zebrafish) has clearly shown that this gene is expressed exclusively in the heart during early development. In zebrafish, bmp10 is expressed in ~ 100 cells. Consistent with this, the read-count data for bmp10 in unexposed controls is just above the detection limit. In mutants with malformed hearts, bmp10 upregulation is always visible by in situ hybridization. Indeed, changes in both the patterns and levels of gene expression – as assessed by in situs – are a fundamental tool in developmental genetics to relate molecular pathways to formation of structure. In situs are regularly used semi-quantitatively, and we predicted similar changes in the RNA-Seq data. We have addressed this concern with additional text in the Discussion section and a new table (Supplementary file 1A) relating raw read-count data to tissue-specific expression.
"Our data demonstrate a transcriptional cascade that is tightly linked to these defects in cardiac function (cardiomyocyte intracellular calcium cycling) and form (heart chamber growth)." – there are 4 genes that are slightly down-regulated at one timepoint preceding the emergence of altered cardiac phenotypes (6 dpf), and none deferentially expressed immediately post-organogenesis (10 dpf) when the proposed E-C uncoupling should be apparent. Does this represent the discovery/demonstration of a "transcriptional cascade"?
This statement is not entirely accurate. The first change observed was a 4-fold upregulation of bmp10, an extremely potent morphogen. Subsequently we observed upregulation of the bmp10 target gene nkx2.5. Defects in E-C coupling at the protein level are likely occurring as soon as PAHs accumulate in the embryo. We are proposing that E-C coupling and excitation-transcription coupling are simultaneously disrupted by PAHs, by processes that involve SR and nuclear calcium handling. The bmp10-initiated transcriptional cascade is the link between SR calcium cycling and morphogenesis.
"Crude oil likely disrupts normal MRC function" – what is the evidence for this?
"MRC function could be impaired by a high metabolic cost of PAH degradation." – why? Evidence or rationale to support this?
The evidence for this is the differential regulation of MRC genes prior to visible defects in circulation. We are suggesting that the high metabolic demand of osmoregulation (the reason MRCs are mitochondria-rich) is compromised by an abnormal metabolic demand for PAH metabolism. We (and others) have shown that PAH exposure leads to massive upregulation of cyp1a throughout the epidermis in the embryos of multiple fish species, including haddock. Thus the epidermis is the major proximal target organ for PAHs, with an associated metabolic cost. The cyp1a response is uniform throughout the epidermis where the MRCs reside. Still unresolved is whether the ATP demands for PAH metabolism competes with the ATP available for ion transport, or whether PAHs or metabolites (or other compounds) directly interfere with channels, transporters or pumps as they do in EC coupling.
"In haddock embryos, fluid moves from the dorsal subdermal space to the yolk sac. At the larval stage, the permeable yolk sac membrane is replaced by the more resistant peritoneal cavity and body wall, causing fluid to move into the dorsal finfold and adjacent tissues." – I assume that the authors mean to write this as a proposal or hypothesis that is consistent with their observations?
This sentence has now been rephrased as a hypothesis that is consistent with our observations.
Discussion section paragraph nine: much of this is highly speculative. There are plenty of transcriptomics studies of oil exposure during fish development, that include edema as an endpoint. If the authors want to make these assertions perhaps they should check those other studies for altered regulation of VEGF-C.
Discussion paragraph ten: This entire paragraph reads as excessive speculation
In fact, there is not an abundance of transcriptomic studies for fish exposed embryonically to crude oil. To our knowledge, the first was just recently published (Xu et al., 2016, Environ. Sci. Technol. 50:7842) while this paper was under initial review, by a team we have previously collaborated with in the context of Deepwater Horizon impacts to mahi mahi in the Gulf of Mexico. There are more studies looking at the transcriptomes of fish larvae (vs. embryos), including several on cod, also by collaborators. A comparison to mahi mahi would be tenuous, as they are much less sensitive than haddock to crude oil exposure, with an edema accumulation phenotype that is mild by comparison (Edmunds et al., 2015, Sci. Rep. 5:17326). The genes discussed in this paragraph are exceptional in that they stand apart from the overall dataset, and an established literature (which we cite) is available to provide some context for interpreting these results. Our intention is not speculation, but rather a suggestion that these genes should be a focus for future studies.
Discussion paragraph eleven: "but tissue localization during craniofacial development in embryonic fish is needed to confirm a role in this adverse outcome pathway" – look up Planchart & Mattingly 2010 TCDD upregulates FOXQ1 in zebrafish jaw primordium
We are aware that TCDD does indeed upregulate foxq1 in zebrafish jaw primordium, but jaw malformation in TCDD-exposed zebrafish is entirely secondary to TCDD-induced heart defects. This has been rigorously demonstrated by a genetic gain-of-function approach (Lanham et al., 2014, Toxicol. Sci. 141:141). In the case of oil-exposed haddock, the craniofacial malformations are uncoupled from cardiac defects/edema formation, as discussed in Sørhus et al., 2016 (Sci. Rep. 6:31058). While Planchart and Mattingly did not demonstrate AHR-dependence of foxq1, this gene was shown to be AHR-responsive in a rat liver progenitor cell line (Faust et al., 2013, Arch. Toxicol. 87:681). Thus we cannot rule in or out a role for foxq1 in craniofacial malformations, or expression in the epidermis or liver with the rest of the AHR battery. This will likely need to be eventually resolved by in situ hybridization.
Discussion section: "First, we used known spatial mRNA distributions in model species (primarily zebrafish) to more accurately phenotypically anchor the transcriptome data for crude oil-exposed haddock" – this sounds fancy but there is no description of this in the methods
The requested description was added to the Methods. We simply utilized zfin.org to identify expression patterns that have been determined in previous studies by the zebrafish community. The search tool on the zfin.org landing page allows for detailed advanced search capabilities, either by gene, tissue, developmental stage, etc. Absent data in the Zfin database, we searched PubMed for papers showing in situ data in either chick or mouse.
Reviewer #2:
[…] Major concerns:
1) Given that haddock is not a model system, better explanation of the developmental time windows will be critical for most readers to understanding the context of the embryonic and larval developmental stages. In zebrafish 50% epiboly occurs at 5.5 hpf while in haddock it occurs at 3 dpf.
This is an excellent suggestion. In response, we have added information on the timing of normal haddock development. We also provide a new figure (Figure 2) that describes the developmental onset of each phenotype investigated.
2) Authors need to add Alcian images of the new facial phenotypes not listed in their previous work, the upper jaw is not shown in previously nor is the basicranium.
The basicranium was shown (labeled as trabeculae crania/ethmoid plate) in our previous paper. In cod and haddock, the upper jaw is not at a sufficient stage of chondrogenesis to be detected by Alcian blue until much later in the larval period (45 days post hatch for cod; von Herbring 2001, J. Fish Biol. 59:767).
3) The authors need to be more explicit in their logic in moving from a broad overview of the transcriptome and to genes that directly impact the phenotypes listed when those genes are not the most highly responsive. Why discuss the broader transcriptional changes (subsection “General patterns of gene regulation in response to crude oil”)?
The addition of results from IPA has significantly expanded the broad overview aspect of this study. This demonstrates that the pathways that contain our individual genes of interest are among the most highly responsive to crude oil exposure.
4) Without the ability to create genetic chimeras or cell-type specific transgenics to directly test the tissue target of PAHs effect on facial development, the authors need to soften their stance that PAH disrupts muscle development thereby affecting skeletal development. Crump and Schilling have shown that in zebrafish edn1 is expressed in and required in NCCs. The authors cannot rule out a direct NCC-PAH impact based only on their data (Discussion section).
This is a fair point, and although we believe this is the most parsimonious explanation, we have qualified our interpretation. We discount a role for NCCs primarily because the phenotype does not appear to be a patterning defect, as would be expected from a disruption of NCC regulators. For example, in zebrafish, a loss of edn1 dorsalizes the ventral jaw structure, and ectopic/excess edn1 ventralizes the dorsal jaw structure.
5) What do the authors mean by the "transcriptional response to disrupted E-C coupling was not a simple pattern of compensation" (subsection “Genes associated with defects in cardiac function and morphogenesis”)?
See response to this same question from reviewer 1.
6) Several the figures need to be reviewed. Two figure supplements are missing from the manuscript, though they are described in the text (subsection “A novel adverse outcome pathway: disruption of cholesterol homeostasis”). Figure 6B shows that myh1 is down regulated at 0 dph not 11 dpf as described in the text (subsection “Genes associated with craniofacial abnormalities”). The supporting datasets need to be organized with the supplemental table 1 for readability. It is difficult to relate these datasets with the table in the current form.
These figures were mislabeled in the previous version of the manuscript. Both text and figure labeling have now been corrected. Our intent was to show which 10 genes are most up-regulated and down-regulated at each several stages. To visualize when some genes were highly expressed at several stages, we collapsed the row were most up or down regulated at several stages during development. We have rewritten the Table legend to clarify this.
Reviewer #3:
[…] Substantive concerns:
1) In the title, Abstract, and manuscript, the authors overreach when they invoke "chemical genetics" and "adverse outcome pathways" (AOPs). Chemical genetics involves high-throughput screening of libraries of individual compounds, a much more precise approach than the exposure to a complex chemical mixture performed here. Although nowhere defined by the authors, an AOP describes the entire sequence of events from a molecular initiating event, across multiple levels of biological organization, to an adverse outcome; it is synonymous with "mechanism of action" (see e.g. Ankley 2010 Envir. Tox. Chem. and Villeneuve 2011 Envir. Tox. Chem.). When the authors refer to AOPs they are actually referring only to adverse outcomes, not the pathways. Their gene expression data may help to inform our understanding of some AOPs associated with oil exposure, but they certainly have not "revealed novel AOPs." And contrary to the claim in the Abstract, it is not clear that they have "identified initiating events"-the specific chemical-protein interactions that lead to the gene expression and phenotypic changes that they report.
Our view of chemical genetics is broader and not necessarily confined to high-throughput screens of small molecule libraries. This is what the field has evolved into, but chemical genetics is simply the premise that small molecules can produce abnormal developmental phenotypes in a manner similar to single-gene mutations. Indeed, the logic of the small molecule screening approach was derived from earlier focal studies of teratogens that caused malformations through specific interactions with developmental signaling molecules, such as cyclopamine-hedgehog signaling (Incardona et al., 1998, Development 125:3553). This history is captured in the Introduction to the first published zebrafish small molecule screen (Peterson et al., 2000, Proc. Nat. Acad. Sci. 97:12965). In traditional forward developmental genetics, molecular pathways are identified by relating mutant phenotypes to the physical location of the mutation on the chromosome. On one hand, chemical genetics is simply the use of compounds to alter gene function at the protein end, thereby working towards the same goal (pathway identification) in the opposite direction. Chemically induced phenotypes can provide a greater level of control over the timing of gene function, for example, and sidestep experimental challenges such as embryonic lethality. The power of genetics – chemical or otherwise – is derived from the anchoring of gene expression to a careful phenotypic quantification and anchoring to gene expression. This was the goal of the present study, and we believe the use of the term is defensible. We also wanted to emphasize the feasibility of this approach in real-world organisms to study traits that may be absent from the lab models typically used in high throughput screening. Finally, as a study focusing on gene-environment interactions, we draw additional comparisons to analyses of natural compounds. Cyclopamine was initially identified from a similarly complex mixture of compounds produced by a range plant (Veratrum viridae) that induced malformations in pregnant sheep.
We agree that a goal in AOP development is to describe a series of related events across biological scales. But this is a process that plays out across many studies conducted over years or decades of research. In this regard, the framework for scaling crude oil-driven adverse health impacts to fish is relatively advanced, as a consequence of > 25 years of science initially spurred by the 1989 Exxon Valdez spill in Prince William Sound. The early life stage-based cardiotoxicity AOP was the focus of a recent review (Incardona and Scholz, 2016, Aquat. Toxicol. 177:555). In the current paper, we refer to population-scale adverse outcomes as the motivation for conducting this research, to inform future fisheries management by the National Marine Fisheries Service in the U.S. and the Institute for Marine Research in Norway. In terms of the actual data generated, our current study clarifies the relationships between initiating events and organ- or organism-level adverse outcomes. In certain cases our findings represent clear initiating events or key events strongly supported by the literature – e.g., bmp10 upregulation, consequent nkx2.5 upregulation, and ventricular malformation. In other cases we propose steps in a possible AOP that will require additional study – e.g., water shifts in the developing neural tube, aquaporin up- or downregulation, and subsequent abnormal CNS development. Although semantic, the distinctions are important, and we have provided additional background on AOPs and their development in the revised text.
2) The authors' approach, which is stated explicitly (Results section), is to interpret changes in haddock gene expression in relation to known zebrafish mutants, i.e. they focused only on specific genes known to be involved in development of the tissues affected by oil. While this is valuable, is there a more objective approach that might be used to identify unexpected associations between changes in gene expression and specific phenotypes? It seems as though they have not taken full advantage of the unbiased RNA-seq dataset in this regard.
We have used the unbiased RNA-seq dataset as a tool. As described above (response to Editor), we conducted an extensive manual curation of all DEGs (Kmeans, pathways analysis, inspecting manually the genes underlying the most regulated pathways, most regulated genes) before focusing on a final set of genes. Moreover, we have since added the results from an unbiased IPA. The reviewer may be more concerned about false negatives, or missing unidentified genes that might play key roles in heart development, for example. Given current state of technology, this is likely an unavoidable limitation of doing functional genomics in a non-model species where it is extremely difficult to manipulate the embryo with the types of tools available to zebrafish. In addition, IPA will only make associations based on previously published data. Nevertheless, given strong positive results with known potent signaling molecules (bmp10, wnt9b, edn1) and key transcription factors, it seems unlikely that we are missing novel key genes. These genes alone are sufficient to explain the phenotypes. In future studies we intend to develop in situ hybridization probes for some of the highly regulated genes for which there are no tissue localization data available, to further define their role in the crude oil injury phenotypes.
3) The authors over-interpret the connections between gene expression patterns and phenotypes, claiming cause-effect relationships in haddock from what are only associations. For example, Discussion section paragraph two claims a "tight linkage" between gene expression and cardiac defects. Whether the gene expression changes are causal or are secondary to the phenotypic changes is not clear. The authors could strengthen their arguments by being more explicit about the temporal and concentration-dependent associations between gene expression patterns and phenotypes, e.g. by adding a measure of phenotypic progression to Figures 4–7.
As discussed in the responses to reviewer 1, we have provided additional data showing how changes in gene expression precede the developmental emergence of each visible phenotype. We have also revised the text to support our conclusions of tight cause-effect linkages based on gene dosage studies from zebrafish. Small changes in the expression of potent morphogens such as bmp10 and master transcriptional regulators such nkx2.5 produce large changes in developmental trajectory. For morphogens in particular, this is how they work in the process of specifying cell fate. For the most potent ligand in the bmp family, bmp10, too much or too little protein leads to defective heart development. The cumulative data from both mouse and zebrafish indicate that it would be impossible to form a normal heart if bmp10 were inappropriately upregulated too early in embryonic development, as we observed here. Morphogen effects are exquisitely dose-dependent, such that even less than 2-fold changes can lead to different target gene activation.
4) The experimental design is not described sufficiently. The methods (paragraph one) refer to a previous paper, but it is not clear whether these samples are from the same experiment described in that paper, or just used similar exposures. Indeed, the oil concentrations in that paper are expressed differently than in the current manuscript. Even if the methods are the same, this paper should be a description of the experimental design, including exposure conditions, numbers of replicates, numbers of pooled embryos in each replicate, etc.
We have clarified the source of the samples (the same as for Sørhus et al., 2016, Sci. Rep. 6:31058). Oil concentrations have also been corrected to match the previous paper.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Reviewer #1:
Details of the experimental design and analysis are still not adequate. The authors refer the reader to the Sorhus 2016b publication for all aspects of experimental design. The authors should provide in the current manuscript at least the basics of experimental design in the methods section, including the number of doses, the specific developmental stages examined within each dose, and the number of biological replicates (replicate pools) within each dose*stage treatment. Furthermore, I have read Sorhus 2016b and it is not obvious WHAT were the exact contrasts without much effort on the part of the reader. The authors should make it easy to understand the basic of the experimental design.
Also, it is still not clear to me what was the unit of replication. I asked for clarification on the nature of replication previously. The revised manuscript reads "cDNA library preparation and sequencing was performed by the Norwegian Sequencing Centre (NSC, Oslo, Norway) on one pool from three replicate tanks per stage for each dose (low, pulse and high) plus control using the Illumina TruSeq RNA Sample Preparation Kit".
A treatment is a dose*stage. Were there replicate pools of embryos assayed per treatment? As far as I can tell, six developmental stages were profiled at each of 3 doses and a control, resulting in 24 treatments. I find 24 "samples" sequenced in the SRA (https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP060012). This would suggest one pool per treatment, unless there are multiple samples embedded within each SRA entry. If there is just one sequenced pool per dose*stage treatment, then I recommend rejection of the manuscript for lack of replication of experimental units. If I misunderstand this, and there is in fact replication of pools of embryos within each dose*stage treatment, then I recommend the following:
The 24 samples the reviewer is referring to is the samples from Sørhus et.al 2016a DevBio. The samples from this study, however, should now be available in the SRA database under this accession ID: PRJNA328092.
Additional information has been added to the Methods section. To make it clear: A total of 132 samples were sequenced, and 126 should be included in SRA (one pulsed group tank appeared to be an outsider and was excluded from further analysis, a matter that is thoroughly explained in the Results section Sørhus 2016b). We have clarified this further in the text. " cDNA library preparation and sequencing was performed by the Norwegian Sequencing Centre (NSC, Oslo, Norway) using the Illumina TruSeq RNA Sample Preparation Kit.A total of 132 samples were sequenced and 126 were subjected for analysis: Three (control, low and high dose) or two (pulse) biological replicates for each dose from 6 stages during and after embryonic exposure and three biological replicates for each dose from 5 stages during larval exposure.”
Subsection “Structure of pelagic larvae and visible phenotypes associated with crude oil exposure”: "with 96% of high dose animals showing abnormal phenotypes, ranging to ~ 60% for pulse dose and ~ 35% for the low dose."
The authors should consider preparing a visual that summarizes this distribution of phenotypes across doses. E.g., stacked bar chart (or pie chart), including proportion of animals showing each phenotype, with dose as series. And perhaps separate plots for different developmental stages.
We understand the benefit of having a visual summary in this paper. However, the Sørhus et.al 2016b paper also includes both figures and tables showing a clear dose dependent effect for abnormal phenotypes, functional and development cardiac abnormalities and edema accumulation. Thus, including the exact same figures and tables should therefore not be necessary.
Subsection “Oil-induced changes in gene expression during embryonic development.”: "p > 0.05" Is it not more standard to notate as p<0.05?
Corrected
Subsection “General patterns of gene expression in response to crude oil”: "At all stages, the subcategory of Organismal Development or Embryonic Development (henceforth combined as Development) was in the top 5 Diseases and Bio Functions category under Physiological System Development and Function with p values ranging from 10-3 to 10-19 " Are categories/subcategories of functions from IPA? Authors should state this.
Yes, these are the classifications provided in the IPA output. We have provided additional detail in the Methods section to highlight this, in addition to stating this in the main text.
In the same subsection: "Pathway enrichment was dose-dependent and clearly associated with the frequencies of abnormal phenotypes (Supplementary file 1C)." It would appear that reviewers do not have access to these supplementary files. This is frustrating.
We have confirmed that these lists have been uploaded, and we have requested them to be available for the reviewers in our letter to editor.
Subsection “Genes associated with defects in cardiac function and morphogenesis”: "overexpression of bmp10, nkx25, or tbx3 is associated with serious heart defects in other vertebrates." This needs a citation
The citations Chen 2006, Ribeiro 2007 and Tu 2009 have been inserted.
Discussion section: "Two major initiating events for crude oil-associated cardiac defects during fish development are chemical blockade of IKr repolarizing potassium currents, (encoded by kcnh2) and disruption of intracellular calcium handling, the latter culminating in sarcoplasmic reticulum (SR) calcium depletion through effects on either RyR or SERCA2 (encoded by ryr2 and at2a2, respectively). In the fully formed heart, these pharmacologic effects impair cardiac function by inducing arrhythmia and reducing contractility." Citations are needed in this section.
The citations Brette 2014 and Incardona 2009 and 2014 have been inserted.
Subsection “Extraction of mRNA, RNA sequencing and bioinformatics”: "150 key genes involved in cardiac and craniofacial development and cardiac function were assembled." Please provide a table of these genes, and relevant citations describing their relationships to key phenotypes. The authors state in their rebuttal "The lists are provided in a new table, Supplementary file 1E." It would appear that I unfortunately do not have access to review these files (or I'm somehow just looking in the wrong place?). If not already done, authors please make sure these files are detailed with the phenotype relationship and citations to all relevant literature (e.g., defend criteria for including a gene in your curated set).
Hopefully the reviewers have now been allowed to view the supplementary files. The table include name and relation to phenotype. A thorough functional description of the genes including citations are provided in Sørhus et al. 2016a. We also adjusted the supplementary file legend to: Supplementary file 1F: Manually curated list of genes involved in cardiac development and function and craniofacial and bone development and maintenance. For references see Sørhus et. al 2016a, DOI: 10.1016/j.ydbio.2016.02.012.
SP, swissprot; GB, genebank.
In the same subsection: "The contigs" is this the reference cod sequence, or the haddock RNA-seq read sequences?
This refers to the cod genes. We have now clarified this in the text: “The cod genes[…]”
The authors state in their rebuttal "Although the effects of oil exposure on salt and water balance in fish embryos have not been examined," This is perhaps true, but it has been examined in adults, and this may be worth noting in the manuscript. E.g.: Kennedy CJ, Farrell AP (2005) Ion homeostasis and interrenal stress responses in juvenile Pacific herring, Clupea pallasi, exposed to the water-soluble fraction of crude oil. Journal of Experimental Marine Biology and Ecology 323, 43-56.
Very good point. This point has now been added to the manuscript text.
Regarding the author's rebuttal "As for direct evidence of disrupted cholesterol homeostasis, we are not sure what is more direct than up-regulation of HMG-CoA-reductase" Up regulation of a gene is not a direct measure of altered cholesterol homeostasis.
Actually, in this case up-regulation of HMG-CoA reductase (hmdh) is unequivocally a direct measure of altered cholesterol homeostasis. This is so well established it can be found in any cell biology text book, and is the subject of the Brown and Goldstein review we cited in the Results section when first mentioning hmdh. However, we have made this more explicit in the relevant Discussion section with the statement “Cellular sterol levels are tightly controlled by membrane-bound transcription factors, which are cleaved and activated when membrane cholesterol levels fall, leading to transcription of hmdh, the rate-limiting enzyme for cholesterol synthesis”, citing Brown and Goldstein again. Like CYP1A induction in response to aromatic xenobiotics, this is a classic example of transcriptional control of metabolism, harking back to Jacob and Monod and the lac operon in E. coli.
Reviewer #3:
I would quibble with a few of the responses to my original comments, but I don't view these issues as serious enough to derail publication of the paper. Nevertheless, I point them out for the authors' consideration.
1) Regarding what is an "initiating event" in an adverse outcome pathway (AOP), I disagree with the authors' claim (response to review 3) that bmp10 upregulation is such an initiating event. It is an early event, certainly. But the true initiating event is the chemical-protein interaction that causes bmp10 to be upregulated. That event has not been identified in this paper.
That is true. We do not in this paper provide the direct evidence for the first initiating event. However, in the paper we suggest that the effect on ETC, resulting in an inappropriate calcium activation of myocardin as the initiating event, leading to the up-regulation of bmp10.
2) With regard to my questioning of their statement claiming that it has not been clear that oil influences mRNA levels as part of a developmental phenotype (Discussion section): The authors pointed out that the several previous papers I noted as showing oil affecting mRNA expression in fish did not involve measurements made in embryos. However, they did not provide this explanation in the revised manuscript, despite the fact that two of the three reviewers raised the same question. In addition, in their response the authors mention a more recent paper, published (27 June) prior to the original submission of this manuscript (30 Aug), that does include transcriptomic analysis of fish embryos exposed to oil (Xu et al. ES&T). It is unfortunate that the authors did not take advantage of the opportunity to better explain their original statement or compare their results with the prior work in embryos.
The paper by Xu et al. in ES&T only looked a gene expression in hatched larvae, not at earlier time point during embryogenesis, and in particular prior to the appearance of abnormal phenotypes. Therefore, none of the papers in question are really comparable. We are uncertain of the value in adding text and additional citations to point this out.
3) The manuscript would be improved by a clearer summary of the experimental design so that the readers don't have to dig out the other paper. And the point about replication is important – this needs clarification.
Additional information has been added to the Materials and methods section.
[Editors' note: further revisions were requested prior to acceptance, as described below.]
Reviewer #1:
The revisions, in general are fine. However, I should note that I have been generally frustrated by the struggle with the authors to include all relevant information, accessibly, within this manuscript. Far too many papers are published these days that have incomplete description of methods. Also, many papers refer the reader to other papers to find out details of the methods. This results in un-necessary additional work for the reader – especially un-necessary since all journals these days have supplemental sections that allow authors to include all relevant information, without cluttering up the main manuscript. The authors are still insisting on sending the review on a hunt through their previous papers to find relevant information ("A thorough functional description of the genes including citations are provided in Sørhus et al. 2016a."). This is not too much to ask. We should all strive to make our research more transparent, and more reproducible.
A reference column has been added to the Supplementary file 1F including one or more reference for every gene.
https://doi.org/10.7554/eLife.20707.025Article and author information
Author details
Funding
Norges Forskningsråd (Project no. 234367)
- Elin Sørhus
- John P Incardona
- Tomasz Furmanek
- Nathaniel L Scholz
- Sonnich Meier
- Rolf B Edvardsen
- Sissel Jentoft
VISTA Foundation (Project no. 6161)
- Elin Sørhus
Institute of Marine Research (Project no. 14236)
- Elin Sørhus
- Tomasz Furmanek
- Sonnich Meier
- Rolf B Edvardsen
The funders had no role in study design, data collection and interpretation, or the decision to submit the work for publication.
Acknowledgements
We acknowledge Michal Rejmer, Domagoj Maksan, Ørjan Karlsen and Terje van der Meeren for providing fish husbandry. Funded by the Research Council of Norway (Project no. 234367, www.forskningsradet.no), the VISTA foundation (Project no. 6161, www.vista.no) and the Institute of Marine Research, Norway. The University of Oslo’s Norwegian Sequencing Centre (NSC; http://www.sequencing.uio.no) created the library and sequenced the transcriptome. Funding organizations had no role in study design, data collection and analysis, or manuscript preparation.
Ethics
Animal experimentation: All animal experiments within the study were approved by NARA, the governmental Norwegian Animal Research Authority (http://www.fdu.no/fdu/, reference number 2012/275334-2). All methods were performed in accordance with approved guidelines.
Reviewing Editor
- Marianne Bronner, California Institute of Technology, United States
Version history
- Received: August 26, 2016
- Accepted: January 20, 2017
- Accepted Manuscript published: January 24, 2017 (version 1)
- Version of Record published: February 10, 2017 (version 2)
Copyright
This is an open-access article, free of all copyright, and may be freely reproduced, distributed, transmitted, modified, built upon, or otherwise used by anyone for any lawful purpose. The work is made available under the Creative Commons CC0 public domain dedication.
Metrics
-
- 1,763
- Page views
-
- 295
- Downloads
-
- 75
- Citations
Article citation count generated by polling the highest count across the following sources: Crossref, Scopus, PubMed Central.
Download links
A two-part list of links to download the article, or parts of the article, in various formats.
Downloads (link to download the article as PDF)
Open citations (links to open the citations from this article in various online reference manager services)
Cite this article (links to download the citations from this article in formats compatible with various reference manager tools)
Novel adverse outcome pathways revealed by chemical genetics in a developing marine fish
eLife 6:e20707.
https://doi.org/10.7554/eLife.20707
Further reading
-
- Ecology
- Plant Biology
Global agro-biodiversity has resulted from processes of plant migration and agricultural adoption. Although critically affecting current diversity, crop diffusion from Classical antiquity to the Middle Ages is poorly researched, overshadowed by studies on that of prehistoric periods. A new archaeobotanical dataset from three Negev Highland desert sites demonstrates the first millennium CE&'s significance for long-term agricultural change in southwest Asia. This enables evaluation of the 'Islamic Green Revolution' (IGR) thesis compared to 'Roman Agricultural Diffusion' (RAD), and both versus crop diffusion during and since the Neolithic. Among the finds, some of the earliest aubergine (Solanum melongena) seeds in the Levant represent the proposed IGR. Several other identified economic plants, including two unprecedented in Levantine archaeobotany-jujube (Ziziphus jujuba/mauritiana) and white lupine (Lupinus albus)-implicate RAD as the greater force for crop migrations. Altogether the evidence supports a gradualist model for Holocene-wide crop diffusion, within which the first millennium CE contributed more to global agricultural diversity than any earlier period.
-
- Ecology
- Evolutionary Biology
Temperature determines the geographical distribution of organisms and affects the outbreak and damage of pests. Insects seasonal polyphenism is a successful strategy adopted by some species to adapt the changeable external environment. Cacopsylla chinensis (Yang & Li) showed two seasonal morphotypes, summer-form and winter-form, with significant differences in morphological characteristics. Low temperature is the key environmental factor to induce its transition from summer-form to winter-form. However, the detailed molecular mechanism remains unknown. Here, we firstly confirmed that low temperature of 10 °C induced the transition from summer-form to winter-form by affecting the cuticle thickness and chitin content. Subsequently, we demonstrated that CcTRPM functions as a temperature receptor to regulate this transition. In addition, miR-252 was identified to mediate the expression of CcTRPM to involve in this morphological transition. Finally, we found CcTre1 and CcCHS1, two rate-limiting enzymes of insect chitin biosyntheis, act as the critical down-stream signal of CcTRPM in mediating this behavioral transition. Taken together, our results revealed that a signal transduction cascade mediates the seasonal polyphenism in C. chinensis. These findings not only lay a solid foundation for fully clarifying the ecological adaptation mechanism of C. chinensis outbreak, but also broaden our understanding about insect polymorphism.
