RNA Sequencing of Lake Charr Epidermal Mucus to Assess Molecular Effects of Diluted Bitumen Exposure in a Boreal Lake

Transport of diluted bitumen (dilbit) from Canada’s oil sands region poses risk for leaks and spills of petroleum-derived contaminants into the environment. Exposure of fish to dilbit is known to cause cardiotoxicity, developmental deformities, and impairment in swim performance. However, previous studies have examined the toxicity of dilbit in laboratory settings which does not account for environmental and biological food-web variables that may alter exposure and/or toxicity of dilbit. Moreover, most methods of assessing organism health following oil exposure require lethal sampling. This work is a part of a larger set of experiments where dilbit spills were simulated within enclosures on a lake; the present study assesses the impacts of residual levels of dilbit that may have entered the surrounding lake environment from the enclosures following model spill cleanup. In order to understand the impacts of residual dilbit in an ecosystem setting without use of lethal sampling, epidermal mucus was collected and sequenced from lake charr (Salvelinus namaycush) exposed to residual dilbit in a boreal lake. While concentrations reached a maximum of 2.29 μg/L total polycyclic aromatic compounds (ΣPAC) within surface waters, surface water ΣPAC concentrations generally remained below 1 μg/L. Results of RNA sequencing were compared to sequencing data from mucus collected prior to dilbit additions. Differential gene expression and pathway analyses indicated dysregulation of genes associated with intermediary and energy metabolism as well as a trend in upregulation of cyp1a3 in epidermal mucus following dilbit exposure. Thus, results of the present study suggest that lake charr undergo consistent biological responses after exposure to residual levels of dilbit following a model spill, and that mRNA-based analysis of mucus may be a viable method for non-lethal oil exposure assessment. Overall, the results provide insight on the response of wild fish to very dilute dilbit exposures after a model spill cleanup.


