Transcriptome Analysis Reveals That Naphthenic Acids Perturb Gene Networks Related to Metabolic Processes, Membrane Integrity, and Gut Function in Silurana (Xenopus) tropicalis Embryos

Naphthenic acids (NAs) are oil-derived mixtures of carboxylic acids and are considered emerging contaminants with the potential to disrupt development of aquatic species. In the Oil Sands Region of Canada, NAs are components of the water released following processing of the bitumen-containing sand. The aim of this research was to identify potential mechanisms of toxicity of NA mixtures. Silurana (Xenopus) tropicalis embryos were raised in water spiked with commercial oil-derived NA extracts (S1 and S2) at a sub-lethal concentration (2 mg/L). The transcriptomic responses of the whole 4-day old embryos following exposure were assessed using a custom oligonucleotide microarray. Both NA mixtures induced embryonic abnormalities that included edema, and cardiac and gut abnormalities. Exposure to NAs also affected morphometric parameters and decreased total length, tail length, and interorbital distance of the embryos. Gene ontology analysis revealed that 18 biological processes, 5 cellular components, and 19 molecular functions were significantly enriched after both S1 and S2 exposures. Sub-network enrichment analysis revealed pathways that were related to phenotypic abnormalities; these included gut function, edema, and cartilage differentiation. Other notable networks affected by NAs included metabolism and cell membrane integrity. In a separate dose-response experiment, the expression of key genes identified by microarray (cyp4b1, abcg2, slc26a6, eprs, and slc5a1) was determined by Real-Time qPCR in S. tropicalis embryos exposed to the commercial NAs and to acid-extractable organics (AEOs) prepared from Oil Sands Process-Affected Water. In general, the RT-qPCR data agreed with the microarray data. In S. tropicalis embryos exposed to the AEOs, the mRNA levels of eprs (bifunctional glutamate/proline-tRNA ligase) and slcs5a1 (sodium/glucose cotransporter 1) were significantly decreased compared to the controls. Such changes are likely indicative of increased edema and disrupted gut function, respectively. These data suggest that NAs have multiple modes of action to induce developmental toxicity in amphibians. Some modes of action may be shared between commercial NAs and AEOs.

Naphthenic acids (NAs) are oil-derived mixtures of carboxylic acids and are considered emerging contaminants with the potential to disrupt development of aquatic species. In the Oil Sands Region of Canada, NAs are components of the water released following processing of the bitumen-containing sand. The aim of this research was to identify potential mechanisms of toxicity of NA mixtures. Silurana (Xenopus) tropicalis embryos were raised in water spiked with commercial oil-derived NA extracts (S1 and S2) at a sub-lethal concentration (2 mg/L). The transcriptomic responses of the whole 4-day old embryos following exposure were assessed using a custom oligonucleotide microarray. Both NA mixtures induced embryonic abnormalities that included edema, and cardiac and gut abnormalities. Exposure to NAs also affected morphometric parameters and decreased total length, tail length, and interorbital distance of the embryos. Gene ontology analysis revealed that 18 biological processes, 5 cellular components, and 19 molecular functions were significantly enriched after both S1 and S2 exposures. Sub-network enrichment analysis revealed pathways that were related to phenotypic abnormalities; these included gut function, edema, and cartilage differentiation. Other notable networks affected by NAs included metabolism and cell membrane integrity. In a separate dose-response experiment, the expression of key genes identified by microarray (cyp4b1, abcg2, slc26a6, eprs, and slc5a1) was determined by Real-Time qPCR in S. tropicalis embryos exposed to the commercial NAs and to acid-extractable organics (AEOs) prepared from Oil Sands Process-Affected Water. In general, the RT-qPCR data agreed with the microarray data. In S. tropicalis embryos exposed to the AEOs, the mRNA levels of eprs (bifunctional glutamate/proline-tRNA ligase) and slcs5a1 (sodium/glucose cotransporter 1) were significantly decreased compared to

