Cardiac remodeling in response to embryonic crude oil exposure involves unconventional NKX family members and innate immunity genes

The developing heart is the first organ to become functional, but it does so at a point very early in its morphogenesis. As a muscular pump driven by physiologically complex excitable cells (modified cardiomyocytes), it is highly susceptible to disruption by a variety of insults, in particular exogenous chemicals. In humans, for example, heart-related toxicities are the most common adverse reaction to pharmaceuticals (Wilke et al., 2007). Because of the intimate links between form and function in the developing heart (Andrés-Delgado and Mercader, 2016; Miquerol and Kelly, 2013), even very low-level cardiotoxicity during embryogenesis can have serious and long-lasting outcomes. This poses a real-world environmental health problem for developing fish, most of which undergo embryogenesis externally without the metabolic protection of a maternal circulation. The exquisite sensitivity of the developing fish heart to chemical insults has been demonstrated through studies of oil spills in the vicinity of spawning fish populations (Incardona and Scholz, 2017, 2018). Certain polycyclic aromatic hydrocarbons (PAHs) from petroleum disrupt K⁺ and Ca²⁺ ion channel conductances controlling cardiac action potentials and excitation–contraction coupling in cardiomyocytes (Brette et al., 2014, 2017). During cardiac organogenesis, acute exposure to PAHs leads to cardiac arrhythmia and reduced contractility, which in turn lead to secondary defects in heart looping and growth of the ventricular chamber (Incardona et al., 2009, 2004, 2014; Jung et al., 2013; Sørhus et al., 2016). At higher exposure levels, severe heart malformation is lethal to developing larvae. By contrast, lower exposure concentrations lead to grossly normal early heart development (Hicken et al., 2011; Incardona et al., 2013). However, physiological studies indicate that surviving fish have reduced swimming performance later in life, coupled with subtle defects in cardiac anatomy and histology that are representative of pathological remodeling (Hicken et al., 2011; Incardona et al., 2015). These defects include abnormal ventricular shape, reduced compact myocardium and an abnormal hypertrophic–hypercellular response in spongy myocardium.

The anatomical and histological changes in the hearts of juvenile fish exposed to oil as embryos link early impacts on cardiac morphogenesis to pathological heart remodeling at much later stages in life. While these effects on heart development are likely to underlie population-level impacts of oil spills (Heintz, 2007; Incardona et al., 2015; Peterson, 2001), they are poorly understood at the molecular level. Largely because of the use of zebrafish as a model for human heart development, early fish heart development is understood at a relatively high level of detail (Sidhwani and Yelon, 2019; Staudt and Stainier, 2012). Much less is known about post-embryonic development (Singleman and Holtzman, 2012). However, several studies indicate that effects on both early signaling pathways and biomechanical forces can influence ventricular structure and function at much later stages. The initial linear heart tube in organogenesis-stage fish embryos must undergo looping of the atrial and ventricular chambers, followed by ventricular ballooning (outgrowth) and subsequent trabeculation to form spongy myocardium. For high-performance species (including zebrafish and salmonids), compact myocardium also forms during the late larval to juvenile period (Gupta and Poss, 2012). Developmental genetics studies in zebrafish demonstrated that normal contractility in the very early heart is required for the trabeculation process (Peshkovsky et al., 2011; Staudt et al., 2014), and embryos with defective trabeculation grow into adults with pathological ventricular hypertrophy (Abdul-Wajid et al., 2018). These observations made in zebrafish mutants are consistent with our observations in oil-exposed zebrafish (Hicken et al., 2011), pink salmon (Oncorhynchus gorbuscha) and Pacific herring (Clupea pallasi) (Incardona et al., 2015), and provide the basis for a mechanistic link between early crude oil cardiotoxicity and latent impacts in juveniles and adults.

Underlying molecular responses associated with more severe heart malformation in oil-exposed embryos and larvae of several species have been examined both through hypothesis-driven gene selection and global transcriptome profiling using RNA sequencing (RNAseq) (Edmunds et al., 2015; Sørhus et al., 2017, 2016; Xu et al., 2016). Here, we used RNAseq in juvenile pink salmon, a key indigenous species affected by the 1989 Exxon Valdez oil spill, to identify molecular changes associated with previously described ventricular defects arising from embryonic exposure to Alaskan crude oil (Incardona et al., 2015). This study has relevance for both environmental health and human medicine. Our findings have important implications for the science of oil spill damage assessment and recovery monitoring (Scholz and Incardona, 2015), as we have identified potential new biomarkers of long-term, delayed injury from oil exposure during fish early life stages. Moreover, because of the unique plasticity of the fish heart for both remodeling and regeneration (Foglia and Poss, 2016; Gamperl and Farrell, 2004), which are linked, our findings could help identify the mechanisms that limit this capacity in the damaged human heart.

A total of 37 hearts from pink salmon juveniles, O. gorbuscha (Walbaum 1792), were transcriptionally profiled following exposure to four different concentrations of crude oil-derived PAHs (here termed ‘trace’, ‘low’, ‘medium’ and ‘high’; detailed exposure concentrations and tissue doses have been published previously; Incardona et al., 2015). From these experiments, eight hearts were selected to maximize statistical differences in ventricular shape from each of the exposure doses and control, with the exception of the trace exposure for which only 5 animals were available for cardiac transcriptome analysis. Total RNA was extracted from each ventricle using TRIZOL^® reagent as recommended by the manufacturer (Thermo Fisher Scientific). RNA quantity was measured using a Qubit assay (Thermo Fisher Scientific). RNA integrity was assessed using a Bioanalyzer assay (Agilent, Santa Clara, CA, USA). Total RNA from a subset of the samples representing each tissue and treatment was pooled in preparation for construction of a de novo assembly. In preparation for sequencing, the pooled sample along with the individual samples (38 in total) were processed and used to generate amplified libraries using the TruSeq RNA library prep kit v2 protocol (38 libraries total). All 38 libraries were sequenced using the Illumina HiSeq 2000 platform and are deposited in the NCBI Gene Expression Omnibus (Edgar et al., 2002), and are accessible through GEO Series accession number GSE102654.