INTRODUCTION
The oil sands region of Alberta, Canada has one of the largest known reserves of bitumen, containing an estimated 50 billion cubic meters of the heavy type crude oil (Canada National Energy Board, 2006). Due to its viscous nature, bitumen is mixed with varying amounts of light oil diluent (e.g., naptha-based condensate) to decrease overall viscosity for more efficient transport, thus forming diluted bitumen (dilbit). Dilbit is primarily exported from the oil sands region via extensive transcontinental pipeline networks, raising concern for potential leaks and spills into lakes and rivers located along transportation routes (Dupuis and Ucan-Marin, 2015;Canada Energy Regulator, 2018). Although export by rail has increased in recent years, volumes remain small compared to pipelines but present an additional potential spill vector (Canada Energy Regulator, 2018). Multiple significant dilbit spills have occurred in North America, with the largest being the 2010 Enbridge Line 6B pipeline rupture which released 1.15 million gallons of dilbit into the Kalamazoo River in Michigan (Crosby et al., 2013;National Academies of Sciences Medicine and Engineering, 2016). The spill impacted 1,560 acres of instream habitats as well as floodplain and upland areas and resulted in ecological injury to fish, benthic invertebrates, birds, mammals, and other wildlife (USFWS, 2015). Although major dilbit spills have occurred in the past and future spills are inevitable, there remain significant research gaps on the behavior of dilbit in the environment and its toxicity to aquatic organisms. Furthermore, oil sands dilbit production and export is projected to continue to expand due to significant economic incentives (Crosby et al., 2013;Boufadel et al., 2015). Continual increases in transport of dilbit and proposed expansions of existing pipeline networks make it critical to fully characterize the effects of dilbit on aquatic ecosystems for adequate risk assessment of future spills and proper guidance in spill response decision making.
As bitumen is heavily biodegraded over geological time, dilbit contains a high proportion of heavy molecular weight compounds such as resins, asphaltenes, and high molecular weight polycyclic aromatic compounds (PACs) and alkylated PACs (Yang et al., 2011;National Academies of Sciences Medicine and Engineering, 2016). Moreover, the chemical fingerprint of dilbit (as quantified by n-alkanes, PACs and their alkylated homologues, terpanes, steranes, bicyclic sesquiterpanes, etc.) is highly variable and dependent on extraction and refining processes, geographic source, and the composition and volume of diluent added by the producer (Adams et al., 2013;Lee et al., 2016). Following a spill, the chemical properties of dilbit are further modified by complex weathering processes such as photo-oxidation, biodegradation, evaporation, emulsification, dispersion, and sedimentation (Dew et al., 2015;National Academies of Sciences Medicine and Engineering, 2016). Due to its complex chemical composition and multiple factors affecting its chemical properties, it may be difficult to understand the fate and behavior of spilled dilbit in freshwater environments. In order to better understand dilbit's behavior in low-energy lake systems, Stoyanovich et al., 2019 conducted experimental releases of dilbit in outdoor microcosms containing natural lake water and sediments (Stoyanovich et al., 2019). As volatile hydrocarbons evaporated, dilbit increased in density and viscosity, resulting in its submergence after 8 days. Furthermore, there was a rapid accumulation of PACs in the water column, which tend to be drivers of oil toxicity in aquatic organisms.
The toxicity of dilbit in freshwater environments has been minimally studied and has only recently begun to receive more extensive research (Dupuis and Ucan-Marin, 2015;National Academies of Sciences Medicine and Engineering, 2016). Similar to toxic endpoints associated with conventional crude oil, dilbit exposure in fish has shown to induce mortality, decreased utilization of energy stores, metabolic alterations, reduced growth, oxidative stress, behavioral changes, developmental malformations, delayed hatching, changes in heart rate, abnormal inflation or complete un-inflation of the swim bladder, among other adverse effects (Madison et al., 2015(Madison et al., , 2017Philibert et al., 2016;Alsaadi F. et al., 2018, Alsaadi et al., 2018Alderman et al., 2018;Avey et al., 2020). Impairments induced by dilbit exposure may have impacts on organismal fitness and, thus, population health. For example, exposure of sockeye salmon parr (Oncorhynchus nerka) to water soluble fractions of dilbit (3.5-66.7 μg/L ΣPAC) resulted in concentration-dependent cardiac remodeling that correlated to reductions in swimming performance (Alderman et al., 2017b). Dilbit may also impede exercise recovery, as dilbit-exposed sockeye salmon experienced increased cellular damage following exercise as suggested by an increase in tissue leakage proteins, such as creatine kinase, in serum compared to controls (Alderman et al., 2017a). As reduced swim performance may impair migratory success, dilbit spills may pose risk to sensitive salmonid populations in freshwater ecosystems.
Following a freshwater oil spill, there is a need for rapid exposure assessment of biota in the impacted area. Currently, many methods used to assess ecological damage to organisms, such as those used following the Kalamazoo River spill, require lethal sampling which may result in further impairment (USFWS, 2015). Use of exposure assessment tools which are nonlethal and minimally invasive is important to reduce further perturbations on animal populations. In fish, molecular analysis of epidermal mucus may be a viable tool for monitoring effects of oil exposure as it is nonlethal, rapid, inexpensive, and sample procurement requires little training. Following exposure to Deepwater Horizon slick oil, RNA sequencing of mucus from juvenile mahi-mahi (Coryphaena hippurus) identified molecular changes in mucus associated with oil toxicity, such as upregulation of cytochrome P450 1A (CYP1A) and dysregulation of transcripts associated with cardiotoxicity and immunosuppression (Greer et al., 2019). Moreover, exposure of dusky splitfin (Goodea gracilis) to crude oil resulted in increased antioxidant defense response in the epidermal mucus (Dzul-Caamal et al., 2016). Results from these previous studies suggest that epidermal mucus may be a useful target for environmental monitoring following oil spills. However, previous studies were laboratory-based and may not reflect changes in chemical properties of oil brought about by environmental conditions such as a evaporation, photooxidation, emulsification, dispersion, biodegradation, and sedimentation (National Academies of Sciences Medicine and Engineering, 2016).
Although research suggests that mucus may be a promising method for nonlethal sampling, the efficacy of mucus as an oil exposure assessment tool has not yet been assessed on wild populations in the field. In the present study, epidermal mucus was collected from lake charr (Salvelinus namaycush) residing in a pristine Canadian boreal shield lake at the International Institute for Sustainable Development Experimental Lakes Area (IISD-ELA) where experimental dilbit spills and cleanups were conducted within enclosures on the lake (detailed in Rodriguez-Gil et al., 2021). While experimental spills were conducted within enclosures, minute levels of dilbit entered the surrounding lake environment leading to potential exposure of resident lake charr. To minimize impacts on the small and vulnerable lake charr population, lethal or invasive sampling of tissue or blood was not used and, instead, epidermal mucus was collected. RNA sequencing was conducted on lake charr mucus collected prior to and after dilbit additions to identify altered molecular pathways associated with dilbit exposure. Lake charr may serve as a model species for the study as many dilbit pipelines are concentrated in areas of Canada where they intersect freshwater ecosystems containing sensitive salmonid species. Results of the present study will identify molecular pathways altered by dilbit in a whole-lake setting and further contribute to the development of mucus as a nonlethal tool for oil exposure assessment. Overall, the present study provides insight on the response of wild fish to very dilute µg/L dilbit exposure after a model spill cleanup.

Study Area and Diluted Bitumen Additions
The study was conducted at IISD-ELA, a freshwater research center located near Kenora, Ontario, Canada (49°40′N, 93°44′W) (https://www.iisd.org/ela/). Being located in a sparsely populated region, the lakes at IISD-ELA are mostly unaffected by anthropogenic influences. IISD-ELA's Lake 260 was selected for the present study (49°42′N, 93°46′W) which has an area of 331,328 m 2 , a volume of 1,975,971 m 3 , an average depth of 5.9 m, and a maximum depth of 15.7 m. Lake 260 is a boreal lake and serves as a model for freshwater water bodies in Canada's Boreal Shield ecozone.
The present study was a part of a larger series of dilbit experiments at IISD-ELA. In order to examine the fate, behavior and potential impacts of spilled dilbit in freshwater lakes, dilbit spills were simulated within 10 m diameter,~100 m 3 volume limnocorrals on Lake 260 beginning in 2018 (Orihel, 2019). The experimental design of the simulated dilbit spills are described in a previous publication (Rodriguez-Gil et al., 2021). Briefly, limnocorrals consist of a floatation collar that supports an impermeable curtain that extends to the bottom of the lake where it is sealed to the sediments using sandbags, providing an enclosed ecosystem to model whole lake processes. A series of 9 limnocorrals were constructed and deployed in water depths of~1.5 m. After installation and pre-spill monitoring, Cold Lake winter Blend dilbit was added to seven limnocorrals using a regression with added volumes ranging from 1: 100,000-1:100 dilbit to water (v:v). The total volume of dilbit added to the enclosures was~382 L. Chemical and physical weathering, as well as potential toxicity, of the dilbit were studied for 70 days, after which residual surface dilbit was removed from the limnocorrals. The limnocorral curtains were then cut at mid water column depth prior to winter freeze up to allow continued study of the dilbit that remained on the lake bottom. As a result, there was potential for dilbit to enter the lake water column throughout the duration of the winter months and prior to final sediment cleaning, which took place in August 2019. Additionally, in June 2019 fourteen shoreline enclosures (5 × 15 m and enclosing volumes of 25,000-30,000 L of lake water) were treated with 1.5 L of weathered dilbit (Palace, 2021). While the environments were enclosed, tracer studies indicated that some leakage to the outside lake environment occurred, potentially contributing to enriched PAC concentrations in the lake. The present study thus examined the potential effects that residual dilbit entering the surrounding lake environment following cleanup had on resident fish.

Water Extractions and Analytical Chemistry
Surface water samples (1 L) were collected from the lake at 10 time points during the study: 3 days before oil addition and on days 1,2,3,5,9,13, 18, 53, and 80 after application. Samples were filtered first using Whatman GF/C filters and then transferred to 2 L separatory funnels. Each funnel was spiked with 20 µl of a 5 ppm recovery internal standard which consisted of d8-naphthalene, d8-acenaphthylene, d10-acenaphthene, d10-fluorene, d10-phenanthrene, d10pyrene, d12-Benz(a)anthracene, d12-chrysene, d12-benzo (b)fluoranthene, d12-benzo (k)fluoranthene, d12-benzo (a) pyrene, d12-indeno (1,2,3-c,d)pyrene, d14-dibenzo (a,h) anthracene, and d14-benzo (g,h,i)perylene. Ten grams of NaCl was dissolved into each sample prior to the liquid: liquid extraction. The samples were extracted by adding 50 ml of dichloromethane (DCM) to each funnel and gently swirling for 1 min. The DCM was drained into a round bottom flask and the water extraction was repeated. The combined 100 ml DCM extract was then evaporated to a final volume of 1 ml and spike with 20 µl of 5 ng/μl d10-anthracene as an instrument performance internal standard. Samples were vortexed and transferred to 2 ml amber glass autosampler vials for storage at 4°C until analysis by GC-MS/MS. An Agilent 7,890 GC coupled to a 7,000°C triple quadrupole mass spectrophotometer fitted with an electron ionization (EI) source used for the native and d-labelled PAC and alkyl PAC quantification. Helium was used as the carrier gas at a flow rate of 1.2 ml/min with an Agilent J&W HP-5ms ultra inert column (30 m × 0.25 mm, 0.25 μm film thickness). Sample volumes of 1 μl were injected by a PAL RSI 85 auto sampler at a temperature of 60°C. GC conditions and details of the multiple reaction monitoring precursor and product ion transitions for the 16 priority PACs and alkylated PACs that were analyzed in the samples are included in Idowu et al., 2018. Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 836640

Samples Collection
The comparisons of interest in the study were 1) 2017 and 2018 to determine yearly gene expression changes in baseline samples and 2) 2018 and dilbit-exposed samples (2019) to determine gene expression changes following dilbit exposure. Thus, prior to dilbit additions into enclosures in September 2017, seven lake charr were collected from the surrounding lake environment and sampled for epidermal mucus. In September 2018, seven additional mucus samples were collected prior to cutting of limnocorral curtains and exposure of the lake water column to sunken dilbit. Finally, mucus was collected in September 2019 following exposure of dilbit to the lake water column. Thus, the experimental design is a temporal comparison pre-and postexposure. Fish were captured from the lake from a boat using hook and line and were placed in separate bins of fresh lake water until processing (held <5 min). Retrieval and post-capture air exposure were minimized to reduce physiological stress. Water was renewed for each fish to prevent cross-contamination. Lake charr were placed in a solution of 0.1 g/L pH buffered (7.0) tricaine methanesulfonate (MS-222) mixed in lake water until fin movement ceased and fish were unresponsive to light pressure on the caudal fin. Total length, fork length, and weight of the individuals were recorded. Sex was determined through external morphological characteristics if possible. Mucus was collected by gently scraping the skin with a stainless-steel spatula, avoiding displacement of scales, to aggregate the mucus along the body axis in the anterior to posterior direction. Mucus was aspirated using a syringe, dispensed into a microcentrifuge tube, placed on dry ice, and stored at −80°C until analyzed. Lake charr were then released back into the lake. As the lake charr population is small and, thus, vulnerable at IISD-ELA, additional tissue types or blood via lethal/invasive sampling were not collected. Tools and work area were cleansed with ethanol and then lake water between fish to prevent crosscontamination. Condition factor (K) was calculated for each individual using formula K = weight/length3 × 100. Statistical differences in average weight, total length, and condition factor among 2017, 2018, and 2019 lake charr were tested using a oneway ANOVA and Tukey multiple comparisons test at a significance threshold of p ≤ 0.05. All capture and handling procedures were conducted in accordance with the Canadian Council on Animal Care approved protocols.
cDNA Library Generation and Illumina Sequencing
The comparisons of interest in the study were 1) 2017 and 2018 to determine yearly gene expression changes in baseline samples and 2) 2018 and dilbit-exposed samples (2019) to determine gene expression changes following potential dilbit exposure. RSEM (Li and Dewey, 2011) was used to quantify abundance of transcripts in each sample. To visualize transcriptome-wide patterns among 2018 baseline and 2019 dilbit-exposed samples, unsupervised hierarchical clustering was used to generate a sample-to-sample Euclidian distance heat map using the RSEM count matrix that underwent a regularized-logarithm transformation. Principal component analysis (PCA) was also performed on the transformed RSEM count matrix for baseline and dilbit-treated samples. Unsupervised hierarchical clustering was used to produce a heat map of the top 500 most-variable transcripts among 2018 and dilbit-treated samples. The top 500 most-variable transcripts were extracted from the transformed count matrix using the genefilter R package (Gentleman et al., 2019). The nontransformed transcript count matrix was input into DESeq2 v.1.23.6 (Love et al., 2014) to determine differentially expressed transcripts among 2018 and dilbit-treated samples at p ≤ 0.05 after Benjamini-Hochberg false discovery rate (FDR) correction. Differential expression analysis was also conducted among 2017 and 2018 baseline data to examine yearly changes in  (Livak and Schmittgen, 2001) with β-actin as the normalizing gene. Primer sequences were designed using IDT PrimerQuest Tool ® with transcript sequences from the brook trout assembly. Supplementary Table S1.