INTRODUCTION
Interest in the potential toxicity of oil sands process-affected water (OSPW) emerged in the 1960s after the initiation of oil extraction activities in the oil sands region of Alberta, Canada (Grewer et al., 2010). Crude bitumen from the oil sands contains petroleum mixed with inorganic material (Government of Canada, 2013). Extraction processes require up to 2.5 barrels of water to produce one barrel of crude oil (Natural Resources Canada, 2019). The OSPW is kept in separate tailings ponds to prevent contamination of natural water sources. As of 2015, there was a total surface area of 98 km 2 of tailings pond water (Environment and Parks Alberta, 2015). Frank et al. (2014) compared the chemical profile of the groundwater of the Athabasca River to OSPW and found important similarities, concluding that tailings pond water was leaching out of the ponds and into the surrounding environment.
The toxicity of the OSPW has mainly been attributed to the presence of naphthenic acids (NAs) (MacKinnon and Boerger, 1986). The O 2 family is represented by classical NAs with just one carboxylic acid group or dihydroxy groups, and published data concur that O 2 species are the main toxic and classical NAs in OSPW (Morandi et al., 2015;Hughes et al., 2017;Huang et al., 2018). The term NAs refers to all the carboxylic acids present in crude oil and these are now classified as emerging water contaminants (Richardson and Ternes, 2018). Naphthenic acids have been also detected in sediments after oil spills (Aeppli et al., 2012;Wan et al., 2014), and are also used commercially in a wide range of products such as paint and ink driers, wood and fabric preservatives, fuel additives, emulsifiers and surfactants, and are used in the production of metal salts (Brient et al., 2000). There are only a few published studies that have attempted to address the mechanisms of NA toxicity in aquatic organisms. Frank et al. (2009) proposed narcosis as a mechanism of NA toxicity in Daphnia magna. Later, Marentette et al. (2017) suggested oxidative stress as a possible toxic mechanism of naphthenic acid fraction components (NAFCs) of OSPW affecting walleye (Sander vitreus) embryos. More detailed studies are needed to assess the effects and potential mechanisms of NA toxicity. Considering that there are more than 100 countries that produce petroleum, the toxicity and biological effects that may be caused by diverse NAs is of international concern.
The interest in frogs as a test organism is due to their importance in the food web, their susceptibility to pollutants and endocrine disruptors in water, and as a model for vertebrate development. Specifically, the African Western clawed frog (Silurana (Xenopus) tropicalis) is a diploid model species commonly used in the field of developmental biology and toxicology (Rosenfeld et al., 2017). A sequenced genome and relatively high number of conserved genes compared to humans, in particular, those controlling development and associated with disease (Hellsten et al., 2010), make S. tropicalis an amenable model organism. Moreover, frogs in early developing stages are susceptible to endocrine disruptors, and studies have shown that the exposure of frogs in larval stages to stressors has effects on future morphology, such as body size, locomotor ability, and survivorship (Crespi and Denver, 2005). There have been a few articles published regarding the toxicity of NAs in frogs (Hersikorn and Smits, 2011;Smits et al., 2012;Melvin and Trudeau, 2012a,b;Melvin et al., 2013;Vu et al., 2015). In this study, we report the developmental and transcriptomic responses of S. tropicalis embryos to NAs in order to determine potential mechanisms of toxicity.