A juvenile pink salmon cardiac transcriptome was assembled de novo from 101 bp paired-end sequences from a pooled cardiac RNA sample. Sequence reads were parsed for quality metrics using Trimmomatic software (Bolger et al., 2014). Reads not passing the quality thresholds were removed, thereby excluding approximately 17 million of 296 million reads from the pooled cardiac library. Remaining reads were inputted into the CLC Genomics Workbench (Qiagen Bioinformatics, Redwood City, CA, USA) de novo assembly algorithm, which uses de Bruijn graph theory to break the reads into shorter sub-sequences.

A total of 222,598 contigs were assembled from the pooled pink salmon cardiac RNA (see Dataset 1, tab 1). Lengths ranged from 64 to 23,422 bp, with N₂₅, N₅₀ and N₇₅ values of 2690, 1132 and 373, respectively. The original sequence reads were mapped back to the assembly using Bowtie2 (Langmead and Salzberg, 2012) to determine accuracy and exclude sequences not correctly assembled (see Dataset 1, tab 2). A total of 244 million reads could be reliably mapped back to the assembly. Read pairs mapping discordantly were identified when mate pairs aligned outside the expected genomic distance range of 200–500 bp. These reads were excluded, leading to a total of 194 million reads mapping concordantly for the de novo assembly.

Sequencing reads (50 bp single-end) from the 37 individual hearts comprising the different experimental treatments were preprocessed for quality control using the same procedure as for the assembly. Following quality filtering, the remaining reads from each sample were aligned to the de novo assembled transcriptome using the Bowtie2 algorithm (Langmead and Salzberg, 2012). All samples aligned similarly to the transcriptome with a total assembly percentage between 80% and 83% of the reads (see Dataset 1, tab 3). Uniquely aligned reads (60–65% of the total) were used for quantitative gene expression analysis. Read count values were obtained using Samtools (Li et al., 2009) and were not normalized as the downstream statistical analysis package used a normalization strategy unique to the analysis.