Soluble Dilbit Constituents in Lake 260
The present study was a part of a larger series of experiments at IISD-ELA where dilbit spills and cleanups were carried out within limnocorrals and shoreline enclosures in a boreal lake (Rodriguez-Gil et al., 2021). There were two potential vectors for dilbit constituents entering the surround lake environment. June 2019, tracer studies indicated that some leakage from the enclosures to the outside lake environment occurred which may have contributed water soluble dilbit compounds into the lake. Thus, the present study assessed the potential impacts of residual dilbit following a spill cleanup on wild fish residing in the lake outside of the enclosures. ΣPAC concentrations of 0.058 μg/L were detected in lake surface waters prior to dilbit additions into shoreline enclosures in June 2019 ( Figure 1); these levels could be linked to residual dilbit in sediments following previous 2018 dilbit experiments, through atmospheric sources, or from background PACs from organic sediment inputs. Following the additions of dilbit into shoreline enclosures on June 22 (Day 0), there was an increase in PACs within surface waters, reaching a maximum concentration of 2.29 μg/ L ΣPAC on Day 5 (Figure 1). ΣPAC concentrations in surface waters declined to 1.05 μg/L on Day 6 and remained under 1 μg/L for the remainder of the exposure period ( Figure 1). By Day 66, ΣPAC concentrations in surface waters declined to 0.082 μg/L ΣPAC concentrations at the time of sampling mucus from dilbit-exposed lake charr was approximately 0.067 μg/L (interpolation from Day 66 and Day 81 data points) (Figure 1). Three-and 4-ring compounds were the most dominant of the total PACs, with phenanthrenes, dibenzothiophenes, pyrenes, retene, and napthalenes being the most abundant compound groups detected in lake water outside of the enclosures (Supplementary Material; additional manuscript forthcoming further detailing these chemistry data). As dilbit was exposed to weathering in the field, it is likely that lighter molecular weight compounds volatilized, leaving an increased proportion of higher molecular weight 3-and 4-ring PACs. Chemistry results are supported by previous studies simulating dilbit spills in outdoor tanks containing lake water and sediments where rapid desorption of PACs into the water column was observed, most notably phenanthrenes, dibenzothiophenes, and fluorenes (Stoyanovich et al., 2019).
Routes of exposure to oil constituents in fish include ventilation through gills, dermal uptake, and ingestion (Nichols et al., 1996;Almeda et al., 2013). As dilbit constituents were detected within the water column, absorption of waterborne compounds through skin and gills was one probable route of exposure. Moreover, lake charr may also have ingested oil-contaminated organisms, such as zooplankton and forage fish which were more likely to have contact with sunken dilbit (Almeda et al., 2013; National Academies of Sciences Medicine and Engineering, 2016). Through consumption of oil-contaminated material, dilbit constituents such as PACs may be transferred to apex predatory species through trophic transfer within the food web. Thus, exposure to PACs may continue even after most constituents have volatilized from the water column. Finally, there was also potential for exposure via direct contact with sunken dilbit that remained on the sediments. Additional studies on the trophic movement of dilbit constituents in lakes are needed to better understand bioaccumulation processes in boreal lake systems.  Table  S2). Average weights for 2017, 2018, and 2019 lake charr were 644.3 ± 87.9 g, 814.3 ± 92.8 g, and 961.9 ± 71.5 g, respectively ( Figure S1). Average total lengths were 42.59 ± 1.9 cm in 2017, 46.8 ± 1.7 cm in 2018, and 48.74 ± 1.3 cm in 2019 (Supplementary Figure S1).  Table S3).
To determine year-year differences in background transcriptional signatures, comparisons were made in fish sampled prior to dilbit exposure in 2017 and 2018. A total of 1,576 transcripts were differentially expressed among 2017 and 2018 baseline years at a FDR adjusted p-value ≤ 0.05 (Supplementary Material). The top GO biological processes were "translation", "ribosomal small unit assembly", and "in utero embryonic development", whereas the top GO molecular functions were "structural constituent of ribosome", "poly(A) RNA binding", and "mRNA binding" (Supplementary Material).
Comparing animals pre-and post-dilbit exposure, the sampleto-sample Euclidian distance heat map (Supplementary Figure  S2) and PCA of transcriptome-wide count data ( Figure 2A) indicated clustering of individuals based on treatment. Unsupervised hierarchical clustering of the top 500 mostvariable transcripts also clustered individuals by treatment ( Figure 2B). Thus, results of clustering analyses indicate that there were clear differences in the mucus transcriptomic profiles of lake charr following additions of dilbit. PCA of transcriptomewide count data did not show distinct clustering patterns based on sex, suggesting that transcriptomic changes among years were driven by treatment (Supplementary Figure S3).
A total of 714 transcripts were differentially expressed at a FDR adjusted p-value ≤ 0.05 among individuals sampled in preand post-dilbit (Supplementary Table S4). The RNA sequencing results were confirmed through qPCR (Supplementary Figure  S4). Among all pathway analyses, there was prevalence of biological pathways related to intermediary metabolism, more specifically relating to carbohydrates and amino acids. For example, of the 12 KEGG pathways identified, three were related to intermediary metabolism including 'arginine biosynthesis' (rank 3), 'alanine, aspartate, and glutamate metabolism' (rank 5), and 'carbon metabolism' (rank 11) (Figure 3; Supplementary Table S5). Intermediary metabolism pathways were also identified in the GO analysis, including the biological process 'urea cycle' (rank 11) as well as the molecular functions 'glucose binding' (rank 9) and 'carbamoyl-phosphate synthase activity' (rank 18) (Figure 3; Supplementary Table S6; Supplementary Table S7). This trend was also seen in the IPA results, where the top canonical pathways (ranked by p value) included 'superpathway of citrulline metabolism' (rank 1), 'arginine degradation 1' (rank 10), 'urea cycle' (rank 18), and 'arginine  (Figure 4). Taken together, results of the pathway analyses suggest molecular changes in metabolic and energetic processes in lake charr following dilbit exposure. It should also be pointed out that the number of differentially expressed genes was greater among baseline year comparisons, suggesting that differences  in environmental variables among years contribute a significant amount of variation. While field-based toxicological studies may be beneficial as they capture effects from multiple stressors, additional laboratory studies may aid in discerning effects from potential confounding environmental variables (discussed further in text). Several studies in fish have shown that exposure to PACs and oil, can alter energy metabolism mechanics and induce energy mobilization causing increased energy demands due to stress (Wendelaar Bonga, 1977;Frasco and Guilhermino, 2002;Tintos et al., 2008;Vijayan et al., 2010;Ings et al., 2012;Hook et al., 2018). Stress-induced changes in metabolism are associated with a surge of cortisol, which results in activation of gluconeogenesis and, consequently, increased plasma glucose concentrations. Moreover, gluconeogenesis relies on non-carbohydrate sources such as lactate, glycerol, and amino acids for the formation of glucose. Due to the multiple pathways associated with amino acid metabolism, it is suggested that lake charr catabolized and degraded protein for energy production (Figure 3; Figure 4). Metabolism of amino acids for energy is further supported by the upregulation of alanine aminotransferase 2 (gpt2) which converts alanine into pyruvate for glucose production and has been shown to be enhanced in response to stress for gluconeogenesis ( Figure 5) (Mommsen et al., 1999;Kumar et al., 2010;Qian et al., 2015). Other gluconeogenic and amino acid catabolism transcripts were also upregulated, including glutamate dehydrogenase (glud1), glutamine-fructose-6-phosphate aminotransferase (gfpt1), and ornithine aminotransferase (oat) ( Figure 5). Furthermore, several of these transcripts are also involved in the urea cycle, likely enhancing the urea cycle to convert excess ammonia following breakdown of amino acids ( Figure 3; Figure 4; Figure 5). Cortisol also elicits decreased potential for glycolysis in fish through decreased activities of pyruvate kinase, hexokinase, fructose 1,6bisphosphate, and glucose-6-phosphate 1-dehydrogenase (Laiz-Carrión et al., 2003). In the present study, glycolytic transcripts such hexokinase-1 (hk-1) and glucose-6-phosphate 1-dehydrogenase (g6pd) were downregulated in lake charr mucus following dilbit exposure ( Figure 5). As hexokinase is a rate-limiting enzyme of glycolysis, decreased transcription of hk-1 and other glycolytic transcripts supports the utilization of gluconeogenesis during stress-induced energy production. Although molecular changes suggested a shift in energetics, there were no differences in condition factor pre-and postdilbit exposure, suggesting that the disruption was not enough to elicit changes in body mass.
Shifts in energy metabolism have been observed in fish following dilbit exposure. Previous studies exposing Atlantic salmon (Salmo salar) smolts and sockeye salmon embryos to dilbit resulted in elevated cardiac and total body triglyceride levels, suggesting decreased lipid utilization for energy compared to controls (Alderman et al., 2018;Avey et al., 2020). In the present study, decreased lipid utilization was suggested by downregulation of the transcript encoding for carnitine palmitoyltransferase I (cpt1a) which is the rate limiting step in long-chain fatty acid oxidation ( Figure 5) (Frøyland et al., 1998;Qu et al., 2016). Furthermore, expression of genes involved in protein degradation was increased in muscle of dilbit-exposed Atlantic salmon smolts, namely ubiquitin which labels protein for proteasomal degradation (Avey et al., 2020). Although ubiquitin was not upregulated in lake charr mucus, several transcripts encoding for enzymes involved in ubiquitination FIGURE 5 | Fold change of select transcripts dysregulated in dilbit-exposed lake charr relating to energy metabolism and stress response from RNA sequencing results. All transcripts were significant at FDR adjusted p-value ≤ 0.05, except cyp1a3 which was significant at a FDR adjusted p-value ≤ 0.1.
Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 836640 of proteins were upregulated, including polyubiquitin, ubiquitin conjugating enzyme E2 D2, and others (Supplementary Table S4). Thus, results from the present and previous studies together suggest that dilbit may induce shifts in energy metabolism, such as increased protein degradation and potentially greater utilization of amino acids for energy. Although dilbit exposure has been shown to increase muscle damage in sockeye salmon (Alderman et al., 2017b), it is not known whether muscle damage is correlated to transcriptional changes involved in muscle protein degradation and catabolism. This hypothesis could be further investigated through analysis of energy stores within muscle in subsequent sampling events via non-lethal biopsy collection. Perturbations in metabolism following dilbit exposure was also suggested by the upregulation of cytochrome c oxidase subunit 1 (mt-co1) and cytochrome c oxidase subunit 2 (mt-co2) which encode for catalytic subunits of cytochrome C oxidase (COX) ( Figure 5) (Capaldi, 1990). COX functions as an aerobic respiratory enzyme and catalyzes the final step of the mitochondrial electron transport chain (Capaldi, 1990;Barrientos et al., 2002). In fish, the enzymatic activity of COX is an indicator of metabolic capacity and a biomarker of metabolic perturbation following exposure to petroleum hydrocarbons (Gagnon, 1999;Cohen et al., 2005). Following waterborne exposures to crude oil wateraccommodated fractions, COX activity in the gills and liver of Australian bass (Macquaria novemaculeata) either increased or decreased compared to controls depending on concentration (14-134 μg/L Σpetroleum hydrocarbons) (Cohen et al., 2005). Similarly, COX activity was altered in the heart and red muscle of Atlantic salmon smolts following subchronic exposure to dilbit, being either increased or decreased relative to controls depending on tissue type and concentration (9.65 and 67.9 μg/L ΣPAC) (Avey et al., 2020). Thus, dysregulation of mt-co1 and mt-co2 in dilbit-exposed lake charr corroborates previous studies suggesting that crude oil induces alterations in COX and, thus, changes in energy metabolism. Although there are many abiotic and biotic factors that may contribute to stress response in the environment, the lack of molecular evidence for energetic changes among baseline years 2017-2018 provides a degree of support that the suggested energy metabolism changes in 2019 were a result of dilbit additions. For example, there were no GO pathways related to intermediary metabolism or energetics among 2017 and 2018, suggesting no molecular changes in energy metabolism year-to-year at baseline conditions. Furthermore, only two genes (g6pd and gpt2) relating to energy metabolism were dysregulated among baseline conditions of the genes identified to be altered following dilbit additions ( Figure 5).