Chemicals
Two commonly available commercial NA mixtures were purchased from Sigma-Aldrich (Cat. # 70340; lot numbers BCBC9959V and BCBK0736V), hereafter called Sigma 1 (S1) and Sigma 2 (S2), respectively. The acid extractable organic mixture (AEO) contains organic compounds from the acid extraction of an OSPW sample from Alberta, Canada. Details on the source, handling and preparation of the extract can be found in Gutierrez-Villagomez et al. (2019). The different mixtures were previously described and characterized by electrospray ionization high-resolution mass spectrometry (Gutierrez-Villagomez et al., 2019) and gas chromatography coupled to mass spectrometry . The S1 and S2 are composed largely of the O 2 and O 4 classes. The AEO extract is composed mainly of O 2 , O 3 , O 4 (89.3%) and other species in a lower proportion (10.7%). The O 2 species proportion in the AEO extract was estimated as 58.6%. Human chorionic gonadotropin (hCG) was purchased from Millipore. The salts used in the preparation of the FETAX solution (NaCl, NaHCO 3 , KCl, CaCl 2 , CaSO 4 , and MgSO 4 ) were purchased from Fisher Scientific. L-Cysteine was obtained from Sigma-Aldrich.

Animals and Exposure
Adult frogs originated from the University of Ottawa's S. tropicalis colony. Adult frogs were fed daily with Nasco pellets and kept under a 12 h light/dark cycle. Embryos were obtained by injection of hCG into the posterior lymph sac of adult S. tropicalis as previously described (Langlois et al., 2010). Briefly, the frogs were injected with a priming dose of 12.5 IU of hCG into the posterior lymph sac and moved into glass tanks with FETAX solution [Standard Guide for Conducting the Frog Embryo Teratogenesis Assay-Xenopus (FETAX)] (ASTM, 2012) with a temperature of 21.4 • C and pH of 5.9. After 24 h, the breeding pairs were injected with a boosting dose of 100 IU of hCG and placed in glass tanks with FETAX solution with a temperature of 21.6 • C and a pH of 5.8. The exposures were carried out according to the Frog Embryo Teratogenesis Assay Xenopus (FETAX). FETAX is a rapid test for identifying teratogenic compounds or mixtures (FETAX) (ASTM, 2012). After spawning, the embryos were collected and the jelly coat was removed by gently mixing for 3 min in a 2% w/v L-cysteine solution prepared with FETAX solution. Cleaned embryos were then moved to a tank with FETAX solution. According to Nieuwkoop and Faber (NF) staging system (Nieuwkoop and Faber, 1994), embryos at stages 9-10 were individually selected under a light microscope. Four replicates with 10 embryos were used for each treatment and were placed in 10 cm Petri dishes containing 50 mL of the NA solutions. Solutions of NAs for the different treatments were prepared daily from a stock that was kept at 4 • C and covered from light. The stock solutions from S1 and S2 were prepared using ethanol as a solvent. The final concentration of ethanol in all of the treatments was 0.0025%. The AEO extract contained 8,050 mg/L NAs in a solution of 0.1 N of NaOH as previously reported (Gutierrez-Villagomez et al., 2019). Glassware preparation also followed (Gutierrez-Villagomez et al., 2019): 15 min soaking with hot water and detergent, rinsed twice with reverse osmosis water, later rinsed in a 10% HCl solution, rinsed twice with reverse osmosis water, rinsed with acetone (Fisher Scientific, Certified ACS grade, 99.8%), rinsed twice with reverse osmosis water and dried in an oven at 120 • C for 3 h. The embryos were exposed to S1 and S2 at 2 mg/L as reported (Gutierrez-Villagomez et al., 2019). This concentration was selected as it induces sublethal effects in S. tropicalis embryos. Control embryos were placed in FETAX solution. Treatment solutions were prepared daily from a respective stock solution that was kept at 4 • C and covered from light. Every 24 h the embryos were observed, dead embryos were removed, and the solutions were replaced to ensure concentration and water quality. The exposure was performed in a controlled light/temperature room with a 12:12 light:dark cycle. The water temperature was 25 ± 1 • C. After 96 h of exposure, tadpoles were sampled to measure gene expression profiles. The tadpoles were anesthetized with a solution of 0.5 g/L tricaine-S [MS-222, tricaine methanesulfonate (TMS); Western, Chemical, Inc.]. Digital images of the embryos at 96 h were taken under a light stereomicroscope (Nikon SMZ 1500), with a Nikon DS-Fi1 camera and NIS Elements version 3.22.00 software. Images were used for morphometric measurements. The abnormalities were classified according to the FETAX Atlas of Abnormalities (Bantle, 1991). The care and treatment of animals used in this study were in accordance with the guidelines of the Animal Care Committee of the University of Ottawa and the Canadian Council on Animal Care.

Microarray Analysis
Four S. tropicalis larval samples for the S1 and S2 treatments were used for microarray analysis. Each sample consisted of a pool of 10 larvae to ensure an adequate quantity of RNA. Total RNA extraction was done using RNeasy Mini Kit (Qiagen) with on-column RNase-free DNase treatment (Qiagen). RNA concentrations were determined using NanoDrop-2000 spectrophotometer (Thermo Scientific). Total RNA integrity was assessed using RNA 6000 Nano Assay Kit with the 2100 Bioanalyzer (Agilent Technologies). The RNA integrity number (RIN) for the S. tropicalis samples ranged from 8.0 to 8.8. A custom 4 × 44 K Agilent microarray developed for S. tropicalis was used to identify transcripts affected by the exposure to NA mixtures. The microarray platform is deposited in Gene Expression Omnibus (GPL15626) at NCBI and comprised 32,899 unique probes and 1417 Agilent control features (Langlois and Martyniuk, 2013). Microarray hybridization was performed according to the Agilent One-Color Microarray-Based Gene Expression Analysis protocol using Cyanine 3 (Cy3). For the production of complementary DNA (cDNA), 1 µg total RNA per sample was used. The labeling protocol followed that provided by the manufacturer. Microarrays were hybridized for 17 h and washed according to the Agilent protocol. The Agilent G2600D Microarray Scanner was used for microarray scanning at 5 µm. The Agilent Feature Extraction Software (v. 11.0.1.1) was used to extract the signal intensities from the TIFF microarray images. Microarray data have been deposited in the Gene Expression Omnibus (accession # GSE118736) and are MIAME compliant 1 . Intensity data were imported into JMP R Genomics (Version 6.0) and data were normalized using quantile normalization. The microarray limit of detection was 2.5, as determined according to the lowest point in the standard curve and the Agilent negative controls. Thus, probes that had a signal intensity below 2.5 were filtered to a value of 2.5. Following the removal of all the control spots, the differentially expressed genes were identified using a one-way ANOVA test with a false-discovery rate of 5% (FDR = 0.05).

Bioinformatics
Any probe that showed a change with respect to the control (p < 0.05) prior to FDR adjustment in one or both treatments of NA mixtures was used to broaden the analysis of the transcriptomic responses. Clustering analysis was performed to determine if there was evidence of separation in transcriptomes among treatments. Two-way hierarchical clustering of differentially expressed probes was performed using the Fast Ward algorithm. Rows were centered to a mean of zero prior to clustering and were also scaled to a variance of one. Gene set enrichment on gene ontology (GO) terms was determined using the parametric analysis of gene set enrichment (PAGE) algorithm, which is a two-sided z-score for GO categories (Kim and Volsky, 2005). Sub-network enrichment analysis (SNEA) was performed using Pathway Studio 9.0 (Ariadne, Rockville, MD, United States) and ResNet 9.0. The SNEA uses known relationships based on expression, binding and common pathways between genes to build networks focused around gene hubs. The GenBank ID was then used for mapping human homologs in Pathway Studio; 4,572 genes were successfully mapped. Enrichment p-value cut-off was set at p < 0.05. The analysis used the function "highest fold change, best p-value" for duplicate probes. This bioinformatics method leverages the entire dataset regardless of p-value and builds a distribution based on fold change to statistically test for enrichment of processes.

RNA Extraction, cDNA Synthesis, and RT-qPCR for Microarray Validation
The samples obtained in a previous dose-response experiment (Gutierrez-Villagomez et al., 2019) were used for RT-qPCR validation of the microarray. The concentrations used for the S1 treatments for RT-qPCR analysis were 0, 2, 4, 6, and 8 mg/L. For S2, the treatments analyzed were 0, 2, 6, 8, and 12 mg/L. RNA was isolated using the RNeasy mini kit (Qiagen) as described in the manufacturer's protocol with on-column RNase-free DNase treatment (Qiagen). Before cDNA synthesis, the concentration and quality of all samples were assessed using NanoDrop-2000 spectrophotometer (Thermo Scientific) and on an agarose gel. The quality of the RNA was indicated by the presence of two defined bands. The top band represents 28S ribosomal RNA (rRNA) and the second band represents the 18S rRNA (Aranda et al., 2012). Total cDNA was prepared using Maxima first strand cDNA synthesis kit for RT-qPCR (Thermo Scientific).
Real-time quantitative polymerase chain reaction (RT-qPCR) with SYBR green dye technology was used to validate relative gene expression. Based on the microarray results, six genes were selected for further analysis (Table 1); the ribosomal protein L8 (rpl8) gene was selected as the reference gene because it did not significantly change with exposure to NAs (Table 1) and is a ribosomal gene often used for RT-qPCR normalization (Dhorne-Pollet et al., 2013). Primers were designed using Primer-BLAST 2 and synthesized by Integrated DNA Technologies (Supplementary Table S1). To confirm the amplification of the regions of interest, the PCR products were cloned into the pGEM R -T Easy vector (Promega, Madison, WI, United States) and sequenced using 3730 DNA Analyzer (Applied Biosystems). The Rotor-Gene SYBR Green kit (Qiagen) and Rotor-Gene Q (Qiagen) were used to amplify and detect the transcripts of interest. The thermal cycling parameters were as suggested by the manufacturer; an activation step at 95 • C for 5 min, followed by 40 cycles of 95 • C denaturation step for 5 s and one primer annealing temperature of 60 • C for 10 s. After 40 cycles, a melt curve was performed over a range of 60-95 • C with increments of 1 • C to ensure a single amplified product. The concentration of each primer in the RT-qPCR reactions was 1 µM. The efficiency of all RT-qPCR reactions was 94-107% (100.7 ± 3.4%) and the coefficient of determination (R 2 ) was ≥ 0.99 (0.995 ± 0.004). Data were analyzed using the Rotor-Gene Q Series software (Qiagen). Analysis of outliers was performed using the ROUT method in GraphPad Prism 6 ( Motulsky and Brown, 2006). The relative standard curve method was used to calculate relative mRNA abundance between samples, normalized using NORMA-GENE algorithm (Heckmann et al., 2011) and then presented as fold change of gene expression from replicates (n = 5-6; assayed in triplicate) for each group.

Statistical Analysis
Frequencies of abnormalities were analyzed as transformed data using the arcsine square root transformation. To assess data normality and homogeneity of variance, Shapiro-Wilk's test and Levene's test were performed, respectively. Data that failed the normality and/or the equal variance tests were transformed (log10). One-way ANOVA and Tukey's post hoc analyses were performed on normally distributed data. Nonparametric Kruskal-Wallis followed by Student-Newman-Keuls (SNK) tests were performed on data that did not pass normality tests. The significance level was set at α = 0.05. Statistical analysis was performed using SPSS V21 IBM and Sigma Plot 12.0.

Effects of Exposure to Commercial NAs
The exposures of S. tropicalis embryos to NAs at 2 mg/L resulted in a significant reduction in total length (TL), snoutvent length (SVL), tail length (TaL), and interorbital distance (IOD) (p < 0.05) (Figures 1A-C and Supplementary Figure S1). Following exposure to NAs, there was a significant increase of all abnormalities assessed. The most common were cranial abnormalities, edema, gut abnormalities, heart abnormalities, and eye abnormalities (Figures 1A-C and Supplementary Figure S1). These types of teratogenic effects have been described in detail previously in S. tropicalis exposed to commercial extracts and AEOs in a dose-response fashion (Gutierrez-Villagomez et al., 2019). Edema and heart abnormalities have also been described in larval zebrafish (Danio rerio) exposed to AEOs extracted from oil sands from the Daqing oil exploring area (Wang et al., 2015) and larval fathead minnows (Pimephales promelas) exposed to NA fraction components (NAFCs) of OSPW from Alberta, Canada (Marentette et al., 2015). The 96 h LC 50 estimates for S. tropicalis were 10.4, 11.7, and 52.3 mg/L for S1, S2, and the AEOs, respectively (Gutierrez-Villagomez et al., 2019). The 96 h EC 50 estimates based on frequencies of developmental abnormalities in S. tropicalis were 2.1, 2.6, and 14.2 mg/L for S1, S2, and the AEOs, respectively. In fathead minnow embryos, an EC 50 for hatching success of 2.5 and 7.5-18.5 mg/L have been reported for a Sigma extract and AEOs, respectively (Marentette et al., 2015). These observations indicate that fish and frogs might have a similar sensitivity with respect to embryotoxicity in response to commercial NAs and AEOs, however, direct comparisons between taxa should be conducted to test this hypothesis.
FIGURE 1 | Photographs of Silurana (Xenopus) tropicalis indicating the main teratogenic effects induced by exposure to NAs, control (A), S1 (B), and S2 (C). Scale bar is equal to 1 mm in (A-C). Some of the observed abnormalities are indicated with arrows. Hierarchical clustering of gene expression data using gene probes that were differentially expressed (p < 0.05) in S. tropicalis embryos exposed to two mixtures of NAs (S1 and S2) during early embryonic development versus control samples (C1-4) in (D). In panel (D) columns represent individual genes and rows represent a treatment replicate. Each cell in the matrix represents the expression level of a gene feature in an individual replicate. The red or green color in cells reflects relative high or low expression levels, respectively, as indicated in the scale bar.
Frontiers in Marine Science | www.frontiersin.org

Gene Expression Profiling Following NA Exposure
The microarray analysis performed on embryos of S. tropicalis exposed to S1 revealed 6,074 differentially expressed genes versus control (p < 0.05), of which 2,873 were up-regulated and 3,201 were down-regulated. A total of 894 genes passed the correction for multiple hypothesis testing (FDR < 0.05). This includes 526 up-regulated and 368 down-regulated genes. For S2-exposed animals, the analysis revealed 6,789 differentially expressed genes versus control (p < 0.05). Of these, 3,395 were up-regulated and 3,394 were down-regulated by the exposure to S2. A total of 1,177 genes passed the FDR test (FDR < 0.05). Of these, 620 were up-regulated and 557 down-regulated. There were 3,656 common genes affected by S1/S2 (p < 0.05). Of these, 564 genes passed the FDR test (FDR < 0.05). The higher number of genes significantly affected (p < 0.05) by the exposure to S2 compared to S1, indicates that S2 has a higher impact at the transcriptomic level of S. tropicalis. These results highlight that the NA composition has a direct impact on the toxicity of a mixture. The S1 and S2 are composed largely of the O2 and O4 classes. However, the composition and proportion of the components in the two mixtures is different (Gutierrez-Villagomez et al., 2019). We observed that the EC 50 and EC 20 for S1 are significantly lower than for S2 (Gutierrez-Villagomez et al., 2019). In this manuscript, we observed that the two mixtures affected the transcriptome of the S. tropicalis embryos in a different manner. Thus, it is clear that subtle differences in chemical composition can have different transcriptional effects linked with differing teratogenic outcomes.

Hierarchical Analysis
The hierarchical clustering analysis of the differential expression (p < 0.05) data showed two main clades that mark the clear separation of the control samples and the two NA treatments (S1 and S2). The second clade is composed only of expression profiles from embryos exposed to the NA extracts. Treatments S1 and S2 are sister groups, indicating that these two NA mixtures affect the transcriptome of S. tropicalis in a similar, but nonidentical manner. The different replicates are located in their respective cluster and there are no replicates mixed between groups, indicating a relatively strong transcriptome response to S1 and S2 ( Figure 1D).

Gene Set Enrichment of Gene Ontology (GO) Terms
In total, 64 biological processes, 25 cellular components, and 52 molecular functions were significantly enriched after S1 exposure (p < 0.05) (Supplementary Table S2). After S2 exposure, 77 biological processes, 29 cellular components, and 57 molecular functions were significantly enriched (p < 0.05) (Supplementary Table S3). The number of significantly enriched processes in common between S1 and S2 for each of biological process (p:), cellular component (c:) and molecular function (f:) were 18, 5, and 19, respectively (p < 0.05) ( Table 2). The differences between the number of GO terms significantly enriched with the exposure to S1 and S2 likely reflects the differences in the chemical composition of the two mixtures (Gutierrez-Villagomez et al., 2019). The significant enrichment of GO terms, such as respiratory electron transport chain, respiratory chain, adenylyl-sulfate kinase activity, and pyruvate kinase activity indicates that NAs affect metabolism and that this could be part of the mechanism of toxicity of NAs ( Table 2). "Protein import into mitochondrial inner membrane" and "mitochondrial intermembrane" GO terms were also significantly enriched (p < 0.05). These results support the hypothesis of Frank et al. (2009) who proposed that narcosis as part of the toxic mechanism of NAs. Furthermore, visual perception, structural constituent of eye lens, actin cytoskeleton, and cytoskeleton organization are GO terms significantly enriched (p < 0.05) after the exposure to S1 and S2. This may be associated with the multiple abnormalities in S. tropicalis induced by S1 and S2, such as eye and morphometric abnormalities (Figures 1A-C).

Sub-Network Enrichment Analysis (SNEA)
Using SNEA, a total of 75 cell process and 42 disease pathways were affected by S1 exposure (p < 0.05; Supplementary  Table S4), and 66 cell process pathways and 46 disease pathways were perturbed in the embryos following S2 exposure (p < 0.05; Supplementary Table S5). There were 38 cell process pathways (p < 0.05; Table 3) and 14 disease pathways (p < 0.05; Table 4) in common as perturbed by S1 and S2. Some pathways that were significantly affected by exposure to S1 (Figure 2A) and S2 ( Figure 2B) included networks related to metabolism, such as liver metabolism, cholesterol metabolism, retinoid metabolism, retinoic acid metabolism, steroid metabolism, arachidonic acid metabolism, drug metabolism, and xenobiotic clearance. This corroborates the GO analysis that also indicated that metabolic processes were significantly enriched following exposure. Other pathways detected as significantly enriched were those related to membrane integrity (Figures 3A,B), such as membrane depolarization, Cl − transport, Na + influx co-transport, ion transport, anion transport, fluid transport, and HCO − 3 transport. This analysis further supports the hypothesis that toxicity of NAs is at least partially through narcosis, defined as the disruption of cell membrane function (Newman and Clements, 2007;Li et al., 2017). Others refer to this as baseline toxicity and have linked narcosis to changes in membrane potentials (Escher et al., 2002). The morphological analysis showed that NAs significantly induce gut abnormalities (Supplementary Figure S1J) in S. tropicalis embryos, and the lack of gut coiling was observed in some of the tadpoles (Figures 1A-C). Importantly, gene networks related to gut function (Figure 4), such as digestion, intestinal absorption, acid secretion, gastrointestinal system digestion, bile acid secretion, and biliary flow were significantly disrupted by exposure to NAs. Supporting this is the identification of disease pathways (Table 4) such as constipation and malabsorption syndromes, both related to gut functions. Other pathways related to the abnormalities that were observed in S. tropicalis (Figures 1A-C and Supplementary Figure S1) included cartilage differentiation and edema (Tables 3, 4).  (biological process (p:), cellular component (c:) and molecular function (f:) altered after the exposure of S. tropicalis to both S1 and S2 using Parametric Analysis of Gene Set Enrichment (PAGE) analysis (p < 0.05).

GO term
Frequency S1 S2

RT-qPCR Validation of Differentially Expressed Genes
Five genes were selected for microarray validation and for further study based on the following characteristics: (1) both exposure to S1 and S2 significantly altered their expression; (2) the transcripts had to be annotated and have a human homolog identified; and (3) the transcripts had to appear in the pathway analysis ( Table 1). The cytochrome P450 4B1b (cyp4b1) gene encodes a 3 | Sub-network enrichment analysis (SNEA) for common cellular processes in S. tropicalis exposed to S1 and S2 mixtures of NAs. All gene set seeds that had > 10% overlap with the annotated pathways are listed (p < 0.05). The total number of neighburs refers to the total number of known entities in a sub-network.
monooxygenase enzyme that catalyzes reactions involved in the metabolism of cholesterol, steroids, and xenobiotics (Baer and Rettie, 2006). This gene forms part of a metabolism network that was significantly disrupted by NAs, specifically arachidonic acid metabolism (Figures 2A,B). Arachidonic acid is a fatty acid involved in eye and brain development (Uauy et al., 2001). The solute carrier family 26 member 6 (slc26a6) gene encodes an anion transporter protein and regulates acid-base homeostasis (Chernova et al., 2005;Freel et al., 2009). This gene is part of networks that were disturbed by the exposure to NAs, in particular, intestinal absorption and Cl − transport (Figures 3A,B). Furthermore, Slc26a6-null mice have defective intestinal oxalate secretion and a higher incidence of formation of stony concretions in the urinary tract (Jiang et al., 2006). The ATP-binding cassette sub-family G member 2 (abcg2) gene encodes a membrane transporter (Allikmets et al., 1998)  All gene set seeds that had >10% overlap with the annotated pathways are listed (p < 0.05). The total number of neighbors refers to the total number of known entities in a sub-network.
involved in the gene network of gut function, specifically in the pathway of intestinal absorption (Figures 4A,B). The bifunctional glutamate/proline-tRNA ligase (eprs) gene encodes for the enzyme that charge tRNA with glutamate and proline (Young et al., 2016). It was selected because the SNEA showed that this gene is related to edema in S. tropicalis embryos (see Table 1). The sodium/glucose cotransporter 1 (slc5a1) gene encodes a sodium-dependent glucose transporter (Wright et al., 1994) and it is involved in several networks affected by NAs, such as liver metabolism, membrane depolarization, Na + influx co-transport, and intestinal absorption (Figures 2-4). Lee et al. (1994) reported that slc5a1 is highly expressed in the intestine and in lower levels in liver in rat, where its function is the uptake of glucose (Turk et al., 1991). Also, the mutation of the slc5a1 gene can induce severe diarrhea and dehydration in humans (Lam et al., 1999). These five genes appear to be associated with the significant increase (p < 0.05) of phenotypic abnormalities in S. tropicalis embryos, such as gut abnormalities, eye abnormalities, and edema ( Figures 1A-C and Supplementary Figure S1). The myod1 gene was selected because it changed significantly with both S1 and S2 treatments and is linked to myogenesis and muscle differentiation (Weintraub et al., 1991). Therefore, it may be associated with the significant decrease of TL, SVL, TaL, and IOD of S. tropicalis exposed to NAs (Figures 1A-C and Supplementary Figure S1). In general, the expression obtained by microarray analysis and RT-qPCR were similar in direction of change (Figures 5, 6 and Table 1). The genes cyp4b1 (Figures 5A,B), abcg2 (Figures 5E,F), slc26a6 (Figures 5G,H), eprs (Figures 5J,K), and slc5a1 (Figures 6A,B) were affected in a concentration-response manner and the fold changes were comparable between the S1 and S2 treatment groups. However, some discrepancies were noted for myod1. The microarray analysis showed a significant down-regulation of myod1 by exposure to S1 and S2 (Table 1), while the RT-qPCR indicated an upregulation (Figures 6G,H).
There are clear differences in the chemical composition of the commercial NA extracts versus the AEOs from OSPW (Gutierrez-Villagomez et al., 2019). However, we had hypothesized that the toxic mechanism would be similar because all extracts contain high levels of O 2 species of NAs, albeit, less in the AEOs. As well, exposures produce the same types of abnormalities. To test this hypothesis, samples of the S. tropicalis exposed to AEO concentrations of 0, 2, 12, 24, and 48 mg/L were analyzed by RT-qPCR for the same genes described for microarray validation. The expression of cyp4b1 (Figure 5C), abcg2 (Figure 5F), and slc26a6 ( Figure 5I) was not affected in a concentration-response manner in AEO-treated animals. After extraction, the AEOs were dissolved in a solution of 0.1N NaOH. The solvent control did not have a significant effect on morphometric measurements and frequency of abnormalities (Gutierrez-Villagomez et al., 2019). Likewise, Marentette et al. (2015) reported no effect on the frequency of abnormalities up to concentrations of 50 mM NaOH in fathead minnow embryos. However, the RT-qPCR results showed that NaOH at 0.001 N reduced the expression of many of the genes analyzed by PCR (Figures 5, 6, AEO column). It is nevertheless important to use and report the effects of NaOH with the toxicity assessment of AEOs. The expression level of eprs ( Figure 5L) and slc5a1 ( Figure 6C). significantly decreased at 48 mg/L AEOs compared to control and NaOH control (p < 0.05). As eprs and slc5a1 are part of networks related to edema, membrane integrity, and metabolism and gut function ( Table 1) and frequency of edema and gut abnormalities was significantly higher in samples FIGURE 2 | Networks associated with metabolism indicating altered gene regulation following exposure of S. tropicalis to 2 mg/L of S1 (A) and S2 (B). Sub-network enrichment analysis (SNEA) was performed using Pathway Studio 9.0. Pink and red indicate that the transcript abundance was increased and green indicates that the transcript abundance was decreased (p < 0.05). The abbreviations are the official gene symbol of the transcripts of the respective pathway (Supplementary  Tables S6, S7). FIGURE 3 | Networks associated with membrane integrity indicating altered gene regulation following exposure of S. tropicalis to 2 mg/L of S1 (A) and S2 (B). Sub-network enrichment analysis (SNEA) was performed using Pathway Studio 9.0. Pink and red indicate that the transcript abundance was increased and green indicates that the transcript abundance was decreased (p < 0.05). The abbreviations are the official gene symbol of the transcripts of the respective pathway (Supplementary Tables S8, S9).
Frontiers in Marine Science | www.frontiersin.org FIGURE 4 | Networks associated with gut function indicating altered gene regulation following exposure of S. tropicalis to 2 mg/L of S1 (A) and S2 (B). Sub-network enrichment analysis (SNEA) was performed using Pathway Studio 9.0. Pink and red indicate that the transcript abundance was increased and green indicates that the transcript abundance was decreased (p < 0.05). The abbreviations are the official gene symbol of the transcripts of the respective pathway (Supplementary Tables S10, S11).
FIGURE 5 | RT-qPCR analysis for cyp4b1 (A-C), abcg2 (D-F), slc26a6 (G-I), and eprs (J-L), mRNA in S. tropicalis exposed to different concentrations of S1 (first column), S2 (second column), and AEOs (third column). Data were normalized and defined as fold change relative to control. Different letters indicate groups that are significantly different (p < 0.05).
Frontiers in Marine Science | www.frontiersin.org 13 August 2019 | Volume 6 | Article 533 FIGURE 6 | RT-qPCR analysis for slc5a1 (A-C), rpl8 (D-F), and (G-I) myod1 mRNA in S. tropicalis exposed to different concentrations of S1 (first column), S2 (second column), and AEOs (third column). Data were normalized and defined as fold change relative to control. Different letters indicate groups that are significantly different (p < 0.05).
exposed to AEOs (Gutierrez-Villagomez et al., 2019), there appears to be an important relationship between the phenotypic abnormalities in S. tropicalis embryos and altered eprs and slc5a1 levels.

CONCLUSION
We report the first comprehensive transcriptomic profiling of the effects of NAs in early S. tropicalis embryos. A diverse number of GO terms and pathways were significantly enriched by exposures to commercial extracts S1 and S2, suggesting that there are several mechanisms involved. The results obtained in the microarray analysis and comparison to morphometric abnormalities in S. tropicalis exposed to S1 and S2 were highly similar. Metabolism and membrane integrity GO terms and pathways were significantly enriched, indicating that exposure to NAs disrupts metabolism and the functionality of cell membranes in S. tropicalis. These data support the hypothesis that narcosis is at least one of the likely mechanisms of toxicity of NAs. Other GO terms and pathways that were also significantly enriched were related to abnormalities observed in the S. tropicalis larvae after the early exposure to NAs. These included gut function, cartilage differentiation, actin cytoskeleton, cytoskeleton organization, and edema. While AEOs contain less of the highly toxic NAs in the O2 class, they still cause similar abnormalities and cause changes in the expression of specific genes. The rank-order potency was S1 > S2 > AEOs (Gutierrez-Villagomez et al., 2019) based on LC50 (classical mortality measurements) and EC50 (incidence of total abnormalities in S. tropicalis embryos). Here, we show that both commercial NA mixtures and AEOs extracted from OSPW at concentrations found in the environment affect some of the same genes. What is new in our study is that we were able to screen the activities of the commercial NAs, then test whether the OSPW extract has the same effects on the consistently regulated genes-stepping toward the very first potential biomarkers for NAs. This is especially the case for eprs and slc5a1, which may be respectively linked to induction of edema and gut abnormalities in exposed embryos. Future research should focus on specific tissues, considering that gut abnormalities were one of the most relevant and unique phenotypic changes observed, and that was strongly supported by the transcriptomic profiling data. It will also be important to determine the relative contributions of AEOs and polyaromatic compounds on the effects of whole OSPW and other effluents associated with petroleum production.

DATA AVAILABILITY
The datasets generated for this study can be found in GEO, GSE118736.

ETHICS STATEMENT
This study was carried out in accordance with the Canadian Council on Animal Care and approved by the University of Ottawa Animal Care and Veterinary Services.

AUTHOR CONTRIBUTIONS
JG-V and VT designed the research project. JG-V conducted the research and wrote the manuscript. CM helped with bioinformatics. LX helped with RT-qPCR analysis. VL and CM designed the microarray platform. BP provided the financial support. All authors read, edited, and approved the final manuscript.

FUNDING
Grant support for JG-V was provided by CONACyT Mexico. The support of NSERC Strategic Grants Program (STPGP 463041-14 to VT and JB), University of Ottawa Research Chair program, and Environment and Climate Change Canada through the Canada-Alberta Joint Oil Sands Monitoring Program for VT are also acknowledged with appreciation.