Identification of differentially expressed genes (DEGs) was performed using the R statistical program (DESeq2 package), within the Bio conductor software project (Love et al., 2014). DESeq2 estimates variance-mean dependence in count data from high-throughput sequencing assays and tests for differential expression based on a model using the negative binomial distribution. To determine unbiased potential relationships, principal components analysis (PCA) and hierarchical clustering were used initially to assess large-scale variation among individual samples. The PCA first component showed high variation from five of the eight low-dose samples compared with all other samples and treatments. The second component showed a loose relationship with exposure to increasing oil concentration. An explanation for this relationship was not obvious, but could reflect an uncontrolled tank effect during rearing for the low-dose treatment group. Pairwise comparisons for each of the treatments were performed using the parametric Wald test with Benjamini–Hochberg adjusted P-values (see Dataset 1, tab 4). Consistent with the initial PCA, the pairwise comparisons with the greatest number of DEGs involved the low-dose treatment. There were no DEGs common among the high- and low-dose treatments, and only a single gene was shared between the medium and high doses. Repeating the PCA without the low-dose samples showed no discernible patterns at the treatment level. In this case, pairwise comparisons (Wald test with Benjamini–Hochberg adjusted P-values) showed a dose-dependent pattern of smaller numbers of differentially expressed genes (see Dataset 1, tab 5). This final list of DEGs was annotated, with ambiguities clarified by additional manual sequence analysis (see below). The remainder of the transcriptome (i.e. genes not differentially expressed) was annotated at a later date. DEGs identified from DESeq2 were aligned to the Swiss-Prot/UniProt protein database using the BLASTx algorithm and analyzed with the pathway analysis software Ingenuity Pathway Analysis (IPA) (https://www.qiagenbioinformatics.com/products/ingenuity-pathway-analysis/).

Plasmid templates for Danio rerio nkx2.3, nkx2.5 and nkx3.3 were cloned with a TOPO^® TA Cloning^® Kit (Thermo Fisher Scientific) from the purified PCR amplicons of the coding region for each gene. Capped/poly(A)-tailed mRNAs were synthesized from linearized plasmids utilizing a mMESSAGE mMACHINE^® T7 Ultra Kit (Thermo Fisher Scientific). Purified mRNAs were diluted in sterile distilled water at concentrations enabling microinjection of 10 or 100 pg in a 4 nl volume into early cleavage stage embryos of AB strain zebrafish using a Harvard Apparatus PL100 injection system (Holliston, MA, USA). The zebrafish broodstock was maintained according to the policies of the US Department of Commerce and the Public Health Service, conforming to the standards of the National Academy of Sciences' Guide for the Care and Use of Laboratory Animals (National Research Council, 2011; Linbo, 2009). At mid-gastrulation, viable embryos were transferred to 100 mm Petri dishes in 40 ml zebrafish system water and incubated at 28.5°C. For experiment 1, 30–40 embryos for each treatment were incubated in a single replicate dish. For experiment 2, 60–90 embryos for each treatment were divided equally among three replicate plates for incubation. At 48 h post-fertilization (hpf), larvae were anesthetized with MS-222 (Sigma-Aldrich) and mounted in 2% methylcellulose for digital videomicroscopy using Nikon SMZ800 stereoscopes, Unibrain Fire-i780c 1394 cameras (unibrain.com) and BTV Pro application (Bensoftware.com). All larvae were imaged for experiment 1; for experiment 2, 10 larvae from each replicate dish were randomly selected for imaging for a total of 30 per treatment. The score for cardiac looping was based on orientation of the atrioventricular canal in left-lateral videos. The pericardial area was measured using the freehand trace tool in ImageJ (https://imagej.nih.gov/ij/).

PCR amplicons were generated for use as DNA templates for digoxigenin-labeled cRNA in situ probes. PCR primers were designed using Primer3 (http://bioinfo.ut.ee/primer3/) and were based on the juvenile pink salmon cardiac transcriptome sequences for ifit, ltbp1, nkx2.3, nkx3.3 and rsad2, or published sequences for zebrafish nkx genes (see Table S2.) For pink salmon genes, each probe DNA sequence was used as a query in a NCBI-BLAST search (https://blast.ncbi.nlm.nih.gov/Blast.cgi?CMD=Web&PAGE_TYPE=BlastDocs&DOC_TYPE=Download) of the pink salmon cardiac transcriptome, and any similar non-target mRNAs that were found had a predicted melting temperature (T_m) at least 12°C lower than that of the target mRNA. Probe lengths for ifit, ltbp1, nkx2.3, rsad2 and zax were 265, 197, 114, 267 and 266 nucleotides, respectively. To generate templates for gene-specific anti-sense probes, the forward PCR primer was matched with a reverse PCR primer that had a T7 RNA polymerase promoter sequence added to the 5′ end. For each gene, the template for the non-specific, sense probe was generated using a forward primer that had an SP6 RNA polymerase promoter sequence added to the 5′ end, matched with the reverse primer. All templates were gel purified, and the DNA sequence was confirmed. Templates were used to generate digoxigenin-labeled cRNA in situ probes as per the manufacturer's instructions for the DIG RNA Labeling Kit (Roche #11175025910).

AB strain zebrafish embryos were fixed at 48 hpf in 4% paraformaldehyde and whole-mount in situ hybridization was carried out as described elsewhere (Edmunds et al., 2015). Wild juvenile pink salmon (fork length 142–176 mm, mass 29–60 g) from Icy Strait, AK, USA, were provided by NOAA Alaska Fisheries Science Center through the July 2015 annual Southeast Coastal Monitoring survey (capture approved by Alaska Department of Fish and Game Scientific Research permit #2015-6). Juvenile pink salmon hearts from wild fish or experimental fish that were oil-exposed as embryos were fixed overnight at 4°C in 10% neutral buffered formalin, then stored in 70% ethanol at 4°C before being embedded in paraffin, sectioned at 5 µm thickness, and placed onto positively charged slides. Slides were de-paraffinized in xylene substitute, hydrated through an alcohol series, and rinsed with dH₂O. In a room temperature humidity chamber with 5× saline sodium citrate (SSC), slides were treated for 2 min with phosphate buffered saline (PBS), 5 min with 0.05 mol l⁻¹ Tris-HCl, pH 7.5, 20 min with 5 µg ml⁻¹ proteinase K (VWR #E195) in Tris, 10 min with 0.1 mol l⁻¹ triethanolamine buffer, pH 8.0 (TEA), 10 min with 0.25% acetic anhydride in TEA, followed by three 5 min treatments with PBS. Slides were then incubated for 2 h at room temperature with hybridization buffer (50% deionized formamide, 5× SSC, 5× Denhardt's solution, 250 µg ml⁻¹ yeast tRNA, 1 mg ml⁻¹ calf thymus DNA and 10% dextran sulfate), followed by hybridization overnight with 0.3 ng µl⁻¹ in situ probe diluted in hybridization buffer (except for the probe for rsad2, which was diluted to 0.15 ng µl⁻¹). Hybridization was performed at the highest temperature at which a significant signal could be obtained, and ranged from 55 to 70°C, depending on the probe. Based on the calculated T_m for each probe, a high-stringency hybridization temperature was optimum for all probes except for the ifit1 probe hybridized to experimental tissue, which was medium stringency. Post-hybridization washes included 30 min at room temperature with 5× SSC, 15 min at the hybridization temperature with 30% deionized formamide in 5× SSC, and two 15 min treatments at the hybridization temperature with 0.2× SSC. Next, slides were treated at room temperature for 5 min with Buffer B1 (0.1 mol l⁻¹ Tris-HCl plus 0.15 mol l⁻¹ NaCl, pH 7.5) and 1 h with Buffer B2 (Buffer B1 plus 1% heat-inactivated goat serum). Overnight incubation at 4°C was performed with anti-digoxigenin-AP, Fab fragments (Roche #11 093 274 910) diluted 1:300 in Buffer B2, followed by room temperature washes: three 5 min washes with Buffer B1 followed by 5 min wash with Buffer B3 (0.1 mol l⁻¹ Tris-HCl, 0.1 mol l⁻¹ NaCl, 0.05 mol l⁻¹ MgCl₂, pH 9.5). Hybridized probe was detected by a coloration reaction performed at room temperature or 4°C using buffer containing 100 mg ml⁻¹ NBT (Roche #11 383 213 001), 50 mg ml⁻¹ BCIP (Roche #11 383 221 001) and 50 mg ml⁻¹ tetramisole (Sigma-Aldrich #L9756) in Buffer B3. Color was allowed to develop for 2 h to 5 days, depending on the probe, and was stopped by adding Buffer TEN (10 mmol l⁻¹ Tris-HCl, 1 mmol l⁻¹ EDTA, 0.9% NaCl, pH 8.0) for 10 min, followed by a 30 s treatment with 100% ethanol, 10 min with Buffer TEN, and mounting in 50% glycerol in Buffer TEN. To unequivocally identify cardiomyocytes following in situ hybridization, immunohistochemistry was performed with MF20, a mouse monoclonal IgG_2b antibody specific for myosin heavy chain (Developmental Studies Hybridoma Bank, University of Iowa). MF20 was used at 3.5 µg ml⁻¹ with Vectastain Elite ABC Kit and DAB Peroxidase Substrate Kit (Vector Laboratories, #PK-6102 and #SK-4100, respectively) for detection, following the manufacturer’s protocol. Antigen unmasking was performed by heating tissue sections at 94–99°C in buffer containing 10 mmol l⁻¹ sodium citrate and 0.05% Tween 20, pH 6.0, for 10 min, followed by cooling at room temperature for 30 min.

Gene-specific primers were designed using Primer3 and were based on the juvenile pink salmon cardiac transcriptome sequences for ifit, ltbp1, nkx2.3, nkx3.3, rsad2 and wdtc1 (see Table S2), except for ef1α (Luckenbach et al., 2010). Where possible, the same primers were used for target gene qPCR and for generation of DNA templates for digoxigenin-labeled cRNA in situ probes. Total RNA from oil-exposed juvenile pink salmon ventricles was treated with DNase (Ambion #AM1906) and reverse transcribed with SuperScript III (Invitrogen #18080-400), following the manufacturer's instructions. qPCR was performed with 80 ng cDNA, 250 nmol l⁻¹ of each primer, and Fast SYBR Green (Applied Biosystems #4385612) on an Applied Biosystems ViiA 7 Real-Time PCR System for 40 cycles at 60°C annealing temperature. The DNA sequence of all amplicons was confirmed. The relative expression of target genes was normalized using the averaged values of two reference genes, wdtc1 and ef1α (Edmunds et al., 2014). For quantification of zebrafish nkx genes, hearts were dissected from 1 month old juveniles and >1 year old adults (AB strain) reared using standard procedures according to the animal care policy of the US Department of Commerce (Linbo, 2009), and was RNA extracted and processed for qPCR as above for salmon cardiac RNA.

In order to confirm the identity of the pink salmon cardiac transcriptome sequences for nkx2.3 and nkx2.5, Primer 3 was used to design gene-specific primers that were positioned as close to the 5′ and 3′ ends of the contig sequences as possible (see Table S2). Total RNA was reverse transcribed as described above, then used in a PCR reaction with the gene-specific primers. DNA sequencing of the resulting PCR amplicons was performed on a 3500 Genetic Analyzer (Applied Biosystems, Life Technologies).

For these molecular experiments, we utilized pink salmon embryos that were exposed from fertilization through late organogenesis to a graded series of environmentally relevant levels of Alaskan crude oil, yielding total PAH concentrations in water ranging from 10 to 45 µg l⁻¹, with tissue concentrations of 200–2500 ng g⁻¹ wet mass (Incardona et al., 2015). These doses are henceforth termed trace, low, medium and high treatment groups, corresponding with the original nominal oiled gravel loadings. This dose series was within the range that reduced juvenile-to-adult marine survival in pink salmon mark–recapture studies following the Exxon Valdez oil spill (Heintz, 2007; Heintz et al., 2000). After 8–10 months of post-exposure rearing in clean water, externally normal juveniles showed dose-dependent reductions in growth and maximum swimming speed, coupled to alterations in heart structure (Incardona et al., 2015). The ventricles of these embryonically exposed salmon had dose-dependent increases in length-to-width ratio (Fig. 1A) and reduced thickness of the compact myocardium (Fig. 1B) (representative examples from quantification previously detailed elsewhere: Incardona et al., 2015). In addition, there were signs of a hypertrophic response in the spongy myocardium, including increased numbers of cardiomyocyte nuclei, consistent with increased proliferation (Incardona et al., 2015).

We isolated RNA from ventricles and sequenced the transcriptomes of eight individuals from each treatment group, except the lowest exposure (trace) dose, for which we sequenced five. We selected individual ventricles to maximize statistical differences in shape across the exposure dose series, with length–width ratios ranging from 1.24 in controls to 1.43 for the high dose group (Incardona et al., 2015). These selections were representative of the peak frequency distributions for each treatment level (Fig. S1). First, we used 101 bp paired-end reads from a pooled sample of ventricular RNA to generate a de novo assembly of the juvenile pink salmon cardiac transcriptome. A total of 222,598 contigs were used in the de novo assembly. For single end reads from the 37 individual samples, 61.9±0.7% (mean 11.8 million) aligned uniquely to the de novo transcriptome and were used to generate read count values. Using several statistical approaches (detailed in Materials and Methods), a relatively small number of genes were identified as differentially expressed in a dose-dependent manner. Among all the pairwise comparisons between treatment groups, a total of 175 unique differentially expressed contigs were aligned to the UniProtKB/Swiss-Prot protein sequence database. Of these, 102 were successfully annotated (e-values <0.0001; see Table S1), with the remainder showing no significant sequence similarities.

Using IPA software, gene lists were analyzed for statistical enrichment of affected pathways. We used IPA for pairwise comparisons between genes from controls versus the genes from the high and medium doses combined, control versus high dose, control versus medium dose and control versus trace dose. As detailed in Materials and Methods, the low-dose group was excluded from this analysis. As a likely consequence of the relatively small number (<200) of DEGs, IPA did not detect alterations in the Top Canonical Pathways category. However, the Top Diseases and Bio Functions category identified DEGs associated with ‘cardiovascular system development and function’, ‘congenital heart anomaly’, ‘cardiac degeneration’, ‘cardiac stenosis’, ‘cardiac hypoplasia’ and ‘cardiac dysfunction’. In addition, DEGs associated with ‘cardiotoxicity’ were identified in the Top Tox Functions category. These included a contig originally annotated as nkx2.5 (upregulated 3.4-fold), latent transforming growth factor β (TGF-β)-binding protein 1 (ltpb1) (downregulated 1.3-fold), notch1 (downregulated 1.4-fold), stat1 (1.4-fold upregulated) and csf1r (1.5-fold upregulated). Moreover, manual curation of the gene lists identified 11 other genes associated with cardiac hypertrophy or cardiovascular development for a total of 14%. These included another downregulated gene that functions in TGF-β signaling, hemicentin1/fibulin6 (hmcn1/fbln6) (Table 1), and genes encoding four collagen isoforms that are targets of TGF-β signaling. In addition, one of the most highly upregulated genes (3.1-fold) was a contig annotated as another NKX family member, zampogna/nkx3.3.

In vertebrates, the nkx gene family encodes transcription factors that are widely involved in morphogenesis of all internal organs (Stanfel et al., 2005). Members of the family are highly similar in sequence, and it can be difficult to distinguish orthologous genes in different species. In fish, Nkx2.5 is the primary transcription factor involved in the determination of cardiomyocytes during early heart development, and is also required for proliferation of cardiomyocytes at later stages (Staudt and Stainier, 2012; Targoff et al., 2013). However, given the relatively low sequence similarity to multiple nkx2.5 orthologs (e.g. 42%), we examined the putative nkx2.5 contig sequence more carefully and used the sequence to amplify a cDNA from pink salmon cardiac RNA. Sequencing this 1716 bp cDNA revealed that the contig more likely represented an ortholog of the closely related family member nkx2.3. We confirmed this by fully annotating the pink salmon cardiac transcriptome (see Dataset 2, tab 1; detailed in Materials and Methods), identifying nkx2.5 contigs that were not differentially expressed, and sequencing additional nkx2.5 cDNA based on these. Comparison of these sequences unambiguously identified the differentially expressed gene as nkx2.3 (Fig. S2).

In mammals, Nkx2.3 functions in development of the intestine and lymphoid tissues (Pabst et al., 1999; Wang et al., 2000), and in lower vertebrate embryos it is also expressed in the pharyngeal arches (Buchberger et al., 1996; Evans et al., 1995; Lee et al., 1996). In chick embryos, transcripts for nkx2.3 were previously detected at later stages of heart development in differentiated cardiomyocytes (Buchberger et al., 1996), and there is evidence that nkx2.3 has some functional redundancy with nkx2.5 in amphibian heart development (Cleaver et al., 1996; Fu et al., 1998). It is expressed in the adult frog heart (Evans et al., 1995). However, nkx2.3 is not expressed during zebrafish embryonic heart development (Lee et al., 1996), and otherwise has not been functionally characterized in fish. Similarly, nkx3.3 is almost completely uncharacterized, likely because it appears to be specific to anamniotes and is absent from mammalian genomes. In Xenopus, nkx3.3 is expressed in the pharyngeal arches and the muscular layer of the developing gut (Newman and Krieg, 1999), but functional roles for the gene remain unknown. Additional sequence comparison differentiated the putative nkx3.3 contig from the closely related nkx3.2 gene (Fig. S3), which is present in mammals and birds.

We further characterized the expression of nkx2.3 and nkx3.3 in zebrafish, relative to nkx2.5. In developing embryos, as in Xenopus, nkx3.3 transcripts were detected in the pharyngeal arches, more strongly in arches 3–7 than in arches 1 and 2, but not the heart (Fig. 2A,B). While neither gene was expressed during early heart development (embryos and early larvae) in zebrafish (Lee et al., 1996), RT-qPCR showed that both nkx2.3 and nkx3.3 mRNAs were expressed in hearts isolated from 1 month old juvenile and mature adult zebrafish, at lower levels than nkx2.5 (Fig. 2C). The expression patterns of these three NKX family members are therefore consistent across fish species, with co-expression beginning during the growth of the juvenile heart and continuing into adulthood. As detailed below, we further characterized the patterns of cellular expression of these genes in ventricles from juvenile pink salmon hearts.

The function of nkx2.5 in early cardiac morphogenesis has been characterized extensively in zebrafish. Normal chamber formation and looping require precise levels of nkx2.5 expression, with abnormal morphology resulting from both loss of function and overexpression (Tu et al., 2009). Although not expressed in the embryonic heart, we determined whether nkx2.3 and nkx3.3 could influence zebrafish cardiac morphogenesis in a manner similar to nkx2.5 by injection of mRNAs for each gene into early cleavage stage embryos (1−2 cell stage). We tested nkx2.3 and nkx3.3 overexpression in two separate experiments, with nkx2.5 included in one experiment. In each case, nkx2.3 and nkx3.3 mRNA injection resulted in dose-dependent defects in cardiac morphogenesis, indicated functionally by the accumulation of pericardial edema (Fig. 3A–H; Fig. S4) and morphologically as defective looping of the chambers (Fig. 3I–L). In both experiments, nkx3.3 produced roughly twice the frequency of looping defects compared with nkx2.3 (Fig. 3M; ∼40% versus ∼20%, respectively at 100 pg mRNA). However, nkx3.3 was comparable in potency to nkx2.5 (Fig. S4; Fig. 3M), and the morphological phenotypes were indistinguishable. These findings are consistent with overlapping functionalities between nkx2.5, nkx2.3 and in particular nkx3.3 within fish cardiomyocytes.

Although not detected by pathway analysis, network analysis in IPA revealed that more than 20% (22/102) of the differentially expressed, annotated genes observed in the juvenile ventricular transcriptome function in innate immunity and inflammation (Table 2). In particular, the most significant network centered on two nodes derived from the stat1 and irf3 transcription factors (Fig. 4). In addition to stat1, this network also includes csf1r also identified in the ‘cardiotoxicity’ category above. Each of these DEGs was upregulated in response to embryonic crude oil exposures. Many are conventionally associated with antiviral responses, including interferon-responsive genes such as radical S-adenosyl methionine domain containing 2 (rsad2)/viperin, interferon-induced proteins with tetratricopeptide repeats (ifit1 and ifit5), interferon response factor (irf)3/7 and interferon induced protein 44 (ifi44). In addition, two pattern recognition receptors (PRRs) were upregulated, the nucleotide-binding oligomerization domain (NOD)-like receptor family member nlrc5 and the retinoid acid-inducible gene I (RIG-I)-like receptor DExH-box helicase 58 (dhx58). A potential PRR ligand was also upregulated, fibrinogen gamma chain (fgg). Although tlr4 expression was not altered in oil-exposed fish, tlr4 contigs were detected in the pink salmon cardiac transcriptome (see Dataset 2, tab 1). In addition, hsp70 (encoding heat shock protein 70), associated with the stress response and a potential PRR ligand (Chen and Nuñez, 2010), was increased in a dose-dependent manner ranging from 1.3- to 1.9-fold higher than in unexposed animals (Table 2; Table S1). The remainder (34) of the DEGs could not be placed directly into known roles for cardiovascular pathology or pathophysiology. We measured expression of five representative DEGs by qPCR, including nkx2.3, nkx3.3, ltbp1, rsad2 and ifit1. Expression levels determined by qPCR generally matched the read count data, with the strongest dose dependency for nkx2.3 and nkx3.3 (Fig. S5).

Neither oil-exposed nor control pink salmon showed any signs of infectious disease during juvenile growth, there were no differences in mortality, and cardiac tissues showed no obvious signs of inflammation (Incardona et al., 2015). To determine whether there may have been a subclinical viral infection, we fully annotated and searched the entire pink salmon cardiac transcriptome for virally encoded mRNAs or other viral response genes by filtering BLAST hit descriptions for mention of ‘virus’ or ‘viral’. Of 148,000 annotated contigs, only 536 contigs (0.4%) met these criteria (see Dataset 2, tabs 2 and 3).

The cardiac ventricle comprises several cell types in addition to cardiomyocytes. Cardiac fibroblasts are interspersed among cardiomyocytes and endothelial cells line the luminal spaces of the chamber. The salmonid heart also has a coronary vasculature, a source of additional cell types such as smooth muscle and nucleated and transcriptionally active fish erythrocytes. Moreover, even in the absence of infection or severe injury (e.g. infarct), resident lymphoid cells (e.g. macrophages) are interspersed throughout the myocardium. To further relate DEGs to the cardiac remodeling phenotype, we examined the expression patterns of representative genes by in situ hybridization in both the experimentally reared fish and juvenile pink salmon captured from the wild.

Both nkx2.3 and nkx3.3 were expressed strongly in cardiomyocytes (Fig. 5), identified by labeling with anti-myosin heavy chain antibody (Fig. 5A,B). nkx2.3⁺ nuclei were distributed sparsely but evenly across both the compact and spongy myocardium (Fig. 5C). Higher magnification images showed these to be cardiomyocyte nuclei in wild fish (Fig. 5C), and in both control and oil-exposed experimental fish (Fig. 5C′,C″). However, many of the nuclei had a more rounded appearance, possibly representing a different cell type such as fibroblasts or myofibroblasts. A slightly different pattern was observed for nkx3.3 (Fig. 5D). Although also expressed predominantly by cardiomyocytes, nkx3.3⁺ cells appeared to be more abundant in the spongy myocardium relative to compact myocardium, with a potential enrichment closer to the compact myocardium–spongy myocardium boundary. There were no qualitative differences between wild (Fig. 5D) and experimentally reared fish (Fig. 5D′,D″). Similarly, ltbp1 was also expressed primarily by cardiomyocytes, with strong expression in both the compact and spongy myocardium, with an apparent higher density of positive nuclei near the compact myocardium–spongy myocardium boundary (Fig. 5E). This pattern was also similar in both control (Fig. 5E′) and oil-exposed (Fig. 5E″) experimental fish. Finally, two of the interferon-responsive genes, ifit1 (Fig. 5F) and rsad2 (Fig. 5G), had similar patterns of expression. The signal for both genes was apparent in cardiomyocytes evenly distributed through both compact and spongy myocardium, but was also evident in endothelial cells and, in the case of rsad2, occasional erythrocytes. Both genes showed a similar expression pattern in wild (Fig. 5F,G) and experimental fish (Fig. 5F′,F″,G′,G″). All genes showed a predominantly nuclear/perinuclear expression pattern, which could reflect either their inherent subcellular localization (Piette et al., 1992; Russell and Dix, 1992) or possibly a re-distribution from a more peripheral localization during fixation due to cessation of contractility (Goldspink et al., 1997). The similar expression patterns in both wild and experimental fish indicate that cardiomyocyte gene expression is not an artifact associated with fish of hatchery origin.

In this study we used transcriptome analyses and in situ hybridization experiments to identify how crude oil exposure results in remodeling in the developing fish heart. The findings reported herein are consistent with ventricular cardiomyocytes as the primary source of crude oil-induced DEGs in pink salmon, several months after embryonic exposure. Oil exposure did not appear to alter the patterns of gene expression across different ventricular cell types. This indicates that altered expression detected by RNAseq and qPCR at the whole-organ level represents upregulation or downregulation at the individual cardiomyocyte level, rather than gain or loss of expression in distinct cell types, e.g. an inflammatory infiltrate. Cardiac remodeling in mammals is most pronounced after serious myocardial injury such as ischemia leading to infarct. Necrotic myocardial cells serve as initiating foci for activation and infiltration of lymphoid cells as well as the activation of fibroblasts, which culminates in a histologically visible inflammatory response and scar formation (Prabhu and Frangogiannis, 2016). Subsequent fibrosis, involving increased collagen deposition, is mediated in part by TGF-β signaling (Edgley et al., 2012). In contrast, cardiac remodeling in fish has been studied primarily in relation to adaptive physiology rather than pathological processes, and the normal fish heart shows considerable plasticity and a distinct capacity for cellular regeneration in comparison to mammals. Exercise-induced remodeling involves coordinated responses in both compact and spongy myocardium (Castro et al., 2013; Gamperl and Farrell, 2004), while remodeling in response to temperature acclimation or stress hormones is associated with differential effects in compact and spongy myocardium. For example, in salmonids, cold acclimation produces a thinner compact myocardium and an increased spongy myocardium with greater collagen deposition, likely mediated by TGF-β signaling (Johnston and Gillis, 2017), while warm acclimation leads to thicker compact myocardium, reduced spongy myocardium and reduced collagen content in both layers (Gamperl and Farrell, 2004; Keen et al., 2016; Klaiman et al., 2011). Normal cold adaptation thus involves collagen deposition (fibrosis), most likely to increase passive ventricular stiffness and force-generating capability at cold temperatures (Keen et al., 2016). Overall molecular responses also differ between cold and warm acclimation, with increased expression of innate immunity and stress-response genes in response to warming (Anttila et al., 2014).

Our results overlap with human pathological hypertrophy, but are also clearly distinct from normal physiological remodeling in fish, in terms of both histological changes and molecular response. First, the cardiac genes identified here such as notch1 and ltbp1 that have been previously linked to hypertrophic responses in mammals are likely a driving force behind the remodeling that arises in response to embryonic oil exposure. For example, notch1 is a key signaling receptor that plays well-known roles in both normal heart development and pathological hypertrophy in humans. Notch1 regulates the proliferative response of fish cardiomyocytes during the regenerative response to injury (Zhao et al., 2014), and reduction of Notch1 signaling in the mammalian heart leads to pathological hypertrophy (Croquelois et al., 2008; Nemir et al., 2014). In addition, development of specialized conductive cells that facilitate transmission of the cardiac pacemaker signal is dependent on Notch1 function (Rentschler et al., 2012). Therefore, the reduced levels of notch1 in oil-exposed fish are consistent with an abnormal hypertrophic response, as well as reduction in contactin 2 (cntn2), a marker for specialized conduction cells (Rentschler et al., 2012). While the laminin 2 (lama2) gene is normally associated with a specific type of muscular dystrophy not involving the heart, both LAMA2 and nidogen1 (NID1) were reduced in cardiomyopathic hearts from mice with Duchenne muscular dystrophy (Holland et al., 2013), and loss of LAMA2 in humans leads to congenital muscular dystrophy with associated cardiomyopathy (Nguyen et al., 2019). Cardiomyopathy is common in humans with Duchenne muscular dystrophy, and NID1 was also reduced in conventional human pathological hypertrophy (Kim et al., 2016). Thus, the findings here support a conserved pathology between fish and mammalian hearts. However, decreased expression of collagen genes (Gil-Cayuela et al., 2016; Goldsmith et al., 2014) and genes involved in augmenting TGF-β signaling (ltbp1 and hmcn1/fbln6) are indicative of an anti-fibrotic rather than pro-fibrotic response (Chowdhury et al., 2017, 2014; Johnston and Gillis, 2017; Spinale et al., 2010; Zavadzkas et al., 2011), suggestive of key distinctions between the regulation of fish pathogenesis and regeneration. Finally, roles for arrestin domain-containing gene 3 (arrdc3) and rnf213 in hypertrophy are not clear, but both are involved in angiogenesis (Nauta et al., 2017; Scholz et al., 2016; Wen et al., 2016), and Arrdc3 protein regulates β-adrenergic receptors (Patwari et al., 2011; Tian et al., 2016).

Under the conditions typically associated with remodeling in the injured human heart, elevated expression of immune and inflammatory genes would be expected to involve the cells of the immune system or cardiac fibroblasts (Diaz-Araya et al., 2015). In contrast, we observed upregulation of these genes intrinsically in cardiomyocytes in the absence of histologically visible cellular injury or death. Moreover, we observed only a very small fraction of transcripts (0.4%) associated with common fish viruses, and these were not differentially regulated. Most of these were related to retroviral proteins and oncogenes, and none were associated with viruses known to cause cardiac disease in fish, such as piscine retrovirus, piscine myocarditis virus or salmonid alphavirus (Haugland et al., 2011; Herath et al., 2016; Palacios et al., 2010). In addition, the transcriptional response of the salmonid heart to active viral infection includes a much larger suite of genes involved in both innate and adaptive immunity, at much higher levels of induction (Johansen et al., 2015; Timmerhaus et al., 2011). These findings indicate that a viral infection is highly unlikely to have caused the differential expression of the identified innate immunity and inflammatory genes. The upregulation of PRRs and PRR ligands such as the clotting factor fibrinogen (fgg) (Hodgkinson et al., 2008; Millien et al., 2013), which is normally expressed only in the liver (Fish et al., 2012), and hsp70 suggest the action of a non-infectious stressor. PRRs were first identified as receptors for non-host constituents of pathogens, e.g. viral DNAs and protein components of bacterial cell walls. However, as originally proposed over 20 years ago (Matzinger, 1994), PRRs recognize ‘hidden-self’ molecules that are released or exposed during cellular stress (Schroder and Tschopp, 2010). Consequently, PRRs and the innate immune response are increasingly associated with non-infectious pathophysiological processes, as well as normal tissue homeostasis (Anders and Schaefer, 2014; Zhang et al., 2015). In particular, recent studies have linked innate immune responses such as interferon signaling to pathological hypertrophic responses in the mammalian heart.

Knockout studies in the mouse indicated that irf genes play a direct role in cardiac hypertrophy (Sun and Wang, 2014), with IRF3/7 playing a protective role (Jiang et al., 2014; Lu et al., 2013). Thus upregulation of IRF3/7 in pink salmon hearts could represent a compensatory response to stem pathological hypertrophy initiated by oil exposure. Similarly, although the physiological functions of interferon-inducible card6 gene are not well characterized, genetic studies in the mouse indicate that this gene also plays a protective role against pathological hypertrophy (Li et al., 2014), and NLRC5 was shown to reduce the fibrotic response to TGF-β in cardiac fibroblasts (Zhou et al., 2017). While RSAD2/viperin has been studied almost exclusively in the context of antiviral properties, there are three previous mechanistic associations with processes potentially acting in cardiac remodeling. The gene was initially identified through its increased expression during bone formation and stimulation by signals causing osteoblast proliferation (Gerwal et al., 2000), and its induction during an ischemic response in the pancreas (Drognitz et al., 2006). More importantly, RSAD2/viperin was upregulated in cardiomyocytes in response to prorenin, a blood pressure-regulating protein secreted by the kidney that causes pathological cardiac hypertrophy when overexpressed (Saris et al., 2006). It is therefore likely that RSAD2/viperin in oil-exposed pink salmon was involved in the proliferative response in spongy myocardium, possibly secondary to a reduced availability of oxygen (ischemia). Finally, many of the same or closely related genes were recently identified as differentially expressed markers of human hypertrophic cardiomyopathy (Chen et al., 2019). Our findings therefore suggest that reduced cardiac function stemming from exposure of embryonic fish to crude oil results in abnormal development and increased stress in heart cells, activating the PRR ligand-mediated signaling and the consequential pathological inflammatory response mediated by the PRRs and interferons. Functional studies will be necessary to determine which of these immune/inflammatory genes reflect a compensatory response to hypertrophy, as opposed to those that might contribute directly to abnormal proliferation and remodeling post-exposure.

While nkx2.3 and nkx3.3 are low-abundance genes in juvenile and adult fish hearts, relative to nkx2.5, both genes were highly upregulated (>3-fold) in salmon hearts with relatively subtle phenotypic changes far downstream of embryonic oil exposure. These oil-exposed hearts were overtly normal and differed only by mild changes in ventricular dimensions and thickness of the compact myocardium (Incardona et al., 2015), which is a very different situation compared with acute cardiac injury and regeneration models. In this case, upregulation of these nkx family members is associated with a process that represents a chronic response to embryonic injury (pathophysiological remodeling), rather than an acute response to tissue necrosis or ablation. While neither of these genes was observed to be upregulated in zebrafish cardiac ablation/regeneration models (Kang et al., 2016; Karra et al., 2018; Lai et al., 2017), upregulation may not be required if nkx3.3 serves as a permissive factor. At the same time, however, nkx3.3 was upregulated in regenerating muscle following tail amputation in the newt Pleurodeles waltl, but this was observed 20–30 days post-injury (Nicolas et al., 1999), a time point beyond the aforementioned zebrafish regeneration studies. The changes in nkx3.3 expression observed here occurred well over 8 months after cardiotoxic injury to the embryo (Incardona et al., 2015).

In summary, in the hearts of juvenile pink salmon embryonically exposed to crude oil, we observed histological changes that were similar to the normal remodeling that occurs in fish in response to cold acclimation. These included a thinner compact myocardium and hypertrophy of the spongy myocardium. However, the molecular response was more aligned with warm acclimation, as evidenced by a large increase in the expression of innate immune and inflammatory genes and reduced fibrotic pathways. Importantly, the expression of innate immune genes was attributable to cardiomyocytes and occurred in the absence of obvious inflammatory infiltrates. How the identified DEGs relate to altered ventricular shape is unclear at present, but identification of downstream targets of nkx2.3 and nkx3.3 could provide insight. Alternatively, altered ventricular dimensions could arise passively as a result of altered mechanical forces related to changes in compliance. The reduced performance associated with embryonic oil exposure (Incardona et al., 2015) thus likely arises from both a reduction in compact myocardium and a lessened ability to transmit mechanical force due to reduced collagen expression and increased compliance.

These findings have important implications for the environmental health of both wild fish species and humans. All of the DEGs identified here have potential as novel molecular biomarkers for adverse impacts on fish populations from oil spills, as well as chronic low-level pollution from other PAH-containing sources such as urban stormwater runoff. Considering both the similarities and differences between the mammalian heart and the hearts of aerobically active fish, such as salmon, these findings may also provide insight into pathological remodeling in mammalian hearts and regenerative potential. Whereas the mammalian heart grows largely by hypertrophy of a limited number of fetal cardiomyocytes, the fish heart continues to increase in size in correspondence with overall growth by continuously adding new cardiomyocytes through proliferation. The unique capacity for the fish heart to regenerate cardiomyocytes after severe injury is linked to this proliferative ability (Foglia and Poss, 2016). The association of nkx2.3 and nkx3.3 with hyperproliferation during fish heart remodeling suggests important roles for these genes in post-embryonic cardiac growth. The role of nkx3.3 in fish heart growth and amphibian tail regeneration and its absence from mammalian genomes are perhaps not coincidental, and its function warrants further investigation. These studies demonstrate how an expanding mechanistic understanding of oil spill and other environmental pollution impacts, including the chemical perturbation of biological systems, can ultimately yield new insights into fundamental biological processes critical to understanding cardiac structure and function.

We thank Maryjean Willis for assistance with preparing sections for in situ hybridizations, and Louisa Harding and anonymous referees for critical review of the manuscript.