Ontogeny-Specific Skeletal Deformities in Atlantic Haddock Caused by Larval Oil Exposure

Bone deformities are one of the main effects of crude oil exposure in marine fish larvae. Craniofacial and jaw deformities, if severe enough, may restrict feeding and ultimately kill the developing larvae. This study aimed to examine the impact of dispersed crude oil on bone development in Atlantic haddock (Melanogrammus aeglefinus) larvae, a fish species spawning in areas approached for oil and gas exploration in the North Atlantic Ocean. Atlantic haddock larvae were exposed to low (60 μg oil/L), high (600 μg oil/L), or pulsed (0–600, average 60 μg oil/L over time) dispersed crude oil from 0 to 18 days post hatch (dph). Endpoints included survival and growth, bone integrity, and transcriptional parameters, which were assessed during (0–18 dph) and after exposure until the fish reached 8 months of age (243 dph). The results showed that the larvae in the high treatment group had reduction in growth at 2–19, 44, 134, and 243 dph. Craniofacial abnormalities were most severe at 8 and 19 dph. These deformities were not present at 44 dph, possibly because the larvae with deformed jaws failed to feed properly and died. Higher prevalence of spinal deformities was observed in haddocks that survived for 243 dph. Three genes encoding proteins critical for osteoblast function, sp7, postn, and col10a1, were downregulated in the high treatment group larvae. We discuss possible mechanisms of action in the developing larvae after oil exposure. In conclusion, this study shows that larval exposure to oil can potentially have long-term effects on growth and bone integrity in Atlantic haddock.