Identification of Biomarkers
The expression patterns of certain genes may be used as biomarkers of exposure and can help researchers quickly assess whether an individual has been exposed to a contaminant of interest. The induction of CYP1A has commonly been used as biomarker of exposure to petroleum hydrocarbons and other organic pollutants (Goksøyr, 1995;Woodin et al., 1997;Whyte et al., 2000) and exposure to dilbit has also been shown to induce CYP1A activity in several fish species (Madison et al., 2015, 2017, Alderman et al., 2017a, 2017bAlsaadi F. M. et al., 2018;McDonnell et al., 2019). In the present study, there was a trend of cyp1a3 upregulation in lake charr following dilbit additions (FDR-adj p = 0.08) ( Figure 5; Supplementary  Table S4). Although not statistically significant at a p = 0.05, a trend in upregulation of cyp1a3 may have biological significance and may be a viable biomarker in mucus, especially at higher PAC concentrations seen following an oil spill. A study exposing Japanese medaka embryos to water accommodated fractions (WAF) of Cold Lake dilbit indicated that cyp1a mRNA synthesis increased monotonically with concentration in hatched fish (Madison et al., 2017). In the present study, ΣPAC levels within surface waters were relatively low at the period of sampling (approximately 0.067 μg/L) and a greater cyp1a response would likely be expected at higher concentrations observed during a spill. Our results are consistent with findings of a previous study in which cyp1a1 and cyp1b1 were the top upregulated transcript in mucus of mahi-mahi following exposure to Deepwater Horizon slick oil (16.55 μg/L ΣPAC) (Greer et al., 2019). In the present study, there was a clear trend of cyp1a3 upregulation which was confirmed by qPCR at relatively low environmental concentrations and confirms that this may be a sensitive biomarker even within mucus of animals.