INTRODUCTION
Research conducted in the aftermath of major accidental oil spills, such as the 1989 grounding of Exxon Valdez in the Prince William Sound, Alaska, and the 2010 Deepwater Horizon blowout in the Gulf of Mexico, has documented developmental abnormalities in fish larvae exposed to crude oil (Peterson et al., 2003;Rice et al., 2007;Incardona et al., 2013). As verified in laboratory studies, exposure to crude oil or its derivatives, such as polycyclic aromatic hydrocarbons (PAHs), induces craniofacial, jaw, and eye deformities, and spinal curvature, as well as pericardial and yolk-sac edema and reduced growth in fish larvae (Carls et al., 1999;Incardona et al., 2004;Ingvarsdottir et al., 2012;Incardona and Scholz, 2016;Sørhus et al., 2017;Hansen et al., 2019;Laurel et al., 2019;Lie et al., 2019;Sun et al., 2019;Johann et al., 2020). Morphological alterations in the head skeleton, including the jaw, have been considered as a secondary effect of reduced blood flow or heart failure after crude oil exposure (Incardona et al., 2004(Incardona et al., , 2013. Developmental abnormalities related to bone formation, therefore, appear to be one of the major impacts of the components in crude oil on fish larvae. With changing climate and melting sea ice, several Arctic nations have planned to expand oil and gas exploration northward into pristine northern waters. This new use of the marine environment raises competing interests and concerns regarding potential impacts on fisheries. For example, the Barents Sea, an Arctic shelf area north of Norway and Russia, is among the richest fishing grounds and home to the largest stock of Atlantic cod in the world (Gadus morhua). Atlantic cod and most of the other major Barents Sea stocks spawn off the coast of Northern Norway, particularly in the area of Lofoten and Vesterålen (Olsen et al., 2010), which constitutes a large spawning area for several species of commercially important fish stocks, such as Atlantic cod, Atlantic haddock (Melanogrammus aeglefinus), and Atlantic herring (Clupea harengus) (Olsen et al., 2010). Because of environmental concerns, the Lofoten and Vesterålen area has, so far, not been opened to oil exploration by the Norwegian Government (MPE, 2011). However, as oil exploration in the Lofoten and Vesterålen area has been a hotly debated issue, there is a need for more knowledge of the potential impact of major oil spills on early life stages of fish species spawning in this region.
Evidence suggests that the sensitivity to dispersed crude oil during early developmental stages may differ considerably among fish larvae from very similar species. For example, Sørhus et al. (2015) showed that Atlantic haddock embryos were more sensitive to oil exposure than Atlantic cod embryos, two closely related gadoid species, because of the accumulation of dispersed oil micro-droplets, possibly resulting from enhanced uptake of toxic components such as PAHs. The authors suggested that the differential sensitivity to oil arises from dissimilar lamella structure of the chorion, with the Atlantic haddock chorion generally being more adherent than the Atlantic cod chorion (Sørensen et al., 2017). This underlines the need of speciesspecific knowledge to be able to generate models predicting the potential impact of oil spills on early life stages of fish in vulnerable areas. Although the mechanisms inducing the observed bone abnormalities following crude oil exposure are poorly understood, PAHs are believed to be possible culprits (Carls et al., 1999;Incardona et al., 2004Incardona et al., , 2006Incardona et al., , 2013. Incardona et al. (2004) demonstrated that an individual PAH might cause skeletal deformities in zebrafish (Danio rerio) in vivo in a similar fashion as crude oil (Incardona et al., 2005(Incardona et al., , 2006. However, crude oil is composed of thousands of other compounds that may also contribute to its toxicity (Meador and Nahrgang, 2019). Many chemicals, such as PAHs, are known to induce bone abnormalities in fish larvae (Kingsford and Gray, 1996;Naruse et al., 2002;He et al., 2011;Corrales et al., 2014;Seemann et al., 2015). Environmental contaminants, such as dioxins, polychlorinated biphenyls (PCBs), and tributyltin can also affect the process of mineralization and cause a weakened skeleton (Hodgson et al., 2008;Carpi et al., 2009;Korkalainen et al., 2009;Herlin et al., 2010;Zhang et al., 2012). While effects on the process of bone metabolism might not lead to deformities directly, a weakened skeleton could make an animal more vulnerable to mechanical damage. In juvenile and adult fish with developed and mineralized vertebrae, there is a close relationship between mineral level and mechanical stiffness and strength (Fjelldal et al., 2006). Bone deformity has, therefore, commonly been considered as an early indicator of adverse effects in environmental studies. Although there are well-documented effects of oil on skeletal malformations in larval fish, there is a lack of studies following fish to later developmental stages.
The goal of this study was to examine the impact of dispersed crude oil on bone development in Atlantic haddock larvae. Atlantic haddock larvae were exposed to crude oil from 0 to 18 days post hatch (dph), and the larvae were analyzed for bone and cartilage integrity at 8 and 19 dph, and fried at 44 dph. Deformities were also assessed in surviving fish after 8 months. Endpoints included bone staining, x-ray of juvenile fish, and transcriptional analysis of genes involved in bone metabolism by RNA-seq. Based on a previous observation on cod larvae exposed to crude oil (Olsvik et al., 2011), we expected to see an upregulation in osteoclast-related genes and a down-regulation in osteoblasts-related genes.

Fish Husbandry and Exposure Experiment
Wild Atlantic haddocks were collected from spawning grounds outside Austevoll, Western Norway, and used as a brood stock at the Austevoll Research Station (IMR). The fish were kept in 7,000-L tanks resembling natural spawning conditions: temperature of 8 • C, light cycle of 12:12 light/dark, water flow of 32 L/h, and about 34.7 ± 0.2‰ salinity. A detailed description of exposure set up, hatching, feeding, and oil exposure regime has been provided previously by Sørhus et al. (2016). Briefly, fertilized eggs were collected from the tanks and placed in egg incubators at 7 • C. At 13 days post fertilization (dpf), embryos were transferred to a total of 16 green PE plastic tanks. Approximately 96,000 eggs were evenly distributed into the 16 tanks, with about 6,000 eggs in each tank. The number of larvae was estimated after hatching to be 5,975 ± 260 (mean ± SD). Starting at 4 days post hatch (dph), the larvae were fed with natural zooplankton, mainly copepod nauplii (Acartia longiremis), collected from the saltwater pond system Svartatjern (Sørhus et al., 2015). Due to the life history strategy of haddock, with females spawning a large number of eggs, high mortality is expected in all the treatment groups, including the control.
The experimental setup comprised three treatment groups that were exposed to mechanically dispersed crude oil, and one control with nominal concentrations: low (60 µg oil/L), pulse (0-600 µg oil/L), high (600 µg oil/L), and control (0 µg oil/L). There were four replicates per treatment. The pulse treatment group was designed to account for the vertical transport of larvae in and out of areas with high oil concentrations, which is relevant after an acute oil spill. Larvae in the pulse treatment group were exposed to high (600 µg oil/L), pulsating concentrations for 2.4 h every 24 h, and the concentration was decreased close to zero before a new exposure. Over time, they were exposed to the same concentration as the larvae in the low treatment group (60 µg oil/L). The larvae were exposed for 18 days, from 0 to 18 dph, and analyzed for gene expression at 1, 2, 4, 9, 14, and 18 dph, and for bone abnormalities at 8, 19, and 44 dph, and at 243 dph (8 months old). For preservation, the larvae were fixed in 4% paraformaldehyde (PFA) for 2 h, washed once in 70% ethyl alcohol (EtOH), and stored in 70% EtOH until further analysis. The juvenile fish (243 dph) were frozen for x-ray.
The larvae were exposed to a weathered blend of crude oil, containing oil from four different formations in the Heidrun oil field in the Norwegian Sea. These formations contain different types of oil, light paraffinic (0.83 g/ml) and heavy biodegraded (0.93 g/ml), and they are exported as a heavy crude oil (0.89 g/ml). The composition of this Heidrun blend oil is thought to be a representative for the types of oil that may be found in the area of Lofoten and Vesterålen (Sørhus et al., 2015). Due to the fast evaporation process that normally occurs after an offshore oil spill, the oil was artificially weathered by distillation at 250 • C to obtain realistic oil for exposure (Sørhus et al., 2015). The weathering procedure is described by Stiver and Mackay (1984). This process corresponds to 2-7 days on the sea surface at 10 • C, and the resulting oil has a decreased concentration of naphthalene and, thus, increased concentrations of PAHs with three rings (Sørhus et al., 2015). The principle of the exposure system and procedures for oil droplet generation are described in detail elsewhere Sørhus et al., 2015). The oil was pumped into the dispersion system using a highperformance liquid chromatography (HPLC) pump (LC-20AD liquid chromatography pump; Shimadzu, Kyoto, Japan) with a flow of 5 µl/min together with a flow of seawater of 180 ml/min. The system generates an oil dispersion with oil droplets in the low µm range with a nominal oil load of 26 mg/L (stock solution). The dose of exposure to the tanks was regulated with a parallel pipeline system with one line of flowing clean seawater and the other line containing a flow of the stock solution. The two pipelines were connected by a three-way magnetic valve, allowing water to be collected from both lines. Different dilutions were made by controlling the relative sampling time from the oil stock solution and clean water. The experimental setup consisted of three treatments and a control, each with four replicates: 60 µg oil/L (low dose), 600 µg oil/L (high dose), 600 µg oil/L for 2.4 h in a 24-h period (pulse dose), and no oil (control).

Analytical Chemistry
Water samples (1 L) were taken from each exposure tank at the beginning of the experiment, preserved by acidification (HCl, pH < 2), and stored at 4 • C in the dark until further handling. The PAHs were analyzed as described by Sørhus et al. (2015). Pooled tissue samples of the larvae (100 individuals, ∼200-250 mg wet weight) were examined for PAH accumulation by gas chromatography/tandem mass spectrometry (GC-MS/MS) as described by Sørhus et al. (2016). The body burden results are given as ng PAH/g wet weight. Newly hatched haddock larvae are a very lean fish containing mainly membrane lipid as phospholipids, and the lipid concentration is ∼0.8 % of the wet weight (Garcia et al., 2008;Sørensen et al., 2016).

Fluorescence Microscope Measurement
The accumulation of PAH metabolites in the bile of the larvae was recorded at 9 dph (after 9 days of exposure) by PAH-associated blue fluorescence staining using a Nikon AZ100 microscope (Nikon, Tokyo, Japan) equipped with a 4 ′ ,6-diamidino-2-phenylindole (DAPI) filter set (peak excitation 330-380 nm; emission > 425 nm) and a Nikon Intensilight C-HGFI fluorescence light source (Nikon, Tokyo, Japan). Pictures were taken using a 2× objective and 1× zoom and a QICAM monochrome camera (Teledyne Photometrics, Tucson, AZ, United States) (1,392 × 1,040 pixels, http://www.qimaging.com), and resulted in images with a resolution of 0.237 pixels/µm. The microscope UV-light, camera exposure time and gain factor settings were adjusted to the high exposure groups to avoid oversaturation in the pictures. These factors were then held constant (2 s exposure time, gain 5) for all the images. Emission was estimated using ImageJ (https://imagej.nih.gov/ij); the image of the gallbladder was manually encircled with an ImageJ region of interest (ROI) tool, and then the level of fluorescence was measured as the mean gray value (scale 0-255; 0 = black, 255 = white) of the ROI. If the gallbladder could not be visually detected because of lack of emission, the larvae were given a value of 0.

Bone and Cartilage Staining
To stain the bone and cartilage of Atlantic haddock larvae exposed to crude oil, slight modifications were made from the acid-free staining procedure described by Walker and Kimmel (2007). Optimal staining was achieved by adjustments in the concentration of magnesium chloride (MgCl 2 ), time of incubation in Alcian blue staining solution, time of incubation in Alizarin red staining solution, and strength of the bleaching solution. Ten larvae (19 dph) from each replicate tank were stained, yielding a total of 150 larvae (one tank was omitted because of increased larval mortality). In addition, a total of 150 juveniles (44 dph) were stained to examine if developmental skeletal abnormalities may appear later in life. The larvae were transferred to 6-well plates and stained with the replicate groups in separate wells. EtOH was removed before the addition of phosphate-buffered saline (PBS), which was replaced with fresh PBS, for 5-8 min (until the larvae sank to the bottom of the wells). The larvae were further washed in 100 mM Tris pH 7.5/100 mM MgCl 2 for 10 min. Cartilage staining was performed by adding to each well 4 ml of the 0.04% Alcian blue staining solution comprising 100 mM MgCl 2 . The larvae were incubated under gentle agitation conditions at room temperature for approximately 2 h. Then, the cartilage staining solution was removed, and the larvae were further washed in three separate steps in a decreasing ethanol concentration (80, 50, and 25%). Bleaching was performed in 0.1% H 2 O 2 /0% KOH for 2 h. To rinse the bleaching solution out of the wells, the larvae were washed once, and then twice 10-30 min in 25% glycerol/.1% KOH. Bone staining was performed in 4 ml 0.01% Alizarin Red under gentle agitation conditions at room temperature for 1 h. The larvae were further destained in 4 ml 50% glycerol/0.1% KOH under constant agitation conditions for 10 min. For further dehydration and preservation, the destaining solution was replaced with a stepwise increase in the concentration of glycerol, 50% glycerol/0.1% KOH for 24 h, 75% glycerol/0.1% for 24 h, and finally 100% glycerol.
The staining protocol was slightly different for 44 dph larvae. The 44 dph larvae were bleached in 0.1% H 2 O 2 /0.9% KOH for ∼5 h until sufficient pigmentation was removed, and 70% dH 2 O/30% saturated sodium borate was added to rinse the bleach out of the wells before the addition of 1% trypsin. Time of trypsinization ranged from 12 to 15 h until sufficient transparency was obtained. If necessary, the solution was replaced with fresh 1% trypsin, and 70% dH 2 O/30% saturated sodium borate was added again to remove the trypsin solution from the wells. Bones were stained in 0.01% Alizarin under gentle agitation conditions for 75 min at room temperature.

Analysis of Morphology, Growth, and Survival
Pictures of each larva were taken using an Olympus SZX Stereo Microscope (Olympus, Tokyo, Japan) with a Digital Slight Ds-Fi1 camera (Nikon, Tokyo, Japan). The larvae were randomized further before they were analyzed. The treatment groups were compared with control based on bone deformities and irregularities, and image analyses of selected parameters were conducted. The following parameters were measured using ImageJ: standard length (SL) of larvae, length from eye to upper jaw (EJ), and length of Meckel's cartilage (LM) (Figure 1). SL was measured along the notochord, from the posterior end and along the dorsal side until the transition to the head, and further in a straight line to the tip of the upper jaw (Skaalsvik et al., 2015). The measurement of SL at 44 dph was performed on fixed larvae. Furthermore, percent mineralization was divided by the SL of the larvae, allowing for comparison among the larvae regardless of size (Kvalheim, 2010). For the 44 dph larvae, an assessment of skeletal abnormalities, such as vertebral fusion, fusion of neural spines, and deformities in the lower and upper jaws, was also conducted. At two days (20 dph) and two months (76 dph) post exposure stop, larval survival was estimated based on multiple pictures obtained from the experimental tanks. A digital single-lens reflex (SLR) camera (Olympus E-3; Olympus, Tokyo, Japan) with a wide-angle lens (Olympus Zuiko Digital 11-22 mm; Olympus, Tokyo, Japan) was used. Images were taken from directly above. Four tanks (one from control, two from the low treatment, and one from the pulse treatment group) were eliminated because of bacterial growth introduced by startfeeding the larvae with living zooplankton. The four tanks were randomly distributed in the room, and the technical issue was not related to crude oil exposure.

Radiology
The Atlantic haddocks were radiographed with a portable Xray apparatus (Porta 100 HF; Eickemeyer Medizintechnikfür Tierärzte KG, Tuttlingen, Germany) and a 35 × 43-cm computed radiography (CR) image plate (Dürr Medical, Bietigheim-Bissingen, Germany) with 40 kV and 10 mAs and 70 cm distance. The image plate was scanned using CR 35 VETwin (Dürr Medical, Bietigheim-Bissingen, Germany), and the resulting image was converted into a TIFF file using Vet-Exam Plus Software version 4.14.0 (Dürr Medical, Bietigheim-Bissingen, Germany). The software program Adobe Photoshop CS2 was used for the evaluation of vertebral deformities according to the description of Witten and Huysseune (2009); vertebra number and deformity type were assessed for each deformed vertebra. Each fish was assessed for vertebral deformities, and the location and type of deformity were recorded.

RNA Isolation and RNA-Seq Analysis
Total RNA was isolated from frozen larvae using a Trizol reagent (Invitrogen, Waltham, MA, United States), followed by the DNase treatment step using a TURBO DNA-free kit (Life Technologies Corporation, Carlsbad, CA, United States). RNA quantity and quality were assessed using a Nanodrop spectrophotometer (NanoDrop Technologies, Wilmington, DE, United States) and the 2100 Bioanalyzer (Agilent Technologies, Sta. Clara, CA, United States). RNA sequencing was performed as described by Sørhus et al. (2017). Briefly, cDNA libraries were generated using the Illumina TruSeq RNA Sample Preparation Kit (Illumina, San Diego, CA, United States). Paired-end libraries were sequenced on an Illumina HiSeq2500 (Illumina, San Diego, CA, United States) platform. The raw data are available from the Sequence Read Archive (SRA) at the National Center for Biotechnology Information (NCBI) (Accession ID: PRJNA328092). In this study, six genes of interest for bone formation (sp7, bglap2, postn, and col101a1) and resorption (ctsk and mmp9) were extracted from the RNA-seq data set and studied in the larvae at 1, 2, 9, 14, and 18 dph.

Statistical Analysis
Statistica version 12 (2013, StatSoft Inc., Tulsa, OK, United States) and Graphpad Prism 9 (2021, GraphPad Software Inc., San Diego, CA, United States) were used for statistical analysis. One-way ANOVA, followed by Dunnett's post-hoc test, was performed for analyses of significant differences. Bartlett's test was performed to check for homogeneity of variance, while normality was checked by the D'Agostino-Pearson test. Data that violated Bartlett's test or the D'Agostino-Pearson normality test were log-transformed before testing for significance. If homogeneity in variance was not achieved, the non-parametric analysis Kruskal Wallis test was performed, followed by Dunn's post-hoc test. In both the parametric and non-parametric tests, the treatment groups were compared with the control group. A significance level of 0.05 was used for all the tests.
exposure water at 18 dph are shown in Supplementary Table 1.
The PAH profile in the exposure tanks was dominated by the C0-C3 naphthalene (∼60%), followed by tricyclic PAHs, which contributed ∼39% of the total PAHs. The PAH pattern in the larvae reflected that of the exposure water. PAH profiles of the larvae at 18 dph are shown in Supplementary Table 2. The uptake of totPAH in the tissue (totPAH body burden) was highest in the high concentration, and the uptake was larger in the low treatment (with continuous exposure) compared with the pulsed high concentration (2.4 h per day) exposure. Larval expression of cyp1a reflected the total PAH tissue concentrations in the treatment groups (Table 1).
Oil compound-associated fluorescence was applied as a measurement for biliary PAH metabolites at 9 dph.
Thirty individuals from each treatment group were scanned for PAH metabolites in the gallbladder. Supplementary Figure 1 shows representative images of four haddock larvae with biliary PAH metabolites in exposed individuals and measured light intensity in control, low, pulse, and high group individuals. There were a large number of individuals that had zero PAH-associated fluorescence in the control group (67%), and the low (54%) and pulse (58%) groups. In the low and pulse groups, 42 and 38%, respectively, of the individuals had medium or higher PAHassociated fluorescence (Supplementary Figure 1). The large majority (97%) of individuals in the high exposure group had high fluorescence, while one larva had intermediate levels.

Growth, Morphology, and Skeletal Development
Reduced SL was observed in all the exposure groups (Figure 2A). The high concentration showed reduced length at all the time points evaluated after 1 dph. At 44 dph, the larvae treated with the highest dose showed a 25% reduction in length compared with those in the control group. Larvae from the low and pulse treatment groups showed a 17.8 and 11.1% length reduction, respectively, at 44 dph. The length reduction percentage at 69, 134, and 243 dph between the larvae in the high treatment group and those in the control group were 6.7, 11, and 9.9%, respectively. At juvenile stages (69, 134, and 243 dph), only the high treatment group displayed reduced growth. At 134 dph, the SL of the high concentration larvae was significantly different from all the other treatments including control ( Figure 2B).
Due to high mortality rate at 44-dph, the number of individual fry analyzed per treatment group differed considerably (Supplementary Table 3). Upper jaw and MC lengths, normalized to total larvae length, at 8 and 19 dph are shown in Figures 2C,D, respectively. These parameters were not recorded in 44-dph fry. The Atlantic haddock larvae exposed to pulse and high concentrations showed significantly reduced upper jaw length at 8 dph ( Figure 2C, p < 0.001 and p < 0.0001, respectively), whereas significantly reduced upper jaw length was only observed for the high exposure concentration at 19 dph ( Figure 2C, p < 0.0001). For the length of MC, oil exposure had a less profound effect ( Figure 2D). Significantly reduced MC length was only observed in the larvae exposed to pulse concentration at 8 dph ( Figure 2D, p < 0.05).
Two days post exposure (at 20 dph), no significant differences in mortality were observed between the treatment groups and control (Supplementary Figure 2A). Two months later, at the juvenile stage (76 dph), the larvae exposed to high concentration had a survival rate of <3% as Length of juveniles at 69-243 dph after exposure termination. Normalized length of (C) upper jaw and of (D) Meckel's cartilage in larvae exposed at 8 and 19 dph. Asterisks indicate statistical differences from the control fish: p < 0.05 = *, p < 0.01 = ** , p < 0.001 = ***, p < 0.0001 = ****. Data are presented as mean ± SD.
FIGURE 3 | Images of (A) control Atlantic haddock larvae compared with (B) high treatment group larvae of approximately equal standard length at 19 dph. The shorter upper jaw of high treatment group larvae is clearly seen. Scale bars: 1 mm. FIGURE 4 | Control (A,C,E,G) and crude oil-exposed (B,D,F,H) 19-dph haddock larvae of equal standard-length stained with Alizarin red and Alcian blue. At 19 dph, ossification of the second neural arch has started in the larger larvae (C,D), 6.7 and 6.6 mm, respectively. Onset of the parasphenoid can be seen in all the larvae. When comparing larvae of similar length, the premaxillary was only observed in the control group (E,F), dorsal view. Branchial arches were also severely compressed in the exposed larvae (H), ventral view in comparison with the control larvae (G), ventral view. Scale bars: 1 mm.
opposed to the 10-14% rate of those in the other treatment groups. However, due to bacterial growth followed by variation in survival in some of the replicate tanks, this observation was not statistically different from control (Supplementary Figure 2B).
At 19 dph, the most prominent effect of crude oil exposure on larvae mineralization was seen in the high treatment group. Figure 3 shows two typical larvae of approximately equal standard length from the control ( Figure 3A) and high treatment (Figure 3B) groups.
Compared with the control larvae, crude oil exposure typically reduced upper jaw size of the high treatment group larvae. Furthermore, the larvae in the high treatment group often showed morphological changes in the overall craniofacial skeleton. The exposed larvae had a stumped/compressed craniofacial morphology and lacked the normal curvature observed in the control larvae (Figure 3). In addition to reduced SL and shorter upper jaw (Figures 3, 4), these larvae also had weaker mineralization of the maxillary in the upper jaw. Compared with most of the control larvae (95%), the premaxilla was not observed in any of the exposed animals from the high treatment group, also when compared with larvae of equal length (Figures 4A-F). Despite the lack of the premaxillary in the high group, other skeletal features, such as the 1st and 2nd neural arch and the parasphenoid, were mineralized in a similar fashion in all the experimental groups. The compressed craniofacial features in the high treatment group also led to severely compressed gill arches (ceratobranchialis), with thinner and overlapping ceratobranchial arches observed from the ventral side (Figures 4G,H). These abnormalities observed in the craniofacial skeleton in the high exposure larvae were not observed in the low and pulse-exposed larvae. The larvae with craniofacial compressions also had a markedly narrow jaw. The overall jaw width was also significantly reduced in all the exposed animals when normalized for standard length (Figure 5). Figure 6 shows the skeleton of a typical set of Atlantic haddock at 44 dph. Fry from all three groups exposed to crude oil showed reduced standard length, an overall trend toward incomplete ossification of the cartilage template in the craniofacial skeleton, as well as incomplete ossification of the notochord. No major skeletal deformities in the craniofacial skeleton were observed at 44 dph. Also, the Alizarin red staining did not disclose any abnormalities in the vertebral column, such as irregular, compressed, or joined vertebral bodies, in 44-dph haddock fry exposed to dispersed crude oil. However, an increased incidence of split neural arches was observed in the exposed animals (Table 2; Supplementary Figure 3). As seen in individuals from the pulse treatment group, which displayed a split, double neural arch on cordacentra numbers 13 and 16 (Figure 6D), minor irregularities were observed in some of the exposed fish. FIGURE 6 | Significantly reduced standard length of 44 dph juvenile Atlantic haddock exposed to crude oil, compared with the control (two control larvae shown). Standard length of haddock juveniles exposed to (C) low, (D) pulsating, and (E) high concentrations were significantly reduced compared with those in the control (A,B). They also had more incomplete mineralization of the vertebral column compared with those in the control group. Scale bars: 1 mm.

Long-term Effect on Vertebra Deformities
To further evaluate the phenotypic impact of early life crude oil exposure, surviving adult Atlantic haddocks were examined for bone deformities by X-ray. In total, 97 fish survived and lived to reach the juvenile stage. Table 3 shows the percentage of bone deformities in the surviving Atlantic haddocks. Of note, relatively few fish from the high treatment group survived this stage, explaining the lower number of examined fish from this group. Typical vertebra deformities in the high treatment group fish are shown in Figure 7, where vertebra compression combined with lordosis was the most prominent phenotype among the deformed individuals.

Transcription of Bone-Associated Genes
Mechanisms associated with bone formations are tightly regulated in animals. In cod larvae exposed to oil, we have previously seen the downregulation of genes encoding osteoblast-associated proteins and the upregulation of genes encoding osteoclast proteins. Four osteoblast-associated and two osteoclast-associated genes were, therefore, examined in detail in this study. Of the osteoblast-genes, three were lower expressed in the oil-exposed larvae in this study. The expression of Sp7, a zinc finger transcription factor expressed in osteoblasts and required for bone formation (osterix), was significantly lower in the high exposure group compared with the control group at 14 dph (ANOVA, with Dunnet's posthoc test, Figure 8A). The expression of Postn, which encodes a secreted extracellular matrix protein expressed in collagenrich tissues (periostin), was significantly lower in the high exposure group compared with the control group at 9, 14, and 18 dph (Figure 8C). The expression of Col10A1 (collagen, type X, alpha 1), a marker for maturing chondrocytes expressed during skeletal development, was significantly lower in the low exposure group compared with the control group at 2 and 9 dph, in the pulse exposure group compared with the control group at 9 dph, and in the high exposure group compared with the control group at 2, 9, 14, and 18 dph ( Figure 8D). In contrast, the expression of bglap2 (ostecalcin), which is among the most abundant bone protein genes and is involved in calcium binding and mineralization, was significantly higher in the low exposure larvae at 18 dph ( Figure 8B). Taken together and except the bglap2 result, oil exposure during early larval development seems to have an inhibitory effect on the transcription of bone formation genes in the Atlantic haddocks.
For the two studied osteoclast-associated genes, the general pattern was less clear. The expression of Ctsk (cathepsin K) (Figure 8E), encoding a lysosomal cysteine proteinase involved in bone remodeling and resorption, was significantly higher in the pulse treatment group compared to the control group at 2 dph, but otherwise did not show any differential expression during development after oil exposure, and neither did mmp9 (Figure 8F), which encodes matrix metallopeptidase 9, a zincdependent protein involved in the breakdown of the extracellular matrix. The transcriptional results from this study, thus, seem to verify a general pattern of oil-induced downregulation of bone formation genes in early life stages of Atlantic haddock, but the impact of oil on bone resorption genes was less pronounced. A typical split neural arch in 44 dph Atlantic haddock fry stained with Alizarin red and Alcian blue can be seen in Supplementary Figure 2. The fish were evaluated for deformities when their weight was about 65-87 g (8 months old).

Main Findings
This study examined the effect of dispersed crude oil on skeletal development in post-hatch larval stages of haddock and in surviving fry and juveniles. The most severe cases of malformations at 8 and 19 dph in the Atlantic haddock larvae exposed to weathered oil from 0 to 18 dph were craniofacial abnormalities. At 44 dph, most of the craniofacial abnormalities were gone. Larvae with major craniofacial abnormalities will have problems with feeding, accordingly, and most of them died before they reached 1.5 months of age. Both during the exposure and after being moved to clean water, the larvae treated with the high concentration showed reduced growth. A slight reduction in the increase in the prevalence of bone deformities was also observed in surviving juveniles exposed to the high concentration at 8 months of age. While high control mortality limits definitive conclusions, craniofacial abnormalities that are induced by oil exposure during the first weeks after hatch cause latent mortality and/or reduced growth due to reduced feeding ability. Further research is needed to confirm these results and translate effects reported under defined lab oil exposures investigated in toxicity tests to simulated field exposures of accidental releases of crude oil near spawning grounds of Atlantic haddock to quantify the ecological implications. Fluorescence detection is a highly sensitive method to document exposure to aromatic compounds like PAHs (Hornung et al., 2004(Hornung et al., , 2007 and oil (Gonzalez-Doncel et al., 2008;Bui et al., 2012;Sørhus et al., 2016) in exposed fish embryos and larvae. The detection is not quantitative; however, it gives a good picture of the distribution of the PAH-associated fluorescence, and when performing image analysis, the light intensity can be used to compare the uptake of fluorescence compounds among exposure groups. The PAH-associated fluorescence was only observed in the gallbladder, suggesting a high capacity for the haddock larvae to metabolize oil-related PAHs. Other oil exposure studies have also detected the accumulation of PAHs in the yolk sac of medaka (Gonzalez-Doncel et al., 2008) and zebrafish (Bui et al., 2012), but in 9-day-old haddock larvae, most of the yolk are utilized, and there is no indication of high fluorescence in the remaining yolk.

Bone Development
Most studies so far have assessed bone deformities in fish larvae after embryonic oil exposure. Numerous studies have reported jaw deformities in fish larvae after pre-hatch exposure to oil (Carls et al., 1999;Pollino and Holdway, 2002;Incardona et al., 2004;Nokame et al., 2008;He et al., 2011;de Soysa et al., 2012;Shi et al., 2012;Sørhus et al., 2016;Hansen et al., 2019;Perrichon et al., 2021). Incardona et al. (2004) observed a reduction in the size of craniofacial skeletal elements, including the jaw, in zebrafish after phenanthrene or dibenzothiophene exposure. These compounds were, after naphthalene, the most abundant PAHs in the Heidrun blend oil used in this study (Sørhus et al., 2016). Haddock embryos are highly susceptible to oil exposure because of their sticky eggshell (Sørhus et al., 2015;Sørensen et al., 2017). Even short exposures during the embryonic phase resulted in more severe functional and developmental defects (Sørhus et al., 2015(Sørhus et al., , 2016 compared with larval stages (Sørhus et al., 2016). Especially, exposure during cardiac development caused detrimental craniofacial and cardiac malformations . Similarly, the exposure of Atlantic halibut (Hippoglossus hippoglossus) during cardiac development led to more severe jaw abnormalities compared with exposure prior to cardiac development (Perrichon et al., 2021). Jaw deformities have also been reported in fish larvae after post-hatch exposure to crude oil (Ingvarsdottir et al., 2012), and in oil-treated fish larvae after pre-to post-hatch exposure (Solberg et al., 1982;Tilseth et al., 1984). Although exposure in late embryonic stages and early larval stages resulted in less abnormalities compared with the early embryonic phase, craniofacial deformities were frequent traits (Sørhus et al., 2016. Here, we document the acute and long-term effects on cartilage formation, bone mineralization, and growth after exposure post hatch. As expected, the craniofacial abnormalities were most pronounced in the high group with reduced upper jaw length, reduced jaw width, compressed craniofacial features, increased incidence of neural arch abnormalities, and lack of premaxilla in 19 dph larvae. We have previously established SL, as opposed to dph, as the most important factor determining the development and sequence of skeletal ossification (Saele et al., 2017). Thus, normalization according to SL, as applied in this study, takes into account different developmental stages at the time of measurement. While the general reduction in skeletal mineralization in the exposed groups was probably related to size and developmental stage, the lack of premaxilla in the high group was not. This indicates an unsynchronized delay of jaw development either through patterning or disruption of mineralization. While endochondral bone structures (lower jaw and gill arches etc.) are ossified cartilages, maxillary structures are dermal bones ossified at onset. The disruption of the ossification process is also supported by the observed reduced expression of the key osteoblast specific transcription factor sp7 involved in osteoblast differentiation and activity (Lian et al., 2006).
Subtle defects on jaw structures might impact the possibility of long-term survival in the wild. A compressed and narrower jaw, in addition to lack of premaxilla as observed in the crude oil-exposed larvae, would also impact the selection of prey size and feeding success. Exposure during early larval stages is less lethal compared with embryonic exposure. The lower sensitivity of larvae compared with embryos reflects their increased ability to metabolize unwanted substances. Nevertheless, the severity of jaw deformities reported in haddock larvae after pre-hatch exposure in this study documents that fish larvae are vulnerable to oil throughout the early life stages. Tilseth et al. (1984) registered reduced food intake along with jaw malformations and reduced SL of Atlantic cod exposed to 0.25 ppm (250 µg oil/L) of water-soluble fractions of crude oil from 10 to 14 dph FIGURE 8 | Transcriptional levels of six genes involved in (A-D) bone formation and (E,F) bone resorption in Atlantic haddock larvae at 1, 2, 7, 9, 14, 18 dph. Data are presented as mean ± SD. RPKM, reads per kilobase million. Asterisks indicate statistical difference from the control fish: p < 0.05 = *, p < 0.01 = **, p < 0.001 = ***, p < 0.0001 = ****.
(21 days). Carls et al. (1999) found that reduced lower jaw size affected swimming ability in Pacific herring, which, again, will affect the ability to capture prey. After spinal defects, reduced lower jaw size was the second most important factor influencing swimming ability. Reduced food intake due to a deformed jaw can have long-term adverse effects. For example, embryonic exposure to oil resulted in altered lipid content and growth impairment in juvenile polar cod (Boreogadus saida) (Laurel et al., 2019). In Atlantic cod, we have seen previously that the exposure of first-feeding larvae (9 to 13 dph) to dispersed oil has resulted in similar transcriptional and metabolic responses as food deprivation (Hansen et al., 2016). Thus, the reduction in SL observed in this study may be a result of reduced food intake due to jaw deformities. However, other factors may also affect SL, such as changes in metabolic rate, or morphological changes, such as yolk sac edema and heart defects, which might affect swimming ability and inhibit prey capture (Carls, 1987;Carls et al., 1999;Incardona et al., 2004Incardona et al., , 2014. Carls et al. (1999) observed that Pacific herring larvae exposed to weathered oil were generally smaller and more immature compared with unexposed larvae. No effect of oil exposure was observed on bone density and bone repartition in adult European bass (Dicentrarchus labrax), but bone mineralization of the vertebrae decreased in a reversible way in exposed fish (Danion et al., 2011). The authors suggested that vertebral bone mineralization could be used as an early stress indicator of PAH pollution in sea bass. Vertebral malformation was also observed in juvenile haddock after a 2-month dietary exposure to produced water, oil, and heavy PAHs (Meier et al., 2020).
Larvae with reduced growth may recover and grow normally, but not many studies have examined recovery after oil exposure. Decreased growth rate due to abnormal jaw development may affect larval survival and, ultimately, population dynamics (Ingvarsdottir et al., 2012). In this study, the larvae were afforded the opportunity to recover from 18 to 44 dph (26 days). At 44 dph, all the exposed larvae had a significantly shorter SL than the larvae in the control group. In addition, at 44 dph, the vertebral column and craniofacial skeleton of haddocks in all the treatment groups were not fully mineralized. This contrasted with the vertebral column of most haddock in the control group, which was very close to or fully mineralized at 44 dph. Considering the differences in SL, incomplete ossification in the three exposure groups was most likely a result of restricted growth rate rather than a direct effect on processes involved in mineralization and bone development. A malformed notochord will result in an abnormal vertebral column (Grotmol et al., 2005). However, no consistent malformations in the notochord were observed in the exposed haddock larvae at 19 dph, and no malformations in the vertebral column were observed at 44 dph other than split spinal arches. Still, a higher prevalence of spinal deformities was observed in the haddocks after 8 months, possibly a longterm effect of slower mineralization in the individuals exposed to crude oil.
In adult fish, skeletal deformities can impact physiological processes such as swimming, reproduction, growth, resistance to disease, and susceptibility to predation (Eissa et al., 2021). In this study, lordosis was observed among the surviving fish from the high exposure group. This type of deformity has been observed earlier, both in wild haddocks (Jawad et al., 2018) and cultured haddocks fed with a phosphorous-deficient diet (Roy et al., 2002). Lordosis development in gadoids starts with vertebra compression, and then more deformed vertebrae aggravate the deformity as the curvature increases, where some of the adjacent compressed vertebrae ultimately fuse (Fjelldal et al., 2009). The currently observed lordosis is mild and involves compressed vertebrae, suggesting an early stage of lordosis. In Atlantic salmon, low vertebral bone mineral content reduces the mechanical stiffness and strength of the vertebra (Fjelldal et al., 2006), and results in compression of the vertebra (Fjelldal et al., 2007). Whether oil exposure during early life stages in fish has long-term negative effects on bone mineralization remains to be elucidated. However, in a previous study, Atlantic cod given a bad nutritional start displayed an increased incidence of vertebral malformations later that were not observed at earlier time points during mineralization of the vertebra (Saele et al., 2017), similar to this study. Other studies have also shown a lack of compensatory growth in cod given a nutritionally insufficient diet (Øie et al., 2015).
Reduced survival, asymmetric head shape, and SL were observed in Atlantic herring 8 weeks after a 12-day larval exposure (15-750 µg oil/L) (Ingvarsdottir et al., 2012). In contrast, the abnormal upper jaw and the overall less organized and defined structure of the craniofacial skeleton observed in the high treatment group at 19 dph were not observed at 44 dph in this study. Furthermore, the split neural arches observed at 44 dph were not visible in the 8-month-old haddocks following x-ray analysis. There may be three possible reasons for this: (i) larvae with the most severe jaw deformities died, (ii) they had arrested growth of the jaw, which they managed to recover from, but the growth continued, (iii) they were more immature compared with the unexposed control (retarded development). The latter reason is implausible because development was adjusted for length. When comparing the 19-dph control larvae with the high treatment larvae of approximately equal length, the shorter upper jaw was still present. Since bone growth is a continuous remodeling process throughout development, larvae with reduced jaws might later recover their ability to feed (Carls et al., 1999), so reason ii is possible. Most likely, larvae with the most severe jaw deformities died; however, further data are needed to confirm the significance and magnitude of long-term latent effects on juvenile survival, growth, and bone-related developmental effects. The Exxon Valdez oil spill overlapped with the Pacific herring spawning period, and field studies that examined the short-term effect on early life stages of herring larvae showed abnormalities such as reduced growth, low larval weight, and craniofacial and jaw deformities (Pollino and Holdway, 2002). However, the level of deformities observed in larval Pacific herring populations decreased over time, suggesting premature death in the most severe cases or alternatively potential reversibility of abnormalities (Carls et al., 1999). Decreased feeding capacity and survival, over time, have also been observed in zebrafish with craniofacial abnormalities after TCDD exposure during the embryonic stage (Chollett et al., 2014). Taken together, the result presented here is in line with the field observation of Pacific herring after the Exxon Valdez oil spill and suggests that fish larvae with severe jaw deformities will not survive during recovery.
The abnormalities observed in the high treatment group larvae were more pronounced than in the low and pulse treatment group larvae. Over time, the larvae in the low and pulse treatment groups were exposed to the same concentrations (60 µg oil/L), but the peak concentration in the pulse dose treatment was 10 times higher (600 µg oil/L). However, at 19 dph there were no significant differences in normalized upper jaw and Meckel's cartilage length between the larvae in the low and pulse treatment groups and those in the control group (although the low and pulse treatment group larvae were shorter). This result indicates that 18-day chronic pulse exposure with a mean measured TPAH exposure of 3.2 ppb exhibits similar effects when compared with 18-day continuous oil exposure corresponding to 0.7 ppb (low concentration, see Table 1). A pulsed exposure scenario may be more environmentally relevant given haddock larval drift in and out of oil-contaminated water plumes following a spill.

Mechanisms of Action
Toxicants affecting bones have been classified into four categories according to the mode of action, activation/inhibition of membrane and/or nuclear receptors, alteration of redox condition, mimicking of bone constituents, and unknown pathways (Fernández et al., 2018). In terms of oil as an osteotoxicant, researchers have mainly focused on highmolecular PAHs and their ability to activate the aryl hydrocarbon receptor (AhR) and impact on downstream target genes, such as cyp1a, encoding enzymes catalyzing the metabolism of these compounds. For some PAHs, embryotoxicity is independent of CYP1A induction but is still mediated by AhR (Hodson, 2017). However, despite a lot of recent research studies, the underlying molecular mechanisms of PAH-mediated effects on skeletal development and bone homeostasis are not fully understood. The maintenance of bone size, shape and integrity is a complex process at the molecular, cellular, and tissue levels. To maintain bone homeostasis, a perfect balance between bone deposition (osteogenesis) and bone resorption (osteoclastogenesis) is needed. Chondrocytes and osteoblasts construct and shape the skeleton for maximal resilience, while osteoclast maintains mineral homeostasis by resorbing the extensive surface of mainly cancellous bones (Zaidi, 2007;Witten and Hall, 2015). In our previous study on Atlantic cod larvae exposed to weathered dispersed crude oil (up to 10.4 µg PAH/L) between 9 and 13 dph, we showed that elements in the oil inhibit bone formation and enhance bone resorption mechanisms (Olsvik et al., 2011). In cod larvae exposed to oil, genes associated with osteoblasts, such as sp7, col10a1, col1a1, col1a2, and bglap2, were downregulated, while genes associated with osteoclasts, such as ctsk, mmp9, and acp5, were upregulated. Accordingly, we show that the osteoblast-associated genes sp7, postn, and col10a1 were downregulated in oil-exposed haddock larvae. In this study, ctsk at 2 dph was the only osteoclastassociated gene with significantly elevated expression. Hansen et al. (2021), however, did not observe any clear patterns of differential expression of osteoblast (sp7, bglap2)-and osteoclast (ctsk)-associated genes in Atlantic cod embryos exposed to watersoluble fractions of fuel oils from 4 to 8 dpf at concentrations ranging from 4.1 to 128.3 µg totPAH/L. Thus, it appears that the expression of genes encoding proteins involved in the ossification process is more distinctly affected by oil after larval exposure than after embryonic exposure. Although we did not see a clear pattern of dose-dependent upregulation of osteoclastassociated genes in the haddocks, activated bone resorption will contribute to skeletal deformities. The reduced expression of osteoblast-associated genes and increased expression of osteoclast-associated genes suggest that components in the oil caused a shift in the balance between formation and degradation favoring the latter.
Numerous studies have documented that various PAHs can inhibit mechanisms associated with bone formation and enhance bone resorption in fish. In Japanese medaka (Oryzias latipes), it has been proposed that PAHs, such as benzo(a)pyrene, typically cause perturbed chondroblast and osteoblast differentiation due to differential expression of key genes regulating osteoblast-mediated bone formation, resulting in a defective notochord sheath repair, and rendering the vertebral column more vulnerable to physical damage (Yamaguchi et al., 2020). In rockfish (Sebastiscus marmoratus) exposed to benzo(a)pyrene during embryogenesis, He et al. (2011) suggested that PAHs act by perturbing the proliferation of chondrocytes and impairing the expression of sonic hedgehog pathway genes, leading to disturbed craniofacial formation. Benzo(a)pyrene reduced osteoclast and osteoblast activity in scales of zebrafish and goldfish (Carassius auratus) after ex-vivo exposure (Torvanger et al., 2018). In rockfish exposed to three PAHs for 6 days after hatch, Li et al. (2011) observed inhibited Na + /K + -ATPase and Ca 2+ -ATPase activity in a dose-dependent manner in larvae with craniofacial defects. Decreased expression of genes in the calcium signaling pathway has also been shown in haddock larvae exposed to oil . Disturbance of calcium homeostasis can impact the proliferation and differentiation of chondrocytes, as shown in rockfish embryos exposed to tributyltin (Zhang et al., 2012). Furthermore, many studies focusing on cigarette smoking and human stem cells developing into osteoblast have shown that PAHs impact the balance between osteogenesis and osteoclastogenesis. For example, Chen et al. (2020) reported that PAH exposure is associated with increased bone resorption among the general adult population in the United States. Similar results have been reported in human mesenchymal stem cells, cells that can differentiate into osteoblasts, exposed to benzo(a)pyrene, in which the PAH exposure suppressed early and late osteogenic differentiation and down-regulated the expression of runt-related transcription factor 2 (Runx2), osteocalcin and osteopontin . The authors suggested that impaired BMP/Smad pathways through AhR regulating BMPRII and Hey1 may be the underlying mechanism for benzo(a)pyrene inhibiting BMP2induced osteogenic differentiation. There is also evidence that benzo(a)pyrene enhances vitamin D3 catabolism, thereby potentially interrupting bone homeostasis (Matsunawa et al., 2009). Research has also shown that endocrinedisrupting chemicals found in oils, such as alkylphenols, have an impact on the function of osteoclasts by nonestrogenic effects (Hagiwara et al., 2008;Prada et al., 2020). Components in oils could also have an impact on bone homeostasis because of oxidative stress (Zancan et al., 2015). Combining these findings will illustrate clearly that PAHs affect mechanisms associated with osteoblasts and osteoclasts in fish larvae.

CONCLUSIONS
This study shows that haddocks exposed to oil for 18 days at the early larval stage had reduced growth and developed skeletal malformations. After a recovery period of 27 days, the exposed larvae did not show abnormality in the upper jaw and less organized and defined structure of the craniofacial skeleton, suggesting that the larvae with jaw deformities most likely failed to feed properly and died. A slight reduction in growth and a higher percentage of vertebra deformities were observed in surviving juveniles from the high treatment group 8 months after the end of the exposure. Reduced mineralization and inhibition of osteoblastcontrolled mechanisms linked to bone formation could, in part, be responsible to the observed negative impacts on skeletal development.

DATA AVAILABILITY STATEMENT
The raw data are available from the Sequence Read Archive (SRA) at NCBI (Accession ID: PRJNA328092).

ETHICS STATEMENT
The animal experiment described in this study was 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.

AUTHOR CONTRIBUTIONS
KL and PO designed the research questions. PO wrote the first draft of the manuscript. KL, IT, ES, SM, and PF contributed to the writing of the manuscript. KL, PO, IT, ES, SM, MT, AT, LS, and PF analyzed the data. KL, PO, IT, ES, and IG conducted sampling of the fish. SM, ES, and ØK conducted the main experiment. SM led the main project, designed the exposure study, and obtained funding. All authors contributed to manuscript revision, read, and approved the submitted version.

FUNDING
The research reported in this study was financially supported by the Norwegian Research Council (Grant Nos: 234367 and 267820), and by the VISTA Foundation (project no: 6161).