Applications and Implications
Overall, results of RNA sequencing suggest stress following dilbit exposure, as indicated by upregulation of transcripts involved in cortisol-induced energy metabolism. More specifically, evidence suggests that there was a potential increased reliance on amino acid metabolism and decreased reliance on fatty acid oxidation for energy production in dilbit-exposed lake charr compared to controls. The changes in energy metabolism and utilization support previous studies which have also shown that crude oil and dilbit specifically may induce changes in energy metabolism (Alderman et al., 2018;Avey et al., 2020). Moreover, cyp1a3 was upregulated in lake charr mucus, suggesting a molecular response to dilbit. Further laboratory studies should be conducted to confirm these molecular changes however, coupled with assessment of phenotypic endpoints. Laboratory studies may also be beneficial as potential confounding variables in the field can be eliminated, such as varying environmental variables among sampling years. Along with this, future studies should utilize greater replication for a more definitive results, as this preliminary study utilized 4 individuals per treatment group. As RNA sequencing identified molecular changes related to energy metabolism and PAC Frontiers in Environmental Science | www.frontiersin.org February 2022 | Volume 10 | Article 836640 exposure in lake charr mucus, results of the present study indicate that transcriptomic analysis of epidermal mucus is a promising vector for development of nonlethal oil exposure assessment methods.
There are some realities that should be considered when using molecular analysis of mucus for exposure assessment following an oil spill. Although RNA sequencing is a holistic approach to capturing molecular changes in an organism, use of targeted qPCR for genes such as cyp1a and ahr may be a more rapid and inexpensive method for monitoring exposure following an oil spill. Along with this, mucus from reference individuals will be needed for either RNA sequencing or qPCR analysis. Although use of baseline mucus samples from the population of interest would be ideal to minimize confounding variables, it may not be feasible. Instead, mucus from individuals from a wellcharacterized reference location may be used. Use of other nonlethal endpoints in conjunction with analysis of mucus for target genes of interest may provide a more complete insight on organism health. Examples include analysis of blood for enzymatic activity or use of muscle plugs for contaminant loading assessment. There are also logistical issues that require further research before the method may be applied for spill monitoring to determine magnitude and geographic extent of effects. For example, further research is needed to identify how rapidly the mucus transcriptome responds to external conditions as well as to determine the turnover rate of the mucus transcriptome. Furthermore, the relationship among gene expression in actively metabolizing tissues, such as liver, and the mucus should be identified in order to better understand the physiological significance of differential expression in the mucus. While additional research is needed, use of epidermal mucus is a promising method for oil exposure assessment without the use of lethal sampling.

DATA AVAILABILITY STATEMENT
The datasets presented in this study can be found in online repositories. The names of the repository/repositories and accession number(s) can be found below: https://www.ncbi. nlm.nih.gov/, PRJNA729809.

ETHICS STATEMENT
The animal study was reviewed and approved by the Canadian Council on Animal Care.

AUTHOR CONTRIBUTIONS
NA contributed to study design, collected samples, conducted bioinformatics and statistical analyses, and wrote the manuscript. VP contributed to study design, collected samples, provided field support, and edited the manuscript. LH collected samples, provided field support, and edited the manuscript. LP conducted chemistry analysis. DS contributed to study design, provided field support, and edited the manuscript. All authors contributed final refinements to the manuscript.